Calculating the primary Lund Jet Plane density

The Lund-jet plane has recently been proposed as a powerful jet substructure tool with a broad range of applications. In this paper, we provide an all-order single logarithmic calculation of the primary Lund-plane density in Quantum Chromodynamics, including contributions from the running of the coupling, collinear effects for the leading parton, and soft logarithms that account for large-angle and clustering effects. We also identify a new source of clustering logarithms close to the boundary of the jet, deferring their resummation to future work. We then match our all-order results to exact next-to-leading order predictions. For phenomenological applications, we supplement our perturbative calculation with a Monte Carlo estimate of non-perturbative corrections. The precision of our final predictions for the Lund-plane density is 5-7% at high transverse momenta, worsening to about 20% at the lower edge of the perturbative region, corresponding to transverse momenta of about 5 GeV. We compare our results to a recent measurement by the ATLAS collaboration at the Large-Hadron Collider, revealing good agreement across the perturbative domain, i.e. down to about 5 GeV.

A huge variety of observables has been explored (e.g. [24][25][26][27][28][29][30][31][32][33]) for studying jet substructure, supplemented in recent years by a range of machine-learning approaches (e.g. [32,[34][35][36][37][38][39][40][41][42][43][44][45]). With such a diverse range of observables, it has become challenging to obtain a detailed understanding of the specific jet features probed by each one. At the same time approaches have emerged in which one designs an infinite set of observables which, taken as a whole, can encode complete information about a jet. Specific examples are energy-flow polynomials [32] and the proposal [33] (see also [14]) to determine a full Lund diagram [46] for each jet. As well as encoding complete information about the radiation in a jet, both of these approaches provide observables that can be directly measured and that also perform well as inputs to machine-learning. Here we concentrate on Lund diagrams.
Lund diagrams [46] are two-dimensional representations of the phase-space for radiation in jets, which have long been used to help understand Monte-Carlo event generators and all-order logarithmic resummations in QCD. The phase-space for a single emission involves three degrees of freedom, and Lund diagrams highlight the logarithmic distribution of two of those degrees of freedom, typically chosen to be the emission's transverse momentum (k t ) and angle (∆). In leading-order QCD, the logarithms of both variables are uniformly distributed.
A core idea introduced in Refs. [33] and [14] and briefly reviewed in section 2, is to use a Cambridge/Aachen declustering sequence to represent a jet's internal structure as a series of points in the two-dimensional Lund plane, with the option of concentrating on "primary" emissions, those that can be viewed as emitted by the jet's main hard prong. The location of a given point immediately indicates whether it is in a perturbative or nonperturbative region, whether it is mainly final-state radiation or a mix with initial-state radiation, underlying event, etc. The set of points obtained for a single jet can be used as an input to multi-variate tagging methods [33], or can be used to construct other specialised observables [47]. Given an ensemble of many jets, one can also determine the average density of points in each region of the (primary) Lund plane, ρ(∆, k t ). This average density is of interest in fundamental measurements of QCD radiation [48][49][50], both in perturbative and non-perturbative regions, and in studies of modifications of jet structure in heavy-ion collisions [14].
The purpose of this article is to carry out a baseline calculation of the all-order perturbative structure of the primary Lund-plane density, identifying the key physical aspects that are relevant for understanding the density, and providing a prediction that accounts for all single-logarithmic corrections α n s ln m ∆ ln n−m k t multiplying the leading-order, O (α s ), result for the density. The relevant contributions are discussed in section 3 and include running-coupling effects, collinear flavour-changing effects and various effects of soft radiation at commensurate angles (which we compute only in the large-N c approximation). 1 In section 4 we match the all-order results to a next-to-leading order (NLO) calculation using the NLOJet++ program [51], in section 5 we address the question of non-perturbative corrections and in section 6 we combine the different results into a set of final predictions that we compare to recent experimental measurements from the ATLAS collaboration [48]. 2 The primary Lund plane density and basic setups Let us assume we have a jet with transverse momentum p ⊥ , obtained from a given jet algorithm such as the anti-k t algorithm [52]. We first re-cluster the constituents of the jet using the Cambridge/Aachen (C/A) algorithm [53,54] as often used in jet substructure techniques. We then iteratively repeat the following steps, starting with j defined as the full (re-clustered) jet: 1. Undo the last step of clustering: j → j 1 + j 2 , taking j 1 to be the harder branch, i.e. p ⊥1 > p ⊥2 .
2. Record the properties of the branching T ≡ {k t , ∆, z, ...} defined as where y and φ denote the rapidity and azimuthal angle of a particle, specifically y = 1 2 ln E+pz E−pz .

(2.2)
Additional variables can be added to each tuple, for example an azimuthal angle. One can also choose to follow softer branchings at each step, which would lead to exploration of secondary, tertiary, etc. Lund planes. Neither of these aspects is relevant for the discussion presented here. The primary Lund plane density is then defined as the density of emissions in the (logarithmic) ∆, k t plane: 2 ρ(∆, k t ) = 1 N jets Alternatively, one can introduce a primary Lund plane density in the ∆, z plane: ρ(∆, z) = 1 N jets dn emissions d ln 1/∆ d ln 1/z , (2.4) as measured by the ATLAS collaboration.
Note that integrating the primary Lund plane density over ln 1/∆ and ln k t (or ln 1/z) gives the average number of primary emissions per jet. If the integration is performed with an additional Soft-Drop [28,29] condition, one obtains the Iterated Soft-Drop multiplicity [11] (see also chapter III of [55]).
In practice, we will focus on two kinematic configurations: • High-p ⊥ setup. This is close to the original proposal from [33]. We cluster jets with the anti-k t algorithm with R = 1, keep all jets with p ⊥ ≥ 2 TeV. The primary Lundplane density ρ(∆, k t ) is then reconstructed according to the procedure described above.
• ATLAS setup. This is similar to the ATLAS measurement presented in [48]. Jets are reconstructed with the anti-k t algorithm with a radius R = 0.4. The two largestp ⊥ jets with |η| < 2.1 are kept (with η the pseudo-rapidity, defined as η = − ln tan θ 2 ). One then imposes that the leading jet has a p ⊥ of at least 675 GeV and that the p ⊥ of the second jet is at least 2 3 of the p ⊥ of the leading jet. For each of the two jets, we construct the Lund planeρ(∆, z) as follows: we take all the particles within a radius R = 0.4 of the jet axis, recluster them with the C/A algorithm with R = 0.4 and apply the de-clustering procedure highlighted above.
In practice the ATLAS measurement only includes charged tracks with p ⊥ above 500 MeV within a distance (η − η jet ) 2 + (φ − φ jet ) 2 = 0.4 of the jet axis 3 in their reconstruction of the Lund plane. We will treat this as a non-perturbative correction (going from a full-particle measurement to a measurement based on tracks above 500 MeV). We note that the use of charged tracks makes this measurement collinear unsafe since, for example, arbitrarily collinear branchings can affect the relative fractions of charged and neutral particles in each branch, and consequently the definition of the harder branch in the de-clustering procedure. Numerically, this effect is small, as we shall verify later.
In all cases, the initial jet clustering is done using FastJet [58,59] and the Lund plane is constructed using the code available with fastjet-contrib [60].

All-order calculation
The average Lund plane density measures an effective intensity of radiation per unit logarithm of k t and of angle. As such at LO, in the simultaneously soft and collinear limit, i.e. 3 This distance uses pseudo-rapidity η instead of rapidity y. The two variables are identical for massless objects. For individual experimental objects, it is the pseudorapidity that is measured. However for any object that is massive, rapidity is to be strongly favoured [56,57] because rapidity differences are invariant under longitudinal boosts, while pseudorapidity differences are not. Every stage of jet clustering creates massive objects, even starting from massless ones. The pseudorapidity of a jet displays various pathologies: for example, a jet consisting of two massless particles with identical pseudorapidities η1 = η2 has a rapidity yjet = η1 = η2, while the jet's pseudorapidity is different ηjet = η1, by an amount that depends non-trivially on the kinematics of the jet. Therefore even if the inputs to the initial jet clustering are selected based on their pseudo-rapidity, we recommend using rapidity for all subsequent operations.
-4 -away from the large-angle and the collinear edges of the plane, it is given by where C i is the Casimir of the hard parton of flavour i initiating the jet, C i = C A for a gluon-initiated jet and C F for a quark-initiated jet. Beyond leading order, each additional factor of α s can be associated with up to one logarithm of either ∆ or of p ⊥ /k t . As a result, at any given order, say α n+1 s , the logarithmically dominant terms have the structure α n+1 s ln m ∆ ln n−m p ⊥ kt with 0 ≤ m ≤ n. Our goal is to calculate this complete set of single-logarithmic contributions to ρ(∆, k t ), i.e. for all n and m, including the full (non-logarithmic) ∆ dependence for terms with m = 0 and the full k t dependence for terms with m = n. In the case ofρ(∆, z), we will equivalently aim to account for all terms α n+1 s ln m ∆ ln n−m z. The logarithms have several physical origins. These are (i) running coupling corrections, enhanced by logarithms of the transverse momentum k t ; (ii) hard-collinear logarithms of the emission angle ∆ which can induce flavour-changing effects and affect the behaviour of ρ close to the z = 1 line; (iii) soft emissions at large angles enhanced by the logarithm of either k t or the emission energy fraction z; and (iv) Cambridge/Aachen clustering effects for emissions with commensurate angles, enhanced by logarithms of k t or z. Each of these effects is discussed separately in the following sections.

Running coupling corrections
This is by far the simplest correction: the scale of the running coupling is simply set by the transverse momentum of the emission. We therefore have We use the 2-loop running coupling in the CMW scheme [61]: The reference α s ≡ α s (p ⊥ R) is taken at the scale p ⊥ R with p ⊥ the transverse momentum of the jet and R the jet radius. We use n f = 5 flavours for k t ≥ m b = 4.78 GeV, n f = 4 for m b > k t ≥ m c = 1.67 GeV and n f = 3 below. Furthermore, we freeze the coupling at k t = 1 GeV. At our accuracy, it would have been sufficient to use the 1-loop running coupling (without the CMW scheme term, i.e. K). We have instead used the 2-loop running for two main reasons. Firstly, several of our result depend on the structure of the hard events -5 -under consideration (dijet events in our case). These events will be obtained using the NLOJet++ program [51] with underlying PDF sets that use at least a 2-loop running. It is therefore more coherent to use a 2-loop running also in the computation of ρ. Secondly, the running coupling is numerically the largest of the logarithmically-enhanced contributions. Including the 2-loop corrections therefore makes sense from a purely phenomenological perspective.

Hard-collinear effects
Eq. (3.2) assumes that the primary branch followed by the declustering procedure keeps the flavour of the initial parton. In this case, all the emissions for a quark-initiated jet (gluon-initiated jet) come with a C F (C A ) colour factor.
In practice, however, hard and collinear branchings have two effects at single-logarithmic accuracy: (i) they can change the flavour of the harder branch via either a q → qg splitting where the daughter gluon carries more than half the parent quark's momentum, or from a g → qq splitting; and (ii) successive collinear branchings can reduce the transverse momentum of the leading parton thereby creating a difference between the z = 1 and the k t = p ⊥ ∆ lines in the Lund plane. These effects are of the form α n+1 s ln n 1/∆, associated with a series of emissions strongly ordered in angle and without soft enhancement.
For an initial parton of flavour i, we use p x, j|i, t coll (∆; ∆ 0 , µ) , (3.5) to denote the probability of having a leading parton of flavour j, carrying a longitudinal fraction x, when the primary Lund declustering procedure has reached an angular scale ∆. The dependence on ∆ is encoded through a "collinear evolution time" where I α is the integration of the coupling between two transverse momentum scales: with α s ≡ α s (µ), λ = 1 + 2α s β 0 ln(k t /µ) and λ 0 = 1 + 2α s β 0 ln(k t0 /µ). The expression on the second line corresponds to a fixed number of flavours. The parameter ∆ 0 in (3.6) is the large-angle scale at which the collinear evolution starts and can be varied to get an estimate of the uncertainties associated with the resummation of the collinear logarithms. The evolution of the probability densities p(x, j|i, t) is given by The kernels P (R) jk (z) and P (V ) jk (z) correspond to real and virtual emissions respectively and are straightforwardly obtained from the DGLAP splitting functions, imposing that the leading parton is defined following the larger-p ⊥ branch: where the P ij (z) are the normal full (real) DGLAP splitting functions, including appropriate symmetry factors. If the leading parton carries a longitudinal momentum xp ⊥ at an angle ∆, the splitting variables (cf. (2.1)) are related through z = kt xp ⊥ ∆ . Therefore, for a jet initiated by a hard parton of flavour i, the primary Lund-plane density including collinear effects takes the form (3.11) This expression includes the effects of the flavour changes and the distribution of the longitudinal momentum fraction of the leading parton. The factor z 2C j P j (z), with P j (z) = k P (R) kj (1 − z), accounts for the fact that close to the z = 1 boundary, one should use the full splitting function instead of its soft limit.
In practice, an exact analytic calculation of p(x, j|i, t coll ) is not possible. It is however straightforward to obtain it numerically using the approach of Ref. [62]. 4 A sample of the resulting distributions is shown in Fig. 1 for both quark-initiated and gluon-initiated jets. The distributions progressively shift towards smaller values of x as expected.
From Eq. (3.8) one can also deduce a few analytic properties of p(x, j|i, t coll ). In particular, if one takes the 0 th and 1 st moments of (3.8) one obtains respectively an evolution equation for the average fraction of quarks and gluons and the average longitudinal momentum of the leading parton. We therefore define   The red (blue) curves correspond to a quark (gluon) initial parton. The solid (dashed) curves correspond to a quark (gluon) leading parton at time t coll . The top axis shows the angle corresponding to t coll for a 2-TeV jet. The t coll axis extends significantly beyond the typical perturbative region so as to help illustrate the asymptotic trends.
Plots of f (j|i, t) andx(j|i, t) are shown in Fig. 2 from which we can make several observations. First, as t coll → ∞, i.e. ∆ → 0 (modulo Landau-pole complications), the quark and gluon fractions tend to constants that are independent of the initial flavour of -8 -the jet. From (A.1) one finds f (q|any, t coll → ∞) = s g s q + s g , f (g|any, t coll → ∞) = s q s q + s g , (3.13) with s q = C F (2 ln 2 − 5 8 ) and s g = 1 3 n f . This means about 62.1% quarks and 37.9% gluons for n f = 5, in agreement with Fig. 2a. Furthermore, for t coll → 0, we find where the w gq and w qg coefficients are given analytically in Appendix A. The results in Eq. (3.14) correspond to the average harder-parton momentum fractions after a single q → qg splitting (with the gluon the harder particle) or a single g → qq splitting. The numerical values are in agreement with Fig. 2b.

Soft emissions at large or commensurate angles
The average primary Lund density is subject to several classes of effect associated with the non-trivial characteristics of soft radiation. At O (α s ), when one goes beyond the collinear limit of Eq. (3.1) and considers radiation at angles ∆ ∼ 1, soft gluon radiation is a coherent sum from all hard coloured partons in the event rather than just a single parton. At higher orders, two further effects arise: the presence of a first soft gluon contributes to the radiation of a subsequent second soft gluon at commensurate angles, and so forth for higher numbers of gluons, contributing effects similar to non-global logarithms; and in the presence of two or more gluons at commensurate angles, one must account for the way in which jet clustering determines whether a given gluon is classified as a primary or a secondary Lund emission. These two effects are present both for large and (perhaps more surprisingly) small ∆. In this section we will consider all of these effects, using the large-N C limit so as to retain simple colour algebra. After considering how we decompose events into separate colour flows, we shall in section 3.3.1 examine the impact of different colour flows on the large-angle part of the Lund plane at O (α s ). Then in section 3.3.2 we shall consider the case of double (energy-ordered) soft gluon emission and derive the structure of the nonglobal and clustering logarithms in the small-angle limit at order α 2 s . Finally in section 3.3.3 we will discuss how we generalise these results to resum all sources of single soft logarithms to all orders.
We start by recalling that in the large-N C limit, a given Born-level process can be expressed as a sum over several partonic channels where each of them is a weighted sum of different colour flows: (3. 15) In this context, for each colour flow, one can view the Born-level process as a superposition of colour dipoles. Let us consider a dijet process with two incoming partons, p 1,2 , and two outgoing partons, a jet parton p j and a recoiling parton p r . The relative weights of the different partonic channels can be obtained from the 2 → 2 squared matrix elements, e.g. using NLOJet++. Then, a q in q in → q out q out channel would have, in the large-N c limit, a single colour flow with weight w = 1 corresponding to dipoles [(q in q out ) + (q in q out )]. Similarly, a g in q in → g out q out channel would have two colour flows: with s, t, u the usual Mandelstam variables. The complete set of colour flows, dipole superpositions and weights can, for example, be deduced from Ref. [63].
In the next sections, it will be helpful to separate ρ soft in different contributions according to the Born-level flavour of the jet: (3.17) where f i denotes the relative fraction of quark and gluon jets. This separation in jet flavours can be straightforwardly obtained from (3.15). We can expand ρ soft,i as a series in α s : In the first subsection below, we will show that ρ (1) soft,i deviates from (3.1) by corrections that are power-suppressed in ∆. At our single-logarithmic accuracy, we have These soft logarithms are either due to non-global configurations or to clustering logarithms associated with the Cambridge/Aachen reclustering used to construct the primary Lundplane density. It is interesting to note that these clustering logarithms are present at arbitrarily small angles, which is somehow uncommon (though addressed also in [64]). We show how they appear at order α 2 s in section 3.3.2 and provide an all-order resummation in section 3.3.3.

Soft emissions at large angles: fixed-order study
For definiteness, let us consider the case of two incoming partons, 1,2 , and two outgoing partons, j,r , with the following kinematics: 5 where Q = p ⊥ cosh y jet (with the jet transverse momentum p ⊥ ≡ ⊥ , while its rapidity is equal to y jet ). We consider the jet of radius R around j , while r corresponds to the recoiling hard jet. In the large-N c approximation we have to consider soft gluon emission from any of 6 possible colour dipoles: one incoming-incoming, two jet-incoming, two recoil-incoming and one jet-recoil. For a given dipole with legs a and b , the contribution from a soft emission of momentum k µ ≡ k ⊥ (cos φ, sin φ, sinh y, cosh y) takes the form For each dipole configuration, we can set y = y jet +∆R cos ψ and φ = ∆R sin ψ. The k ⊥ and ∆R integrations can then be evaluated trivially, leaving the integration over ψ. This integration usually cannot be computed exactly so we instead perform a series expansion in ∆ 2 : soft,1j and ρ (1) soft,2r is obtained from ρ (1) soft,1r via the replacement y jet → −y jet .
In the collinear limit, ∆ 1, ρ soft,2j , ρ soft,1j and ρ (1) soft,jr all tend to a constant, with corrections taking the form of power corrections in ∆ 2 , as expected. The other dipoles are suppressed by a factor ∆ 2 .
While running-coupling and collinear flavour-changing effects only depend on the flavour of the jet, the corrections due to soft emissions at large angles involve the structure of the whole event. The relative weight of each dipole depends on the channel and colour flow under consideration, cf. (3.15).

Soft emissions and Cambridge/Aachen clustering: fixed-order study
Say we want to extend the calculation from section 3.3.1 to order α 2 s . The same calculation as above would have to be repeated with two soft emissions, k 1 and k 2 , strongly ordered in energy (k ⊥1 k ⊥2 ). Measuring the emission k 2 and integrating out k 1 yields a contribution to the primary Lund plane density of the form In the limit ∆ → 0 the prefactor is simple, and we calculate it here.
-11 -For concreteness, we illustrate the case of a quark-induced jet. We denote by θ i the angle between k i and the quark and by θ 12 the angle between k 1 and k 2 and work in a limit where all angles are small. In contrast to section 3.3.1, we now use a frame where the jet is perpendicular to the beam. In conjunction with the small-angle limit, this ensures that angles and rapidity-azimuth distances are equivalent, as are energies and transverse momenta (with respect to the beam).
Let us first consider three simple nested-collinear limits. When θ 1 θ 2 , gluon k 2 is emitted with colour factor C F and declustered as a primary emission, i.e. the LO 2α s C R /π emission intensity for gluon k 2 is unaffected by the presence of gluon k 1 . The situation is similar when θ 1 θ 2 . When θ 12 θ 1 , gluon k 2 is emitted with colour factor C A and declustered as a secondary emission, i.e. on gluon k 1 's Lund leaf. The only non-trivial situation is when the angles θ 1 , θ 2 and θ 12 are commensurate, θ 1 ∼ θ 2 ∼ θ 12 1. In this region, one needs to account for the non-trivial matrix element for the emission of two gluons at commensurate angles and for the effects of the Cambridge/Aachen clustering used to construct the Lund plane. Together, these induce an O(α 2 s ln k t /(p ⊥ ∆)) correction to the 2α s C R /π behaviour, where the logarithm is associated with the integral over the transverse momentum of gluon k 1 .
This contribution to the primary Lund density for a quark-induced jet can be written where k ⊥i is the transverse momentum of k i relative to the beam. In this expression the first (second) term of the curly bracket corresponds to a real (virtual) emission k 1 . For the real emission, we have two factors: the first square bracket corresponds to the matrix element for the emission of the soft gluon k 2 and the second square bracket imposes that the gluon k 2 is reconstructed as a primary emission, i.e. is not clustered with the emission k 1 . Eq. (3.24) genuinely encodes clustering effects. If we first consider the C 2 F contribution, naively associated with two emissions from the hard quark, the virtual term partially cancels the real contribution, leaving a negative contribution with a factor Θ(θ 12 < θ 1 )Θ(θ 12 < θ 2 ), i.e. where emission k 2 clusters with emission k 1 . Obviously, this contribution disappears in the collinear limit θ 2 θ 1 , θ 12 as expected. Focusing now on the C F C A term, naively associated with secondary k 2 emission, the only contribution comes from the situation where the emission is not clustered with its emitter, k 1 , which vanishes in the collinear limit θ 21 θ 1 . The k ⊥1 integration in Eq. (3.24) has a logarithmic enhancement from strong energy ordering k ⊥2 k ⊥1 p ⊥ , leading to a contribution proportional to α 2 s ln kt p ⊥ ∆ , as anticipated in Eq. (3.23), and no collinear divergence. The integrand is suppressed in the limits θ 1 → 0 and θ 1 → ∞ and only receives a contribution from θ 1 ∼ θ 2 . 6 6 A consequence of this is that, at our logarithmic accuracy, we can safely set the upper bound on the -12 - The integration over k ⊥1 , k ⊥2 , θ 2 and one of the azimuthal angles (ϕ 1 or ϕ 2 ) can be trivially performed, leaving an integration over θ 1 and an azimuthal angle ϕ. One finds 7 The calculation for a gluon jet can be obtained by replacing C F → C A in Eq. (3.24). That replacement carries through directly to Eq. (3.25), giving Thus for a purely gluonic theory, the energy-ordered double-soft emission pattern and the C/A clustering combine in such a way that there is no α 2 s ln k t /(p ⊥ ∆) correction to the Lund density when ∆ is small.
The above results are valid at small angles. Two additional classes of effect arise at large angles. Firstly, the clustering effects become sensitive to the coherent structure of the radiation from the complete hard event. This relates to the discussion in section 3.3.1. Secondly, if one identifies the jet with the anti-k t algorithm and reclusters its constituents with the C/A algorithm, there is an interplay between the two clusterings. This leads to another source of logarithmic enhancement The ln( R R−∆ ) structure appears when a first emission, close to but outside the anti-k t jet boundary, splits collinearly such that one of its offspring is inside the boundary. 8 The allorder resummation of these boundary logarithms is beyond the scope of this paper. They are however briefly discussed in Appendix B. Note that if the original jet is identified with the C/A algorithm, these boundary logarithms are absent.

Soft emissions: all-order treatment
The treatment of soft single logarithms to all orders requires us to consider configurations with arbitrarily many energy-ordered gluons at commensurate angles [66], for which analytic approaches exist only in specific limits [67]. The technique we adopt is similar to that originally proposed for the resummation of non-global logarithms in [66] and the related clustering logarithms [68,69]. We rely on a large-N C approximation (with subleading colour corrections up to order α 2 s in the collinear limit). Techniques that exist to resum non-global logarithms at full N C are so far applicable only for a limited set of observables [70,71].
Compared to the typical treatment of non-global and clustering logarithms we have one extra difficulty and one simplification. The difficulty has the following origin. Since θ1 integration to infinity. 7 The numerical pre-factor is analytically found to be −iπ . 8 It is related to the ∆η ln ∆η term observed in Eq. (3.13) of [65].
-13 -non-global logarithms stem from emissions at commensurate angles one can usually impose an angular cut-off at an angle θ min that is small compared to the physical angle one probes (π/2, a jet radius, the rapidity width of a slice, ...). This helps limit the particle multiplicity and associated computational cost of the calculations. In the case of the primary Lund plane, we instead have to probe a large range of angles, meaning that potential angular cut-offs have to be taken small/large enough to cover this extended phase-space, resulting in increased computational demands. 9 The simplification relative to normal non-global logarithm calculations relates to the fact that the Lund plane density doesn't involve any Sudakov suppression (unlike say a hemisphere mass). That Sudakov suppression, specifically the part associated with primary emissions, can lead to low computational efficiency unless dedicated subtraction techniques are applies. In the case of the Lund plane density one can simply generate all the emissions without separating primary emissions from the other ones.
The basic approach to simulating soft emissions is to directly order them in energy, as done in [66]. In this case we just need to generate the full angular structure of the emissions and only retain their energy ordering. If gluon k j is emitted after k i , then it has a much smaller energy. In this case clustering with the anti-k t algorithm is equivalent to keeping the particles in a radius R around the hard parton. For C/A clustering, all the necessary distances are available from the angular structure of the event and the recombination of two particles is equivalent to replacing them with the harder one.
For our ultimate primary Lund plane predictions, we want to focus on the phasespace above a certain k t cut, below which the non-perturbative effects dominate. For a given minimum relative transverse momentum k t,min the minimum accessible angle is ∆ min = k t,min /p ⊥ . We therefore have to take an angular cut-off for the event simulations that is sufficiently smaller that ∆ min . If we generate emissions down to an energy E min , many of these emissions will have a k t much smaller than k t,min . For example for ∆ = ∆ min , emissions would be generated down to k t = E min ∆ min k t,min . This is not a problem per se, except for the fact that this approach generates many more emissions than absolutely necessary, which ends up being computationally challenging.
For this reason, we have adopted a different approach, more traditional in partonshower event generators, namely we generate the emissions ordered in k t . Say we work with a fixed coupling. The event is described as a collection of dipoles, each with a hard scale corresponding to their invariant mass Q. For the initial condition, we decompose the Born-level event as a sum over all possible dipole configurations. At any given stage of the event generation, corresponding to a given scale k t = k ti , we should be able to generate the next emission at a scale k t,i+1 < k ti . If a dipole (p 1 , p 2 ) of invariant mass Q 12 splits, this is done by first generating k t,i+1 according to the following Sudakov factor (which includes 9 Similar non-global and clustering logarithms have been studied down to small angles in Ref. [64] for the study of the SoftDrop grooming radius at NLL accuracy. There, the authors relied on the fact that the behaviour of the clustering logarithms becomes independent of ∆ at small-enough ∆. In this paper, we decided not to rely on this behaviour so as to also reach a good level of numerical precision for the approach to this asymptotic regime within our single-logarithmic accuracy, and to do so over the relatively large energy range needed to cover the full Lund plane. -14 -collinear radiation at each end of the dipole): One then decides the rapidity η of the emission, uniformly distributed between ln k t,i+1 /Q 12 and ln Q 12 /k t,i+1 as well as an azimuthal angle φ. The 4-momentum of the new emission is thus reconstructed as and n 1,2 two unit vectors orthogonal to p 1 and p 2 . Note that when a dipole (p 1 , p 2 ) splits into two new dipoles (p 1 , k i+1 ), (k i+1 , p 2 ), the energy scales of the two new dipoles can be straightforwardly obtained using In order to reduce the inclusion of uncontrolled corrections beyond our intended resummation of single logarithms, we can perform another simplification. Recall that we are aiming to resum the energy logarithms due to clustering effects. Instead of computing the energy of each particle explicitly from its 4-momentum, we can directly project this momentum along the direction of the initial dipole. Let us denote by (p,p) the hard-scattering momenta of a given initial colour dipole. For each emission, we want to find the contributions z andz of its momentum fractions along the p andp directions respectively. For the emission of a new particle k from a (p 1 , p 2 ) dipole, we have The new emission k therefore has a projection z k ,z k along the p,p directions given by where we have used Eq. (3.30) and replaced the sum over the two contributions by its maximum at our accuracy. Iterating the above procedure for emissions ordered in k t produces an event where each particle has a 4-momentum as well as longitudinal fractions z andz along their initial (p,p) dipole. To reconstruct the primary Lund plane density we then proceed as follows: the anti-k t jet of radius R is made of all the emissions within a radius R of an initial hard parton. These particles can then be clustered using the C/A jet algorithm. This clustering uses the exact 4-momenta of the jets and a winner-takes-all-like [72] recombination scheme where the recombined particle is taken as the one of the two recombining particles with the largest z ( -) momentum along the jet direction. 10 When we reconstruct the primary Lund plane density, the ∆ coordinate is again taken from the exact angular kinematics, the z variable from the event z ( -) and hence k t is obtained as p t,hard z∆ with p t,hard the jet transverse momentum (w.r.t. the colliding beams) of the initial hard parton.
The above discussion is strictly speaking valid only in a fixed-coupling approximation. To account for running-coupling effects, we use the following procedure. For given Bornlevel kinematics (see e.g. (3.20)) and a given colour flow corresponding to a given initial set of dipoles, we generate a Monte Carlo event using a fixed α s . 11 Following the procedure outlined in the previous paragraph we obtain the coordinates ∆ and z of the primary Lund declusterings. From the z coordinate, we then determine an emission "time" t soft,fc defined as t soft,fc = α s ln 1/z. This procedure yields a resummed density ρ soft (∆, t soft,fc ). To include running-coupling corrections at a given k t and ∆, we simply use ρ soft (∆, t soft ) with a t soft defined to include running-coupling effects: Additionally, this approach can be straightforwardly extended to generate results at fixedorder. This will be useful to compare to the results derived in section 3.3.1 and 3.3.2 as well as for matching with exact fixed-order results in section 4. As a validation of our numerical approach, we first compare the output of the numerical approach to the analytic results for soft gluon radiation at fixed order. The predictions for soft radiation at large angle at O(α s ), obtained in section 3.3.1, are compared to our numerical results in Fig. 3a. The comparison is done assuming the large-N c limit and is independent of ln k t (modulo the soft approximation of the kinematic constraint, k t /(p ⊥ ∆) ≤ 1). The figure shows excellent agreement with the analytic results from Eq. (3.22).
With Fig. 3b, we study the numerical results for C/A clustering effects at O(α 2 s ), Eqs. (3.25) and (3.26). For a gluon jet, we need to consider two dipoles. Since our calculation is done in the collinear limit, we have considered a range of small values of ∆. The linear rise with ln z, with the expected analytic coefficient, is clearly visible for quark jets, together with no effects at this order for gluon jets.
All-order results are shown in Fig. 4 for two different regions in angle. We see that apart from the region of very small z (large t soft ), the resummation has a relatively small 10 The usage of the z ( -) momentum fractions guarantees that in the collinear limit only the logarithmicallyenhanced contributions are kept. At large angles, and, in particular, for initial dipoles which do not involve the jet momentum, one could generate subleading corrections as well. In practice however, our fixed-order tests indicate that these subleading corrections are very small, if present at all. 11 In practice, we take the coupling at the scale p ⊥ R.  For the quark configuration we have used a (1j) dipole connecting the incoming particle 1 and the outgoing "jet" j , while for the gluon configuration, we have considered the superposition of 2 dipoles (1j, 2j), i.e. one dipole connecting the incoming 1 with j and a second dipole connecting the other incoming parton 2 with j (cf. (3.20)). effect compared to the NLO result. A feature that is particularly intriguing is that in the collinear region, ∆ 1, Fig. 4b, the result appears to be independent of t soft . Recall that for gluon jets, at order α 2 s the soft logarithmic term was identically zero for small ∆, Eq. (3.26). Fig. 4b leads us to wonder whether the soft single logarithmic terms remain zero for gluon jets at all orders, or whether they are non-zero but simply too small to observe in our calculation. Note however that at large angles, Fig. 4a, there is a clear t soft dependence both at α 2 s and beyond, i.e. the soft single logarithmic coefficients are non-zero.

Full resummed result
Our final resummed predictions include all the effects discussed in this section: the running of the strong coupling, collinear effects -flavour changes, splitting functions and the momentum of the leading parton -as well as soft-gluon emissions to all orders including large-angle contributions and clustering effects for emissions at commensurate angles: In this expression, the factor α s /π includes the 2-loop running coupling discussed in section 3.1. The scales µ R = ξ R p ⊥ R, µ F = ξ F p ⊥ R and the factors ξ K and ξ Z probe the scale uncertainties and are discussed below. The factor p(x, j|i, t coll ) -computed numerically by solving Eq. (3.8) with an approach similar to Ref. [62] -encodes the probability for the leading parton to have a momentum fraction x and a flavour j, starting from a jet of flavour i (with initial fraction f i ) computed in the collinear limit as in section 3.2. Similarly, the factor zP j /(2C j ) accounts for the collinear structure associated with an observed Lund-plane emission at finite z (cf. e.g. Eq. (3.11)). Finally, the factor ρ soft resums the soft logarithms at large angles as well as C/A clustering logarithms, as described in section 3.3.3. In practice, ρ soft depends on the full colour structure of an event. We have computed it by interfacing Born-level events obtained with the NLOJet++ program to the numerical code from section 3.3.3. Each Born-level event is separated (at large-N c ) into different (weighted) dipole configurations. The result is binned as a function of t soft , ∆, and the jet p ⊥ and flavour. The different dipole configurations contributing to a given jet flavour are summed, as only the sum is needed to combine ρ soft with the collinear effects in writing (3.39). The structure of Eq. (3.39) is illustrated in Fig. 5, which shows that to obtain the density at a point (∆, k t ) in the Lund plane (the black dot), one first resums collinear effects down to the angle ∆ (the solid blue line) then resums the soft emissions at commensurate angles between xp ⊥ ∆ and k t (the solid red line). In particular, one sees that at large angles, where the details of the dipole configuration matter, collinear effects can be neglected in (3.39) and the sum over dipole configurations can be performed trivially. At small angles, the clustering logarithms resummed in ρ soft depend only on the jet flavour.
Our results for the resummation of the soft gluons are strictly-speaking obtained in the large-N c limit. It is however possible to restore the full-N c behaviour up to and including O(α 2 s ) in the collinear limit, i.e. in the limits that have been discussed in sections 3.3.1 ln R/∆ ln k t Figure 5: Schematic representation of Eq. (3.39) where the density at a point (∆, k t ) in the Lund plane is obtained by first resumming collinear effects at angles larger than ∆ and then soft gluons (at commensurate angles) down to a scale k t . and 3.3.2. First, we multiply the soft density for quark jets, ρ soft,q by a factor 2C F /C A to guarantee the proper result from Section 3.3.1 at order α s . Then, we multiply t soft by 2(C A − C F )/C A for quark jets, to guarantee the proper expansion, Eq. (3.25), at order α 2 s in the collinear limit. At large angles the structure of subleading-N C corrections is more complicated, and we will rely on matching with fixed-order calculations to address these terms up to order α 2 s . To obtain our final predictions integrated over p ⊥ , we have again used NLOJet++ to obtain the jet cross-section, the quark and gluon fractions f q,g , the average jet p ⊥ and the average α s (µ R ) (as well as ρ soft ) in a series of bins in p ⊥ . The contribution of each bin is evaluated using (3.39) at the average p ⊥ in the bin and summed with weight proportional to the bin cross-section.
Compared to section 3.3.3, the definition of t soft from Eq. (3.38) has to be adjusted to ensure t soft → 0 when k t → 1 2 xp ⊥ ∆ (i.e. z → 1 2 ). This is simply done by writing where we have introduced a parameter ξ Z that allows us, by a standard variation of ξ Z between 1/2 and 2, to probe the uncertainties associated with the resummation of soft gluons. Similarly, we estimate the renormalisation (µ R = ξ R p ⊥ R) and factorisation (µ F = ξ F p ⊥ R) scale uncertainties using the 7-point rule [73] around µ R = µ F = p ⊥ R (ξ R = ξ F = 1). The factorisation scale only influences the Born-level spectrum and the quark/gluon fractions f q,g . The choice of µ R should also be reflected in the factor α s /π in (3.39) as well as in the definition of t coll and t soft , via the reference scale µ R = ξ R p ⊥ R for α s in (3.3). Additionally, the uncertainty of the choice of scale for the argument of α s in (3.2) is taken into account by setting the scale to ξ K k t and varying ξ K between 1/2 and 2. This is the dominant source of uncertainty in our calculation. The uncertainty on the collinear resummation could be estimated by varying ∆ 0 in (3.39). However, since the effect of the collinear resummation is small (see e.g. Figs. 6 and 7), we have neglected this and set -19 - ∆ 0 = R. 12 To be conservative, the final perturbative uncertainty is obtained by summing in quadrature the three individual sources of uncertainties: the 7-point variation of µ R (or ξ R ) and µ F (or ξ F ), the variation of ξ K and the variation of ξ Z . 13 We present some representative results obtained with Eq. (3.39) in Figs. 6 (for slices of the Lund plane in a narrow bin of ∆) and 7 (for slices in a narrow bin of k t ). In each plot, we 12 Varying ∆0 would come with the additional complication that, for ∆0 > R, collinear radiation at angles larger than the jet radius would cause the Born-level p ⊥ and the jet p ⊥ to differ. Since, in our case, ln(R) is not large, we can neglect this effect. 13 Recall that t coll , Eq. (3.6) and t soft , Eq. (3.40) are written terms of Iα, Eq. (3.7) and that they all have a structure α n s L n , where each factor of L can be one of ln ∆ or ln ptR/kt. To probe uncertainties, we should examine variations that generate terms α n s L n−1 . The variation of µ in Eq. (3.7) does not generate such terms, but only terms α n s L n−2 . One approach to generating terms α n s L n−1 is to change the argument of αs within the integral in Eq. (3.7), i.e. replacing αs(qt) with αs(ξqt), where ξ is the scale variation factor. This is equivalent to replacing Iα(kt, kt0; µ) → Iα(ξkt, ξkt0; µ), i.e. changing both integration boundaries. A second approach is to change just one boundary by a factor ξ, which can be thought of as a replacement L → L±ln ξ. The prescription that we have adopted for t soft corresponds to the second approach, specifically varying the lower boundary (which has a larger numerical impact than varying the upper boundary). Ultimately, the choice we make here is not especially critical, because the overall perturbative uncertainty is dominated by the ξK variations. Those soft-gluon effects are most significant at small k t and at large angle. It is worth noting though that their effect is also visible at large k t in Fig. 6a. This is due to the power corrections in ∆ 2 starting at order α s , as discussed in section 3.3.1.
Collinear effects are small except close to the k t = 1 2 p ⊥ ∆ endpoint where the use of the full splitting function and the probability distribution for the momentum fraction of the leading parton have a clearly visible effect (see Fig. 7 in particular). Flavour-changing collinear effects are small but are still visible in Fig. 7a, reflected in the difference between the green and blue lines for ∆ 0.02. In particular, as one goes to smaller values of ∆, there is an increase in the fraction of jets whose leading parton is a gluon. This flavourchanging effect is modest in size, in part because the initial Born-level spectrum has a quark fraction of about 77.5%, relatively close to the asymptotic fraction of 62% that is visible in Fig. 2a (cf. Eq. (3.13)).
The perturbative scale uncertainties are about 10% at large k t , slowly growing to ∼ 15% at k t ∼ 20 GeV and to ∼ 30% at k t = 2 GeV (averaging the upper and lower uncertainties). They are dominated by the scale variation, ξ K , in the argument of α s with -21 -an additional small contribution from the variation of ξ Z at small k t . 14 While all the above expressions are given for the primary Lund-plane density ρ(∆, k t ), they can almost straightforwardly be adapted toρ(∆, z), as measured e.g. by the ATLAS collaboration. Specifically, Eq. (3.39) becomes We just note that, while keeping k t large enough in Eq. (3.39) guarantees that we stay in the perturbative region, the integration over x in (3.41) potentially extends to arbitrarily small xzp ⊥ ∆ momentum scales. This is regulated by our freezing of the running coupling at 1 GeV. In practice, this only affects the small values of z in a region where the nonperturbative corrections dominate anyway.
In anticipation of the matching of our resummed predictions to exact fixed-order results for ρ(∆, k t ), we note that our all-order equations (3.39) and (3.41) can be expanded to fixed-order. For Eq. (3.39), at NLO we have with α s ≡ α s (p ⊥ R). We have explicitly written α 2 s ρ

Matching with fixed-order
In order to get a full coverage of the primary Lund-plane density, including regions which are not dominated by large logarithms, it is useful to supplement our resummation with as many orders of the α s series expansion of ρ(∆, k t ) as are known exactly.
In this paper we focus on dijet events, for which we can obtain the primary Lund-plane density using the NLOJet++ program, The first (LO) and second (NLO) contributions are accessible using respectively LO and NLO 3-jet calculations [51].
Compared to the all-order calculation discussed in section 3, the LO contribution includes the first-order soft gluon radiation at large angles. The NLO contribution includes the first non-trivial running-coupling, flavour-changing and clustering corrections. 15 We have checked numerically that there was an agreement between NLOJet++ and our analytic calculations for the soft-and-collinear behaviour at O(α s ) and for the logarithmic dependence at O(α 2 s ), although small deviations expected from our large-N c approximationused to calculate dipole decompositions and soft logarithms beyond the collinear limitare observed at large angles. We show some explicit examples in Appendix C.
Knowing both the all-order resummation and the exact fixed-order results, we obtain a matched prediction using where ρ resum,NLO is the expansion to O(α 2 s ) of the resummed result (3.39). This expression is such that it reproduces the resummed calculation in the region where large logarithms are present, and the exact NLO result when expanded to second order in α s .
Explicit examples of matched predictions, including different levels of approximations for the resummation, are presented in Fig. 8 for the k t dependence at fixed ∆. First, we see that the exact NLO results (red) are close to what is obtained using the expansion of our resummed calculation (green). Next, the resummed result (blue) shows a strong enhancement at small k t , primarily due to the running coupling, and to soft-gluon clustering effects. Finally, the matched result (black) smoothly interpolates between the fixed-order result at large angle and large k t and the resummed result at smaller angle or k t .
The bands in the upper panel of Fig. 8 as well as the curves in the lower panel show our theoretical uncertainties. One of the striking features is that the matching with NLO reduces the uncertainties compared to the resummed result. This is valid across the whole kinematic range and especially visible at larger k t . Final uncertainties after matching are ∼ 6% at k t = 200 GeV, increasing to ∼ 12% at 20 GeV and ∼ 30% at 2 GeV.

Non-perturbative effects
Before making our final predictions it is interesting to estimate non-perturbative corrections to the calculation we have provided so far. We do so using a Monte-Carlo approach. We have studied the primary Lund-plane density using 5 different Monte-Carlo generators/tunes: Pythia8 (v8.230) [74] with the Monash 2013 tune [75], tune 4C [76] and the 15 In these fixed-order calculations the central renormalisation and factorisation scales have been set to p ⊥ R with p ⊥ the jet transverse momentum.  [82]. For each generator/tune we first study the primary Lund-plane density at parton level. We can then switch to hadron level to study the effect of hadronisation, include multi-parton interactions (MPI) to study the effects of the Underlying Event, and examine the impact of using only charged tracks as done in the ATLAS measurement [83].
For the central value of the non-perturbative corrections, we take the average of the Monte-Carlo generators, excluding Herwig7. The reason behind this exclusion is that our perturbative results are in the same ballpark as parton-level results from Pythia8 and Sherpa but differ significantly from parton-level Herwig7 results (see Appendix D). We obtain the (upper and lower) uncertainties on the non-perturbative corrections from the envelope of the Lund-plane density ratios for the 5 Monte-Carlo generators/tunes. To remain conservative, we keep the Herwig7 results in our non-perturbative uncertainty estimates.
Our results are presented in Fig. 9, for our high-p ⊥ setup separately for hadronisation and Underlying Event corrections. It is clearly visible that hadronisation corrections become sizeable at low k t , with a negative effect above ∼ 3 GeV and a positive effect below. Their effect is almost invisible for k t 10−20 GeV. Underlying-Event corrections are instead important (and positive) at low-to-moderate k t and large angles. The   between the initial anti-k t clustering and the C/A re-clustering as already discussed in [33] and related boundary logarithms discussed in section 3.3.2 and Appendix B.
The diagonal dashed line in Fig. 9 corresponds to k t = 1 2 p ⊥ ∆ for p ⊥ = 2 TeV. This is the kinematic limit for the lowest-energy selected jets. The Lund plane density quickly decreases above that line. The large fluctuations and uncertainties observed in Fig. 9 around the dotted line are a trace of the statistical fluctuations in our Monte Carlo samples.

Final predictions
Our final predictions include both the matched perturbative predictions, discussed in section 4, multiplied by the non-perturbative corrections obtained in section 5.
We show in Fig. 10a the resulting two-dimensional average primary Lund-plane density ρ(∆, k t ), and in Figs. 10b and 10c the associated relative uncertainty at perturbative level and at the non-perturbative level respectively. Fig. 11 shows slices at fixed angle ∆, which help to better visualise certain features. The density plot, Fig. 10a, shows all the expected features: the gradual increase towards small k t due to the running of α s ; the extra enhancement due to soft-gluon emissions, both at large angles and at small k t /∆ (or equivalently z); the reduction close to the kinematic limit associated with the "energy loss" of the leading branch; and the increase at low k t and in the bottom-left corner of the Lund plane due to non-perturbative effects.
The uncertainties are dominated by the perturbative component for k t 3−5 GeV, except at large angles where non-perturbative effects can have a sizeable impact up to k t ∼ 10−20 GeV. The total uncertainty is found to be about 20% at k t ∼ 5 GeV (away from the large-angle region), and decreases to 5−7% for k t in the 200−500 GeV range. Relative to the LO+resum results, visible in Fig. 11, the inclusion of NLO corrections reduces the uncertainties mainly at high k t . Even if the non-perturbative corrections have -26 - High-p setup: 0.074 < < 0.091 (b) small angles: 0.074 < ∆ < 0.091 Figure 11: Slices of the primary Lund-plane density ρ(∆, k t ) at constant ∆. The upper panels correspond to the density ρ(∆, k t ) itself while the lower panels show the relative scale uncertainties δρ ρ (∆, k t ). We show results for the matched result at LO (red) and NLO (blue), as well as NLO results including non-perturbative corrections (black). a negligible impact on the uncertainty above ∼ 3−5 GeV (10−20 GeV) at small (large) ∆, they result in a (small-but-visible) shift of the central value up to larger values of k t .
Finally, we discuss our analytic calculations supplemented with non-perturbative corrections forρ(∆, z), corresponding to the ATLAS setup. Besides the differences discussed in section 2, we follow the same strategy as for the high-p ⊥ setup: the resummed prediction is obtained using Eq. (3.41), matched to NLOJet++ fixed-order results using Eq. (4.2) and supplemented with non-perturbative corrections -this time correcting so as to correspond to a measurement performed using charged-tracks above 500 MeV -following the procedure outlined in section 5. Details of the non-perturbative corrections are given in Appendix E.
We compare our results to the ATLAS data from Ref. [83] for slices in ∆ in Fig. 12 and slices in z in Fig. 13. The vertical dashed lines correspond to the k t scales estimated using z = k t /(p ⊥ ∆), i.e. assuming a jet at the lower p ⊥ cut of 675 GeV and a leading parton/subjet carrying a fraction x = 1 of the initial jet transverse momentum. The shaded grey bands indicate regions where the uncertainty on the non-perturbative corrections is larger than 10%. Shaded red bands correspond to the regions sensitive to the boundary logarithms discussed in section 3.3.2. We recall that we have not resummed these terms, so our calculation should be considered incomplete in the red shaded regions. A rough estimate of their potential size is given in Appendix B.
-27 -  Figure 12: Comparison between our calculations and the ATLAS measurement from Ref. [83], for different bins of ∆. The dashed vertical lines, corresponding to z = kt p ⊥ ∆ for p ⊥ = 675 GeV and several k t values, are meant to indicate the transverse scales one is typically sensitive to. The shaded grey bands indicate bins where the relative uncertainty on the non-perturbative corrections is larger than 10%. The shaded red regions indicate that our calculation is incomplete because of the missing resummation of the boundary logarithms.
For all unshaded bins in Figs. 12 and 13, we see agreement between our predictions and the data to within the experimental and theoretical uncertainties. Generally speaking, the theoretical uncertainties are larger than the experimental ones, though they are comparable at values of z and ∆ that correspond to large k t values. Recall that the theoretical uncertainties are to a large extent dominated by the choice of scale of α s in the resummation and a higher-order resummation would therefore be beneficial to reduce the uncertainties.
If we consider the grey shaded regions, i.e. those where non-perturbative uncertainties are larger than 10%, the agreement between data and theory remains good to within the total uncertainties in most of the bins, almost all the way down to 2 GeV. In practice this -28 - agreement is facilitated by the non-perturbative blow-up of the uncertainties at low k t and our predictions' central values are systematically above the data points. Recall, however, that our estimates of non-perturbative corrections rely on the assumption that the partonlevel event-generator results are structurally similar to a full perturbative calculation. This assumption is questionable at low k t : for example, a parton shower may contain a low-k t cut, with the phase-space below that k t value being filled up by hadronisation (there is a hint of this occurring in Appendix D, Fig. 16); in contrast our perturbative calculation has no such cut, and so the hadronisation contribution to that region, supplemented with our perturbative contribution, could effectively lead to double counting and so an overestimate relative to the data. In this respect it might be interesting to develop a more analytic understanding of expected hadronisation effects on the Lund plane density.
One region where there is clear disagreement between our predictions and the data is in the (red-shaded) largest angle bin 0.287 < ∆ < 0.400. This disagreement is only mildly alleviated by our estimate of the potential size of boundary logarithms, cf. Fig. 14 in Appendix B. Several avenues could be of interest for further exploring this region, for example a full resummation of the boundary logarithms, or a measurement with jets whose original clustering was with the C/A algorithm (rather than anti-k t ), so as to remove -29 -these boundary logarithms altogether. Note also that this region is potentially sensitive to underlying-event effects, and if they are incompletely modelled in event generators, this could also contribute to the disagreement. 16

Conclusions and outlook
In this paper we have carried out the first calculation of all-order logarithmic contributions to the average primary Lund plane density. We have resummed three classes of single logarithmic terms: (i) running-coupling effects, which are relatively straightforward and the numerically dominant contribution over most of the Lund plane; (ii) soft effects, which involve large-angle contributions and clustering logarithms, both evaluated in the large-N C limit; and (iii) collinear effects at large momentum fractions, which include contributions that can change both the momentum and the flavour of the leading parton. We have also discovered a new class of logarithmic effects in jets that arise when reclustering an anti-k t jet's constituents with the Cambridge/Aachen algorithm. The corresponding terms are relevant close to the large-angle boundary of the Lund plane. We defer their full single logarithmic resummation to future work.
For the purposes of making phenomenological predictions, we have matched our allorder, resummed, calculation to an exact (3-jet) calculation at next-to-leading order with the NLOJet++ program. We then supplemented the perturbative predictions with nonperturbative effects extracted from Monte Carlo simulations with Herwig, Pythia and Sherpa.
The theoretical uncertainty on our perturbative predictions ranges from 5−7% at large k t to ∼ 20% at k t ≈ 5 GeV. Hadronisation and underlying-event corrections are relevant below 20−30 GeV, but in most of the Lund plane dominate the overall uncertainty only below k t ∼ 3−5 GeV (15 GeV at large angles, where the underlying event is a significant contributor), cf. Figs. 10 and 11.
We have made our predictions for two variants of the Lund plane definition, one using angle and absolute transverse momentum (the default for most of this paper), and the other using angle and relative momentum fraction in a given branching. The latter corresponds to the choice made by the ATLAS collaboration in their recent pioneering unfolded measurement of the primary Lund plane density with charged tracks [48]. 17 We have compared our results to the ATLAS data, including an additional non-perturbative correction to account for the use of charged tracks, and found good agreement in all regions where we have confidence in our predictions, i.e. the non-shaded regions of Figs. 12 and 13. This includes a broad swathe of the Lund plane, down to scales corresponding to transverse momenta of about 5 GeV.
Our work opens a series of questions that should be kept in mind for future work. First, it would be interesting to extend our calculation beyond single-logarithmic accuracy. We expect that this would give a considerable reduction in the uncertainty, notably from control over the effective scale to be used in the coupling. Such a calculation, however, remains challenging. Other effects like the subleading-logarithmic and subleading-N c corrections to the clustering logarithms, or NNLO fixed-order corrections (requiring a NNLO pp → 3-jet calculation) would also be expected to bring significant improvements in certain specific regions of the Lund plane.
It would also be of interest to understand the resummation of the boundary logarithms that originate from the interplay between the initial anti-k t clustering and the C/A reclustering. Practically, however, these logarithms could be avoided by using the C/A algorithm for both the initial clustering and the reclustering. One observation that would also deserve better understanding is the apparent absence of any resummation effect from clustering logarithms in the soft-collinear part of the Lund plane for gluon-induced jets Keeping the above theoretical limitations in mind (and possible future improvements), one might wish to investigate whether a measurement of the Lund plane density, which intrinsically covers a wide range of transverse-momentum scales, could be helpful to make an extraction of the strong coupling constant, α s , extending existing work on strong coupling determinations from soft-drop measurements [84,85]. In a similar spirit, one could perhaps extend the approach of Ref. [86] to develop an analytic approach to non-perturbative corrections at small k t and potentially even use Lund-plane measurements to determine an effective coupling constant down to small transverse momenta.
Finally, it would be interesting to compare both analytical predictions and measurements of the primary Lund-plane density to recent efforts to develop parton showers with perturbative control beyond leading double logarithmic accuracy (e.g. [47,87]) and leading colour (e.g. [88,89]).

A Analytic results for collinear resummation
In this Appendix we give the explicit analytic solutions for the average quark/gluon fractions f (j|i, t coll ) and momentum fractionx(j|i, t coll ) defined in Eq. (3.12). We find and Note that the coefficients s q and s g are in agreement with the flavour-changing effects calculated in [62].

B Boundary logarithms for ∆ ∼ R
In this Appendix, we show how new logarithms of R − ∆ arise from secondary emissions at order α 2 s . We show that these logarithms are a consequence of the interplay between the initial anti-k t clustering used to obtain the initial jets and the C/A clustering used to construct the primary Lund plane.
Say we start from a dipole ( i , j ) and have 2 emissions, k 1 and k 2 , strongly ordered in transverse momentum as discussed in Section 3.3.2. Emission 1 (real or virtual) is integrated over and the softer emission 2 (real) is measured as a contributing to ρ (2) soft (k t , ∆). We denote by (k|ij) the geometrical pattern associated with the radiation of gluon k from the dipole ( i , j ) (i.e. the transverse momenta with respect to the beam are factored -32 -out). We also denote by R i the distance of i to the jet axis (in rapidity-azimuth) and R ij the distance between i and j. For a parton with Casimir C R , we have If one combines the C R contributions, performs the k ⊥i integrations and switches to polar coordinates for the y 2 , φ 2 integration and uses the δ(∆ − R 2 ) constraint to simplify, we get We can evaluate this numerically, separating the C R C A term in an "inside" contribution where k 1 is inside the jet (integrated in polar coordinates around the jet axis) and an "outside" contribution where k 1 is outside the jet (integrated directly in y 1 and φ 1 ). We have done this explicitly as a check of the Monte-Carlo implementation introduced in section 3.3.3 and found perfect agreement (in the large-N c limit). Note that the combination of dipoles in the first square bracket of Eq. (B.2) vanishes when y 1 → ±∞, showing explicitly that there are no divergences collinear with the beam. The main purpose of this Appendix is to show that the "out" C R C A contribution has a collinear divergence when ∆ → R. To see this, we set ∆ = R − with → 0 (or take ln(R/∆) → 0). The collinear divergence comes from the situation where emission k 1 is close to emission k 2 (with k 1 outside the jet and k 2 inside), where the combination of dipoles can be simplified to 4/θ 2 12 . After a few straightforward manipulations, we reach soft (R, k t ) can be taken from Eq. (3.22). This exhibits a logarithmic behaviour when ∆ → R (which is integrable if one considers a bin in ∆ between some lower bound and R). 18 We have checked that this behaviour is reproduced by the Monte-Carlo described in section 3.3.3.
The physical origin of the collinear enhancement in (B.4) is the interplay between the anti-k t clustering used to obtain the initial jet and the C/A clustering used to construct 18 A similar enhancement was observed for narrow slices in Ref. [65].
-33 -the primary Lund plane. For a jet initially clustered with the C/A algorithm, emissions k 1 and k 2 would be clustered together and emission k 2 would then not be seen as a primary emission.
Equation (B.4) exhibits a double logarithmic behaviour. One should also expect singlelogarithmic corrections, proportional to ln(R/(R − ∆)) without the soft enhancement. In principle, these single-logarithmic terms should be resummed to all orders. We have not, so far, found a simple prescription to achieve this resummation, which involves an interplay between the complex structure of soft emissions at commensurate angles and (potentially hard) collinear splittings at the boundary of the jet.
We however give in this Appendix a simple (incomplete) prescription from which one can gauge the potential impact of this resummation. Coming back to our calculation at O(α 2 s ), we see that the boundary logarithm ln(R/(R − ∆)) comes from the fact that emissions k 1 and k 2 are collinear to each other and the logarithm ln(p ⊥ R/∆) comes from the energy ordering between the two emissions. One would obviously get a single-logarithmic contribution if the two emissions were still collinear but no longer strongly ordered in energy. This contribution, where the first emission is soft and just outside the jet, and the second emission is collinear to the first one and inside the jet, can be straightforwardly computed. A calculation similar to the previous one shows that the energy logarithm is replaced by an integration over the Altarelli-Parisi splitting function P (z) with z the momentum fraction of the collinear branching. One then gets with B g = − 11 12 + n f 6C A the standard gluon hard-collinear branching contribution obtained from integrating the finite part of the gluon (to anything) splitting function. This is but the first of a tower of terms enhanced by logarithms of R/(R − ∆). We will examine its magnitude shortly. The full structure of the series involves other potentially complicated effects: (a) an interplay between non-trivial clustering logarithms and these new purely collinear effects; and (b) the way in which the anti-k t jet clustering affects the jet axis and subsequent identification of the set of particles (or tracks) that gets reclustered with the C/A algorithm, specifically in presence of hard splittings at angles comparable to the jet radius. In addition to these subtleties, one might want to consider a number of combinations of jet clustering: e.g. reclustering a full anti-k t jet with R C/A = ∞, reclustering it with R C/A = R anti-kt , reclustering only the particles within a distance R anti-kt of the anti-k t jet axis, etc. Given that these effects concern only a single bin in ∆, and that their treatment brings many complications, we postpone their study to future work.
We do nevertheless wish to investigate the size of the one contribution we have outlined in Eq. (B.5). It can be included in the all-order resummation via the following redefinition of t soft   Figure 15: Differences between the exact NLO results and the NLO expansion of our resummation for different slices of the Lund plane. Control of single logarithmic terms in our resummation implies that these differences should tend to a constant both for small k t and for small ∆, as observed, except as concerns subleading-N C terms in the small-k t , large ∆ region. The exact NLO results are normalised here to the Born-level jet cross-section.
effect is relatively small (in particular, relative to our uncertainties and to the discrepancy with the data), we see that our results move in the right direction. Pending a full treatment of these boundary logarithms -left for future work -the bin closest to the jet edge should be treated with caution. We signal this limitation by shading the corresponding region in red in our overall comparisons with the ATLAS data, Figs. 12 and 13.

C Validation of the resummation at NLO
As with any resummed calculation, it is important to check that its expansion to fixed order reproduces the behaviour seen in the exact fixed-order calculation, to within the expected accuracy of the resummation. In our case, this means that, at NLO, one should reproduce all contributions of the form α 2 s ln, where the argument of the logarithm is any variable in the Lund plane. In practice, one therefore expects the difference ρ NLO − ρ resum,NLO to tend to a constant when k t becomes small at a fixed ∆, or when ∆ 1 at a fixed k t . Fig. 15 shows that this is indeed the case for both limits. We note that, while in the main text of the paper, the NLO Lund-plane density has been normalised to the NLO inclusive jet cross-section, for the purpose of Fig. 15 both ρ NLO and ρ resum,NLO have been normalised using the Born-level jet cross-section.

D Comparison between our calculation and Monte Carlo simulations
We show in Fig. 16 a comparison between our analytic calculations and Monte-Carlo simulations for a slice of the Lund plane at constant angle. In Fig. 16a we compare our perturbative predictions to parton-level simulations and the (blue) uncertainty band corresponds to our perturbative scale uncertainty. In Fig. 16b the comparison is made for the full prediction, including non-perturbative corrections. 19 At hadron+MPI level, we see a globally-decent agreement between our results and those from each Monte Carlo event generator. At parton-level however, the Herwig7 results are systematically much smaller that our analytic results for k t below ∼ 10 GeV. This is

E Non-perturbative corrections for the ATLAS setup
The set of non-perturbative corrections included in our calculation ofρ(∆, z) for the "ATLAS setup" (cf. section 2) differs from those included in our default "high-p ⊥ setup." The main differences are (i) the use of charged tracks instead of all particles, (ii) a slightly different clustering procedure using a radius of 0.4 for the C/A reclustering, and (iii) the selection of tracks above 500 MeV and within a distance to the jet axis calculated using pseudo-rapidity instead of rapidity.
In Fig. 17, we split the full non-perturbative correction into two separate factors: the corrections due to hadronisation and multi-parton interactions computed on all particles, -37 - Fig. 17a, and the extra corrections associated with the use of charged tracks, the 500 MeV transverse-momentum cut and the track selection based on pseudo-rapidity, Fig. 17b. 20 The final set of corrections is shown in Fig. 17c. We clearly see that the use of just charged tracks with momenta greater than 500 MeV induces a strong reduction ofρ(∆, z) beyond the perturbative domain (k t < 2 GeV) and a positive correction for k t > 2 GeV. This effect partially cancels the original effect of hadronisation in the region where hadronisation depleted the Lund plane density, i.e. k t 5 GeV. The corresponding uncertainties on the non-perturbative corrections are shown in Fig. 17(d)-(f). While the additional corrections associated with the selection of charged tracks add a little to the uncertainty at large z, the final pattern of uncertainties is largely unmodified compared to that obtained solely from the hadronisation and MPI uncertainties.
[5] ATLAS Collaboration, G. Aad et al., Measurement of the jet mass in high transverse momentum Z(→ bb)γ production at √ s = 13 TeV using the ATLAS detector, arXiv:1907.07093.  20 The effect of reclustering the jet constituents with a finite jet radius R = 0.4 is included together with the hadronisation and MPI corrections. We have checked that this effect itself is very small, below 0.1%.