Resummation of non-global logarithms in cross sections with massive particles

Abstract
 A factorization formalism for jet processes involving massive colored particles such as the top quark is developed, extending earlier results for the massless case. The factorization of soft emissions from the underlying hard process is implemented in an effective field theory framework, which forms the basis for the resummation of large logarithms. The renormalization group evolution giving rise to non-global logarithms is implemented into a parton shower code in the large-Nc limit. After a comparison of the massive and massless radiation patterns, the cross section for $$ t\overline{t} $$
 t
 
 t
 ¯
 
 production with a veto on additional central jet activity is computed, taking into account radiation both from the production and the decay of the top quarks. The resummation of the leading logarithms leads to an improved description of ATLAS measurements at $$ \sqrt{s} $$
 
 s
 
 = 7 TeV.


Introduction
The study of jet cross sections plays a crucial role in high-energy physics. While theoretical calculations are carried out in terms of interactions at the field level, detectors are only able to measure properties of outgoing particles after they have fully hadronized, i.e. transformed from colored quarks and gluons to color-neutral final states such as mesons. Consequently, it is impossible to measure the underlying hard scattering process directly, but one needs to reconstruct it by measuring jets and analyzing their properties.
While the total energy of the particles inside the jets is typically of the same order as the partonic center-of-mass energy of the collision, the total energy of the particles not ending up in a jet is considerably lower. Due to this scale separation effective field theory methods, in particular Soft-Collinear Effective Theory (SCET) [1][2][3] (see [4][5][6] for reviews), are useful in the study of jet cross sections. In the effective theory the cross sections factor into hard, collinear and soft functions, each of which can be safely evaluated in fixed-order perturbation theory at their characteristic energy scale. To connect these factors, it is then necessary to evolve one of the factors from its characteristic scale to the scale of the other factor by using the Renormalization Group (RG) equation. This procedure was first applied to jet cross sections in [7,8].
Because of their multi-scale nature, jet cross sections are sensitive to potentially large logarithmic corrections. When evaluating the phase-space integrals of matrix elements and applying cuts to the allowed energies, logarithms of the ratios of the energy scales involved in the process appear in the calculations. For example, when the energy of particles inside the jets (denoted by Q) is unconstrained and of the same order as the partonic centerof-mass energy, i.e. Q 2 ∼ŝ, but the energy outside the jets (denoted by Q 0 ) is required to be small, the phase-space integrals produce terms proportional to ln (Q/Q 0 ). These logarithms become large if Q 0 Q. The factorization formula studied in [7,8], derived within the effective field theory approach, can be used to resum these corrections, in principle to all logarithmic orders. Based on this theoretical framework, a dedicated parton shower code was developed and applied to resum the large logarithms appearing in jet processes and isolation-cone cross sections up to leading logarithmic (LL) order in [9]. Subsequently, higher-order matching corrections in both the hard and the soft function were added. This led to the resummation of the interjet energy flow up to LL accuracy and to the resummation of the jet mass up to next-to-leading-logarithmic (NLL ) accuracy in [10]. As usual, the prime in LL and NLL indicates that the matching corrections are included one order higher than what it would be required in RG improved perturbation theory. In the present case, this means that NLO hard and soft functions were used. By supplementing these calculations with the two-loop corrections to the anomalous dimension matrix one would achieve full NLL and next-to-next-to-leading logarithmic (NNLL) accuracy, respectively.
The work done so far was carried out in the high-energy limit where all partons can be considered massless. The purpose of this paper is to extend the approach of [7][8][9][10] to processes involving heavy colored particles and to develop and validate a parton shower code for the resummation of jet cross sections in top-quark production processes. Soft radiation is obtained from matrix elements of Wilson line operators along the directions of the emitting particles, independently of the mass of the emitting parton. Because of this fact, the factorization theorem has the same general form as in the purely massless case. However, the soft radiation pattern and its generation by the parton shower code differ significantly in the two cases. At one-loop order, the angular dependence of the radiation of a soft parton with momentum k µ = E n µ k between legs carrying momenta p i and p j is given by the usual product of eikonal factors This factor is the same in both cases, but massless particles are traveling along lightlike directions, while massive particles travel along time-like directions. This difference in kinematics must be accounted for in the shower code. Furthermore, in contrast to what happens in the high-energy limit, the radiation factor in (1.1) does not vanish when i = j, if p i is a time-like momentum. Therefore, in addition to the usual dipole emission pattern, φ y y min −y min y max −y max b-jetb -jet 2π Figure 1. Sketch of the veto region as defined by ATLAS in [15]. The gap, in which additional radiation is vetoed, is represented by the shaded red area with rapidity y min < |y| < y max . Radiation inside the b-tagged jets is not vetoed. For y min = 0, this setup reduces to the usual central jet veto.
it is necessary to include monopole contributions in the massive case. The latter describes radiation that is emitted and absorbed by the same Wilson line rather than exchanged between two color-connected Wilson lines. This difference in the massive and massless radiation pattern is of course well known, in particular the different collinear behavior, which is often referred to as the dead cone effect [11][12][13][14].
As an application of the new parton shower code described in this work, we consider tt production with a veto on additional central jet energy. This process was measured by ATLAS at the Large Hadron Collider (LHC) with the goal of testing the description of soft radiation in parton showers [15]. The top pair production process involves two initial-state partons producing a tt-pair in the final state. The top quarks then decay into bottom quarks and W bosons. The measurement is performed using events in which the W 's decay leptonically and in which two b-jets are detected. The veto on central jets is imposed by requiring that, with the exception of the two bottom-tagged jets, no additional jets above a given transverse momentum Q 0 are allowed to be present in the rapidity range y min < |y| < y max (see Figure 1). With the veto, only particles of low energy are allowed inside this rapidity range, while the energy is unconstrained anywhere else. This is a typical situation in which large non-global logarithms appear. In this work these logarithms are resummed at LL accuracy and the results of the resummation are matched to NLO predictions in fixed-order perturbation theory.
In addition to radiation effects associated with the production process, one should also include radiation emerging from the decay products of the top quarks. We work in the narrow-width approximation for the top quarks, in which they are treated as stable particles and the process factorizes into a production cross section multiplied by the decay of the top quarks. It is well known that radiation from the b-quarks that would contribute to non-factorizable corrections in fixed-order perturbation theory is suppressed by factors of O(Γ t /m t ) [16][17][18][19][20][21]. To account for the factorizable contributions, we run a separate shower for the top decay to also account for the b-quark radiation. Numerically, the effect of this  Figure 2. Diagram for the process gg → tt → bb l + l − νν. In the large-N c limit, the radiation can be split into a set of color dipoles. The color dipoles associated to the production of the tt pair are shown in blue, the ones associated to the decay in green. The full LL cross section will include the emissions from all five dipoles. radiation is smaller than the one from the production of the top pair since the radiation inside the b-jet is not constrained. However, the radiation from the decay is large enough that it must be taken into account. Figure 2 shows one of the several tree-level diagrams contributing to the tt-pair production process measured by ATLAS in [15]. We also depict the color dipoles, which are the source of the emissions in the large-N c limit.
The remainder of this paper is organized as follows. In Section 2 the factorization theorem [7,8] is reviewed and the changes needed in presence of massive partons are discussed. Section 3 describes in detail how the relevant phase-space integrals can be evaluated in the parton shower code. In Section 4 we assess the impact of massive partons in the resummation of non-global logarithms. An explicit example of the resummation of non-global logarithms for a cross section involving top quarks at LL accuracy is presented in Section 5. As indicated above, the observable we consider is top-pair production with a veto on central jet energy. The predictions for this observable are then matched to the NLO result and compared to experimental measurements carried out by ATLAS. Section 6 contains our conclusions and an outlook. In Appendix A, we use a sample event to illustrate our parton shower code step by step. In Appendix B we explain how to use the shower to also compute the first two orders of the fixed-order expansion of the resummed result.

Factorization for cross sections involving massive quarks
Before discussing the factorization of the cross section, we should determine which scales are present and which scale hierarchies can arise in the observable under study. Throughout this paper, we consider scattering processes at a large center-of-mass energy Q and impose a veto on radiation in a certain phase-space region. We are interested in a regime where the energy scale Q 0 of the soft radiation in the veto region is much smaller than Q. The presence of the massive particles introduces additional scales in the process. On top of the masses themselves, which we denote generically with M , the most important new scale is the production threshold: when massive particles are produced, only part of the energy Q is available for additional radiation. We denote by Q 1 the energy above threshold that can be radiated. For tt production, the threshold is at Q T = 2m t and the maximal energy available to radiation is Kinematically, this value corresponds to a configuration where a collinear tt pair recoils against a gluon. This is a corner of phase space and the typical gluon energy will be much lower. However, the scale of the hardest possible emission plays an important role since it corresponds to the large scale in the emission process which should be compared to the veto scale Q 0 . Since we are interested in non-global logarithms associated with soft radiation, we only consider Q 1 Q 0 , but even under this assumption, one should consider two different hierarchies, namely a) Q ∼ Q 1 and b) Q Q 1 . The simpler of the two cases is Q ∼ Q 1 , which implies that the process energy is far larger than the threshold energy, and that the masses are smaller than the maximum emission energy, Q 1 M . It is then interesting to ask what role the masses themselves play and whether we encounter logarithms of the masses. If the heavy partons are not in the veto region, the vetoed cross section is collinear safe and mass effects are power suppressed in the limit M → 0; the massless limit is smooth. On the other hand, if the massive partons are inside the veto region, the limit M → 0 becomes complicated. Of course, in top-pair production, several additional complications arise when considering the limit m t → 0. In this paper, we only consider the case where M is larger or of the same order as Q 0 .
In the case in which Q Q 1 instead, the process occurs near threshold and the emissions are always soft compared to the particle masses. At the same time, we want to have Q 1 Q 0 , therefore the distance of Q from the threshold must still be large compared to the scale of soft radiation. Phenomenologically, this situation can only be relevant for top quarks and quite stringent vetoes. Since the radiation is always soft compared to the heavy particles, we should describe the entire process in Heavy Quark Effective Theory (HQET) (for reviews see [22,23]). One would first match QCD onto HQET at the scale Q and evolve to Q 1 before computing the soft emissions. It should be noted that this first matching will also have to be performed for the total cross section in the same kinematic regime. The effect will therefore largely cancel in ratios of cross sections such as the gap fraction. Furthermore, if one gets very close to the threshold Q 1 ∼ M α s , the heavy quarks become non-relativistic, but in view of Q 1 Q 0 this regime is not important phenomenologically. When we apply our formalism to top-pair production at the LHC at √ s = 7 TeV, we find that the average Q of the partonic collisions is Q ≈ 520 GeV ∼ 3m t and Q 1 ≈ 150 GeV. Therefore, in the phenomenological application considered in this work the scale hierarchy lies in between cases a) and b).
We discussed the two scenarios a) and b), but together with the scale M , also combinations of scenarios can arise. For example, for Q 1 M , it is possible to emit additional massive partons (at leading logarithmic accuracy only gluons are emitted). Then, for Q 1 M Q 0 , one could imagine a two step procedure, where one would start in scenario a), but after a number of emissions, only a small energy is left and one would switch over to scenario b). Massive theories have a much richer set of kinematic configurations than massless ones and can involve complicated interplays of different scales.
Here we first describe factorization for the simple case a), restricting ourselves to e + e − cross sections for the moment. After presenting the result, we discuss how it should be modified to account for the case b). The factorization formula for a jet production with a veto on radiation in part of the phase space takes the form . . , n m }. As discussed above, the hard functions H m also depend on the particle masses and derived quantities such as Q 1 . In order not to overburden the notation, we suppress this dependence. The symbol ⊗ indicates an integral over the directions of the m particles and . . . denotes the color trace which is taken after combining the functions. Up to the fact that some reference vectors are time-like, this formula is identical to the one studied in [8][9][10].
The hard functions H m are free of large logarithms if one chooses a value µ ∼ Q for the renormalization scale. The same is achieved for the soft functions S m for µ ∼ Q 0 . For Q Q 0 , at least one of these two functions will involve large logarithms, irrespective of the scale choice. These large logarithms can be resummed by solving the RG equation of the hard function and evolving it from its characteristic scale µ h ∼ Q down to a soft scale µ s ∼ Q 0 , leading to where the evolution factor is just the path-ordered exponential of the anomalous dimension RG evolution generates additional massless particles and⊗ denotes the integration over their directions before integrating over the hard directions.
Up to now we worked under the assumptions that Q 1 ∼ Q. Alternatively, if Q Q 1 , the hard functions involve large logarithms of the ratio Q 1 /Q which are not resummed by the above treatment. In order to factorize the two scales, one must first match onto HQET.
For e + e − → γ * → tt, it is necessary to match the electromagnetic current operator which induces the process onto the corresponding HQET operator where h v and h v are the two HQET fields describing the top and the anti-top quarks. One would then derive an expression analogous to (2.3) in HQET. The hard functions arising would be related to the ones in (2.3) by To resum the logarithms of Q 1 /Q, one will first solve the RG of C V to run from the scale µ ≈ Q down to µ ≈ Q 1 . When computing the gap fraction, one will also compute the total cross section in HQET using (2.5). The anomalous dimension driving the running of C V is the massive cusp anomalous dimension with the cusp angle defined by the two vectors v and v . In the ratio defining the gap fraction, the Wilson coefficient C V and its running will drop out. The situation is more complicated for hadron colliders, which involve sums of different partonic channels with different running so that the cancellation between numerator and denominator will not be complete. The general form of the anomalous dimension for a process with massive partons was given in [24] and the explicit forms relevant for top production can be found in [25], but we will not study the small effect of this running in this work. However, an important lesson from the above discussion is that one should set the scale µ h ∼ Q 1 in observables such as the gap fraction, since most of the running above this scale will drop out in the ratio of cross sections. In (2.4) we have presented the formal solution to the evolution equation. We will now discuss how the general solution simplifies at LL accuracy and how it can be implemented as a parton shower. In dijet processes at lepton colliders, one only needs to consider the case l = m 0 = 2 at LL, as the contribution of additional partons to the hard function would be suppressed by additional powers of α s for µ h ∼ Q. On the other side of the energy spectrum, the LL soft function is simply the unit matrix in the color space of the m final-state partons, since any soft correction would again be suppressed by a factor α s at the low scale µ s ∼ Q 0 . When computing tt production at a future electron-positron collider with a sufficiently high center-of-mass energy at LL accuracy, the general result (2.3) therefore simplifies to The situation is more complicated at hadron colliders such as the LHC, where the initial state contains two additional colored hard partons, which give rise to non-perturbative Parton-Distribution Functions (PDFs). In addition, Glauber gluons can induce interactions between soft and collinear partons. This complication is absent in the large-N c limit in which we perform our computations. In this limit, the only difference to the e + e − case is that there are two additional Wilson lines which describe the soft initial-state radiation.
For LL resummation one needs the anomalous dimension only at one-loop accuracy. Consequently, the exponent of the evolution matrix in (2.4) reduces to The "evolution time" t measures the separation of the scales µ s and µ h : one finds t = 0 for µ s = µ h and a growing t for increasing separation µ s < µ h . As the soft scale approaches the Landau pole, one finds t → ∞. If the scale µ h is kept fixed the function t ≡ t(µ s ) is bijective.
The discussion so far applies both to massive and to massless partons. The difference between the two cases becomes evident when one considers the one-loop anomalous dimension matrix where the matrix elements R m and V m (which are themselves matrices in color space) are associated with the emission of a real or virtual soft gluon The color matrices T i,L act on the hard function from the left, i.e. on the amplitude, while T i,R act on the conjugate amplitude. The function Θ in enforces that the hard emission is inside the allowed region. The factor Π ij is equal to +1 if i and j are both incoming or outgoing legs, and equal to 0 otherwise. When considering both massless and massive partons, the dipole radiator takes one of the following forms: mixed: . (2.14) In the special case of i = j (which can not occur in the mixed case (2.13), as it implies that the two legs are the same), the radiator (2.12) vanishes for massless legs, as n i · n i = 0, but is non-zero for massive quarks (2.14). The different kinematics and the presence of the monopoles distinguish the massive from the massless case. As mentioned above, we work in the large-N c limit in which the color structure becomes trivial and reduces to factors of N c . This is a huge simplification over the general case in which the m-parton terms act in the color-space of the m-partons. There is currently a large effort by several groups aiming to extend parton showers beyond the large-N c case, but we restrict ourselves to this limit. The fact that the color structure becomes trivial implies that the Glauber phases in V m in (2.10) vanish. Furthermore, all interference effects are suppressed and exchanges are only possible between neighbouring legs. However, the monopole contributions are present and need to be included, as is obvious from the diagrams shown in Figure 3.
The full corrections in the large-N c limit read The sum includes all dipoles i consisting of the legs i and i + 1 and we have absorbed the monopole contributions into the dipoles by defining In the rest of this work, the framework discussed here is applied to top-pair production. In this case the massive legs are always chosen to be the the first and the last in the list of Wilson-line directions, so that monopole radiation can only occur at i = 1 and i = m − 1, as the monopole radiator W k ii is manifestly zero for the massless gluonic legs in between. In Figure 3, we have depicted all possible real emissions for one dipole of two massive Wilson lines. The relative sign of the dipole and monopole contributions in (2.17) can be understood intuitively by looking at the figure: the partons in the dipole have opposite charge, in contrast the monopole. The factor of two of the dipole term compared to the monopole ones arises because one has to add the identical contribution of the two dipoles (ij) and (ji).
The details on how one gets from the RG equation to a parton shower are thoroughly explained in [9], but for completeness we briefly review the derivation here. The parton shower is based on the RG equation of the hard function which reads By changing variable from the scale µ to the evolution time t and by making use of the fact that the one-loop anomalous dimension matrix has the simple form (2.9), the evolution equation at LL accuracy takes the form This differential equation (2.19) can also be rewritten as an integral equation: Starting from (2.20), one can generate the hard functions in an iterative way as since H k (0) = 0 for k > 2. The cross-section at LL finally reads The iterative structure of (2.21) is well suited for implementation into a Monte Carlo code which generates successive emissions and thereby also performs the angular phase-space integrals of (2.22). For later convenience, we introduced the quantity R(t) given by the ratio of the resummed cross section with a veto to the inclusive cross section σ tot . At LL accuracy, one can replace σ tot by the Born-level result σ 0 . The inclusion of the massive Wilson lines into the Monte Carlo code is achieved in a straightforward way. The change compared to the massless case boils down to implementing the angular integrations in (2.15), where the modified dipole emitter W k ij replaces the massless one. A general algorithm for the evaluation of the angular integrals is discussed in the next section. The details of the Monte Carlo algorithm, which showers tree-level event files obtained by means of MadGraph5_aMC@NLO [26], are presented in Appendix A.

Evaluation of the massive angular phase space integrals
The goal of this section is to evaluate the integral for arbitrary Wilson lines u i and u j , which are either both massless (u i = n i and u j = n j ), both massive (u i = v i and u j = v j ) or one massive and one massless (u i = v i and u j = n j or vice versa). For the discussion below, it is convenient to normalize all reference vectors in such a way that the zero component of the four-vector is equal to one, i.e.
where E i is the energy component of the vector p µ i . With this convention, one finds that This differs from the definition adopted in heavy-quark effective theory, where one usually normalizes to the mass, i.e. with our convention (3.2) v 2 = 1.
The integral in (3.1) is evaluated by first boosting the vectors into the center-of-mass frame of the dipole and by subsequently changing the angular integration variables to (an appropriate generalization of) rapidity. The reader who is not interested in the technical details of the calculation of (3.1) can skip the following discussion and move directly to Section 4.

Boost to the center-of-mass frame
In order to calculate the integral in (3.1), it is convenient to first boost the dipole momenta from the laboratory frame into the center-of-mass frame of the dipole. To construct the relevant boost we use a form of the Lorentz transformation introduced by Householder [ where ∆ µ = n µ −ñ µ is the difference of two light-like vectors n µ andñ µ . One immediately verifies that Λ µ ν (∆)n ν =ñ µ , so the transformation maps n µ →ñ µ . In addition, it is straightforward to check that Λ µ ρ Λ ρ ν = δ µ ν and det(Λ) = −1. The transformation (3.3) is easily implemented into a computer code and here we use it to construct a boost of two arbitrary time-like or light-like directions u i and u j into a frame where these momenta are back-to-back alongside the z-axis. The transformation is carried out in three steps, denoted by X µ ν , B µ ν and Z µ ν . We denote lab frame vectors p µ i in the three frames reached by each of the transformations as The transformation X rotates the total dipole three momentum such that it points along the x-axis. Then B boosts into the center-of-mass frame and the last step Z rotates the system so that the back-to-back vectors lie along the z-axis.
Let us now discuss the three transformations in turn. The sum of the momenta associated to the the two vectors p i = E i u i and p j = E j u j is (3.5) By using the transformation (3.3) one can find the rotation to a frame where the spatial component of the light-like vector n P ≡ (1, n P ) points along the x-axis. This rotation (more precisely a rotation with parity inversion since det(Λ) = −1) is defined as Consequently, by applying the rotation X µ ν to the total momentum one findš The Lorentz transformation needed to obtain two back-to-back vectorsũ i andũ j from the original vectors u i and u j is now a boost along the x-axis. The corresponding transformation in matrix form is where β was introduced in (3.5) and γ = 1/ 1 − β 2 . Consequently, the two vectors are in a back-to-back configuration, i.e.
Finally, it is convenient to apply a last rotation in order to align the vectors along the z-axis. This can be achieved by employing again a Lorentz transformation of the form described in (3.3). In particular one can define In conclusion, the complete Lorentz transformation of any vector from the lab frame into a frame where the vectors u i = p i /E i and u j = p j /E j are back to back and aligned along the z-axis is (3.14) One finds that det(L) = 1, since L is the product of one proper and two improper transformations.

Evaluation of the angular integral
After applying the Lorentz transformation L and normalizing the vectors according to (3.2), one finds In this frame, one can write the integral over the dipole as where the light-like momentum n k in the center-of-mass system is parameterized as and y indicates the rapidity-like quantity The boundaries of the rapidity integration are In the massless limits, β i and/or β j become equal to 1 and y max and/or y min go to infinity. In that case, the collinear divergence in the integral (3.16) needs to be regularized. To this end, we apply a hard cutoff |y| < y cut in numerical computations, in addition to the constraints (3.19). We then verify that the physical cross sections are cutoff independent. The specific form of the collinear cutoff we use in our code is given in Appendix A, see (A.14). A discussion of different cutoffs can be found in Appendix A of [9]. The integral over the monopoles gives Combining the monopole and dipole contributions, the final result for the virtual correction reads with δ v i = 1 if v i is a time-like direction and zero otherwise. Note that the integration boundaries y max and y min depend on β i and β j , see (3.19).

The real emission corrections
are evaluated using Monte Carlo methods by randomly choosing a value of y and φ in the integrand of (3.16). The factor inside the bracket in (3.22) is a positive weight factor, as shown below. To see whether a given real-emission vector is inside the jet region, one transforms the vector n k back to the laboratory frame by using the inverse transformation to L given in (3.14).

Positive definiteness of W k ij
We now show that the weight factor in (3.22) is positive. This is done most conveniently in the center-of-mass frame. When written in terms of scalar products, the factor W k ij reads To see that this expression is indeed non-negative, one replaces the scalar products by By inserting the relations in (3.24) in (3.23) one finds (3.25) Consequently, the factor W k ij in (3.23) is always larger than or equal to zero.

Emissions from massive partons and non-global logarithms
In Section 2 we have shown that in the large-N c limit the monopole contributions can be absorbed into the dipole terms by replacing the usual dipole emitter W k ij given in (1.1) by the modified emitter W k ij introduced in (2.17). It is interesting to compare the massless and massive cases to illustrate the dead cone effect [11][12][13][14] mentioned in the introduction. In Figure 4 we plot the real-emission integrand (3.25) multiplied by the measure sin θ as a function of θ. The plot shows the collinear divergences at θ = 0 and θ = π, which are present in the massless case β i = β j = 1, while the massive integrand vanishes at the end points. One also observes that the monopole contribution significantly reduces the radiation, compared to the pure dipole contribution shown by the dotted line in the plot.
To see what effect the mass has on the size of non-global corrections, we consider the gap fraction in e + e − collisions. To define a gap region, we fix a direction n for each event and impose a veto E tot < Q 0 on radiation outside a cone around this direction. We then define the rapidity of an emission with momentum k as An emission is outside the cone, i.e. inside the gap region, if |y| < y max and we define the gap fraction as For massless final-state quarks e + e − → qq, one has to ensure that the reference vector n is chosen such that radiation collinear to the original partons is included to obtain a collinear safe cross section. To do so, one uses for n the thrust axis or the direction obtained from running a jet algorithm on the events. In the massive case e + e − → tt, on the other hand, we are completely free to choose the reference vector and we compare results obtained when choosing n collinear or perpendicular to the top-quark direction.
To study the contribution from the first two emissions, we expand in the evolution time t, which is directly related to Q 0 , see (2.8). The coefficients of the expansion can be obtained by iterating the one-loop anomalous dimension which determines the evolution factor (2.4) at LL. Following the steps outlined in Section 5.2 of [8] for the massless case, one finds The real-emission parts R m of the anomalous dimension in (2.16) generate an additional parton and the symbol⊗ indicates the integral over its direction. The angular brackets denote the normalized color trace, which in the large-N c limit reduces to the trivial trace GL are shown with dashed lines, the non-global parts S NGL using solid lines. The black lines in both panels are identical and correspond to radiation from a massive dipole with β i = β j = 0.5 and a reference vector n along the direction of the massive quarks. Left panel: Comparison to the massless case (blue lines). Note that the massless coefficients have been divided by 10 to make their size similar to the massive ones. Right panel: Comparison to the same β i = β j = 0.5 dipole with n perpendicular to the massive quarks (red lines). 1 = 1. Let us first discuss S (1) . We label the initial hard partons as 1 and 2 and the newly emitted gluon as 3. Then where we introduced the short-hand notation The virtual correction V m given in (2.15) has opposite sign and includes an integral over the entire solid angle. Combining it with the real-emission part, one finds that where 3 out = 1 − 3 in . The dipole structure after the first emission is (q, g, q) = (1, 3, 2). To be consistent with the notation in the anomalous dimensions (2.15) and (2.16) one should relabel the particles after the emission as (1, 2, 3), but we prefer to keep the original labels so that the neighboring dipoles in the second step are (1, 3) and (3,2), and We can rewrite all terms appearing in (4.4) in terms of angular integrals and combine real and virtual parts as we did in (4.7). This leads to the two-loop result NGL + S in agreement with the results given in [8] for the massless case. The global part of S (2) is just one half of the squared one-loop contribution, while the non-global piece has a more complicated structure that arises from the emission of a second gluon from the one produced in the first emission. Figure 5 shows the two-loop coefficients in different situations. In the left plot, where the cone is chosen along the direction of the original dipole, we compare the massive case with β = 1/2 to the massless one. In the massless case (shown in blue) S (1) ∝ y max so that the global part increases quadratically as y max is increased. In the massive case (shown in black), on the other hand, the radiation stops as the gluon becomes collinear to the quark so that the global part of the gap fraction goes to a constant as y max becomes large. Interestingly, the non-global part becomes constant as y max → ∞ in the massless case, while it vanishes for a non-zero mass. The radiation from a massless dipole is much larger than the one in the massive case; indeed it was necessary to divide the massless two-loop coefficients by a factor of ten to make them similar in size to the massive ones in the figure. In the right panel of Figure 5 we check how much of a difference the choice of the cone vector n makes. The red curves show the result when n is chosen perpendicular to the direction of the massive quarks instead of collinear to them. In this case, the massive quarks lie in the middle of the gap region. We observe that the size of the two-loop coefficients for the two choices of n is quite similar.
Having discussed the two-loop corrections, it is interesting to see how the fixed-order expansions compare to the full resummed result. In Figure 6 we show the result of the LL resummation of the gap fraction starting with a single dipole in the center of mass along the cone axis n for a gap with maximal rapidity y max = 0.8. The left plot shows the result for a massive dipole (β = 1/2) while the one on the right starts with a massless one (β = 1). Along with the full LL result, we also plot its NLO and NNLO expansion. The point made above is fully confirmed; the radiation from a massless dipole is much stronger than from a massive dipole. In fact, both the one-loop and two-loop coefficients are an order of magnitude larger for the massless case than for the massive one. The figure shows the gap fraction as a function of t. The relation among t and Q 0 depends on the value of Q. We stress that the larger values of t in the figure correspond to very small values of Q 0 . Indeed, for Q 1 = 1 TeV, t ≥ 0.1 corresponds to Q 0 1 GeV.

Resummation of tt production with veto on central jets
In this section, the formalism is applied to the resummation of non-global logarithms in a cross section involving soft radiation from top quarks. We consider tt production at the LHC with a veto on additional central jet activity as measured by ATLAS [15]. This measurement was performed to test the modeling of soft radiation from top quarks in parton shower Monte Carlo codes and is therefore well suited to study resummation effects.
In the measurement ATLAS considers events with at least two energetic b-jets, oppositesign leptons and missing energy, subject to a set of selection requirements designed to enhance the tt signal and reject background. In detail, the imposed cuts are as follows: Two of the b-jets must have p T > 25 GeV, |y| < 2.4 and ∆R(j, l) > 0.4, where ∆R(x, y) = (∆φ(x, y)) 2 + (∆η(x, y)) 2 with ∆φ(x, y) and ∆η(x, y) being the difference of the azimuthal angle and the rapidity of particles x and y. The opposite charged leptons must fulfill the usual ATLAS cuts: for muons p T > 20 GeV, |η| < 2.5 and for electrons p T > 25 GeV, |η| < 2.47. If the two leptons are of the same flavor, one imposes that their invariant mass is not too small, m > 15 GeV, and not near the Z-resonance, |m − m Z | > 10 GeV. In addition one requires missing E miss T > 40 GeV. In the mixedflavor µe-channel, one instead imposes that H T > 130 GeV, where H T is the scalar sum of the visible transverse momenta.
Starting with this event sample, ATLAS then defines a gap region as depicted in Figure 1. The gap consists of rapidity intervals y min < |y| < y max , but the bottom-tagged jets are removed from the gap region. In [15], four rapidity regions with various y min and y max are measured. We will focus on the two regions with gap regions |y| < 0.8 and |y| < 2.1.
For a given region, the gap fraction is defined as the fraction of events which do not involve a jet with transverse momentum above Q 0 in the gap. The luminosity drops out in the ratio so that the gap fraction is the ratio of the corresponding cross sections, which are both computed in the presence of the selection cuts discussed above, as defined in (4.2).
For our fixed-order predictions we use MadGraph5_aMC@NLO [26] and the Les-Houches Event (LHE) files produced by this code are taken as an input for our resummation code. We use NNPDF2.3 leading-order PDF sets, with α s (M Z ) = 0.130 [28]. For the fixedorder prediction of the gap fraction, we use the relation Up to corrections of O(α 2 s ), we can use lowest-order cross sections in this formula. We obtain these by generating tree-level events for tt and ttg with MadGraph5_aMC@NLO.
To be able to impose the selection cuts, specifically the exclusion of the bottom-tagged jets, we let the tt pair decay into leptons and a bb pair. The b andb are then acting as centers of a jet with size R = 0.4 in the plane of azimuthal angle and rapidity. Therefore, a particle q belongs to the gap region, if In the plots of this section the fixed-order predictions for the cross section are shown in green. As in any multi-scale problem, it is not clear what default value one should use for the renormalization and factorization scales. The average partonic center-of-mass energy √ŝ for tree-level tt events at the LHC with √ s = 7 TeV in the presence of the ATLAS selection cuts is about 520 GeV, which is about three times bigger than the top-quark mass and significantly larger than Q 0 , the lowest scale in the problem. We use an intermediate value µ r = µ f = 2m t as the default choice for the renormalization and factorization scales, but one could argue that the relevant scale for α s in the ratio in (5.1) is a lower value µ r ∼ µ s ∼ Q 0 since the factors of α s associated with the production cross section drop out in the ratio and only the coupling constant associated with the soft gluon emission remains. Indeed, choosing a lower scale would somewhat improve the agreement of the fixed-order prediction with data. The fixed-order uncertainty bands in the plots come from varying the scales µ r and µ f by factors of two around their default values while imposing 1/2 ≤ µ r /µ f ≤ 2, i.e. we are using the 7-point method to get the scale bands. Looking at the predictions for different scale values, we observe that the largest variations arise when both scales are simultaneously varied up or down. The fixed-order scale bands are fairly narrow. While the cross sections themselves have a relatively large scale uncertainty, most of it drops out in the ratio in (5.1).
Let us now turn to the resummation, which is performed on the basis of the LHE files for tt production. The shower code reads out the momenta of the top quarks and the initial-state particles in order to obtain the directions of the initial Wilson lines. These, together with the large-N c color dipole structure provided in the event file are the starting point for the shower, which then emits gluons until an emission goes into the gap. The value of the evolution time t ≡ t(µ h , µ s ) in (2.8) is later translated into a value of Q 0 , the scale associated with the emission. The shower also calculates the angular integrals in (2.22) as it evolves from the hard to the soft scale.
As the default hard scale we use µ h = Q 1 = 150 GeV which was calculated using (2.1) with an average Q = √ŝ ≈ 520 GeV. The soft scale µ s should be chosen to be of the order of Q 0 . However, we want to switch off the resummation at larger Q 0 values where we enter the fixed-order regime. To this end, we use a profile function which switches off the resummation for Q 0 → Q max . We choose Q max = µ h = 150 GeV and use the same functional form as in [10], namely whereQ 0 = Q 0 /Q max . The profile function is constructed such that µ s → x s Q 0 for Q 0 → 0 and that µ s → µ h forQ 0 → 1. The higher-power terms in the denominator are chosen such that the first few derivatives atQ 0 = 1 vanish and the parameter x s = { 1 2 , 1, 2} is used for scale variation. Beyond Q 0 = Q max all resummation effects are switched off and only the fixed-order prediction remains.
The results of the resummation (blue) along with its fixed-order expansion to the second order (brown and pink) as well as the fixed-order calculation results (green) are given in Figure 7. From this plot, one can clearly see that the difference between the LL result (blue) and its first-order expansion (brown) is moderate for the gap region with |y| < 0.8, while this difference is very large when the gap covers the interval |y| < 2.1. The effect of the radiation from the top decays is not negligible and reduces the gap fraction. By leaving it out (dashed blue line) one obtains a gap fraction that is sizeably larger especially at low Q 0 .
The ultimate goal of this section is to match the NLO calculation to the LL resummation in order to obtain LL+NLO predictions. The size of the difference between the LL result and its first order expansion, which we denote by LL@NLO, is relevant for the matching procedure, as discussed below.
When expanding the LL resummed result, we also expand the evolution time the envelope. It turns out that the variation of the soft scale is dominant throughout the plot.
There are two common schemes to combine the resummed results and fixed-order predictions, namely additive and multiplicative matching. In the additive matching scheme, one simply adds the LL gap fraction to the NLO prediction and subtracts the one-loop expansion of the LL gap fraction to avoid double counting Predictions obtained with the additive matching scheme (red) are shown in Figure 8 together with the NLO fixed-order results and data from the ATLAS measurement. In the multiplicative matching scheme, one exponentiates the matching corrections and multiplies the exponential by the LL result The results obtained by means of multiplicative matching are shown in red in Figure  9. Multiplicative matching exponentiates the entire first emission, which is similar in spirit to what is done in the POWHEG method [29,30]. The ATLAS paper compared their measurements to NLO results matched to parton showers using POWHEG and also MC@NLO [31,32]. Both schemes reproduce the data to better than 5%, POWHEG is typically even within 1%-2% of the measurement. The ATLAS paper does not provide the uncertainties of the theory prediction to which they compare, but we would expect them to be similar in size to the NLO uncertainty bands in our plots. In Section 4 it was shown that the radiation from massless legs is numerically much larger than from massive ones. Consequently, we expect that, in order to get a good description of the gap fraction, the |y|<2.1 Figure 9. Same as Figure 8, except that the multiplicative matching scheme was adopted.
modeling of the initial-state radiation is the most important effect. For this reason, it is not clear to us if a comparison to the ATLAS data provides a sufficiently stringent test of the description of soft radiation from massive quarks in a parton shower.
One observes that the additive matching scheme works well for the gap region |y| < 0.8 and actually mildly improves the agreement of central value with the data. However, for the case in which the gap region is |y| < 2.1, the predictions obtained with additive matching become unphysical for small values of Q 0 . This is not surprising, since the higher-order emissions are enhanced by factors of the gap size ∆y. If these rapidity logarithms become larger, they must be resummed. The formalism to carry out this resummation exists [7,8] but we do not implement it in the present work.
The multiplicative matching leads to better results since the matched gap fraction correctly vanishes for Q 0 → 0, as the resummed result does. Predictions obtained by means of multiplicative matching are shown in Figure 9, which shows that they are in good agreement with the experimental data, within the large scale uncertainty bands. To reduce these, it would be important to go to higher logarithmic accuracy, or to at least include higher-order corrections to the hard and soft functions, as it was done in the massless case [10].
In order to compare predictions to the Run I ATLAS measurement [15], all calculations were carried out at √ s = 7 TeV. For the tree-level top production process at √ s = 13 TeV, one finds that the average partonic center-of-mass energy is Q ≈ 550 GeV, which translates into Q 1 ≈ 170 GeV, only 20 GeV higher than at 7 TeV. Consequently, we conclude that the result for the gap fraction at √ s = 13 TeV would be quite similar to the ones at Run I.

Conclusion
In this paper, we have developed the necessary formalism to carry out the resummation of non-global logarithms for processes involving massive quarks. More specifically, we discussed how the parton shower approach needs to be modified to go beyond the high-energy limit, implemented those changes and then compared the radiation patterns of massive and massless partons. As an application, we have performed the leading logarithmic resummation of the cross section for tt production with a veto on central jet activity, an observable measured at Run I of the LHC. Soft radiation has a well-known eikonal form, independent of the mass of the emitting parton. However, in the massive case, the velocity vector of the emitting parton is time-like rather than light-like. This fact makes the kinematics of the process more complicated. A second important difference to the massless case is that for massive emitters one needs to account for monopole radiation. The radiator W k ij describing an emission between legs i and j is nonzero for i = j for massive legs; therefore the parton shower must also include radiation from a single leg, despite the fact that it is a purely soft shower. We have shown that in the large-N c limit this radiation can be absorbed into the dipoles by replacing the usual radiator with a modified one, indicated with W k ij . The monopole radiation has a negative relative sign with respect to the dipole contribution, but the total contribution W k ij remains positive. These properties make it straightforward to implement massive partons in the shower code that was previously developed for the emission from massless quarks [9,10].
Comparing the two cases, we observe that the massive dipole radiator is numerically significantly smaller than the massless one and that the radiation is further reduced by the monopole terms. For example, when analyzing the fixed-order expansion of the leading logarithmic resummation for a gap in the central rapidity region of size ∆y = 1.6, both the one-loop and the two-loop coefficients are an order of magnitude larger for a massless dipole compared to one with two massive legs with β = 0.5 each.
ATLAS measured the gap fraction in tt production with a veto on central jet activity [15]. This provides an interesting test case for the computational framework developed here. However, to compare to experiment we also need to account for radiation from the top-quark decay. To do so, we work in the narrow-width approximation in which the process factors into production and decay and then apply the parton shower to all color dipoles associated to the tt production as well as to the dipoles associated with the decay of the tt pair. The predictions that we compare to the ATLAS measurements are obtained by matching the LL resummed result to the NLO fixed-order computation of the gap fraction. There are two schemes commonly used to combine resummed and fixed order results: additive and multiplicative matching. For small gap sizes ∆y both schemes give similar results while for larger gaps the additive matching yields unphysical gap fractions. The problems for large gap sizes are not unexpected since the higher order corrections (and also the power-suppressed terms added in the matching) are enhanced by ∆y, i.e. by collinear logarithms. If these logarithms become large they must be resummed as well. The formalism necessary to implement this resummation exists [7,8], but the corresponding calculation is beyond the scope of the present paper.
In the present work, we resummed the leading non-global logarithms. In order to go to higher logarithmic accuracy, one needs to include the one-loop corrections to both the hard and the soft function, the tree-level result of the hard function with one additional emission, and to evolve with the two-loop anomalous dimension matrix. In the massless case, calculations including the first three ingredients listed above were recently presented [10]. Work on the final ingredient, the two-loop anomalous dimension matrix, is ongoing.
In this paper, we have extended the resummation of non-global observables to processes involving massive partons in the large-N c limit. Obviously, it would be desirable to extend the formalism to include logarithmic corrections beyond the large-N c limit. This would be especially interesting since Glauber phase effects then start to play a role in hadronic collisions. There are a few first finite-N c results in [33,34] based on a different formalism [35] and there is a considerable amount of ongoing work focused on the inclusion of subleading color effects into parton showers [36][37][38][39][40][41][42][43][44], but a full implementation of all subleading-color effects, in particular Glauber phases, is not yet available.
An understanding of non-global logarithms could prove useful in the context of the topquark mass determination. Given the complicated structure of these types of logarithms and our limited ability to perform all-order resummations, it is of course desirable to avoid them in the context of precision physics. On the other hand, to maximize sensitivity to the top-quark mass, jet observables are preferable to inclusive cross sections. It has been proposed to use jet substructure techniques such as grooming to reduce the sensitivity to soft radiation [45] and a factorization theorem implementing grooming has been put forward [46]. These techniques can reduce the size of non-global logarithms, and our approach could be used to get a better understanding of the remaining effects and their uncertainty.

A Details of the Monte Carlo algorithm
This appendix describes in detail the algorithm used to obtain the results presented in Section 5. We will start with a sample tree-level file and will then show how it is processed step by step by our code. This level of detail is not necessary for most readers, but should be useful for someone implementing a similar shower. It can also serve as a documentation of our code (written in Python), which we plan to make public in the future.
The starting point for the LL resummation algorithm is a Les Houches Event File (LHEF) [47] for the hard process produced using MadGraph5_aMC@NLO [26]. The generated process is the collision of two protons with center-of-mass energy √ s = 7 TeV producing a tt pair. Each top quark in the pair decays in a bottom quark and a Wboson. The latter is required to decay leptonically. In this way, the final state includes a bb pair, two leptons and two neutrinos. In MadGraph5_aMC@NLO syntax, the process is generated by the following command: generate p p > t t~> vl l+ vl~l-b bĨ n tt production at leading order, the partonic initial state includes either two incoming gluons (gluon fusion channel) or two incoming quarks (quark annihilation channel). Mad-Graph5_aMC@NLO computes the cross section for N c = 3 and then randomly assigns one of the possible large-N c color structures to each tree-level event so that it can be analyzed by a parton shower. The large-N c color structure is given by a set of dipoles, as illustrated in Figure 2. An event consists of four or five dipoles in total, namely two (quark annihilation) or three (gluon fusion) dipoles associated to the production of the top pair and two dipoles from the radiation in the decay to bottom quarks.
In narrow-width approximation, the amplitudes squared factorize into production and decay. We will separately compute the emissions from production and decay and obtain the cross section as a product where σ Event 0 is the Born-level event weight supplied by MadGraph5_aMC@NLO. The factors R tt (t), R t→b (t) and Rt ,→b (t) are computed by showering the color dipoles arising in the production and decay process. The product form (A.1) holds on the level of the squared amplitudes in the large-N c limit, but the observable Q 0 , the energy inside the veto region, is additive and the cross section will therefore be a convolution of the different pieces, not simply a product. However, at LL accuracy, the cross section is Q 0 -independent as it only depends on t ≡ t(µ h , µ s ) and the convolution then reduces to the product in (A.1).
Below, we illustrate the parton shower using the production process R tt (t), but we run exactly the same shower for the dipoles in the decay. We consider an event in the gluon fusion channel to discuss the showering process in detail in the following, since this is the most involved case. The dipole structure of this event is shown in Figure 10. We could separately shower each of the three dipoles, but it is more efficient to treat the event as one dipole with two intermediate gluons, see below. The form of the shower for R tt (t) is then similar to (2.22) except that the process starts with four partons where σ tt 0 is the Born-level production cross section, and H i ≡ H i . In the quark annihilation channel, the two dipoles need to be showered separately.
The remainder of this appendix is organized as follows. The set up of the shower is discussed in Section A.1 by looking at an explicit example. The shower procedure is discussed in A.2, also in this case an explicit event is considered as an example. Finally, a brief outline of the algorithm is provided in Section A.3.
The final state leptons in the event must satisfy the cuts listed in Table 1 of [15]. The momenta of the bb pair are needed, because they define the direction of the b-jets, which are cut out of the gap region (or veto region), see (5.2) and Figure 1.
In the following we illustrate the shower algorithm with the tt production process. The additional dipoles arising from the top-quark decay can be showered exactly in the same way. In the end, it is necessary to multiply the results of each shower to get the complete result of each event, see (A.1). The showering process of the dipoles associated to the top-pair production starts by selecting the momenta of the initial-state partons and the momenta of the tt pair. These momenta are stored in dipoles according to the large-N c color information assigned by MadGraph5_aMC@NLO. In the sample event depicted in Figure 10, the color index associated to the top quark is c1: 501, while the color indices associated to the gluon in the second line of the list above (g 2 ) are c1: 501, c2: 503. The color indices of the gluon in the first line of the list (g 1 ) are c1: 503, c2: 502 and the color index of the anti-quark is c2: 502. The color indices indicate which lines are color connected and the shower algorithm orders the particles in such a way that equal color indices are adjacent to each other; so that the list of the color indices is 501,501,503,503,502,502. Therefore the ordering of the particles is (t,g 2 ,g 1 ,t) which represents a dipole with two intermediate gluons.
Since the algorithm only requires information about the direction of the particles, we normalize the components of the momenta to their energy. These normalized momenta are stored in an array such that each adjacent pair of vectors represents a color dipole.
We then calculate the virtual correction of each dipole using (3.21) by boosting the two vectors in each dipole into a frame where they are back-to-back, as explained in Section 3, and by subsequently evaluating the velocities β i , β j . Each dipole contributes to the virtual corrections a factor where N c = 3.
Let us illustrate the process by explicitly calculating the virtual corrections for the first dipole in the list, V 12 , according to the method put forward in Section 3. The first step consists in boosting the two normalized momenta in the rest frame. The sum of the two normalized momenta is This vector is then aligned to the x-axis by means of the matrix withǓ having vanishing y and z components by construction. with g = diag(1, −1, −1, −1). Let us briefly explain the value of y cut which cuts off the collinear divergence arising for massless partons. Following (A.1) of [9], we obtain the value of the cut-off by imposing a rapidity cut η cut in the lab frame and then computing the value of y cut that this corresponds to in the center-of-mass frame. This leads to the somewhat complicated expression y cut = ln cos θ 2 + cos 2 θ 2 + sin 2 θ 2 e 2ηcut , (A.14) where θ is the angle between the two vectors forming the dipole in the lab frame. One immediately sees that y cut = η cut for back-to-back vectors θ = π. The value y cut for the dipole under consideration is obtained after setting η cut = 6. Since the top quark is massive and the gluon is massless, one has δ v 1 = 1 and δ v 2 = 0 in calculating the contribution of this dipole to the virtual corrections (A.4). The virtual correction associated to the first dipole in the event is given by Let us now analyze the individual terms in (A.16). The first termĤ 4 (t) denotes the evolution from the hard scale to the low scale without any radiation from any of the dipoles (compare to the first line of (2.21)) and is given bŷ since by definitionĤ 4 (0) = R tt (0) = 1. At this stage it is convenient to define the probability distribution such thatĤ 4 (t) = P(V 4 , t).
The functionĤ 5 (t) consists of the initial four hard partons plus one additional parton emitted by any of the dipoles at an evolution time earlier than t. In the large-N c limit, the new emission occurs from any one of the three dipoles in the list so that we get three termsĤ The quantity R 5 ij corresponds to the real correction factor as given in (3.22) when the fifth parton is emitted in the direction n 5 from the dipole of legs i and j To bring (A.20) into a form suitable for Monte Carlo implementation, we now strategically insert factors of one. We rewrite the first term on the right-hand side of the first line in the formĤ The integral dΩ 5 /4π, over the direction of the emission in (A.16), is evaluated by Monte Carlo methods. The factor V 12 /V 4 is the weight of the dipole (12) in the total virtual correction V 4 and is interpreted as the probability of having an emission from the dipole (12). For the case ofĤ The contribution of terms arising fromĤ       This iterative procedure can be repeated to calculate all of the termsĤ i (t).
In the parton shower algorithm, this procedure is implemented as follows. For the numerical example, we again consider the event already set up for the showering in Appendix A.1. At first, the shower is initialized as follows: The initial weight w of an individual event is the inverse of the number of showerings of a tree-level event in the LHE file, n sh ; additional weight factors arising from integrands such as the one in (A.23) are discussed below. At first it is necessary to randomly generate a time step ∆t according to the probability density P(V tot , ∆t). We generate random time steps ∆t according to this distribution by taking the cumulant u ≡ P/V tot ∈ [0, 1], inverting it ∆t = − ln u V tot , (A. 30) and using an equally distributed random variable u ∈ [0, 1]. For our sample event, we choose u = 0.5 for illustration, which yields ∆t = 0.00252. To account forĤ 4 (∆t) = P(V 4 , ∆t), we add a weight w into a histogram in a bin corresponding to the time t new = ∆t = 0.00252. Once the shower will be finished, this histogram will provide the gap fractionR tt (t). Next, we use (A.20) to iteratively computeĤ 5 (t) at a time t > t new . To do so, we interpret t new as the time at which one of the three original dipoles emits a parton. Looking at (A.20), one sees thatĤ 5 (t) has three terms. The first of these terms is given in (A.23); it involves a product of several factors. The first factor is the probability density P(V 4 , t ) which was taken into account when generating the time step ∆t. Showering multiple times, one gets a Monte Carlo approximation of the integral over dt . We then have the three factors R 5 12 /V 12 , V 12 /V 4 and V 4 /V 5 . The last two of these factors can be treated as probabilities since the value of these ratios is always in the interval [0, 1]. The factor R 5 12 /V 12 corresponds to the phase-space integral and will be treated as a weight. The very last factor P(V  To implement this, one can draw a random value u ∈ [0, 1]. Then, the i-th dipole is assumed to emit, if the cumulative sum of the virtual corrections of all the dipoles from 1 through i divided by the total virtual corrections is smaller than u. In the example under consideration, this means that if u < 0.190 = V 12 /V tot it is the first dipole that emits, if 0.190 < u < 0.715 = (V 12 + V 23 )/V tot the emission arises from the second dipole, and if 0.715 < u the third dipole emits. For the purposes of this discussion, let us assume that u = 0.1, such that first dipole emits. Now that the algorithm determined which dipole radiates, one can boost into the back-to-back frame of the selected dipole. The procedure to do this for the first dipole was already illustrated in Appendix A.1. In this frame, the algorithm draws two more random numbers, namely φ = 2πu φ , y = y min + u y (y max − y min ) , (A. 31) with u i ∈ [0, 1] and the integration boundaries y min and y max as given previously in (A.12). For illustration, assume that u φ = u y = 0.5, which yields φ = π and y = −1.754. With this input one then obtains the four-vector of the newly emitted parton as with n z = (e 2y − 1)/(β 1 e 2y + β 2 ) and n T = 1 − n 2 z .
The new vector is then boosted back to the lab frame by using the inverse boost matrix L −1 given in (A.13). Subsequently, the vector is normalized in such a way that the energy component is 1: (A.33) With the new vector n 5 , the algorithm evaluates the factor The quantity in (A.34) is strictly positive (as shown in Section 3.3), but it can not be treated as a probability, since it can exceed unity. 1 Therefore the algorithm accounts for it by modifying the weight factor We have to ensure that the real emission is not in the veto region, as imposed by the Θ in (n 5 )-function in (A.22). This is done by checking the conditions (5.2) with q = n 5 . We obtain ∆R(b, n 5 ) = 2.860 , ∆R(b, n 5 ) = 2.561 , y(n 5 ) = −3.496 , (A. 36) implying that the new emission fulfills all the conditions Θ in (n 5 ) and one can proceed to the next step of the algorithm. If the condition would have not been fulfilled, the shower would have been stopped at this point and the algorithm would have restarted at the beginning, after erasing all information on this showering other than the histogram entry forĤ 4 (t) at t = ∆t. Since the generated vector n 5 was not in the gap region, the algorithm continues by adding the new vector to the list of vectors in between the first and second vector of the list in (A.3). In addition, the algorithm updates the virtual correction list by replacing the virtual correction V 12 by V 15 and V 52 : 1 This arises because we include the monopoles as a weight factor into the dipole integral. Alternatively, one could integrate the full, modified dipole W k ij given in (3.25). The integral can be done and leads to a more complicated version of the rapidity variable (3.18). However, the inversion to the angle θ can then not be done analytically, in contrast to (3.18).
There is one last factor inĤ (1) 5 (t) that was not accounted for so far, namely the factor V 4 /V (1) 5 . This last factor can again be treated as a probability, as the new virtual corrections are always larger than the old ones. 2 This means that instead of multiplying the weight by this factor, one can again draw a random variable u ∈ (0, 1) and if u < 0.810 = V 4 /V (1) 5 = V old tot /V new tot , the algorithms continues with the generation of a time step starting at t = t new with weight w = w new . In the opposite case, the shower is stopped. One could also treat V 4 /V (1) 5 as a weight factor, but, due to the iterative nature of the shower, the weights of the individual steps multiply each other, leading to events with small weight which would render the shower inefficient.
The new time step is generated in exactly in the same way as before according to the new probability density P(V new tot , ∆t ) and completes the calculation ofĤ 5 (t) by writing the weight into the histogram at t = ∆t + ∆t . After this is done, the algorithm proceeds to calculateĤ 6 (t). Looking at (A.26), one sees that the same procedure described for the calculation ofĤ 5 (t) can be used, since it involves the same type of ingredients: a) An emitting dipole is chosen, each with probability of V ij /V tot .
b) The emission is generated and the factor R 6 ij /V ij is calculated and multiplied to the weight. c) If the emission is not in the veto region, the algorithm proceeds with probability V tot /V new tot and and generates a new time step using P(V new tot , ∆t ), which givesĤ 6 (t) with t = ∆t + ∆t + ∆t .
The iterative calculation of all theĤ i (t) with i > 6 can be carried out in the same way until one reaches the necessary maximal value for t determined by the lowest value of µ s ∼ Q 0 in the problem under consideration. Each showering generates several hard functions at successively larger times until it terminates. In the calculations presented in this work, we used an upper limit of t max = 0.1, which corresponds to µ s ≈ 0.75 GeV after applying the profile (5.3) and µ h = 150 GeV.

A.3 Parton shower algorithm
This section summarizes the different steps in the shower algorithm, which were discussed in detail in Appendix A.2 in the context of the showering of a particular event. The shower algorithm described below is applied n sh times to each tree-level event. In the following, we describe one such shower event. For the results presented in our paper, we used about 10 5 tree-level events and worked with n sh = 10 4 .
Step 0. Set up the shower Store all Wilson-line directions according to their color information into an array, and calculate the virtual corrections of each dipole {u} = {u 1 , u 2 , . . . , u m } , 2 Whether this is true depends on the form of the angular cutoff used in the shower, see [9]. where An expression for the rapidity values can be found in (3.19). Initiate the shower algorithm with the initial settings Step 1. Generate time step Step 2. Insert weight into histogram At this new time t, insert w into the histogram.
Step 3. Choose emitting dipole Randomly choose a dipole which emits the next emission, where each dipole with legs i and j emits with probability Step 4. Create emission Boost to the frame where the emitting legs are back-to-back along the z-axis and in that frame choose an emission direction n k by generating a random angle φ ∈ [0, 2π] and rapidity y ∈ [y min , y max ]. Boost n k back into the lab frame and normalize it as n k = (1, n k ). Update the event weight according to Step 5a. Emission not in veto region Go back to step 1 with a probability V tot V new tot := V 12 + · · · + V ij + · · · + V (m−1)m V 12 + · · · + V ik + V kj + · · · + V (m−1)m . (A.47) The algorithm restarts from step 1 with the new arrays {u}, {V } and V tot = V new tot . With the probability of 1 − (V tot /V new tot ) the shower is stopped. Of course one can also set an upper limit t max on t after which the shower stops.
Step 5b. Emission into veto region If the emission n k lands in the gap region, the shower stops.

B Fixed-order expansion of the LL result
In this appendix, we detail how one can obtain the coefficients S (1) and S (2) in the fixedorder expansion of the leading logarithmic resummation (4.3) from the parton shower. When extracting the fixed-order coefficients for a given tree-level event, we look at each dipole at t = 0 individually, calculate its expansion coefficients and then combine the results. Averaging the expansion coefficients of the individual events then gives the final result for the two-loop expansion of the resummed cross section. In the following we describe the computation for a single dipole.

B.1 One-loop coefficient
The one-loop coefficient of the fixed-order expansion may be extracted easily from the shower. To do so we write (4.7) as S (1) = −V 12 dΩ(n 3 ) 4π where 1 and 2 are the legs of the dipole and R 3 12 = 4N c W 3 12 . Please note that throughout this appendix we write out the appropriate Θ in (n k ) and Θ out (n k ) angular constraints and we do not include the factors Θ in (n k ) into the definition of R k ij as we did in (A.22). The factor R 3 12 /V 12 is produced by the shower, see Appendix A, and gives S (1) after multiplication by −V 12 . All that needs to be done is to account for the constraint Θ out (n 3 ) which ensures that the emission is in the veto region. From the first step of the parton shower one obtains n sh , if Θ out (n 3 ) = 1 , 0, otherwise. (B.2)

B.2 Two-loop coefficient
While the global part of the two-loop coefficient is just one half of the one-loop coefficient squared, the non-global part is much more involved mainly due to collinear divergences in individual terms in the integrand in (4.9). Rewriting the non-global part as one notices that collinear singularities arise in R 3 12 since the light-like direction n 3 is in the jet region (not in the veto region) and can become collinear to n 1 or n 2 . However, the full expression S (2) NGL is collinear finite since the terms multiplying R 3 12 vanish in both collinear limits: R 4 13 + R 4 23 − R 4 12 → 0 for n 3 → n 1 or n 3 → n 2 . (B.4) The two terms R 3 12 R 4 13 and R 3 12 R 4 23 have a simple interpretation in the parton shower. The first emission of the shower produces R 3 12 and results in the dipole configuration (1, 3, 2). The term R 4 13 arises when the second emission occurs in the dipole (1, 3), while R 4 23 corresponds to the emission from (3,2). However, the subtraction term R 4 12 does not arise in the parton shower. To include it, the factor R 4 12 can be split in two parts according to R 4 12 = R 4 12 θ(n 2 · n 3 − n 1 · n 3 ) + R 4 12 θ(n 1 · n 3 − n 2 · n 3 ) .

(B.5)
In this way, the first term is evaluated if the spatial angle between n 3 and direction n 1 is smaller than the one between n 3 and direction n 2 and is therefore used to cure the collinear singularity when n 3 → n 1 .The same argument holds for 1 ↔ 2 and removes the other divergence. One can thus write The terms in the second and third line of this expression are separately collinear finite and the factors multiplying the dipoles R 3 12 R 4 13 and R 3 12 R 4 23 in the two lines can be implemented as weight factors in the shower.
To implement (B.6) in the shower, we store the weight R 3 12 /V 12 of the first emission. We then go on to the second emission and check whether it is emitted by the dipole (n 1 , n 3 ) or (n 3 , n 2 ). We also check if n 3 , the direction of the first emission, is closer to n 1 or n 2 . With this information, one can then calculate the two-loop coefficient as a sum of weights (B.7) The weights for the two cases (i, j) = (1, 2) and (i, j) = (2, 1), corresponding to the second emission R 4 i3 , are obtained as follows , if Θ out (n 4 ) = 1 and n i · n 3 > n j · n 3 , n sh , if Θ out (n 4 ) = 1 and n i · n 3 < n j · n 3 , 0 , otherwise . (B.8) We have written s 2 in terms of factors V 12 /V 3 and V i3 /V 3 , with V 3 = V 13 + V 32 , which arise in the shower algorithm, analogous to (A.23). They represent the probability to continue the shower after the first emission and the probability to choose the dipole (n i , n 3 ) rather than (n j , n 3 ) for the second emission.