# Lund multiplicity in QCD jets

Rok Medves, Alba Soto-Ontoso, Gregory Soyez

► To cite this version:

Rok Medves, Alba Soto-Ontoso, Gregory Soyez. Lund multiplicity in QCD jets. Journal of High Energy Physics, 2023, 04, pp.104. 10.1007/JHEP04(2023)104. hal-03921030

HAL Id: hal-03921030

https://hal.science/hal-03921030v1

Submitted on 18 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RECEIVED: December 23, 2022

REVISED: March 10, 2023

ACCEPTED: April 3, 2023

PUBLISHED: April 24, 2023

# Lund multiplicity in QCD jets

Rok Medves,ª Alba Soto-Ontosoᵇ and Gregory Soyezᵇ,ᶜ

E-mail: rok.medves@physics.ox.ac.uk, alba.soto.ontoso@cern.ch, gregory.soyez@ipht.fr

**ABSTRACT:** We compute the average Lund multiplicity of high-energy QCD jets. This extends an earlier calculation, done for event-wide multiplicity in $e^{+}e^{-}$ collisions [1], to the large energy range available at the LHC. Our calculation achieves next-to-next-to-double logarithmic (NNDL) accuracy. Our results are split into a universal collinear piece, common to the $e^{+}e^{-}$ calculation, and a non-universal large-angle contribution. The latter amounts to 10–15% of the total multiplicity. We provide accurate LHC predictions by matching our resummed calculation to fixed-order NLO results and by incorporating non-perturbative corrections via Monte Carlo simulations. Including NNDL terms leads to a 50% reduction of the theoretical uncertainty, with non-perturbative corrections remaining below 5% down to transverse momentum scales of a few GeV. This proves the suitability of Lund multiplicities for robust theory-to-data comparisons at the LHC.

**KEYWORDS:** Jets and Jet Substructure, Resummation

ARXIV ePRINT: 2212.05076

<footnote>ªRudolf Peierls Centre for Theoretical Physics, Clarendon Laboratory, Parks Road, University of Oxford, Oxford OX1 3PU, U.K.</footnote>

<footnote>ᵇUniversité Paris-Saclay, CNRS, CEA, Institut de physique théorique, 91191, Gif-sur-Yvette, France</footnote>

<footnote>ᶜTheoretical Physics Department, CERN, 1211 Geneva 23, Switzerland</footnote>
# Contents

1 Introduction 1

2 Recap of average Lund multiplicity in $e^{+}e^{-}$ collisions 4

3 Calculation of the average Lund multiplicity in jets 9

3.1 DL accuracy 9

3.2 NDL accuracy 10

3.3 NNDL accuracy 13

3.3.1 Hard matrix-element correction (top of the Lund plane) 14

3.3.2 Interplay between the top and bottom of the Lund plane 17

3.3.3 Two large-angle emissions 18

3.3.4 A large-angle emission and a hard-collinear emission 20

3.3.5 A large-angle emission and a running-coupling correction 21

3.4 Final result 21

3.5 From fixed kinematics to fiducial cuts 22

4 Numerical results 23

4.1 Resummation 24

4.2 Phenomenological predictions 25

5 Conclusions 29

A Kinematic weights for dijets 31

B Extraction of the top-of-the-Lund-plane coefficients $D_{hme,R}^{(bc)}$ 32

# 1 Introduction

Jet substructure observables are instrumental in advancing our understanding of Quantum Chromodynamics (QCD) at hadron colliders. Today, infrared-and-collinear (IRC) safe jet substructure observables are successfully being utilised for a variety of experimental applications including quark-gluon jet discrimination [2–5], jet tagging [6–9], and precision measurements [10–22].

Theoretically, the description of jet substructure observables, especially when accounting for multiple splittings, hinges on resummation techniques. These techniques typically organise the perturbative expansion of the observable in terms of single logarithmic $\alpha_s L$ terms, with $L$ a large logarithm of the considered jet substructure observable. In the literature, the accuracy of such calculations is regarded as leading logarithmic (LL), next-to-leading logarithmic (NLL), or generally $N^kLL$. NLL accuracy has now become the standard
in jet substructure calculations for proton-proton collisions and several predictions are available in the literature for both ungroomed [23–31] and groomed [32–38] observables. Quite recently, the first next-to-next-to-leading logarithmic (NNLL) calculations for the SoftDrop jet mass have been put forward [39–41]. This increased level of accuracy has opened the path to extractions of $\alpha_s$ at hadron colliders, potentially competitive with fits to lepton collider data [42].

However, some observables, including the one studied in this paper, do not obviously follow the (N$^k$LL) exponentiating resummation pattern. Instead one relies on an expansion starting from double-logarithms (DL) and resumming terms proportional to $\alpha_s^n L^{2n-k}$ at N$^k$DL accuracy. With respect to this, calculations related to subjet multiplicity have achieved next-to-double logarithmic (NDL) accuracy [1, 43, 44]. This approach has also proven useful for resumming advanced substructure techniques like dynamical grooming observables [45] whose resummation structure does not exponentiate, in particular for Sudakov-safe [46] observables. These calculations are known to next-to-next-to-double logarithmic (NNDL) accuracy. These first-principle calculations of jet substructure observables can then be directly compared to both experimental data [10–22] and/or full-fledged Monte Carlo simulations.<fnref n="1" />

In recent years, the so-called Lund-plane techniques [47] have proven to be a powerful tool for addressing a wide range of jet substructure questions [1, 5, 18, 26, 48]. The most recent observable in this framework has been the *Lund multiplicity* recently introduced in ref. [1]. This new observable counts the number of jets in a hemisphere above a certain momentum cutoff $k_{t,cut}$. In general, this number fluctuates from event to event. We focus here on its average which can be more easily brought to higher-order calculations in perturbation theory. The presence of two disparate scales, i.e. the hard scale $Q$ and the infrared cutoff, leads to the appearance of large logarithms $L \equiv \ln(Q/k_{t,cut})$ that must be resummed to all-orders in order to guarantee the convergence of the perturbative expansion. The resummation structure of this observable can be organised in terms of double logarithms (N$^k$DL) as follows:

$$ \langle N^{(\text{Lund})}(\alpha_s; L) \rangle = \left[ \underbrace{h_1(\alpha_s L^2)}_{\text{DL}} + \underbrace{\sqrt{\alpha_s} h_2(\alpha_s L^2)}_{\text{NDL}} + \underbrace{\alpha_s h_3(\alpha_s L^2)}_{\text{NNDL}} + \dots \right] + \mathcal{O}(e^{-|L|}) , \quad (1.1) $$

with $L \equiv \ln(Q/k_{t,cut})$. Eq. (1.1) is valid in the perturbative limit in which $\alpha_s \ll 1$, $\alpha_s L^2 \sim 1$. Note that this expression differs from the more standard exponentiating (N$^k$LL) form $\ln\langle N \rangle = g_1(\alpha_s L)/\alpha_s + g_2(\alpha_s L) + \alpha_s g_3(\alpha_s L) + \dots$ which instead would correspond to the limit $\alpha_s \ll 1$, $\alpha_s L \sim 1$.<fnref n="2" /> Since in the N$^k$DL limit we have $\alpha_s L \ll 1$, the expansion in eq. (1.1) only requires a finite number of subleading logarithmic corrections at each

<footnote n="1">¹Note that in some phase-space regions, non-perturbative effects such as hadronisation or multi-parton interactions cannot be neglected and the perturbative results have to be complemented with a non-perturbative factor, determined by either Monte Carlo or phenomenological modelling.</footnote>

<footnote n="2">²To the best of our knowledge, there is no full NLL calculation of the average multiplicity. The generating function approach of ref. [49] is close to achieving NLL accuracy in the sense that it resums an arbitrary number of single-logarithmic corrections. However, in order to obtain a closed analytic expression, approximations which invalidate NLL accuracy have to be taken.</footnote>
logarithmic order. This allows one to obtain analytic expressions for the average Lund multiplicity with a well-defined and systematically-improvable logarithmic accuracy.

In ref. [1], $\langle N^{(\text{Lund})} \rangle$ was calculated at NNDL accuracy for $e^+e^-$ collisions, i.e. up to and including $h_3$ in eq. (1.1), and at NDL accuracy for colour singlet production at hadron colliders. In the present manuscript, we extend the calculation to jets produced in high-energy proton-proton collisions, achieving NNDL accuracy. To that end, we first generalise the Lund-multiplicity definition so as to count the number of subjets inside high-energy jets. As any other Lund-plane-based observable, the procedure begins with reclustering a jet of interest at a high-energy collider with the Cambridge-Aachen algorithm [50, 51].<fnref n="3" /> Then, its Lund multiplicity is defined as follows:

1. Set $N^{(\text{Lund})} = 1$.

2. Undo the last clustering step to get two subjets $j_1$ and $j_2$, with $j_1$ the harder branch, i.e. $p_{t1} > p_{t2}$.

3. Calculate the relative transverse momentum of the splitting as

$$k_t = \min(p_{t1}, p_{t2})\Delta_{12} \quad \text{where} \quad \Delta_{12}^2 = (y_1 - y_2)^2 + (\phi_1 - \phi_2)^2, \quad (1.2)$$

with $y_{1,2}$ ($\phi_{1,2}$) the rapidities (azimuthal angles) of the subjets with respect to the beam.

4. If $k_t \ge k_{t,\text{cut}}$ the splitting contributes to the Lund multiplicity, i.e. $N^{(\text{Lund})}$ is increased by one, and we go back to step 2 for each of the two subjets.

5. Otherwise, if $k_t < k_{t,\text{cut}}$, repeat from step 2 following only the harder subjet $j_1$.

The procedure terminates when there is nothing left to decluster. The result of the above algorithm is one plus the total number of emissions in both the primary Lund plane and subsidiary Lund leaves whose transverse momentum is above $k_{t,\text{cut}}$. Naturally, the average Lund multiplicity in a high-energy jet follows the resummation structure of eq. (1.1). The main complication compared to the earlier calculation of the Lund multiplicity in $e^+e^-$ events comes from a richer structure of large-angle gluon radiation. In this context, the treatment of gluon radiation at large angles is known to be prohibitively complex if one wants to retain full-colour accuracy. In this work, we will target full-colour accuracy for the first large-angle emission but limit ourselves to a leading-colour prescription for subsequent large-angle emissions.

The paper is organised as follows. In section 2 we present a brief overview of the NNDL calculation of $\langle N^{(\text{Lund})} \rangle$ in $e^+e^-$ collisions that was extensively discussed in ref. [1]. Then, we move on to the novel ingredients of the pp calculation in section 3 which start appearing at NDL accuracy. We supplement our resummed calculation with a matching to the exact $O(\alpha_s^2)$ result and a study of its sensitivity to non-perturbative corrections in section 4. The final phenomenological predictions for LHC energies can be found in section 4.2. We conclude with a summary of our main findings in section 5.

<footnote n="3">³We note that using a different recombination scheme for the Cambridge-Aachen reclustering procedure, such as the winner-takes-all [52] recombination scheme, only brings differences beyond our target NNDL accuracy (cf. appendix A.2 of ref. [1]).</footnote>
## 2 Recap of average Lund multiplicity in $e^+e^-$ collisions

The purpose of this section is to summarise the working principles of the resummation strategy that we pursue for Lund multiplicity. An extended discussion and full details can be found in ref. [1]. We begin by considering the Lund multiplicity of a $q/g$-initiated hemisphere in an $e^+e^-$ collision at a centre-of-mass energy $Q$. Note that in this case, we use energies and angles instead of $p_\perp$ and $\Delta$, with the relative transverse momentum defined as $k_t = \min(E_1, E_2) \sin\theta_{12}$, which matches the pp definition in the collinear limit. Our aim is to compute the average number of clusterings, within this hemisphere, that satisfy the $k_t > k_{t,\text{cut}}$ condition at NNDL accuracy. To this end, we calculate the $h_1, h_2$ and $h_3$ functions entering into eq. (1.1).

**DL accuracy.** We begin by discussing the case in which a parton with energy $Q$ radiates a set of $n$ soft-and-collinear gluons, that is, the double-logarithmic (DL) limit. In this case, emissions are strongly ordered both in energy, i.e. $Q \gg E_1 \gg E_2 \gg \dots \gg E_n$, and in angle. Further, we can take the coupling constant to be fixed. The first important observation is that, at DL accuracy, the only non-zero contribution comes from the kinematic configuration where the second gluon is emitted from the first gluon (i.e. with $\theta_{21} \ll \theta_{10}$ in the collinear limit), the third gluon is emitted from the second ($\theta_{32} \ll \theta_{21}$), etc... All the other configurations cancel between real and virtual contributions. In terms of colour factors, it means that the only non-zero contribution comes with a colour factor $C_i C_A^{n-1}$ with $C_i = C_F$ ($C_i = C_A$) for a quark-initiated (gluon-initiated) jet. In what follows we will refer to this configuration as a nested chain of double-logarithmic emissions.

To compute the average DL Lund multiplicity, we first write the angular ordering condition as $0 < \eta_1 \ll \eta_2 \ll \dots \ll \eta_n$, where $\eta_j = -\ln\tan(\theta_{j,j-1}/2)$, $\theta_{j,j-1}$ is the angle between the $j$th and $(j-1)$th emission and $\theta_{1,0} \equiv \theta_1$ is the angle between the first emission and the initial hard parton. At DL accuracy, the differential probability for each emission can be written as $\bar{\alpha} \frac{dk_t}{k_t} d\eta = \bar{\alpha} \frac{dx}{x} d\eta$, with $\bar{\alpha} = 2\alpha_s C_A / \pi$, $x$ the energy fraction of the emission and $k_t$ its relative transverse momentum with respect to the emitter. By iterating this branching probability for $n$ emissions, and imposing angular and energy ordering, we obtain the average Lund multiplicity at DL accuracy

$$ 
\begin{aligned}
N_i^{(\text{DL})} &= h_{1,e^+e^-}^{(i)} = 1 + \frac{C_i}{C_A} \sum_{n=1}^{\infty} \bar{\alpha}^n \int_0^{\infty} d\eta_1 \int_{\eta_1}^{\infty} d\eta_2 \dots \int_{\eta_{n-1}}^{\infty} d\eta_n \\
& \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \int_0^1 \frac{dx_1}{x_1} \int_0^{x_1} \frac{dx_2}{x_2} \dots \int_0^{x_{n-1}} \frac{dx_n}{x_n} \Theta(x_n e^{-\eta_n} > e^{-L}) \\
&= 1 + \frac{C_i}{C_A} (\cosh \nu - 1), \quad \text{with } \nu = \sqrt{\bar{\alpha} L^2}, \bar{\alpha} = \frac{2\alpha_s C_A}{\pi}, \qquad (2.1)
\end{aligned}
$$

where we have omitted the angular brackets on the left hand side for conciseness. In this expression, the ‘1’ accounts for the jet initiator. Similarly, the differential distribution of Lund declusterings at a given $k_t$ is given by

$$ n_i^{(\text{DL})} = \frac{dN_i^{(\text{DL})}}{dL} = \frac{C_i}{C_A} \sqrt{\bar{\alpha}} \sinh \nu. \qquad (2.2) $$
<figure><img src="image_2.png" /><figcaption>(a) Hard collinear splitting.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(b) Running coupling.</figcaption></figure>

<caption>**Figure 1.** Lund representation of the NDL configurations that are identical in $e^{+}e^{-}$ collisions and pp jets.</caption>

**NDL accuracy.** A key aspect of our resummation strategy is that in order to achieve NDL accuracy one only has to lift the soft-and-collinear constraint for a single emission in the nested chain of double-logarithmic emissions. This NDL-like emission occurs at a scale $\ell = \ln(Q/k_t)$ and can be either real or virtual and either primary or secondary. Let us denote by $K_{\text{NDL}}$ the genuine NDL corrections to the matrix element for this emission. All other emissions preceding and following this NDL emission down to the infrared cutoff $k_{t,\text{cut}}$ are then treated in the DL, soft-and-collinear, approximation and resummed using eq. (2.1) or (2.2). The master formula that embodies the above ideas and gives the pure NDL correction to the Lund multiplicity reads

$$ \delta N_i^{(\text{NDL})} = \int_0^L d\ell \left\{ K_{\text{NDL}}^R \left[ N_{\text{hard}}^{(\text{DL})}(L;\ell) + N_{\text{soft}}^{(\text{DL})}(L;\ell) \right] - K_{\text{NDL}}^V N_i^{(\text{DL})}(L;\ell) \right\} + \int_0^L d\ell_1 n_i^{(\text{DL})}(\ell_1) \int_{\ell_1}^L d\ell_2 \left\{ K_{\text{NDL}}^R \left[ N_{\text{hard}}^{(\text{DL})}(L;\ell_2) + N_{\text{soft}}^{(\text{DL})}(L;\ell_2) \right] - K_{\text{NDL}}^V N_g^{(\text{DL})}(L;\ell_2) \right\}, \quad (2.3) $$

where $N^{(\text{DL})}(L;\ell) \equiv N^{(\text{DL})}(L-\ell)$ was defined in eq. (2.1). The first line of eq. (2.3) accounts for the case in which the NDL emission takes place along the primary branch, while the second line describes a chain of DL emissions between scales $Q$ and $Qe^{-\ell_1}$, followed by an NDL correction in some secondary branch between scales $Qe^{-\ell_1}$ and $Qe^{-\ell_2}$, followed again by DL emissions between scales $Qe^{-\ell_2}$ and $k_{t,\text{cut}} = Qe^{-L}$. For corrections associated with a real emission (terms involving $K_{\text{NDL}}^R$) one includes subsequent branchings on both the hard and soft branches, as indicated by the subscripts in the squared brackets of eq. (2.3).

In order to systematically identify all NDL corrections we first split them into two categories: corrections related to the running of the strong coupling and kinematic corrections. The former corresponds to 1-loop running coupling corrections, i.e. $\alpha_s \to \alpha_s(k_t) \simeq \alpha_s(1+2\alpha_s\beta_0 \ln Q/k_t)$. A powerful tool for listing all corrections of kinematic origin is the Lund plane representation [47], in which the location of a given emission in the Lund plane reveals at which accuracy it contributes. For example, a DL emission lives in the bulk of the Lund plane, while NDL ones populate regions of the Lund plane where the integrated matrix element is single-logarithmic, i.e. scales as $\alpha_s L$. This is the case when the emission is either hard-collinear, wide-angle or has $k_t \sim k_{t,\text{cut}}$. Of these, ref. [1] showed that only
the hard-collinear yields a NDL correction. The NDL Lund multiplicity, including the running-coupling and hard-collinear corrections, displayed in figure 1, is found to be

$$h_{2,e^+e^-}^{(q)} = h_{2,\text{coll}}^{(q)} = \frac{C_F}{\sqrt{2\pi C_A}} \left\{ \frac{\pi\beta_0}{2C_A} \left[ (\nu^2 - 1) \sinh \nu + \nu \cosh \nu \right] + (B_{gg} + c_\delta B_{gq})\nu \cosh \nu \right. \\ \left. + [2B_q - B_{gg} + (2 - 3c_\delta)B_{gq}] \sinh \nu + 2(c_\delta - 1)B_{gq}\nu \right\}, \qquad (2.4a)$$

$$h_{2,e^+e^-}^{(g)} = h_{2,\text{coll}}^{(g)} = \sqrt{\frac{C_A}{2\pi}} \left\{ \frac{\pi\beta_0}{2C_A} \left[ (\nu^2 - 1) \sinh \nu + \nu \cosh \nu \right] \right. \\ \left. + (B_{gg} + c_\delta B_{gq})\nu \cosh \nu + [B_{gg} + (2 - c_\delta)B_{gq}] \sinh \nu \right\}, \qquad (2.4b)$$

for quark and gluon hemispheres respectively. Here, $c_\delta = (2C_F - C_A)/C_A$, and the values of the hard-collinear $B_{ij}$ coefficients can be found in table 1. The subscript “coll” highlights that these corrections are purely in the collinear region, which will be useful later.

**NNDL accuracy.** The same logic applies when going one order higher in logarithmic accuracy, where only a finite number of emissions should be treated beyond the soft-and-collinear approximation. Again, it is useful to distinguish kinematic corrections from running coupling effects. Let us first focus on the former. An NNDL contribution of the form $\alpha_s^n L^{2n-2}$ can be obtained either through one genuine “NNDL-like” emission contributing a factor $\alpha_s$ together with $(n-1)$ DL, soft-and-collinear, emissions, i.e. $\alpha_s^n L^{2n-2} = (\alpha_s)(\alpha_s L^2)^{n-1}$, or through two NDL-like emissions dressed with $(n-2)$ DL emissions, i.e. $\alpha_s^n L^{2n-2} = (\alpha_s L)^2(\alpha_s L^2)^{n-2}$. NNDL running coupling corrections arise through terms proportional to $\beta_0^2$, generated by either one or two emissions, and mixed cases in which one emission carried a $\beta_0$-term and another pertains to the hard-collinear regime. Note that two-loop running-coupling corrections, proportional to $\beta_1$, would only start appearing at N$^3$DL.

Hence, in order to resum the average Lund multiplicity to NNDL accuracy, one must: (i) systematically identify the relevant contributions, (ii) for NNDL-like contributions, compute the corresponding kernel $K_{\text{NNDL}}$ by means of a fixed-order calculation, and (iii) perform the resummation using eq. (2.3) with $K_{\text{NDL}}$ replaced by $K_{\text{NNDL}}$, or use a similar equation with two NDL insertions for contributions of the form $(\alpha_s L)^2(\alpha_s L^2)^{n-2}$. In total, 8 non-zero corrections contribute to the average Lund multiplicity at NNDL accuracy. The result can be written as

$$h_{3,e^+e^-}^{(q,g)} = h_{3,\text{coll}}^{(q,g)} + h_{3,e^+e^-,\text{la}}^{(q,g)}, \qquad (2.5)$$

with

$$2\pi h_{3,e^+e^-,\text{la}}^{(q)} = D_{\text{hme}}^{qgg} \cosh \nu, \qquad (2.6)$$
<table><caption>**Table 1.** Coefficients entering the NDL $h_{2,e^+e^-}$ (top entries) and NNDL $h_{3,e^+e^-}$ functions (bottom entries) for the average Lund multiplicity in $e^+e^-$ events.</caption><thead><tr><th>Contribution</th><th>Coefficient</th></tr></thead><tbody><tr><td>running coupling</td><td>β<sub>0</sub> = <sup>11C<sub>A</sub> - 4n<sub>f</sub>T<sub>R</sub></sup>&frasl;<sub>12π</sub></td></tr><tr><td>hard-collinear correction</td><td>B<sub>q</sub> = -<sup>3</sup>&frasl;<sub>4</sub>, B<sub>gg</sub> = -<sup>11</sup>&frasl;<sub>12</sub>, B<sub>gq</sub> = <sup>n<sub>f</sub>T<sub>R</sub></sup>&frasl;<sub>3C<sub>A</sub></sub></td></tr><tr><td>hard matrix-element (e<sup>+</sup>e<sup>-</sup>)</td><td>D<sub>hme</sub><sup>qgg</sup> = C<sub>F</sub>(<sup>π<sup>2</sup></sup>&frasl;<sub>6</sub> - <sup>7</sup>&frasl;<sub>4</sub>),<br>D<sub>hme</sub><sup>ggg</sup> = C<sub>A</sub>(<sup>π<sup>2</sup></sup>&frasl;<sub>6</sub> - <sup>49</sup>&frasl;<sub>36</sub>), D<sub>hme</sub><sup>gqq</sup> = n<sub>f</sub>T<sub>R</sub><sup>2</sup><sup>2</sup></td></tr><tr><td>collinear endpoint</td><td>D<sub>end</sub><sup>q→gg</sup> = C<sub>F</sub>(3 + 3 ln 2 - <sup>π<sup>2</sup></sup>&frasl;<sub>3</sub>)<br>D<sub>end</sub><sup>g→gg</sup> = C<sub>A</sub>(<sup>137</sup>&frasl;<sub>36</sub> + <sup>11</sup>&frasl;<sub>3</sub> ln 2 - <sup>π<sup>2</sup></sup>&frasl;<sub>3</sub>)<br>D<sub>end</sub><sup>g→qq</sup> = n<sub>f</sub>T<sub>R</sub>(-<sup>29</sup>&frasl;<sub>18</sub> - <sup>4</sup>&frasl;<sub>3</sub> ln 2)</td></tr><tr><td>commensurate k<sub>t</sub> and angle</td><td>D<sub>pair</sub><sup>qq</sup> = <sup>13</sup>&frasl;<sub>9</sub>n<sub>f</sub>T<sub>R</sub>, D<sub>pair</sub><sup>gg</sup> = (<sup>π<sup>2</sup></sup>&frasl;<sub>6</sub> - <sup>67</sup>&frasl;<sub>18</sub>)C<sub>A</sub><br>K = (<sup>67</sup>&frasl;<sub>18</sub> - <sup>π<sup>2</sup></sup>&frasl;<sub>6</sub>)C<sub>A</sub> - <sup>10</sup>&frasl;<sub>9</sub>n<sub>f</sub>T<sub>R</sub></td></tr><tr><td>clustering</td><td>D<sub>clust</sub><sup>(prim)</sup> = -<sup>5π<sup>2</sup></sup>&frasl;<sub>54</sub>, D<sub>clust</sub><sup>(sec)</sup> = <sup>π<sup>2</sup></sup>&frasl;<sub>27</sub></td></tr><tr><td>energy loss</td><td>D<sub>e-loss</sub><sup>q</sup> = <sup>7</sup>&frasl;<sub>2</sub>C<sub>A</sub> + (<sup>5</sup>&frasl;<sub>2</sub> - <sup>2π<sup>2</sup></sup>&frasl;<sub>3</sub>)C<sub>F</sub><br>D<sub>e-loss</sub><sup>g</sup> = (<sup>67</sup>&frasl;<sub>9</sub> - <sup>2π<sup>2</sup></sup>&frasl;<sub>3</sub>)C<sub>A</sub> - <sup>26</sup>&frasl;<sub>9</sub>C<sub>A</sub><sup>2</sup>n<sub>f</sub>T<sub>R</sub></td></tr></tbody></table>

$$2\pi h_{3,\text{coll}}^{(q)} = D_{\text{end}}^{q \to gg} + (D_{\text{end}}^{g \to gg} + D_{\text{end}}^{g \to q\bar{q}}) \frac{C_F}{C_A} (\cosh \nu - 1) \quad (2.7) \\ + \frac{C_F}{C_A} \left[ (1 - c_\delta) D_{\text{pair}}^{q\bar{q}} (\cosh \nu - 1) + \left(K + D_{\text{pair}}^{gg} + c_\delta D_{\text{pair}}^{q\bar{q}}\right) \frac{\nu}{2} \sinh \nu \right] \\ + C_F \left[ \left(\cosh \nu - 1 - \frac{1 - c_\delta}{4} \nu^2\right) D_{\text{clust}}^{(\text{prim})} + \left(\cosh \nu - 1\right) D_{\text{clust}}^{(\text{sec})} \right] \\ + \frac{C_F}{C_A} \left[ D_{\text{e-loss}}^g \frac{\nu}{2} \sinh \nu + \left(D_{\text{e-loss}}^q - D_{\text{e-loss}}^g\right) (\cosh \nu - 1) \right] \\ + \frac{C_F}{2} \left\{ (B_{gg} + c_\delta B_{gq})^2 \nu^2 \cosh \nu + 8 \left[ 2c_\delta B_{gg} - 2c_\delta B_q - (1 - 3c_\delta^2) B_{gq} \right] B_{gq} \cosh \nu \right. \\ \left. + [4B_q(B_{gg} + (2c_\delta + 1)B_{gq}) - (B_{gg} + c_\delta B_{gq})(B_{gg} + 9c_\delta B_{gq})] \nu \sinh \nu \right. \\ \left. + 4(1 - c_\delta^2) B_{gq}^2 \nu^2 + 8 \left[ 2c_\delta B_q - 2c_\delta B_{gg} + (1 - 3c_\delta^2) B_{gq} \right] B_{gq} \right\} \\ + \frac{C_F}{C_A} \frac{\pi \beta_0}{2} \left\{ (B_{gg} + c_\delta B_{gq}) \nu^3 \sinh \nu + [2B_q - 2B_{gg} + (6 - 8c_\delta)B_{gq}] \nu \sinh \nu \right. \\ \left. + 2(B_q + B_{gg} + B_{gq}) \nu^2 \cosh \nu - 4(1 - c_\delta)B_{gq}(2\cosh \nu - 2 + \nu^2) \right\} \\ \left. + \frac{C_F}{C_A} \frac{\pi^2 \beta_0^2}{8C_A} \left[ 3\nu(2\nu^2 - 1) \sinh \nu + (\nu^4 + 3\nu^2) \cosh \nu \right], \right.$$

for quark-initiated jets, and

$$2\pi h_{3,e^+e^-, \text{la}}^{(g)} = D_{\text{hme}}^{ggg} \cosh \nu + D_{\text{hme}}^{gq\bar{q}} (c_\delta \cosh \nu + 1 - c_\delta), \quad (2.8)$$
<figure><img src="image_4.png" /><figcaption>Figure 2. Lund representation of the NNDL configurations that are identical in $e^{+}e^{-}$ and pp collisions. Black dots indicate emissions which contribute a factor $\alpha_s L^2$, blue dots emissions which contribute a factor $\alpha_s L$, and red dots emissions which contribute as $\alpha_s$. Running coupling corrections are represented by open circles.</figcaption></figure>

$$2\pi h_{3,\text{coll}}^{(g)} = (D_{\text{end}}^{g \to gg} + D_{\text{end}}^{g \to q\bar{q}}) \cosh \nu \qquad (2.9) \\ + \left[ (1-c_{\delta})D_{\text{pair}}^{q\bar{q}}(\cosh \nu - 1) + (K + D_{\text{pair}}^{gg} + c_{\delta}D_{\text{pair}}^{q\bar{q}}) \frac{\nu}{2} \sinh \nu \right] \\ + C_A (D_{\text{clust}}^{(\text{prim})} + D_{\text{clust}}^{(\text{sec})}) (\cosh \nu - 1) + D_{\text{e-loss}}^{g} \frac{\nu}{2} \sinh \nu \\ + \frac{C_A}{2} \left\{ (B_{gg} + c_{\delta}B_{gq})^2 \nu^2 \cosh \nu - 8(1-c_{\delta}^2)B_{gq}^2(\cosh \nu - 1) \right. \\ \left. + [(B_{gg} + c_{\delta}B_{gq})(3B_{gg} - 5c_{\delta}B_{gq}) + 4(1+c_{\delta})B_{gq}B_{q}] \nu \sinh \nu \right\} \\ + \frac{\pi\beta_0}{2} \left\{ (B_{gg} + c_{\delta}B_{gq})\nu^3 \sinh \nu + 6(1-c_{\delta})B_{gq}\nu \sinh \nu \right. \\ \left. + 2[2B_{gg} + (1+c_{\delta})B_{gq}]\nu^2 \cosh \nu - 8B_{gq}(1-c_{\delta})(\cosh \nu - 1) \right\} \\ + \frac{\pi^2\beta_0^2}{8C_A} [3\nu(2\nu^2 - 1)\sinh \nu + (\nu^4 + 3\nu^2)\cosh \nu],$$

for gluon-initiated jets.

The values of the coefficients entering in eq. (2.5) can be found in table 1. Note that the results in the above expressions have been explicitly separated in two parts: (i) $h_{3,\text{coll}}^{(q,g)}$ that involves the NNDL contributions sensitive to collinear physics and running-coupling corrections, see figure 2, and (ii) $h_{3,\text{la}}^{(q,g)}$ that includes all terms involving large-angle kinematics. As we will argue in the next section, the collinear contributions can be directly
recycled for the case of the Lund multiplicity in a high-energy jet, while the large-angle contributions have to be recomputed.

## 3 Calculation of the average Lund multiplicity in jets

In this section, we derive the main result of this paper, namely we extend the calculation of the Lund multiplicity to the case of the subjet multiplicity within a high-energy jet.

We therefore start with a jet of transverse momentum $p_{\perp}$, and rapidity $y_{jet}$, assuming that it has been reconstructed using the anti-$k_t$ jet algorithm [53] with a jet radius $R$. The average Lund multiplicity $\langle N_{jet}^{(Lund)}(\alpha_s; L) \rangle$ above a relative transverse momentum cut, $k_{t,cut}$, follows the logarithmic expansion of eq. (1.1) with the logarithm $L$ defined as $L \equiv L_{jet} = \ln(p_{\perp}R/k_{t,cut})$. In practical applications, the transverse momentum and rapidity of the jet will not be fixed but integrated over a fiducial phase space. We discuss how to translate our results for a fixed jet $p_{\perp}$ and $y_{jet}$ to a realistic scenario in section 3.5.

A large fraction of the calculation can be simply recycled from the $e^{+}e^{-}$ results presented in section 2 with $L = L_{jet}$. This includes all the contributions involving only collinear physics (including running-coupling corrections), as a consequence of the universality of the collinear limit of QCD. We are therefore left with the computation of the contributions depending on emissions at large angles (commensurate with the jet radius). In general, these corrections depend on the overall structure of the hard process. In practice, we will consider jets in either $Z$+jet or inclusive jet samples. However, we keep the discussion below as generic as possible so as to facilitate possible future extensions to other processes.

The $e^{+}e^{-}$ results in section 2 have been presented separately for quark and gluon hemispheres. In a pp context, jets can also be considered as quark- or gluon-initiated (as given by the Born-level hard process) as long as we are only interested in the contributions involving either soft gluon emissions and/or collinear branchings. More concretely, these include the collinear contributions recycled from our $e^{+}e^{-}$ results, as well as new large-angle contributions which only depend on soft-gluon emissions. For these contributions we can consider that we work with a jet of a given flavour, and ultimately get the total multiplicity by weighting the quark- and gluon-initiated jet multiplicities by the corresponding Born-level fractions. Up to our target NNDL accuracy, the only contribution where the quark/gluon separation is not obvious is the case of hard matrix-element corrections, sensitive to emissions which are both hard and large-angle. We address this case in section 3.3.1.

As with the $e^{+}e^{-}$ recap of section 2, we start with the DL expressions and then proceed with the calculation of the NDL and NNDL corrections.

### 3.1 DL accuracy

Since the DL Lund multiplicity only involves soft-and-collinear emissions, the results for $h_1^{(q,g)}$ for quark and gluon jets can be taken directly from the corresponding results in $e^{+}e^{-}$ events, i.e. eq. (2.1). The DL average Lund multiplicity for a given hard process is thus given by

$$h_1(\xi) = f_q h_1^{(q)}(\xi) + f_g h_1^{(g)}(\xi), \qquad (3.1)$$
<figure><img src="image_5.png" /><figcaption>Figure 3. Lund representation of the only new NDL contribution for high-energy jets: a large-angle emission.</figcaption></figure>

with

$$f_i = \frac{\sigma_i}{\sigma_q + \sigma_g}, \quad \xi = \alpha_s L^2 \quad \text{and} \quad h_1^{(q,g)}(\xi) = h_{1,e^+e^-}^{(q,g)}(\xi), \qquad (3.2)$$

where the quark/gluon fractions $f_{q,g}$ are computed from the Born-level cross-sections $\sigma_{q,g}$.

## 3.2 NDL accuracy

NDL corrections arise from a series of branchings where a single emission gives a single-logarithmic contribution (proportional to $\alpha_s L$) and all other emissions are soft-and-collinear (proportional to $\alpha_s L^2$). The two contributions which appear in the $e^+e^-$ case, namely the hard-collinear and running coupling correction, cf. figure 1, are unchanged in the case of the Lund multiplicity in a jet. The only potential new source of NDL corrections is the case of an emission at large angles, close to the jet boundary as depicted in figure 3. Including the proper relative weight of quark- and gluon-initiated jets, we therefore write

$$h_2(\xi) = f_q h_{2,\text{coll}}^{(q)}(\xi) + f_g h_{2,\text{coll}}^{(g)}(\xi) + h_{2,\text{la}}(\xi). \qquad (3.3)$$

While $h_{2,la}$ vanishes in $e^+e^-$ collisions, it is no longer the case for high-energy jets and we compute this contribution below.

Generically, the radiation of a large-angle gluon depends on the overall structure of the Born event. Since the gluon is soft, it is however still possible to separate the contribution from quark- and gluon- initiated jets based on the underlying Born hard process. We thus write $h_{2,la} = f_q h_{2,la}^{(q)} + f_g h_{2,la}^{(g)}$. These contributions would nevertheless be hard-process-dependent, i.e. a “gluon” jet in a Z+jet event would be different from a “gluon” jet in an inclusive jet sample. We proceed by writing the Born hard process as a sum over different partonic flavour channels:

$$\sigma^{(i)} = \sum_{C \in \text{flavour channels}} \sigma_C^{(i)}, \qquad (3.4)$$

where the (optional) flavour superscript, i, limits the sum to channels contributing to a specific jet flavour.

We start by considering the $O(\alpha_s)$ contribution to the Lund multiplicity where a single soft gluon is emitted. For a given channel, the matrix element describing such an emission
can be further decomposed into a sum over coloured dipoles (ab) with a and b two (incoming or outgoing) partons in the Born process:

$$ \sum_{a,b \in \text{legs}} \frac{\alpha_s}{2\pi} \omega_{ab}^C \int k_{\perp} \, dk_{\perp} \int dy \int \frac{d\phi}{2\pi} (k|ab), \quad \text{with} \quad (k|ab) = \frac{(p_a \cdot p_b)}{(p_a \cdot k)(k \cdot p_b)}, \qquad (3.5) $$

and $\omega_{ab}^C = (-2\mathbf{T}_a \cdot \mathbf{T}_b)$ is a colour factor or kinematic weight which, in general, depends on the kinematics of the full hard event. Therefore, the NDL correction to the Lund multiplicity associated with the radiation of a soft large-angle gluon takes the form

$$ \delta N_{\text{la}}^{(\text{NDL},i)} = \sqrt{\alpha_s} h_{2,\text{la}}^{(i)}(\xi) = \sum_{C \in \text{flavour channels}} \frac{\sigma_C^{(i)}}{\sigma^{(i)}} \sum_{(ab) \in \text{event}} \delta N_{C,(ab)}^{(\text{NDL})}, \qquad (3.6) $$

where $\delta N_{C,(ab)}^{(\text{NDL})}$ is the correction to the multiplicity due to a soft gluon emitted close to the jet boundary from a dipole (ab). At the first non-trivial order in $\alpha_s$, it is given by

$$ \delta N_{C,(ab),\mathcal{O}(\alpha_s)}^{(\text{NDL})} = \frac{\alpha_s}{\pi} \omega_{ab}^C \int_0^{p_{\perp}} \frac{dk_{\perp}}{k_{\perp}} dy \frac{d\phi}{2\pi} \left[ \frac{k_{\perp}^2}{2} (k|ab) - \frac{1}{\Delta_k^2} \delta_a \text{ or } b=\text{jet} \right] \Theta(k_t > k_{t,\text{cut}}) \Theta(\Delta_k < R). \qquad (3.7) $$

In this expression, the second term in the square bracket subtracts the collinear (DL) contribution and the “$\delta$” selects the cases where one of the two dipole legs corresponds to the jet under consideration. This is done so as to isolate the NDL correction. Also, $\Delta_k^2 = (y - y_{\text{jet}})^2 + (\phi - \phi_{\text{jet}})^2$ is the (squared) distance between the jet axis and the emission. To evaluate (3.7), we first make a change of variable from $k_{\perp}$ to $k_t = k_{\perp} \Delta_k$. Including the constraint that $k_t$ has to be above $k_{t,\text{cut}}$ we can then write

$$ \int_{k_{t,\text{cut}}}^{p_{\perp} \Delta_k} \frac{dk_t}{k_t} = L + \mathcal{O}(\text{NNDL}), \qquad (3.8) $$

where only the logarithmically-enhanced term contributes at NDL. This is allowed because the angular integration does not contain any collinearly-enhanced contribution so that $\Delta_k$ is of order $R$ up to a constant factor which would not generate a logarithmic enhancement. We are then left with the integration over the geometric variables $y$ and $\phi$ so that the NDL correction takes the form

$$ \delta N_{C,(ab),\mathcal{O}(\alpha_s)}^{(\text{NDL})} = \frac{\alpha_s}{\pi} L \omega_{ab}^C D_{ab}^{\text{la}}. \qquad (3.9) $$

We now compute the coefficients $D_{ab}^{\text{la}}$ for a generic $2 \to 2$ process. Let us denote the four-momenta of the incoming partons as $p_1$ and $p_2$, and those of the outgoing jet (for which one measures the multiplicity) and recoiling parton as $p_j$ and $p_r$.<fnref n="4" /> Without loss of generality, we can parametrise the momenta of the hard legs and the radiated gluon, $k$, as

<footnote n="4">⁴In the case of $Z+\text{jet}$ events the recoiling object is the colour singlet $Z$ so that $p_r$ does not contribute to the coloured dipoles.</footnote>
follows (using a $(p_x, p_y, p_z, E)$ convention):

$$p_1 \equiv p_{\perp} \cosh \hat{y} (0, 0, 1, 1), \qquad (3.10a)$$

$$p_2 \equiv p_{\perp} \cosh \hat{y} (0, 0, -1, 1), \qquad (3.10b)$$

$$p_j \equiv p_{\perp} (1, 0, \sinh \hat{y}, \cosh \hat{y}), \qquad (3.10c)$$

$$p_r \equiv p_{\perp} (-1, 0, -\sinh \hat{y}, \cosh \hat{y}), \qquad (3.10d)$$

$$k \equiv k_{\perp} (\cos \phi, \sin \phi, \sinh y, \cosh y), \qquad (3.10e)$$

where $p_{\perp}$ is the transverse momentum of the jet, $\hat{y}$ is the rapidity of the jet in the partonic centre-of-mass frame,<fnref n="5" /> $k_{\perp} \ll p_{\perp}$, and we have picked the jet azimuthal angle to be zero. The soft gluon can be emitted from any of the 6 possible colour dipoles: $(p_1, p_2)$, $(p_1, p_j)$, $(p_2, p_j)$, $(p_1, p_r)$, $(p_2, p_r)$ and $(p_j, p_r)$. Except for the (12) dipole, the integrals in eq. (3.7) cannot be computed exactly analytically. We thus make a series expansion in the jet radius $R$ and find [26]

$$D_{12}^{\text{la}} = \frac{R^2}{2}, \qquad (3.11a)$$

$$D_{1j}^{\text{la}} = \frac{R^2}{8} + \frac{R^4}{576} + \mathcal{O}(R^8), \qquad (3.11b)$$

$$D_{jr}^{\text{la}} = \frac{\tanh^2 \hat{y}}{8} R^2 + \frac{(\cosh^2 \hat{y} - 3)^2}{576 \cosh^4 \hat{y}} R^4 + \frac{\tanh^2 \hat{y}}{384} R^6 + \mathcal{O}(R^8), \qquad (3.11c)$$

$$D_{1r}^{\text{la}} = \frac{e^{2\hat{y}}}{8 \cosh^2 \hat{y}} R^2 + \frac{1}{64 \cosh^4 \hat{y}} R^4 + \frac{\tanh^2 \hat{y}}{384 \cosh^4 \hat{y}} R^6 + \mathcal{O}(R^8). \qquad (3.11d)$$

Then, $D_{2j}^{\text{la}} = D_{1j}^{\text{la}}$ and $D_{2r}^{\text{la}}$ is obtained from $D_{1r}^{\text{la}}$ taking $\hat{y} \to -\hat{y}$. One notices that the (1r) dipole tends to the (12) dipole when $\hat{y} \to \infty$ and to 0 when $\hat{y} \to -\infty$. Similarly, the (jr) dipole is symmetric under $\hat{y} \to -\hat{y}$ and converges to the (1j) dipole at large $\hat{y}$. Dipoles involving the recoiling parton depend on the kinematics of the $2 \to 2$ process via $\hat{y}$. This will have to be taken into account when considering jets within a rapidity range (see section 3.5). In practice, we have checked numerically that the small-$R$ expansion used in eq. (3.11) reproduces the exact result to well within 1% up to $R=1$, and up to deviations of at most 0.3% for $R=0.8$ which we use as a default below.

At all orders, one should include the extra contribution to the multiplicity due to subsequent soft-and-collinear radiation from the large-angle gluon. Using eq. (2.3), we find that $(\nu = \sqrt{\bar{\alpha}L^2})$

$$\delta N_{C,(ab)}^{\text{(NDL)}} = \frac{\alpha_s}{\pi} \omega_{ab}^C D_{ab}^{\text{la}} \int_0^L d\ell N_g^{\text{(DL)}}(L-\ell) = \sqrt{\frac{\alpha_s}{2C_A\pi}} \omega_{ab}^C D_{ab}^{\text{la}} \sinh \nu. \qquad (3.12)$$

To compute the NDL multiplicity for a given hard process, one should also specify its corresponding weights $\omega_{ab}^C$. In the case of $Z$+jet events there are only two partonic channels, one corresponding to quark jets and the other to gluon jets. Since there are only three coloured partons (two incoming and one outgoing), the colour algebra for each channel is

<footnote n="5">⁵In the lab frame, $y_{jet} = \hat{y} + \frac{1}{2} \ln \frac{x_1}{x_2}$ with $x_{1,2}$ the longitudinal momentum fractions of the incoming partons.</footnote>
<table><caption>**Table 2.** Colour coefficients $\omega_{ab}^C \equiv (-2T_a \cdot T_b)$ for different flavour channels in Z+jet events.</caption><thead><tr><th>Channel</th><th>ω<sub>12</sub></th><th>ω<sub>1j</sub></th><th>ω<sub>2j</sub></th></tr></thead><tbody><tr><td>q<span style="text-decoration: overline;">q</span> → Zg</td><td>2C<sub>F</sub> − C<sub>A</sub></td><td>C<sub>A</sub></td><td>C<sub>A</sub></td></tr><tr><td>qg → Zq</td><td>C<sub>A</sub></td><td>C<sub>A</sub></td><td>2C<sub>F</sub> − C<sub>A</sub></td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Lund representation of the new configurations with respect to $e^+e^-$ collisions that contribute to the multiplicity at NNDL.</figcaption></figure>

straightforward and the $\omega_{ab}^C$ weights are listed in table 2. The situation is more involved for dijet events due to both the increased number of partonic channels and the more intricate colour structure of the four-leg process. The corresponding $\omega_{ab}^C$ weights can be extracted from ref. [54] and are given in appendix A for completeness.

### 3.3 NNDL accuracy

At NNDL accuracy, we again have to include two broad families of contributions: collinear contributions which can be obtained from the corresponding calculation in $e^+e^-$, namely $h_{3,\text{coll}}$ given by eq. (2.5), as well as extra corrections involving large-angle emissions. The latter are depicted in figure 4 and we briefly discuss them here before proceeding with their explicit calculation in the following sections.

NNDL corrections can be obtained in two ways (see the discussion in section 2): either from one “NNDL-like” emission contributing a factor $\alpha_s$ (with no logarithmic enhancement) accompanied by any number of soft-and-collinear emissions, or from two “NDL-like” emissions, each contributing a factor $\alpha_s L$, plus any number of soft-and-collinear emissions. In
the first case, contributions at large angles can arise from two kinematic configurations:

* an emission at the top of the Lund plane, figure 4a, corresponding to a hard matrix-element correction, discussed in section 3.3.1,

* or an emission at the bottom-left corner of the Lund plane, figure 4b, corresponding to an emission at large angle and close to the $k_{t,\text{cut}}$ kinematic limit. We have structured the calculation so that this contribution is zero, as discussed in section 3.3.2.

In the second case, the first NDL-like emission must happen at large angle, and the second one can correspond to any NDL correction. The second emission can therefore be

* a collinear emission close the kinematic $k_{t,\text{cut}}$, figure 4c, which also yields a vanishing contribution as discussed in section 3.3.2,

* another large-angle emission, figure 4d, discussed in section 3.3.3,

* a hard-collinear branching, figure 4e, discussed in section 3.3.4,

* or a running-coupling correction, figure 4f, discussed in section 3.3.5.

Altogether, the NNDL function $h_3(\xi)$ for the Lund multiplicity in high-energy jets is written as<fnref n="6" />

$$h_3(\xi) = f_q h_{3,\text{coll}}^{(g)}(\xi) + f_g h_{3,\text{coll}}^{(q)}(\xi) + h_{3,\text{hme}}(\xi) + h_{3,\text{la}^2}(\xi) + h_{3,\text{la}\times\text{hc}}(\xi) + h_{3,\text{la}\times\text{rc}}(\xi), \quad (3.13)$$

where the collinear $h_{3,\text{coll}}^{(i)}(\xi)$ were defined in eqs. (2.7), (2.9), and we discuss the calculation of the final four pieces in the following sections.

### 3.3.1 Hard matrix-element correction (top of the Lund plane)

We first consider the NNDL correction to the Lund multiplicity of a hard, wide-angle emission that arises due to the fact that the soft eikonal approximation is no longer valid and the exact matrix element has to be used instead. Our starting point is the exact expression for the Lund multiplicity at NLO, i.e. $O(\alpha_s)$ relative to the Born-level process

$$\langle N^{(\text{Lund})} \rangle_{O(\alpha_s)} = \frac{1}{\sigma_0 + \sigma_1} \left( \sigma_0 + \int d\Phi |\mathcal{M}_R|^2 [1 + \Theta(k_t > k_{t,\text{cut}})] + \int d\Phi |\mathcal{M}_V|^2 \right), \quad (3.14)$$

where $\sigma_0$ and $\sigma_1$ are respectively the Born-level cross section and its NLO correction, and $|\mathcal{M}_{R,(V)}|$ is the real (virtual) NLO matrix-element for the process being studied. For the real contribution, we have assumed that the two partons are inside the jet under consideration. The generalisation to the case where only one parton is inside the jet is straightforward and discussed below.

The extension of eq. (3.14) to all orders requires further dressing of each term with towers of DL emissions. Since this depends on the partonic flavour, one has to separate each term in (3.14) into flavour channels. Let us therefore assume that one can associate a suitably-defined flavour index $a$ with the virtual contribution, and two flavour indices $b$ and

<footnote n="6">⁶In this expression, the quark and gluon fractions are taken from the Born-level process as in the DL and NDL sections. At NNDL accuracy, NLO corrections to the Born-level cross-sections have to be taken into account and we include them in $h_{3,\text{hme}}(\xi)$ as detailed in section 3.3.1.</footnote>
c with the two partons in the real contribution, corresponding to an $a \to bc$ splitting. For the sake of the argument, let us further assume that we can separate the real and virtual contributions into a “hard” and a “soft” piece:

$$ |\mathcal{M}_R^{(bc)}|^2 = |\mathcal{M}_R^{(bc)}|_{\text{hard}}^2 + |\mathcal{M}_R^{(bc)}|_{\text{soft}}^2, \qquad (3.15) $$

$$ |\mathcal{M}_V^{(a)}|^2 = |\mathcal{M}_V^{(a)}|_{\text{hard}}^2 + |\mathcal{M}_V^{(a)}|_{\text{soft}}^2. \qquad (3.16) $$

For the real contribution, one can define $|\mathcal{M}_R^{(bc)}|_{\text{soft}}^2$ so that it contains the soft and/or collinear divergences which reproduce the contributions to the multiplicity not associated with hard-wide-angle physics. This includes the DL and NDL contributions, the NNDL collinear-endpoint correction (which can be recycled from the $e^+e^-$ result and corresponds to the terms involving the $D_{\text{end}}^{a \to bc}$ coefficients in eqs. (2.7) and (2.9)) and potential NNDL soft-large-angle corrections (discussed in section 3.3.2). This leaves only a non-zero $|\mathcal{M}_R^{(bc)}|_{\text{hard}}^2$ contribution at the top of the Lund plane. Such a separation is less obvious for the virtual correction, but, due to infrared safety, one can take $|\mathcal{M}_V^{(a)}|_{\text{soft}}^2 = -|\mathcal{M}_R^{(bc)}|_{\text{soft}}^2$ for an $a \to bc$ splitting. In practice, the NNDL hard matrix-element contribution to the multiplicity only requires the integrated cross-sections $\sigma_{R(V),\text{hard}(\text{soft})} = \int d\Phi |\mathcal{M}_{R(V)}|_{\text{hard}(\text{soft})}^2$.

The all-orders average multiplicity including all hard matrix-element effects can then be written as

$$ \langle N^{(\text{Lund})} \rangle = \frac{1}{\sigma_0 + \sigma_1} \left\{ \sum_a \sigma_0^{(a)} N_a(L) + \sum_{bc} \sigma_{R,\text{hard}}^{(bc)} [N_b(L) + N_c(L)] \right. \\ \left. + \sum_{bc} \sigma_{R,\text{soft}}^{(bc)} \delta_{a=b+c} N_a(L) + \sum_a \sigma_V^{(a)} N_a(L) \right\}. \qquad (3.17) $$

In writing this expression, we have taken into account the fact that for the DL dressing of the real “soft” contribution, one has to define the flavour sum $a = b + c$ as a proxy for the flavour of the parent parton corresponding to the (bc) flavour pair. The exact procedure for the reconstruction the flavour of a from the daughter partons is irrelevant as long as it is correct in the soft limit.

Adding and subtracting $\sigma_{R,\text{hard}}^{(bc)} \delta_{a=b+c} N_a(L)$ to eq. (3.17) we obtain

$$ \langle N^{(\text{Lund})} \rangle = \frac{1}{\sigma_0 + \sigma_1} \left\{ \sum_a \sigma_0^{(a)} N_a(L) + \sum_{bc} \sigma_{R,\text{hard}}^{(bc)} [N_b(L) + N_c(L) - N_a(L)] \right. \\ \left. + \sum_a \left[ \sum_{bc} \sigma_R^{(bc)} \delta_{a=b+c} + \sigma_V^{(a)} \right] N_a(L) \right\}, \qquad (3.18) $$

where each term can be determined by means of fixed-order Monte Carlo simulations which provide access to the flavours of the partons. The first term is obtained from the quark and gluon contributions to the Born-level cross section. As we have already mentioned, the second term is calculated from the real NLO correction to the cross-section. To isolate the “hard” component, we run with a small $k_{t,\text{cut}}$ and subtract the $O(\alpha_s)$ contributions not associated with hard-wide-angle physics, i.e. the DL and NDL contributions as well as the collinear endpoint NNDL contribution, figure 2a. The DL-dressing factor $N_b(L) + N_c(L) - N_a(L)$ can be interpreted as if a hard parton of flavour a — which would have been dressed with
a DL multiplicity $N_a(L)$ — is replaced by two hard partons b and c, dressed with DL multiplicities $N_b(L)$ and $N_c(L)$ respectively. Finally, the third term in eq. (3.18) can simply be interpreted as the change in quark and gluon fractions at NLO, including both real and virtual corrections. While the individual contributions from the second and third terms depend on the precise recipe to define the combined flavour of a flavour pair (bc), their sum is well-defined.

Let us briefly return to the question of what happens if in a $2 \to 3$ real event each of the three partons are in separate jets. In this case, the second term in eq. (3.18) disappears (imagine b is in the jet, we have $a = b$ and the square bracket becomes $N_b - N_a = 0$). The last term would then just have one of b or c in the jet which corresponds to the NLO correction to the jet cross-section.

We conclude this discussion by summarising the recipe for obtaining the NNDL hard matrix-element corrections to the Lund multiplicity.

1. Our prescription to define the jet flavour is as follows. If the jet has a single parton, the jet flavour is taken as the flavour of that parton. If instead the jet has 2 partons of flavour indices b, c, the sum $b+c$ is defined as follows: for 2 gluons the sum is a gluon, for a quark and a gluon the sum is a quark, for a same-flavour $q\bar{q}$ pair the sum is a gluon, and for all the other cases the sum is taken as a quark.

2. The Born-level cross-sections $\sigma_0^{(a)}$ for quark and gluon flavours can be obtained from a Monte Carlo simulation including all the required fiducial cuts.

3. With the above flavour prescription, the NLO corrections $\sigma_1^{(a)}$ for quark and gluon jets can be obtained, including fiducial cuts, from any Monte Carlo generator providing information about the flavour of the final-state partons. This defines

$$ \frac{\alpha_s}{2\pi} D_{hme,NLO}^{(a)} = \frac{\sigma_1^{(a)}}{\sigma_0^{(a)}}. \qquad (3.19) $$

4. We compute the second term in (3.18) from a simulation of the real NLO corrections to the cross-section. We select jets with two partons (of flavour b and c) and study the cumulative distribution of their relative $\ln k_t$. For each flavour channel, after subtracting the known DL and NDL contributions, the cumulative distribution tends to a constant at small $\ln k_t$. Further subtracting the NNDL contribution from the collinear endpoint from this constant gives $\sigma_{R,hard}^{(bc)}$. We then define

$$ \frac{\alpha_s}{2\pi} D_{hme,R}^{(bc)} = \frac{\sigma_{R,hard}^{(bc)}}{\sigma_0^{(b+c)}}. \qquad (3.20) $$

A practical example of the extraction of $D_{hme,R}^{(bc)}$ from Monte Carlo simulations is given in appendix B.

The last hard matrix-element correction that has to be taken into account arises by expanding the normalisation in (3.18) as $\frac{1}{\sigma_0+\sigma_1} \sim \frac{1}{\sigma_0}(1-\frac{\sigma_1}{\sigma_0})$. The NNDL hard matrix-
element correction to the Lund multiplicity thus takes the form

$$ \delta N_{\text{hme}}^{(\text{NNDL})} = \sum_{a} \left( \frac{\alpha_s}{2\pi} D_{\text{hme,NLO}}^{(a)} - \frac{\sigma_1}{\sigma_0} \right) f_a N_a^{(\text{DL})}(L) \\ + \sum_{bc} f_{b+c} \frac{\alpha_s}{2\pi} D_{\text{hme,R}}^{(bc)} [N_b^{(\text{DL})}(L) + N_c^{(\text{DL})}(L) - N_a^{(\text{DL})}(L)], \quad (3.21) $$

with $f_a = \frac{\sigma_0^{(a)}}{\sigma_0}$ the Born-level partonic fractions. One can gather the quark and gluon contributions and write

$$ \alpha_s h_{3,\text{hme}}(\xi) = \delta N_{\text{hme}}^{(\text{NNDL})} = \frac{\alpha_s}{2\pi} \left[ D_{\text{hme}}^q N_q^{(\text{DL})}(L) + D_{\text{hme}}^g N_g^{(\text{DL})}(L) \right], \quad (3.22) $$

with

$$ D_{\text{hme}}^{q} = \left( D_{\text{hme,NLO}}^{(q)} + D_{\text{hme,R}}^{(qg)} + D_{\text{hme,R}}^{(qq)} - \frac{2\pi\sigma_1}{\alpha_s\sigma_0} \right) f_q + 2D_{\text{hme,R}}^{(q\bar{q})} f_g, \quad (3.23a) $$

$$ D_{\text{hme}}^{g} = \left( D_{\text{hme,NLO}}^{(g)} + D_{\text{hme,R}}^{(gg)} - D_{\text{hme,R}}^{(q\bar{q})} - \frac{2\pi\sigma_1}{\alpha_s\sigma_0} \right) f_g. \quad (3.23b) $$

### 3.3.2 Interplay between the top and bottom of the Lund plane

Let us now focus on the contribution associated with diagram 4b, corresponding to an endpoint contribution for the radiation of a soft gluon at large angles. Going back to the derivation of the NDL large-angle contribution from section 3.2, we note that, after we changed variable from $k_\perp$ to $k_t = k_\perp \Delta$, the integration over $k_t$ gave both the NDL contribution and an NNDL correction. More precisely, eq. (3.8) can be rewritten as

$$ \int_{k_{t,\text{cut}}}^{p_{\perp}\Delta_k} \frac{dk_t}{k_t} = \int_{k_{t,\text{cut}}}^{p_{\perp}R} \frac{dk_t}{k_t} - \int_{p_{\perp}\Delta_k}^{p_{\perp}R} \frac{dk_t}{k_t} = L + (\text{NNDL}_{\text{top}}), \quad (3.24) $$

where we see that the NNDL correction is associated with the hard $k_t$ scale and can therefore be absorbed into the hard matrix-element correction computed in the previous section. As a corollary, the contribution at the bottom of the Lund plane, diagram 4b, vanishes.

The argument can be extended to the case of diagram 4c, where a first large-angle correction would be followed by a secondary emission close to the $k_{t,\text{cut}}$ threshold. For this, we consider the all-order expression obtained by adding a tower of soft-collinear gluons to a first large-angle emission as given by eq. (3.7) (assuming for simplicity that the (ab) dipole does not involve the jet leading parton):

$$ \delta N_{C,(ab),\mathcal{O}(\alpha_s)} = \frac{\alpha_s}{2\pi} \omega_{ab}^c \int_0^{p_{\perp}} \frac{dk_{\perp}}{k_{\perp}} \int_{\Delta_k < R} dy \frac{d\phi}{2\pi} k_{\perp}^2 (k|ab) \Theta(k_t > k_{t,\text{cut}}) N^{(\text{DL})} \left( L; \ln \frac{p_{\perp}R}{k_{\perp}\Delta_k} \right), \quad (3.25) $$

where the last factor includes soft-collinear gluon radiation from the $k_t$ scale of the first emission, $k_\perp \Delta_k$, down to the cut-off scale $k_{t,\text{cut}}$. Changing variables from $k_\perp$ to $k_t = k_\perp \Delta_k$, one recovers the same picture as above, namely the NDL contribution plus an NNDL correction when the first emission is at the top of the Lund plane (the second term in eq. (3.24)) dressed with DL soft-collinear radiation. In particular, written in terms of $k_t$,
the dressing factor $N^{(\text{DL})}(L; \ln \frac{p_{\perp}R}{k_{\perp}\Delta_k})$ is simply $N^{(\text{DL})}(L; \ln \frac{p_{\perp}R}{k_t})$, i.e. the DL result with no NDL corrections. This means that the contribution from diagram 4c also vanishes.

This vanishing of diagrams 4b and 4c carries some degree of arbitrariness and indicates a freedom in the calculation to reshuffle contributions across different corrections. More concretely, a fraction of the contribution associated with the top of the Lund plane can be redistributed into a contribution at the bottom of the Lund plane (starting at $\mathcal{O}(\alpha_s)$) and a contribution associated with diagram 4c (starting at $\mathcal{O}(\alpha_s^2)$). This reorganisation is done without changing the all-order NNDL final result.<fnref n="7" /> To illustrate this, let us take a different approach where we instead change variable from $k_{\perp}$ to, e.g., $\tilde{k}_{\perp} = k_{\perp}R$. In this case, the integration with a single large-angle gluon, i.e. the analogue of eq. (3.24), would give

$$ \int_0^{p_{\perp}R} \frac{d\tilde{k}_{\perp}}{\tilde{k}_{\perp}} \Theta(k_t > k_{t,\text{cut}}) = \int_{k_{t,\text{cut}}}^{p_{\perp}R} \frac{dk_t}{k_t} = \int_{k_{t,\text{cut}}}^{p_{\perp}R} \frac{dk_t}{k_t} + \int_{k_{t,\text{cut}}}^{k_{t,\text{cut}}} \frac{dk_t}{k_t} = L + (\text{NNDL}_{\tilde{k}_{\perp},\text{bottom}}), \quad (3.26) $$

with

$$ (\text{NNDL}_{\tilde{k}_{\perp},\text{bottom}}) = (\text{NNDL}_{\text{top}}) = \ln \frac{\Delta_k}{R}, \quad (3.27) $$

with $(\text{NNDL}_{\text{top}})$ the contribution from the top of the Lund plane from eq. (3.24). After integrating over $\Delta_k$, this would give a constant. At order $\alpha_s$, this is simply equivalent to reshuffling an NNDL correction from the top to the bottom of the Lund plane. At all orders, the contributions from the top and bottom of the Lund plane are dressed with different DL contributions. In practice, one can show that since the contribution from diagram 4c is also non-zero, the all-order result remains the same independently of the use of $\tilde{k}_{\perp}$ or $k_{\perp}$, as expected.

### 3.3.3 Two large-angle emissions

We now investigate the case in which a Born-level configuration $(p_1p_2 \rightarrow p_jp_r)$ emits a soft-and-large-angle gluon, $k_1$, which then radiates a soft emission $k_2$ at a commensurate angle, $\Delta_{k_1k_2} \sim \Delta_{jk_1}$. Due to the purely soft nature of this correction, we can again split it in a quark- and a gluon-jet contribution

$$ h_{3,\text{la}^2}(\xi) = f_q h_{3,\text{la}^2}^{(q)}(\xi) + f_g h_{3,\text{la}^2}^{(g)}(\xi) \quad \text{with} \quad \alpha_s h_{3,\text{la}^2}^{(i)}(\xi) = \sum_C \frac{\sigma_C^{(i)}}{\sigma^{(i)}} \sum_{(ab)\in\text{event}} \delta N_{C,(ab),\text{la}^2}^{(\text{NNDL})}, \quad (3.28) $$

where we decomposed the hard Born-level process into a sum over partonic channels, $C$, and coloured dipoles $(ab)$, as in section 3.2. While the methods in section 3.2 allow us to account for the primary emission at full colour, we here adopt a simplified approach and describe the second emission in the large-$N_c$ limit. If the gluon $k_1$ is radiated from dipole $(ab)$, we thus only have to consider, for a real $k_1$, the radiation of gluon $k_2$ from either the $(ak_1)$ or the $(k_1b)$ dipole and, for a virtual $k_1$, the radiation of $k_2$ from the $(ab)$ dipole. At

<footnote n="7">⁷A similar effect is present in the collinear limit where there is some freedom in reshuffling contributions between the top of the Lund plane, the collinear endpoint and a hard-collinear emission with a second emission at commensurate angles, see section 4.3.3 of ref. [1].</footnote>
$O(\alpha_s^2)$, the Lund multiplicity is thus given by

$$ \delta N_{C,(ab),\text{la}^2,\mathcal{O}(\alpha_s^2)} = \left(\frac{\alpha_s}{2\pi}\right)^2 \omega_{ab}^C N_c \int_0^{p_{\perp}} \frac{dk_{\perp 1}}{k_{\perp 1}} \int_0^{k_{\perp 1}} \frac{dk_{\perp 2}}{k_{\perp 2}} \int dy_1 dy_2 \frac{d\phi_1}{2\pi} \frac{d\phi_2}{2\pi} \\ \times k_{\perp 1}(k_1|ab) k_{\perp 2}[(k_2|ak_1)+(k_2|bk_1)-(k_2|ab)] \Theta(\Delta_2 < R) \Theta(k_{t2} > k_{t,cut}), \quad (3.29) $$

with $\Delta_i$ the distance between gluon $k_i$ and the jet axis. Here, the eikonal factors were defined in eq. (3.5), $\Theta(\Delta_2 < R)$ ensures that $k_2$ is emitted inside of the jet, and we made explicit the use of the large-$N_c$ approximation for the second emission. Note that the gluon $k_1$ can be emitted either inside or outside of the jet. In writing eq. (3.29), we have implicitly assumed that the jet had been clustered with the anti-$k_t$ algorithm, allowing us to write the condition that the second emission is in the jet as $\Theta(\Delta_2 < R)$. For a different jet algorithm, this would have to be adapted, hence potentially yielding a different result at NNDL accuracy.

In order to isolate the NNDL terms in eq. (3.29), we first subtract the DL and NDL contributions. For dipoles not involving the jet ((12), (1r) and (2r)), $(k_1|ab)$ has no collinear singularity when $\Delta_1 \to 0$. The only contribution up to NDL is thus that of $k_1$ being emitted inside the jet at large angles with $k_2$ being a DL soft-and-collinear dressing. One therefore has the following subtracted integrand:

$$ I_{ab} = k_{\perp 1}^2(k_1|ab) \left\{ k_{\perp 2}^2[(k_2|ak_1) + (k_2|k_1b) - (k_2|ab)]\Theta(\Delta_2 < R) - \frac{4}{\Delta_{12}^2}\Theta(\Delta_1, \Delta_{12} < R) \right\}. \quad (3.30) $$

For dipoles involving the jet ((1j), (2j) and (rj)), there are two contributions to subtract: a DL one where both gluons are soft and collinear with the second radiated as a secondary emission from the first one, and a NDL one where $k_1$ occurs at large angles and $k_2$ is a DL soft-and-collinear dressing. This gives

$$ I_{aj} = k_{\perp 1}^2(k_1|1j)k_{\perp 2}^2[(k_2|ak_1) + (k_2|jk_1) - (k_2|aj)]\Theta(\Delta_2 < R) \\ - \frac{2}{\Delta_1^2} \frac{4}{\Delta_{12}^2} \Theta(\Delta_{12} < \Delta_1 < R) - \left[k_{\perp 1}^2(k_1|aj) - \frac{2}{\Delta_1^2}\right] \frac{4}{\Delta_{12}^2} \Theta(\Delta_1, \Delta_{12} < R). \quad (3.31) $$

In writing the subtraction terms, we have split the radiation of the first soft gluon into a collinear contribution $2/\Delta_1^2$ and a large-angle contribution $k_{\perp 1}^2(k_1|aj) - 2/\Delta_1^2$. After performing these subtractions, the angular integrations are finite and do not show any collinear enhancement. At $O(\alpha_s^2)$ we can thus trivially perform the $k_{\perp}$ integrals by replacing $k_{t2} \to k_{\perp 2}R$ (see eq. (3.8)), which simply gives a factor $L^2/2$. The NNDL correction to the multiplicity is then

$$ \delta N_{C,(ab),\text{la}^2,\mathcal{O}(\alpha_s^2)}^{(\text{NNDL})} = \left(\frac{\alpha_s}{2\pi}\right)^2 \omega_{ab}^C N_c \frac{L^2}{2} D_{ab}^{\text{la}^2}, \quad (3.32) $$

with

$$ D_{ab}^{\text{la}^2} = \int_{-\infty}^{\infty} dy_1 \int_{-\infty}^{\infty} dy_2 \int_{-\pi}^{\pi} \frac{d\phi_1}{2\pi} \int_{-\pi}^{\pi} \frac{d\phi_2}{2\pi} I_{ab}. \quad (3.33) $$
<figure><img src="image_7.png" /><figcaption>**Figure 5.** NNDL coefficients $D_{ab}^{\text{la}^2}$ for each dipole coming from two emissions at large angles. The results are provided for both $R = 0.4$ and $R = 0.8$ jets. The table on the right gives numerical values for $\hat{y} = 0$.</figcaption></figure>

We could not find analytic expressions for the integrals in eq. (3.33) so we evaluated them numerically.<fnref n="8" /> As for the NDL $D_{ab}^{\text{la}}$ coefficients, the $D_{ab}^{\text{la}^2}$ coefficients involving the recoiling jet depend on the rapidity difference between the jet and the recoiling jet, $2\hat{y}$ in the notations of section 3.2. For reproducibility purposes, we show the coefficients $D_{ab}^{\text{la}^2}$ in figure 5 together with the numerical values obtained for $\hat{y} = 0$ for two different jet radii. Furthermore, $D_{2j}^{\text{la}^2} = D_{1j}^{\text{la}^2}$ and $D_{2r}^{\text{la}^2}$ is obtained from $D_{1r}^{\text{la}^2}$ taking $\hat{y} \to -\hat{y}$. As expected, $D_{1r}^{\text{la}^2}$ goes to 0 at large negative $\hat{y}$ and to $D_{12}^{\text{la}^2}$ at large positive $\hat{y}$. Similarly, $D_{jr}^{\text{la}^2}$ tends to $D_{1j}^{\text{la}^2}$ at large $|\hat{y}|$.

At all orders, since this correction affects just the first two emissions, we only need to account for subsidiary DL radiation from the second gluon. That is, the NNDL correction when $k_1$ is emitted from the (ab) dipole is

$$ \delta N_{C,(ab),\text{la}^2}^{(\text{NNDL})} = \left(\frac{\alpha_s}{2\pi}\right)^2 \omega_{ab}^C N_c D_{ab}^{\text{la}^2} \int_0^L d\ell_1 \int_{\ell_1}^L d\ell_2 N_g^{(\text{DL})}(L-\ell_2) = \frac{\alpha_s}{2\pi} \omega_{ab}^C D_{ab}^{\text{la}^2} \frac{1}{4}(\cosh \nu - 1). \quad (3.34) $$

The total result is obtained after summing over all possible coloured dipoles and partonic channels as written in eq. (3.28). We remind the reader that the kinematic weights $\omega_{ab}^C$ for inclusive jet (Z+jet) events are given in appendix A (table 2).

### 3.3.4 A large-angle emission and a hard-collinear emission

This corresponds to the case of a soft-and-large-angle gluon emission together with a subsidiary hard-collinear emission (either a $g \to gg$ or a $g \to q\bar{q}$ splitting) as displayed in

Values for $\hat{y} = 0$

<table><thead><tr><th>R</th><th>0.4</th><th>0.8</th></tr></thead><tbody><tr><td>$D_{12}^{\text{la}^2}$</td><td>0.5864261(12)</td><td>0.5712475(12)</td></tr><tr><td>$D_{1j}^{\text{la}^2} = D_{2j}^{\text{la}^2}$</td><td>-0.24098(15)</td><td>-0.97564(17)</td></tr><tr><td>$D_{1r}^{\text{la}^2} = D_{2r}^{\text{la}^2}$</td><td>0.259460(11)</td><td>0.582097(13)</td></tr><tr><td>$D_{jr}^{\text{la}^2}$</td><td>-0.0017776(17)</td><td>-0.0285608(17)</td></tr></tbody></table>

<footnote n="8">⁸We could tentatively do an expansion in the jet radius $R$ as in section 3.2 but the coefficients of this expansion would also have to be evaluated numerically so we preferred the fully numerical solution. We note however that, at small $R$, the $D_{ab}^{\text{la}^2}$ coefficients behave like $R^2 \ln R^2$. This behaviour comes from configurations where the first emission is outside the jet but very close to the jet boundary with the second emission inside the jet and also close to the jet boundary so that the $\ln R^2$ enhancement is of collinear origin. See e.g. [23, 55] for additional discussions.</footnote>
figure 4e. As in the above section we split the contribution with respect to the jet flavour

$$h_{3,\text{la}\times\text{hc}}(\xi) = f_q h_{3,\text{la}\times\text{hc}}^{(q)}(\xi) + f_g h_{3,\text{la}\times\text{hc}}^{(g)}(\xi) \quad \text{with} \quad (3.35)$$
$$\alpha_s h_{3,\text{la}\times\text{hc}}^{(i)}(\xi) = \sum_C \frac{\sigma_C^{(i)}}{\sigma^{(i)}} \sum_{(ab)\in\text{event}} \delta N_{C,(ab),\text{la}\times\text{hc}}^{(\text{NNDL})}.$$

The starting point for computing $\delta N_{C,(ab),\text{la}\times\text{hc}}^{(\text{NNDL})}$ is the result of section 3.2, where we calculated the correction to the Lund multiplicity due to a primary, large-angle emission. Instead of dressing this soft, large-angle, gluon with a tower of DL emissions, we now have to account for a single hard-collinear correction anywhere in this chain of gluons. For this, we can recycle the NDL hard-collinear correction to the multiplicity, $\delta N_{g,\text{hc}}^{(\text{NNDL})}(L; \ell)$, computed between the relative transverse momentum scale $e^{-\ell}p_{\perp}R$ at which the first large-angle gluon is emitted and the cut-off scale $k_{t,\text{cut}} = e^{-L}p_{\perp}R$ (see section 3.2.2 of ref. [1] for explicit expressions). We therefore have

$$\begin{aligned}
\delta N_{C,(ab),\text{la}\times\text{hc}}^{(\text{NNDL})} &= \frac{\alpha_s}{\pi} \omega_{ab}^{C} D_{ab}^{\text{la}} \int_0^L d\ell \, \delta N_{g,\text{hc}}^{(\text{NNDL})}(L; \ell) && (3.36) \\
&= \frac{\alpha_s}{2\pi} \omega_{ab}^{C} D_{ab}^{\text{la}} [(B_{gg} + c_{\delta} B_{gq})\nu \sinh \nu + 2B_{gq}(1 - c_{\delta})(\cosh \nu - 1)],
\end{aligned}$$

where $c_{\delta} = (2C_F - C_A)/C_A$. The total result is obtained after summing over all possible coloured dipoles and partonic channels as written in eq. (3.35). Again, the kinematic weights $\omega_{ab}^C$ can be found in table 2 for Z+jet events and in appendix A for an inclusive jet sample.

### 3.3.5 A large-angle emission and a running-coupling correction

The structure of this correction is similar to the one in the previous section with the main difference being that the running-coupling correction can affect either the soft large-angle emission itself, or any subsidiary soft-and-collinear emission. We therefore write

$$\begin{aligned}
\delta N_{C,(ab),\text{la}\times\text{rc}}^{(\text{NNDL})} &= \frac{\alpha_s}{\pi} \omega_{ab}^{C} D_{ab}^{\text{la}} \int_0^L d\ell \left[ 2\alpha_s \beta_0 \ell N_g^{(\text{DL})}(L-\ell) + \delta N_{g,\text{rc}}^{(\text{NNDL})}(L;\ell) \right] && (3.37) \\
&= \frac{\alpha_s}{2C_A} \omega_{ab}^{C} D_{ab}^{\text{la}} \frac{\beta_0}{2} \nu [\nu \cosh \nu + \sinh \nu],
\end{aligned}$$

where $\delta N_{g,\text{rc}}^{(\text{NNDL})}(L; \ell)$ is the NDL resummed running-coupling correction of a gluon radiating from scale $k_t = p_{\perp}Re^{-\ell}$ to $k_{t,\text{cut}}$ (see eq. (4.60) of ref. [1] for an explicit expression). The first (second) term in the integrand corresponds to dressing the large-angle emission (any subsidiary emission) with a running-coupling correction, respectively. As before, the full NNDL correction $h_{3,\text{la}\times\text{rc}}(\xi)$ is obtained by summing the contributions in eq. (3.37) for each dipole and partonic channel, as in eq. (3.35).

### 3.4 Final result

For completeness, we summarise here the necessary expressions needed to compute the average Lund multiplicity within a high-energy anti-$k_t$ jet up to NNDL accuracy, i.e. up to $h_3$ in eq. (1.1). We define $\alpha_s \equiv \alpha_s(p_{\perp}R)$, $L \equiv \ln(p_{\perp}R/k_{t,\text{cut}})$, $\xi \equiv \alpha_s L^2$, $\nu \equiv \sqrt{2C_A\xi/\pi}$, $c_{\delta} \equiv$
$(2C_F - C_A)/C_A$ and $f_{q,g} \equiv \sigma_{q,g}/(\sigma_q + \sigma_g)$ the Born-level fractions of quark and gluon jets. The DL function is

$$h_1(\xi) = f_q h_{1,\text{coll}}^{(q)}(\xi) + f_g h_{1,\text{coll}}^{(g)}(\xi) \quad \text{with} \quad h_{1,\text{coll}}^{(i)}(\xi) = 1 + \frac{C_i}{C_A} (\cosh \nu - 1). \quad (3.38)$$

The NDL function is

$$h_2(\xi) = f_q h_{2,\text{coll}}^{(q)}(\xi) + f_g h_{2,\text{coll}}^{(g)}(\xi) + \frac{1}{\sqrt{2\pi C_A}} \sum_{i \in (q,g)} \sum_{C \in \text{flavour channels}} \frac{\sigma_C^{(i)}}{\sigma^{(i)}} \sum_{(ab) \in \text{event}} \omega_{ab}^C D_{ab}^{\text{la}} \sinh \nu, \quad (3.39)$$

and the NNDL function is

$$2\pi h_3(\xi) = 2\pi [f_q h_{3,\text{coll}}^{(q)}(\xi) + f_g h_{3,\text{coll}}^{(g)}(\xi)] + D_{\text{hme}}^q \left[ \frac{C_F}{C_A} (\cosh \nu - 1) + 1 \right] + D_{\text{hme}}^g \cosh \nu \\ + \sum_{i \in (q,g)} \sum_{C \in \text{flavour channels}} \frac{\sigma_C^{(i)}}{\sigma^{(i)}} \sum_{(ab) \in \text{event}} \omega_{ab}^C \left\{ D_{ab}^{\text{la}2} \frac{1}{4} (\cosh \nu - 1) + D_{ab}^{\text{la}} \frac{\pi \beta_0}{2C_A} \nu [\nu \cosh \nu + \sinh \nu] \right. \\ \left. + D_{ab}^{\text{la}} [(B_{gg} + c_\delta B_{gq})\nu \sinh \nu + 2B_{qg}(1 - c_\delta)(\cosh \nu - 1)] \right\}. \quad (3.40)$$

In the above expressions, we have recycled the collinear pieces of the $e^+e^-$ results, i.e. $h_{2,\text{coll}}^{(i)}(\xi)$ is defined eq. (2.4), and $h_{3,\text{coll}}^{(q)}(\xi)$ and $h_{3,\text{coll}}^{(g)}(\xi)$ are defined in eq. (2.7) and (2.9), respectively.

## 3.5 From fixed kinematics to fiducial cuts

So far, we have considered contributions to the average Lund multiplicity stemming from jets with a fixed transverse momentum and rapidity. In order to provide realistic predictions for a collider environment, we have to relax that assumption and average our analytic result with the cross-section for producing jets within a given acceptance, i.e.

$$\langle \overline{N^{(\text{Lund})}} \rangle (\alpha_s; \ln(p_{\perp,\text{min}} R/k_{t,\text{cut}})) = \frac{1}{\sigma} \int_{p_{\perp,\text{min}}} \frac{d^2\sigma}{dp_{\perp} dy} \langle N^{(\text{Lund})} (\alpha_s; L_{\text{jet}}) \rangle dp_{\perp} dy, \quad (3.41)$$

where $p_{\perp,\text{min}}$ is the minimum jet transverse momentum selected, $d^2\sigma/dp_{\perp} dy$ is the double-differential distribution in both the transverse momentum $p_{\perp}$ and rapidity $y$ of the jet. Finally, $\sigma$ is the total jet cross-section, meaning that we are effectively computing the average multiplicity per jet passing the cuts. At NNDL accuracy, our expressions for the multiplicity already include, through the contribution from the top of the Lund plane derived in section 3.3.1, the non-trivial NLO corrections from the hard matrix element. The convolution in eq. (3.41) should therefore be performed using the Born-level cross-section.

In practice, one could compute the average multiplicity in jets as a function of $p_{\perp}$ and $y$, and perform the integration in eq. (3.41) explicitly. However, this integration can be simplified by realising that, due to the steeply falling nature of the jet spectrum, $p_{\perp}$ and $p_{\perp,\text{min}}$ are commensurate. One can thus replace $L_{\text{jet}}$ by $L \equiv L_{\text{min}} = \ln(p_{\perp,\text{min}} R/k_{t,\text{cut}})$ as
the resummation scale in eq. (1.1). Since $L_{jet} - L \equiv \delta_L = \ln(p_{\perp}/p_{\perp,\min})$ does not bring new large logarithms, one can expand the right hand side of eq. (3.41), up to NNDL accuracy, as

$$ \langle N^{\text{(Lund)}}(\alpha_s; L_{\text{jet}}) \rangle = h_1(\xi) + \sqrt{\alpha_s} [h_2(\xi) + 2\delta_L \sqrt{\xi} h'_1(\xi)] \quad (3.42) \\ + \alpha_s [h_3(\xi) + (\delta_L^2 - 2\beta_0 \delta_L \xi) h'_1(\xi) + 2\delta_L^2 \xi h''_1(\xi) + 2\delta_L \sqrt{\xi} h'_2(\xi)], $$

with $\xi = \alpha_s L^2$, $\alpha_s = \alpha_s(p_{\perp,\min}R)$ and $h'_i(\xi)$ ($h''_i(\xi)$) the first (second) derivative of $h_i(\xi)$ with respect to $\xi$. The main advantage of this approach is that the resummation functions $h_i(\xi)$, given by eqs. (3.38)–(3.40) are independent of the jet $p_{\perp}$ and $y$. The integration over the jet spectrum in eq. (3.41) thus reduces to computing a few specific quantities such as

$$ \bar{\delta}_L = \frac{1}{\sigma} \int \frac{d^2\sigma}{dp_{\perp}dy} \ln \frac{p_{\perp}}{p_{\perp,\min}} dp_{\perp}dy, \quad (3.43a) $$

$$ \bar{\delta}_L^2 = \frac{1}{\sigma} \int \frac{d^2\sigma}{dp_{\perp}dy} \ln^2 \frac{p_{\perp}}{p_{\perp,\min}} dp_{\perp}dy, \quad (3.43b) $$

as well as quantities related to the radiation of soft gluons at large angles entering at NDL and NNDL, i.e.<fnref n="9" />

$$ \overline{\omega D_i^{\text{fl}}} = \sum_{C \in \text{flavour channels}} \frac{1}{\sigma^{(i)}} \int \frac{d^2 \sigma_C^{(i)}}{dp_{\perp} dy} \sum_{(ab) \in \text{event}} \omega_{ab}^C D_{ab}^{\text{fla}} dp_{\perp} dy, \quad (3.44a) $$

$$ \overline{\omega D_i^{\text{fla}^2}} = \sum_{C \in \text{flavour channels}} \frac{1}{\sigma^{(i)}} \int \frac{d^2 \sigma_C^{(i)}}{dp_{\perp} dy} \sum_{(ab) \in \text{event}} \omega_{ab}^C D_{ab}^{\text{fla}^2} dp_{\perp} dy. \quad (3.44b) $$

In eqs. (3.43) and (3.44), we can again perform the integration using the Born-level matrix element. All the coefficients entering eq. (3.42) depend on the underlying hard process as well as the imposed fiducial cuts.

## 4 Numerical results

In this section we discuss some phenomenological considerations for the average Lund jet multiplicity, including the matching of our resummed predictions to exact fixed-order results, as well as non-perturbative corrections. For definiteness, we specialise our discussion to the case of $Z+jet$ events and an inclusive jet sample in proton-proton collisions at a centre-of-mass energy $\sqrt{s} = 13.6$ TeV. In the $Z+jet$ case we impose that $p_{\perp Z} > p_{\perp,\min} = 500$ GeV and $|y_Z| < 2.5$ (for simplicity, we have disregarded the $Z$ decay in the fixed-order simulations). If the event passes these cuts, we then calculate the Lund multiplicity within the leading $R = 0.8$ anti-$k_t$ jet, as obtained from FastJet [56], with $p_{\perp,\text{jet}} > p_{\perp Z}/2$ and $|y_{\text{jet}}| < 1.7$. In the inclusive jet sample, we select all $R = 0.8$ anti-$k_t$ jets with $p_{\perp,\text{jet}} > p_{\perp,\min} = 750$ GeV and $|y_{\text{jet}}| < 1.7$.

As discussed in section 3.5, the NNDL-accurate resummed predictions involve a series of coefficients which depend on the underlying Born-level process. We obtain the $Z+jet$ ones

<footnote n="9">⁹In practice, one needs to integrate over the full kinematics of the $2 \to 2$ Born event since the $\omega_{ab}^C D_{ab}^{\text{fla}}$ and $\omega_{ab}^C D_{ab}^{\text{fla}^2}$ coefficients depend on the relative $\hat{y}$ between the measured jet and the recoiling Born-level parton.</footnote>
<table><caption>**Table 3.** Coefficients entering the NDL $h_2$ and NNDL $h_3$ functions for our two jet samples, inclusive jets and Z+jet events with $R = 0.8$ in both cases.</caption><thead><tr><th></th><th colspan="2">Coefficient</th></tr><tr><th></th><th>Incl. jets</th><th>Z+jet</th></tr></thead><tbody><tr><td>kinematic cuts</td><td>p<sub>&perp;,jet</sub> &gt; 750 GeV<br>(p<sub>&perp;,min</sub> = 750 GeV)<br>|y<sub>jet</sub>| &lt; 1.7</td><td>p<sub>&perp;,Z</sub> &gt; 500 GeV, p<sub>&perp;,jet</sub> &gt; 1/2 p<sub>&perp;,Z</sub><br>(p<sub>&perp;,min</sub> = 500 GeV)<br>|y<sub>Z</sub>| &lt; 2.5, |y<sub>jet</sub>| &lt; 1.7</td></tr><tr><td>strong coupling</td><td>&alpha;<sub>s</sub>(Rp<sub>&perp;,min</sub>) = 0.091</td><td>&alpha;<sub>s</sub>(Rp<sub>&perp;,min</sub>) = 0.097</td></tr><tr><td>quark, gluon jet fractions</td><td>f<sub>q</sub> = 0.553<br>f<sub>g</sub> = 0.447</td><td>f<sub>q</sub> = 0.818<br>f<sub>g</sub> = 0.182</td></tr><tr><td>spectrum convolutions</td><td>&delta;&#773;<sub>L</sub> = -0.166<br>&delta;&#773;<sup>2</sup><sub>L</sub> = 0.052</td><td>&delta;&#773;<sub>L</sub> = -0.195<br>&delta;&#773;<sup>2</sup><sub>L</sub> = 0.071</td></tr><tr><td>hard matrix-element</td><td>D<sup>q</sup><sub>hme</sub> = -9.331<br>D<sup>g</sup><sub>hme</sub> = 5.508</td><td>D<sup>q</sup><sub>hme</sub> = -11.545<br>D<sup>g</sup><sub>hme</sub> = 8.285</td></tr><tr><td>two large-angle emissions</td><td>&omega;D<sup>la<sup>2</sup></sup><sub>q</sub> = 0.030<br>&omega;D<sup>la<sup>2</sup></sup><sub>g</sub> = -0.246</td><td>&omega;D<sup>la<sup>2</sup></sup><sub>q</sub> = -0.074<br>&omega;D<sup>la<sup>2</sup></sup><sub>g</sub> = -0.504</td></tr></tbody></table>

from MadGraph [57] and the inclusive jet ones from a version of NLOJet++ [58] augmented so as to provide parton flavour information [59]. All the results reported here use the central PDF of the PDF4LHC21 set [60]. We set the central renormalisation scale to $p_{\perp,\min}R$ and the central factorisation scale to $\hat{H}_T$, the scalar sum of the transverse momentum of the final-state partons. The values for the resummation coefficients with this set of choices are listed in table 3.

## 4.1 Resummation

Throughout this paper, we have organised the calculation of the average Lund multiplicity as the sum of a collinear piece, identical to the $e^+e^-$ result, and a large-angle piece. The collinear contribution is universal in the sense that it is process-independent and only depends on the Born-level jet flavour. Conversely, the large-angle contribution depends on the full structure of the hard process. Up to NNDL accuracy, it can be further split into a contribution associated with soft emissions (close to the edge of the jet), and genuine fixed-order corrections associated with an emission at the top of the Lund plane. The former is still associated with the Born-level jet flavour, but for the latter the breakdown into a quark and a gluon contribution is ambiguous.

To gain insight into the relative sizes of each of these contributions we plot in figure 6 the relative weights of the different contributions for both the inclusive jet sample and Z+jet
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Relative weight of each term to the total Lund multiplicity at NNDL accuracy for an inclusive sample of jets (left) and Z+jet events (right). We distinguish universal, collinear terms (blue) from those contributions which depend on the event structure, i.e. wide-angle radiation (red) and hard matrix-element corrections (grey).</figcaption></figure>

events. Both the soft and the collinear contributions have explicitly been split according to flavour, and the contribution labelled as “non-collinear” includes both the contributions from soft-large-angle radiation as well as the contributions from the top of the Lund plane. For both hard processes, almost 90% of the Lund multiplicity is dominated by contributions of collinear origin. Focusing on the large-angle corrections, one sees that the contribution from hard emissions, i.e. at the top of the Lund plane, is further suppressed (roughly by a factor 2) compared to the soft-and-wide-angle contributions. The hard contribution is however larger for the Z+jet case than for inclusive jets, and, in both cases, changes sign as a function of $k_{t,cut}$. Globally speaking, this picture is in agreement with the fact that collinear effect start at DL accuracy, soft wide-angle effects start at NDL and top-of-the-Lund-plane corrections only contribute at NNDL. The fact that the number of emissions inside a jet is dominated by universal collinear physics depending only on the Born-level flavour of the jet is in quantitative agreement with the study of the dependence of quark and gluon jets on the topology of the hard process from ref. [61].

Regarding the flavour decomposition of the result, the relative sizes of quark and gluon contributions differ drastically between the inclusive and the Z+jet sample. As expected, we find a larger relative contribution of quark-initiated jets in Z+jet events than in the inclusive sample. In the latter, the slightly higher initial fraction of quark-initiated jets (see table 3) is compensated by Casimir scaling for small values of $k_{t,cut}$, i.e. gluons give a larger multiplicity than quarks at small $k_{t,cut}$. However, quark jets dominate down to $k_{t,cut} \sim 30$ GeV.

## 4.2 Phenomenological predictions

To produce realistic predictions for phenomenological applications at the LHC we match our resummed calculation, computed in section 3, to the exact fixed-order result at NLO,
i.e. $O(\alpha_s^2)$ compared to the underlying $2 \to 2$ hard process.<fnref n="10" /> We further account for theoretical uncertainties and include non-perturbative corrections. We perform these three steps in this section, focusing only on our inclusive jet sample for simplicity.

**Theoretical uncertainties.** We evaluate the uncertainties associated with missing higher orders by varying the renormalisation ($p_{\perp,\min}R$), factorisation ($\hat{H}_T$) and resummation scales ($p_{\perp,\min}R/k_{t,\text{cut}}$). To this end we multiply each of the scales by a dimensionless parameter, $x_R$, $x_F$ and $x_L$ respectively. Variations of $x_F$ only affect quantities that depend on the averaging over the Born-level jet spectrum, i.e. quantities listed in table 3. Conversely, when varying $x_R$ and $x_L$ one introduces spurious terms that spoil the perturbative accuracy of the resummation. To maintain NNDL accuracy, we introduce counter-terms so that our resummed expression becomes

$$ \langle \overline{N^{\text{(Lund)}}} \rangle_{\text{resum}} = \bar{h}_1(\xi) + \sqrt{\alpha_s} [h_2(\xi) + 2(\bar{\delta}_L - \ln x_L)\sqrt{\xi} \bar{h}'_1(\xi)] \quad (4.1) \\ + \alpha_s \{ \bar{h}_3(\xi) + [2\beta_0(\ln x_R - \bar{\delta}_L)\xi + \ln^2 x_L - 2\ln x_L \bar{\delta}_L + \bar{\delta}_L^2] \bar{h}'_1(\xi) \\ + 2(\ln^2 x_L - 2\ln x_L \bar{\delta}_L + \bar{\delta}_L^2)\xi \bar{h}''_1(\xi) + 2\sqrt{\xi} [\bar{\delta}_L \bar{h}'_2(\xi) - \ln x_L \bar{h}'_2(\xi)] \}, $$

where the bar indicates an average over the jet spectrum as explained in section 3.5. Note that, in the last term, we have used the short-hand notation $\bar{\delta}_L h'_2(\xi)$ to emphasise that, in the large-angle NDL contribution, one has to compute explicitly the overall $\bar{\delta}_L \omega_C D_i^{\text{la}}$ coefficient averaged over Born-level events. All the other, collinear, NDL terms simply involve $\bar{\delta}_L$. In practice, we vary separately $x_R$, $x_F$ and $x_L$ between 0.5 and 2. Our total perturbative uncertainty is obtained by taking the envelope of the 7-point scale variations [62] in $x_{R,F}$, and adding this uncertainty in quadrature with the resummation uncertainty obtained by varying $x_L$.

The numerical evaluation of eq. (4.1) is displayed in figure 7a at DL, NDL and NNDL accuracy. Let us first focus on the central values. We observe that the NDL correction enhances the Lund multiplicity with respect to the DL result when $k_{t,\text{cut}} \lesssim 6$ GeV. We note that this was not the case in the $e^+e^-$ result, where the NDL correction reduced the multiplicity for all $k_{t,\text{cut}}$ values. This increase in the multiplicity of high-energy jets is solely due to the large-angle emission correction calculated in section 3.2. The NNDL corrections further increase the Lund multiplicity for $k_{t,\text{cut}} \lesssim 15$ GeV, and have only moderate impact for all larger values of $k_{t,\text{cut}}$. Turning to the uncertainty band, we find that it is remarkably reduced beyond DL accuracy. For example, at $k_{t,\text{cut}} = 5$ GeV (10 GeV), the error band shrinks from a ~28% (~24%) at DL to ~10% (~7%) at NDL accuracy and further to ~5% (~3.4%) at NNDL accuracy, i.e. an additional reduction around 50%. We note that these small relative uncertainties at NNDL accuracy, even for very small $k_{t,\text{cut}}$, can be traced back to the fact that we are computing the cumulative distribution of the average Lund multiplicity. For the differential distribution, we find uncertainties larger by about a factor of 2 at small $k_{t,\text{cut}}$, e.g. 10% (7.5%) at $k_{t,\text{cut}} = 5$ GeV (10 GeV). For

<footnote n="10">¹⁰The first non-trivial order for the Lund multiplicity requires at least one emission on top of the underlying $2 \to 2$ process. This means that the LO and NLO fixed-order multiplicities we refer to in this section correspond respectively to a relative $O(\alpha_s)$ and $O(\alpha_s^2)$ compared to the hard $2 \to 2$ process.</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 7. Left: resummed Lund average multiplicity at DL (dashed, green), NDL (dash-dotted, red) and NNDL (solid, blue) as a function of $k_{t,cut}$ for LHC kinematics. Right: average Lund multiplicity after matching with NLO. The pure NLO result is shown in dotted (grey) for reference, together with matched NLO+NDL (dash-dotted, red) and NLO+NNDL (solid, blue). The bottom panels display the ratio to the most accurate result (NNDL on the left and NLO+NNDL on the right).</figcaption></figure>

$k_{t,cut} \gtrsim p_{\perp,\min}R/2 = 300$ GeV, well beyond the regime of validity of the resummed result, we note the irregular behaviour of the uncertainties associated with the NDL and NNDL accurate results. This unphysical behaviour disappears after matching as we discuss next.

**Matching.** We use an additive scheme for the Lund multiplicity:

$$ \langle \overline{N}^{(\text{Lund})} \rangle_{\text{match}} = \langle \overline{N}^{(\text{Lund})} \rangle_{\text{fo}} \qquad (4.2) \\ + \left( \langle \overline{N}^{(\text{Lund})} \rangle_{\text{resum}} - \langle \overline{N}^{(\text{Lund})} \rangle_{\text{resum,fo}} \right) \left( 1 - \frac{2k_{t,cut}}{p_{\perp,\min}R} \right) \Theta \left( k_{t,cut} < \frac{p_{\perp,\min}R}{2} \right), $$

where $\langle \overline{N}^{(\text{Lund})} \rangle_{\text{fo}}$ is the exact fixed-order result obtained from NLOJet++ taken here at NLO, i.e. at $O(\alpha_s^2)$ relative to the Born-level event, $\langle \overline{N}^{(\text{Lund})} \rangle_{\text{resum}}$ is the NNDL-accurate resummed result, eq. (3.41), and we subtract $\langle \overline{N}^{(\text{Lund})} \rangle_{\text{resum,fo}}$, the $O(\alpha_s^2)$ expansion of the resummed result, to avoid double counting. We also introduce a damping factor to smooth the transition between the resummed and fixed-order results that switches off the resummation at $k_{t,cut}^{\text{max}} = p_{\perp,\min}R/2$ where it should anyway not be trusted. This ensures that the matched results reduce to the fixed-order at large values of $k_{t,cut}$, while not spoiling the NNDL accuracy at small $k_{t,cut}$ values.<fnref n="11" /> We note that scale uncertainties are obtained

<footnote n="11">¹¹One could alternatively have redefined $L$ so that $L \to \ln(p_{\perp,\min}R/k_t)$ at small $k_t$ and $L \to 0$ when $k_t$ is large. Note also that, compared to our earlier $e^+e^-$ study, the $k_t$ spectrum no longer has a well-defined fixed-order endpoint but has instead a tail extending far beyond $p_{\perp,\min}R$.</footnote>
using eq. (4.1) for the resummed part. A similar expression is used for the fixed-order expansion $\langle N^{\text{(Lund)}} \rangle_{\text{resum,fo}}$ where the $h_i(\xi)$ functions and their derivatives are expanded in $\alpha_s$. In the exact fixed-order results we use the $x_R, x_F$ scale variations as provided by NLOJet++.

We show in figure 7b the effect of matching our NDL and NNDL accurate resummed result with the exact NLO result, as described in eq. (4.2). We observe that the resummation slightly increases the value of the Lund multiplicity with an effect more pronounced as $k_{t,cut}$ decreases. As an example, taking the central value of Lund multiplicity curves at $k_{t,cut} = 5$ GeV we find that $\langle N^{\text{(Lund)}} \rangle$ is $\approx 2.6$ at NLO and 2.7 at NLO+NNDL. The uncertainty band after matching increases due to the impact of $x_R$ and $x_F$ variations on the fixed-order result. Nevertheless, the error band on the NLO+NNDL result remains small, with a value of $\sim 6\%$ at 5 GeV and below 15% down to 1 GeV, as observed in the bottom panel of figure 7b. We also note that the matching procedure indeed washes out the unphysical behaviour observed at large $k_{t,cut}$ values for the resummed predictions in figure 7a.

**Non-perturbative corrections.** The last ingredient that we incorporate into our theoretical prediction alongside matching and scale variations is non-perturbative corrections due to hadronisation and multi-parton interactions (MPI). To estimate these we run the Pythia8.3 [63] (with Monash13 [64], 4C [65], and Atlas14 [66] tunes), Herwig7.20 [67] and Sherpa2.2.11 [68] event generators at parton level and at hadron level with MPI turned on.

In figure 8a we show the average Lund multiplicity distributions obtained from our Monte Carlo simulations as well as the ratio between the hadron-level (including MPI) and parton-level results. For readability purposes, we only show the Pythia8 results for the 4C tune, which typically drives the upper edge of the non-perturbative corrections (with the other two Pythia8 tunes showing intermediate non-perturbative corrections). Overall the non-perturbative corrections are small, staying within $\sim 5\%$ down to $k_{t,cut} \sim 3$ GeV.

Below $k_{t,cut} \sim 2$ GeV the non-perturbative corrections start to become significant. However, as noted in ref. [1], our resummation becomes questionable for values of $L$ of order $1/(\beta_0\alpha_s)$. This can be traced to the fact that running-coupling corrections dominate the resummation functions $h_k(\xi)$ at large $\xi$ with contributions proportional to $\xi^{k-1}h_1(\xi)$. In this regime, which corresponds to $k_{t,cut} \lesssim 4$ GeV with our kinematic cuts, one should also consider resumming these running-coupling corrections beyond their expansion to NNDL.

In practice we implement the average of the hadron/parton ratios from the five Monte Carlo setups as a multiplicative factor into the perturbative results. The envelope of these five hadron/parton ratios (the shaded region in the bottom panel of figure 8a) is used as an estimate of the uncertainties associated to non-perturbative effects. When applying non-perturbative corrections to our resummed (or NLO-matched) results, the non-perturbative uncertainties are added in quadrature with the perturbative uncertainties.

The main phenomenological result of this paper is displayed in figure 8b where we show the average Lund multiplicity for our inclusive jet sample. We present results at NNDL accuracy matched to NLO with (black) and without (blue) non-perturbative corrections.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Left: Lund multiplicity distributions obtained from Monte Carlo simulations (upper panel) together with the hadron+MPI/parton ratio (lower panel) used as our estimate of the non-perturbative corrections. Results are shown for Pythia8 (with the 4C tune), Sherpa2 and Herwig7 (with their default tune). The (green) shaded area in the bottom panel corresponds to the envelope of the individual non-perturbative corrections. Right: prediction for the average Lund multiplicity in LHC high-energy jets. The matched NNDL+NLO perturbative results are showed both with (black) and without (blue) non-perturbative corrections. The bottom panel displays the breakdown of uncertainties in their perturbative and non-perturbative contributions.</figcaption></figure>

As expected from figure 8a, the effect of the non-perturbative corrections starts being clearly visible for $k_{t,cut}$ values of a few GeV and below. The theoretical uncertainties are displayed in the lower panel of figure 8b. They are dominated by the perturbative uncertainties across the full $k_{t,cut}$ range, with non-perturbative corrections being negligible for $k_{t,cut} \gtrsim 4$ GeV. This suggests that there is room for further theoretical developments on the perturbative part of the calculation. This includes subleading all-order contributions ($N^{k \ge 3} \text{DL}$), subleading fixed-order results ($N^{k \ge 2} \text{LO}$) [69, 70] or multiple running-coupling correction effects discussed above and relevant for $k_{t,cut} \lesssim 4$ GeV.

## 5 Conclusions

A new jet observable, dubbed Lund multiplicity, has recently been introduced in ref. [1]. It provides a procedure for evaluating subjet multiplicity with the use of an IRC-safe cutoff $k_{t,cut}$ in the spirit of modern Lund-plane declustering techniques [47]. In ref. [1], the event-wide average Lund multiplicity was resummed up to NNDL accuracy in $e^+e^-$ collisions and up to NDL accuracy in colour singlet production in pp collisions. The present
manuscript extends the definition and calculation of the average Lund multiplicity to the case of high-energy jets at hadron colliders by counting the mean number of subjets per anti-$k_t$ jet with relative $k_t$ above $k_{t,cut}$. The main motivation behind this work is to provide a multiplicity-based observable which is perturbatively well-defined in QCD, and can benefit from the large energy range accessible at the LHC. Indeed, for a jet of transverse momentum $p_\perp$ and radius $R$, the Lund multiplicity probes energy scales from $k_{t,cut}$ values around 1 GeV up to $p_\perp R$ values of several hundreds of GeV.

From a theoretical standpoint, the main achievement of this paper is the all-order computation of the average Lund multiplicity up to NNDL accuracy, which is particularly relevant given the large range of transverse momenta probed by this observable. There are two main differences with respect to the $e^+e^-$ event-wide result: the notion of a jet radius, which impacts the large-angle components of the resummation starting at NDL accuracy (see sections 3.2 and 3.3), and the presence of experimental fiducial cuts used for the jet analysis in a collider environment (see section 3.5). For a jet of fixed $p_\perp$ and rapidity, NNDL accuracy is reached by recycling universal contributions sensitive to collinear emissions and running coupling corrections (see section 2) from the results in $e^+e^-$ collisions, and supplementing these with jet-radius effects specific to the case of high-energy jets in pp collisions. The resulting NNDL distribution for a generic $2 \to 2$ hard process, eqs. (3.38)–(3.40), is the main theoretical result of this paper. Reaching NNDL accuracy for a realistic situation with fiducial cuts is achieved via a simple procedure (detailed in section 3.5).

We have performed phenomenological studies specialising our discussion to $Z$+jet and inclusive jet samples at the LHC. Our study has shown (figure 6) that universal, collinear, terms contribute towards $\sim 90\%$ of the average Lund multiplicity in both hard processes. Decomposing contributions in terms of jet flavour, we have found that the $Z$+jet result is overall governed by quark-initiated jets, while for inclusive jets both quark- and gluon-initiated jets contribute commensurately for all values of the transverse momentum cut $k_{t,cut}$.

Furthermore, we have presented (section 4.2) predictions for the average Lund multiplicity of high-energy inclusive jets at the LHC. In order to cover the full range of $k_{t,cut}$ values, we have matched our NNDL calculation to exact NLO fixed-order results. Non-perturbative corrections have been included via a multiplicative factor extracted from general-purpose Monte Carlo event generators. The result of this procedure is shown in figure 8. Theoretical uncertainties remain below 15% down to $k_{t,cut}$ values of 1 GeV and are dominated by their perturbative component. In particular, non-perturbative uncertainties are remarkably small ($< 3\%$) for $k_{t,cut} \gtrsim 4$ GeV, with the overall uncertainty remaining under $\sim 7\%$. This leaves some margin for further theoretical improvements like the inclusion of higher-order effects such as $O(\alpha_s^3)$ matching, N$^3$DL resummation, or joint resummations of running-coupling effects at small $k_{t,cut}$.

In a more generic context, experimental measurements of the number of charged particles within a jet [71, 72] and of the number of subjets inside groomed jets [17] have revealed sizeable discrepancies between Monte Carlo event generators. A high-precision calculation of a multiplicity-based observable can therefore help clarifying this situation. Resolving these discrepancies can then result in a better determination of the jet energy
scale [73], a key quantity in any jet analysis at the LHC. Therefore, besides their purely theoretical interest, the results presented in this paper have the potential to serve as benchmarks to test and develop Monte Carlo event generators and, more concretely, to pin down the logarithmic accuracy of their parton shower in the spirit of [74–80].

Acknowledgments

We are grateful to our PanScales collaborators (Melissa van Beekveld, Mrinal Dasgupta, Frédéric Dreyer, Basem El-Menoufi, Silvia Ferrario Ravasio, Keith Hamilton, Jack Helliwell, Alexander Karlberg, Pier Monni, Gavin Salam, Ludovic Scyboz and Rob Verheyen) for discussions and comments on this manuscript. This work has been supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 788223, PanScales).

A Kinematic weights for dijets

Here we provide the kinematic weights entering into the formulae of section 3.2. We have extracted them from ref. [54] and s, t, u refer to the usual Mandelstam variables.

$$
\begin{align}
\omega_{12}^{gg \to gg}(s, t, u) &= \omega_{34}^{gg \to gg}(s, t, u) = \omega_{13}^{gg \to gg}(t, s, u) = \omega_{24}^{gg \to gg}(t, s, u) = \nonumber \\
&= \omega_{14}^{gg \to gg}(u, t, s) = \omega_{23}^{gg \to gg}(u, t, s) = C_A \left(\frac{2}{3} + \frac{H_{1d}(s, t, u)}{H_d(s, t, u)}\right), \tag{A.1a}
\end{align}
$$

$$
\omega_{12}^{gg \to q\bar{q}}(s, t, u) = \omega_{24}^{q\bar{q} \to qg}(t, s, u) = \omega_{34}^{q\bar{q} \to gg}(s, t, u) = 2N_c \left[ 1 - \left(1 - \frac{1}{N_c^2}\right) \frac{H_{1c}(s, t, u)}{H_c(s, t, u)} \right] \quad (A.1b)
$$

$$
\omega_{34}^{gg \to q\bar{q}}(s, t, u) = \omega_{12}^{q\bar{q} \to gg}(s, t, u) = \omega_{13}^{q\bar{q} \to qg}(t, s, u) = 2C_F \left(1 - 2\frac{H_{1c}(s, t, u)}{H_c(s, t, u)}\right), \quad (A.1c)
$$

$$
\begin{align}
\omega_{13}^{gg \to q\bar{q}}(s, t, u) &= \omega_{24}^{gg \to q\bar{q}}(s, t, u) = \omega_{13}^{q\bar{q} \to gg}(s, t, u) = \omega_{24}^{q\bar{q} \to gg}(s, t, u) = \nonumber \\
&= \omega_{12}^{q\bar{q} \to qg}(t, s, u) = \omega_{34}^{q\bar{q} \to qg}(t, s, u) = 2C_F \frac{H_{1c}(s, t, u) - H_{2c}(s, t, u)}{H_c(s, t, u)}, \tag{A.1d}
\end{align}
$$

$$
\begin{align}
\omega_{14}^{gg \to q\bar{q}}(s, t, u) &= \omega_{23}^{gg \to q\bar{q}}(s, t, u) = \omega_{14}^{q\bar{q} \to gg}(s, t, u) = \omega_{23}^{q\bar{q} \to gg}(s, t, u) \nonumber \\
&= \omega_{14}^{q\bar{q} \to qg}(t, s, u) = \omega_{23}^{q\bar{q} \to qg}(t, s, u) = 2C_F \frac{H_{1c}(s, t, u) + H_{2c}(s, t, u)}{H_c(s, t, u)}, \tag{A.1e}
\end{align}
$$

$$
\begin{align}
\omega_{12}^{qq' \to qq'}(s, t, u) &= \omega_{34}^{qq' \to qq'}(s, t, u) = \omega_{14}^{q\bar{q} \to q'\bar{q}'}(s, t, u) = \omega_{23}^{q\bar{q} \to q'\bar{q}'}(s, t, u) = \nonumber \\
&= \omega_{14}^{q\bar{q}' \to q\bar{q}'}(s, t, u) = \omega_{23}^{q\bar{q}' \to q\bar{q}'}(s, t, u) = 2/N_c, \tag{A.1f}
\end{align}
$$

$$
\begin{align}
\omega_{13}^{qq' \to qq'}(s, t, u) &= \omega_{24}^{qq' \to qq'}(s, t, u) = \omega_{12}^{q\bar{q} \to q'\bar{q}'}(s, t, u) = \omega_{34}^{q\bar{q} \to q'\bar{q}'}(s, t, u) = \nonumber \\
&= \omega_{13}^{q\bar{q}' \to q\bar{q}'}(s, t, u) = \omega_{24}^{q\bar{q}' \to q\bar{q}'}(s, t, u) = -1/N_c, \tag{A.1g}
\end{align}
$$

$$
\begin{align}
\omega_{14}^{qq' \to qq'}(s, t, u) &= \omega_{23}^{qq' \to qq'}(s, t, u) = \omega_{13}^{q\bar{q} \to q'\bar{q}'}(s, t, u) = \omega_{24}^{q\bar{q} \to q'\bar{q}'}(s, t, u) = \nonumber \\
&= \omega_{12}^{q\bar{q}' \to q\bar{q}'}(s, t, u) = \omega_{34}^{q\bar{q}' \to q\bar{q}'}(s, t, u) = 2C_F - 1/N_c, \tag{A.1h}
\end{align}
$$
$$ \begin{aligned} \omega_{12}^{qq\to qq}(s,t,u) &= \omega_{34}^{qq\to qq}(s,t,u) = \omega_{14}^{q\bar{q}\to q\bar{q}}(u,t,s) = \omega_{23}^{q\bar{q}\to q\bar{q}}(u,t,s) = \\ &= \frac{2 N_c H_a(s,t,u) + N_c H_a(s,u,t) - (N_c^2 + 1) H_{1b}(s,t,u)}{N_c^2 H_b(s,t,u)}, \end{aligned} \quad (A.1i) $$

$$ \begin{aligned} \omega_{13}^{qq\to qq}(s,t,u) &= \omega_{24}^{qq\to qq}(s,t,u) = \omega_{13}^{q\bar{q}\to q\bar{q}}(u,t,s) = \omega_{24}^{q\bar{q}\to q\bar{q}}(u,t,s) = \\ &= \frac{1 - N_c H_a(s,t,u) + N_c(N_c^2 - 2)H_a(s,u,t) + 2H_{1b}(s,t,u)}{N_c^2 H_b(s,t,u)}, \end{aligned} \quad (A.1j) $$

$$ \begin{aligned} \omega_{14}^{qq\to qq}(s,t,u) &= \omega_{23}^{qq\to qq}(s,t,u) = \omega_{12}^{q\bar{q}\to q\bar{q}}(u,t,s) = \omega_{34}^{q\bar{q}\to q\bar{q}}(u,t,s) = \\ &= \frac{1}{N_c^2} \frac{N_c(N_c^2 - 2)H_a(s,t,u) - N_c H_a(s,u,t) + 2H_{1b}(s,t,u)}{H_b(s,t,u)}. \end{aligned} \quad (A.1k) $$

In the previous equations, we have defined:

$$ H_a = \frac{s^2 + u^2}{t^2}, \qquad H_b = \frac{s^2 + u^2}{t^2} + \frac{s^2 + t^2}{u^2} - \frac{2}{N_c} \frac{s^2}{tu}, \quad (A.2a) $$

$$ H_c = (t^2 + u^2) \left[ \left(1 - \frac{1}{N_c^2}\right) \frac{1}{tu} - \frac{2}{s^2} \right], \qquad H_{1b} = \frac{s^2}{tu}, \quad (A.2b) $$

$$ H_{1c} = \frac{N_c}{4C_F} (t^2 + u^2) \left[ \left(1 - \frac{2}{N_c^2}\right) \frac{1}{tu} - \frac{2}{s^2} \right], \qquad H_d = 3 - \frac{ut}{s^2} - \frac{us}{t^2} - \frac{st}{u^2}, \quad (A.2c) $$

$$ H_{2c} = \frac{N_c}{4C_F} (t^2 - u^2) \left[ \frac{1}{tu} - \frac{2}{s^2} \right], \qquad H_{1d} = \frac{1}{6} \left( \frac{st}{u^2} + \frac{su}{t^2} - \frac{2tu}{s^2} - \frac{3s^2}{tu} + 3 \right). \quad (A.2d) $$

## B Extraction of the top-of-the-Lund-plane coefficients $D_{hme,R}^{(bc)}$

In this appendix we illustrate how we extract the coefficients $D_{hme,R}^{(bc)}$ from fixed-order Monte Carlo simulations. These coefficients are part of the NNDL resummation arising from hard matrix-element corrections at the top of the Lund plane, cf. section 3.3.1. We take the example of the inclusive jet sample studied in section 4.2 as a reference. An analogous procedure, not shown here, has been performed for the Z+jet sample.

The first step is to run a Monte Carlo simulation for the real part of the NLO corrections, i.e. Born-level 2 → 3 events in NLOJet++ in the present case. We build the Lund multiplicity at $O(\alpha_s)$ from these events split in flavour channels according to the procedure described at the end of section 3.3.1. The resulting distributions, differential in $\ln k_t$, are shown in figure 9a together with the analytic NDL expectations (dashed curves). The fact that the two converge at large $\ln(p_\perp R/k_t)$ validates our NDL calculation.

Next, we subtract the NDL contribution from the Monte Carlo distribution and build the cumulative distribution of the remainder. This is expected to asymptote to the total $O(\alpha_s)$ NNDL contribution. To isolate the piece originating from the top of the Lund plane, we further subtract the collinear endpoint correction, $D_{end}^{a\to bc}$ (cf. table 1). The resulting (cumulative) distributions are plotted in figure 9b. The results have been normalised by $\frac{\alpha_s}{2\pi}$ so as to directly correspond to $D_{hme,R}^{(bc)}$. We indeed see a convergence to a constant value which we extract for each flavour channel by taking the value at $\ln(p_\perp R/k_t) = 7$. The first uncertainty quoted in the legend of the plot corresponds to the Monte Carlo statistical
<figure><img src="image_11.png" /></figure>

(a) Differential multiplicity distributions, as obtained from NLOJet++ for the real NLO correction to our inclusive jet sample (solid lines). The dashed lines are the NDL expectations. The results are split in flavour channels.

<figure><img src="image_12.png" /></figure>

(b) Extraction of the $D_{h\text{me},R}^{(bc)}$ asymptotic values. These are obtained by integrating the results on the left and subtracting the collinear endpoint contribution ($D_{\text{end}}$ in table 1). The results are again split in flavour channels.

<caption>**Figure 9.** Illustration of the extraction of the $D_{h\text{me},R}^{(bc)}$ coefficients for our inclusive jet sample.</caption>

uncertainty, and the second number is an estimate of the systematic uncertainties obtained by extracting $D_{h\text{me},R}^{(bc)}$ at $\ln(p_{\perp}R/k_t) = 6.8$ and 7.2 and taking their difference. It is delicate to comment on the relative size of the various flavour channels as only the full sum in eq. (3.18) is independent of the prescription used to assign flavour.

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited. SCOAP³ supports the goals of the International Year of Basic Sciences for Sustainable Development.

## References

[1] R. Medves, A. Soto-Ontoso and G. Soyez, *Lund and Cambridge multiplicities for precision physics*, JHEP 10 (2022) 156 [arXiv:2205.02861] [INSPIRE].

[2] J. Gallicchio and M.D. Schwartz, *Quark and Gluon Tagging at the LHC*, Phys. Rev. Lett. 107 (2011) 172001 [arXiv:1106.3076] [INSPIRE].

[3] ATLAS collaboration, *Light-quark and gluon jet discrimination in pp collisions at $\sqrt{s} = 7$ TeV* with the ATLAS detector, Eur. Phys. J. C 74 (2014) 3023 [arXiv:1405.6583] [INSPIRE].

[4] C. Frye, A.J. Larkoski, J. Thaler and K. Zhou, *Casimir Meets Poisson: Improved Quark/Gluon Discrimination with Counting Observables*, JHEP 09 (2017) 083 [arXiv:1704.06266] [INSPIRE].

[5] F.A. Dreyer and H. Qu, *Jet tagging in the Lund plane with graph networks*, JHEP 03 (2021) 052 [arXiv:2012.08526] [INSPIRE].
[6] R. Kogler et al., Jet Substructure at the Large Hadron Collider: Experimental Review, Rev. Mod. Phys. 91 (2019) 045003 [arXiv:1803.06991] [INSPIRE].

[7] A.J. Larkoski, I. Moult and B. Nachman, Jet Substructure at the Large Hadron Collider: A Review of Recent Advances in Theory and Machine Learning, Phys. Rept. 841 (2020) 1 [arXiv:1709.04464] [INSPIRE].

[8] T. Lapsien, R. Kogler and J. Haller, A new tagger for hadronically decaying heavy particles at the LHC, Eur. Phys. J. C 76 (2016) 600 [arXiv:1606.04961] [INSPIRE].

[9] D. Krohn, J. Thaler and L.-T. Wang, Jets with Variable R, JHEP 06 (2009) 059 [arXiv:0903.0392] [INSPIRE].

[10] A. Tripathee et al., Jet Substructure Studies with CMS Open Data, Phys. Rev. D 96 (2017) 074003 [arXiv:1704.05842] [INSPIRE].

[11] A. Larkoski et al., Exposing the QCD Splitting Function with CMS Open Data, Phys. Rev. Lett. 119 (2017) 132003 [arXiv:1704.05066] [INSPIRE].

[12] ATLAS collaboration, Measurement of the Soft-Drop Jet Mass in pp Collisions at $\sqrt{s}$ = 13 TeV with the ATLAS Detector, Phys. Rev. Lett. 121 (2018) 092001 [arXiv:1711.08341] [INSPIRE].

[13] CMS collaboration, Measurement of the Splitting Function in pp and Pb-Pb Collisions at $\sqrt{s_{NN}}$ = 5.02 TeV, Phys. Rev. Lett. 120 (2018) 142302 [arXiv:1708.09429] [INSPIRE].

[14] CMS collaboration, Measurement of jet substructure observables in $t\bar{t}$ events from proton-proton collisions at $\sqrt{s}$ = 13 TeV, Phys. Rev. D 98 (2018) 092014 [arXiv:1808.07340] [INSPIRE].

[15] CMS collaboration, Measurement of the groomed jet mass in PbPb and pp collisions at $\sqrt{s_{NN}}$ = 5.02 TeV, JHEP 10 (2018) 161 [arXiv:1805.05145] [INSPIRE].

[16] ATLAS collaboration, Measurement of soft-drop jet observables in pp collisions with the ATLAS detector at $\sqrt{s}$ = 13 TeV, Phys. Rev. D 101 (2020) 052007 [arXiv:1912.09837] [INSPIRE].

[17] ATLAS collaboration, Measurement of jet-substructure observables in top quark, W boson and light jet production in proton-proton collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector, JHEP 08 (2019) 033 [arXiv:1903.02942] [INSPIRE].

[18] ATLAS collaboration, Measurement of the Lund Jet Plane Using Charged Particles in 13 TeV Proton-Proton Collisions with the ATLAS Detector, Phys. Rev. Lett. 124 (2020) 222002 [arXiv:2004.03540] [INSPIRE].

[19] STAR collaboration, Measurement of groomed jet substructure observables in $p+p$ collisions at $\sqrt{s}$ = 200 GeV with STAR, Phys. Lett. B 811 (2020) 135846 [arXiv:2003.02114] [INSPIRE].

[20] A LARGE ION COLLIDER EXPERIMENT and ALICE collaborations, Measurement of the groomed jet radius and momentum splitting fraction in pp and Pb-Pb collisions at $\sqrt{s_{NN}}$ = 5.02 TeV, Phys. Rev. Lett. 128 (2022) 102001 [arXiv:2107.12984] [INSPIRE].

[21] ALICE collaboration, Measurements of the groomed jet radius and momentum splitting fraction with the soft drop and dynamical grooming algorithms in pp collisions at $\sqrt{s}$ = 5.02 TeV, arXiv:2204.10246 [INSPIRE].

[22] ALICE collaboration, Measurements of groomed-jet substructure of charm jets tagged by D⁰ mesons in proton-proton collisions at $\sqrt{s}$ = 13 TeV, arXiv:2208.04857 [INSPIRE].
[23] M. Dasgupta, K. Khelifa-Kerfa, S. Marzani and M. Spannowsky, On jet mass distributions in $Z^+$ jet and dijet processes at the LHC, JHEP 10 (2012) 126 [arXiv:1207.1640] [INSPIRE].

[24] D. Napoletano and G. Soyez, Computing $N$-subjettiness for boosted jets, JHEP 12 (2018) 031 [arXiv:1809.04602] [INSPIRE].

[25] P. Cal, F. Ringer and W.J. Waalewijn, The jet shape at NLL', JHEP 05 (2019) 143 [arXiv:1901.06389] [INSPIRE].

[26] A. Lifson, G.P. Salam and G. Soyez, Calculating the primary Lund Jet Plane density, JHEP 10 (2020) 170 [arXiv:2007.06578] [INSPIRE].

[27] D. Reichelt et al., Phenomenology of jet angularities at the LHC, JHEP 03 (2022) 131 [arXiv:2112.09545] [INSPIRE].

[28] N. Ziani, K. Khelifa-Kerfa and Y. Delenda, Jet mass distribution in Higgs/vector boson + jet events at hadron colliders with $k_t$ clustering, Eur. Phys. J. C 81 (2021) 570 [arXiv:2104.11060] [INSPIRE].

[29] S. Caletti et al., Jet angularities in $Z^+$ jet production at the LHC, JHEP 07 (2021) 076 [arXiv:2104.06920] [INSPIRE].

[30] K. Lee, B. Meçaj and I. Moult, Conformal Colliders Meet the LHC, arXiv:2205.03414 [INSPIRE].

[31] E. Craft, K. Lee, B. Meçaj and I. Moult, Beautiful and Charming Energy Correlators, arXiv:2210.09311 [INSPIRE].

[32] S. Marzani, L. Schunk and G. Soyez, The jet mass distribution after Soft Drop, Eur. Phys. J. C 78 (2018) 96 [arXiv:1712.05105] [INSPIRE].

[33] Z.-B. Kang, K. Lee, X. Liu and F. Ringer, The groomed and ungroomed jet mass distribution for inclusive jet production at the LHC, JHEP 10 (2018) 137 [arXiv:1803.03645] [INSPIRE].

[34] Z.-B. Kang et al., The soft drop groomed jet radius at NLL, JHEP 02 (2020) 054 [arXiv:1908.01783] [INSPIRE].

[35] P. Cal, D. Neill, F. Ringer and W.J. Waalewijn, Calculating the angle between jet axes, JHEP 04 (2020) 211 [arXiv:1911.06840] [INSPIRE].

[36] D. Anderle et al., Groomed jet mass as a direct probe of collinear parton dynamics, Eur. Phys. J. C 80 (2020) 827 [arXiv:2007.10355] [INSPIRE].

[37] P. Cal, K. Lee, F. Ringer and W.J. Waalewijn, Jet energy drop, JHEP 11 (2020) 012 [arXiv:2007.12187] [INSPIRE].

[38] P. Cal, K. Lee, F. Ringer and W.J. Waalewijn, The soft drop momentum sharing fraction $\zg$ beyond leading-logarithmic accuracy, Phys. Lett. B 833 (2022) 137390 [arXiv:2106.04589] [INSPIRE].

[39] C. Frye, A.J. Larkoski, M.D. Schwartz and K. Yan, Factorization for groomed jet substructure beyond the next-to-leading logarithm, JHEP 07 (2016) 064 [arXiv:1603.09338] [INSPIRE].

[40] Z.-B. Kang, K. Lee, X. Liu and F. Ringer, Soft drop groomed jet angularities at the LHC, Phys. Lett. B 793 (2019) 41 [arXiv:1811.06983] [INSPIRE].

[41] M. Dasgupta, B.K. El-Menoufi and J. Helliwell, QCD resummation for groomed jet observables at NNLL + NLO, JHEP 01 (2023) 045 [arXiv:2211.03820] [INSPIRE].

[42] H.S. Hannesdottir, A. Pathak, M.D. Schwartz and I.W. Stewart, Prospects for strong coupling measurement at hadron colliders using soft-drop jet mass, arXiv:2210.04901 [INSPIRE].
[43] E. Gerwick, S. Schumann, B. Gripaios and B. Webber, QCD Jet Rates with the Inclusive Generalized kt Algorithms, JHEP 04 (2013) 089 [arXiv:1212.5235] [INSPIRE].

[44] B. Bhattachacherjee et al., Associated jet and subjet rates in light-quark and gluon jet discrimination, JHEP 04 (2015) 131 [arXiv:1501.04794] [INSPIRE].

[45] P. Caucal, A. Soto-Ontoso and A. Takacs, Dynamical Grooming meets LHC data, JHEP 07 (2021) 020 [arXiv:2103.06566] [INSPIRE].

[46] A.J. Larkoski, S. Marzani and J. Thaler, Sudakov Safety in Perturbative QCD, Phys. Rev. D 91 (2015) 111501 [arXiv:1502.01719] [INSPIRE].

[47] F.A. Dreyer, G.P. Salam and G. Soyez, The Lund Jet Plane, JHEP 12 (2018) 064 [arXiv:1807.04758] [INSPIRE].

[48] F.A. Dreyer, G. Soyez and A. Takacs, Quarks and gluons in the Lund plane, JHEP 08 (2022) 177 [arXiv:2112.09140] [INSPIRE].

[49] S. Catani, Y.L. Dokshitzer, F. Fiorani and B.R. Webber, Average number of jets in $e^{+}e^{-}$ annihilation, Nucl. Phys. B 377 (1992) 445 [INSPIRE].

[50] Y.L. Dokshitzer, G.D. Leder, S. Moretti and B.R. Webber, Better jet clustering algorithms, JHEP 08 (1997) 001 [hep-ph/9707323] [INSPIRE].

[51] M. Wobisch and T. Wengler, Hadronization corrections to jet cross-sections in deep inelastic scattering, in the proceedings of the Workshop on Monte Carlo Generators for HERA Physics (Plenary Starting Meeting), Hamburg Germany, April 27–30 April (1998), pp. 270–279 [hep-ph/9907280] [INSPIRE].

[52] A.J. Larkoski, D. Neill and J. Thaler, Jet Shapes with the Broadening Axis, JHEP 04 (2014) 017 [arXiv:1401.2158] [INSPIRE].

[53] M. Cacciari, G.P. Salam and G. Soyez, The anti-$k_t$ jet clustering algorithm, JHEP 04 (2008) 063 [arXiv:0802.1189] [INSPIRE].

[54] R.K. Ellis, G. Marchesini and B.R. Webber, Soft Radiation in Parton Parton Scattering, Nucl. Phys. B 286 (1987) 643 [Erratum ibid. 294 (1987) 1180] [INSPIRE].

[55] Y. Hatta, E. Iancu, A.H. Mueller and D.N. Triantafyllopoulos, Resumming double non-global logarithms in the evolution of a jet, JHEP 02 (2018) 075 [arXiv:1710.06722] [INSPIRE].

[56] M. Cacciari, G.P. Salam and G. Soyez, FastJet User Manual, Eur. Phys. J. C 72 (2012) 1896 [arXiv:1111.6097] [INSPIRE].

[57] J. Alwall et al., The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations, JHEP 07 (2014) 079 [arXiv:1405.0301] [INSPIRE].

[58] Z. Nagy, Next-to-leading order calculation of three jet observables in hadron hadron collision, Phys. Rev. D 68 (2003) 094002 [hep-ph/0307268] [INSPIRE].

[59] A. Banfi, G.P. Salam and G. Zanderighi, Accurate QCD predictions for heavy-quark jets at the Tevatron and LHC, JHEP 07 (2007) 026 [arXiv:0704.2999] [INSPIRE].

[60] PDF4LHC WORKING GROUP collaboration, The PDF4LHC21 combination of global PDF fits for the LHC Run III, J. Phys. G 49 (2022) 080501 [arXiv:2203.05506] [INSPIRE].

[61] S. Bright-Thonney and B. Nachman, Investigating the Topology Dependence of Quark and Gluon Jets, JHEP 03 (2019) 098 [arXiv:1810.05653] [INSPIRE].
[62] M. Cacciari et al., The t anti-t cross-section at 1.8 TeV and 1.96 TeV: A Study of the systematics due to parton densities and scale dependence, JHEP 04 (2004) 068 [hep-ph/0303085] [INSPIRE].

[63] C. Bierlich et al., A comprehensive guide to the physics and usage of PYTHIA 8.3, arXiv:2203.11601 [DOI:10.21468/SciPostPhysCodeb.8] [INSPIRE].

[64] P. Skands, S. Carrazza and J. Rojo, Tuning PYTHIA 8.1: the Monash 2013 Tune, Eur. Phys. J. C 74 (2014) 3024 [arXiv:1404.5630] [INSPIRE].

[65] R. Corke and T. Sjostrand, Interleaved Parton Showers and Tuning Prospects, JHEP 03 (2011) 032 [arXiv:1011.1759] [INSPIRE].

[66] ATLAS collaboration, ATLAS Pythia 8 tunes to 7 TeV data, ATL-PHYS-PUB-2014-021, CERN, Geneva (2014) [INSPIRE].

[67] J. Bellm et al., Herwig 7.0/Herwig++ 3.0 release note, Eur. Phys. J. C 76 (2016) 196 [arXiv:1512.01178] [INSPIRE].

[68] SHERPA collaboration, Event Generation with Sherpa 2.2, SciPost Phys. 7 (2019) 034 [arXiv:1905.09127] [INSPIRE].

[69] S. Abreu et al., Leading-color two-loop QCD corrections for three-jet production at hadron colliders, JHEP 07 (2021) 095 [arXiv:2102.13609] [INSPIRE].

[70] M. Czakon, A. Mitov and R. Poncelet, Next-to-Next-to-Leading Order Study of Three-Jet Production at the LHC, Phys. Rev. Lett. 127 (2021) 152001 [Erratum ibid. 129 (2022) 119901] [arXiv:2106.05331] [INSPIRE].

[71] ATLAS collaboration, Measurement of the charged-particle multiplicity inside jets from $\sqrt{s} = 8$ TeV pp collisions with the ATLAS detector, Eur. Phys. J. C 76 (2016) 322 [arXiv:1602.00988] [INSPIRE].

[72] CMS collaboration, Study of quark and gluon jet substructure in $Z$+jet and dijet events from pp collisions, JHEP 01 (2022) 188 [arXiv:2109.03340] [INSPIRE].

[73] ATLAS collaboration, Jet energy scale measurements and their systematic uncertainties in proton-proton collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector, Phys. Rev. D 96 (2017) 072002 [arXiv:1703.09665] [INSPIRE].

[74] M. Dasgupta et al., Parton showers beyond leading logarithmic accuracy, Phys. Rev. Lett. 125 (2020) 052002 [arXiv:2002.11114] [INSPIRE].

[75] K. Hamilton et al., Colour and logarithmic accuracy in final-state parton showers, JHEP 03 (2021) 041 [arXiv:2011.10054] [INSPIRE].

[76] A. Karlberg, G.P. Salam, L. Scyboz and R. Verheyen, Spin correlations in final-state parton showers and jet observables, Eur. Phys. J. C 81 (2021) 681 [arXiv:2103.16526] [INSPIRE].

[77] K. Hamilton et al., Soft spin correlations in final-state parton showers, JHEP 03 (2022) 193 [arXiv:2111.01161] [INSPIRE].

[78] M. van Beekveld et al., PanScales showers for hadron collisions: all-order validation, JHEP 11 (2022) 020 [arXiv:2207.09467] [INSPIRE].

[79] M. van Beekveld et al., PanScales parton showers for hadron collisions: formulation and fixed-order studies, JHEP 11 (2022) 019 [arXiv:2205.02237] [INSPIRE].

[80] F. Herren et al., A new approach to color-coherent parton evolution, arXiv:2208.06057 [INSPIRE].