Analytic and Monte Carlo studies of jets with heavy mesons and quarkonia

We study jets with identified hadrons in which a family of jet-shape variables called angularities are measured, extending the concept of fragmenting jet functions (FJFs) to these observables. FJFs determine the fraction of energy, z, carried by an identified hadron in a jet with angularity, τa. The FJFs are convolutions of fragmentation functions (FFs), evolved to the jet energy scale, with perturbatively calculable matching coefficients. Renormalization group equations are used to provide resummed calculations with next-to-leading logarithm prime (NLL’) accuracy. We apply this formalism to two-jet events in e+e− collisions with B mesons in the jets, and three-jet events in which a J/ψ is produced in the gluon jet. In the case of B mesons, we use a phenomenological FF extracted from e+e− collisions at the Z0 pole evaluated at the scale μ = mb. For events with J/ψ, the FF can be evaluated in terms of Non-Relativistic QCD (NRQCD) matrix elements at the scale μ = 2mc. The z and τa distributions from our NLL’ calculations are compared with predictions from monte carlo event generators. While we find consistency between the predictions for B mesons and the J/ψ distributions in τa, we find the z distributions for J/ψ differ significantly. We describe an attempt to merge PYTHIA showers with NRQCD FFs that gives good agreement with NLL’ calculations of the z distributions.


Introduction
The study of jets and heavy flavor continues to play an important role at the Large Hadron Collider (LHC) and many other high energy and nuclear experiments. Such studies are essential for testing our understanding of Quantum Chromodynamics (QCD) and for calculating backgrounds in searches for new physics. In this paper we calculate cross sections for e + e − to jets, where one of the jets contains a hadron with either open or hidden heavy flavor. In particular, we will derive factorization theorems and perform analytical Nextto-Leading-Log prime (NLL') resummation 1 for these cross sections using renormalization group (RG) techniques. We will also compare our results with monte carlo simulations of the same cross sections.
Recently, there has been considerable interest in cross sections of this type [2][3][4][5][6][7][8][9][10][11]. Ref. [2] demonstrated that the cross section for producing a jet with an identified hadron can be determined using a distribution function called the fragmenting jet function (FJF). FJFs are in turn related to the more commonly studied fragmentation functions (FFs) by a matching calculation at the jet energy scale. This implies that cross sections for jets with an identified hadron provide a new arena to measure FFs, which are more commonly extracted from the semi-inclusive cross section e + e − → H + X. Especially important is that this provides an opportunity to extract gluon FFs [10,11], since quark FFs are more j 1 0 dzzJ ij (s, z, µ) . (1. 3) The properties of FJFs were further studied in refs. [3][4][5][6][7]. These papers focused on the FJFs for light hadrons such as pions. FJFs for particles with a single heavy quark were studied in ref. [8] and FJFs for quarkonia were calculated in ref. [9]. One important goal of this work is to generalize FJFs to jets in which the angularity is measured. The angularity, denoted τ a , is defined as [13] where the sum is over all the particles in the jet, and ω = i p − i is the large light-like momentum of the jet. The angularity should be viewed as a generalization of the invariant mass squared of the jet since s = ω 2 τ 0 . We calculate the matching coefficients appropriate for jets in which the angularity has been measured, denoted J ij (τ a , z, µ), and verify the s → τ a generalization of the sum rules in eq. (1.3) in appendix B of this paper. The other goal of this work is to study the z and τ a dependence of the cross section for jets with identified heavy hadrons in e + e − collisions and compare our analytical results to monte carlo simulations. We will do this for two-jet events in which e + e − → bb is followed by JHEP06(2016)121 fragmentation to B mesons. We will also study three-jet events with e + e − → bbg followed by the gluon fragmenting to a jet with a J/ψ. At the LHC we expect high energy gluons fragmenting to a jet with J/ψ to be an important production mechanism of J/ψ at high p T and ref. [9] showed this process is sensitive to the mechanisms underlying quarkonium production. The study in this paper will allow comparison of analytic calculations with monte carlo simulations of gluons fragmenting to J/ψ in jets. In order for this cross section to be physically observable one would either include quarks and antiquarks fragmenting to jets with J/ψ or one would have to ensure experimentally that the J/ψ came from the gluon jet in the three-jet event, which could be possible if the other jets are b-tagged.
In section 2, we discuss the basics of FJFs for events containing jets where the angularity of the one of the jets is probed. We review various properties of FJFs and their relationship with the more commonly studied FFs. We also present our results for the matching coefficients J ij (τ a , z, µ) for jets with measured angularities. Further details of that calculation can be found in appendix B. In section 3, we present our results for the NLL' cross section for e + e − → 2 jets where one of the jets contains a B meson and the angularity of that jet is measured. We find reasonable agreement in both z and τ a distributions between our analytic calculations and monte carlo simulations performed using Madgraph [14] + PYTHIA [15,16] and Madgraph + HERWIG [17]. In section 4, we show similar comparisons of analytic versus monte carlo calculations for the cross section for e + e − → 3 jets where one of the jets contains a J/ψ created via gluon fragmentation. In this case the τ a distributions for the jet are in good agreement, but the monte carlo predictions for the z distributions are inconsistent. We believe that this is due to PYTHIA's modeling of radiation from color-octet states that produces a harder z distribution than the analytic calculations. In an effort to improve the consistency between NLL' and monte carlo calculations, we turn off hadronization in PYTHIA and then convolve the distribution of momenta of the gluons within a jet with the NRQCD color-octet FF at the scale 2m c . This ad-hoc procedure brings monte carlo calculations into much better agreement with analytic NLL' calculations. This suggests that if NRQCD fragmentation could be properly implemented in PYTHIA, consistency with NLL' calculations would be obtained, though more work needs to be done on this problem. In section 5 we give our conclusions. Appendix A summarizes the renormalization group evolution (RGE) needed for NLL' calculations and also gives the profile functions that are used when computing the scale variation in the NLL' calculations. Appendix B describes the calculation of the matching coefficients and checks that they satisfy the required sum rules that relate them to the jet function.

Fragmenting jet functions with angularities
In this section we extend the calculation of ref. [4] to FJFs with measured angularities. We will follow the terminology of ref. [18], in which a jet whose angularity is measured is referred to as a "measured" jet, while a jet for whom only the total energy is measured but the angularity is not is called an "unmeasured" jet. Here we consider the case of two particles as this is the most that will appear in a one-loop calculation. In ref. [4] the JHEP06(2016)121 measurement operator in the definition of FJFs forces the mass squared of the jet to be s. The measurement operator takes the form δ(ω(k + − l + − p + )) = δ(s − ω(l + + p + )), (2.1) where k µ is the parent parton's momentum and l µ and p µ are the momenta of the partons carrying large lightcone components l − = (1 − z)k − and p − = zk − of the parent's momentum, respectively. The operator definition of the FJF with measured angularities is given by where at O(α s ) the operatorτ a takes the form (cf. eq. (1.4)) Other than replacing eq. (2.1) with eq. (2.3), the integrals of all diagrams are the same as in ref. [4]. However, rather than using the δ-regulator and a gluon mass, we will use pure dimensional regularization to regulate all divergences. In this limit, it is possible to show that the one-loop evaluation of the FF yields where T ij are the color structures, T qq = C F , T gg = C A , T qg = C F , T gq = T R . Additionally, we have verified that the same 1/ IR poles appear in the calculation of FJFs and appropriately cancel in the matching between the FJFs and FFs for all values of a < 1. This justifies the formula which is the analog of eq. (1.1) for FJFs that depend on the angularities. Since the matching coefficients J ij (τ a , z, µ) are free of IR divergences, we can simplify the matching calculation by using pure dimensional regularization, setting all scaleless integrals to zero and interpreting all 1/ poles as UV. A detailed calculation of the renormalized finite terms of J ij (τ a , z, µ) can be found in appendix B, the results of which are shown below. We parametrize the matching coefficients J ij (τ a , z, µ) as and where theP ij are the splitting functions of ref. [4] except for the case i = j = q, (2.9) Notice that our results for the matching coefficients J ij (τ a , z, µ) are independent of the jet algorithm and the jet size parameter R. To include modifications of the J ij (τ a , z, µ) that come from these effects, one would have to multiply the measurement operator in eq. (2.3) by an additional Θ-function that imposes the phase space constraints required by the jet algorithm. However, for jets with measured angularities, it was shown in ref. [18] that jet-algorithm dependent terms for cone and k T -type algorithms are suppressed by powers of τ a /R 2 . Inuitively, this is because as τ a → 0 all the particles in the jet lie along the jet axis so the result must be insensitive to which algorithm is used and to the value of R in this limit. For the values of τ a and R considered in this paper, τ a /R 2 is negligible and we will drop these corrections. As a non-trivial check of our results we show in appendix B that our J ij (τ a , z, µ) satisfy the following identities and sum rules, lim a→0 J ij (τ a , z, µ) = ω 2 J ij (s, z, µ) , (2.10)
3 e + e − → 2 jets with a B meson In this section we present an analytic calculation of the cross section for e + e − to two b jets in which the B meson is identified in a measured jet. Following the analysis of ref. [18], the factorization theorem for the cross section for one measured b jet and one unmeasuredb jet is where H 2 (µ) is the hard function, S unmeas (µ) and S meas (τ a , µ) are the unmeasured and measured soft functions, J n (τ a , µ) is the measured jet function containing the b quark. These describe the short-distance process, surrounding soft radiation, and radiation collinear to unmeasured and measured jets, respectively. At NLO the τ a -independent functions are given by where Λ is a veto on out-of-jet energy, r = tan (R/2) and J q alg (µ) is a function that depends on the algorithm used (and we will use the cone algorithm below) and is given in eq. (A.18) of ref. [18]. We note that unlike measured jets, algorithm dependent contributions to the unmeasured jet are not power suppressed. We also note that, beginning at O(α 2 ), nonglobal logarithms of the ratio Qτ a /(2Λr 2 ) begin to appear in the cross-section [19]. For the values of the parameters we consider, these ratios are such that we can treat these logarithms as O(1) and thus these would enter as fixed order corrections needed at NNLL' accuracy, which is beyond the scope of this work.
We suppress the dependence of all these functions on scales other than the renormalization scale µ. Measured functions are convolved according to To calculate the differential cross section for a measured jet with an identified B hadron, we apply the analogous replacement rule in eq. (1.1) to eq. (3.1) and use the expression for

JHEP06(2016)121
the FJF in eq. (2.5) to obtain To obtain an NLL' resummed formula for the cross section, we evaluate each function in the factorization theorem in eq. (3.4) at its "characteristic" scale (where potentially large logarithms are minimized) and, using renormalization group techniques, evolve each function to a common scale, µ, which we will choose to be equal to the hard scale. The details of this evolution are discussed in appendix A. The convolutions in eq. (3.4) must be performed over angularity over S meas , J ij , and factors arising from RG equations. Since such RG factors are distributions (δ or plusdistributions) in the angularity our final answer is written in terms of distributions that can be computed analytically using eqs. (A.18)-(A.20). Upon performing convolutions and resummation to NLL' accuracy we find for the cross section where the '+' distribution is defined in eq. (A.15) (and acts on all τ a -dependent quantities, including any implicit dependencies arising from the choice of scales µ F ) and Ω(µ Jn , µ S meas ) = ω Jn (µ, µ Jn ) + ω S meas (µ , µ S meas ), the functions ω Jn and ω S meas are given in appendix A, the function f S is given by [18] and f ij J are written in terms of the coefficients c ij 0 , c ij 1 and c 2 presented in eq. (2.7) as (3.8) Table 1. Characteristic scales of the different functions in the factorization theorem of eq. (3.1).

JHEP06(2016)121
The evolution kernel Π is given in terms of K F (µ, µ 0 ) and ω F (µ, µ 0 ) (cf. appendix A), where µ F , m F and j F are given in table 1. Because they involve FFs (cf. appendix B), the z convolutions must be evaluated numerically. For the fragmentation of the b quark we use a two-parameter power model FF introduced in ref. [20], Values for the parameters α = 16.87 and β = 2.628 with χ 2 d.o.f. = 1.495 were determined using a fit to LEP data in ref. [21] for the inclusive process e + e − → B + X. Errors in these parameters were not quoted in ref. [21], so we cannot quantify errors associated with the extracted FF in our calculation. Additionally, we neglect the contribution from the fragmentation of other partons for our e + e − collider studies as in ref. [21]. In proton-proton collisions at the LHC, gluon FJFs must also be included since the dijet channel gg → gg gives a significant contribution to the production of jets with heavy flavor [11]. For the evolution of the FF up to the jet scale we solve the DGLAP equation using an inverse Mellin transformation as done in ref. [9]. Figure 1 shows the z distributions from dσ(τ 0 , z) for τ 0 = (1.5, 2.0, 2.5) × 10 −3 of our analytic NLL' calculation (green) and monte-carlo simulations using Madgraph + PYTHIA (black) and Madgraph + HERWIG (red). For each monte carlo and for each NLL' calculation, the graphs are independently normalized to unit area. For plots with fixed τ a we use a z-bin of ± 0.1 and for plots with fixed z we use a τ a bin of size ± 2 × 10 −4 . Jets are reconstructed in PYTHIA using the Seedless-Infared-Safe Cone (SISCONE) algorithm in the FastJets package [22] with R = 0.6, which will be used throughout this work. We produced simulated dijet events at E cm = 250 GeV in which each jet has an energy of at least (E cm − Λ)/2 where Λ = 30 GeV. 2 The central green line corresponds to the NLL' calculation with the various functions in the factorization theorem evaluated at their characteristic values shown in table 1, and the green band corresponds to the estimate of theoretical uncertainty obtained by varying the scales of the unmeasured functions by ±50%, and using profile functions [23][24][25]  breaks down and at a fixed scale for small values of τ a where we reach the non-perturbative regime. This method for estimating theoretical uncertainties is used throughout this work. Additional details on the profile functions we use can be found in appendix A.
The blue curves in figure 2 show the differential cross section as a function of z for fixed τ 0 where µ J (τ ) = µ J (τ, z = 0) = ωτ 1/(2−a) is chosen as the characteristic scale of the measured jet function, and the error band is obtained the same way as for figure 1. As in figure 1, the green curves show the cross section for a measured jet scale µ J (τ, . The reorganization of logarithms of (1 − z) shown in eq. (A.24) suggests that we can improve the accuracy of our calculations for z → 1 by choosing the characteristic value of the measured jet scale to be µ J (τ, z). This improvement is clearly seen in figure 2 which shows the scale variation for the choices µ J (τ ) and µ J (τ, z), the latter choice gives smaller scale variation near the peak in the z distribution.
In figure 3 we present the results for the τ 0 distributions of the differential cross section dσ(τ a , z) for z = 0.4, 0.6, and 0.8. The color and normalization schemes match those in figure 1. We see that for higher values of z the distributions of τ 0 are shifted towards smaller values. This is expected since the majority of the energy of the jet is carried by the B meson which results in narrower jets. Figures 1 and 3 show that our results are consistent within the monte carlo uncertainty that is suggested by the difference between PYTHIA and HERWIG predictions. This gives us confidence that the FJF formalism combined with NLL' resummation can be used to correctly calculate both the substructure and the identified hadron's energy fraction within a jet.
4 e + e − → 3 jets with the gluon jet fragmenting to J/ψ We can also use the FJF formalism to calculate the cross section for e + e − → 3 jets with a J/ψ. As we expect gluon fragmentation to be the dominant production channel at the LHC, we focus on the case where J/ψ is found within a gluon jet. In addition, we assume that the angularity of this jet is also measured. To obtain a physical observable, one must also include contributions from all jets fragmenting to J/ψ, however, we expect the contribution from quark jets to be smaller. It is theoretically possible to isolate the J/ψ coming from gluon jets in experiments by b-tagging the other two jets in the event, so we will focus on the process e + e − → bbg followed by gluon fragmentation to J/ψ. The analytic expression for this cross section is where Ω ≡ Ω(µ Jn 3 , µ S meas ) = ω Jn (µ, µ Jn 3 ) + ω S meas (µ , µ S meas ), the b-quark initiated jets J n 2 are unmeasured, the expression for f S is the same as eq. (3.7) with C F replaced by C A , and our expressions for f ij J are given in terms of the coefficients c ij 0 , c ij 1 and c 2 given in eq. (2.7). Here σ 0 is the LO cross section for e + e − → bbg. We will focus on the Mercedes Benz configuration in which all three jets have (approximately) the same energy, and consider jets with energies large enough that the mass of b-quark can be neglected. Here, H 3 (µ) is 1 + O(α s ) where the O(α s ) comes from the NLO virtual corrections to e + e − → bbg. We do not include this correction. The primary effect of its omission will be on the normalization of the cross section, which is not important for our discussion of the distributions we show below, and to increase the scale uncertainty associated with varying µ H ; however this is not a very important source of uncertainty in our calculations.
While the calculation for B mesons requires a phenomenological FF, the FFs for J/ψ production can be calculated in NRQCD [12]. Refs. [26][27][28][29] showed that a J/ψ FF can be calculated in terms of analytically calculable functions of α s (2m c ) and z multiplied by nonperturbative NRQCD long-distance matrix-elements (LDMEs). In J/ψ production, the most important production mechanisms are the color-singlet mechanism, in which the cc is produced perturbatively in a 3 S 1 state, and the color-octet mechanisms, in which the cc is produced perturbatively in a 1 S   corresponding LDMEs are taken to be the central values from the global fits performed in refs. [30,31], and are shown in table 2. The color-singlet LDME scales as v 3 , where v is the typical relative velocity of the cc in the J/ψ, while the color-octet LDMEs scale as v 7 [12].
This Before discussing the comparison of our results with monte carlo, we briefly review how the Madgraph + PYTHIA monte carlo handles color-singlet and color-octet quarkonium production. We produce quarkonia states in Madgraph from the following processes: e + e − → bbggcc[ 3 S are for the cc produced in the event. We only include diagrams in which the virtual photon couples to the bb so in all cases the cc plus any additional gluons come from the decay of a virtual gluon. We did not simulate production in the 3 P (8) J channel in e + e − → bbg → bbccg because IR divergences in the matrix elements require much longer running times to get the same number of events. We then perform showering and hadronization on these hard processes using PYTHIA. Analysis is done using RIVET [32]. During PYTHIA's showering phase, color-singlet J/ψ do not radiate gluons. Thus if these J/ψ are produced within a jet, all surrounding radiation is due to the other colored particles in the event [15,16]. We require that after showering there are only three jets in the event, two from the b-quarks and one from a gluon that contains the J/ψ. We simulate three-jet events at E cm = 250 GeV in the Mercedes-Benz configuration by requiring the jets each have energies E jet > (E cm − Λ)/3 with Λ = 30 GeV, analagous to what was done in section 3.
For cc produced in a color-octet state PYTHIA allows the color-octet cc to emit gluons with a splitting function 2P qq (z). Since P qq (z) is peaked at z = 1, the color-octet cc pair typically retains most of its energy after these emissions. This model of the production mechanism is very different than the physical process implied by the NLL' calculation. In

JHEP06(2016)121
the NLL' calculation, the FF is calculated at the scale 2m c , then evolved up to the jet energy scale using Altarelli-Parisi evolution equations. Since this is a gluon FF, the most important splitting kernel in this evolution is P gg (z). We find that the FFs obtained at the jet energy scale are not significantly changed if we use only this evolution kernel and ignore mixing with quarks. Thus the production process implied by the NLL' calculation is that of a highly energetic gluon produced in the hard process with virtuallity of order the jet energy scale, which then showers by emitting gluons until one of the gluons with virtuality of order 2m c hadronizes into the J/ψ. Because P gg (z) is peaked at z = 0 and z = 1 the resulting J/ψ distribution in z is much softer than the model employed by PYTHIA. PYTHIA does not allow one to change the actual splitting function, only to modify the color-factor. Therefore, in order to get a softer z distribution we changed the coefficient of PYTHIA's splitting kernel for a gluon radiating off a color-octet cc pair from 2P qq to C A P qq = 3P qq . This results in a slighter softer z distribution than default PYTHIA, but is still inconsistent with the NLL' calculation. This change does not have significant impact on the τ a distributions. The τ a distributions are generally in better agreement. The variable τ a depends on all of the hadrons in the jet and is therefore less sensitive to the behavior of the J/ψ, especially when the J/ψ carries a small fraction of the jet energy. In that case, τ a distributions in the NLL' calculation look similar for all color-octet mechanisms.
In an attempt to see if PYTHIA can be modified to reproduce the z distributions obtained in our NLL' calculations, and confirm the physical picture of the NLL' calculation described above, we generate e + e − → bbg events in Madgraph and allow PYTHIA to shower but not hadronize the events. If we allow the shower to evolve to a scale where the typical invariant mass of a gluon is 2m c and then convolve the gluon distribution with the NRQCD FFs at this scale, we expect that the resulting z distributions should mimic our NLL' calculation. The lower cutoff scale in PYTHIA's parton shower is set by the parameter TimeShower:pTmin, which is related to the minimal virtuality of the particles in the shower, and whose default value is 0.4 GeV. We change this parameter to 1.6 GeV, which corresponds to a virtuality of ∼ 2m c , then obtain a z distribution for the gluons by randomly choosing a gluon from the gluon initiated jet. We then numerically convolve this z distribution with the analytic expression for the NRQCD FF. This procedure, which we will refer to as Gluon Fragmentation Improved PYTHIA (GFIP), yields z distributions that are consistent with our NLL' result, as we will see below. We tested an analogous procedure for two-jet events with B mesons by showering e + e − → bb with PYTHIA with hadronization turned off. We then convolved the resulting b quark distribution with the b-quark FF at the scale 2m b , and found results for B mesons that are consistent with our NLL' calculations. Note that PYTHIA treats the radiation coming from the octet cc pair the same regardless of the angular momentum quantum numbers. In contrast, GFIP like the NLL' calculation gives different results for all three channels by applying different FFs at the end of the parton shower phase. Also GFIP can be applied to all four NRQCD production mechanisms, since convergence issues for the 3 P     for smaller values of z and notice some qualitative differences in the tail regions, especially for the 1 S (8) 0 channel. At higher values of z where the number of final state particles is small, differences in the τ 0 distributions could be attributed to the increasing influence of Pythia's unrealistic model of quarkonium production. As z → 0, we also see similar τ 0 dependence for the two color-octet channels in our analytic results. This suggests that in the small z region, the jet substructure is independent of the production mechanism. Thus, attempts to use angularity distributions to extract the various LDMEs should focus on the range 0.3 < z < 0.7.
In figure 5, we show the angularity distributions (without uncertainties) for the 1 S mechanisms for a = +1/2, 0, −1/2, −1. These are computed analytically and using monte carlo and we again see reasonable agreement. As a is decreased, we see less discrimination between the two production mechanisms. Thus extraction of LDMEs should ideally be done with larger values of a, for a < 1 where factorization in SCET I holds, with the caveat that there is a trade-off since the predictability of the analytical results is limited for a too close to 1 since power corrections grow as 1/(1 − a) [33].
In contrast to the angularity distributions, figure 6 shows that analytic and monte carlo calculations of the z distributions using Madgraph + PYTHIA yield strikingly different results, with Madgraph+PYTHIA yielding a much harder z-distribution.  is far from a proper modification of PYTHIA, it shows us that implementing the missing g → J/ψ fragmentation yields encouraging similarities to our analytical calculations using the FJF formalism with NRQCD FFs. This also suggests that if monte carlo is modified to properly include NRQCD FFs at the scale 2m c it will yield results that are consistent with FJFs combined with NLL' resummation. Correct monte carlo implementation of the NRQCD FFs is important because the GFIP modification can only be used to calculate the z distribution. There are many other jet shape observables, such as N -subjettiness or ∆R (where ∆R is the angle between the J/ψ and the jet axis), that should be able to discriminate between NRQCD production mechanisms, and many of these are most easily predicted using monte carlo.

Conclusion
The study of hadrons within jets provides new tests of perturbative QCD dynamics. The distribution in z (the fraction of jet energy carried by the identified hadron) can be calculated as a convolution of the well-known fragmentation functions (FFs) for that hadron with perturbative matching coefficients that are calculable at the jet energy scale, which is typically well above Λ QCD . At hadron colliders this provides a new way to extract FFs and will be especially important for pinning down gluon FFs, which are of subleading importance JHEP06(2016)121  in e + e − colliders where FFs are usually measured. The production of heavy quarkonia within high energy jets in collider experiments also provides new tests of NRQCD.
In this paper, we studied cross sections for jets with heavy mesons as a function of z and the substructure variable angularity, τ a . We provided for the first time the NLO matching coefficients for jets with measured τ a , and used these along with the known RGE for the hard, jet, and soft functions to obtain NLL' accuracy calculations of cross sections for jets with heavy mesons. We considered the production of B mesons in twojet events in e + e − collisions at E cm = 250 GeV as well as J/ψ production in three-jet events at the same energies. Though not relevant to any experiment, this is useful for comparing NLL' calculations with monte carlo simulations of fragmenting jets whose energy is comparable to those measured at the LHC. In the simulations of quarkonia production, the underlying hard process was generated using Madgraph and then PYTHIA was used to shower and hadronize the events. In the simulations involving B meson production we also used HERWIG.

JHEP06(2016)121
For B mesons, we find that the z and τ a distributions computed using monte carlo and NLL' are in excellent agreement, giving us confidence in our analytic approach. In the case of J/ψ, we considered three-jet events in which the jets all had the same energy and the J/ψ in both simulation and NLL' calculations was required to come from the gluon jet. This allowed us to study J/ψ production via the fragmentation of high energy gluon initiated jets, which we expect to be an important mechanism at the LHC. Earlier studies of gluon FJFs in ref. [9] indicated that the z and E dependence of these jets could discriminate between various NRQCD production mechanisms. The analytic NLL' studies of this paper are consistent with ref. [9]; we also find that the τ a and z distributions can discriminate between different various NRQCD production mechanisms.
For monte carlo simulations, we used Madgraph to calculate e + e − → bbg followed by the gluon fragmenting into a a cc pair in either a 3 S We attribute this to a naive model that PYTHIA uses for simulating the radiation of gluons from color-octet cc pairs.
We then considered an alternative simulation approach where e + e − → bbg events are generated using Madgraph, then PYTHIA is used to shower the event to a low scale near 2m c without hadronization. The resulting gluon distribution is then convolved with the analytically calculated NRQCD FFs calculated at the scale 2m c . This procedure yields z distributions that are in much better agreement with our NLL' calculations.
Future work will focus on extending the NLL' calculations to hadron colliders, where the unmeasured jet and soft function recently calculated in ref. [25] must be combined with the FJFs of this paper. It would be of great interest to compare the results of these calculations with data from the LHC on high energy jets with heavy mesons and quarkonia. Finally, there needs to be more work on improving the understanding of the differences between NLL' and monte carlo simulations. Monte carlo simulations that can properly simulate the production of quarkonia within jets will be essential for calculating other jet observables for which NLL' calculations are either unavailable or impractical.

A.1 Evolution of measured and unmeasured functions
The RGEs satisfied by the elements of the factorization theorem are separated into two categories; terms that do depend on the variable τ a and terms that do not. The latter satisfy the following RGE where m F is related to the characteristic scale for the particular function, and Z F (µ) is the renormalization function for F (µ). The coefficient Γ F (α s ) is proportional to the cusp anomalous dimension, Γ cusp (α s ), which can be expanded in α s The non-cusp part, γ F (α s ), has a similar expansion The solution to RGE is given by where the exponents K F and ω F are given in terms of the anomalous dimension, and for up to NLL and NLL' accuracy are given by where r = α(µ)/α(µ 0 ) and β n are the coefficients of the QCD β-function, The RGEs for functions that depend on the variable τ a are of the form 12) and the solution to this equation is given by

A.2 Plus-distribution identities
We begin with the equation (A.14) which can be easily proven using Laplace transforms and the defining equation of the plus distribution, where F (τ ) is defined as which yields The following equations can be derived by setting τ → 0 in eq. (A.14), expanding in ω 2 both sides and matching powers: where we used [18] Θ(τ ) τ 1+ω The first expression is singular as τ a → 0 the second is singular as z → 1 and τ a → 0, but the singularities are regulated by dimensional regularization. Employing the distributional identity 1 and similarly for τ a we find for the divergent terms where P qq is defined in eq. (2.9). The first four terms in this expression are the expected UV poles for the angularity jet function (multiplied by δ(1 − z)), see eq. (3.37) of ref. [34].
In order to simplify this expression we have redefined 4πe −γ E µ 2 → µ 2 , i.e., we are working in the M S scheme. The last term is the expected UV pole in the perturbative evaluation of