Fragmentation of jets containing a prompt J$/\psi$ meson in PbPb and pp collisions at $\sqrt{s_\mathrm{NN}} =$ 5.02 TeV

Jets containing a prompt J$/\psi$ meson are studied in lead-lead collisions at a nucleon-nucleon center-of-mass energy of 5.02 TeV, using the CMS detector at the LHC. Jets are selected to be in the transverse momentum range of 30 $\lt$ $p_\mathrm{T}$ $\lt$ 40 GeV. The J$/\psi$ yield in these jets is evaluated as a function of the jet fragmentation variable $z$, the ratio of the J$/\psi$ $p_\mathrm{T} $ to the jet $p_\mathrm{T}$. The nuclear modification factor, $R_\mathrm{AA}$, is then derived by comparing the yield in lead-lead collisions to the corresponding expectation based on proton-proton data, at the same nucleon-nucleon center-of-mass energy. The suppression of the J$/\psi$ yield shows a dependence on $z$, indicating that the interaction of the J$/\psi$ with the quark-gluon plasma formed in heavy ion collisions depends on the fragmentation that gives rise to the J$/\psi$ meson.


Introduction
Dissociation of quarkonium states in nucleus-nucleus collisions is one of the best studied signatures of the formation of the quark-gluon plasma (QGP), a deconfined state of quarks and gluons. Although other nuclear effects have been identified, it is widely accepted that at least part of the suppression of the various quarkonium states in central lead-lead (PbPb) collisions, at the collision energies probed at the CERN LHC, is indeed coming from Debye-like screening of heavy-quark pairs in the QGP, as anticipated in the seminal paper by Matsui and Satz [1]. This picture may be probed with quarkonia produced at rest. With the emergence of competing mechanisms, however, it becomes interesting to study the momentum dependence of the nuclear suppression. In particular, this is true for regeneration [2], wherein quarkonia may be created from heavy quarks produced independently. This effect is expected to become more relevant with increased collision energy, as more heavy-quark pairs are produced. This explains the relatively low nuclear suppression at low transverse momentum (p T ) that is observed at the LHC [3], as compared to data at lower collision energies [4,5]. Independent of these nuclear modifications, the interpretation of quarkonium results, and heavy-flavor results in general, is typically based on the assumption that quarkonia are formed at early times compared to the formation time of the QGP (on the order of 1 fm/c).
However, the formation time estimate of quarkonia is based on general arguments rather than on a detailed calculation. Despite decades of theoretical developments, models generally are not able to describe the entirety of the quarkonium data. In particular, they are not able to simultaneously describe the polarization [6] and the p T -differential cross section [7]. A recent measurement by the LHCb Collaboration that looked at hadrons produced at small angles with respect to prompt J/ψ mesons (those not produced in b-hadron decays) in proton-proton (pp) collisions [8] gives some new insight into this puzzle. The observable is the J/ψ-jet fragmentation function, which corresponds to the distribution of z, the ratio of the J/ψ p T to the p T of the jet into which it is clustered. For jets in pseudorapidity 2.5 < η < 4.0 and p T > 20 GeV, LHCb observes that prompt J/ψ mesons are produced with far more in-jet associated hadroproduction than predicted by models, i.e., J/ψ mesons tend towards lower values of z. Models of J/ψ production typically couple fixed-order perturbative quantum chromodynamics calculations with nonperturbative matrix elements that describe hadronization of the charm quark pair into a color-neutral state. A solution to this discrepancy was proposed in Ref. [9], where the evolution of the parton shower prior to the formation of the J/ψ is accounted for. By including this parton shower contribution, which is not adequately described in hadronization generators, such as PYTHIA [10], the authors were able to successfully describe the data. Although the LHCb measurement only concerns the subset of J/ψ mesons found inside a relatively high-p T jet, a recent measurement by CMS indicates that, for J/ψ mesons with energy larger than 15 GeV, nearly all of them are produced in association with a significant degree of jet activity [11].
Assuming this explanation of the LHCb data is correct, this paradigm shift in our understanding of J/ψ production has important implications for the interpretation of J/ψ data in nucleusnucleus collisions. It implies that J/ψ are not exclusively produced at short times, but may also be produced in the course of the interaction of a hard-scattered parton with the QGP. Hence, the suppression of the yield, as quantified by the nuclear modification factor (R AA ), may be sensitive to parton energy loss, the same phenomenon that gives rise to jet quenching. There are already some hints in this direction. First, as observed in Ref. [12], the J/ψ R AA in PbPb collisions [13,14] appears to exhibit the same rise with p T as light hadrons show [15], which for the case of light hadrons is well described by parton energy loss models [16][17][18][19][20][21]. Second, in mid-central collisions the J/ψ show a significant magnitude of elliptic anisotropy in their azimuthal angle (v 2 ) [22,23] at large p T , a region where the hydrodynamical effects that produce such an anisotropy are expected to die out. We know of no explanation for this high-p T v 2 feature other than path-length-dependent parton energy loss.
The goal of the current measurement is to investigate the z dependence of the nuclear modification factor of jets containing a J/ψ meson, i.e., This ratio compares the per-event yield in PbPb collisions (N AA ) to the expectation based on pp collisions, by scaling the cross section of the latter (σ pp ) by T AA . The factor T AA is the average effective nucleon-nucleon luminosity delivered by a single heavy ion collision for a given centrality selection (a quantity related to the collision impact parameter) [24]. In the absence of nuclear effects, R AA = 1.
This Letter constitutes the first direct study of the nuclear modification of J/ψ mesons inside jets. Jets with a constituent J/ψ meson are studied in the jet p T range of 30 < p T < 40 GeV for 0-90%, as well as 0-20% and 20-90% collision centralities. The jets are required to be in the pseudorapidity range |η| < 2, such that they are completely contained in the tracker acceptance, with no explicit selection on the rapidity of the J/ψ. The p T of the J/ψ is measured down to a threshold of 6.5 GeV. This gives a range of z from 0.22 to 1. We investigate to what extent the R AA varies with z, and indirectly, the formation time of the J/ψ. These data potentially constrain the roles of the different QGP interaction mechanisms that may be at play, in particular parton energy loss and Debye screening.

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Forward calorimeters extend the η coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid. Muons are measured in the range |η| < 2.4, with detection planes made using three technologies: drift tubes, cathode strip chambers, and resistive plate chambers. The efficiency to reconstruct and identify muons is greater than 96% over the full range of η. Matching muons to tracks measured in the silicon tracker results in a relative transverse momentum resolution, for muons with p T up to 100 GeV, of 1% in the barrel and 3% in the endcaps [25]. The forward hadron (HF) calorimeter uses steel as an absorber and quartz fibers as the sensitive material. The two halves of the HF are located 11.2 m from the interaction region, one on each end, and together they provide coverage in the range 3.0 < |η| < 5.2.
Events of interest are selected using a two-tiered trigger system [26]. The first level (L1), composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz for high luminosity pp collisions and 30 kHz for PbPb collisions. The second level, known as the high-level trigger (HLT), consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage, for both pp and PbPb collisions.
The particle-flow algorithm [27] aims to reconstruct and identify each individual particle in an event, with an optimized combination of information from the various elements of the CMS detector. The energy of photons is obtained from the ECAL measurement. The energy of electrons is determined from a combination of the electron momentum at the primary interaction vertex as determined by the tracker, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with originating from the electron track. The energy of muons is obtained from the curvature of the corresponding track. The energy of charged hadrons is determined from a combination of their momentum measured in the tracker and the matching ECAL and HCAL energy deposits, corrected for the response function of the calorimeters to hadronic showers. Finally, the energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energies.
Collision centrality is determined from the total transverse energy deposited in both of the HF calorimeters. The centrality is expressed as a percentage of the total inelastic hadronic cross section, with 0% representing the most head-on (central) collisions and 100% the most peripheral collisions. Hadronic events are selected by requiring at least three towers in each half of the HF calorimeter with an energy larger than 4 GeV. In this analysis, we restrict to the 90% most central events, where the hadronic event selection is fully efficient.
A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [28].

Simulation
Monte Carlo simulations are used as a baseline for the J/ψ efficiency and acceptance, as well as the detector response to jets. Correction factors to account for data-to-simulation discrepancies are discussed in the following section. Events are generated with PYTHIA 8 (version 8.230) [10], using the CP5 underlying event tune [29]. Prompt J/ψ are produced using all color-singlet and color-octet modes with the default matrix elements [30]. Simulation of PbPb collisions is done by embedding PYTHIA 8 events into PbPb collision events produced with the HYDJET generator (version 1.9) [31]. The event activity at mid-rapidity is matched to that of minimum bias collision data via an analysis of the energy deposited in randomly distributed cones. The response of the CMS detector is simulated using the GEANT4 package [32].

Prompt J/ψ meson signal extraction
J/ψ mesons are measured via their decays into oppositely charged muon pairs. Events are selected with a trigger that requires that at least two muon candidates are reconstructed in the muon subsystem, first at L1, and then using refined information at HLT. In PbPb collisions, the rate is reduced at HLT by requiring one of the two muons to match a track from the silicon tracking subsystem. The invariant mass of the pair (m µµ ) is required to be within the range of 1 to 5 GeV.
A set of offline muon selection criteria that is optimized for low p T muons is applied [13]. The J/ψ yield is obtained from a fit to the m µµ distribution in bins of jet p T and z. An example of a m µµ fit is shown in the left panel of Fig. 1. The signal is modeled as the sum of two Crystal Ball functions [33] with different widths, but common mean and tail parameters. The tail parameters are obtained from simulation. The dimuon background is modeled by a Chebyshev polynomial.
Nonprompt J/ψ mesons, i.e., those that are produced in b hadron decays, are separated from the prompt component by exploiting the long lifetime of these decays. A fit is performed to the distribution of the pseudo-proper decay length, l J/ψ = Lm/|p|, where L is the distance along the beam axis between the J/ψ vertex and the nearest primary collision vertex, and m and p are the world-average J/ψ mass (3.097 GeV) [34] and the J/ψ candidate momentum, respectively. An example of a l J/ψ fit is shown in the right panel of Fig. 1. In order to parameterize the l J/ψ distribution, the sPlot technique [35] is used to decompose it into the J/ψ signal and non-J/ψ background components, using m µµ as the discriminating variable. Prompt J/ψ mesons can have nonzero values of l J/ψ , both positive and negative, because of the finite detector resolution. The resolution is modeled as a sum of two Gaussian functions, which are constrained by fitting the negative part of the l J/ψ distribution. The nonprompt component is modeled as an exponential function, with the lifetime treated as a free parameter, convolved with the same double-Gaussian resolution function. The combinatorial background is itself well described by prompt-and nonprompt-like components, again using the same resolution function. The nonprompt component of the background is described by an empirical sum of exponential functions. To extract the yield of prompt J/ψ mesons a two-dimensional fit to the joint m µµ and l J/ψ distribution is performed. All parameters are fixed by the one-dimensional fits, aside from the prompt and nonprompt J/ψ signal normalizations, and the background normalizations, as well as the nonprompt signal lifetime parameter. A more detailed description of the fitting procedure can be found in Ref. [13].
The J/ψ meson acceptance, and the efficiency of reconstruction and identification, are determined in simulation in finely binned histograms of the p T and η of the J/ψ. For PbPb collisions, the efficiency is also determined in bins of collision centrality. The correction for acceptance and efficiency is applied as a weight factor to each J/ψ candidate prior to the signal extraction. Differences between data and simulation are evaluated in-situ, from events collected with a single-muon trigger, using the tag-and-probe technique [36]. These additional corrections are applied to the simulation in relatively coarse bins. Tag-and-probe corrections are derived separately for trigger, identification, and reconstruction efficiency.

Jet p T determination
Jets are clustered using the anti-k T algorithm [37,38] with a distance parameter of R = 0.3. The constituents of the jets are the output of the particle-flow algorithm, described in Section 2, with the exception of the J/ψ candidate. Muon pairs with an invariant mass consistent with a J/ψ decay are replaced by the reconstructed J/ψ candidate. Jet energy corrections are derived from simulation as a function of p T and η using the framework described in Ref. [39], in which datato-simulation scale factors are obtained from Z+jet, γ+jet and dijet p T balancing studies. In the case of PbPb collisions, these corrections are derived from peripheral events to avoid additional p T imbalance from jet quenching. The jet energy correction procedure is only applied to the non-J/ψ component of the jet, as the momentum of the J/ψ is determined with high precision from its dimuon decay. In PbPb collisions, the constituent subtraction method [40] is employed to subtract the contribution to the jet momentum from the underlying event. As the signal particle of interest, the J/ψ meson is excluded from the subtraction, as it would be incorrect to attribute any fraction of its energy to the underlying event.
In addition to the energy scale, the energy resolution of jets is somewhat degraded in data compared to simulation. The corresponding scale factors are derived from dijet balancing studies performed in pp collisions at √ s = 13 TeV in 2017 and 2018 [39]. The resolution is found to be between 10 and 20% worse than the simulation, depending on η. In order to incorporate this effect, a Gaussian smearing is applied to the measured jet p T values in simulation to match the resolution in data. This smeared mapping is then applied in the unfolding procedure, as described below.
The finite p T resolution of jets causes migration of jets across z bins. It also causes migration of jets into and out of the nominal p T selection of 30-40 GeV. In order to capture the full jet response, we also perform the measurement in underflow and overflow bins of 10-30 GeV and 40-50 GeV, respectively. With the current data, a reliable yield cannot be extracted for p T intervals larger than the selected overflow bin, which motivates the choice of the nominal p T interval. To correct for these bin migration effects, we unfold the data simultaneously in these two dimensions (jet p T and z) using the iterative D'Agostini method [41], as implemented in the ROOUNFOLD package [42]. The detector response matrix, which defines the relationship between the true and measured values, is taken from the PYTHIA 8 simulation, embedded into HYDJET for the case of PbPb collisions, aside from the aforementioned data-to-simulation scale factors for the jet energy scale and resolution. The response matrices for pp and PbPb collisions are shown in Fig. 2. The coarse bins show the jet p T dependence, whereas the finer inset bins show the dependence on z. Each column of measured bins is normalized to unity, such that each bin represents the fractional contribution of the given true p T and z values to the measured values for that column. The PbPb data exhibit substantially larger off-diagonal contributions than for pp, as expected from the larger underlying event, which drives the bin migration in PbPb.
In this method, the unfolding is initialized with a "prior" distribution in the two variables, which is taken from simulation. To avoid bias from the presumed shape of the z distribution in simulation, which is known to be inaccurate [8], the prior distribution is flattened in z. After a tunable number of iterations, which corresponds to the degree of regularization, the procedure is truncated. To improve the performance of the unfolding, we run the full set of iterations multiple times. The prior of each such "super-iteration" is obtained using the z distribution that is output from the previous super-iteration. The numbers of iterations and super-iterations are optimized using simulated prompt J/ψ events, by applying a random Gaussian smearing to the measured z distributions according to the relative statistical uncertainties in data, in order to emulate the effect of statistical fluctuations. The optimization criterion is the minimum χ 2 between the unfolded and true z distributions. The use of three iterations is found to give the best performance in both pp and PbPb simulated samples. However, while only one superiteration is sufficient for pp, 25 super-iterations are required for PbPb collisions. Since the optimal settings depend not only on the statistical uncertainties of the measured z distribution, but also on its shape, we additionally evaluated them in nonprompt J/ψ meson simulation, where the measured distribution is more similar to the data than in the prompt J/ψ meson simulation, by unfolding the nonprompt z distributions with the nonprompt response matrices and using the settings with the best performance as a systematic variation in the regularization procedure.

Systematic uncertainties
The systematic uncertainties may be divided into three categories: J/ψ signal extraction, jet energy scale and resolution, and normalization uncertainties.
J/ψ signal extraction: Uncertainties in the extraction of the J/ψ yields arise from the signal and background shape modeling, as well as from the acceptance and efficiency of muon reconstruction and identification. The signal shape of the m µµ distribution is varied from the two Crystal Ball parameterization to a single Crystal Ball function convoluted with a Gaussian function. The radiative tail of the signal model is also varied by treating the tail parameters of the Crystal Ball function as free parameters, rather than taking their values from simulation. The relative uncertainty in the J/ψ yield coming from the signal modeling is less than 1% in pp and less than 2% in PbPb collisions. The dependence of the signal shape on l J/ψ is evaluated by switching from a parameterization of the lifetime distribution of nonprompt J/ψ signal to a template built from simulation. The relative uncertainty from this variation is less than 1.5% in both pp and PbPb collisions, except at low J/ψ p T (<10 GeV) in pp collisions, where it reaches 4.5%, and where it dominates the uncertainties in the J/ψ signal extraction. The Gaussian function describing the lifetime resolution is varied to use the value obtained from simulation, rather than obtaining it from the data. This variation affects both the signal and background modeling, and results in a relative uncertainty of less than 2% in pp and PbPb collisions.
The uncertainty in the description of the background is obtained as follows. Instead of a Chebyshev polynomial, an exponential of a Chebyshev polynomial is used to describe the m µµ distribution. The range of m µµ used in the fit is varied as well, to reduce the contamination from the ψ(2S) state. This results in an uncertainty of less than 1% in pp and PbPb collisions, except at low J/ψ p T (< 10 GeV) in PbPb collisions, where it reaches 6% and dominates the uncertainties in the J/ψ signal extraction. As was done for the signal, a template is used to describe the l J/ψ distribution, rather than a parameterization. In this case, the template is obtained from a background-like event sample using sPlot. The source results in a relative uncertainty of < 2% in the J/ψ yields in both pp and PbPb collisions.
The acceptance and efficiency corrections lead to an associated relative uncertainty in the J/ψ yields in pp collisions of 1.5-2.5%, and 3-4% in PbPb collisions, depending on the value of J/ψ p T , and comprise the dominant source of uncertainty in the J/ψ signal extraction for J/ψ p T > 10 GeV. The largest component of the uncertainty is statistical in nature, coming from the finite size of the single-muon trigger samples used in the tag-and-probe method. The statistical uncertainty of the simulation samples is also taken into account, but is generally subdominant. A systematic component is evaluated by variations of the signal and background modeling in the extraction of the J/ψ yield from fits to the m µµ distribution, similar to the procedure for the J/ψ yield extraction in the main analysis, described above. The relative uncertainty from this source is around 1% in both pp and PbPb collisions, with no strong dependence on p T . Further details on the procedure used to evaluate systematic uncertainties for muons from the tag-and-probe method are reported in Refs. [25,36].
Jet energy scale and resolution: The uncertainty in the jet energy scale is evaluated from dijet and γ+jet balancing methods, as described in Ref. [39]. The uncertainty in the jet energy scale in pp collisions is around 3-4%, increasing as a function of |η|. In PbPb collisions, the uncertainty is around 4%, except for the barrel-endcap transition region (1.3 < |η| < 1.6), where it can become as large as 10%. An additional source of uncertainty in the jet energy scale is considered in PbPb collisions to take into account the differences in the particle mixture in simulation and data due to jet quenching, as described in Ref. [43]. This uncertainty is around 5% for the 30% most central events and 2.5% for all other centralities. The uncertainty in the jet energy resolution is evaluated from a dijet balancing method, as also described in Ref. [39]. In pp collisions, the uncertainty in the resolution is in the range of 2-4% in the barrel region, but is larger in the endcap and transition regions, where it varies in the range of 10-20%, depending on η. For PbPb collisions, the uncertainty is evaluated from peripheral data, as well as from pp data recorded in the same yearly running period, and varies from 3% in the barrel to 7% in the endcap and transition regions. In PbPb, there is an additional contribution to the uncertainty because of the modeling of the underlying event in simulation with HYDJET. This uncertainty is evaluated by comparing the energy in randomly distributed cones in data and simulation. The difference between the random cone distributions in data and simulation is used to estimate the effect on the jet resolution.
Unfolding: The uncertainty in the regularization applied is estimated by varying the number of iterations used in the unfolding, from the settings found to be optimal for prompt J/ψ to those found to be optimal for the nonprompt J/ψ signal, which has a very different signal shape as a function of z. For nonprompt simulations in pp and PbPb collisions, the use of ten iterations is found to give the best performance, in contrast to three iterations for the prompt simulation. In both cases, the statistical precision of the data is emulated in simulation. The assumption of a prior that is flat in z is also relaxed, by instead initializing the prior to match the shape of the truth z distributions in simulated nonprompt J/ψ signal samples, which are more similar in shape to the measured z distributions. Finally, the statistical uncertainty of the simulated samples used in the unfolding is taken into account, and is generally subdominant. Systematic uncertainties related to the jet energy scale and resolution are propagated by producing variations of the detector response matrix and repeating the unfolding procedure for each variation.
Normalization uncertainties: Cross section measurements in pp collisions have an overall uncertainty from the integrated luminosity of 1.9% [44] that is obtained from an analysis of data from a van der Meer scan [45]. The PbPb data are normalized by the equivalent number of hadronic nucleon-nucleon interactions in the data sample, which has an uncertainty of 1.3% coming from the selection of such events, taking into account possible contamination from electromagnetic interactions and beam backgrounds. To compare to data from pp collisions, the PbPb data are normalized by T AA , which is determined from the Monte Carlo implementation of the Glauber model described in Ref. [46]. For the centrality intervals used in this analysis, the corresponding values of T AA are 6.27 ± 0.14, 18.79 ± 0.36, and 2.717 ± 0.098 mb −1 , for the 0-90, 0-20 and 20-90% centrality intervals, respectively.
The various sources of systematic uncertainty are shown in Fig. 3. Sources related to the J/ψ yield extraction are combined into a single component, for visibility. The total uncertainty is the quadrature sum of all the sources. Some sources of systematic uncertainty are largely correlated bin-to-bin, notably the uncertainty in the jet energy scale, jet energy resolution, and prior. No cancellation of systematic uncertainties is assumed, however, when comparing pp to PbPb collisions.   Figure 4 shows the distribution in pp data of the fragmentation variable z for prompt J/ψ mesons. Its shape is compared to generator-level predictions from PYTHIA 8 for prompt and nonprompt J/ψ signals. In contrast to the PYTHIA 8 simulation, where prompt J/ψ are produced directly in the matrix element partonic scattering, the data show a relatively large degree of surrounding jet activity, indicative of J/ψ production inside of parton showers. The z distribution in data more closely resembles that of the nonprompt J/ψ PYTHIA 8 simulation, which contains a larger jet-like component from fragmentation, as well as other products of the b-hadron decay. The data confirm the trends observed in Ref.

Results
[8], but in a different rapidity range and nucleon-nucleon center-of-mass energy.   The PbPb data show a suppression level that is generally comparable to that observed for "inclusive" prompt J/ψ production, i.e., without an explicit jet requirement [13]. This is quantified by the ratio of these two distributions, R AA , shown in Fig. 5 (right). The data show a slight rising trend as a function of z, with a significance of around two standard deviations. Fig. 6 shows the R AA for two centrality selections, 0-20 and 20-90%. A larger degree of suppression is suggested for the more central selection, as expected for final-state effects related to the QGP. The rising trend with increasing z is somewhat more pronounced in central events.
In the largest z bin, where the J/ψ is produced with fewer associated particles, the suppression is significantly reduced as compared to lower values of z, most importantly in the centralityintegrated results. Such a reduction of suppression at large z has a natural interpretation in terms of the jet quenching phenomenon. Lower values of z should be populated with jets with a J/ψ produced late in the parton shower. Such a parton cascade is expected to have a large degree of interaction with the QGP in the form of subsequent medium-induced emissions, as compared to a jet with a small partonic multiplicity [47]. In this picture, the rising trend observed for inclusive prompt J/ψ production would be explained by the same mechanism, as z tends to increase with increasing p T .

Summary
Jets containing a prompt J/ψ meson were studied in proton-proton (pp) and lead-lead (PbPb) collisions at √ s NN = 5.02 TeV, for jets with transverse momentum 30 < p T < 40 GeV and pseudorapidity |η| < 2. The distribution of the fragmentation variable z, the ratio of the J/ψ p T to that of the jet, was measured in both systems. In pp collisions, prompt J/ψ mesons were found to have more surrounding jet activity, i.e., to populate lower values of z than predicted by PYTHIA 8 simulations, suggesting that J/ψ production late in the parton shower is underestimated. The pp and PbPb distributions were compared by calculating the nuclear modification factor, R AA , the ratio of PbPb data to the expectation based on pp data. The value of R AA as a function of z shows a rising trend. The suppression at low z is found to be larger in the 20% most central events (i.e. "head-on" collisions), as compared to the less central selection. The results show explicitly that the J/ψ produced with a large degree of surrounding jet activity are more highly suppressed than those produced in association with fewer particles. This finding emphasizes the importance of incorporating the jet quenching mechanism in models of J/ψ suppression.   [6] LHCb Collaboration, "Measurement of J/ψ polarization in pp collisions at √ s = 7 TeV", Eur. Phys. J. C 73 (2013) 2631, doi:10.1140/epjc/s10052-013-2631-3, arXiv:1307.6379.