Jet Substructure from Dark Sector Showers

We examine the robustness of collider phenomenology predictions for a dark sector scenario with QCD-like properties. Pair production of dark quarks at the LHC can result in a wide variety of signatures, depending on the details of the new physics model. A particularly challenging signal results when prompt production induces a parton shower that yields a high multiplicity of collimated dark hadrons with subsequent decays to Standard Model hadrons. The final states contain jets whose substructure encodes their non-QCD origin. This is a relatively subtle signature of strongly coupled beyond the Standard Model dynamics, and thus it is crucial that analyses incorporate systematic errors to account for the approximations that are being made when modeling the signal. We estimate theoretical uncertainties for a canonical substructure observable designed to be sensitive to the gauge structure of the underlying object, the two-point energy correlator $e_2^{(\beta)}$, by computing envelopes between resummed analytic distributions and numerical results from Pythia. We explore the separability against the QCD background as the confinement scale, number of colors, number of flavors, and dark quark masses are varied. Additionally, we investigate the uncertainties inherent to modeling dark sector hadronization. Simple estimates are provided that quantify one's ability to distinguish these dark sector jets from the overwhelming QCD background. Such a search would benefit from theory advances to improve the predictions, and the increase in statistics using the data to be collected at the high luminosity LHC.


Introduction
The physics program at the Large Hadron Collider (LHC) has reached a very mature stage. Run II is now completed, and ATLAS and CMS each have ∼ 150 fb −1 of 13 TeV data to explore. This data has already taught us a variety of lessons regarding the Standard Model and beyond, but detection of new physics has thus far remained elusive. Given the strong theory motivations provided by, e.g. supersymmetry and/or WIMP dark matter, most signal regions have been developed to target perturbative extensions of the Standard Model, which yield relatively clean, easily interpretable observables. This is made sharp by the notion of Simplified Models [1][2][3], which typically introduce one or two new physics states whose dynamics and interactions can be fully captured via a few additional terms that one adds to the Standard Model Lagrangian. However, not all Standard Model extensions have collider signatures that can be captured in the weakly-coupled Simplified Model framework. A good understanding of the novel signal regions associated with more out-of-the-box ideas is crucial to achieving full coverage when searching for new physics potentially being produced at the LHC.
Of particular relevance here is the idea that the dark matter could be a stable remnant of some new strong dynamics that resides in a hidden sector . It is then reasonable to assume the presence of some non-gravitational connection to the visible sector, such that the hidden sector was in thermal contact with the Standard Model at some point in the early Universe. This could result from a renormalizable interaction involving the Higgs, Neutrino, and/or Hypercharge Portals [29][30][31] or could be due to the exchange of some new mediator. Depending on the properties of the portal, it could be possible to access the hidden sector at the LHC. Furthermore, the dark strong dynamics could obfuscate the resultant signatures, as has been demonstrated concretely through many examples, e.g. lepton jets [32][33][34][35][36][37], emerging jets [38][39][40], semivisible jets [41][42][43][44], and soft bombs [6,14].
All of these examples share a common characteristic: a hard collision can generate a dark sector parton that subsequently undergoes a dark sector parton shower. This often yields a high multiplicity of soft final state particles, smearing out the kinematics of the underlying partons and making it difficult to distinguish the associated signal against large backgrounds. There is a further practical complication due to the fact that these signatures rely on the presence of dark strong dynamics -the theoretical predictions are not nearly as well understood as in the Simplified Model case. As a result, searches for this class of models are usually designed to be very inclusive, avoiding over reliance on details of the modeling. The resulting trade-off between signal significance and systematic error mitigation motivates the work presented here: our goal is to understand the systematic uncertainties associated with making predictions that rely on dark sector strong dynamics. An appreciation of which aspects of the observable can be reliably considered is crucial for the optimization of resulting search strategies.
Specifically, we focus on scenarios where the dark hadrons that result from a dark sector shower promptly decay back to Standard Model hadrons. Our goal is to explore the properties of the resulting jets' substructure, and to quantify the uncertainty inherent to making such predictions. Since substructure is sensitive to a variety of IR effects, such as the dark hadron mass spectrum and hadronization model, our work provides an observable-driven window into the systematic issues associated with making predictions for these strongly coupled dark sector scenarios.
As the use of jet substructure has become routine (see Refs. [45][46][47][48][49][50][51][52][53] for some reviews), many observables have been proposed to distinguish quark and gluons, or to tag boosted objects, and applications to dark sector showers have also been previously explored [43]. Detailed comparisons of parton and hadron level predictions for substructure observables have been performed in the context of the Standard Model, e.g. see the Les Houches 2017 report [54]. Of particular interest here are variables that were designed to be sensitive to the showering history of a jet, since our goal is to find ways to distinguish QCD jets from those that resulted from showering within a dark sector. We are also interested in taking advantage of advancements in analytic calculations that rely on resummation techniques to capture the showering contribution to substructure. To this end, our benchmark observable will be the energy correlation function e (β) 2 [55], where β controls the sensitivity to wide-angle radiation; see Eq. (2.1) below for details. We choose to focus on e (β) 2 since this family of observables is primarily sensitive to the gauge charge of the associated parton in the underlying hard process, which could be our only handle for uncovering dark shower signatures. 1 There is potential concern when predicting the efficiency of jet-substructure assisted searches. The discriminating power of nearly all substructure observables only becomes calculable if large logarithms that can appear in perturbation theory are resummed to all orders. If this calculation is performed using a Monte Carlo generator such as Pythia, only the leading logarithms (LL, defined in Sec. 2) are correctly captured, 1 A number of observables has been considered for the problem for quark/gluon discrimination that are expected to provide superior discirminiation to e (β) 2 . These include both intrinsically IRC unsafe variables like track multiplicity or N 95 [56], and also more complicated IRC safe observables that try to exploit correlations between multiple particles to approximate the behavior of these multiplicity variables [57][58][59]. The distributions of these observables are dominated by non-perturbative corrections, as discussed in more detail in Sec. 2.1. For QCD, this information can be extracted from suitably chosen control regions, while in the case of a new hidden sector, our only recourse is to appeal to phenomenological models whose systematics are challenging to quantify. This implies that the ability to extract meaningful limits using such observables is significantly reduced, and so we will not consider them in this paper. resulting in large expected theory uncertainties, which cannot be quantified by running the generator alone. 2 For QCD studies, such concerns are partially ameliorated by the fact that the parameters of generators are tuned to real data, allowing them to often match the real world better than their formal accuracy would suggest. When looking for physics beyond the Standard Model that we have not yet observed, we have no such recourse. To better address this state of affairs, we take advantage of theoretical technology developed to resum the soft and collinear QCD logarithms that contribute to e (β) 2 at leading and next-to-leading logarithmic order along with modern numerical implementations within Pythia. Sensibly enveloping across the spread of associated predictions will allow us to quantify the systematic error band that is the main result of this work. These error bands can then be utilized to consistently include substructure information into LHC searches for dark sector physics.
Throughout this paper, we assume the dark sector includes n F families of dark quarks which bind into dark hadrons at energies below some dark confinement scale Λ due to a non-Abelian dark SU ( N C ) gauge group. Dark quarks will be produced with large transverse momentum p T Λ such that they shower and hadronize, yielding jets of dark hadrons. We assume that these dark hadrons decay promptly back to Standard Model quarks, 3 yielding QCD-like jets. We then explore the impact on the e (β) 2 observable as we vary the dark sector parameters Λ, n F , N C , and the effect of making the dark quarks massive. In addition, we provide an approximate characterization of the non-perturbative uncertainties associated with dark hadronization by exploring the impact of varying the phenomenological parameters associated with the Lund string model [66]. We then use our error bars to estimate the extent to which dark sector showers can be distinguished from QCD when including the impact of substructure.
The rest of this paper is organized as follows. In Sec. 2, we introduce the twopoint energy correlation function, which will be used as our benchmark substructure variable. We then review how to calculate this observable to next-to-leading-logarithmic accuracy utilizing traditional resummation techniques. Our enveloping procedure that combines the analytic predictions with numerical results derived from Pythia is then introduced, and provides a proxy for the systematic error associated with making a dark substructure prediction. In Sec. 3, we present the extent to which the substructure 2 Automating parton showers beyond leading log and leading color is extremely challenging. Some progress towards formalizing the problem was made in Ref. [60], followed by a numerical approach to address aspects of subleading color in Ref. [61]. For recent progress in automating aspects of nextto-leading-logarithm accurate parton showers, see Refs. [62][63][64], with a recent candidate full proposal in Ref. [65]. 3 Decays to gluons are also in principle possible but to have them dominate the decay rate would require more involved model building.
changes as a function of some of the dark sector parameters: the dark confinement scale Λ, the number of dark colors N C , dark flavors n F , and the dark quark mass m q . In Sec. 4, we explore the effect of varying the parameters that model the dark sector hadronization. In Sec. 5, we estimate our ability to experimentally probe a dark sector jet against the QCD background. We present our conclusions in Sec. 6. In App. A, we detail the expressions that are used to derive the analytic contributions to our systematic error envelopes.

Substructure Observables with Error Envelopes
A large array of jet substructure observables and algorithms have been developed, and are being combined in analyses in increasingly complicated ways. However, the majority of substructure techniques are designed to find evidence of hard processes buried within boosted hadronic events, 4 and as such, most observables are optimized for the identification of distinct multi-prong structures within a jet. A dark sector has no guarantee that it will produce such structure. Instead, we are interested in observables that are sensitive to the structure of the color charge and gauge group of the radiation making up the parton shower. This problem is closely analogous to the problem of quark/gluon discrimination in QCD, and we may look to prior work in this context for guidance [55,56,59,[69][70][71][72][73][74][75][76][77][78][79][80][81][82][83][84][85]. Additionally, we would like to work with infrared and collinear (IRC) safe observables, so that they are perturbatively calculable. This is particularly important for a dark sector search since, unlike the situation for QCD, we have no data from which to extract any of the non-perturbative parameters which are required to make predictions. Thus, there is no way to estimate their uncertainties without resorting to ad hoc empirical models. These two considerations almost uniquely limit us to considering observables which characterize the angular spread of radiation within the jet. A representative choice is the two-point energy correlation function [55], defined as where β is the angular dependence parameter that determines how sensitive the variable is to the angular distribution of the radiation. The jet algorithm determines the constituent particles in jet J that are summed over in Eq. (2.1). In the context of a hadron collider like the LHC, it is most useful to define z i ≡ p T i /p T J and θ ij ≡ R ij /R 0 , where p T J is the total p T of the jet, R ij is the Euclidean distance between the i th and j th partons in the η-φ plane, and R 0 is the jet radius. 5 For brevity, we will usually drop the (β) superscript below when making general statements, and will also refer to the two-point energy correlation function as the energy correlator when appropriate from context. Note that e (β) 2 is equivalent to the C (β) 1 variable introduced in Ref. [55] and widely used in experimental studies [86][87][88].
To build intuition, one can consider a jet with two constituents; in the infrared and collinear limit, the jet mass is given by can be seen as a generalization of jet mass that incorporates arbitrary angular dependence. It is also closely related to the family of jet angularities [89,90], without the need to define a jet axis.
Our essential idea is to calculate the distributions of interest analytically and numerically assuming various approximations, and then use these to determine an error bar such that it spans the range of predictions. First, we review the analytic calculation of the resummed substructure distributions at leading and next-to-leading log order, followed by a brief discussion of the numerical implementation using Pythia. Then, we explain how we combine the various approximations into an error envelope in the context of a QCD calculation. This will set the stage for Sec. 3, where we explore the range of predictions for the substructure distributions resulting from a dark sector shower.

Analytics Using Traditional Resummation Techniques
To understand the robustness of the e 2 distributions, it is useful to explore the range of predictions that result from analytic techniques for calculating the normalized differential cross section. These formulas were derived in Ref. [91], and we present a summary of the main steps for the calculations in App. A. The collinear limit of the leading order e 2 distribution generates a collinear logarithm from the integral over the splitting angle θ and a soft logarithm from the integral over the momentum fraction z. Enforcing the kinematics of two-body momentum conservation with a delta function, we can write down the differential distribution for e 2 by appealing to the definition in Eq. (2.1): 2) 5 For an e + e − collider, a more convenient choice would be z i ≡ E i /E J and θ ij ≡ 2p i · p j /E i E j or the actual Euclidean angle between the i th and j th partons. In the strict collinear limit, all these definitions collapse to be equivalent, and thus only differ in terms that are non-singular in the small e (β) 2 limit. We choose to normalize θ ij by the jet radius R 0 to eliminate the leading dependence on R 0 .
where R 0 is the jet radius 6 and p i (z) is the appropriate parton splitting function for a quark-initiated jet or a gluon-initiated jet, which are given by where T R = 1 2 is the index of the quark representation, i.e., the fundamental representation. These splitting functions encode the divergences associated with a shower that is initiated by the emission of a soft gluon.
In the limit where e 2 1, we can simplify z(1−z)(θ/R 0 ) β z(θ/R 0 ) β by assuming z 1. It is then straightforward to evaluate Eq. (2.2), which yields 2N C and C g = C A = N C are the color factors associated with the jet, and B q = − 3 4 and B g = − 11 12 + n F T R 3C A encode the subleading terms in the splitting functions that arise from hard collinear emissions. Identifying the characteristic logarithm L ≡ ln 1/e 2 , the cumulative distribution at leading order exhibits a characteristic double logarithm in the limit of small e 2 : This shows that perturbation theory breaks down in the limit of small e 2 , so we would like to resum this double logarithm to derive a convergent prediction. The authors of Refs. [94,95] derived a concise expression for the next-to-leading logarithmic (NLL) resummation of the cumulative distribution for recursively IRC safe observables such as e 2 : where the "radiator" R i is given by with R i ≡ dR i dL and κ = zθp T J , and R 1,i is the fixed-order (FO) correction at nextto-leading order, which allows one to match (in the Log-R scheme [96]) between the resummed and perturbative regimes, ensuring the appropriate kinematic endpoint is respected. As such, G 2,i L 2 and G 1,i L are the logarithms appearing in the fixed-order expression (in the collinear limit) that must to be subtracted to avoid double counting the resummed logarithms. Simplifying z(1 − z) to z in Eq. (2.2) is justified by the identical structure of the two collinear limits and is compensated by a suitable combinatoric factor, as further discussed in App. A.
In the context of quark/gluon discrimination, a number of observables have been proposed that seemingly satisfy our property of being perturbatively calculable while claiming to offer improved discrimination over the energy correlation function above [57][58][59]. This comes at a price. Instead of contributions from individual emissions contributing linearly to the observable, each emission's weight depends on the entire shower history. However, this feature also increases the resulting sensitivity to non-perturbative corrections by reducing the parametric suppression of these effects, and until more detailed understanding of these features is available, it is difficult to recommend the use of such substructure variables in situations where these effects cannot be constrained by data. Note that even in the case of the better understood e (β) 2 , the β dependence of quark/gluon discrimination has been measured, and it noticeably deviates from that of the perturbative predictions [97].
An analytic evaluation of R i is possible, although challenging, e.g. see Ref. [98]. The calculation of the resulting efficiencies at NLL due to a cut on e 2 requires evaluating the gauge coupling α s at two-loop order using the CMW scheme [99], such that efficiencies still need to be computed numerically. Another issue is related to α s becoming non-perturbative as the integral is evaluated at low enough scales. To mitigate these complications, we follow the procedure outlined in Ref. [91]: the coupling is only run at one-loop order and is frozen at a "non-perturbative scale" µ NP = 7Λ, where the factor of 7 is an arbitrary choice. This allows us to find a closed-form solution to Eq. (2.7) at the expense of limiting its logarithmic accuracy. We will call this approximate evaluation of Eq. (2.6) the "modified leading logarithmic" (MLL) resummed cumulative distribution with FO corrections. All analytic distributions presented will be MLL+FO accuracy (with the exception of Fig. 1).

Numerics From Pythia
Our analytic expressions have the benefit that they are transparent, in that we can precisely identify the approximations that go into the calculations. However, they do not account for important corrections from, e.g. hadronization or finite quark masses. They also do not provide any way for us to assess the impact of dark sector hadroniza-tion on our prediction. To address these shortcomings, we compare our results for the e 2 observables to those of a Monte Carlo parton shower that models a new confining gauge group. Although all parton showers in common use are formally accurate to leading log, they include various corrections with the goal of modeling certain higher order effects, e.g. see the Monte Carlo Event Generators review in Ref. [100]. It is worth emphasizing here that all such corrections assume QCD, and as such should be revisited in the context of more general confining theories. Specifically here, we simulate events using Pythia 8.240 [101]. We simulate pp collisions at √ s = 14 TeV including initial-and final-state radiation (without multiple parton interactions) for all our events. The signal is generated via a direct portal fromq q pairs to dark sector quarks, and the evolution of the dark sector is implemented in Pythia's Hidden Valley module [10][11][12], including a dark parton shower, hadronization, and decay back to Standard Model states. Events are clustered into anti-k t jets [102] with radius R 0 = 1.0 and e 2 computed for each jet using FastJet 3.3.2 [103], subject to a jet-level cut of p T > 1 TeV.
We will briefly comment on the implementation of the parton shower in Pythia's Hidden Valley module. 7 The underlying physics model is the same as that used for the time-like QCD shower. Showering proceeds via the emission of a dark gluons from both dark quarks and gluons. The dark quarks may be duplicated up to eight flavors, n F , with identical masses and integer spin by default (we set the dark quark spin to be 1/2 in this paper). Running of the dark gauge coupling is included up to one-loop for an arbitrary SU ( N C ) gauge group, assuming massless quarks. Although the functionality to include an arbitrary dark quark mass spectrum is available, we take the masses m q to be degenerate throughout this study. We do not include any states that are charged under both the Standard Model and dark sector symmetry groups, although such states may be considered to extend the range of phenomenological handles in the resulting signal.
A number of aspects of our analytic calculation make its perturbative accuracy greater than that of the Pythia parton shower. Dark gluon splitting into quark pairs is not currently implemented in Pythia; the P q←g (z) splitting function is not singular in the soft limit, and therefore provides contributions beyond LL accuracy. A choice of minimum allowed p T for emissions controls the termination of the shower at low scales. This threshold may be tuned to data in the case of QCD, but for a dark parton shower, this is a parameter that should not be much larger than the confinement scale. Matrix element corrections ensuring the accuracy of parton splitting to one-loop order

Model Showering
Hadronization are included in the QCD parton shower of Pythia but, being model-dependent, not for the Hidden Valley module. Comparing the analytic results to the Pythia predictions will estimate the resulting uncertainties, which are either included or (in the case of the p T cut) have no impact on our analytic results. The dark sector is assumed to confine, and hadronization is implemented via the Lund string model [66], which has some associated parameters whose values are unknown a priori. We explore the consequences of this fact below in Sec. 4. Hadronization proceeds exclusively to dark pions and dark rho mesons, which all decay back to the Standard Model using flat matrix elements (assuming no flavor symmetries leading to stable dark mesons, see e.g. Refs. [42,44]). Table 1 enumerates the relevant parameters discussed along with their default settings.

Error Envelopes
In this section, we describe the procedure used to compute the error envelopes presented in Sec. 3. To capture the "perturbative" theoretical uncertainty associated with these distributions, we combine a number of variations that probe the systematic uncertainties inherent to making dark shower predictions. First, to incorporate uncertainties in the showering step, we capture the range of parton level predictions by comparing the LL order and the MLL + FO order analytics (which we refer to as MLL in the figures). Next, we compare the MLL order analytics and the parton level numerics, i.e., turning off hadronization. Finally, we compare the parton level and the hadron level numerics to account for the effects of hadronization. For events originating from a dark sector shower, we also compare the dark hadron level and the visible hadron level numerics to capture the effects of decaying dark hadrons and their subsequent recombination into Standard Model hadrons. To construct our error bands, we sum the widths of these comparison sub-envelopes in quadrature to produce an averaged final envelope. The results of this procedure when applied to QCD are presented in Fig. 1. Note that for the later plots we show the central value of the envelope merely to guide the eye; this curve does not simply follow from our analytic results. Then in Sec. 4 below, we investigate the uncertainty due to hadronization modeling. The total error band that includes the perturbative and hadronization errors is then used as the input to our search sensitivity estimates for the LHC presented in Sec. 5.
We note that a common approach to calculating a theory uncertainty is to vary factorization, resummation, and (when considering exclusive observables) fragmentation scale parameters by a factor of two away from their canonical choices. This is a way of estimating higher-order terms that have not been explicitly computed by assuming they are dominated by their logarithmically enhanced pieces. The logarithms dominating our distributions are not due to a running effect so that uncertainties in the resummation procedure will not be captured by such an approach. Theoretical uncertainties for resummed calculations typically require more involved multi-scale variational schemes using effective field theory frameworks. 10 The enveloping approach advocated here is designed to incorporate this uncertainty, while also accounting for unknown details of the hadronization and decay properties of the dark sector. Our estimates of perturbative errors are comparable to those of the effective field theory scale variational approaches for QCD, where similar calculations have been done [106]. Depending on the precise treatment of the normalization when taking scale variations, it is possible to find significantly larger errors below the Sudakov peak (e.g. see Fig. 5 in Ref. [107]), where the interplay of constraints from the integrated cross section calculation and breakdown of resummation convergence makes uncertainties particularly sensitive to choice of scheme [108]. However, the resulting effect on signal yields is minimal, since such large uncertainties occur in a vary rapidly falling part of the distribution.
Before showing the results from varying parameters in the dark sector, we note that the analytic approximation for the radiator R i used in our calculations is not continuously differentiable, see Eqs. (A.21) to (A.23). This is a consequence of sharply cutting off the integrals using the non-perturbative scale µ NP introduced in Sec. 2.1 above, which leads to a kink in the second derivative of the radiator R i . To avoid this We show the predictions derived using the MLL analytic calculation, along with the parton and hadron level numerical results from Pythia. Larger (smaller) angular dependence emphasizes the contribution from pairs of partons with larger (smaller) angular distance. The analytic calculations begin to break down for angular dependence values β < 0.5, which is reflected here in the fact that the β = 0.2 curve does not appropriately terminate at the kinematic endpoint.
issue, we follow Ref. [91] and replace this derivative with a discrete approximation: where the choice δ = 1 is an additional source of theoretical uncertainty that is negligible to single logarithmic accuracy. Figure 2 shows the analytic and numeric e 2 distributions for QCD jets across various angular dependence values β. We note the agreement between the analytic and the numeric distributions begin to diverge for low angular dependences β < 0.5. Furthermore, the β = 0.2 analytic distribution does not appropriately terminate at the kinematic endpoint. We conclude that even though we are working in parameter space where the resummation techniques should be a good approximation, the low angular dependence regime of e 2 is not well modeled. For this reason, we will focus our analysis on the behavior of the e to explore the impact of varying β. From the definition of e 2 in Eq. (2.1), we see that increasing β gives greater weight to emissions at larger angular distances in the distribution. Since emissions at large angle within a jet are preferentially softer at large angles, giving lower weight to large-angle emissions leads to e 2 distributions closer to their kinematic endpoint, behavior that is clearly reflected in Fig. 2. Simultaneously, the distribution of e 2 is dominated by emissions in singular regions of phase space, so that lower values of β provide more sensitivity to the structure of the collinear singularity of partonic splitting functions. This comes at the cost of loss of perturbative control. Sec. 2.1 makes clear that the effective coupling in the calculation of e (β) 2 is α s /β and that for values of β 1, perturbative control of the e 2 distribution is lost throughout phase space.

Applying the Predictions
In addition to plotting the normalized e 2 distributions, we will provide a few different ways of presenting the predictions. We will show the cumulative cross section, which is derived by taking the differential distribution and numerically evaluating the following integral, see the bottom row of Figs. 3, 5, and 7: where e 2,max = 1 4 R β 0 . To incorporate the error envelopes, we assume they are fully correlated. In practice, this simply means we compute the upper (lower) error envelope of the cumulative distribution by integrating the upper (lower) edge of the differential distribution. The choice of x cut will be optimized below when we discuss the discovery potential of dark substructure in Sec. 5.
We also provide some quantitative insight into how different the signal and background distributions are using the MLL analytic predictions directly, see Figs. 4, 6, and 8. The left and middle panels of these figures provide two different figures of merit, which give a quantitative sense of how well one could distinguish signal from background, in this case approximated by quark initiated jets. Specifically, on the left we show ROC curves, which are the parametric curve that traces the background rejection 1 − B as a function of signal acceptance S , due to varying a cut on e 2 . The middle panels show the parametric curve for discovery significance S / √ B as a function of signal acceptance S , again due to varying a cut on e 2 . The right panels show the change in the signal rate as a function of the dark sector parameter that is being varied, for a benchmark fixed background factor, which is taken to be 1 − B = 90%. As we will explore in the next section, these various ways of presenting the predictions provide additional insight into the behavior of the e 2 observable across the dark sector parameter space.

Distinguishing Dark Substructure from QCD
Now that we have established a method to estimate the theoretical uncertainties inherent to calculating substructure distributions, we will apply this technology to explore the range of predictions one can expect from a dark sector including error bars. This demonstrates the behavior of the dark sector as a function of its parameters. In particular, we will highlight how the uncertainties depend on the parameters. While we incorporate the effects of hadronization in this section, we set the hadronization parameters to their default values. The results presented here will be combined with an estimate of hadronic uncertainties in Sec. 4, which are computed by varying the non-perturbative parameters. These are then used as the inputs to the estimates performed in Sec. 5, where we study to what extent it is possible to distinguish dark sector showers from QCD via substructure measurements. Note that we have made the simplifying assumption that the QCD background is entirely composed of quark jets in what follows, since the signal will dominate in the central region of the detector. A more realistic study should of course incorporate a more sophisticated modeling of the background. However, since our uncertainties are dominated by the signal modeling, a more careful accounting of the quark/gluon composition of the background should be a subdominant effect.

Λ Dependence
In this section, we explore the dependence on the dark sector confinement scale Λ. The plots shown in Fig. 3 compare the e 2 distribution for a dark-quark-initiated jet against a QCD-quark-initiated jet for a range of confinement scales Λ = Λ QCD compared to the QCD quark background, for two choices of β. As the confinement scale increases, the dark sector distribution shifts toward larger values of e 2 . The larger confinement scale implies that the dark sector coupling is larger than the QCD coupling at the energy scale of the jet. This implies that the peak of the differential distribution occurs at a larger value of e 2 , or equivalently, that the resummation approximation αL 2 ∼ 1 becomes relevant for larger values of e 2 . Therefore, the distribution peaks closer to the kinematic endpoint.
In the bottom row of Fig. 3, we provide the cumulative distribution Σ(x cut ) for the various choices of Λ. For β = 2, the envelope saturates at x cut = 10 −3 for large values of Λ and shifts toward x cut = 10 −4 as Λ decreases. The range of this envelope is 0.22 and insensitive to the size of Λ. Similarly, for β = 0.5, the envelope saturates at x cut = 10 −2 . The envelope range increases as Λ decreases, from a minimum of 0.26 and a maximum of 0.40. As Fig. 4 shows, the discriminatory power of a dark sector signal against a QCD background increases as the dark sector's confinement scale Λ increases. However, this increased discrimination power saturates for large confinement scales Λ 50 GeV. This saturation is caused by freezing the running coupling at the "non-perturbative scale" µ NP = 7 Λ, which we emphasize is a nonphysical prescription designed to obtain a closed-form solution to (2.7). Using the explicit dependence of µ NP on Λ, we can derive a naïve small coupling expansion for the discriminator, This provides a reasonable estimate of the scaling until Λ 5 GeV, when the approximation begins to fail. This can be traced back to the behavior of Eqs. (A.21) to (A.23), from which we infer that as the confinement scale increases, the non-perturbative effects become more relevant for larger values of e 2 .

N C Dependence
In this section, we explore how the substructure depends on the number of dark colors, N C . The set of plots shown in Fig. 5 compare the e 2 distribution for a QCD-quarkinitiated jet against a dark sector-quark-initiated jet for various choices of the number of dark colors N C > 3. As the number of dark colors increases, the β-function for the dark sector gauge coupling becomes more negative, so the scale evolution is faster for the dark sector than for the QCD background. This faster running of the coupling shifts the dark sector distribution toward smaller values of e 2 since α < α s at the scale set by the jet p T . In the bottom row of Fig. 5, we provide the cumulative distribution Σ(x cut ) for the various choices of N C . For β = 2, the envelope saturates at x cut = 10 −4 , regardless of the value of N C . The range of this envelope is 0.22 and insensitive to the size of N C . Similarly, for β = 0.5, the envelope saturates at x cut = 10 −2 . The envelope range increases as N C decreases, from a minimum of 0.36 and a maximum of 0.40.
As Fig. 6 shows, the discriminatory power of a dark sector signal against the QCD background decreases as the dark sector's number of dark colors N C increases. However, this decrease is rather marginal, and saturates for N C ∼ 10. We can understand this behavior analytically, by expanding the LL resummed cumulative distribution Eq. (A.6) to leading order in α. We find that the N C dependence is well approximated by This makes it clear that the discriminator quickly asymptotes as one increases N C , thereby explaining the qualitative behavior in the figures, i.e., that the sensitivity of the observables studied here to the number of dark colors is minimal.

n F Dependence
In this section, we explore how the substructure depends on the number of dark flavors n F . The plots shown in Fig. 7 compare the e 2 distribution for a dark-quark-initiated jet against a QCD-quark-initiated jet for a range of dark flavors with n F > 5; note that we take the number of flavors for QCD to be n F = 5. As the number of dark flavors increases, the β-function for the dark sector coupling α decreases, and in particular the dark sector is no longer asymptotically free when n F > 11 N C 4T R . This implies that the renormalization group evolution is slower for the dark sector than for the QCD background. This impacts the dark sector distribution by shifting it towards larger values of e 2 , since α > α s at the characteristic hard scale of the jet.
In the bottom row of Fig. 7, we provide the cumulative distribution Σ(x cut ) for the various choices of n F . For β = 2, the envelope saturates at x cut = 5 × 10 −4 for large values of n F and shifts toward x cut = 10 −4 as n F decreases. The range of the envelope is 0.24 and insensitive to the size of n F . Similarly, for β = 0.5, the envelope saturates at x cut = 10 −2 , regardless of the value of n F . The envelope range increases as n F increases, from a minimum of 0.34 and a maximum of 0.50. While we are limited by how many flavors we allow the dark sector to have if we want a confining dark sector, the differential distribution shifts toward larger values of e 2 as n F increases.
As Fig. 8 shows, the ability to discriminate a dark sector signal against a QCD background increases as the number of dark flavors increases. Furthermore, this effect increases rapidly as n −1 F . The dark flavor dependence can be estimated by expanding the LL resummed cumulative distribution given in Eq. (A.6) to leading order in the coupling. This yields While naively this implies that we should be able to find regions of parameter space that are very non-QCD-like, the framework breaks down for n F > 11 N C 4T R , because the dark sector does not confine as mentioned above. Practically, Pythia has limited the number of dark flavors one can include to be eight at most. Therefore, we are not able to numerically probe the discriminator beyond this point in parameter space. However, the trend agrees between the numeric and analytic calculations, and follows the analytic estimate in Eq. (3.3) to a good approximation.

m q Dependence
Finally, we explore the impact of varying m q on the e 2 distribution. Since the analytic calculations assume massless partons, we are not in a position to include the analytic contributions to our error envelopes. However, for IRC-safe observables such as e (β) 2 , the mass dependence of our distributions is suppressed as a power of m q / Λ when m q Λ, with the resulting effect on our results being negligible. This is not the case when quark masses exceed the confinement scale since m q then set the scale where the parton shower terminates. In the latter case, an accurate analytic treatment of finite quark masses is challenging, due to the presence of multiple overlapping logarithms of both e 2 and ratios of quark masses and energy scales. As a result, the resummation of differential distributions becomes a more involved procedure, and we will content ourselves with simply providing the results of a numerical study, and will not estimate the error band for different choices of m q .
With a degenerate spectrum, the impact of finite dark quark masses within Pythia is limited at stopping the parton shower from emitting at scales below m q , since the resulting partons would not be able to subsequently hadronize, and treating the color strings as having massive endpoints in the evolution of the Lund string during the hadronization step. Since gluon splitting to quark pairs is not included, potential finite  Table 1. The associated dark hadron masses are 2 m q . Only a numerical study using Pythia is presented. We provide a cubic fit to these distributions to guide the eye. Larger dark quark masses move the peak to higher values due to the resulting cutoff imposed on collinear divergences for emissions from massive quarks. We show the results for two choices of the angular dependence: β = 2 [left] and β = 0.5 [right].
mass effects due to radiation dead cones around additional massive quarks from g → qq splitting play no role. Matrix element corrections in emission, which induce additional mass-dependence in analogous QCD showers, are not included. When the quark masses are above the confinement scale, the hadrons are more akin to quarkonia like the J/ψ or Υ than they are to light mesons like the π or ρ. Hadronization still occurs, since individual dark quarks cannot decay and can only annihilate once they become bound into hadrons. While the properties of these states may be well approximated by perturbative methods, as long as the multiplicity of quarkonia produced is larger than a few, a parton shower is still expected to provide a good approximation of the final state. 11 The result is displayed in Fig. 9, where we compare the e 2 distribution for a quarkinitiated QCD jet against the Pythia distributions for different choices of the dark quark mass m q ; we assume that the dark quarks are degenerate and that the dark hadron masses are 2 m q for simplicity. The other dark sector parameters are set to the default given in Table 1. We see that the peak of the distributions moves to higher values of e 2 as m q is increased. Additionally, we note that the impact on the distributions is not as dramatic as when we varied Λ above (cf., Fig. 3). This can be understood due to the fact that increasing the quark masses for fixed gauge coupling simply acts to cut out more of the IR region of the shower phase space where the sector is becoming strongly coupled. While this has an impact on the resultant multiplicity of dark hadrons that are produced in a shower, their subsequent decay from a higher rest mass to nearly massless QCD hadrons obscured the impact of the specific mass scale set by m q on the observable distribution.

Quantifying Hadronization Uncertainties
The enveloping procedure includes variations among predictions that result from either an analytic or a numerical approach to capture the dominant IR logs that result from showering. When considering sources of systematic uncertainties, it is critical to investigate the irreducible error on predictions due to incalculable strong coupling effects. Specifically, the numerical results rely on a phenomenological model of hadronization. In the case of Pythia, the hadronization step uses the Lund string model [66], which models the physics of confinement by iteratively connecting partons to each other with color strings, and breaking these strings by pair producing quarks from the vacuum when energetically favorable until an equilibrium configuration is achieved. 12 This approach introduces incalculable parameters, which can be tuned to data in the case of real QCD, but must simply be set by hand for the dark sector. It is therefore critical to our goals here to include the uncertainty associated with these choices. As we will show here, hadronization systematics are of the same size as the perturbative ones included in the error envelopes thus far. Clearly, they should additionally be included for searches performed by the experimental collaborations.
The results of varying the hadronization parameters is given in Fig. 10, where all other dark sector model parameters are set to the benchmark values given in Table 1. We then explored the hadronization parameter space to find a choice that resulted in the least (most) number of dark hadrons, which corresponds to the parameter choices aLund = 0, bmqv2 = 2, and rFactqv = 0 (aLund = 2, bmqv2 = 0.2, and rFactqv = 2). The hadronization band in Fig. 10 is then computed by taking the envelope across the result of the default hadronization parameters and these two extreme choices. For reference, we also plot the perturbative prediction, and provide the error envelope as computed above with default hadronization parameters, and we also show the combination of the two envelopes by adding them in quadrature. The largest impact is that the peak of these distributions do shift; this is expected since the position of the turn over is not under robust theoretical control. We see that the variation from hadronization is of the same order as the perturbative uncertainty. 13 We will use the total envelope in the next section where we estimate the impact of non-trivial error envelopes on a mock search for dark sector substructure.

Discovering Dark Substructure
Having quantified the perturbative and hadronization theory errors on the prediction for substructure that results from dark sector showering, we will briefly turn to estimating the impact of including our error envelopes for a search. Our goal here is to simply estimate the discovery potential. Unsurprisingly given existing limits, the subtle nature of the signature and the overwhelming size of the QCD background will imply that additional handles are required to reduce the background by a factor of O(10 5 ) if there is any hope of seeing evidence for dark substructure signals. For example, in models where some of the dark hadrons are stable, a cut on missing energy could play this role. In this case, ignoring the effect of jet to jet fluctuations in the number of unstable mesons, the predictions made above are unchanged, except that the statistics are reduced due to the fact that some particles are missing. We expect the associated theory uncertainty to be a subleading effect. One important mitigating factor is that stringent limits on new physics contributions to QCD distributions already exist from ATLAS [110] and CMS [111]. Since these searches simply look for high p T jets in the final state, the dark jets would fall in the signal region with essentially equal efficiency to QCD jets. Therefore, our first step to quantifying the discovery reach for models that yield substructure from dark showers is to interpret these bounds as a limit on the dark quark production cross section.
We assume the portal to the dark sector can be modeled by a contact interaction: By hunting for deviations in the tails of jet distributions, ATLAS [110] and CMS [111] have derived comparable limits Λ CI 22 TeV. We emphasize that this limit is essentially unchanged for our model, since the searches do not make any cuts on substructure. We convert this limit on the new physics scale into a bound on the production cross section using an implementation of a B − L extension of the Standard Model [112,113] publicly available in the FeynRules [114] model database. We take the Z mass to be large so that the production processq q → Z →¯ q q is well approximated by Eq. (5.1). Events are simulated using MadGraph5_aMC@NLO [115], taking the model parameters to correspond to the lower bound on Λ CI . This allows us to compute the cross section for pp → q q, and we then simply interpret the result as the rate for dark quark production. We implement generator level cuts on rapidity η < 2 and transverse jet momentum p T > 1 TeV. Our dijet backround is produced by all 2 → 2 QCD processes applying the same cuts. This results in a signal cross section σ S = 5 × 10 −5 pb, which can be compared to the enormous QCD background σ B = 13 pb. 14 These cross sections are used to compute the expected number of events for two choices of integrated luminosity; the final Run III data set of 300 fb −1 and the complete high luminosity data set of 3000 fb −1 . These values should be interpreted as the number of events that survive a loose "pre-selection" for the search.
Next we approximate the discovery significance including the impact of both statistical and systematic uncertainties using where S is the number of signal events, B is the number of background events, and δ i are their respective systematic uncertainties. Given the already stringent limits on the production of the dark quarks, it is easy to check that using dark substructure alone will not provide enough discriminating power to beat down the QCD background. Therefore, we will reframe the question in terms of a background reduction factor , which provides an estimate of what one must be able to achieve by incorporating other handles into the search, e.g. missing energy, resonances, and/or displaced objects. 15 To compute , we solve Eq. (5.2) using the substitution B → B: Larger values of correspond to improved discrimination. First, we estimate how large would need to be in order to see a 2σ excess of signal events without a cut on substructure and assuming no uncertainty on the signal production rate and assuming the cut has no impact on signal statistics, see the left panel of Fig. 11. This provides a baseline against which we can compare how much improvement can be obtained using substructure. Next, we include the substructure cut, using the models with varying Λ as a concrete example. We assume the theory error bands on the dark sector distributions are fully correlated, just as we did above when computing the cumulative distributions, e.g. Fig. 3. For the QCD background, there is a wealth of data that is used for tuning and calibration, and as such the systematic error bars can be controlled by leveraging a variety of inputs. For the results presented in Fig. 11, we use the background uncertainty δ B = 0.1 as determined by a recent NNLO calculation [116]. We additionally assume that δ B does not depend on the substructure cut. As a point of comparison, data driven approaches currently yield δ B ∼ 20% [110].
In Fig. 11, we plot the background rejection factor required to achieve a 2σ exclusion as a function of the dark confinement scale Λ, by optimizing a substructure cut for each choice of model parameters. In order to explore the impact of the error envelopes, we provide the result with δ S = 0 in black and δ S = 0 in red, and we also provide the results for β = 2 and 0.5 to investigate varying the angular dependence parameter. We assume either 300 fb −1 or 3000 fb −1 of integrated luminosity, which allows us to explore the scaling as the data set size is increased. 16 Most importantly, we see that a cut on substructure improves one's ability to discover these models, even when the systematic error on the signal shape is included. In particular, taking β = 2 and Λ = 20 GeV the relative change ∆ = 0.9(0.6) for no the range 1.3 − 1.5 [116]. 15 These changes to the model would obviously also impact the limits on signal production rates, i.e., the limit on Λ CI in Eq. (5.1). 16 It is worth noting that the bound on the scale for the contact operator Λ CI will also improve with more data, which is not being taken into account here. error (with error) for 300 fb −1 ; the relative change for 3000 fb −1 ∆ = 1.5 (1.4) for no error (with error). 17 This motivates future work quantifying the error envelopes for a wider variety of substructure distributions that could result from dark sector showers, so that cuts on these variables can be properly incorporated into searches. In particular, it is important to include such systematics when deriving limits on signal parameter space, since the non-trivial error bands can result in more realistic exclusion regions. Finally, we note that for the 300 fb −1 data set, the optimized value of the cut yields a signal region that is statistics dominated. Then when we increase the data set size to 3000 fb −1 , we find that optimal signal region has comparable statistical and systematic errors. We conclude that this subtle signature of dark sector physics is an interesting target scenario for the physics program at the high luminosity LHC.

Conclusions
In this paper, we explored the theory uncertainties associated with making predictions for a scenario where the presence of a new strongly coupled dark sector leaves its imprint on the substructure of QCD-like jets. We focused on the two-point energy correlation function, e (β) 2 . In particular, we quantified the error resulting from perturbative uncertainties associated with truncating to finite order in the logarithmic and gauge coupling expansions. We also explored the uncertainty due to incalculable non-perturbative hadronization effects. Varying the dark confinement scale Λ had the most pronounced impact on the shape of the resulting distributions. We showed e (β) 2 to be relatively insensitive to the number of dark colors N C but observed more striking variations when varying the number of dark flavors n F . We also briefly explored the dependence on the dark quark mass, although we did not provide an error envelope for these distributions due to the technical limitations discussed above.
We then used these error estimates to quantify one's ability to distinguish dark sector jets from the QCD background. We assumed that current bounds on four-quark contact operators apply, which was used to set the production rate for the dark sector. Achieving sensitivity to this subtle signal requires introducing additional handles for the search strategy that could reduce the QCD background by a factor of O(10 5 ) assuming little impact on the signal efficiency. Depending on the model, one could implement a cut on missing energy, a requirement of one or more b-tagged jets, or identification of displaced vertices or resonances -these additional uncorrelated fea-tures could additionally impact the interpretation of the limit on the production cross section, a full exploration of the open parameter space for variations of the base model is an interesting topic for future work. 18 This signature may also provide an interesting target opportunity for model agnostic approaches to new physics searches that rely on machine learning, e.g. [117][118][119][120][121][122][123][124]. While such approaches could mitigate the impact of theory uncertainties on the discovery potential of searches using substructure, the importance of uncertainties in setting accurate limits or extracting model parameters in the case of discovery cannot be ignored. Regardless of these details, this study makes clear that a dedicated search that relies on subtle features in substructure will benefit from the full data set collected at the high luminosity LHC, thereby providing a compelling physics target for future experimental efforts.
Moving forward, we acknowledge the practical need for the generalization of the error envelopes presented here to additional substructure variables. It is important to note that properly accounting for the impact of theory errors for a different observable of interest would require a similar study to what we have presented above. In particular, comparable analytic calculations are necessary to characterize theory uncertainties. We do expect that for a class of mass-like observables, i.e., those that display Casimir scaling at LL [55], one would find conclusions broadly similar to the case of e 2 . However, there are cases, e.g. those briefly mentioned at the beginning of Sec. 2, with a sufficiently different structure, such that a dedicated study would be necessary to determine the size and scaling of the errors. In the case of uncertainties which can be reliably characterized via Monte Carlo alone, e.g. hadronization modeling, modern machine learning methods similar to those of Ref. [125] might prove helpful in reducing the effort involved. However, we emphasize that a proper analytic accounting of expected theory errors in a resummed calculation has no true substitute. The work presented here makes the case that a comprehensive characterization of how substructure observables can be most useful for LHC applications should be performed.

A Analytic Calculation
In this appendix, we review the analytic calculation of the e 2 distribution to NLL order. Our discussion closely follows those of Refs. [55,91], which in turn are based on the framework developed in Refs. [94,95,126]. Our primary goal here is to provide some additional clarification on technical points that may be less familiar to reader not as versed in the details of QCD resummation. For a recent introduction to the general principles of final state resummation accessible to non-experts, see Ref. [127].
We begin with the collinear limit of the e 2 distribution, which is doubly divergent due to a collinear logarithm from the angular integral and a soft logarithm from the integral over the so-called splitting functions. These splitting functions p i (z), which depend on the momentum fraction z can be used to derived resummed distributions. The leading order (LO) contribution is due to a single emission. This can be simply modeled by integrating the splitting function against a delta function that enforces the 2-body momentum conservation as applied to Eq. (2.1). To this order, the differential distribution is where p i (z) is the appropriate parton splitting function for a quark-initiated jet or a gluon-initiated jet, which are given by For quark-initiated jets, only P g←q is included, since the function P q←q is not divergent in the soft limit and would effectively double count the jet core. Likewise, for gluoninitiated jets, the factor of 1 2 multiplying P g←g accounts for a double counting that results from there being the two gluons emerging from a single gluon, while the factor of n F multiplying P q←g provides the proper counting statistics for the gluon to split into n F different quark pairs. In the limit where e 2 1, we can simplify z(1−z)(θ/R 0 ) β z(θ/R 0 ) β by assuming z 1. It is then straightforward to evaluate Eq. (A.1), which yields 2N C and C g = C A = N C are the color factors associated with the jet and B q = − 3 4 and B g = − 11 12 + n F T R 3C A encode the subleading terms in the splitting functions and arise from hard collinear emissions. At LO, the cumulative distribution exhibits a characteristic double logarithm in the limit of small e 2 . Denoting the logarithm as L ≡ ln 1 e 2 , one finds Note that the first integral is divergent, since we have not accounted for virtual corrections. However, we can sidestep this issue by assuming that the probability to emit anywhere is finite. Instead of computing the missing O(α s ) corrections to the total rate, we instead invoke unitary to write the integral in the second finite form which implicitly includes the virtual corrections. Due to presence of the logarithm in Eq. (A.4), perturbative control over the differential distribution is lost for small values of e 2 . Particles with different color charges are going to give qualitatively different behavior in precisely this limit, and so it is necessary to resum the resulting logarithms to all orders to explore how the distributions differ. To leading-log (LL) accuracy, one can consider the emission of n collinear partons within the jet as independent, with the scale of the (one-loop) coupling for each splitting m chosen at the relative transverse momentum scale κ m = z m θ m p T J . Virtual corrections do not change the kinematics, so they will contribute to the distribution for any value of the observable, whereas real emissions will only contribute if the kinematic configuration is such that the emission angle is smaller than the jet radius. At LO, virtual corrections only yield a divergent correction to the tree-level value of e 2 = 0. Thus, to LL accuracy, the resummed cumulative distribution can be computed by simply summing over all emissions off the initial parton while treating them as uncorrelated. In the small z limit, and taking the second form of the integral in Eq. (A.4) to work with finite quantities, the resummed cumulative distribution is given by where the second line sums over virtual emissions, which have the same matrix element as real emissions by unitarity (modulo a sign difference) [94,95]. The series above is readily resummed into a single term, correct to double logarithmic accuracy: The function R i is called the radiator for the jet, and it captures the Sudakov double logarithms associated with the IR divergences that result from soft or collinear emissions from the hard parton. In the fixed coupling approximation, the radiator has the form so that expanding Σ LL to leading order in the radiator recovers the LO behavior in Eq. (A.4). At NLL order a number of new effects appear: multiple emissions, the two-loop running coupling, and non-global logarithms that arise from out-of-jet-emissions falling within the cone. The resummed cumulative distribution can be improved to single logarithmic accuracy by explicitly summing over uncorrelated emissions: 19 The angular ordering condition comes from the fact that when inserting an eikonal emission factor i T i k i · /k i · q into an existing matrix element M, the squared matrix element picks up a kinematic factor of . (A.9) Each such term can be rewritten as W ij = W The benefit of this rewriting is that every such term satisfies an angular ordering property, where R i here is the Laplace space version of the expression in Eq. (A.6). Logarithmic accuracy in ν tracks the logarithmic accuracy in e 2 , since they are Laplace conjugates of each other. Therefore, to derive the NLL cumulative distribution, one must compute the radiator to single logarithmic accuracy in ν. Expanding about ν −1 = e 2 gives 13) where N = 1 + O(α s ) is a matching coefficient that can be determined by comparing with the fixed-order cumulative distribution, γ E is the Euler-Mascheroni constant, and the radiator R i is given in Eq. (A.6).
Note that improving predictability to NLL order requires matching the resummed calculation to the fixed-order distribution. To this end, we implement the Log-R matching scheme [96] by first considering the LO cumulative distribution, i.e., the properly integrated form of Eq. (A.1), where R 1,q = C F β − 4 Li 2 1 + u 2 + 3u + ln 2 (1 − u) − 2 ln(1 + u) ln(1 − u) + 4 ln 2 − ln(1 + u) ln(1 + u) − 3 tanh −1 u + π 2 3 − 2 ln 2 2 , and u ≡ √ 1 − e 2 . Here u takes values between u = 1 − e 2,max = 1 − 1 4 R β 0 and 1. With the Log-R matching scheme, it is straightforward to match the resummed and fixed-order results, where G 2,i L 2 and G 1,i L are the logarithms appearing in the fixed-order expression which must to be subtracted from R 1,i to avoid double counting the resummed logarithms. From Eq. (A.4), these logarithms are explicitly Using this analytic form in Eq. (A.16) requires evaluating the radiator R i , which is given in Eq. (A. 6). An analytic evaluation of R i is possible, although challenging, e.g. see Ref. [98]. The calculation of the resulting efficiencies at NLL due to a cut on e 2 requires evaluating the gauge coupling α s at two-loop order using the CMW scheme [99], such that efficiencies still need to be computed numerically. Another issue is related to α s becoming non-perturbative as the integral is evaluated at low enough scales. Following the procedure in Ref. [91], the coupling is only run at one-loop order and is frozen at the non-perturbative scale µ NP ≡ 7Λ. These choices result in a closed-form solution for R i while limiting its logarithmic accuracy, so Eq. (A.16) provides a modified leading logarithmic (MLL) resummed cumulative distribution with FO corrections. All analytic distributions presented above are to LL or MLL+FO accuracy, but strictly not accurate to full NLL order.
The prescription of freezing the coupling at the non-perturbative scale µ NP ≡ 7Λ leads to the explicit form where µ = µ NP p T R 0 is the relevant scale associated with the non-perturbative transition. Finally, we will write down the explicit expressions for the radiator functions that are used here. Their form depends on the choice of the angular dependence β. For β > 1, while for β < 1, 2 > µ