1 Introduction

Jets, as sprays of collimated hadrons resulting from fragmentation of high-energy partons produced in hard scatterings, are multipurpose tools to explore various properties of the strong interaction, via measurements of the strong coupling constant [1, 2], heavy-flavour production [3], and the global fit analysis of parton distribution functions (PDF) [4]. Jet processes are also pivotal to address fundamental questions such as the validity of factorisation theorems [5,6,7], or existence of gluon saturation due to nonlinear QCD dynamics [8]. Jets, produced abundantly at LHC energies due to the large centre-of-mass energy and high luminosity, provide high precision tests of QCD [9,10,11,12]. Moreover, measuring the “jet quenching” phenomenon [13] observed for jets produced in heavy-ion collisions offers an insight into how high-momentum partons interact with the medium created in the collisions. Through this interaction one probes QCD at high energy densities and temperatures, where the strongly-interacting matter enters the guark–gluon plasma (QGP) phase. The observed properties of this matter are consistent with a strongly-coupled, low-viscosity fluid of quarks and gluons [14, 15].

Recent measurements in high-multiplicity pp and pA collisions [16] exhibit several collective effects qualitatively similar to the ones observed in AA collisions for various observables. Such collective behaviour encompasses long-range (\(|\Delta \eta |\!\!\!>\!\!\!2\)) near-side (\(|\Delta \varphi |\!\!\!\approx \!\!\!0\)) two-particle angular correlations, known as the “ridge” [17,18,19,20,21]; enhanced yield of charged- or identified-particle production in high-multiplicity events with respect to the reference using minimum-bias (MB) charged- or pion-particle production [22,23,24,25,26]; elliptic flow of heavy-flavour hadrons [27]; and enhanced baryon production at intermediate transverse momentum (\(p_\mathrm {T}{\,\sim \,} 3{\mathrm {GeV}}/c\)) [28]. With the lack of experimental observations of jet quenching effect with present accuracy in small collision systems [29, 30], the measurements aforementioned raise the difficult but intriguing question of whether these observations arise similarly to heavy-ion collisions, namely from the formation of a hot and dense fluid-like medium, or rather involve other physical mechanisms [31, 32]. Several theoretical approaches and models have been put forward to explain these QGP-like effects in small systems while accounting for the absence of jet quenching [33,34,35,36,37], such as multiple parton interactions (MPI) [38], string shoving [39], or rope hadronisation [40]. However, these models cannot explain the measured non-zero elliptic flow at high \(p_\mathrm {T}\) from two-particle correlations [19], which is usually attributed to in-medium path-length dependent energy loss [41,42,43].

To deepen our understanding of the mechanisms that are at play in high-multiplicity collisions of small systems, the multiplicity dependence of the charged-jet production has been studied in pp collisions at a centre-of-mass energy \(\sqrt{s}=13\ {\mathrm {TeV}}\) by ALICE. Charged-particle jets are reconstructed from tracks measured at midrapidity using the anti-\(k_\mathrm {T}\) clustering algorithm [44] with jet resolution parameters R ranging from 0.2 to 0.7. The measured inclusive jet cross sections are compared to model calculations, allowing us to test the relative importance of various mechanisms at play in particular hadronisation and MPI to which the lowest transverse momentum jets are most sensitive.

The event activity is quantified by the charged-particle multiplicity measured by the ALICE V0 detector at forward rapidity, in order to minimise the autocorrelations between the event selection and the measured observable. By taking advantage of the largest integrated luminosity collected so far by the ALICE experiment, this work extends previous measurements [9, 12, 45] to higher collision energy, broader jet kinematic range, larger jet radii (up to \(R{=}0.7\)), and higher event activities. A complementary insight is thereby provided to similar CMS [46] and ALICE measurements focusing on the soft sector based on transverse spherocity [47] and two-particle correlations [48]. Thanks to these new measurements, stronger constraints are placed on models that describe fundamental mechanisms of particle production in hadronic collisions.

The results presented in this article test whether the seemingly universal pattern of the self-normalised production of hard probes as a function of event activity, emerging from the ALICE measurements of \(\mathrm {J/\psi }\) [49,50,51], \(\mathrm D\) [52], and \(\mathrm B\)-meson production [53], holds also for jets. While this self-normalised hard probe production pattern could be ascribed to mere collision geometry up to a certain multiplicity, a new regime of high-density gluon configurations is expected to set in at the highest multiplicities [54].

The article is organised as follows: Sect. 2 describes the ALICE detectors, the experimental data samples, and Monte-Carlo simulations used for this analysis; Sect. 3 discusses the multiplicity selection and the jet reconstruction methods; Sect. 4 outlines the unfolding corrections and the estimation of the systematic uncertainties; Sect. 5 presents the results; and finally, Sect. 6 gives a summary and outlook.

2 Experimental setup and data sample

ALICE is the dedicated heavy-ion experiment at the LHC. A detailed description of the ALICE apparatus can be found in Ref. [55]. In the following, only the detector components used in the data analysis presented in this article are described.

The ALICE apparatus comprises a central barrel (pseudorapidity coverage \(|\eta |<0.9\) over the full azimuth) situated in a uniform \(0.5\ \mathrm {T}\) magnetic field along the beam axis (z), which is supplied by the large L3 solenoid magnet [56]. The central barrel contains a set of tracking detectors: a six-layer silicon inner tracking system (ITS) surrounding the beam pipe, and a large-volume (\(5\ \mathrm {m}\) length, \(0.85\ \mathrm {m}\) inner radius and \(2.8\ \mathrm {m}\) outer radius) cylindrical Time Projection Chamber (TPC). The first two layers of the ITS are instrumented with high-granularity Silicon Pixel Detectors (SPD), followed by two layers composed of Silicon Drift Detectors (SDD), and finally, the two outer layers are made of double-sided Silicon micro-Strip Detectors (SSD). In the forward and backward rapidity regions, a pair of plastic scintillator counters called V0A and V0C are positioned on each side of the interaction point, covering pseudorapidity ranges \({2.8<\eta <5.1}\) and \({-3.7<\eta <-1.7}\), respectively. The V0 system provides the interaction trigger for the whole experiment, and is further used to suppress machine-induced background events [57].

The measurement was based on the data from pp collisions at a centre-of-mass energy of \(\sqrt{s}=13~\,{\mathrm {TeV}}\) collected between 2016 and 2018. During this period, MB events were selected online using the high purity V0-based MB trigger [58], which required a charged-particle signal coincidence in the V0A and V0C arrays. The visible cross section satisfying the MB trigger requirement was determined in a van der Meer scan [59, 60]. The integrated luminosity of the used sample, measured with V0, is \(8.12 \pm 0.16~\mathrm nb^{-1}\) for 2016, \(10.67 \pm 0.29~\mathrm nb^{-1}\) for 2017, and \(13.14 \pm 0.27~\mathrm nb^{-1}\) for 2018, respectively. The luminosity uncertainty was evaluated to be 1.6 % by taking into account the correlations during the combination of the samples [61]. For the offline analysis, further event selection was made by requiring a primary vertex position within \(\pm 10\ \mathrm {cm}\) in the longitudinal direction around the nominal interaction point to ensure full geometrical acceptance in the ITS for \({|\eta |<0.9}\). Pile up, i.e. the average number of simultaneous interactions per bunch crossing, was maintained well below unity through beam separation in the horizontal plane. Residual pile-up events were rejected based on a multiple vertex finding algorithm using tracking information from the SPD. The final corresponding data sample consists of \(2.2\times 10^9\) events after the trigger and offline selection.

Reconstructed tracks with transverse momenta larger than \(0.15\ {\mathrm {GeV}}/c\) were selected in the analysis. The track selection criteria was optimised for good momentum resolution and minimal contamination from secondary particles, as described in Refs. [11, 12]. To ensure a uniform \((\eta ,\varphi )\) distribution in the regions where the SPD was inactive, tracks with no hit in either of the two SPD layers were constrained to the primary vertex. The tracking efficiency estimated from a full detector simulation amounts to \(80\%\) for \(p_\mathrm {T}>0.4\ {\mathrm {GeV}}/c\), decreasing to \(60\%\) at \(0.15\ {\mathrm {GeV}}/c\). The transverse momentum resolution is better than \(3\%\) for tracks with \(p_{\mathrm T}\) below \(1\ {\mathrm {GeV}}/c\), and increases linearly up to \(10\%\) at \(p_\mathrm {T}=100\ {\mathrm {GeV}}/c\).

The response of the ALICE detector to produced particles was evaluated using GEANT3 [62]. Based on this response, measured distributions were corrected for instrumental effects, see Sect. 4.1. In this paper, the default Monte Carlo (MC) event generator used for comparison with measurements is the PYTHIA 8.125 general-purpose Leading Order (LO) MC generator (hereafter referred to as PYTHIA8) with the Monash-2013 set of tuned parameters (tune) [63] for the underlying event (UE) and NNPDF2.3 LO PDF set [64], while MPI and Colour Reconnection (CR) models being enabled. Furthermore, in order to reduce the large theoretical uncertainties affecting the computations at LO in perturbative QCD, like the residual dependence of the unphysical factorisation and renormalisation scales, jet production at next-to-leading order (NLO) accuracy was obtained within the POWHEG framework [65]. Unlike pure fixed order calculations, POWHEG interfaces NLO calculations with PYTHIA8 parton showers to generate exclusive final states. The particle level outputs of such simulations were then directly compared with the experimental data, which were corrected for instrumental effects.

Fig. 1
figure 1

Scaled V0M distribution which is used to determine the forward multiplicity classes in pp collisions at \(\sqrt{s} = 13\) TeV. The colour shaded areas represent V0M multiplicity classes obtained from real data. The PYTHIA8 distribution is shown with the open black markers

Table 1 Average charged-particle pseudorapidity densities at midrapidity \(\left\langle \mathrm {d}N_{\mathrm{ch}}/\mathrm {d}\eta \right\rangle \) from data for inclusive events and different V0M multiplicity classes [66]

3 Multiplicity selection and jet reconstruction

3.1 Multiplicity selection

In order to study the multiplicity dependence of inclusive charged-particle jet production, the MB sample was divided into event classes based on the “V0M amplitude” that is proportional to the total number of charged particles passing through the V0A and V0C detectors. The distribution of the self-normalised V0M amplitude from data and the PYTHIA8 event generator is shown in Fig. 1. The distribution is normalised to its average value, \(\left\langle \mathrm V0M~amplitude\right\rangle \), to reduce the sensitivity of the multiplicity percentile determination on the amplitude. PYTHIA8 MC does not reproduce the measured multiplicity distribution, as was already reported in Ref. [66]. To reduce the potential model-dependent bias, corrections of the multiplicity dependent jet yields were done using a data-driven method instead of pure MC samples, which is discussed in Sect. 4.1.

The event classes used in the analysis and the corresponding midrapidity charged-particle densities for experimental data are summarised in Table 1. The multiplicity classes were defined in terms of percentile intervals of experimental V0M amplitude/\(\left\langle \mathrm {V0M~amplitude}\right\rangle \) as shown in Fig. 1. The average charged-particle multiplicity densities in MB pp collisions and for events of a given multiplicity class were obtained by integrating the corresponding fully corrected \(p_{\mathrm T}\) spectra given in Ref. [66]. When comparing the data to MC predictions, the multiplicity percentile was calculated from data and MC using their respective self-normalised distribution accordingly in order to minimise the difference observed in the V0M amplitude distribution. The 0–1% range corresponds to the highest multiplicity class (I), while the 60–100% interval corresponds to the lowest multiplicity class (VIII).

3.2 Jet reconstruction

Jets were reconstructed from tracks with \(p_{\mathrm{T,track}} > 0.15\ {\mathrm {GeV}}/c\) and \(\left| \eta _{\mathrm{track}} \right| <0.9\), using the anti-\(k_{\mathrm{T}}\) sequential recombination algorithm [44] from the FastJet package [67]. The jet transverse momenta were calculated using boost-invariant \(p_{\mathrm{T}}\) recombination scheme as a scalar sum of their constituent transverse momenta. Jet resolution parameters were varied in the range from \(R=0.2\) to 0.7. The pseudorapidity of the reconstructed jets was limited to \(\left| \eta _{\mathrm{jet}} \right| <0.9 - R\) to ensure they remain in the fiducial acceptance [68]. The transverse momentum range of the inclusive charged-particle jet spectra spans from 5 to \(140\ {\mathrm {GeV}}/c\). The spectra measured in V0M multiplicity classes have the upper limit \(100\ {\mathrm {GeV}}/c\). The cross section of inclusive jet production was measured as a function of \(p_{\mathrm T}\) considering the van der Meer minimum bias visible cross section mentioned in Sect. 2. For the multiplicity dependence study, the per-event jet yield was measured as a function of multiplicity classes defined by the V0M percentile intervals. In addition, the integrated jet yield was calculated as a function of the charged-particle density \(\mathrm {d}N_{\mathrm{ch}}/\mathrm {d}\eta \), self-normalised similarly to other ALICE measurements [49,50,51,52,53].

The measured jets are inevitably affected by the UE activity originating from MPI, fragmentation of beam remnants, as well as initial- and final-state radiation [69]. In pp collisions, the UE effect on jet measurements is rather limited [12]. However, since the UE contribution depends on event multiplicity, the measured jets were not affected in the same way for events falling in different multiplicity classes. In order to perform fair comparisons between different multiplicity intervals, the results presented in this paper include the UE subtraction.

The UE contribution to the charged-particle jet \(p_{\mathrm T}\) was estimated event-by-event using the same approach as in previous measurements in pp [70] and p–Pb collisions [29, 30]. The background density \(\rho _{\mathrm{ch}}\) was determined using the \(k_{\mathrm{T}}\) algorithm [71] with a fixed radius of 0.2, taking into account only jets containing at least one physical particle, while removing the two \(k_{\mathrm{T}}\) clusters of highest transverse momentum to limit the impact of the jet signal on the underlying event estimation. The background density \(\rho _{\mathrm{ch}}\) is then used to subtract the average background from each jet in the same event: \(p_{\mathrm{T}}^{\mathrm{corr}} = p_{\mathrm{T}}^{\mathrm{raw}} - \rho _{\mathrm{ch}} \cdot A\), where A is the jet area.

During the subtraction of the average UE background of each jet, the local background fluctuations smear the reconstructed jet transverse momentum. To study jet-by-jet fluctuations of the background, the random cone (RC) method was used [72]. In this method, cones of radius R positioned at random (\(\eta , \varphi \)) coordinates within the detector acceptance (fiducial region) were generated in each event. The sum of the charged-particle \(p_{\mathrm{T}}\) in a given cone was then compared to the expected average background obtained from \(\rho _{\mathrm{ch}}\) as follows:

$$\begin{aligned} \delta p_{\mathrm{T}}^{\mathrm{RC}} = \sum ^{\mathrm{RC}}p_{\mathrm{T,track}}-\rho _{\mathrm{ch}}\pi R^{2}, \end{aligned}$$
(1)

where the sum runs over track \(p_{\mathrm{T}}\) inside the random cone, and \(\pi R^{2}\) is the area of the random cone. The width of \(\delta p_{\mathrm{T}}^{\mathrm{RC}}\) is a measure of the momentum smearing due to local background fluctuations [73]. To minimise the influence of signal jets on the \(\delta p_{\mathrm{T}}^{\mathrm{RC}}\) distribution, a minimum distance from the random cone to the two highest momentum jets (leading jets) in the event was required. The \(\delta p_{\mathrm{T}}^{\mathrm{RC}}\) distribution, obtained for different cone radii R in events with excluded leading jets, are shown in Fig. 2 a). It clearly shows stronger background fluctuations with increasing jet radius, as expected due to the larger number of particles within the jet cone.

Fig. 2
figure 2

a Comparison of the \(\delta p_{\mathrm{T}}\) distribution obtained for different random cone radii (\(R=0.2, 0.4, 0.7\)). b Comparison of the \(\delta p_{\mathrm{T}}\) distribution with the RC method (including and excluding two leading jets) and the track embedding approach for \(R=0.4\). c Comparison of measured \(\delta p_{\mathrm{T}}\) distribution using RC method without leading jets for \(R=0.4\) in different multiplicity classes

An alternative embedding method was used to quantify the background fluctuations. In this procedure, a probe track was embedded into an event [74]. The azimuthal angle of the probe track was required to be perpendicular to the jet (\(\varphi _{\mathrm{track}}^{\mathrm{emb}} = \varphi _{\mathrm{ch\,jet}} + \pi /2\)) while retaining its \(\eta \) value (\(\eta _{\mathrm{track}}^{\mathrm{emb}} = \eta _{\mathrm{ch\,jet}}\)). The transverse momentum of the probe track (\(p_{\mathrm{T}}^{\mathrm{emb}}\)) was uniformly chosen between 0 and \(100\ {\mathrm {GeV}}/c\). After embedding the probe track into the event, the jet finding algorithm was relaunched with the same background subtraction method as described above. The embedded \(\delta p_{\mathrm{T}}^{\mathrm{emb}}\) was evaluated in a similar way to Eq. 1 after removing the momentum of the embedded probe track:

$$\begin{aligned} \delta p_{\mathrm{T}}^{\mathrm{emb}} = p_{\mathrm{T,ch\,jet}}^{\mathrm{raw,emb}}-\rho _{\mathrm{ch}}A_{\mathrm{ch\,jet}}^{\mathrm{emb}} -p _{\mathrm{T}}^{\mathrm{emb}}, \end{aligned}$$
(2)

where \(p_{\mathrm{T,ch\,jet}}^{\mathrm{raw,emb}}\) and \(A_{\mathrm{ch\,jet}}^{\mathrm{emb}}\) are the transverse momentum and area of the reconstructed jet with the embedded probe track.

The background fluctuations determined by the random cone and the track embedding method are presented in Fig. 2b. While the distributions obtained from the different methods show very similar negative tails, the tail on the right-hand side of the distribution caused by real jets is much less pronounced when a minimum distance to two leading jets is required. Therefore, the \(\delta p_{\mathrm{T}}\) distribution from the RC method without two leading jets was used as default to build up the background fluctuation response, while the track embedding method was used for the assessment of the systematic uncertainty on the background fluctuation estimate. Figure 2c compares the \(\delta p_{\mathrm{T}}\) distribution in different multiplicity intervals with the RC method when avoiding that the cone overlaps with the two leading jets. The figure suggests that the magnitude of local background fluctuations grows with imposed multiplicity bias, as expected.

4 Corrections and systematic uncertainties

The differential charged-particle jet cross section after UE subtraction was corrected for jet energy scale smearing due to local background fluctuations and detector effects. The resulting cross section unfolded to charged particle level allows for a direct comparison with the theoretical predictions. These unfolding correction approaches are explained in Sect. 4.1. The corresponding systematic uncertainties are presented in Sect. 4.2.

4.1 Unfolding corrections

The measured jet momentum is affected by background fluctuations and a variety of instrumental effects, including finite momentum resolution and track reconstruction efficiency. The jet spectrum was corrected for these effects to obtain the particle-level jet spectrum through an unfolding procedure [72]. First, the particle-level true jets were constructed from the PYTHIA8 (Monash 2013) event generator [63], by selecting only those stable charged particles defined as all particles with a mean proper lifetime larger than \(1\ \mathrm {cm}/c\), and excluding the decay products of these particles [75]. Next, jets were reconstructed at detector level from tracks coming from MC particles propagated through the GEANT3 model of the ALICE setup. The corresponding jet energy scale residual and resolution were estimated to be about 20% for jet transverse momentum larger than \(10~{\mathrm {GeV}}/c\) with jet resolution parameter \(R=0.4\).

Based on the geometrical matching between the corresponding jets at the particle level and detector level, the detector level response matrix was constructed. The response matrix employed in the unfolding procedure was a combination of the detector response matrix and the background fluctuation response matrix. The background fluctuation response matrix utilises the \(\delta p_{\mathrm{T}}\) distribution obtained from data, based on distributions shown in Fig. 2. This data-driven approach ensures that the response matrix used for unfolding reflects the accurate multiplicity dependence. Finally, the measured jet spectrum was unfolded using the combined response matrix which corrects for background fluctuations as well as detector effects.

In this analysis, the default unfolding method was based on the singular value decomposition (SVD) approach as implemented in the RooUnfold package [76]. The regularisation parameter k, which suppresses high-frequency variations in the unfolded result, was selected by examining the d-vector distribution [77]. In addition to the SVD unfolding, Bayesian unfolding [78] was also used for systematic uncertainty evaluation. Consistent results were obtained between both unfolding methods. To validate the unfolding process and identify potential biases, closure tests that compare the unfolded distributions to the particle-level true distributions were performed. The consistency of the unfolding procedure was also checked by folding the solution to detector level and comparing it to the measured raw spectrum. In both cases, no significant difference was found.

4.2 Systematic uncertainties

The sources of systematic uncertainty that affect the jet measurement were divided into the following categories: tracking efficiency and \(p_{\mathrm{T}}\) resolution, unfolding corrections, background fluctuations, contamination from secondary particles and normalisation. Additionally, in the case of multiplicity dependent measurements, a systematic uncertainty due to multiplicity estimation was also considered. All uncertainties from these sources were considered as uncorrelated, except the uncertainties of the unfolding corrections that come from several correlated sources. Therefore, the uncertainty of the unfolding category is calculated by varying each source and calculating the RMS of all variations. Then, all systematic uncertainty categories were treated separately and their respective contributions we added in quadrature.

Table 2 summarises the contributions to the systematic uncertainties for the inclusive jet cross section in a few selected jet \(p_{\mathrm{T}}\) and R bins. Table 3 summarises the systematic uncertainties for multiplicity-dependent jet production in a few selected jet \(p_{\mathrm{T}}\) and multiplicity intervals for a jet reconstruction resolution parameter of \(R = 0.4\).

The total systematic uncertainty increases with jet \(p_{\mathrm{T}}\) and R, and its evolution is similar for all jet radii and multiplicity intervals that are not listed here. In the following, the individual sources of systematic uncertainties listed in Tables 2 and 3 are briefly described.

The dominant source of systematic uncertainty comes from track reconstruction efficiency since it directly affects the jet energy scale and resolution. The systematic uncertainty on tracking efficiency was estimated to be 3% based on variations of the criteria used in the track selection [12]. To evaluate the effect of this uncertainty on the measured jet spectra, a new detector response matrix was computed by generating a PYTHIA simulation that accounts for a 3% reduction of tracking efficiency, and then used in the unfolding procedure. The difference between the spectra corrected with the default response matrix and with the response matrix obtained with the decreased tracking efficiency was taken as a systematic uncertainty. The relative uncertainty on the jet spectra caused by tracking efficiency increases slowly with increasing jet \(p_{\mathrm{T}}\) and has a weak multiplicity dependence.

The relative systematic uncertainty on track momentum resolution was estimated from the study of the invariant mass distributions of \(\Lambda \) and \(\mathrm {K^0_s}\) as a function of \(p_{\mathrm T}\) and amounts to 20% [79]. This track \(p_{\mathrm T}\) resolution uncertainty was then propagated to the corrected jet spectra with a similar method as used for the tracking efficiency uncertainty evaluation. The resulting uncertainty from track \(p_{\mathrm T}\) resolution is about 2%.

Table 2 Summary of systematic uncertainties on inclusive jet cross section for three selected jet transverse momentum bins with different resolution parameters
Table 3 Summary of systematic uncertainties on jet production for three jet transverse momentum bins with different multiplicity percentile intervals for jets with resolution parameter \(R = 0.4\)

In order to assign the uncertainty arising from the unfolding corrections, several variations are considered. By default, the reconstructed jet spectra were unfolded using a detector response matrix and prior spectrum obtained based on events generated by the PYTHIA8 generator with the Monash 2013 tune [63]. The prior spectrum is used as initial guess of the true spectrum. The dependence on MC event generator was quantified by comparing the spectra unfolded using the response matrix and prior from the default generator with those unfolded with response from the EPOS generator with LHC tune [80]. The resulting uncertainty is of the order of 3%. Second, the SVD unfolding method [77], as the default approach used in this paper, was regularised by the choice of parameter k for each cone radius parameter. To estimate the related systematic uncertainty, the regularisation parameter was varied by \(\pm 2\) around the optimised value. The unfolded results were stable against regularisation parameter variations, and the corresponding systematic uncertainty is negligible. To validate the unfolding procedure, the SVD unfolded spectra were compared with the results obtained with the Bayesian unfolding method and the remaining differences were taken as an uncertainty. In addition to the above variations, the bin boundary migration uncertainty was also evaluated by changing the boundaries of the input spectra and the response matrix during the unfolding process. The uncertainties discussed above were then with the RMS calculated for all variations and referred to as the unfolding systematic uncertainty in Tables 2 and 3.

The systematic uncertainty due to the background fluctuation estimation was quantified by comparing unfolded spectra obtained with \(\delta p_{\mathrm{T}}\) distributions using the method of RC without two leading jets (default) and the track embedding method as discussed in Sect. 3.2.

The difference between both corrected jet spectra was assigned as background fluctuation uncertainty. A systematic uncertainty on the background density \(\rho _{\mathrm{ch}}\) measurement was estimated to be 5%, resulting in a 2% uncertainty on the UE-subtracted jet cross section at \(p_{\mathrm T}= 5\ {\mathrm {GeV}}/c\), and smaller for higher jet transverse momentum. This uncertainty is highest at low \(p_{\mathrm {T,jet}}\) in high multiplicity events.

By default the multiplicity percentiles were determined from the measured distribution of V0M amplitude in data as listed in Table 1. As shown in Fig. 1, PYTHIA8 MC cannot reproduce such multiplicity distribution, which is mainly attributed to a limited description of the UE [69, 70, 81] and inaccurate description of secondary particles and magnetic field during particle transport in the detector simulation. To prevent such multiplicity differences from being propagated as multiplicity dependent results, the unfolding corrections use the data-driven approach with the background fluctuation response matrices taken from the corresponding event multiplicity class based on data directly, see Fig. 2. This matrix was multiplied with the instrumental matrix obtained from PYTHIA minimum bias events. To account for the multiplicity estimation uncertainty, a response matrix obtained from pure MC simulation was also used, where the multiplicity intervals and the background fluctuations were both extracted from PYTHIA8 generator. The systematic uncertainty was assigned based on the comparison of the unfolded spectra obtained from the default analysis and this variation. Such uncertainty reaches 5.7% for low-\(p_{\mathrm T}\) jets in the highest multiplicity class, and decreases in the lower multiplicity percentile intervals. The uncertainty is independent of the jet R since the multiplicity estimation is at the event level, and the ratio of jet yields of different R is independent of multiplicity.

Fig. 3
figure 3

Inclusive charged-particle jet cross sections in pp collisions at \(\sqrt{s} = 13\) TeV using the anti-\(k_{\mathrm{T}}\) algorithm for different jet resolution parameters R from 0.2 to 0.7, with UE subtraction. Statistical uncertainties are displayed as vertical error bars. The total systematic uncertainties are shown as solid boxes around the data points

Secondary charged particles are mostly produced by weak decays of strange particles (\(\mathrm {K^0_s}\), \(\Lambda \), etc.), photon conversions, hadronic interactions in the detector material, and decays of charged pions. Contamination from such secondary charged particles was significantly reduced by a requirement on the maximum distance of closest approach (DCA) of the tracks to the primary vertex. Therefore, the systematic uncertainties due to secondaries were estimated by varying the DCA threshold of track selection, resulting in a jet \(p_{\mathrm T}\) scale uncertainty of 0.5% [9, 45], which turns into a jet cross section uncertainty of about 3%.

A systematic uncertainty on the integrated luminosity measurement of 2% [61] was assigned to the inclusive jet cross section measurement as a normalisation uncertainty which consequently does not affect the ratios of the measured cross section spectra.

When calculating the systematic uncertainties on the ratios of jet spectra, each uncertainty source was varied simultaneously both for numerator and denominator, and the ratio was calculated using the varied spectra. The resulting difference between the varied spectra and the nominal one is taken as the uncertainty on the ratio for that given uncertainty source. This results in a significant reduction of the correlated uncertainties from the cancellation between the numerator and the denominator [72]. The remaining relative difference from each source was added in quadrature.

The statistical uncertainties on the jet production ratio were also treated carefully between numerator and denominator. To avoid the statistical correlations, the total event sample was divided into two parts for the calculation of the numerator and denominator, respectively. The resulting statistical uncertainty on the ratio remains smaller than the systematic uncertainty.

5 Results

5.1 Inclusive jet cross section measurements

The fully-corrected inclusive charged-particle jet cross sections after UE subtraction in pp collisions at \(\sqrt{s}=13~{\mathrm {TeV}}\) are shown in Fig. 3 as a function of jet \(p_{\mathrm T}\) for jet resolution parameters ranging from \(R = 0.2\) to \(R = 0.7\) and pseudorapidity ranges \(|\eta _{\mathrm{jet}}|<0.9 - R\). The choice of R changes the relative strength of perturbative and non-perturbative (hadronisation and UE) effects on the jet transverse momentum distribution [82]. To be consistent with the multiplicity dependent results, all figures presented hereafter are obtained with UE subtraction, while the same measurements without UE subtraction are listed in Appendix A.1.

Fig. 4
figure 4

Inclusive charged-particle jet cross sections in pp collisions at \(\sqrt{s} = 13\) TeV with UE subtraction. Data for different jet resolution parameters R varied from 0.2 to 0.7 are compared to LO and NLO MC predictions. The statistical uncertainties are displayed as vertical error bars. The systematic uncertainties on the data are indicated by shaded boxes in the top panels and shaded bands drawn around unity in the bottom panels. The red lines in the ratio correspond to unity

Fig. 5
figure 5

Ratio of charged-particle jet cross section for resolution parameter \(R = 0.2\) to other radii \(R = X \), with X ranging from 0.3 to 0.7, after UE subtraction. Data are compared with LO (PYTHIA) and NLO (POWHEG+PYTHIA8) predictions as shown in the bottom panels. The systematic uncertainties of the cross section ratios from data are indicated by solid boxes around data points in the upper panel and shaded bands around unity in the mid and lower panels. No uncertainties are shown for theoretical predictions for better visibility

Fig. 6
figure 6

Comparison of charged-particle jet cross section ratio with UE subtraction in pp collisions at \(\sqrt{s} = 5.02\) [12], 7 [45], and 13 \({\mathrm {TeV}}\) and in p–Pb collisions at \(\sqrt{s_\mathrm {NN}} = 5.02 \ {\mathrm {TeV}}\) [29]. Results are \(\sigma (R=0.2)/\sigma (R=0.4)\) (a), and \(\sigma (R=0.2)/\sigma (R=0.6)\) (b)

Fig. 7
figure 7

Charged-particle jet yields in different V0M multiplicity percentile intervals for resolution parameters R varied from 0.2 to 0.7 in pp collisions at \(\sqrt{s} = 13\ {\mathrm {TeV}}\). Statistical and total systematic uncertainties are shown as vertical error bars and boxes around the data points, respectively

Fig. 8
figure 8

Ratio of charged-particle jet yield measured in different multiplicity classes with respect to that in MB events as a function of \(p_{\mathrm{T}}\) for different resolution parameters R from 0.2 to 0.7. Statistical and total systematic uncertainties are shown as vertical error bars and boxes around the data points, respectively

Fig. 9
figure 9

Ratios of charged-particle jet spectra with \(R = 0.2\) to that with other jet resolution parameters R from 0.3 to 0.7, shown in different V0M multiplicity classes. Statistical and systematic uncertainties are shown as vertical error bars and boxes around the data points, respectively

Fig. 10
figure 10

Comparison of jet spectra ratios of \(R = 0.2\) to other radii \(R = 0.3\) (a), 0.5 (b), 0.7 (c) in MB events and in three multiplicity intervals (0–1%, 10–15% and 60–100%). Statistical and systematic uncertainties are shown as vertical error bars and boxes around the data points, respectively. Results for other radii can be found in Appendix Fig. 19

Fig. 11
figure 11

Comparison of jet spectra ratios of \(R = 0.2\) to \(R =0.3\) (a), 0.5 (b), 0.7 (c) in three multiplicity intervals (0–1%, 10–15% and 60–100%) and compared with PYTHIA8 simulations. Statistical and systematic uncertainties are shown as vertical error bars and boxes around the data points, respectively. Results for other radii can be found in Appendix Fig. 20

Figure 4 compares the inclusive charged-particle jet cross sections with predictions from the PYTHIA8 and POWHEG MC event generators after UE subtraction, with the same selections and background subtraction procedure applied. The ratios of the MC simulations to ALICE data are shown in the bottom panels of Fig. 4. The POWHEG MC provides a better description of the data within uncertainties for \(p_\mathrm {T,jet}^{\mathrm{ch}}\gtrsim 20\ {\mathrm {GeV}}/c\). Nevertheless, large deviations occur for jet transverse momenta below \(20\ {\mathrm {GeV}}/c\) where POWHEG overestimates the data. Such deviation increases with the jet R. A similar enhancement of POWHEG with respect to the data is also observed at 7 TeV [9], where the implementation of MPI in PYTHIA shows a significant effect on the low \(p_{\mathrm T}\) jet yield when coupled with POWHEG.

Figure 5 shows the inclusive jet cross section ratios for jets reconstructed with a resolution parameter of \(R = 0.2\) to other resolution parameters \(R = 0.3\) to 0.7. The observable defined by the ratio of inclusive jet cross sections relates directly to the relative difference between jet \(p_{\mathrm T}\) distributions when using different resolution parameters and therefore provides insights into the angular dependence of jet fragmentation.

This observable is also less sensitive to experimental systematic uncertainties since the correlated uncertainty on the numerator and denominator spectra are largely cancelled in the ratio. Consequently, the comparisons between data and model predictions provide better precision than those for inclusive spectra. In order to compare the ratios within the same jet pseudorapidity range, the ratios were studied for jet \(|\eta _{\mathrm{jet}}| < 0.2\), which coincides with the fiducial jet acceptance for the largest resolution parameter studied (\(R = 0.7\)). Statistical correlations between the numerator and denominator of the jet cross section ratios were removed by using exclusive subsets events for their respective assessments. The measured ratios were compared with PYTHIA and POWHEG calculations in Fig. 5. Both predictions give a reasonable description of the data for high-\(p_{\mathrm T}\) jets within 10%, although they fail to describe the low-\(p_{\mathrm T}\) region, especially for large resolution parameters, where non-perturbative and UE contributions become large. Even though PYTHIA8 overestimates the jet yields (see Fig. 4), the jet production ratio can still be well described by PYTHIA8 MC.

Figure 6 shows the ratio of the charged-particle jet cross section with different R values for (a) \(R=0.2 / R=0.4\) and (b) \(R=0.2 / R=0.6\) in pp collisions at \(\sqrt{s} = 5.02\) [12], 7 [45], 13 \({\mathrm {TeV}}\), and p–Pb collisions at \(\sqrt{s_\mathrm {NN}} = 5.02 \ {\mathrm {TeV}}\) [29]. These results, which are in good agreement within uncertainties, show a similar increase of the jet cross section ratio as a function of jet \(p_{\mathrm T}\), as expected from the stronger collimation of high-\(p_{\mathrm T}\) jets. No significant energy nor collision species dependence is observed within uncertainties.

5.2 Multiplicity dependence of jet production

The jet production yields measured in different V0M multiplicity intervals as a function of jet \(p_\mathrm {T}\) for different resolution parameters R varied from 0.2 to 0.7 in pp collisions at \(\sqrt{s} = 13\ {\mathrm {TeV}}\) are shown in Fig. 7. A higher (lower) jet yield is observed in higher (lower) multiplicity classes. To better investigate this multiplicity dependence, the ratios of jet spectra from multiplicity classes and with MB events (Fig. 3) are presented in Fig. 8. The charged-particle jet yield ratio in the highest event multiplicity class (0–1%) is about 10 times higher than in the MB case, independent of the jet resolution parameter R, while in the lowest event class (60–100%), it amounts to only about 10% of the MB yield. Furthermore, such ratio has a weak \(p_\mathrm {T}\) dependence, except for the very low \(p_\mathrm {T}\) region. This indicates that jet production changes with event activity, but the slope of the spectrum stays similar to the one measured in MB events.

Figure 9 shows the ratios of the \(R = 0.2\) jet spectrum to other radii for different multiplicity classes. To better understand the multiplicity dependence of the jet spectra ratios, Fig. 10 compares these ratios observed in three selected multiplicity intervals (0–1%, 10–15% and 60–100%) to the ones measured in MB events for (a) \(R = 0.2/0.3\), (b) 0.2/0.5, and (c) 0.2/0.7. The ratios are consistent with the ones obtained in the MB case (Fig. 5) for small jet radii. At larger jet radii, a hint of ordering of the jet production ratios with event multiplicity is observed. It is more pronounced for large radii (\(R = 0.2/0.7\)) and low \(p_{\mathrm{T}}\). However, with the current systematic uncertainty on data, it is difficult to draw final conclusions on such dependence. Similar behaviour is observed in MC simulations as shown in Fig. 11. The MC predictions tend to underestimate the data and this discrepancy increases with jet radius. However, within the current experimental systematic uncertainties, there is no clear indication of multiplicity dependence for jet yield ratios.

The \(p_{\mathrm T}\)-integrated (\(5\le p_{\mathrm{T}}< 100\ {\mathrm {GeV}}/c\)) jet yields and the average transverse momentum of charged-particle jets as a function of the self-normalised charged-particle multiplicity are shown in Fig. 12 for different resolution parameters R from 0.2 to 0.7. Both jet yields and the average jet \(p_{\mathrm{T}}\) increase with multiplicity, the increase is more evident at larger R. Jets with \(R = 0.2\) exhibit very weak dependence of their \(\left\langle p_{\mathrm T}\right\rangle \) on multiplicity, indicating that jets reconstructed with small radii are dominated by the leading particle inside in the jet.

Fig. 12
figure 12

Integrated jet yields (a), and \(\left\langle p_{\mathrm{T}} \right\rangle \) (b) of jets with \(5\le p_{\mathrm{T,jet}}^{\mathrm{ch}}< 100\ \mathrm{GeV}/c\) as a function of self-normalised charged-particle multiplicity for different resolution parameters R varied from 0.2 to 0.7, with the charged-particle multiplicities provided in Ref. [66]. Statistical and systematic uncertainties are shown as vertical error bars and boxes around the data points, respectively

Fig. 13
figure 13

Self-normalised jet yields as a function of the self-normalised charged-particle multiplicity for different resolution parameters R varied from 0.2 to 0.7 in different jet \(p_{\mathrm{T}}\) intervals: \(5\le p_{\mathrm{T,jet}}^{\mathrm{ch}}< 7\ \mathrm{GeV}/c\) (a), \(9\le p_{\mathrm{T,jet}}^{\mathrm{ch}}< 12\ \mathrm{GeV}/c\) (b), \(30\le p_{\mathrm{T,jet}}^{\mathrm{ch}}< 50\ \mathrm{GeV}/c\) (c), and \(70\le p_{\mathrm{T,jet}}^{\mathrm{ch}}< 100\ \mathrm{GeV}/c\) (d). The charged-particle multiplicities are taken from Ref. [66]. Statistical and systematic uncertainties are shown as vertical error bars and boxes around the data points, respectively

Fig. 14
figure 14

Comparison of self-normalised jet yields as a function of the self-normalised charged-particle multiplicity in four selected jet \(p_{\mathrm T}\) intervals (\(5\le p_{\mathrm{T,jet}}^{\mathrm{ch}}< 7\ \mathrm{GeV}/c\), \(9\le p_{\mathrm{T,jet}}^{\mathrm{ch}}< 12\ \mathrm{GeV}/c, 30\le p_{\mathrm{T,jet}}^{\mathrm{ch}}< 50\ \mathrm{GeV}/c\), and \(70\le p_{\mathrm{T,jet}}^{\mathrm{ch}}< 100\ \mathrm{GeV}/c\)) for a given jet radii: \(R = 0.3\) (a), \(R = 0.5\) (b), \(R = 0.7\) (c) between data and PYTHIA8 predictions, with the charged-particle multiplicities provided in Ref. [66]. Statistical and systematic uncertainties are shown as vertical error bars and boxes around the data points, respectively. Results for other radii can be found in Appendix Fig. 21

Figure 13 presents the jet yield ratios in different multiplicity percentiles with respect to MB events as a function of self-normalised charged-particle multiplicity. The ratios are shown for four selected jet \(p_{\mathrm T}\) bins (\(5\le p_{\mathrm{T,jet}}^{\mathrm{ch}}< 7\ \mathrm{GeV}/c\), \(9\le p_{\mathrm{T,jet}}^{\mathrm{ch}}< 12\ \mathrm{GeV}/c\), \(30\le p_{\mathrm{T,jet}}^{\mathrm{ch}}< 50\ \mathrm{GeV}/c\), and \(70\le p_{\mathrm{T,jet}}^{\mathrm{ch}}< 100\ \mathrm{GeV}/c\)), and for resolution parameters \(R = 0.2\) – 0.7. The jet yield ratios increase with multiplicity for all resolution parameters and \(p_{\mathrm T}\) intervals. No significant dependence of the jet yields with the jet resolution parameter R is seen.

To explore the \(p_{\mathrm T}\) dependence of the normalised jet production as a function of self-normalised charged-particle multiplicity, Fig. 14 shows the self-normalised jet yields as a function of the self-normalised multiplicity in four selected jet \(p_{\mathrm T}\) intervals for resolution parameter \(R = 0.3\), 0.5 and 0.7. The PYTHIA8 predictions are also compared against data.

The jet production ratios measured at midrapidity increase with multiplicity in a similar way to the results presented in earlier publications for identified particles when using forward multiplicity V0 estimator [49,50,51]. The increase is weaker for the lowest jet \(p_{\mathrm T}\) in the highest multiplicity interval. In general, PYTHIA8 simulations predict the overall increasing trend, however, the absolute magnitude is overestimated by the PYTHIA8 MC, especially in the highest multiplicity interval.

6 Summary

The inclusive charged-particle jet production cross sections measured with transverse momentum from \(5\ {\mathrm {GeV}}/c\) to \(140\ {\mathrm {GeV}}/c\) in pp collisions at \(\sqrt{s} = 13\ {\mathrm {TeV}}\) have been reported. The measurements were performed using the anti-\(k_\mathrm {T}\) jet finding algorithm with different resolution parameters R varied from 0.2 to 0.7 and the pseudorapidity range \(|\eta _{\mathrm{jet}}|<0.9 - R\). The inclusive charged-particle jet cross sections were compared to LO PYTHIA and NLO POWHEG pQCD calculations. As expected, a better agreement between data and MC is observed for the NLO predictions, although the NLO prediction overestimates the jet yield below \(20\ {\mathrm {GeV}}/c\). The cross section ratios for different resolution parameters were also studied. These ratios increase with jet \(p_\mathrm {T}\) and saturate at the high end of the jet \(p_\mathrm {T}\) range, indicating a stronger collimation for high-momentum jets.

The multiplicity dependence of jet production using different resolution parameters has also been studied. A higher (lower) jet yield is observed in higher (lower) multiplicity classes. Jet production in different multiplicity intervals compared to MB has a weak \(p_\mathrm {T}\) and jet resolution parameter dependence. Furthermore, the self-normalised jet production yields and average jet \(p_\mathrm {T}\) as a function of the self-normalised charged-particle multiplicity have been measured. The integrated jet yields and \(\left\langle p_{\mathrm T}\right\rangle \) in the integrated \(p_{\mathrm T}\) interval between 5 and 100 \(\mathrm{GeV}/c\) increase with the self-normalised charged-particle multiplicity. No strong dependence of jet \(p_\mathrm {T}\) and the resolution parameter R are observed except at low transverse momentum in the highest multiplicity percentile interval. A similar multiplicity dependence has also been reported for prompt D mesons in p-Pb collisions at \(\sqrt{s_{\mathrm{NN}}} = 5.02\ {\mathrm {TeV}}\) and non-prompt J/\(\psi \) (from B hadron decays) production in pp collisions at \(\sqrt{s} = 7\ {\mathrm {TeV}}\) when using a forward multiplicity estimator. Current MC event generators can only predict the rising trend but cannot describe the absolute yields, especially in the highest multiplicity class.

Fig. 15
figure 15

Inclusive charged-particle jet cross sections in pp collisions at \(\sqrt{s} =\) 13 TeV using the anti-\(k_{\mathrm{T}}\) algorithm for different resolution parameters R varied from 0.2 to 0.7, without UE subtraction. Statistical uncertainties are displayed as vertical error bars. The total systematic uncertainties are shown as solid boxes around the data points

Fig. 16
figure 16

Inclusive charged-particle jet cross sections in pp collisions at \(\sqrt{s} = 13\) TeV without UE subtraction and compared to LO and NLO MC predictions with different resolution parameters R varied from 0.2 to 0.7. The statistical uncertainties are displayed as vertical error bars. The systematic uncertainties on the data are indicated by shaded boxes in the top panels and shaded bands drawn around unity in the bottom panels. The red dashed lines in the ratio correspond to unity

Fig. 17
figure 17

Ratio of charged-particle jet cross section for resolution parameter \(R = 0.2\) to other radii \(R = X \), with X ranging from 0.3 to 0.7, without UE subtraction, and the comparison of calculations from LO (PYTHIA) and NLO event generators (POWHEG+PYTHIA8). The systematic uncertainties of the cross section ratios from data are indicated by solid boxes around data points in the upper panels, and shaded bands around unity in the lower panels. No uncertainties are shown for theoretical predictions for better visibility

The measurements presented in this paper provide further insight into the interplay between soft particle production and hard processes. Detailed comparisons of models with data will help to elucidate the relationship between jet production mechanisms and high-multiplicity events in small systems, particularly at LHC energies.