Predictions for energy correlators probing substructure of groomed heavy quark jets

We develop an effective field theory (EFT) framework to perform an analytic calculation for energy correlator observables computed on groomed heavy-quark jets. A soft-drop grooming algorithm is applied to a jet initiated by a massive quark to minimize soft contamination effects such as pile-up and multi-parton interactions. We specifically consider the two particle energy correlator as an initial application of this EFT framework to compute heavy quark jet substructure. We find that there are different regimes for the event shapes, depending on the size of the measured correlator observable, that require the use of different EFT formulations, in which the quark mass and grooming parameters may be relevant or not. We use the EFT to resum large logarithms in the energy correlator observable in terms of the momentum of a reconstructed heavy hadron to NLL$'$ accuracy and subsequently match it to a full QCD $\mathcal{O}(\alpha_s)$ cross section, which we also compute. We compare our predictions to simulations in \textsc{Pythia} for $e^+e^-$ collisions. We find a good agreement with partonic simulations, as well as hadronic ones with an appropriate shape function used to describe nonperturbative effects and the heavy quark hadron decay turned off. We also predict the scaling behavior for the leading nonperturbative power correction due to hadronization. Consequently, we can give a prediction for the energy correlator distribution at the level of the reconstructed heavy hadron. This work provides a general framework for the analysis of heavy quark jet substructure observables.


Introduction
Jet substructure observables are playing a significant role in numerous experiments at various energies (e.g. LHC, RHIC). The focus is on precision standard model measurements [1][2][3] as well as searches for new physics [4][5][6][7]. At the same time, jet substructure measurements are increasingly being used as a probe of the Quark-Gluon Plasma (QGP) medium [8]. QCD jets produced from early stage collisions of beam quarks and gluons from two nuclei play a central role in studying the transport properties of QGP. During their propagation through the hot and dense medium, the interaction between the hard jets and the colored medium will lead to parton energy loss (jet quenching) [9][10][11]. There have been several experimental signatures of jet energy loss observed at RHIC and the LHC such as modification of reconstructed jets [12][13][14][15] and jet substructure [16][17][18] as compared to the expectations from proton proton (pp) collisions. Continued progress relies on achieving a deeper understanding of the dynamics of jets, allowing for subtle features in a jet to be exploited. This understanding has progressed rapidly in recent years, both due to advances in explicit calculations [19][20][21][22][23][24][25] as well as due to the development of techniques for understanding dominant properties of substructure observables using analytic [26][27][28] and machine learning [29][30][31][32][33][34][35] approaches. Developments in jet substructure (see [36] for a recent review.) have shown that the modified mass drop tagging algorithm (mMDT) or soft-drop grooming procedure robustly removes contamination from both underlying event and non-global color-correlations, see Refs. [21,22,37], and have been applied to study a wide variety of QCD phenomenology within jets [38][39][40][41][42][43][44][45][46].
Heavy b quark jets play a particularly prominent role in QGP tomography. It was pointed out in [47] that heavy quarks have fundamentally different radiation patterns from light quarks and thus heavy quark jets are expected to be more greatly affected by plasma interactions. Experiments at RHIC [48] and LHC take advantage of these properties to make extensive studies of b jets, and the planned sPHENIX experiment in particular promises improved tracking and flavor tagging capabilities for higher precision measurements on b jets [49]. Thus improved theory calculations of b jet properties are now very timely.
In this paper we take steps to build and apply a theoretical framework to compute the substructure of jets initiated by heavy quarks, starting in vacuum. Much of the necessary framework has already been developed in the context of computations of jet mass for ungroomed [50,51] and groomed [39] top quark jets as well as transverse momentum spectra for B hadrons [52]. We can also build upon the EFT framework for groomed light jets in [43]. Meanwhile [53] built an extensive EFT framework for massive quark effects on measurements in Drell-Yan. We draw upon many of these analyses, organizing their pieces to be applied to a large class of Infrared and Collinear (IRC) safe event shapes known as energy correlation functions [54][55][56] computed for groomed massive quark jets. We find the appropriate EFT to describe these observables for groomed massive quark jets is an elegant combination of those for ungroomed massive jets in [50,51] and groomed light-parton jets in [43], resulting in an EFT very similar to that used in [39] for groomed top jets. We use this EFT framework to smoothly describe the whole spectrum in the 2-point energy correlator called e (α) 2 : • For very large e (α) 2 , the spectrum is predicted by fixed-order perturbation theory in full QCD, to which the predictions of resummed singular terms in e (α) 2 for smaller values will have to matched.
• For intermediate e (α) 2 , the spectrum is sensitive to soft drop grooming effects, but the radiation is sufficiently energetic or at wide enough angles not to be sensitive to quark mass. The e (α) 2 dependence in the spectrum factorizes into a collinear-soft mode [57] and a collinear (jet) mode equivalent to that for massless parton jets. 1 • For smaller e (α) 2 , the spectrum is sensitive to radiation at small enough angles to be affected by the nonzero quark mass. The e (α) 2 dependence in the spectrum factorizes into the aforementioned collinear-soft mode, a matching coefficient to HQET for a heavy quark of mass m b , and an ultracollinear mode [50,51] describing soft radiation from the b quark in its rest frame that is highly boosted into the lab frame, where it looks collinear.
• For very small e (α) 2 , we encounter a nonzero minimum cutoff e (α) 2 ≥ e (α) 2,min imposed by grooming and the finite quark mass. Near this scale, we find the collinear-soft and ultracollinear modes actually merge back into each other, recombining into a single groomed massive quark jet function describing the full e (α) 2 dependence here, and requires a fixedorder matching as well, something not usually required at small values of a jet shape observable. This is a novel feature of applying the EFTs developed in [39,43] that we find for this particular class of observables with massive quark jets.
In the main part of the paper we will explain each of these regions in detail. See in particular Fig. 1 for an intuitive illustration of how these modes arise in each region. We see that it is primarily in the last two regions that the heavy quark mass affects the distribution, causing it to differ from those for light parton jets.
The EFT developed here will thus reveal general properties of the radiation pattern in a massive quark jet as well as the impact of the quark mass in determining the jet substructure. In this paper we apply the EFT to compute energy correlation functions for heavy quark jets produced in e + e − collisions, as a first step towards computations for pp collisions, and ultimately, heavy-ion collisions. The act of grooming, however, should help isolate properties intrinsic to the heavy quark jet itself, reducing differences between e + e − and pp cases, at least in the shape of the energy correlator distributions. An analytic understanding of the jet substructure properties for heavy quark-initiated jets in collisions in vacuum will serve as a comparison baseline for understanding the modification of the properties by the medium. Here, we focus just on two-point energy correlation functions, but the EFT that we develop here is quite general and can be easily extended to other jet substructure observables computed on heavy quark jets.
At the same time, the hope is that by understanding the decay kinematics of heavy quark hadrons, jet substructure observables can ultimately serve as a complementary technique to heavy quark jet tagging. Generalized energy correlators have been shown to be excellent discriminants of light quark vs. gluon jets [58], and we have performed preliminary studies showing promising performance for heavy vs. light jet discrimination [59]. This application will require the computation of energy correlation functions in terms of the light decay products of heavy hadrons, which goes beyond the scope of this paper. Here we focus here on two-point energy correlator functions in terms of the momentum of a reconstructed heavy hadron, assuming it has been identified by other means, such as displaced vertex tracking. We leave the extension to include the effects of B hadron decay to a future publication [59].
In Sec. 2, we introduce the specific class of observables that we wish to use to study heavy quark jet substructure. In Sec. 3, we give details of the factorization theorem (formulated within SCET and HQET) for resumming large logarithms in our jet observable. Sec. 4 gives details about the one loop EFT calculation of the factorization ingredients and anomalous dimensions used for resummation. Sec. 5 discusses the numerical implementation of the full theory fixed order calculation at O(α s ) used to match in the far-tail region of the energy correlator distributions, along with an analytical calculation of the singular limit. Sec. 6 provides the analytical expression for the resummed result to NLL accuracy. We then compare our resummed result matched to the O(α s ) full theory fixed order cross section with Pythia at both parton level in Sec. 6 and hadron level in Sec. 7. In Sec. 7 we also discuss the anticipated impact of B hadron decay on the shape of the distribution. Here we also give a prediction for the scaling of the leading nonperturbative corrections using energy flow operator technique. We conclude in Sec. 8.

Probing substructure of heavy quark jets
In this section we describe the jet substructure observables that are our focus in this paper. We also provide a review of the soft drop grooming algorithm that we use throughout this paper. The general class of substructure observables that we consider fall in the class known as IRC safe observables. This means that these observables are insensitive to arbitrarily soft or collinear radiation or splitting. Of course, in case of heavy quark initiated jets, the mass of the heavy quark provides a cutoff on the collinearity of the radiation. This, combined with grooming (which restricts the softness of the radiation), restricts the minimum value achievable for any jet substructure observable. This is one of the general features that distinguishes massive quark jets from massless ones.

Energy Correlation Functions
We focus on the measurement of the two-point energy correlation function [55,56] on a jet initiated by a heavy quark. The jet is groomed using a soft drop grooming algorithm. For e + e − collisions this observable is defined as where we sum over pairs of particles (i, j) inside the groomed jet. We define the angle (θ ij ) between the two particles with momentum p i and p j as while the energy fractions are z i = E i /E J , E J being the energy of the jet. The energy correlation functions are insensitive to recoil effects [54,60]. At the same time, they do not include explicit axes in their definition.
In the small angle and massless limits, Eq. (2.2) coincides with the geometric angle between the three-vectors p i and p j , which we will denote ϑ ij when we need to refer to it: where we took one of the particles (i) to have a nonzero mass. In the small mass and small angle regimes, we can drop the corrrections, as appropriate. Eq. (2.3), does, however, imply a nonzero minimum for the θ between a heavy quark and soft gluon radiation it emits, for which we can let ϑ go all the way to zero: an effect we will find important to keep. For jets produced in pp collisions, the definition of the energy correlation function is modified to be boost invariant along the beam direction.
where p T J is the transverse momentum of the jet, R ij is defined as, and p T i , φ i and y i are the transverse momentum, azimuthal angle and rapidity of particle i, respectively. In this paper, we focus on the e + e − calculation. However, it was shown in [43], that for jets at central rapidities and in the limit that all the particles in the jet are collinear, the two definitions of the energy correlator function in Eq. (2.1) and Eq. (2.5) are equivalent.

Soft-Drop Grooming algorithm
The modified mass-drop procedure (mMDT) [21,22] or its generalization known as soft-drop [37] removes contaminating soft radiation from a jet, which is originally identified by any suitable algorithm such as (anti-)k t [61][62][63][64][65], by reconstructing an angular ordered tree of the jet constituents by using the Cambridge/Aachen (C/A) clustering algorithm [62-64, 66, 67], and removing the branches at the widest angles which fail an energy requirement. As soon as a branch is found that passes the test, it is declared the groomed jet, and all the constituents of the branch are the groomed constituents. One simply finds the branch whose daughters are sufficiently energetic. Formally the daughters could have any opening angle, though their most likely configuration is collinear. The strict definition of the algorithm is as follows. Given an ungroomed jet, first we build the clustering history by starting with a list of particles in the jet. At each stage we merge the two particles within the list that are closest in angle 2 . This gives a pseudo-particle, and we remove the two daughters from the current list of particles, replacing them with the merged pseudo-particle. This is repeated until all particles are merged into a single parent. Then we open the tree back up. At each stage of the declustering, we have two branches available, label them i and j. We require: where z cut is the modified mass drop parameter, β is the parameter which controls the angularities, θ ij is the angle between i th and j th particle defined in Eq. (2.2), R is the jet radius and E i is the energy of the branch i. (In this paper we will actually just stick to β = 0, for which soft drop coincides mMDT.) If the two branches fail this requirement, the softer branch is removed from the jet, and we decluster the harder branch, once again testing Eq. (2.7) within the hard branch. The pruning continues until we have a branch that when declustered passes the condition Eq. (2.7). All particles contained within this branch whose daughters are sufficiently energetic constitute the groomed jet. Intuitively we have identified the first genuine collinear splitting. For a hadron-hadron collision, one uses the transverse momentum (p T ) with respect to the beam for the condition of Eq. (2.7), We formally adopt the power counting z cut 1, though typically one chooses z cut ∼ 0.1. See [40] for a study on the magnitude of the power corrections with respect to z cut for jet mass distributions.

Factorization in SCET + and bHQET
We work within the formalism of SCET (Soft Collinear Effective Theory) [70][71][72][73][74], which provides an EFT framework for studying IR modes in jet physics. Due to the presence of a heavy quark, we will also need HQET (Heavy Quark Effective Theory) [75][76][77], which is an EFT with an expansion parameter 1/m q , the inverse of heavy quark mass, or more specificaly, boosted HQET (bHQET) [50,51] due to large energy of the b quark jet. Due to additional scales induced by the jet grooming we will also need to use the extension of SCET known as SCET + [57,78].
We will predict the cross section in e (α) 2 defined in Eq. (2.1) measured on jets initiated by a heavy quark to which soft drop grooming has been applied. Schematically, the prediction for this spectrum takes the form where σ sing contains logs of e (α) 2 , which become large for e α 2 1 and which we will resum to give an accurate prediction in this region. It is to σ sing that the factorization in this section applies. Meanwhile, σ ns is the nonsingular (that is, integrable as e (α) 2 → 0) remainder function that matches the singular prediction to the prediction of fixed-order perturbation theory in full QCD, which is accurate for larger e (α) 2 . We will compute σ ns in Sec. 5 and discuss how we smoothly interpolate between the two terms in Eq. (3.1) in Sec. 6.
The precise set of EFT modes that are needed to factorize σ sing depends, as we will explain, on the relative hierarchy of scales amongst e (α) 2 , z cut and m q . In the rest of this section we will consider these various possible hierarchies, identify the appropriate modes, and provide a factorized form of σ sing .

Power counting and modes
An efficient approach for studying jet substructure in a systematic fashion is power counting [26,27], which allows us to determine the parametric scaling of observables. This is determined by the soft and collinear limit of QCD and is a very powerful technique to determine the structure of factorization. We want to calculate the two point energy correlation function Eq. (2.1) on a massive quark jet, which has been groomed with a soft-drop grooming algorithm described in Sec. 2.2. The grooming parameters are z cut and β and we will exclusively work with β = 0 for which soft drop coincides with mMDT. A typical value of z cut used is ∼ 0.1.
For e + e − collisions, we consider events in which we have identified a jet using an appropriate jet algorithm such as anti-k t or Cambridge/Aachen. This jet is then groomed with a soft drop grooming algorithm, which makes it relatively insensitive to MPI (multi-parton interactions) and recoil effects. At the same time, the effect of non-global logarithms is negligible [43].
We consider a jet with energy E J and radius R ∼ O(1) in power counting, so we do not consider resumming logs of R, which could be done using additional modes as in [79][80][81][82]. For 2 ∼ z cut , we transition to the fixed-order σ ns in Eq. (3.1) and computed in Sec. 5.) We also have a hierarchy m q E J , so that the heavy quark is highly boosted. The jet axis is defined by the light-like n direction, defined as n = (1,ẑ) in Minkowski coordinates, whereẑ is a unit 3-vector in the direction of the jet.
Let us see what the measurement of e (α) 2 on the constituents of a groomed jet implies for the momentum scaling of relevant degrees of freedom in the final state. We can divide the modes within the groomed jet into two categories, those which are at a wide angle θ s ∼ 1 relative to the jet axis and those which are collinear to it, θ c 1. In the following discussion it will be useful to make reference to Fig. 1 to see where the relevant modes appear. This figure is analogous to Fig. 2 in [43] for massless parton jets, to which the new ingredients due to the heavy quark mass appear along the line at θ = θ min in Fig. 1. The plot is in log(1/z) and log(1/θ), where z and θ are the fraction of the jet energy z = E/E J and angle from the heavy quark θ ib defined by Eq. (2.2) carried by the radiation from the b quark. In these variables the phase space constraints imposed by soft drop, the jet radius, and the minimum angle due to the nonzero quark mass, are all simple straight lines, and contours of constant e (α) 2 are also straight lines, 2) The precise hierarchies of modes that appears depends on the size of α; the case α < 1, which will be the main focus of our paper, is shown in Fig. 1. We will include the discussion of α > 1 in what follows; the analogous diagram of modes is in Fig. 2.

Wide-angle soft modes
First we consider the wide angle radiation, θ s ∼ 1. Assuming that the heavy quark in the jet carries most of the energy (i.e. z ∼ 1), the contribution of this wide-angle radiation with energy fraction z s to the measurement of e z cut , which holds in the singular region we are considering, this implies z s z cut , so any such radiation that could contribute to the measurement of e (α) 2 would be groomed away. The grooming acts effectively as a veto on soft radiation with z s < z cut and contributes to the normalization of the cross section through a soft function S(E J z cut , R), but does not affect the shape of the e (α) 2 distribution. This soft function describes the contribution of modes whose momentum (p s ) scales as, in light-cone coordinates p = (n · p, n · p, p ⊥ ), The energy of this mode indicates that it is sensitive to the z cut parameter but is groomed away since it fails soft-drop. This mode is unable to resolve any smaller angle structure induced by a measurement of e (α) 2 z cut (e.g. e 2(iii) in Fig. 1), the wide angle radiation could contribute to the measurement; however, in this region, the distribution must be calculated in fixed-order perturbation theory, i.e. σ ns in Eq. (3.1), which we will calculate in Sec. 5.
The factorization of the cross-section at this stage then yields the same form as for massless groomed jets in [43], but with a replacement of the massless jet function by a massive jet function: where J qz is a massive quark jet function that contains all the radiation in the jet that may contribute to the measurement and S(E J z cut ) is the soft function. The jet function J qz is still sensitive to multiple scales, m q , the mass of the heavy quark, the jet energy E J and the grooming parameter z cut . For the range of values e (α) 2,min e (α) 2 and z cut that we will be interested in, there is still a wide scale separation within this function, as we are about to discuss below, and is already illustrated along the e (α) 2(i) contour in Fig. 1. It thus requires further factorization. We will give a formal definition at these later steps of factorization.
The factor σ 0 in Eq. (3.4) contains the Born cross section, the hard matching function H(Q 2 , µ), and an unmeasured jet function Jn(QR, µ) [79,83,84] for the jet in the opposite direction on which e (α) 2 is not measured, as well as another soft factor for radiation between the two jets (which disappears in the case of hemisphere jets R → 1). All of these factors describe effects at much higher energy scales that we integrate out of our EFT, represented by the bottom left corner of the phase space in Fig. 1, and just contribute to the normalization but not the shape of the e (α) 2 distribution; they will all be divided out later.

Hierarchy of collinear modes
We next consider the collinear radiation. Due to the grooming and the quark mass, there are multiple "collinear" scales. The first important observation to make is that there is minimum angle θ ij as defined in Eq. (2.2) that can exist between collinear radiation and the massive quark that initiated it, due to nonzero m q . For the case of massless jets there is no such lower bound. Hence, in such a jet, for a mode with a given energy, the angle is set by the measurement of e (α) 2 alone. On the other hand, for a jet initiated by a massive quark, the angle θ ij defined in Eq. (2.2) between the heavy quark and a light parton has a minimum value as we observed in Eq. (2.4): which comes from the condition ϑ ij ≥ 0 for the geometric angle, and E i ≤ E J for the quark energy, in Eq. (2.3). Since the minimum energy fraction of a light parton in the jet to pass soft-drop needs to be z cut , this automatically sets a lower limit on the values of e Collinear-soft modes Consider first, the widest angle radiation which we could have in the groomed jet. To pass grooming, it must have an energy fraction z z cut . Since the jet is measured to have a value of e (α) 2 , this radiation must have an angle scaling as Given the lower limit Eq. (3.6), on e (α) 2 , we immediately see that this angle satisfies Hence the angular scaling will be determined by the measured e (α) 2 as in Eq. (3.7), and the light-cone components of the momentum of this mode must scale as which, as we have labeled, is the scaling of a collinear-soft mode [57]. It has a soft enough energy to be sensitive to grooming at z cut , and wider angle than the "ordinary" collinear modes we look at next, but still has a degree of collinearity in its angular scaling, due to measurement of small e (α) 2 . It can resolve the phase space boundary around the region labeled S C in Fig. 1, but none others. This mode may or may not pass soft-drop and contributes to the measurement only when it passes soft-drop. This mode is then identical to the one in the case of light quark jets [43] and is insensitive to the mass of the heavy quark. The cut , because z cut < 1 and α < 1, is actually larger than the "ordinary" collinear scale we will consider below. Thus, at this point, we can then further factorize the cross section as in [43], where as explained in [43], the collinear-soft function S c , depends only on the single scale shown, related to its virtuality. Dependence on all lower, more collinear, scales is still in the function J q , which we need to factorize further.
Collinear and ultracollinear modes In Eq. (3.10), J q is the massive quark jet function, defined by the matrix element where n is the light-like direction of the jet and ω = 2E J , where E J is the energy of the jet. The field χ n represents a quark with mass m q moving in the n direction. The measurement operatorê is defined by its action on a collinear state |X n , taking the collinear limit of the full measurement function: The jet function Eq. (3.11) contains all the modes that potentially pass soft drop. Let us consider the modes that make up this jet function. The relevant modes depend on how the angle of the collinear radiation θ c from the initiating heavy quark compares to θ min in Eq. (3.5). There are two cases: Figure 2. Phase space in z, θ and associated modes, for the case α > 1 (in this plot, α = 2), to be contrasted with Fig. 1. In this case, the e itself. In this region, the factorization is the same as for massless parton groomed jets, giving rise to the massless jet function J q . For smaller e

• θ c > θ min
In this case, the scaling of the angle for this mode is set by the measurement and not by the mass of the quark. Its energy fraction is not limited, so z c ∼ 1. This corresponds to measurement values e (α) The light-cone momenta of the mode (p c ) then scale as (3.13) which is insensitive to the quark mass. This region actually exists as a separate region only for α > 1, illustrated in Fig. 2, in the lower region labeled J q . The jet function thus becomes independent of the quark mass and is the same as for a massless jet. In this regime, the factorized cross section then becomes (3.14) which is the same form obtained in [43] for massless quarks.
• θ c ∼ θ min Once the angle θ c hits θ min , it can go no lower, and thus the angular scaling of collinear modes in this regime is fixed. For z c ∼ 1, it contributes e (α) 2 ∼ (θ min ) α . To go to any lower values of e (α) 2 , the only way to do so is to lower the energy fraction of the emitted collinear radiation, down to the lower limit e (α) 2,min in Eq. (3.6), moving along the right-most red line in Fig. 1 or Fig. 2. The measurement of such e (α) 2 then constrains that (3.15) and the mode's momentum scales as where we have defined, This way (Eq. (3.16)) of writing this mode's momentum suggests that it is the momentum of a small fluctuation around a boosted heavy quark at nonzero e Here v µ is the four velocity of the boosted heavy quark, given again in light-cone coordinates. This mode is the ultra-collinear mode of boosted HQET (bHQET) [39,51]. This tells us that we need to match the massive jet function onto a boosted HQET jet function with the heavy quark mass m q (which is now like a hard scale) integrated out. The scale Γ serves as the IR scale for this EFT. It is analogous to the top quark width which served as Γ in [39,51], but here is due simply to the radiation from the heavy quark. We can immediately see that z uc > z cut so that, as expected, this mode automatically passes soft-drop. In this region of e (α) 2 , the cross section Eq. (3.10) factorizes further as where H(m q ) is a Wilson coefficient from integrating out m q to match onto bHQET, and B q is the bHQET jet function [39,51].
Finally, we note that when e 2,min defined in Eq. (3.6), and as illustrated in Fig. 1 and Fig. 2, the collinear-soft and ultracollinear modes actually merge: We compute each of these functions to one loop order in Sec. 4. There is also a fixed-order matching required between the leading-order combination 2,min in the former region but not the latter. We perform this matching in Sec. 6.

Regions of EFTs
Hence the applicable EFT now depends on the value of e (α) 2 relative to the scale (m q /E J ) α , at the threshold where the "ordinary" collinear modes above hit the angle θ min and, below this scale, have to be matched onto bHQET ultracollinear modes.
In this range, radiation collinear to the heavy quark is restricted to have an angle ∼ θ min , and thus the e (α) 2 cross section is sensitive to the mass of the quark. The applicable factorization is Eq. (3.18).
2 < z cut In this range, radiation collinear to the heavy quark has an angle larger than θ min , and is instead restricted by the measurement of e (α) 2 itself. Then the cross section is insensitive to the mass and essentially behaves as for a massless groomed jet, and we have the factorization Eq. (3.14).
• "Region 0": e In this range, radiation collinear to the heavy quark is restricted to have an angle ∼ θ min . At the same time it is sensitive to the grooming parameter z cut . There is no independent collinear-soft mode and we have the factorization Eq. (3.21). These regions are indicated by the bottom right columns in Fig. 3.
A natural question is whether the transition between regimes happens in a smooth manner, especially at the boundary e First of all, we observe that at the value e (α) 2 = (m q /E J ) α , the ultra-collinear and collinear mode have the same scaling At the same time, not surprisingly, the massive and massless jet function defined by Eq. (3.11) will turn out to be identical at this value of e (α) 2 , as we will see in Sec. 4.1. These properties ensure that the distribution in e (α) 2 is indeed continuous at the transition point, as we will see explicitly in Sec. 6. Now, it is possible, depending on the relative sizes of z cut and (m q /E J ) α , that region II simply doesn't exist. This will happen if (m q /E J ) α z cut . For typical values of these parameters, m q ∼ 5 GeV, E J ∼ 100 GeV, and z cut ∼ 0.1, this can only happen if α < 1. And indeed, we choose to consider α ∼ 0.5 in what follows. Then the scales satisfy (m q /E J ) α ∼ 0.2 z cut . In such a case, only region I above exists. This is the situation illustrated in Fig. 1. When e (α) 2 reaches the upper limit of the region, we simply match our massive EFT onto the full theory fixed order cross section at the common scale (m q /E J ) α , without an explicit transition into the massless EFT regime.
In the case that (m q /E J ) α z cut , we would first transition from the massive (region I) to the massless (region II) EFT at e (α) 2 ∼ (m/E J ) α and then match the massless EFT onto the full theory fixed order cross section at the scale z cut . This is illustrated in Fig. 2, where the value of e (α) 2 determines whether J q gets factored further into S C and B + . 3 Fig. 1 and Fig. 2 illustrate clearly that the value of the angular exponent α in the definition of the correlator e (α) 2 determines whether or not one passes through a region II EFT before transitioning to a region I EFT. For us, the more interesting scenario is actually the first, in Fig. 1, with only a region I EFT, since it means that the e (α) 2 distribution is sensitive to the mass over most of its range. So henceforth we will work with α < 1, and thus only have a region I EFT setup. We also note that this choice of α is what made the virtuality of the collinear-soft mode above to be greater than that of the HQET ultra-collinear mode, which is reflected in Fig. 3 make notation in figure consistent with text and the order of our discussion and factorization above. We will see later that this has important consequences for the nonperturbative corrections to the e (α) 2 distribution.

One-loop EFT results
In this section we compute to one-loop order the functions in the factorized cross section Eq. (3.10). All results are computed using dimensional regularization in the MS scheme.

Massive quark jet function
The quark jet function is defined by the matrix element given in Eqs. (3.11) and (3.12). We use the Feynman rules for massive SCET [85].
We have contributions from two real diagrams shown in Fig. 5.
and Here and below we define δ + (p 2 ) ≡ δ(p 2 )θ(p 0 ) for a four-vector p. It is easier to compute these integrals in Laplace space, in which the jet function is defined for which the renormalized(finite) result for the diagrams R a,b yield the results There are two contributions from virtual diagrams (Fig. 5): with Λ 2 being the off-shellness of the quark. To extract the residue of the propagator, we Taylor expand the integrand in Λ 2 and retain the Λ 0 term which gives us the renormalized(finite) result, Combining the results for real and virtual emissions together with the tree-level result, we get the renormalized result These results give us the anomalous dimension for the jet function, Unsurprisingly, since the global soft, collinear soft and the hard function remain unchanged compared to the massless groomed jet, the anomalous dimension is the same as that of the massless jet function. Comparing with [43], we also observe that the massive and the massless jet function have identical values at 1/s ∼ (m/E J ) α , which marks the transition point between the massive and massless regime.

Global Soft
The global soft function is defined as the following matrix element of Wilson lines so the natural scale for this function is µ gs = 2E J z cut . 4 This leads to the anomalous dimension 4 For a finite R, the scale scale would be modified to µgs = 2EJ zcutR. However, since we are working in a regime R ∼ 1, we set R = 1 for this calculation, i.e. we choose not to resum any logarithms in R since they are small, essentially giving the result for a hemisphere jet. The R dependence becomes important in the tail region where we match onto the full theory fixed order cross section and is implemented numerically.

Collinear Soft function
As explained in Section 3.1, the collinear soft function mode is not affected by the mass of the heavy quark. Hence, it remains exactly the same as that for a massless jet. We reproduce the result ( [43]) for completeness. The collinear soft function is defined as the matrix element The W t Wilson line is the same one that appears in the massive quark jet function but is composed of Collinear Soft fields so that The collinear soft modes only contribute to the measurement if they pass soft-drop, which is implemented by theΘ SD term. In Laplace space, we have the result for the renormalized (finite) function is and this leads to the anomalous dimension (4.20)

Boosted HQET jet function
We evaluate the boosted HQET jet function at one loop in two regimes, one in which the ultra-collinear mode automatically passes soft drop, which happens when e We define the jet function in boosted HQET as, v + is the velocity of the boosted heavy quark v + = (m/ω, ω/m, 0 ⊥ ) and the residual momentum that make up the modes of this jet function (known as ultra collinear modes) scales as where Γ = me (α) 2 /(∆) α/2 is the IR scale for this EFT (m being the hard scale). We can also write down the result of the measurement functionê (α) 2 acting on a state with one gluon of momentum k emitted from the heavy quark as, At one loop, the diagrams again are the same as in Fig. 5. We have two contributions which is exactly the same as the integral R a evaluated in the previous section. We have one more diagram which evaluates to which turns out to be the same as R b in Eq.
The matching betweenJ in Eq. (4.9) and B + in Eq. (4.26) then tells us that the Wilson coefficient H in Eq. (3.18) is i.e., all the virtual corrections of the SCET jet function. This also makes sense since it gives us a clean separation of the scales m and me (α) 2 /(4∆) α/2 . The anomalous dimensions for these functions are now given as The matching coefficient agrees with the one obtained in [51]. In fact it is possible to obtain all the results for the observable which is sensitive to the heavy quark mass simply from knowing the collinear soft function and the matching coefficient from massive SCET to boosted HQET (which is independent of the specific observable considered). So an easy way to extend the result to NNLL would be the computation of the two-loop collinear-soft function.

e
(α) We define the jet function in boosted HQET with an explicit soft-drop constraint as, There are two diagrams as before, and for each of them, we can explicitly divide up the phase space in terms of regions that pass or fail soft-drop.
Due to the phase space and measurement constraints, R 1,a is finite and can be directly evaluated in 4 dimensions. R 1,b contains a divergence. The renormalized results for these diagrams gives us α s C F π ln 2 µ mz cut + π 2 24 (4.31) We have one more diagram which can likewise be split up into two pieces.
The renormalized (finite) result is then The result at nonzero e (α) 2 reproduces the one loop fixed order result we will compute later in the singular limit Eq. (5.19). The anomalous dimension contribution (in Laplace space) is entirely from the piece that fails soft-drop which we can verify is the combined anomalous dimension of the unconstrained HQET jet function B + and the collinear-soft function S c as desired by RG invariance.

Consistency of RG equations
At each step of factorization we can verify the consistency of RG equations where the hard anomalous dimension is given as [57,84] (4.38)

Full range of e (α) 2
To predict the distribution over the full range of e (α) 2 , we will need to match the resummed result to a fixed order cross section. The massive EFT is valid in the region e 2 . The usual procedure is to turn off the resummation and match the resummed result to a fixed order result using a profile function. This enables the distribution to make a smooth transition to the fixed order result around e (α) 2 ∼ z cut . We now compute the fixed order cross section at O(α s ). We implement a k t type jet algorithm to isolate two jets (back-to-back) of radius R ∼ 1. At the same time, we implement the soft drop jet grooming algorithm described in Sec. 2.2 to remove soft radiation from the jet. The grooming parameter is z cut which is the fraction of the energy of the hard scale (E J ). Any radiation with E < z cut E J is removed from the jet and does not contribute to the measurement of e (α) 2 . The cross section for this process at one loop, with three final-state particles, is given by the formula: where q = p 1 + p 2 is the total incoming momentum from the e + and e − , with q = (Q, 0, 0, 0) in the CM frame, k 1,2 are the momentum of the outgoing b andb quarks, respectively, and k 3 the momentum of the outgoing gluon. As usual we average over incoming spins and sum over final spins and colors. In Eq. (5.1), δ J is a delta function imposing the restrictions due to the measurement on the groomed jet in the final state, whose complete form we will work out below. Without loss of generality, we can assume the unmeasured jet (in this case k 2 ) to be along the −z direction. We then need the other two final state particles to have k z i > 0 so that they are in the right hemisphere. We need both of these particles to have energy greater than z cut Q/2. The angle between k 1 and k 3 should be less than R. Altogether these constraints determine δ J to take the form: Meanwhile, the phase space integration measure in Eq. (5.1) is given by: Finally, the spin-averaged squared amplitude in Eq. (5.1) is given by where Q b = −1/3 is the b-quark electric charge, and we have defined x i = 2E i /Q. Momentum conservation imposes x 1 + x 2 + x 3 = 2. This expression Eq. (5.4) is exact in the quark mass m b . Though the terms on the second line may appear to be power suppressed by powers of m b /Q, they actually contribute to the singular cross section. Performing as many of the phase space integrals in Eq. (5.1) analytically as possible, we obtain the cross section in the form: . The first term in Eq. (5.5) is exactly the same as that of the massless jet. The second and third terms may appear as power corrections in m 2 /Q 2 but they also contribute to the singular cross-section. As indicated, x 2 in the integrand is function of x 1,3 , and θ J is the set of further constraints imposed by the jet algorithm and jet grooming, which we will now write out more explicitly.
The soft drop condition is implemented as Meanwhile, the (geometric) angle between k 1 and k 3 is given by For the two particles to form a jet, we need the angle between them to be less than R. The lower limit is set by ϑ 13 = 0, which is the collinear limit and we need the mass to be nonzero to regulate this collinear divergence. Since cos ϑ 13 < 1 so, where we have defined ∆ = m 2 b /Q 2 . This equation tells that the ∆ factor under the square root prevents the singularity x 1 + x 3 = 1 or x 2 = 1. We can rewrite Eq. (5.8)) as, The other limit is set by ϑ = R which is not a region of any singularity, hence we can drop ∆ entirely in this case: would fall entirely out of the allowed phase space, and is sensitive to both quark mass and grooming; (ii) an intermediate value where it is sensitive to quark mass but not grooming; and finally (iii) a large value where it is also insensitive to quark mass. We use fixed-order perturbation theory to compute in regions (ii) and (iii). The region where resummation is important is very close to the corner of the quark mass and grooming boundaries in region (i), near the collinear and soft divergences. This can be visualized more clearly in Fig. 1. The only other restriction which regulates soft singularities is the soft drop condition which gives us These four conditions in Eqs. (5.9), (5.10), and (5.11), define the theta function θ J in Eq. (5.5) and define the boundaries of our phase space, which are illustrated in Fig. 6. We compute the phase space integral numerically to give the full theory fixed order cross section. We do this by integrating over the whole allowed area of phase space illustrated in Fig. 6 above each line of constant e (α) 2 , and differentiating the result with respect to e (α) 2 to obtain the differential cross section Eq. (5.5). Looking at the result in Fig. 7, we see that there are several turning points/kinks in the fixed order cross section and it is instructive to try and understand how these points arise and what scales they correspond to. Given the physical constraints, we can get an intuitive understanding about why these kinks exist. Its clear that due to the energy cut-off imposed by soft drop and the angular cut-off imposed by distribution. The massive full theory fixed order cross section (green) is computed numerically by integration over the phase space in Fig. 6 and has kinks in the shape which corresponding to physical scales where quark mass or grooming constraints are crossed. The massive singular distribution is in blue (which is regulated by the quark mass itself) while the massless full theory fixed-order prediction is shown in orange. the mass of the heavy quark, we have a minimum value of e In the singular limit, we keep the leading term up to power correction in z cut , which gives us the scale e (α) 2,min ≈ z cut (m/E J ) α . As the value of e (α) 2 increases, the constraint imposed by the angle due to the mass of the quark will remain relevant as long as the accompanying radiation is soft. When this radiation becomes collinear (i.e. energy fraction z ∼ 1), beyond the value e (α) 2 ∼ (m/E J ) α , the angle constraint will no longer be relevant. This then is the scale at which we expect a kink in the cross section. Since the mass is no longer relevant, the jet behaves as a massless groomed jet upto power corrections in mass, which is indeed what we observe in Fig. 7. The exact scale corresponding to this kink can be calculated by considering the collinear limit which we describe in the next section.
Finally we consider a configuration where the angle θ ij ∼ 1 and the accompanying radiation just passes soft drop (z ∼ z cut ), e (α) 2 ∼ z cut . Beyond this value of e (α) 2 , we expect that the soft drop condition will no longer be relevant and the jet will behave as a massless ungroomed jet. This gives the second kink in the cross section at e (α) 2 ∼ z cut .

Fixed-order collinear limit
We can compute the cross section at one loop analytically in the collinear (low e (α) 2 ) limit. By comparing the singular terms from the expansion of the resummed cross-section and the fixed order collinear limit, we can fix the normalization for our resummed result in order to do the matching to the full theory one loop fixed order result.
Computing the cross section Eq. (5.5) in this limit is most transparent in terms of the variables z, θ: Solving these relations for x 1,3 , we obtain (5.14) In the collinear limit, we only need these relations in the small θ limit, which gives the Jacobian in the collinear limit, Then, the cross section Eq. (5.5) in this limit becomes (5.17) up to terms that will be suppressed by additional powers of m 2 /Q 2 , z cut , or e (α) 2 . The phase space limits come directly from Fig. 1. We perform the θ integral using the delta function, which in order to have a solution, places a tighter upper limit on the z integral, and also imposes the bound e (α) 2 > e (α) 2,min , as is clear from Fig. 1:

Results at NLL
To compute the cross section at next-to-leading-logarithmic (NLL) accuracy, i.e. up to and including terms of order α n s ln n s in Laplace space, see e.g. [86], we need the one and two loop cusp anomalous dimensions as well as the the one loop non-cusp anomalous dimensions. The required one loop results have been given in Sec. 4. The two loop cusp anomalous dimension can be obtained from the literature [87]. This gives us all the ingredients for NLL resummation. We will concentrate on the regime 2 . This depends on the value of α and ω, but is typically true for α < 1. Thus we will only need to make use of the "region I" and "region 0" EFT described in Sec. 3.2.
As we move towards larger values of e (α) 2 , we reach a transition region in which we must smoothly match to the full theory fixed order cross section which includes power corrections in e 2,min , we need to smoothly transition to the "region 0" EFT described in Sec. 3.2, again using a profile function for the renormalization scale.

Resummation at NLL
Consider a function with an anomalous dimension γ which can be written as, where α s is the strong coupling constant, Γ F is the cusp anomalous dimension and γ F is the non-cusp anomalous dimension. We expand these anomalous dimensions as a series in α s . where C i is the appropriate color factor. The one loop (Γ 0 ) and two loop (Γ 1 ) cusp anomalous dimension are given by [87], The resummed result for a function F at the scale µ to NLL accuracy can be written down as where we have evolved the function from µ 0 to µ with which have expansions up to NLL accuracy, and r = α s (µ)/α s (µ 0 ).

Region I
We have to run three functions (S c , B + , H + ) defined in Sec. 4, from their natural scales to the hard scale E J . The path for resummation is shown in Fig. 8. The collinear soft function S c is the same as for the massless case and so its resummation is the same as in Ref. [43]. We do the scale setting in e (α) 2 space. For the three functions that we consider, we write down their natural scales and anomalous dimensions.
• Collinear Soft • HQET matching function • HQET jet function Then the resummed HQET jet function can be written as The argument of B + on the right-hand side indicates that its natural argument, L C defined in Eq. (4.10), is to be replaced by the derivative operator ∂ ω B + , which generates the same logarithms when acting on the evolution kernel on the right (see [86,88,89]). A similar formula will hold for the collinear soft function. Then the final result for the cross section in Laplace space can be written as, Since the resummation of B + and S c involves the Laplace variable s, and we need to go back to e (α) 2 space. To that end we can use We can now write the resummed cross section in e (α) 2 space Since we are working at NLL accuracy, we keep all the terms of O(α s L ∼ 1). We wish to minimize the logs so we make the choice with the scale µ is set to E J . 5 We will vary around these default choices later to estimate theoretical uncertainty.

Region 0
We have to run two functions (B SD + , H + ) defined in Sec. 4, from their natural scales to the hard scale E J . The path for resummation is shown in Fig. 8. The resummation for the H + function is identical to that of Region I. For the soft drop constrained HQET jet function, the natural scale is simply mz cut which is independent of the measurement scale e (α) 2 (or s in Laplace space) so that Following the same process as for Region I, we write the resummed cross section in e (α) 2 space as 2 ) , (6.16) where B SD + is evaluated in a fixed order expansion, given to 1-loop in Eq. (4.34) As before µ ∼ E J and to minimize the logarithms, we set It is clear from this form that the resummation in this region is only acting as an overall normalization factor. The shape of the e (α) 2 distribution is given entirely by the fixed order result. However, it is still necessary to keep the normalization factor so as to smoothly match to the Region I result. Moreover, the resummation exponent in Region 0 can be obtained from that in Region I simply by setting µ

Partonic e (α) 2
spectrum To obtain the distribution over the full range of e (α) 2 , we need to smoothly match Region 0, Region I and the full theory fixed order cross section in the tail. In order to do that, we need to have profile scales that alter or turn off resummation factors while making a transition from one region to another. . This can be added to the cross section to ensure a complete matching between the two regions. This is like fixed-order matching to the full QCD cross section in the tail later.

• Region I -Tail
For e (α) 2 ≥ z cut , (m/E J ) α , we need to turn off the resummation and make a transition to the full theory fixed order cross section. To do that, we turn off all resummation by smoothly taking all the IR scales µ to E J . This then leaves behind the singular part of the fixed order cross section. The matching is completed by adding in the power correction, i.e., the difference between the singular and the full theory fixed order cross section.
The profile scales for the HQET jet function B + and the collinear-soft function S C , which enable us to implement this matching between various regions are shown in Fig. 9. We implement these profile scales piecewise for the transition between the various regions. For the region 0-I transition we can employ a suitably modified hyperbolic tangent function. For example for µ where e . Once again r and t choose the transition value and the rate of transition between regions. The profile for the collinear-soft function can be tailored in exactly the same way.

Comparison with Partonic Pythia
We present the comparison of our calculations, resummed to NLL accuracy and matched to the one loop (O(α s )) full theory fixed order cross section, with Pythia v. 8.219. To make the comparison, we consider parton level Pythia results with hadronization and decays turned off. Fig. 11 shows the comparison for e + e − collisions at a center-of-mass energy of 300 GeV (left) and 700 GeV (right) with α = 0.5 and a z cut = 0.1. Although these are not typical √ s values at any past or existing e + e − collider, these choices in Pythia provide us with enough statistics for jets with energies of 50-350 GeV for which we will present results, and which are typical energies of jets at current colliders like LHC or RHIC.
To do the matching with the full theory fixed order result, we use profile functions for the collinear-soft and ultracollinear scales (Fig. 9) to match the regions I and 0 and also turn off the resummation to smoothly transition to the fixed order cross section. At the same time, we make an estimate of the errors due to a truncation of the perturbative series by varying the  renormalization scale by a factor of 2 above and below the central value as shown in Fig. 9. Each scale is varied independently and the error band generated is added in quadrature to give the final uncertainty estimate. As we see, the error band is of the order of ±10%. To obtain the result at NNLL accuracy, we would need to compute the collinear-soft function at two loops.
We get a good agreement over the whole range of e (α) 2 . As expected, the cross section goes to zero at a finite non-zero value of e (α) 2 . We also show the comparison of the fixed order and resummed cross sections in Fig. 12. As expected, the impact of resummation is significant in the low e (α) 2 region. We also show a comparison with the resummed, purely perturbative result of the massless quark case at NLL accuracy in Fig. 13, as a rough estimate of the effect of the quark mass m b in the resummation region. Since the value of α is small, the non-perturbative corrections set in at a value e (α) 2 ∼ 0.05 for the case of the massless groomed jet. To get a correct picture of the massless spectrum below this value, it is necessary to include non-perturbative  corrections. This plot is only meant as a rough illustration of the size of the impact of nonzero m b in the peak region.

Hadronization and decay effects and comparison with Pythia
We now turn on hadronization in Pythia while still keeping the B hadron decays off. This enables us to gauge the impact of hadronization corrections. Fig. 14 reveals that the dominant effect of hadronization for large enough e (α) 2 is to produce a simple shift in the e (α) 2 distribution, much like event shapes and other jet shapes (see, e.g., [91][92][93][94]). This can be understood as the impact of soft radiation (E ∼ Λ QCD ) from outside the jet, hadronizing the colored partons inside the jet to color singlet bound states. Hence, we expect the energy of the partons to increase by an order Λ QCD , thus shifting the distribution towards higher values of e (α) 2 . The way to incorporate this shift and extract out the nonperturbative physics in terms of measurable parameters is explored in Sec. 7.2.

B decay effects
Next we also turn on the B hadron decay. Comparing it with the purely partonic cross section in Fig. 15, shows that the impact of B decay is significant. While the peak of the distribution appears to be at a similar value of e (α) 2 , there is a five-fold increase in the cross section. From Fig. 15, we can see that there are two effects of decay. One simply changes the shape of the original distribution and the other is an addition of a large number of extra events in each e (α) 2 bin.
The change in the shape is due to change in the radiation pattern for a given event due to the decay of the B hadron. This will change the value of e (α) 2 that is being contributed by the partonic radiation configuration. On the other hand, the large number of extra events can be understood as coming from the zero bin of the partonic distribution. In this bin, at the partonic level, we have a single b quark in the jet, with all the accompanying radiation groomed away since it fails soft drop. Since this event has only one particle (at the partonic level), it will not contribute to a non-zero value of e (α) 2 . As the b quark hadronizes and subsequently decays (into 3 or more particles), the consequent radiation pattern now has the possibility of contributing to a non-trivial e  light-jet discriminant. (See [39] for inclusion of effects of decay products in top quark jet mass distributions.) Even though at this stage we only include effects of hadronization but not the decay of the B hadron, we can still make meaningful comparisons with the experiments where the B hadron four momentum can be reconstructed from its decay products. This can be done using standard techniques for b-tagging such as vertex displacement. Then at the level of the reconstructed B hadron, we can measure e (α) 2 on the jet, which gives a valuable tool for analyzing the radiation pattern in a jet initiated by a massive quark.

Hadronization corrections
To determine the sector which will give the dominant nonperturbative corrections, we look at the scaling of these corrections for each of our factorized matrix elements. We first look at the Region I which has both collinear soft and bHQET functions. The measurement function for the collinear soft function can be written as (Eq. (4.17)), To probe for nonperturbative corrections, we reduce the virtuality of this mode to Λ QCD while maintaining the angular scaling θ cs ∼ (e 2,min ), it is clear that the nonperturbative corrections from the boosted HQET jet function will dominate in the region of measurement. Following the discussion in [94,95], using the definition of rapidity as η = 1 2 ln(n · k/n · k), we can re-write the measurement delta function for the bHQET jet function as Let us come back to the definition of the bHQET jet function, If we consider a regime where k ⊥ ∼ Λ QCD , then we can consider the exponent as a power series in Λ QCD /Γ where Γ is defined in Eq. (3.17). The dominant power correction will then be captured by the leading non-trivial term of this series, Following [94,95], we can write this in terms of the energy flow operator [96][97][98][99][100][101] as where the operator E T (η) gives us (7.11) which then defines f (η) as Unlike for soft power corrections in [94,95], we cannot completely disentangle the η dependence of f and the matrix element in Eq. (7.9), which otherwise would determine P 1 in terms of a single universal nonperturbative matrix element. However we can still extract the scaling of P 1 in m and ∆ using similar manipulations as in [94,95].
In order to do this, we can boost all the operators to the rest frame of the heavy quark. This is basically a boost in the −z direction. We know that the collinear Wilson line W n remains invariant under this boost as a result of the RPI III symmetry of SCET [102]. The vacuum state is invariant under Lorentz transformations. Under our boost, the field h v + transforms to hṽ whereṽ = (1, 1, 0), i.e., the rest frame of the heavy quark. Finally, the energy flow operator transforms as E T (η) → E T (η + η ) (7.13) where e η = (m/E J ). We then have P 1 (s) = dηf (η) 0|hṽW n E T (η + η )|B + X B + X|W † n hṽ|0 . (7.14) Redefining η as η + η , we are now left with P 1 (s) = dηf (η − η ) 0|hṽW n E T (η)|B + X B + X|W † n hṽ|0 , (7.15) which simplifies to P 1 (s) = − s(4∆) α/2 m dηF (η, α) 0|hṽW n E T (η)|B + X B + X|W † n hṽ|0 , (7.16) where F (η, α) = e η 1 + e −2η α/2 , (7.17) which is now independent of both the mass and energy, though it still depends on α and z cut . In Laplace space we can write the leading power correction as  which is an α dependent (but independent of the mass and energy) nonperturbative parameter which we need to extract from simulation/experiment. The extraction can be done at a single arbitrary energy, mass and then can be used for any other energy and mass. In order to extract this parameter, we consider the first moment of the cross section. Assuming the perturbative distribution is normalized so that the area under the perturbative differential distribution is 1, we can then immediately write with the mass of the heavy quark and the jet energy. For the case of the b quark at α = 0.5, we get a value of A(α) = 0.5 GeV.
In the Region III of factorization, there is a single function which is the bHQET jet function with a Soft-Drop constraint (B SD + ). The analysis for the hadronization correction in this region follows in the same manner as the before. However, since this region covers only as small range of e (α) 2 and exists near e (α) 2,min , the contribution to the average value of e (α) 2 and hence A(α) is very small. Hence, to a very good approximation, the value of A(α) is entirely determined by from the region I.

Conclusion
In this paper, we have presented the first analytic EFT resummation for energy correlator distrbutions computed on groomed massive jets, subsequently matched to numerical fixedorder O(α s ) results. We begin with a calculation of the full theory fixed order cross section for the correlator e where grooming and quark mass effects kick in. We provide an intuitive understanding of these features in terms of specific physical configuration of the radiation in the jet. This analysis helps us understand the relevant physical scales in this process and hence the scaling of the radiation in a given range of e is enforced by the mass of the quark and the soft-drop energy cut-off. We conclude that the mass of the quark is relevant in the regime e (α) 2,min < e (α) 2 < (m/E J ) α ("Region I"), while the system behaves as a massless quark groomed jet in the region (m/E J ) α < e (α) 2 < z cut ("Region II") up to power corrections in m 2 /E 2 J . In the tail region, e 2 > z cut , the distribution is essentially that of a massless quark ungroomed jet. This tells us that, depending on the size of α, we need either two or three stages of matching onto EFTs to describe the singular part of the cross section in the low e (α) 2 regime. To that end we utilize an EFT framework combining the SCET + and bHQET formalisms for the heavy quark and the radiation it emits, which is valid for e (α) 2,min < e (α) 2 < (m/E J ) α . We then need to make a transition to the massless quark EFT regime which has been developed in [43], again within the SCET + formalism. We observe that the two EFT descriptions are identical at the boundary value e (α) 2 = (m/E J ) α , so that one EFT naturally transitions into the other at this value, if both regions exist (α > 1). At the same time, we observe that the EFT for the region e to NLL accuracy. The resummed cross section is matched onto the one loop O(α s ) full theory fixed order cross section. The result shows good agreement with partonic Pythia simulation results. Comparing with hadronic Pythia (with B hadron decays turned off) shows a simple shift in the distribution. We then use the energy flow operator formalism to calculate the scaling of the leading nonperturbative correction to the cross section, which can describe this shift in the distribution. An unusual feature of this distribution is that the dominant nonperturbative corrections appear in the bHQET jet function (as opposed to the soft function for a massless jet). Using the Lorentz transformation properties of the bHQET jet function, we can extract out a universal α dependent nonperturbative parameter (independent of mass and energy), which then can be extracted from simulations/experiment. The distribution we have computed gives an accurate prediction at the level of B hadron production, before its decay. In other words, to compare with experiment, it is necessary to reconstruct the momentum of the B hadron from its decay products and then compute the observable using the B hadron four-momentum. This can, in principle be achieved using b tagging techniques such as vertex displacement. Turning on B decay produces a significant deviation from the distribution that we computed in this way, primarily due to a large number of extra events which populate nonzero e (α) However, this is something we leave for the future.
In order to extend the accuracy of our results to NNLL accuracy, the only additional ingredient needed is the two loop collinear soft function. At present, this function is only known for the special case of α = 2, which corresponds to jet mass. Even though we have considered only the two particle correlator, it is clear that the general form of the EFT developed in this paper and the observations made about the relevance of the quark mass will be applicable for other jet shape observables computed on heavy quark initiated jets.
We believe that this is a significant step into a much more detailed study of jet substructure for heavy quark jets. This analysis will serve as a template for comparison with the corresponding cross section for heavy quark jets in a medium (QGP), which should be able to reveal the effects of medium on jet substructure (see e.g. [25]).