Dynamical grooming meets LHC data

In this work, we analyse the all-orders resummation structure of the momentum sharing fraction, $z_g$, opening angle, $\theta_g$, and relative transverse momentum, $k_{t,g}$, of the splitting tagged by the dynamical grooming procedure in hadronic collisions. We demonstrate that their resummation does non-exponentiate and it is free of clustering logarithms. Then, we analytically compute the probability distributions of ($z_g$, $\theta_g$, $k_{t,g}$) up to next-to next-to-double logarithm accuracy (N2DL) in the narrow jet limit, including a matching to leading order in $\alpha_s$. On the phenomenological side, we perform an analytic-to-parton level comparison with Pythia and Herwig. We find that differences between the analytic and the Monte-Carlo results are dominated by the infra-red regulator of the parton shower. Further, we present the first analytic comparison to preliminary ALICE data and highlight the role of non-perturbative corrections in such low-$p_t$ regime. Once the analytic result is corrected by a phenomenologically determined non-perturbative factor, we find very good agreement with the data.


Introduction
Jet physics aims at pinning down the microscopic properties of Quantum Chromodynamics (QCD [1]. In the context of heavy-ion physics, the modification of jets with respect to their vacuum counterparts is regarded as an experimental evidence for the formation of a dense, thermal medium, namely the Quark-Gluon plasma [2]. Nowadays, most of the efforts in the field from a theoretical point of view, both from an analytic perspective and with machine learning tools (see Ref. [3] for a review), are directed towards studying the space-time structure of a jet by characterising its radiation pattern through jet substructure observables, i.e. constructed from one (e.g. z g [4, 5]), or a few branchings (e.g. N-subjettiness [6] or the Lund jet plane [7,8]) at most. In this paper we focus on the former category where typically the one branching that defines the observable is selected in a region of phase space where perturbative QCD calculations are applicable, that is, far away from the soft and wide angle sector. This tagging task is handled by so-called 'grooming methods' through which the hard and collinear core of the jet is isolated. The way this general goal is achieved differs from one groomer to the other, e.g. Modified Mass Drop Tagger(mMDT) [9] or its extension Soft Drop (SD) [4] selects the splitting whose momentum sharing fraction obeys z > z cut θ β , while trimming [10] first reclusters a jet into subjets with a smaller radius R sub and then keep only those subjets whose p subjet t > z cut p jet t . In the previous expressions (z cut , β and R sub ) are free parameters that need to be tuned with Monte-Carlo simulations to achieve an optimal performance [11]. These methodological differences leave their imprint into the analytic behavior of the observables that they define [12,13]. For example, as we shall see in more detail in what follows, due to the presence of an explicit z cut in the Soft Drop grooming condition this method is free of non-global logarithms in the resummation function. This fact has enable to push the accuracy of the calculation of Soft Drop groomed observables up to next-to-leading log accuracy in p + p [14][15][16][17][18][19][20][21] and even next-to-next-to-leading log in e + + e − [22][23][24]. It is then clear that the usefulness of a given grooming method should not be judged only on the basis of its resilience to non-perturbative physics, but also on its analytic structure from a pQCD point of view. This paper aims at deepening our analytic understanding of jet substructure observables as defined by a novel grooming technique that has been recently introduced and dubbed 'dynamical grooming' (DyG) [25][26][27].
The dynamical grooming method consists in identifying the 'hardest' branching in a jet tree as a proxy for the physical jet scale. The hardness measure is given by where a is a continuous free parameter that has to be larger than zero in order to guarantee collinear safety. For certain values of a in Eq. (1.1), the hardness measure translates into familiar kinematical quantities, e.g κ (1) = k t , where k t is the transverse momentum of the splitting, or κ (2) = m 2 , with m being the branching mass. In addition, we define θ = ∆R/R where ∆R is the angular separation between the sub-jets and R corresponds to the cone size. The hardest splitting is obtained after re-clustering the jet sample with Cambridge/Aachen algorithm [28] and finding the node with the largest κ in the clustering sequence. In fact, it can be proven through analytical arguments [29] that it is sufficient to look for the hardest splitting along the primary Lund plane of the jet, i.e. following the branch with the larger transverse momentum at each de-clustering step. 1 First steps towards the calculation from first-principles in perturbative QCD of the probability distribution of the momentum sharing fraction, z g , the mass and the relative transverse momentum, k t,g , of the hardest splitting were presented in the original dynamical grooming paper in the resummation region, i.e. when z g (k t,g ) 1 [25]. Interestingly, it was found that similarly to the Soft Drop case, the z g distribution pertains to a special class of jet observables known as Sudakov safe [5,30]. Together with the modified leading-log calculation of DyG jet substructure observables, a Monte-Carlo study of the impact of non-perturbative physics was presented in Ref. [25]. An overall similar performance than Soft Drop was shown, but with a remarkable resilience to hadronization in some cases like the z g distribution as tagged by a = 1. This novel idea has triggered the interest of the ALICE collaboration that has recently conducted some preliminary measurements on the z g , θ g [31], and k t,g [32] distributions at √ s = 5.02 TeV in the jet transverse momentum bin of 60 < p ch t < 80 GeV. As we will see, the low p t reach of the ALICE detector challenges the analytic description of such data set given that non-perturbative effects are sizeable. In addition, first steps towards the experimental use of DyG in heavy-ion collisions were reported in Ref. [33].
From an analytic point of view the purposes of this paper are multifold: (i) understand the resummation structure of dynamical grooming observables and propose a definition for their logarithmic accuracy which circumvents their non-exponentiating nature (double log, next-to-double log, etc) (ii) advance the resummation of z g and k t,g from modified-leading logarithm [25] to next-to-next-to double logarithmic accuracy, 2 as well as presenting for the first time the resummation of θ g , (iii) highlight the absence of clustering logarithms in dynamically groomed observables, (iv) perform a fixed-order matching for all three dynamically groomed observables. This last point is not trivial for the pair of Sudakov safe observables and we propose a novel method to match the resummed and fixed-order distributions. All these ingredients are contained in Section 2. After a few sanity checks on the analytic side, in Sec. 3.1 we compare our results to Monte-Carlo simulations at parton level in a high-p t setup, where non-perturbative effects are mild. Next, in Sec. 3.2, we present the first comparison between an analytic calculation and the preliminary ALICE data for (z g , θ g and k t,g ). In addition, Monte-Carlo studies with different general purpose event generators are performed showing the impact of different details in the definition of the dynamical grooming method in Appendices B, C. The discriminating power of these type of jet substructure observables with respect to different hadronization models and parton showers are shown in App. E. 1 We have numerically checked that our results are robust if we look for the hardest splitting in the whole tree and not only on the primary branch.
2 For a precise definition of N p DL accuracy, see Eq. (2.26) and the discussion below.

Theoretical analysis of dynamically groomed observables
In this section, we present the all-order perturbative calculation of dynamically groomed observables in the κ (a) 1 region 3 and their matching to fixed order results applicable when κ (a) ∼ 1. In the soft limit, z 1 and the (1−z) factor can be removed from Eq. (1.1). Furthermore, as the hardest splitting takes place along the primary branch, we neglect momentum degradation such that p t = p t,jet . Therefore, instead of Eq. (1.1), the definition κ (a) = zθ a is adopted throughout this section, and we sometimes omit the a superscript to lighter the notation. 4

Double-logarithmic estimation and basic properties
We shortly revisit the baseline calculation performed in Ref. [25]. In the κ 1 limit, the two-dimensional probability distribution of a splitting, with kinematic variables (z, θ), to be hardest in the clustering sequence is given by where i indicates the flavor of the jet initiating parton. The two ingredients entering the right-hand side of the previous equation are actually connected through In physical terms, the branching kernel, P (z, θ), represents the probability of a splitting with (z, θ) to occur, while ∆(κ|a) is the so-called Sudakov form factor and vetoes all harder emissions, i.e. those with κ > κ. From Eq. (2.2), it is easy to see that a > 0 is required to regulate the collinear singularity. The normalised probability distribution to measure an observable κ (b,c) = z b θ c on the κ (a) tagged splitting is given by where a sum over flavors including the proper quark/gluon fraction is implicit. The observables that we focus on are obtained from Eq. (2.3) by setting: (b = 1, c = 0) for z g , (b = 0, c = 1) for θ g , and (b = 1, c = 1) for k t,g .
We start by considering branchings in the soft-collinear limit (z 1 and θ 1) that generate terms with powers of α s ln 2 (κ (b,c) ) in Eq. (2.1). That is, we achieve double logarithmic accuracy (DLA) in the language of logarithmic resummation, as we will see below. The soft-collinear limit of the branching kernel reads where P i is the leading-order Altarelli-Parisi splitting function that, in this approximation, is given by with C i being the color factor of the jet initiator parton; C i = C A for gluons, and C F for quarks. The running of the strong coupling is beyond DLA and therefore we fix to its value at the jet scale, i.e. α s ≡ α s (p t,jet R). In this limit, the Sudakov reduces to its opening angle and its relative transverse momentum Location of the peak. An important feature of Eqs. (2.7)-(2.9) is the value at which they are cut off. Its location can be obtained by taking the derivative of e.g. Eq. (2.8) where x ≡ᾱa ln 2 (1/θ g ). Then, the maximum value of the distribution, θ max , satisfies the implicit equation If θ max 1, the left hand side can be approximated by its asymptotic behaviour (x → ∞) such that This equation indicates that the smaller the value of a, the deeper the tagged splitting is on the angular ordered shower, i.e. at smaller angles. In other words, at fixed θ g , larger values of a lead to a bigger Sudakov suppression. Therefore, the distribution shifts to larger θ g for larger a. Thus, Eq. (2.13) confirms and provides an analytic explanation for the result reported in Ref. [25] on the location of the tagged branching in the jet tree using Pythia [34] simulations. Notice that, in order to solve the implicit equation for the peak position, we have assumed that θ max 1. This approximation holds for not too large values of a. Otherwise, the smallness ofᾱ can be compensated by a in the product aᾱ appearing in Eq. (2.13) and θ max ∼ 1.
Following similar steps for the maximum of the momentum sharing fraction, we obtain 5 Again, the previous expression is valid as long as a is not too small.
Finally, in the case of k t,g , we find that This analytic estimate confirms the ordering observed numerically in Figure 9 of Ref. [25], i.e. the a = 0.1 curve is peaked at a smaller k t than the a = 2 case is, being a = 1 the curve peaking at the largest value.
Infra-red and collinear safety. The first step towards boosting the accuracy of our calculation is to analyse the IRC (un)safety of the observables that we are dealing with. As we have already mentioned, and was shown in Ref. [25], dynamically groomed observables are collinear unsafe for a ≤ 0. For a > 0, while k t,g is a standard IRC safe observable, both z g and θ g are Sudakov safe only [5] . This means that the all-order resummation encompassed in the Sudakov form factor regulates the singularities that appear at each order in α s when θ g → 0 or z g → 0. Notice that for θ g , this behavior represents a stark difference with respect to Soft Drop grooming, where this observable is, in fact, IRC safe [4]. This can be understood as a result of the z cut that appears in the Soft Drop condition and regulates the soft singularity. In turn, Dynamical Grooming does not introduce any sharp cut-off on the radiation phase-space and thus nothing forbids the hardest splitting to be in the soft (z ∼ 0) region.
A well-known consequence of Sudakov safety [5,30] is that the z g and θ g -distributions have an ill-defined expansion in (integer) powers of α s . To illustrate this fact, we introduce the cumulative distribution, that defines the probability to measure an observable below a certain value ν, i.e.
(2. 16) 5 Notice that this equation can be also obtained by applying the a → 1/a transformation in Eq. (2.13) The k t,g cumulative distribution at DLA reads 17) and its expansion in α s (or equivalently inᾱ) From the previous expression, it is clear that Σ(k t,g ) admits an analytic expansion inᾱ, as it is expected for an IRC safe observable.
In contrast, the θ g -cumulative distribution is 19) and its expansion in powers ofᾱ is given by In this case, the resumming function is not analytic as the second term in Eq. (2.20) is of order √ᾱ . The non-analyticity on the dynamically groomed θ g is caused uniquely by the √ᾱ term. That is, all other powers ofᾱ appearing in Eq. (2.20) are integer and the function is, in fact, analytic. The same arguments apply to z g where again one can utilise the a → 1/a transformation, to confirm that has an analytic dependence onᾱ at any perturbative order. Interestingly, this α s -expansion of z g is remarkably different from its Soft Drop counterpart when β > 0. For Soft Drop, the expansion is driven by α n/2 s terms where the integer n ≥ 1 [5]. Whether this is a purely mathematical statement, or an explanation in physical terms exists, is beyond our degree of understanding and further work is required to clarify it.
To sum up, in this section we have shown that the opening angle and momentum sharing fraction of the splitting tagged by dynamical grooming are unconventional observables from a pQCD point of view. The Sudakov safety of the z g and θ g distributions leads to an ambiguous definition of the logarithmic accuracy in their resummation, as was noted in Ref. [30]. Furthermore, the standard matching to fixed-order calculations is not trivial due to the non-analyticity of the resummed result. In this context, a careful definition of logarithmic accuracy in the resummation is required and will be provided next.

Revisiting the meaning of accuracy: from IRC to Sudakov safe observables
We start by considering a general resummed formula for an IRC safe distribution obtained with dynamical grooming. Following our previous notation, we denote κ (b,c) the observable that we measure on the splitting whose hardness, κ (a) = zθ a , is the largest in the shower. The cumulative distribution to measure κ (b,c) 1 reads where we have omitted the flavour index for simplicity. An important comment regarding the values of (a, b, c) is in order. Only when b ≤ 1 and c ≤ a, the hierarchy κ (a) ≤ κ (b,c) ≤ 1 is satisfied and thus ∆(κ|a) can be accurately computed through resummation techniques. Any other combination of (a, b, c) leads to a situation in which κ (b,c) 1 does not necessarily imply κ a 1 such that fixed order contributions to ∆(κ|a) become relevant.
Having these constraints in mind, we derive the Sudakov safe distributions of z g ≡ κ (1,0) and θ g ≡ κ (0,1) as two limits of Eq. (2.23), i.e. and The key point is that we define the accuracy of z g and θ g through the accuracy of the IRC safe distribution Σ(κ (b,c) ). For instance, we shall state that Σ(z g ) is known at DLA, if Σ(κ (b,c) ) is known at the same degree of accuracy for all c > 0, or at least in the neighbourhood of c = 0. Our prescription to define the accuracy of Sudakov safe observables follows the spirit of Ref. [30]. However, instead of defining the accuracy of the Sudakov-safe observable by marginalization of an IRC safe double differential distribution, we exploit the IRC safety of the κ (b,c) observable itself. It's important to realise that the perturbative expansion of the Sudakov safe observables is only defined after taking first the appropriate limit on Σ(κ (b,c) ) as given by Eqs. (2.24), (2.25). If these steps are taken in reverse order, i.e. expanding Σ(κ (b,c) ) in powers of α s first and subsequently taking the limit of b(c) → 0, one can show that the correct α s -expansion, given by Eqs. (2.20)-(2.22) at DLA, is not recovered. In short, these two operations do not commute.
The perturbative expansion of Σ(κ (b,c) ) can be written as where the c nm coefficients have to be determined. Then, we adopt the following convention [35,36]: the logarithmic accuracy of Σ(κ (b,c) ) is said to be N p DL if the c nm coefficients are known for all n and 2n − p ≤ m ≤ 2n. Notice that in many other jet substructure calculations it is customary to define the logarithmic accuracy at the level of ln Σ instead of on the cumulative distribution itself. The reason why we use Σ(κ (b,c) ) is because, in general, due to the marginalisation procedure stated in Eq. (2.23) the resummation of DyG observables does not exponentiate [37,38] as it clear from Eq. (2.17). This no exponentiation property is part of other jet substructure observables such as subjet multiplicities. Yet, there is a specific case for which it does: when b = 1 and c = a. That is, when the kinematic variable used for tagging coincides with the measured observable. For instance, select the splitting with the largest k t in the shower, and compute its k t -distribution. In this case, the cumulative distribution is simply the Sudakov form factor, i.e.
A natural question at this point is how does one relate the c nm coefficients with the accuracy of P (z, θ) and ∆(κ|a). In other words, which are the relevant terms that one needs to include in the branching kernel and in the Sudakov form factor in order to reach a given accuracy? To answer this question we rely on the exponentiate property of ∆(κ|a), to write its logarithmic structure in the traditional form 6 [38] with C n being constant coefficients, g i analytic functions and x ≡ α s ln κ. In the N p LL type of counting, the resuming function g 1 would be referred as LL, g 2 as NLL and so on. Our targeted accuracy is N 2 DL in the rest of the paper, with the possibility of keeping sub-leading terms. After expanding Eq. (2.28) in powers of α s we realize that one has to account for the following g nm coefficients at the corresponding level of accuracy DL(p = 0) : g 11 , (2.30) NDL(p = 1) : g 11 , g 12 , g 21 , (2.31) N 2 DL(p = 2) : g 11 , g 12 , g 13 , g 21 , g 22 , C 1 . (2.32) The g 11 was already computed in Sec. 2.1 where we accounted for soft and collinear emissions only The other coefficients and their physical interpretation are provided in the following section up to N 2 DL. Given that the constant C 1 term is related to the interplay between the resummation and fixed-order calculations, we postpone its discussion to Sec. 2.3.2 and neglect it in the resummation-related part.
Turning to the terms that are needed in P (z, θ), we start by working out the plain case (b = 1 and c = a). The exponentiation property of the resummation, in this particular case, leads to a one-to-one mapping between the terms in the Sudakov and in the branching kernel. More concretely, following Eq. (2.23) one gets (2.34) that reduces in the N 2 DL case to where, again, x ≡ α s ln κ. The previous equation, derived exploiting the exponentiation property of the plain case, is sufficient to reach N 2 DL for all values of (b, c). Using Eq. (2.35) with any b and c may produce power suppressed or sub-leading (p ≥ 3) logarithmic corrections in front of α n s ln 2n−2 (κ) terms, that are nevertheless negligible in the resummation region. The physical insight behind the 'universality' of Eq. (2.35) relates to the fact that P (z, θ) is just a probability to have a splitting with a given z and θ. Thus, the branching kernel should be a priori independent of both a and the observable we measure on this branching.

k t,g at LO+N 2 DL accuracy
After this rather formal discussion, we would like to shed light on our statements through an explicit calculation. Namely, we compute the IRC safe k t,g distribution in the small jet radius limit at N 2 DL accuracy on the resummation side and include its matching to a fixed-order calculation at leading order, thus achieving a solid analytic description for all values of k t,g .

Resummation
From the general formula given by Eq. (2.23) it is straightforward to calculate the cumulative k t,g distribution by setting b = c = 1. It reads, such that the differential cross section is where σ 0 represents the Born level total cross-section. In what follows, we calculate the necessary g nm coefficients that enter in the Sudakov form factor and the branching kernel, see Eqs. (2.28), (2.35), and organise them according to the underlying physical effect.
Hard-collinear emissions. Due to its simplicity, the first term that we add to our calculation is the one arising from including hard-collinear corrections (z ∼ 1, θ 1) in the splitting function. This amounts to take into account the finite part of the splitting functions as follows: where B q = 2/C F , B g = 11/12 − n f T r /(3C A ), T r = 1/2, and we fix the number of flavors to n f = 5. The analytic integration of the new finite piece that appears both in the Sudakov and the branching kernel is useful to illustrate the point about sub-leading terms that appear naturally in the calculation. In fact, From the previous equation one can easily read off the g 21 coefficient while the term proportional to B i and no ln(κ) dependence is sub-leading, although might be large when a 1. Strictly speaking, this latter contribution is not required to reach N 2 DL in our calculation, but we will check its numerical impact by the end of this section. Notice that since there is no soft singularity for flavor switching splittings, they contribute as a power correction to κ (a) in our Sudakov form factor and we do not include them here. This argument is valid as long as κ 1 along the lines of the role played by y cut in App. B of Ref. [12].
Running coupling. Up to now, we have fixed the coupling in order to achieve compact, fully analytic expression. However beyond DLA, the running of the coupling has to be taken into account. At 1-loop in perturbation theory it is given by: Next, we integrate analytically the branching kernel with the 1-loop running coupling Note that in the previous expression we have only kept the relevant terms up to N 2 DL, as indicated by the O(N 3 DL) notation. Now, we can identify the terms corresponding to the soft and collinear piece of the splitting function to be while the hard-collinear correction results into In the last equation, the upper subscript in the coefficient indicates that this is not the only term that contributes to the g 22 coefficient, i.e. g 22 = i g i 22 . To achieve N 2 DL accuracy, we need to go to the next order in the running coupling. We work in the CMW scheme [39] which enables to include also the 2-loop contribution of the splitting functions in the soft limit. Then, the running of the coupling at two loops is given by . Again, we can integrate the branching kernel with the 2-loop running coupling to identify another contribution to the g 22 coefficient: An important aspect is that the domain of applicability of Eqs. (2.41) and (2.47) is restricted to perturbative scales, i.e. above µ fr ∼ 1 GeV. Hence, in order to avoid the divergence appearing at such scale, a freezing of the coupling, i.e.
is implemented. We would like to point out that this choice is completely ad-hoc and one could think of replacing Eq. (2.49) by a more general functional form and systematically study its impact on jet substructure observables. We will investigate this possibility in future studies.
Soft emissions at large angles. The dynamically groomed k t,g pertains to the category of so-called non-global observables, 7 i.e. it is sensitive to a certain region of the radiation phase space. As such it is affected by a particularly complex class of logarithms known as non-global logs [40][41][42]. A typical configuration that can give rise to these contributions is a collection of large angle gluons outside the jet which subsequently radiate softer gluons inside. In order to understand how this topology contributes to the k t,g distribution, we first calculate the lowest O(α 2 s ) term coming from such configurations. For illustrative purposes, we start with the calculation of the leading non-global logarithm in e + e − annihilation, in which the color structure of the event is simpler. We discuss the straightforward generalization to p + p collisions in the following paragraph. Once again, we rely on the small-R limit and sketch how to lift this approximation in the next section.
The calculation of the non-global contribution at O(α 2 s ) is standard: one calculates the cross-section for two correlated gluon emissions strongly ordered in energy, with the first emission outside the jet and the second inside. For dynamical grooming, one can rely on the fact that the gluon inside the jet is necessarily the hardest, since it is the only one at this order. Thus, the double differential distribution for having a dynamically groomed z g and θ g value from a non-global configuration initiated by a qq dipole is: where the first Θ-function in the second line enforces the first gluon to be outside the jet, while the second Θ-function constrains the tagged emission to be inside. Note that this is the real term only. 8 The function Ω is the azimuthal average of the real cross-section for correlated double gluon emission from a quark [3]: . (2.51) The R g = θ g R integral in Eq. (2.50) is non-singular in the collinear limit, so that one can perform the two angular integrals exactly to get the leading term in the soft and R → 0 limit: In the previous equation, the soft singularity when z g → 0 induces a single log contribution which has to be taken into account at N 2 DL as part of the g 2 function in Eq. (2.28).
In p + p collisions, the situation is a priori more involved. Each Born level partonic configuration needs to be broken into distinct hard dipoles. However, as shown in Ref. [17], only the dipoles involving the measured jet matter in the small R limit (i.e. neglecting terms proportional to θ n ), and all such contributions are enhanced by the same π 2 /3 factor as in the e + e − result in Eq. (2.52). Consequently, the non-global contribution to the resummed distributions factorize according to the flavour of the jet, in the same way as the collinear piece calculated above. In other words, by imposing the small R limit we can use the e + e − result from Eq. (2.52) for the p + p case. By doing so, we can extract the last piece of the g 22 coefficient, namely At this point let us summarize the main ingredients obtained so far in a compact way that shall facilitate the reproducibility of our results. At N 2 DL accuracy, and in the small-R limit, the Sudakov form factor given by Eq. (2.28) involves the coefficients provided in Table 1.
Physical origin Soft and collinear Hard and collinear Equivalently, using Eq. (2.35) we arrive to the following, non unique expression for the branching kernel In order to estimate the uncertainty of the resummation, we have introduced the dimensionless multiplicative factor µ K that will be varied between 0.5 and 2. A subtle issue 9 concerning this µ K variation is that the first non-trivial correction that arises after integrating over (z, θ) Eq (2.54) is given by ∝ α 2 s ln 2 (κ) ln(µ K ). This is of the same order as the corresponding K-term, i.e. ∝ α 2 s ln 2 (κ)K. To overcome the non-desirable variation of a g nm coefficient, we vary µ K under the condition that the K term is constant, i.e. K is shift to K + 4πβ 0 ln(µ K ) in the calculation (and similarly when considering varition of the renormalization scale Q).
Note that the approximations made to derive the coefficients inside the Sudakov factor ∆(κ|a) and the branching kernel P lead to a differential cross-section which is not necessarily normalized to the Born jet cross-section. To restore the correct normalization, one can simply divide the cumulative distribution Eq. (2.36) by Σ(1). This overall normalization factor is a non-logarithmic correction which does not spoil our targeted accuracy.
Lastly, we would like to draw the reader's attention to the fact that multiple gluon emissions were not considered in this calculation. Generically speaking, at leading logarithmic accuracy, a single emission dominates jet substructure observables. This strong ordering might be broken beyond leading-log, like in the jet mass case, such that an arbitrary number of emissions give comparable contributions to the final measured value. The region of phase space for which this happens has to be determined on an observable basis thus increasing the complexity of analytic calculations. In the dynamical grooming case, multiple emissions do not have to be considered for the observables computed in this paper. This property is a direct consequence of how the method is built. That is, dynamically groomed observables in tagging mode are not additive but defined on the hardest emission and thus it is the only one that contributes to all orders in the resummation. N 2 DL and N 2 DL'. Insofar, we have provided the minimal set of g nm coefficients that lead us to N 2 DL accuracy. For that purpose we have neglected all terms that are not logarithmically enhanced. In order to gauge the impact of these sub-leading contributions, we will also provide results with the 'complete' branching kernel, i.e. 55) and the Sudakov ∆(κ|a) calculated exactly from this complete branching kernel whose explicit expression can be found in App. A. Note that the resulting differential crosssection is then normalized by construction. The running of the coupling is neglected in the non-global term for simplicity, adding it would enable to account for part of the full resummation of the non-global soft function [40]. We will refer to this resummation as N 2 DL', where the prime indicates that the resummation actually includes some of the subleading logarithmic corrections with p ≥ 3. That said, we emphasize that in all rigour, both ways of doing the resummation -either 'minimally' using Eq. (2.54) and the coefficients in Tab. 1 in the Sudakov or with the complete branching kernel given by Eq. (2.55)reach the same N 2 DL logarithmic level accuracy, and not more. 10 Therefore, we will use this resummation scheme freedom to leverage our uncertainty.
Beyond N 2 DL and the small-R limit Before we move on to the fixed-order section, we would like to sketch which steps have to be taken in order to extend the calculation that we have just presented.
In the first place, if the small-R constraint is lifted, one has to account for process dependent terms that enter the calculation as a power series in the jet radius. Physical scenarios that lead to such contributions involve soft and large angle emissions that end up being clustered in the reconstructed jet. For example, a splitting originated from the initial state partons can be tagged by Dynamical Grooming and induce single logarithmic terms suppressed by powers of the jet radius R in the resummation. The difficulty with soft emissions at large angles comes from the fact that such emissions have a complicated color structure which depend on the full Born level event and not only on the Casimir factor of the measured jets. In order to handle such corrections, which are expected to be important for R ∼ 1, one could decide to rely on the large N c limit and decompose each Born processes into different colour flows, as done in Ref. [8,40] in the context of the Lund plane density. Then, each color flow corresponds to a superposition of hard dipoles, which can radiate a soft large angle gluon into the measured jet. In practical terms, adding these contributions would promote the jet flavour dependence of P (z, θ) and ∆ i (κ|a) to a color flow one. Once these new terms are taken into account N 2 DL accuracy is reached beyond the small jet radius limit, but in the large N c approximation. If one does not resort to the large N c limit, one has to deal with matrix formulae in color space, as in Ref. [17]. It is however unclear if the simple structure of Eq. (2.36) remains when the exponentiation has a matrix form and it deserves a dedicated study.
From the non-global logarithms side, their full resummation is required if a higher accuracy in the resummation is intended. This is a complicated task for the κ (b,c) observable in p + p collisions, even in the large N c limit. With the latter approximation, one could resort to the same numerical method as in Ref. [17] (see also Ref. [18] for the θ g distribution defined with Soft Drop).

Matching to fixed-order
In order to produce reliable predictions when k t,g = O(1) and to achieve N 2 DL via the C 1 term, the resummed distribution obtained in the previous section needs to be matched with a fixed-order calculation. Several matching schemes are available in the literature. For our purposes, it is clearly desirable to have a matching scheme satisfying the two following conditions: (i) the matching scheme should produce 'for free' the C 1 term, (ii) the matching scheme should preserve the fixed order endpoint of the distribution at k tg,max = 0.5. Two possible matching schemes that satisfy these requirements are the multiplicative and the log(R) matching [37,43]. In what follows, we shall use multiplicative matching at leading order (O(α s )) and discuss how to extend it to next-to-leading order (O(α 2 s )). As the colour structure of the resummation is tremendously simplified within our targeted accuracy, i.e. it only depends on the jet flavor. The matching formula can be decomposed accordingly as follows [43]: We proceed to describe the ingredients entering the previous equation, except the meaning of the last term that will become clear later on. First, σ 0 and σ 1 are the inclusive dijet cross-section at leading order and next-to-leading order respectively. TheΣ i is the resummed cumulative k t,g cross-section for i-jets, that shares the same endpoint, k tg,max , as the fixed order distribution. At a given accuracy, this is achieved through the following transformation: 11Σ Notice that, besides the shift in the endpoint, we have multiplied the resummed cumulative distribution by σ 0,i in order to ensure thatΣ N 2 DL and Σ LO have the same units. Following up with the pieces entering Eq.
One can check that the limiting behavior of the matched distribution is correct. Indeed, Eq. (2.56) gives back the LO distribution for k t,g ∼ k tg,max . In turn, when k t,g 1 the distribution behaves like where the C 1 term is given by its standard definition: Thus, Eq. (2.59) shows that the matched distribution reproduces the resummed result in its regime of validity. In our calculation, C 1 is a constant up to R 2 -suppressed single logarithmic contributions.
In practice, the LO k t,g differential cross-section and the LO and NLO jet cross-sections σ 0 and σ 1 are obtained using MadGraph5 [44] (in fixed order mode) with CT10nlo PDF set [45]. The factorisation scale for the PDF convolution is set to µ F Q, with µ F a dimensionless factor introduced to estimate the uncertainty relative to this prescription. For a given jet selection [p t,min , p t,max ], a unique generation cut is imposed in the fixed order calculation. Namely, the sum of the transverse momenta of the partons is required to be larger than p t,min and one asks for at least one jet with p t > p t,min /4. We have checked that the resulting cross-sections are insensitive to the precise value of these cuts. The reference value of the strong coupling at the jet scale, α s (µ R Q), is evaluated in the MS scheme. The µ R is a dimensionless factor used to gauge the uncertainty with respect to the renormalization scale. When the p t selection is broad, such as in the ATLAS set-up detailed in the following section, the p t range is divided into smaller bins in which the inclusive jet and k tg cross-sections are calculated. The extension of Eq. (2.56) in this case is straightforward.
The last ingredient in Eq. (2.56) involves the decomposition according to the flavour of the jet. This is done in an IRC safe way for both the k t,g differential and inclusive jet cross-section. More concretely, at LO the jets have at most two constituents. Then, when the jet has zero or one net flavour, the jet is tagged as a gluon or quark jet, respectively. Otherwise, whenever the jet is multi-flavored, i.e. it contains two (anti)-quarks of different flavor, it pertains to what we call the 'else' category. The LO k t,g differential cross-section for these multi-flavored jets goes to zero at small k t,g and contributes to the full match result via the Σ LO else (k t,g ) term in Eq. (2.56). Finally, as for the resummation part, we would like to comment on how to further extend the matching procedure to higher accuracy. In this case, the equivalent of the multiplicative matching formula Eq. (2.56) at NLO can be found in Ref. [43] and it involves the NLO k t,g differential cross-section. The latter can be obtained by generating 3-jet events at NLO with MadGraph, or any other code dedicated to matrix element calculations. Even if there is no conceptual difficulty in promoting our matching to NLO, we postpone it for further studies given that its quantitative impact on the resulting distributions could be as sizeable as the missing power of R suppressed terms on the resummation. Therefore, we believe that these two endeavors should be pursued simultaneously and must be included for refining our phenomenological studies presented in Sec. 3.

Results
Once the analytic framework has been presented, we proceed to show some numeric results for high-p t jets (800 < p t < 1000 GeV) at top LHC energy √ s = 13 TeV with cone size R = 0.4. The central value of the following curves is obtained with µ F = µ R = 1. Further, the error bars are obtained by varying a factor of two the following parameters in the calculation: factorization and renormalization scales through the 7-point rule [46], the parameter µ K that controls the scale at which the strong coupling runs (see Eq. (2.54)) and the freezing scale µ fr used to avoid the Landau pole in Eq. (2.49). Then, we combine the various uncertainties by taking the envelope of all distributions.
In Fig. 1 we present four distributions of k t,g : (i) the fixed-order result, (ii) the resummed result at N 2 DL, (iii) the O(α s ) expansion of the latter and (iv) the matched distribution at LO+N 2 DL. It is clear from this figure that the matching procedure works as expected, i.e. the LO+N 2 DL recovers the N 2 DL result at small k t,g , while it tends towards the leading-order curve in the opposite regime. In addition, the endpoint of the resummation is shifted by the matching procedure to the fixed-order one at k t,g = 0.5. By comparing the fixed-order result and the first term in the α s -expansion of the resummed result, we can get a hint on the size of the O(R n ) logarithmically enhanced terms that we have so far neglected. In fact, the difference between the O(α s ) term of the N 2 DL curve and the exact leading order result converges towards a constant at small k t,g . This indicates that these power suppressed terms enter with a small coefficient in the cumulative distribution and can be safely neglected for the setup studied in this work. All the previous statements hold for both values of a. In particular, the fixed order result is independent of a because there is only one splitting tagged. Next, we compare in Fig. 2 the two prescriptions to perform the resummation that we have discussed above, i.e. keeping uniquely the logarithmically enhanced terms at N 2 DL or including sub-leading corrections (N 2 DL'). In the large k t,g regime, we observe no difference between the LO+N 2 DL and the LO+N 2 DL' as it is expected since in this limit the fixedorder contribution dominates the matched result. This is no longer the case for k t,g 1, where details of the resummation structure do matter. An important remark is that the discrepancy between the two curves diminishes when increasing the parameter a that determines the hardness condition in the grooming algorithm. We attribute this to the a−scaling of the g nm parameters that, as one can see in Tab. 1, satisfies g nm ∼ 1/a. Hence, the larger the value of a is, the smaller the coefficients in front of the higher order terms are and the narrower the difference between N 2 DL and N 2 DL' becomes. In the phenomenological section, we will include the differences between N 2 DL and N 2 DL' as part of our uncertainty band given that, from a logarithmic counting point of view, there is no preferred option.
2.4 z g , θ g at LO+N 2 DL accuracy As discussed at length in Sec. 2.1, the momentum sharing fraction, z g , and opening angle, θ g , of the splitting tagged by dynamical grooming are Sudakov safe observables. In Eqs. (2.24)-(2.25), we defined their distribution as the limit of the IRC safe k t,g distribution. This allows us to follow the same logic as in the previous section to obtain the resummation part of their distribution. In turn, the fixed-order result is not even well defined, as shown in Eqs. (2.22) and (2.20), and thus the matching strategy differs to that presented in Sec. 2.3.2. In what follows, we provide the necessary ingredients to reach LO+N 2 DL accuracy in the small-R limit.

Boundary logarithms for the θ g distribution
In the case of z g , the resummation proceeds in exactly the same fashion as for k t,g . In turn, for θ g , another source of logarithmic enhancement appears, caused by the interplay between the anti-k ⊥ algorithm [47] used to cluster the jet, and the C/A algorithm to decluster it in the dynamical grooming procedure. These so-called clustering boundary logarithms [8] are of the form: where the double logarithmic enhancement becomes important when R g → R (θ g ≡ R g /R close to 1).
Boundary logarithms arise from a non-global configuration where the first emission is outside the anti-k ⊥ jet, while the second is inside and almost collinear to the first one. In this situation, C/A algorithm would cluster these two emissions together, as part of the jet. As such, their leading contribution to the branching kernel can be obtained from the same calculation as the one done in Eq. (2.50), but focusing in the regime where θ 1 θ 2 R 1 [8,48], such that: and 1 σ 0 The upper boundary in the θ 1 integral can be safely sent to ∞ since the integral is dominated by the region θ 1 θ 2 R. To get the last line, we have kept the dominant contribution at θ g ∼ 1.
As stated above, this O(α 2 s )-contribution is enhanced by two soft logarithms of the type ln(z g )/z g and one boundary logarithm ln(1−θ g ). However, the logarithmic divergence associated with θ g ∼ 1 is integrable in a neighbourhood of 1. Consequently, as long as one deals with an observable in which the angle R g ∼ R is integrated out between some lower bound and 1 (such as z g or k t,g distributions), the boundary divergence is harmless. More precisely, its integral over θ g should give back the soft single logarithmic divergence that was part of our treatment of non-global configurations. This argument also applies to the Sudakov factor, since vetoing all emissions with hardness larger than κ translates into the following integral where it is clear that θ is always marginalized in the neighbourhood of 1. In other words, the dynamical grooming procedure lowers the singularity associated with boundary logarithms from double-to single-log, and this single-log term is already taken into account by the coefficient g 3 22 in the Sudakov. From that perspective, the θ g distribution is peculiar since the R g → R logarithmic divergence from Eq. (2.64) is not integrated out. To effectively include boundary logarithms for the θ g distribution, we replace the last term in the branching kernel P (z, θ) given in Eq. (2.54) by withθ defined such that (2.67) Numerically, one findsθ 0.66. The step function guarantees that the single logarithmic term from soft non-global configurations is correctly accounted for within our targeted accuracy and without double counting. 12 Such a constraint is also physically expected since boundary logarithms come from the region where θ g ∼ 1 by definition.
Finally, we would like to discuss how this new logarithmic divergence affects the logarithmic counting provided in Sec. 2.2. For the θ g distribution, we have found two sources of logarithmic enhancement that are either of the form ln(θ g ) or ln(1−θ g ). Since the veto factor in ∆(κ|a) suppresses boundary logarithms, there is only one power of ln(1−θ g ) that appears in the α s expansion of the θ g distribution and it comes from the α 2 s result given by Eq. (2.66). In order to have the correct logarithms at N 2 DL in front of this single 12 The spirit of the step function is essentially the same as in our treatment of hard collinear emissions via the effective splitting function given by Eq. (2.38) Another way of including boundary logarithms without double counting is to use directly Eq. (2.63) (without step function) since − 1 0 dθ ln(1−θ 2 )/θ = π 2 /12. power of ln(1−θ g ), it is enough to solely include the first one-loop correction in the running coupling α 2 s → α 2 s (1−4α s β 0 log(z)) and the hard-collinear correction at fixed coupling 1/z → Θ(e −B i −z)/z in Eq. (2.66).

Comparison between the resummation structure of Soft Drop and Dynamical Grooming
The idea of studying the momentum sharing fraction and opening angle of a given splitting in the shower was originally proposed in Ref. [4]. In this work, the splitting at issue was selected through the Soft Drop procedure, that is, the first branching in the de-clustering sequence that satisfies z > z cut θ β . These observables, (z g , θ g ), have been measured experimentally [31,49] and resummed to modified-leading log [50] and next-to-leading log accuracy [18], respectively. An important comment at this point is that Soft Drop observables do exponentiate and, therefore, a N p LL counting applies. Hence, strictly speaking, an apples-to-apples comparison on the resummation structure for Soft Drop and Dynamical Grooming does not exist.
We have identified one major simplification in the resummation structure of θ g when it is defined through the dynamical grooming procedure instead of with Soft Drop: dynamical grooming is free of clustering logarithms. Let us briefly recap how these contributions arise for two correlated emissions [42]. Consider the emission of a gluon, p 1 , off a hard quark p 0 together with a secondary emission, p 2 , off p 1 . These two emissions have commensurate angles θ 01 ∼ θ 02 , while their energies (and thus transverse momentum) are strongly ordered z 2 z 1 (k t,2 k t,1 ). The C/A algorithm will miss-cluster the secondary gluon as a primary if θ 02 < θ 12 , with θ 12 being the relative distance between the two emissions. Then, if p 2 is a real emission it will trigger the Soft Drop condition, even though z 1 z 2 , and consequently θ 02 = θ g . In turn, if p 2 is virtual, the tagged splitting would be p 1 and θ 1 = θ g . This mismatch between the real and virtual contributions lead to a tower of logarithms at NLL that were numerically computed, in the large N c -limit, in Ref. [18] for θ g , as defined by Soft Drop, and also discussed in the context of the Lund plane in Ref. [8]. In the Dynamical grooming case, even if some secondary emissions can be 'wrongly' pushed by the C/A algorithm into the primary Lund plane, it will never be the hardest given that z 1 z 2 and, therefore, κ 1 = z 1 θ a 1 will be larger than κ 2 = z 2 θ a 2 even if the angles are commensurate. The only effect of these emissions on DyG, beyond N 2 DL, would be a small contribution to the p t degradation of the primary branch.
From the non-global logarithms side, that affect not only z g but also k t,g , we have shown in Sec. 2.3, that they are proportional to ln(z g ) for Dynamical grooming. In the Soft Drop case, the soft singularity is cured by the definition of the grooming condition. That is, non-global logs enter in the Soft Drop calculation as ∝ ln(z cut θ β ) and thus have a smaller impact that in the DyG option. Therefore, the cleanest grooming procedure from a theoretical point of view in order to avoid the resummation of both non-global and clustering logarithms at NLL would be to combine the two methods. Then, the grooming procedure would be a two-step process: first, one removes all emissions with z < z cut and then one looks for the hardest one in the Dynamical Grooming sense. This possibility will be further studied in an upcoming publication [51].

Matching to fixed-order
The first leading order matching scheme for Sudakov safe observables was proposed in Ref. [30]. It is based on constructing a n-dimensional IRC safe distribution, that we dub 'IRC safe companion', and re-defining the Sudakov safe observable by an appropriate marginalization. In our case, the 2-dimensional IRC safe distribution would be d 2 σ/dz g dθ g , that can be interpreted as the joint probability distribution for having a tagged branching with momentum sharing fraction z g and (normalised) opening angle θ g . Although z g and θ g are Sudakov safe observables by themselves, measuring them simultaneously, i.e. z g in a given bin of θ g or vice versa, restores IRC safety.
To define a matching scheme for a Sudakov safe observable, one then rely on the matching of the IRC safe companion. Such matching can be done for instance in a multiplicative way, This formula guarantees that d 2 σ LO+N 2 DL i has exactly the same O(α s ) coefficient as the LO result and reproduces the resumed calculation in the kinematic region enhanced by large logarithms. Notice also that at LO, d 2 σ LO i coincides with the primary Lund plane density. We would like to point out that this matching scheme applies to jets with a given flavour i. As we shall see, this simplifies the resulting formula as our resummed result depends on the jet's flavour via the Casimir factor of the jet. At LO, such a decomposition is trivial, as explained in Sec. 2.3.2. However, beyond LO, this requires to determine in an IRC safe manner the jet's flavour. This can be done using, for instance, the flavour-k t clustering algorithm [52].
Once d 2 σ LO+N 2 DL i is known, the z g or θ g distributions are computed by marginalization. In the case of z g , for example, this amounts to where σ is the inclusive jet cross-section. An important feature of Sudakov safe observables is hidden behind the apparent simplicity of Eq. (2.69). In fact, not all matching schemes for the IRC safe companion lead to a well-defined integral after marginalization. For instance, choosing an additive matching, induces a collinear divergent term (the one inside the bracket) which is not cured by the Sudakov. On the contrary, the integral in Eq. (2.69) is well defined because the Sudakov factor in d 2 σ N 2 DL i shields the θ = 0 logarithmic divergence.
We now turn the concrete implementation of Eq. (2.69) used in this paper. In principle we could compute d 2 σ LO i using MadGraph as in our matched calculation of k t,g . However, we decide here to take another path that we find more enlightening from the physics point of view and, at the same time, easier to implement numerically. Namely, in the small jet radius limit that we consider throughout this paper, it is possible to provide an explicit analytic expression for d 2 σ LO i . Up to powers of θ g corrections, it reads where P i is the symmetrized splitting function of a parton i: P i (z) = P i (z) + P i (1 − z), summed over all decay channels. Matching our resummed distribution to this form of the LO result is then straightforward as it amounts to replace 2C i Θ(e −B i − z)/z in P (z, θ) (Eq. (2.55)) by 2C i P i (z). 13 In Fig. 3, we show a comparison between the exact result of d 2 σ LO i /dz g dθ g obtained through MadGraph and the approximation given by Eq. (2.71). In addition, we compare these two options with the one that we get after replacing the full symmetrized splitting function by its soft limit in Eq. (2.71). More concretely, we have computed d 2 σ LO i for the gluon channel in the three ways that we have just mentioned and show the z g and θ g -projections for two bins of θ g and z g , respectively. We see that Eq. (2.71) matches the exact leading order result while the soft limit of the splitting function is not enough to accurately reproduce the MadGraph output throughout the whole range of z g . Similar conclusions can be drawn by analysing the θ g -projection of d 2 σ LO i . Both for z g and θ g , the deviation of Eq. (2.71) to the exact result remains below 5%. We have deliberately chosen a low p t bin to ensure that using P i as a proxy for d 2 σ LO i is valid in the regime in which the ALICE measurement has been recorded.
We decide to normalize the z g and θ g distributions to the Born level cross-section, σ 0 , in contrast to the k t,g case where the NLO correction to the inclusive jet cross section, σ 1 , was taken into account. For k tg , the reason we included this correction is to account for the C 1 term, but in principle, a LO matched cross-section can safely be normalized by the Born cross-section without spoiling the targeted accuracy. For Sudakov safe cross-sections such as z g and θ g , the question of finding a matching scheme which makes possible the inclusion of the NLO correction to σ in a consistent way with respect to the resummation counterpart is in fact closely related to the C 1 problem that we proceed to tackle.
The C 1 term in Sudakov safe observables. The impossibility to perturbatively expand Sudakov safe observables in powers of α s invalidates the definition of the C 1 term given by Eq. (2.60). Up to now, the question on whether a 'C 1 -like' contribution to the resummation exists for this type of observables has not even been addressed in the literature. In this paper, we would like to outline a new, dedicated matching scheme for Sudakov safe observables. The main difference with respect to the original proposal by the authors in Ref. [30] is to rely on an IRC safe cross-section which is one-dimensional. This IRC safe companion is built from the Sudakov safe distribution with an additional cut on the kinematic variable that is integrated out. More explicitly, for the dynamically groomed z g distribution, one defines the IRC safe cumulative distribution Σ(z g |θ cut ) using the same grooming procedure, but with an additional cut-off on the angle of the splitting, θ g , that is denoted θ cut . Then, our matching formula is: with Σ LO+N 2 DL (z g |θ cut ) defined using multiplicative matching as in Eq. (2.56). 14 Notice that the normalization of Σ LO+N 2 DL (z g ) to σ 0 +σ 1 is ensured. At first sight, this formula seems to depend on the value of θ cut . Yet, it is not the case provided that θ cut is low enough. To understand this, recall that the Sudakov form factor in Σ N 2 DL (z g ) provides a natural cut-off, θ c , on the angular integration of the branching kernel that scales at DLA with θ c exp(−1/ √ᾱ a) see Eq. (2.13) (also in Ref. [25]). Thus, if θ cut is chosen smaller than θ c , we expect that Σ N 2 DL (z g ) = Σ N 2 DL (z g |θ cut ) (2.73) for z g θ a cut . Consequently, far from the resummation region, the dominant term in Eq. (2.72) is Σ LO+N 2 DL (z g |θ cut ) which correctly captures the large z g ∼ 0.5 domain. On the contrary, in the small z g limit (but not smaller than θ a cut ), we obtain: The second term in the previous equation is a correction to our resummed formula, which looks like a C 1 term. It depends on the value of θ cut as a reminiscence of the non IRC safety of the z g distribution. To see that, one notices that up to a constant factor, C 1 (θ cut ) ∝ ln(θ cut ). Since θ cut cannot be larger than θ c , the C 1 correction is actually of order O(α s / √ α s a) = O( α s /a), at least. There are two interesting features in this scaling behavior. First, the appearance of the square root of α s is characteristic of Sudakov safe quantities. Second, we observe how C 1 can become sizeable for a 1. The latter point reminds us that introducing and ad-hoc parameter, θ cut , in the matching scheme comes with some associated difficulties. In short, from the resummation point of view, one would like θ cut to be as small as possible such that Eq. (2.73) holds. However, the smallness of θ cut can lead to a sizeable C 1 correction in Eq. (2.74), thus spoiling the correct asymptotic limit. A clear trade-off exists and the concrete value of θ cut in the proposed matching scheme and its applicability to phenomenological applications deserve further investigation.

Results
Following the reasoning of the k t,g section, we would like to highlight some features of the z g and θ g analytic distributions before moving on to the comparison against Monte-Carlo simulations and ALICE preliminary data. In the left panel of Fig. 4, we quantify the difference between the double-log calculation of the z g distribution and the LO+N 2 DL'. The purpose of this figure is to highlight the deviation of the z g -distribution from the 1/z behavior when higher orders in the resummation are included. Indeed, we have shown in Eq. (2.14), the z g -distribution has a dynamically generated cut-off at z cut ∼ e −a/ᾱ . For z > z cut , it was shown in Ref. [25] that the distribution falls off like the soft limit of the Altarelli-Parisi splitting function. We observe that NDL and N 2 DL contributions such as the running of the strong coupling or the presence of non-global logarithms induce an almost 50% difference with respect to the DLA result. This should be taken into account when interpreting the experimental data specially when searching for modifications in heavy-ion measurements [53,54].
In the right panel of Fig. 4, we asses the impact of the boundary logarithms in the θ g distribution that were discussed in Sec. 2.4.1. As expected, they only matter at large angles and diverge when θ g → 1. Their contribution amounts to a 10−20% and is therefore mandatory to include them if a theory-to-data comparison is aimed.

Phenomenology at LHC energies
The analytic calculations that we have presented above rely mainly on two approximations: the narrow jet limit and the use of the Altarelli-Parisi splitting function in the matching as a proxy for the leading order result in the case of z g and θ g . In order to evaluate the goodness of such approximations, we test our analytic results against parton level simulations in realistic experimental conditions together with the available experimental data. The results are presented for two values of a in the dynamical grooming condition: a = 1 and a = 2. The reason why we do not consider smaller values of a and, in particular, a = 0.1 as done in the ALICE measurement, is because non-perturbative phenomena, beyond the reach of our analytic pQCD calculation, notably affect dynamically groomed observables when a < 1. In addition, we utilise the N 2 DL' prescription on the resummation side. Ratio Figure 4. Left: z g -distribution for a = 1 (red) and a = 2 (blue) at two different accuracies, DLA and LO+N 2 DL' and their ratio. Right: θ g -distribution for a = 1 (red) and a = 2 (blue) with and without boundary logarithms in the LO+N 2 DL' result and their ratio.

Analytics vs. Monte-Carlo parton level
In this section, we compare our analytic calculation for (k t,g , z g , θ g ) to parton level simulations of dijet events with Pythia8.235 [34] and Herwig7.1.2 [55]. For the latter we use both the default angular-ordered shower that we denote 'Herwig7-AO' [56] and the dipole-type shower, 'Herwig7-Dip', based on Ref. [57]. Given that non-perturbative effects are reduced when going to larger p t , we study an experimental setup, that lies within the ATLAS capabilities [49], where the comparison to pQCD calculations are deemed to be cleaner. The centre of mass energy is set to √ s = 13 TeV, jets are clustered with the anti-k ⊥ algorithm with R = 0.4 and re-clustered with Cambridge/Aachen using FastJet3.3.1 [58]. The analysis is performed on those jets that satisfy: 800 < p t < 1000 GeV and |η| < 1.5. For the Monte-Carlo studies, we used the DyG condition given by Eq. (1.1).
We show the comparison between our analytic result and parton level Monte-Carlo simulations with Pythia8 and Herwig7 for the relative transverse momentum of the dynamically groomed splitting in Fig. 5. A crucial point to understand the fixed-order dominated regime, i.e. k t,g 10 −1 , is that in the default setting of both Monte-Carlos, the parton shower starts off a leading-order 2 → 2 matrix element. 15 Therefore, the fixed-order contribution to the k t,g distribution is exactly zero for these event generators. Hence, at large k t,g , an exact agreement between our analytic result, dominated by the exact NLO matrix element, and the Monte-Carlos, where k t,g is exclusively generated by the parton shower, is not expected. Nevertheless, both event generators use at the very least the leading order Altarelli-Parisi splitting functions. As we have discussed in Sec. 2.4.3, the use of the full splitting function in the resummation (or, similarly, in the parton shower) effectively reproduces the fixed-order result in the narrow jet limit. Then, part of the higher-order corrections to the Born-level process are incorporated through the splitting function in the 15 The αs counting might be misleading at this point. Notice that what we refer to as LO in the analytic result is actually a NLO contribution in the sense that it enters at order αs, i.e. we consider p + p → jj at NLO. parton shower. This can explain the nice agreement between the analytic result and the Monte-Carlo curves for k t,g 10 −1 .
On the resummation side, both showers in Herwig are in relatively good agreement and notably differ from Pythia. This is, a priori, rather counterintuitive based on the nature of the three parton showers that we are evaluating. The dipole-style Herwig shower and the Pythia one use a Catani-Seymour like [59] dipole map, transverse momentum ordering and implement a local recoil scheme. In turn, Herwig7-AO evolves through 1 → 2 splittings by means of a generalised angular variable and employs a global recoil scheme. Based on these general arguments, one would expect Herwig7-Dip and Pythia showers to deliver somewhat similar results. The opposite behavior observed in Fig. 5 points out to a more general Pythia-to-Herwig difference rather than to the showers themselves. We identify the scale at which the QCD shower is stopped to be the source of this discrepancy. In fact, in the default setting, Pythia imposes a relatively low infra-red cut-off of 0.5 GeV, while Herwig uses a more conservative scale of ∼ 1 GeV that is common to both showers [60]. Then, more phase-space is available for radiation in the Pythia case and this leads to the differences observed on the low k t,g side in Fig. 5. Thus, we conclude that the small-k t,g part of the differential distribution is sensitive to the way the infra-red is handled and thus to hadronisation. This point will be further emphasized in the following section. Moreover, any higher order term contained in the Monte-Carlo and not present at N 2 DL in the resummation, e.g. energy-momentum conservation, would affect the low k t,g regime.
In what concerns the comparison between MCs and the analytic result, an enhancement at low k t,g values appears. A very similar trend is observed in Fig. 11 of App. B where we evaluate the impact of removing the 1−z in the hardness variable κ (see Eq.(1.1)) for the Monte-Carlo results. We remind the reader that this factor is a sub-leading, nonlogarithmic correction in our analytic calculation at N 2 DL accuracy. However, this mismatch in the κ definition on the analytics and the Monte-Carlos amounts to a ∼ 10% difference on the low k t,g regime and is, therefore, partly responsible for the bump at k t,g ∼ 10 −3 .
The distribution of the opening angle θ g is displayed in Fig. 6. On the Monte-Carlo side, we again observe a strong sensitivity to the momentum scale at which the shower is cut-off. In fact, no significant differences are observed between the angle tagged by Herwig7-AO and Herwig7-Dip, thus indicating a strong dominance of the choice of IRscale. In particular, we have checked that the small angle bump for a = 2 disappears if the infra-red scale is lowered down in Herwig. We have pinned down two other sources for the analytic-to-Monte-Carlo discrepancy. The first one concerns again the 1−z factor in the definition of κ. In Fig. 11 in App. B, we quantify this effect and observe that the small θ g region can be distorted by ∼ 20−40%, depending on the value of a. The enhancement of θ g at large angles in the MCs with respect to the analytic curve could be explained by the O(θ n g ) terms that we have neglected all along our calculation, both on the resummation side and also on the matching procedure where power suppressed terms in the fixed order result were ignored. On the other hand, it is not guaranteed that the branching kernels implemented in the Monte-Carlos recover the exact soft and large angle limit. Then, we conclude that the disagreement between the analytic calculation of the θ g -distribution and the parton shower results can be understood as a result of the choice of the infra-red scale, the finite z corrections in the κ definition and the jet clustering procedure, being the first one the strongest effect.
Finally, the momentum sharing splitting fraction z g is presented in Fig. 7. We clearly observe the presence of the dynamically generated cut-off that separates the fall-off of the distributions from the flattening. The latter starts earlier for a = 2 given that z cut is smaller in this case, see Eq.(2.14). The agreement between the theory calculation and the Monte-Carlos is reasonable in the intermediate regime of 10 −2 < z g < 10 −1 . Outside this interval, the recoil factor in the hardness definition is responsible for both the depletion at large z g in the MC's with respect to the analytic as well as for the excess at small-z, as can be seen in App. B. In addition, the reduced phase space for emissions at infra-red scales in Herwig as compared to Pythia is manifest and further studied in App. D.

Comparison to preliminary ALICE data
In this section, we scrutinise our analytic calculation against preliminary measurements of DyG observables [31,32]. The experimental analysis is performed at √ s = 5.02 TeV on jets clustered with the anti-k ⊥ algorithm with R = 0.4. An important aspect is that only charged jets that satisfy 60 < p ch t < 80 GeV and |η| < 0.5 are considered. In the analytic calculation, no distinction is made between charged and neutral particles, thus a method to translate the charged p t bin of the data into its full counterpart has to be designed. A few possibilities exist to tackle this problem. One could be to identify via Monte-Carlo simulations the transverse momentum bin that, after subtracting the neutral component, yields most of the jets in the 60 < p t < 80 GeV interval. We have carried out this exercise and found that 80% of the jets that fulfil 64.5 < p t < 102.5 GeV, fall in the 60 < p ch t < 80 GeV category. Another option is to absorb this p t -shift from charged to full jets into a non-perturbative factor that also accounts for the effect of hadronization, initial state radiation and multi-parton interactions. This latter is the approach followed in this paper (see also Ref. [8]) and works as follows. We perform the analytic calculation in the same p t -bin as where the experimental measurement is carried out, i.e. 60 < p t < 80 GeV. Then, to construct the non-perturbative factor two samples have to be generated with Monte-Carlo. The first one includes all non-perturbative effects and only charged particles are clustered. Then, the dynamical grooming analysis is performed on the charged jets that satisfy 60 < p ch t < 80 GeV. The second sample is generated as parton level events. Again, we select jets in the same p t -bin as the theoretical calculation, i.e. 60 < p t < 80 GeV, without any charge selection. Then, our non-perturbative factor is defined as the ratio of the charged hadron level sample and the parton level one, both computed in the same transverse momentum bin. These results are plotted in App. D, where further details on the role of the infra-red cutoff in the parton shower are provided. Finally, the theoretical 60 < p ch T < 80 GeV |η| < 0.5, anti-k ⊥ (R = 0.4) DyGa = 2 Figure 8. Comparison between the analytic result obtained in this paper and the preliminary ALICE data [32] of k t,g . results are multiplied by this phenomenological parameter. Then, the theoretical error band includes both the uncertainty of the non-perturbative factor and the analytic uncertainties characterised in the end of Sec. 2.3.3. We label these results as 'LO+N 2 DL'+NP'.
Like in the previous section, we start the discussion with the k t,g distribution shown in Fig. 8. To begin with, an important remark is that a mismatch exists between how the k t,g is defined in the analytic calculation, i.e. k t,g = z g R g , 16 and in ALICE's measurement, where k t,g = z g sin(R g )p t . As we have already mentioned, p t degradation is ignored in our calculation because the hardest branching is located on the primary Lund plane at our degree of accuracy. Then, to accommodate the p t dependence of the experimental definition we simply multiply our analytic result by the lower bound of the p t bin, 60 GeV in this case. We have checked that changing this factor by any other value within the explored p t -bin leads to variations that are well covered by our uncertainty bands. The functional form of the angular dependence of the two k t,g definitions is a bit more delicate. This is so because considering sin(R g ) instead of R g brings additional power-corrections in the calculation that we have so far neglected based on our narrow jet approximation. Besides this fact, we observe in Fig. 8 that the data points are only described by the theoretical calculation if the non-perturbative factor, displayed in Fig. 13, is included. In particular, its role is most prominent for the first bin and generates a large uncertainty. This is yet another manifestation of the different methods that Pythia and Herwig employ to regulate the infra-red sector in the shower. We also provide the Monte-Carlo to data comparison in App. E and find that all three explored setups result into 10−20% deviations with respect to the data both for a = 1 and a = 2. Therefore, we conclude that the agreement between the analytic result presented in this paper and ALICE data is satisfactory in spite of the low p t selection where hadronization effects are very large. Turning to θ g , represented in Fig. 9, we observe that the ad-hoc non-perturbative factor completely dominates the result in the first bin for a = 1. The a-dependence of these results is also interesting from the point of view of missing terms in the analytic calculation. Indeed, we see a deficit in the analytic result for splittings with angles θ g > 0.6 that disappears for the a = 2. We have already mentioned that the g nm coefficients are inversely proportional to a and thus higher-order terms would impact less the a = 2 result than the a = 1 one. In addition, the discrepancy appears in the region where the powersuppressed terms that we have neglected all along the calculation based on the narrow jet approximation, i.e. contributions of O(θ n g ), may matter. An obvious way to confirm this hypothesis would be either to include them or, alternatively, to make a jet radius scan of this observable on the experimental side.
To end up this phenomenological section, we present the theory-to-data comparison for the momentum sharing fraction in Fig. 10. In this scenario, the non-perturbative factor is prominent at small z g , but has a relatively mild effect for z g > 0.2. In fact, in this interval both the LO+N 2 DL' and LO+N 2 DL'+NP results agree with the data within uncertainties. Clearly, this observable together with k t,g are the ones for which our theoretical calculation results provides the best description of the experimental measurements. This is a remarkable result given that the Sudakov nature of this observable complicates its theoretical analysis in many different ways, as we have discussed throughout the paper. The Monte-Carlo results are also consistent with the experimental measurement (see App. E).

Conclusions and outlook
The work presented in this paper follows the current global effort towards a precise theoretical description of jet substructure observables that will help us to deepen our understanding of the space-time evolution of QCD jets, both in vacuum and in heavy-ion collisions. In particular, we have focused on three substructure observables defined on the splitting selected by the dynamical grooming method, that is, the hardest one in the jet tree. These three observables are the momentum sharing fraction z g , the opening angle θ g and the relative transverse momentum k t,g . Out of the three, k t,g is the only infra-red and collinearly safe observable. Then, the definition of logarithmic accuracy for z g and θ g is far from trivial and we extensively discuss a possible approach to tackle the problem that consists in defining the accuracy of the Sudakov safe observable in terms of the cumulative distribution of an IRC safe companion. In this way, we demonstrate that the resummation of dynamically groomed observables does not exponentiate, in general, and that a logarithmic counting at the level of the cumulative distribution is therefore better suited. Further, we present all the necessary ingredients to reach next-to-next-to-double logarithmic accuracy in the narrow jet limit. This includes: (i) the resummation of collinear logarithms arising from the running of the coupling and the hard-collinear correction to the splitting function, (ii) the contribution of non-global logarithms at O(α 2 s ) and (iii) the O(α 2 s ) contribution of boundary logarithms in the case of θ g . Remarkably, neither clustering logarithms nor multiple emissions affect these dynamically groomed distributions. We make use of a matching scheme that naturally includes the C 1 term and allows us to recover the exact leading-order result, computed with MadGraph, for large values of k t,g . On the Sudakov safe cases, we employ the full splitting function as a proxy for the fixed-order result and propose a dedicated matching scheme that depends on an ad-hoc cut-off. All in all, we achieve LO+N 2 DL accuracy in our analytic computation.
The analytic framework is tested against three different parton-level Monte-Carlo simulations at high-p t : Pythia8, angular-ordered, and dipole-style showers in Herwig. In the resummation dominated regime, we find a strong dependence of the Monte-Carlo results on the transverse momentum cut-off at which the parton shower stops. This value is smaller for Pythia (∼ 0.5 GeV) than for Herwig (∼ 1 GeV) by default and, as such, more radiation is allowed in the former case. The analytic result regulates the infra-red singularity through a freezing of the running coupling below 1 GeV and, therefore, allows for emissions with all possible transverse momenta. Due to this fact, it is reasonable that the analytic result is closer to Pythia than to Herwig. Even in the region in which the analytic calculation reduces to the fixed-order contribution, the parton shower is fully responsible of the Monte-Carlo results, given that a leading-order matrix element is implemented by default. Despite this mismatch, an overall good agreement is found for the three jet substructure observables that we attribute to the use of the full splitting function in the parton showers that, as we have stated, generates the correct matrix element in the narrow jet limit.
Our last step is to compare the analytic predictions against the preliminary ALICE data. To do so, we supplement the perturbative results with a non-perturbative factor extracted from Monte-Carlo simulations with Pythia and Herwig that accounts for the use of charged tracks, hadronisation and underlying event. This ingredient is particularly crucial in this experimental setup given that only low-p t jets, i.e. p ch t 200 GeV are measurable by the ALICE detector. In fact, it dominates the theoretical prediction in the lower bins of the k t,g , z g and θ g distributions. A quantitative description of k t,g and z g is achieved up to 5−10% deviations in some bins. In the case of θ g , we find deviations of up to 15% in the moderately large angle region. This is precisely the regime in which we have less confidence in our result considering that we have neglected all power suppressed terms of the type O(θ n g ). The natural extension of this work is to go beyond the small-R limit and develop a numerical routine to account for the resummation of non-global and boundary logarithms. Notice that, as far as we are aware, the latter has yet never been achieved in the literature. On the collinear side of the resummation, we could include the recoil of the hard branch, make use of the NLO splitting function together with higher orders in the running coupling. Further, promoting our leading-order matching to NLO is straightforward from a conceptual point of view in the case of k t,g and would improve the agreement with data/parton level MC simulations both at low and high p t and reduce the theoretical uncertainties. Regarding the Sudakov safe observables, we would like to understand the feasibility of the dedicated matching scheme proposed in this paper and, specifically, quantify its dependence on the ad-hoc cut off. Finally, it would be insightful to compare these analytic results for dynamically groomed observables to the newly developed parton showers that aim at achieving perturbative control beyond leading double logarithmic accuracy and leading color [35,36,61]. From an experimental perspective, these theoretical efforts would highly benefit from a high-p t measurement where non-perturbative corrections are deemed to be milder.
Beyond the possible improvements of the p + p calculation, we would like to discuss two further extensions in terms of collision systems: e + p (and eventually e + A), relevant for the future Electron Ion Collider and heavy-ion collisions. The former, despite the relatively low number of constituents per jet [62], provides a cleaner environment with respect to p + p given that both multi-parton interactions and the underlying event will have a residual effect. On the heavy-ion side, dynamically groomed observables can be used to characterise the properties of an in-medium parton shower. In particular, the θ g distribution can be used to experimentally measure the critical resolution angle of the Quark-Gluon Plasma [63], while deviations at large k t,g from respect to the vacuum baseline could suggest rare, hard scatterings between the propagating parton and the medium.
Cambridge/Aachen. This two-step process has advantages from an experimental point of view. However, from a theoretical perspective we have seen in Sec. 2.4.1 that this two-step process induces boundary logarithms in the calculation. In this Appendix, we would like to investigate at the Monte-Carlo level if the observables are modified when defining the jets Impact of different clustering strategies in k t,g as a function of a for dijet events at parton level in PYTHIA with √ s = 13 TeV for k t,g (top), θ g (center) and z g (bottom).
only with C/A. This is shown in Fig. 12. Interestingly, we observe how the bump at large angles whose origin we have discussed in the main text disappears when clustering with C/A. Nevertheless, the impact of these two jet clustering strategies is mild for all cases.

D Non-perturbative corrections with Pythia and Herwig
In order to compare our analytic predictions with ALICE's experimental data, we add nonperturbative effects through a single parameter extracted from Monte-Carlo simulations. This factor, thoroughly explained in Sec. 3.2, is provided in Fig. 13, where we took the ratio of MCs before and after hadronisation. Besides the default settings of Pythia and Herwig, we show two additional curves in which the parton shower cutoff, denoted as µ parton NP , in Herwig is changed from its default value of 1 GeV to the number used in Pythia where µ parton NP = 0.5 GeV. Note that changing this factor does not necessarily imply a one-to-one correspondence between the two event generators. The value of this factor, like any other hadronisation-related parameter, is tuned to data. Then, one cannot vary it when running the Monte-Carlo at hadron level because its predictive power would be negatively affected. Therefore, we only vary this factor for the parton level result, that is, for the denominator of our non-perturbative factor.
The point of the variation of the parton shower stopping is to demonstrate the sensitivity of the dynamically groomed observables to that scale, and the limitations of this method for incorporating hadronisation corrections into analytic calculations. This is manifest in Fig. 13, where the hadron-to-parton ratio varies from 0.5 to 2.5 for those settings that share the same value of µ parton NP , while it explodes for the default Herwig-AO and Herwig-Dip in the limit of non-perturbative values of (k t,g , z g , θ g ). In the latter case, the reason for the rapid growth of the non-perturbative factor, e.g. in the low k t regime, is rooted in the fact that the parton-level shower does not generate splittings below µ parton NP , while hadronization and underlying event populate this part of the phase-space. In terms of Lund planes, the area covered by the parton-level result and the hadron level one are clearly distinct in Herwig. This effect is less pronounced whenever µ parton NP is low, as in default Pythia. This is explicitly shown in Fig. 14.
As we have already mentioned, there is no preferred value of µ parton NP when running parton level simulations and the large variations encountered in the non-perturbative factor simply indicate that the parton-level results are out of their regime of applicability. Then, we decide to use the average of the Monte-Carlo generators with the same value of µ parton NP as the central value of the non-perturbative factor. The uncertainty band is obtained from the envelope of the five MC settings.
E Monte-Carlo description of (z g , θ g , k t,g ) data In this appendix we compare the three Monte-Carlo settings that we explore through this paper, i.e. Pythia8, Herwig7-AO and Herwig7-Dip, to the preliminary ALICE data. The results are shown in Fig. 15. Notice that through these comparisons we are testing simultaneously the parton shower, i.e. dipole-style or angular-ordered, and the hadronization mechanism, i.e. Lund string or cluster models. In the case of k t,g , no significant differences are observed among all Monte-Carlos. For θ g , Herwig7-Dip provides the best description of the data from small to large angles. All three Monte-Carlo settings are able to capture the data in the intermediate range of this measurement 0.4 < θ g < 0.7 and differences only appear in the tails of Fig. 15, where the hadronization mechanism seems to dominate for θ g < 0.4. Finally, all Monte-Carlos show a significant depletion at 0.2 < z g < 0.3 that is ameliorated for a = 1. Pythia achieves the best theory-to-data ratio, but its not obvious for this observable to disentangle between parton-shower dominated differences and hadronization mechanisms. Monte-Carlo to data comparison of k t,g (top), θ g (middle) and z g (bottom) for a = 1 (left) and a = 2 (right) in the dynamical grooming condition, see Eq.(1.1). In the bottom panels, the theory-to-data ratios, computed using the same binning as the data, are presented.