A new approach to QCD final-state evolution in processes with massive partons

Abstract

We present an algorithm for massive parton evolution which is based on the differentially accurate simulation of soft-gluon radiation by means of a non-trivial azimuthal angle dependence of the splitting functions. The kinematics mapping is chosen such as to to reflect the symmetry of the final state in soft-gluon radiation and collinear splitting processes. We compute the counterterms needed for a fully differential NLO matching and discuss the analytic structure of the parton shower in the NLL limit. We implement the new algorithm in the numerical code Alaric and present a first comparison to experimental data.

I Introduction

The production and evolution of massive partons are an important aspect of collider physics, and they play a particularly prominent role at the Large Hadron Collider at CERN. Key measurements and searches, such as $t\bar{t}H$ and double Higgs boson production, involve final states with many $b$-jets. The success of the LHC physics program therefore depends crucially on the modeling of heavy quark processes in the Monte-Carlo event generators used to link theory and experiment. With the high-luminosity phase of the LHC approaching fast, it is important to increase the precision of these tools in simulations involving massive partons.

Heavy quark and heavy-quark associated processes have been investigated in great detail, both from the perspective of fixed-order perturbative QCD and using resummation, see for example Mangano et al. (1992); Bonciani et al. (1998); Cacciari et al. (2004, 2012). Various proposals were made for the fully differential simulation in the context of particle-level Monte Carlo event generators Norrbin and Sjöstrand (2001); Gieseke et al. (2003); Schumann and Krauss (2008); Gehrmann-De Ridder et al. (2012). Recently, a new scheme was devised for including the evolution of massive quarks in the initial state of hadron-hadron and lepton-hadron collisions Krauss and Napoletano (2018). In this manuscript, we will introduce an algorithm for the final-state evolution and matching in heavy-quark processes, inspired by the recently proposed parton-shower model Alaric Herren et al. (2022). The soft components of the splitting functions are derived from the massive eikonal and are matched to the quasi-collinear limit using a partial fractioning technique. In contrast to the method of Catani et al. (2002); Schumann and Krauss (2008), we partial fraction the complete soft eikonal, leading to strictly positive splitting functions and thus keeping the numerical efficiency of the Monte-Carlo algorithm at a maximum. We also propose to use a kinematic mapping for the collinear splitting of gluons into quarks that treats the outgoing particles democratically. This algorithm can be extended to any purely collinear splitting (i.e.,  after subtracting any soft enhanced part of the splitting functions) while retaining the NLL precision of the evolution. Our algorithm does not account for the effect of spin correlations.

Multi-jet merging and matching of parton-shower simulations to NLO calculations in the context of heavy-quark production were discussed, for example, in Mangano et al. (2002); Frixione et al. (2003, 2007); Höche et al. (2013). The NLO matching is typically fairly involved, because of the complex structure and partly ambiguous definition of the infrared counterterms. In this publication, we compute the integrated counterterms for our new parton-shower model, making use of recent results for angular integrals in dimensional regularization Lyubovitskij et al. (2021). This calculation provides the remaining counterterms needed for the matching of the Alaric parton-shower model at NLO QCD. We will discuss the extension to initial-state radiation in a future publication.

This paper is structured as follows: Section II briefly reviews the construction of the Alaric parton-shower model and generalizes the discussion to massive particles. Section III introduces the different kinematics mappings. Section IV discusses the general form of the phase-space factorization and provides explicit results for processes with soft radiation and collinear splitting. The computation of integrated infrared counterterms is presented in Sec. V. Section VI discusses the impact of the kinematics mapping on sub-leading logarithms, and Sec. VII provides first numerical predictions for $e^{+}e^{-}\to$hadrons. Section VIII contains an outlook.

II Soft-collinear matching

We start the discussion by revisiting the singularity structure of $n$-parton QCD amplitudes in the infrared limits. If two partons, $i$ and $j$, become quasi-collinear, the squared amplitude factorizes as

$${}_{n}\langle 1,\ldots,n|1,\ldots,n\rangle_{n}=\sum_{\lambda,\lambda^{\prime}=\pm}\,{}_{n-1}{\Big{<}}1,\ldots,i\!\!\backslash(ij),\ldots,j\!\!\!\backslash,\ldots,n\Big{|}\frac{8\pi\alpha_{s}\,P^{\lambda\lambda^{\prime}}_{(ij)i}(z)}{(p_{i}+p_{j})^{2}-m_{ij}^{2}}\Big{|}1,\ldots,i\!\!\backslash(ij),\ldots,j\!\!\!\backslash,\ldots,n\Big{>}_{n-1}\;,$$ (1)

where the notation $i\!\!\backslash$ indicates that parton $i$ is removed from the original amplitude, and where $(ij)$ is the progenitor of partons $i$ and $j$. The functions $P^{\lambda\lambda^{\prime}}_{ab}(z)$ are the spin-dependent, massive DGLAP splitting functions, which depend on the momentum fraction $z$ of parton $i$ with respect to the mother parton, $(ij)$, and on the helicities $\lambda$ Dokshitzer (1977); Gribov and Lipatov (1972); Lipatov (1975); Altarelli and Parisi (1977); Catani et al. (2001, 2002). For the remainder of this work, we will use only the spin-averaged splitting functions.

In the limit that gluon $j$ becomes soft, the squared amplitude factorizes as Bassetto et al. (1983)

$${}_{n}\langle 1,\ldots,n|1,\ldots,n\rangle_{n}=-8\pi\alpha_{s}\sum_{i,k\neq i}\,{}_{n-1}\big{<}1,\ldots,j\!\!\!\backslash,\ldots,n\big{|}{\bf T}_{i}{\bf T}_{k}\,w_{ik,j}\big{|}1,\ldots,j\!\!\!\backslash,\ldots,n\big{>}_{n-1}\;,$$ (2)

where ${\bf T}_{i}$ and ${\bf T}_{k}$ are the color insertion operators defined in Bassetto et al. (1983); Catani and Seymour (1997). In the remainder of this section, we will focus on the eikonal factor, $w_{ik,j}$, for massive partons and how it can be re-written in a suitable form to match the spin-averaged splitting functions in the soft-collinear limit. The eikonal factor is given by

$$\begin{split}w_{ik,j}=&\;\frac{p_{i}p_{k}}{(p_{i}p_{j})(p_{j}p_{k})}-\frac{p_{i}^{2}/2}{(p_{i}p_{j})^{2}}-\frac{p_{k}^{2}/2}{(p_{k}p_{j})^{2}}\\ \end{split}$$ (3)

Following Refs. Webber (1986); Herren et al. (2022), we split Eq. (3) into an angular radiator function, and the gluon energy

$$w_{ik,j}=\frac{W_{ik,j}}{E_{j}^{2}}\;,\quad\text{where}\quad W_{ik,j}=\frac{1-v_{i}v_{k}\cos{\theta_{ik}}}{(1-v_{i}\cos{\theta_{ij}})(1-v_{k}\cos{\theta_{jk}})}-\frac{(1-v_{i}^{2})/2}{(1-v_{i}\cos{\theta_{ij}})^{2}}-\frac{(1-v_{k}^{2})/2}{(1-v_{k}\cos{\theta_{jk}})^{2}}\;.$$ (4)

The parton velocity is defined as $v=\sqrt{1-m^{2}/E^{2}}$ and is frame dependent. When we label velocities by particle indices in the following, it is implicit that they are computed in the frame of a time-like momentum $n^{\mu}$. In this reference frame they reduce to the relative velocity $v_{i}=v_{p_{i}n}$, where

$$v_{pq}=\sqrt{1-\frac{p^{2}q^{2}}{(pq)^{2}}}\;.$$ (5)

For convenience, we also introduce the velocity four-vector

$$l_{p,q}^{\mu}=\sqrt{q^{2}}\,\frac{p^{\mu}}{pq}\;.$$ (6)

When this vector is labeled by a single particle index only, it again refers to the four-velocity of that particle in the frame of $n^{\mu}$, for example $l_{i}^{\mu}=l_{p_{i},n}^{\mu}$. When partial fractioning Eq. (4), we aim at a positive definite result in order to maintain the interpretation of a probability density and, with it, the high efficiency of the unweighted event generation in a parton shower. Following the approach of Ref. Herren et al. (2022), we obtain

$$W_{ik,j}=\bar{W}_{ik,j}^{i}+\bar{W}_{ki,j}^{k}\;,\qquad\text{where}\qquad\bar{W}_{ik,j}^{i}=\frac{1-v_{k}\cos\theta_{jk}}{2-v_{i}\cos\theta_{ij}-v_{k}\cos\theta_{jk}}\,W_{ik,j}\;.$$ (7)

The partial fractioned angular radiator function can be written in a more convenient form using the velocity four-vectors. We find an expression that makes the matching to the $ij$-collinear sector manifest

$$\bar{W}_{ik,j}^{i}=\frac{1}{2l_{i}l_{j}}\left(\frac{l_{ik}^{2}}{l_{ik}l_{j}}-\frac{l_{i}^{2}}{l_{i}l_{j}}-\frac{l_{k}^{2}}{l_{k}l_{j}}\right)\;,\qquad\text{where}\qquad l_{ik}^{\mu}=l_{i}^{\mu}+l_{k}^{\mu}\;.$$ (8)

In the quasi-collinear limit for partons $i$ and $j$, we can write the eikonal factor in Eq. (3) as Catani et al. (2001, 2002)

$$w_{ik,j}\underset{m_{i}\varpropto p_{i}p_{j}}{\overset{i||j}{\longrightarrow}}w_{ik,j}^{\rm(coll)}(z)=\frac{1}{2p_{i}p_{j}}\left(\frac{2z}{1-z}-\frac{m_{i}^{2}}{p_{i}p_{j}}\right)\;,\qquad\text{where}\qquad z\underset{m_{i}\varpropto p_{i}p_{j}}{\overset{i||j}{\longrightarrow}}\frac{E_{i}}{E_{i}+E_{j}}\;.$$ (9)

This can be identified with the leading term (in $1-z$) of the DGLAP splitting functions $P_{aa}(z,\varepsilon)$, where

$$\begin{split}P_{qq}(z,\varepsilon)&=C_{F}\left(\frac{2z}{1-z}-\frac{m_{i}^{2}}{p_{i}p_{j}}+(1-\varepsilon)(1-z)\right)\;,\\ P_{gg}(z,\varepsilon)&=C_{A}\left(\frac{2z}{1-z}+z(1-z)\right)\;,\\ P_{gq}(z,\varepsilon)&=T_{R}\left(1-\frac{2z(1-z)}{1-\varepsilon}\right)\;.\end{split}$$ (10)

To match the soft to the collinear splitting functions, we therefore replace

$$\frac{P_{(ij)i}(z,\varepsilon)}{(p_{i}+p_{j})^{2}-m_{ij}^{2}}\to\frac{P_{(ij)i}(z,\varepsilon)}{(p_{i}+p_{j})^{2}-m_{ij}^{2}}+\delta_{(ij)i}\,{\bf T}_{i}^{2}\Bigg{[}\frac{\bar{W}_{ik,j}^{i}}{E_{j}^{2}}-w_{ik,j}^{\rm(coll)}(z)\Bigg{]}\;,$$ (11)

where the two contributions to the gluon splitting function are treated as two different radiators Höche and Prestel (2015). As in the massless case, this substitution introduces a dependence on a color spectator, $k$, whose momentum defines a direction independent of the direction of the collinear splitting Herren et al. (2022). In general, this implies that the splitting functions which were formerly dependent only on a momentum fraction along this direction, now acquire a dependence on the remaining two phase-space variables of the new parton. It is convenient to define the purely collinear remainder functions

$$\begin{split}C_{qq}(z,\varepsilon)&=C_{F}(1-\varepsilon)\left(1-z\right)\;,\\ C_{gg}(z,\varepsilon)&=C_{A}\,z(1-z)\;,\\ C_{gq}(z,\varepsilon)&=T_{R}\left(1-\frac{2z(1-z)}{1-\varepsilon}\right)\;.\end{split}$$ (12)

They can be implemented in the parton shower using collinear kinematics, while maintaining the overall NLL precision of the simulation. This will be discussed further in Sec. VI.

III Momentum mapping

Every parton shower algorithm requires a method to map the momenta of the Born process to a kinematical configuration after additional parton emission. This mapping is linked to the factorization of the differential phase-space element for a multi-parton configuration. Collinear safety is a basic requirement for every momentum mapping. In addition, a mapping is NLL safe if it preserves the topological features of radiation simulated previously Dasgupta et al. (2018, 2020). We will make this statement more precise in Sec. VI.

This section provides a generic momentum mapping for massive partons that is both collinear and NLL-safe. It is constructed for both initial- and final-state radiation, inspired by the identified particle dipole subtraction algorithm of Catani and Seymour (1997). In the massless limit, we reproduce this algorithm and thus agree with the existing parton-shower model Alaric Herren et al. (2022).

III.1 Radiation kinematics

We will first describe the kinematics mapping needed for soft evolution. This is sketched in Fig. 1. The momenta $\tilde{K}$ and $\tilde{p}_{ij}$ are mapped to the momenta $K$, $p_{i}$ and $p_{j}$, while $\tilde{p}_{k}$ acts as a spectator that defines the azimuthal direction. We assume that any of the particles $i$, $j$ and $k$ can be massive. The momentum $\tilde{K}$ can be chosen in a suitable way to reflect the dynamics of the process Herren et al. (2022) and absorbs the recoil in the splitting. We define the variables

$$\mu_{ij}^{2}=\frac{\tilde{p}_{ij}^{2}}{2\tilde{p}_{ij}\tilde{K}}\;,\qquad\bar{\mu}_{ij}^{2}=\frac{2\mu_{ij}^{2}}{1+v_{\tilde{p}_{ij}\tilde{K}}}\;,\qquad\kappa=\frac{\tilde{K}^{2}}{2\tilde{p}_{ij}\tilde{K}}\;,\qquad\bar{\kappa}=\frac{2\kappa}{1+v_{\tilde{p}_{ij}\tilde{K}}}\;,$$ (13)

and the analogous final-state variables $\mu_{i/j}^{2}=m_{i/j}^{2}/(2\tilde{p}_{ij}\tilde{K})$ and $\bar{\mu}_{i/j}^{2}=2\mu_{i/j}^{2}/(1+v_{\tilde{p}_{ij}\tilde{K}})$.
The scalar invariants after the splitting are defined in terms of the energy fraction, $z$ Catani and Seymour (1997)

$$2p_{i}n=z\,2\tilde{p}_{ij}\tilde{K}\;,\qquad n^{2}=\big{(}1-z+\kappa+\mu_{ij}^{2}-\mu_{i}^{2}\big{)}\,2\tilde{p}_{ij}\tilde{K}\;.$$ (14)

The momentum of the radiator after the emission is

$$p_{i}^{\mu}=\bar{z}\,\tilde{p}_{ij}^{\mu}+\frac{\mu_{i}^{2}-\bar{z}^{2}\mu_{ij}^{2}}{\bar{z}\,v_{\tilde{p}_{ij}\tilde{K}}}\left(\tilde{K}^{\mu}-\bar{\kappa}\,\tilde{p}_{ij}^{\mu}\right)\;,$$ (15)

where

$$\bar{z}=\frac{z+2\mu_{i}^{2}}{1+v_{\tilde{p}_{ij}\tilde{K}}+2\mu_{ij}^{2}}+\sqrt{\bigg{(}\frac{z+2\mu_{i}^{2}}{1+v_{\tilde{p}_{ij}\tilde{K}}+2\mu_{ij}^{2}}\bigg{)}^{2}-\frac{2\mu_{i}^{2}(1+\bar{\kappa})}{1+v_{\tilde{p}_{ij}\tilde{K}}+2\mu_{ij}^{2}}}\;.$$ (16)

We define the variable

$$\bar{v}=\frac{\displaystyle\frac{2\,v\,z}{1+v_{\tilde{p}_{ij}\tilde{K}}}-\frac{\bar{\mu}_{i}^{2}}{\bar{z}}\bigg{(}1-\bar{z}+\bar{\kappa}-\frac{\bar{\kappa}-\bar{\mu}_{j}^{2}}{\zeta}\bigg{)}}{\displaystyle\bar{z}\,-\frac{\bar{\mu}_{i}^{2}}{\bar{z}}\frac{1-\bar{z}+\bar{\kappa}}{\zeta}}\;,\qquad\text{where}\qquad\zeta=\frac{2}{1+v_{\tilde{p}_{ij}\tilde{K}}}\,\frac{1-z+\kappa+\mu_{ij}^{2}-\mu_{i}^{2}}{1-\bar{z}+\bar{\kappa}}\;,$$ (17)

and where $v=(p_{i}p_{j})/(p_{i}n)$ is defined as in the massless case Catani and Seymour (1997); Herren et al. (2022). In terms of these quantities, the gluon momentum is given by

$$p_{j}^{\mu}=\frac{1-z+\mu_{ij}^{2}-\mu_{i}^{2}+\mu_{j}^{2}}{v_{\tilde{p}_{ij}\tilde{K}}\zeta}\left(\tilde{p}_{ij}^{\mu}-\bar{\mu}_{ij}^{2}\tilde{K}^{\mu}\right)+\bar{v}\,\frac{1+v_{\tilde{p}_{ij}\tilde{K}}}{2v_{\tilde{p}_{ij}\tilde{K}}}\left[\left(\tilde{K}^{\mu}-\bar{\kappa}\tilde{p}_{ij}^{\mu}\right)-\frac{1-\bar{z}+\bar{\kappa}}{\zeta}\left(\tilde{p}_{ij}^{\mu}-\bar{\mu}_{ij}^{2}\tilde{K}^{\mu}\right)\right]+k_{\perp}^{\mu}\;.$$ (18)

where

$${\rm k}_{\perp}^{2}=\tilde{p}_{ij}\tilde{K}(1+v_{\tilde{p}_{ij}\tilde{K}})\,\bar{v}\left[\left(1-\frac{\bar{v}}{\zeta}\right)(1-\bar{z})-\frac{1-\zeta+\bar{v}}{\zeta}\,\bar{\kappa}+\frac{\bar{\mu}_{j}^{2}}{\zeta}\,\right]-m_{j}^{2}\;.$$ (19)

Figure 1: Sketch of the radiation kinematics in final-state evolution. See the main text for details. Note that $p_{k}$ is unaltered by the mapping and only acts as a reference for the azimuthal angle $\phi$ (cf. Eqs. (21) and (22)).

In order to determine a reference direction for the azimuthal angle $\phi=\arctan({\rm k}_{y}/{\rm k}_{x})$, we note that the soft radiation pattern must be correctly generated. To achieve this, we compose the transverse momentum as

$$k_{\perp}^{\mu}={\rm k}_{\perp}\left(\cos\phi\,\frac{n_{\perp}^{\mu}}{|n_{\perp}|}+\sin\phi\,\frac{l_{\perp}^{\,\mu}}{|l_{\perp}|}\right)\;,$$ (20)

where the reference axes $n_{\perp}$ and $l_{\perp}$ are given by the transverse projections

$$n_{\perp}^{\mu}=p_{k}^{\mu}-\frac{\big{(}p_{k}(\tilde{K}-\bar{\kappa}\tilde{p}_{ij})\big{)}\,\big{(}\tilde{p}_{ij}^{\mu}-\bar{\mu}_{ij}^{2}\tilde{K}^{\mu}\big{)}}{2\tilde{p}_{ij}\tilde{K}\,v_{\tilde{p}_{ij}\tilde{K}}^{2}/(1+v_{\tilde{p}_{ij}\tilde{K}})}-\frac{\big{(}p_{k}(\tilde{p}_{ij}-\bar{\mu}_{ij}^{2}\tilde{K})\big{)}\,\big{(}\tilde{K}^{\mu}-\bar{\kappa}\tilde{p}_{ij}^{\mu}\big{)}}{2\tilde{p}_{ij}\tilde{K}\,v_{\tilde{p}_{ij}\tilde{K}}^{2}/(1+v_{\tilde{p}_{ij}\tilde{K}})}\;,$$ (21)

and

$$l_{\perp}^{\,\mu}=\varepsilon^{\mu}_{\;\nu\rho\sigma}\,\tilde{p}_{ij}^{\,\nu}\,\big{(}\tilde{K}^{\rho}-\bar{\kappa}\tilde{p}_{ij}^{\,\rho}\big{)}\,n_{\perp}^{\sigma}\;.$$ (22)

Since the differential emission phase-space element is a Lorentz-invariant quantity, the azimuthal angle $\phi$ is Lorentz invariant Herren et al. (2022). This allows us to write the emission phase-space in a frame-independent way. Upon construction of the momenta $p_{i}$, $p_{j}$ and $K$, the momenta $\{p_{l}\}$ which are used to define $\tilde{K}$ are subjected to a Lorentz transformation,

$$p^{\mu}_{l}\rightarrow\Lambda^{\mu}_{\nu}(K,\tilde{K})p^{\nu}_{l},\quad{\rm where}\quad\Lambda^{\mu}_{\nu}=g^{\mu}_{\nu}-2\frac{(\tilde{K}+K)^{\mu}(\tilde{K}+K)_{\nu}}{(\tilde{K}+K)^{2}}+2\frac{\tilde{K}^{\mu}K_{\nu}}{K^{2}}$$ (23)

III.2 Splitting kinematics

Figure 2: Sketch of the splitting kinematics in final-state evolution. See the main text for details. Note that the unlabeled momentum is unaltered by the mapping, and only acts as a reference direction to define the azimuthal angle $\phi$.

In this section, we describe the kinematics for the implementation of the purely collinear components of the splitting functions in Eq. (12). This is sketched in Fig. 2. We make use of some of the notation in Catani et al. (2002), in particular

$$y=\frac{p_{i}p_{j}}{p_{i}p_{j}+p_{i}K+p_{j}K}\;,\qquad z=\frac{p_{i}K}{p_{i}K+p_{j}K}\;,$$ (24)

and we define the scaled masses

$$\mu_{i}^{2}=\frac{m_{i}^{2}}{2\tilde{p}_{ij}\tilde{K}}\;,\qquad\mu_{j}^{2}=\frac{m_{j}^{2}}{2\tilde{p}_{ij}\tilde{K}}\;,\qquad\kappa=\frac{\tilde{K}^{2}}{2\tilde{p}_{ij}\tilde{K}}\;,\qquad\bar{\kappa}=\frac{2\kappa}{1+v_{\tilde{p}_{ij},\tilde{K}}}\;.$$ (25)

We also introduce the following variables

$$\alpha_{ij}=y+(1-y)(\mu_{i}^{2}+\mu_{j}^{2})\;,\qquad\bar{\alpha}_{ij}=\frac{2\alpha_{ij}}{1+v_{\tilde{p}_{ij},\tilde{K}}}\;,\qquad\mu_{ij}^{2}=\frac{m_{ij}^{2}}{2\tilde{p}_{ij}\tilde{K}}\;,\qquad\bar{\mu}_{ij}^{2}=\frac{2\mu_{ij}^{2}}{1+v_{\tilde{p}_{ij},\tilde{K}}}\;.$$ (26)

In terms of the additional variables

$$\begin{split}z_{ij}=&\;\frac{1}{2\bar{\alpha}_{ij}}\Bigg{(}\frac{1+\bar{\alpha}_{ij}}{1+\bar{\kappa}}+\bar{\mu}_{ij}^{2}-\sqrt{\left(\frac{1-\bar{\alpha}_{ij}}{1+\bar{\kappa}}+\bar{\mu}_{ij}^{2}\right)^{2}-\frac{4\bar{\alpha}_{ij}\bar{\kappa}}{(1+\bar{\kappa})^{2}}}\;\Bigg{)}\;,\qquad\bar{z}_{ij}=\frac{2z_{ij}}{1+v_{\tilde{p}_{ij},\tilde{K}}}\;,\\ \bar{z}=&\;\frac{\displaystyle z\big{(}1-\alpha_{ij}+\mu_{ij}^{2}\big{)}-\big{(}\alpha_{ij}+\mu_{i}^{2}-\mu_{j}^{2}\big{)}\frac{\bar{z}_{ij}\kappa}{1-\bar{z}_{ij}\alpha_{ij}+\bar{\mu}_{ij}^{2}}}{\displaystyle\frac{1-\bar{z}_{ij}\alpha_{ij}+\bar{\mu}_{ij}^{2}}{\bar{z}_{ij}}-\frac{\bar{z}_{ij}\alpha_{ij}\kappa}{1-\bar{z}_{ij}\alpha_{ij}+\bar{\mu}_{ij}^{2}}}\;,\end{split}$$ (27)

the momenta after the splitting are given by

$$\begin{split}p_{i}^{\mu}=&\;\bar{z}\;\frac{\tilde{p}_{ij}^{\mu}-\bar{\mu}_{ij}^{2}\tilde{K}^{\mu}}{\bar{z}_{ij}\,v_{\tilde{p}_{ij},\tilde{K}}}\,+\Big{(}y(1-\bar{z})(1+\mu_{ij}^{2}-\mu_{i}^{2}-\mu_{j}^{2})-\bar{z}(\mu_{i}^{2}+\mu_{j}^{2})+2\mu_{i}^{2}\Big{)}\,\frac{\tilde{K}^{\mu}-\bar{\kappa}\,\tilde{p}_{ij}^{\mu}}{v_{\tilde{p}_{ij},\tilde{K}}/z_{ij}}+k_{\perp}^{\mu}\;,\\ p_{j}^{\mu}=&\;(1-\bar{z})\;\frac{\tilde{p}_{ij}^{\mu}-\bar{\mu}_{ij}^{2}\tilde{K}^{\mu}}{\bar{z}_{ij}\,v_{\tilde{p}_{ij},\tilde{K}}}\,+\Big{(}y\,\bar{z}\,(1+\mu_{ij}^{2}-\mu_{i}^{2}-\mu_{j}^{2})-(1-\bar{z})(\mu_{i}^{2}+\mu_{j}^{2})+2\mu_{j}^{2}\Big{)}\,\frac{\tilde{K}^{\mu}-\bar{\kappa}\,\tilde{p}_{ij}^{\mu}}{v_{\tilde{p}_{ij},\tilde{K}}/z_{ij}}-k_{\perp}^{\mu}\;.\end{split}$$ (28)

The transverse momentum squared is given by

$${\rm k}_{\perp}^{2}=2\tilde{p}_{ij}\tilde{K}\,\Big{[}\,\bar{z}(1-\bar{z})\alpha_{ij}-(1-\bar{z})\mu_{i}^{2}-\bar{z}\mu_{j}^{2}\,\Big{]}\;.$$ (29)

The construction of the transverse momentum vector proceeds as described in Sec. III.1. While the choice of the reference vector defining $n_{\perp}^{\mu}$ is in principle arbitrary, it can be made conveniently, e.g.,  to aid the implementation of spin correlations in collinear gluon splittings.

IV Phase space factorization

In this section, we discuss the factorization of the differential $n+1$ particle phase-space element into a differential $n$ particle phase-space element and the radiative phase space. We start from the generic four-dimensional expression for the initial-state momenta $p_{a}^{\mu}$ and $p_{b}^{\mu}$ and final-state momenta, $\{p_{1}^{\mu},\ldots,p_{i}^{\mu},\ldots,p_{j}^{\mu},\ldots,p_{n}^{\mu}\}$,

$$\begin{split}{\rm d}\Phi_{n}(p_{a},p_{b};p_{1},\ldots,p_{i},\ldots,p_{j},\ldots,p_{n})=&\;\left[\prod_{i=1}^{n}\frac{1}{(2\pi)^{3}}\frac{{\rm d}^{3}p_{i}}{2p_{i}^{0}}\right](2\pi)^{4}\delta^{(4)}\Big{(}p_{a}+p_{b}-\sum p_{i}\Big{)}\;.\end{split}$$ (30)

Replacing particles $i$ and $j$ with a combined “mother” particle $\widetilde{\imath\jmath}$, we can write the expression for the underlying Born phase space as

$$\begin{split}&{};\tilde{p}_{1},\ldots,\tilde{p}_{ij},\ldots,\backslash\!\!\!\!\tilde{p}_{j},\ldots,\tilde{p}_{n})=\\ &\qquad\;\Bigg{[}\prod_{\begin{subarray}{c}k=1\\ k\neq i,j\end{subarray}}^{n}\frac{1}{(2\pi)^{3}}\frac{{\rm d}^{3}\tilde{p}_{i}}{2\tilde{p}_{i}^{0}}\Bigg{]}\frac{1}{(2\pi)^{3}}\frac{{\rm d}^{3}\tilde{p}_{ij}}{2\tilde{p}_{ij}^{0}}(2\pi)^{4}\delta^{(4)}\Big{(}\tilde{p}_{a}+\tilde{p}_{b}-\sum_{k\neq i,j}\tilde{p}_{k}-\tilde{p}_{ij}\Big{)}\;.\end{split}$$ (31)

The ratio of differential phase-space elements after and before the mapping is defined as the differential phase-space element for the one-particle emission process

$${\rm d}\Phi_{+1}(\tilde{p}_{a},\tilde{p}_{b};\tilde{p}_{1},\ldots,\tilde{p}_{ij},\ldots,\tilde{p}_{n};p_{i},p_{j})=\frac{{\rm d}\Phi_{n}(p_{a},p_{b};p_{1},\ldots,p_{i},\ldots,p_{j},\ldots,p_{n})}{{\rm d}\Phi_{n-1}(\tilde{p}_{a},\tilde{p}_{b};\tilde{p}_{1},\ldots,\tilde{p}_{ij},\ldots,\backslash\!\!\!\!\tilde{p}_{j},\ldots,\tilde{p}_{n})}\;.$$ (32)

This expression naturally generalizes to $D$ dimensions. It can be computed using the lowest final-state multiplicity, i.e.,  $n=3$. In order to do so, we first derive a suitable expression for the $D$-dimensional two-particle phase space in the frame of an arbitrary, time-like momentum $n$. It reads

$$\begin{split}{\rm d}\Phi_{2}(P;p,q)=&\;\frac{1}{(2\pi)^{2D-2}}\frac{{\rm d}^{D-1}p}{2p^{0}}\frac{{\rm d}^{D-1}q}{2q^{0}}\,(2\pi)^{D}\delta^{D}(P-p-q)=\frac{1}{4(2\pi)^{D-2}}\,\frac{|\vec{p}\,|^{D-2}}{P^{0}|\vec{p}\,|-p^{0}|\vec{P}|\cos{\theta}_{pn}}\,{\rm d}\Omega_{p,n}^{D-2}\;.\end{split}$$ (33)

where ${\rm d}\Omega_{p,n}$ is the integral over the $D-1$ dimensional solid angle of $p^{\mu}$ in the frame of $n^{\mu}$. All energies and momenta are computed in this frame, and we have defined $P=p+q$. We can re-write the momentum-dependent denominator factor on the right-hand side of Eq. (33) as

$$\frac{p^{0}}{|\vec{p}|}\left(P^{0}p^{0}-|\vec{p}||\vec{P}|\cos{\theta}_{pn}\right)-P^{0}|\vec{p}|\left(\frac{(p^{0})^{2}}{|\vec{p}|^{2}}-1\right)=\frac{(pn)(Pp)-p^{2}(Pn)}{\sqrt{(pn)^{2}-p^{2}n^{2}}}\;.$$ (34)

Inserting this into Eq. (33), and using $D=4-2\varepsilon$, we obtain a manifestly covariant form of the phase-space element,

$${\rm d}{\Phi}_{2}(P;p,q)=\frac{(4\pi^{2}n^{2})^{\varepsilon}}{16\pi^{2}}\,\frac{((pn)^{2}-p^{2}n^{2})^{3/2-\varepsilon}}{((pn)(Pp)-p^{2}(Pn))n^{2}}\,{\rm d}\Omega_{p,n}^{2-2\varepsilon}\;.$$ (35)

IV.1 Radiation kinematics

To derive the emission phase space for final-state splittings in the radiation kinematics of Sec. III.1, we make use of the standard factorization formula

$${\rm d}\Phi_{3}(-K;p_{i},p_{j},q)={\rm d}\Phi_{2}(-K;p_{j},-n)\,\frac{{\rm d}n^{2}}{2\pi}\,{\rm d}\Phi_{2}(-n;p_{i},q)\;,$$ (36)

where $q=\sum_{k\neq i,j}p_{k}$ is the sum of all final-state momenta except $p_{i}$ and $p_{j}$. We can use Eq. (35) to relate the phase-space factor ${\rm d}\Phi_{2}(-n;p_{i},q)$ to the underlying Born phase space as follows

$$\begin{split}\frac{{\rm d}\Phi_{2}(-n;{p}_{i},q)}{{\rm d}\Phi_{2}(-\tilde{K};\tilde{p}_{ij},\tilde{q})}=&\;\bigg{(}\frac{(p_{i}n)^{2}-p_{i}^{2}n^{2}}{(\tilde{p}_{ij}n)^{2}-\tilde{p}_{ij}^{2}n^{2}}\bigg{)}^{3/2-\varepsilon}\;\frac{(\tilde{p}_{ij}n)(\tilde{p}_{ij}\tilde{K})-\tilde{p}_{ij}^{2}(\tilde{K}n)}{({p}_{i}n)^{2}-{p}_{i}^{2}n^{2}}\\ =&\;\frac{\big{(}1-2\mu_{ij}^{2}(2(\kappa-\mu_{i}^{2})-z+\sigma_{i,ij}(1+2\mu_{ij}^{2}))\big{)}\big{[}zv_{p_{i},n}\big{]}^{1-2\varepsilon}}{\big{[}\big{(}1+2\mu_{ij}^{2}(1-\sigma_{i,ij})\big{)}^{2}-4\mu_{ij}^{2}(1-z+\kappa+\mu_{ij}^{2}-\mu_{i}^{2})\big{]}^{3/2-\varepsilon}}\;,\end{split}$$ (37)

where we have defined

$$\sigma_{i,ij}=\bar{z}+\frac{\mu_{i}^{2}/\mu_{ij}^{2}-\bar{z}^{2}}{\bar{z}v_{\tilde{p}_{ij},\tilde{K}}}\frac{1-2\bar{\kappa}\mu_{ij}^{2}}{2}\;.$$ (38)

The differential two-body phase-space element for the production of $n^{\mu}$ and $p_{j}^{\mu}$ is given by

$$\begin{split}{\rm d}{\Phi}_{2}(-K;n,p_{j})=&\;\frac{(4\pi^{2})^{\varepsilon}}{16\pi^{2}}\,\frac{\big{(}(p_{j}n)^{2}-p_{j}^{2}n^{2}\big{)}^{1/2-\varepsilon}}{(n^{2})^{1-\varepsilon}}\,{\rm d}\Omega_{j,n}^{2-2\varepsilon}\\ =&\;\frac{(2\tilde{p}_{ij}\tilde{K})^{-\varepsilon}}{(16\pi^{2})^{1-\varepsilon}}\,\frac{\big{[}(1-z+\mu_{ij}^{2}-\mu_{i}^{2}+\mu_{j}^{2})v_{p_{j},n}\big{]}^{1-2\varepsilon}}{2(1-z+\kappa+\mu_{ij}^{2}-\mu_{i}^{2})^{1-\varepsilon}}\,{\rm d}\Omega_{j,n}^{2-2\varepsilon}\;,\end{split}$$ (39)

Finally, we combine Eqs. (37) and (39) to obtain the single-emission phase space element in Eq. (32)

$$\begin{split}&{};\ldots,\tilde{p}_{ij},\ldots;p_{i},p_{j})\\ &\quad=\bigg{(}\frac{2\tilde{p}_{ij}\tilde{K}}{16\pi^{2}}\bigg{)}^{1-\varepsilon}\,\frac{\big{(}1-2\mu_{ij}^{2}(2(\kappa-\mu_{i}^{2})-z+\sigma_{i,ij}(1+2\mu_{ij}^{2}))\big{)}\big{[}v_{p_{i},n}v_{p_{j},n}\big{]}^{1-2\varepsilon}}{\big{[}\big{(}1+2\mu_{ij}^{2}(1-\sigma_{i,ij})\big{)}^{2}-4\mu_{ij}^{2}(1-z+\kappa+\mu_{ij}^{2}-\mu_{i}^{2})\big{]}^{3/2-\varepsilon}}\\ &\qquad\times\frac{\big{(}z(1-z+\mu_{ij}^{2}-\mu_{i}^{2}+\mu_{j}^{2})\big{)}^{1-2\varepsilon}}{(1-z+\kappa+\mu_{ij}^{2}-\mu_{i}^{2})^{1-\varepsilon}}\,{\rm d}z\,\frac{{\rm d}\Omega_{j,n}^{2-2\varepsilon}}{4\pi}\;.\end{split}$$ (40)

In the massless limit, this simplifies to

$${\rm d}\Phi_{+1}(\tilde{p}_{a},\tilde{p}_{b};\ldots,\tilde{p}_{ij},\ldots;p_{i},p_{j})=\bigg{(}\frac{2\tilde{p}_{ij}\tilde{K}}{16\pi^{2}}\bigg{)}^{1-\varepsilon}\,\frac{\big{(}z(1-z)\big{)}^{1-2\varepsilon}}{(1-z+\kappa)^{1-\varepsilon}}\,{\rm d}z\,\frac{{\rm d}\Omega_{j,n}^{2-2\varepsilon}}{4\pi}\;.$$ (41)

We demonstrate in App. A.1 that this expression agrees with the result derived in Ref. Catani and Seymour (1997).

IV.2 Splitting kinematics

To derive the emission phase space for final-state radiation in the splitting kinematics of Sec. III.2, the recoiler is chosen to be the sum of the remaining final-state partons. We make use of the standard factorization formula

$${\rm d}\Phi_{3}(q;p_{i},p_{j},K)={\rm d}\Phi_{2}(q;K,p_{ij})\,\frac{{\rm d}p_{ij}^{2}}{2\pi}\,{\rm d}\Phi_{2}(p_{ij};p_{i},p_{j})\;,$$ (42)

where $q=\sum_{k\neq i,j}p_{k}$ is the sum of final-state momenta except $p_{i}$ and $p_{j}$. Working in the frame of $q=\tilde{K}+\tilde{p}_{ij}$, we can use Eq. (35) to relate the phase-space factor ${\rm d}\Phi_{2}(q;K,p_{ij})$ to the underlying Born differential phase space element ${\rm d}\Phi_{2}(q;\tilde{K},\tilde{p}_{ij})$ as follows

$$\begin{split}\frac{{\rm d}\Phi_{2}(q;K,p_{ij})}{{\rm d}\Phi_{2}(q;\tilde{K},\tilde{p}_{ij})}=&\;\bigg{(}\frac{(Kp_{ij})^{2}-K^{2}p_{ij}^{2}}{(\tilde{K}\tilde{p}_{ij})^{2}-\tilde{K}^{2}\tilde{p}_{ij}^{2}}\bigg{)}^{1/2-\varepsilon}=(1-y)^{1-2\varepsilon}(1+\mu_{ij}^{2}-\mu_{i}^{2}-\mu_{j}^{2})^{1-2\varepsilon}\bigg{(}\frac{v_{p_{ij},K}}{v_{\tilde{p}_{ij},\tilde{K}}}\bigg{)}^{1-2\varepsilon}\;,\end{split}$$ (43)

The decay of the intermediate off-shell parton is associated with the differential phase-space element

$$\begin{split}{\rm d}{\Phi}_{2}(p_{ij};p_{i},p_{j})=&\;\frac{(4\pi^{2})^{\varepsilon}}{16\pi^{2}}\,\frac{((p_{i}p_{j})^{2}-p_{i}^{2}p_{j}^{2})^{1/2-\varepsilon}}{(p_{ij}^{2})^{1-\varepsilon}}\,{\rm d}\Omega_{i,ij}^{2-2\varepsilon}\\ =&\;\frac{(2\tilde{p}_{ij}\tilde{K})^{-\varepsilon}}{(16\pi^{2})^{1-\varepsilon}}\,\frac{\big{(}y^{2}(1+\mu_{ij}^{2}-\mu_{i}^{2}-\mu_{j}^{2})^{2}-4\mu_{i}^{2}\mu_{j}^{2}\big{)}^{1/2-\varepsilon}}{2\big{(}y(1+\mu_{ij}^{2}-\mu_{i}^{2}-\mu_{j}^{2})+\mu_{i}^{2}+\mu_{j}^{2}\big{)}^{1-\varepsilon}}\,{\rm d}\Omega_{i,ij}^{2-2\varepsilon}\;.\end{split}$$ (44)

Finally, we combine Eqs. (43) and (44) to obtain

$$\begin{split}&{};\tilde{p}_{1},\ldots,\tilde{p}_{ij},\ldots,\tilde{p}_{n};p_{i},p_{j})\\ &\qquad=\bigg{(}\frac{2\tilde{p}_{ij}\tilde{K}}{16\pi^{2}}\bigg{)}^{1-\varepsilon}\,(1-y)^{1-2\varepsilon}(1+\mu_{ij}^{2}-\mu_{i}^{2}-\mu_{j}^{2})^{2-2\varepsilon}\bigg{(}\frac{v_{p_{ij},K}}{v_{\tilde{p}_{ij},\tilde{K}}}\bigg{)}^{1-2\varepsilon}\\ &\qquad\qquad\times\frac{\big{(}y^{2}(1+\mu_{ij}^{2}-\mu_{i}^{2}-\mu_{j}^{2})^{2}-4\mu_{i}^{2}\mu_{j}^{2}\big{)}^{1/2-\varepsilon}}{\big{(}y(1+\mu_{ij}^{2}-\mu_{i}^{2}-\mu_{j}^{2})+\mu_{i}^{2}+\mu_{j}^{2}\big{)}^{1-\varepsilon}}\,{\rm d}y\,\frac{{\rm d}\Omega_{i,ij}^{2-2\varepsilon}}{4\pi}\;.\end{split}$$ (45)

In the massless limit, this simplifies to

$$\begin{split}{\rm d}\Phi_{+1}(\tilde{p}_{a},\tilde{p}_{b};\tilde{p}_{1},\ldots,\tilde{p}_{ij},\ldots,\tilde{p}_{n};p_{i},p_{j})=\bigg{(}\frac{2\tilde{p}_{ij}\tilde{K}}{16\pi^{2}}\bigg{)}^{1-\varepsilon}\,(1-y)^{1-2\varepsilon}\,y^{-\varepsilon}\,{\rm d}y\,\frac{{\rm d}\Omega_{i,ij}^{2-2\varepsilon}}{4\pi}\;.\end{split}$$ (46)

In App. A.2, we will show the equivalence of Eq. (45) to the single-emission differential phase-space element of Catani et al. (2002).

V Infrared subtraction terms at next-to-leading order

The matching of parton showers to fixed-order NLO calculations in dimensional regularization based on the MC@NLO algorithm Frixione and Webber (2002) requires the knowledge of integrated splitting functions in $D=4-2\varepsilon$ dimensions. Since our technique for massive parton evolution is modeled on the Catani-Seymour identified particle subtraction and the Alaric parton-shower, we can use the methods developed in Höche et al. (2019); Liebschner (2020). We will limit the discussion to the main changes needed to implement the algorithm for massive partons. The results of this section provide an extension of the subtraction method for identified hadrons first introduced in Catani and Seymour (1997).

V.1 Soft angular integrals

In this section, we compute the angular integrals of the partial fractioned soft eikonal. While we focus on pure final-state radiation, the results of this integration are generic and apply to the case of final-state and initial-state emission of massless vector bosons. The integrand is given by Eq. (8),

$$\begin{split}\bar{W}_{ik,j}^{i}=&\;\frac{l_{ik}^{2}}{2(l_{i}l_{j})(l_{ik}l_{j})}-\frac{l_{i}^{2}}{2(l_{i}l_{j})^{2}}-\frac{l_{k}^{2}}{2(l_{i}l_{j})(l_{k}l_{j})}\\ \end{split}$$ (47)

We can write the angular phase space for the emission of the massless particle $j$ as

$$\int{\rm d}\Omega_{j}^{2-2\varepsilon}=\Omega(1-2\varepsilon)\int{\rm d}\cos\theta(\sin^{2}\theta)^{-\varepsilon}\int{\rm d}\phi(\sin^{2}\phi)^{-\varepsilon}\;,$$ (48)

where $\Omega(n)=2\pi^{n/2}/\Gamma(n/2)$ is the $d$-dimensional line element, in particular $\Omega(1-2\varepsilon)=2(4\pi)^{-\varepsilon}\Gamma(1-\varepsilon)/\Gamma(1-2\varepsilon)$. We finally find the following expression for the cases with massive emitter

$$\begin{split}&\int\frac{{\rm d}\Omega_{j}^{2-2\varepsilon}}{\Omega(1-2\varepsilon)}\;\bar{W}_{ik,j}^{i}=\frac{l_{ik}^{2}}{4}I_{1,1}^{(2)}\left(\frac{l_{i}l_{ik}}{2},l_{i}^{2},\frac{l_{ik}^{2}}{4}\right)-\frac{l_{i}^{2}}{2}I_{2}^{(1)}(l_{i}^{2})-\frac{l_{k}^{2}}{2}I_{1,1}^{(2)}\left(l_{i}l_{k},l_{i}^{2},l_{k}^{2}\right)\;.\end{split}$$ (49)

For massless emitter, we obtain (see also Catani and Seymour (1997); Herren et al. (2022))

$$\begin{split}&\int\frac{{\rm d}\Omega_{j}^{2-2\varepsilon}}{\Omega(1-2\varepsilon)}\;\bar{W}_{ik,j}^{i}=\frac{l_{ik}^{2}}{4}I_{1,1}^{(1)}\left(\frac{l_{i}l_{ik}}{2},\frac{l_{ik}^{2}}{4}\right)-\frac{l_{k}^{2}}{2}I_{1,1}^{(1)}(l_{i}l_{k},l_{k}^{2})\;,\end{split}$$ (50)

The angular integrals $I_{1,1}$ and $I_{2}$ have been computed in van Neerven (1986); Beenakker et al. (1989); Somogyi (2011); Isidori et al. (2020); Lyubovitskij et al. (2021). They read Lyubovitskij et al. (2021)

$$\begin{split}I_{1,1}^{(2)}(v_{12},v_{11},v_{22})={}&\frac{\pi}{\sqrt{v_{12}^{2}-v_{11}v_{22}}}\left\{\log{\frac{v_{12}+\sqrt{v_{12}^{2}-v_{11}v_{22}}}{v_{12}-\sqrt{v_{12}^{2}-v_{11}v_{22}}}}+\varepsilon\left(\frac{1}{2}\log^{2}{\frac{v_{11}}{v_{13}^{2}}}-\frac{1}{2}\log^{2}\frac{v_{22}}{v_{23}^{2}}\right.\right.\\ &\left.\left.+2{\rm Li}_{2}\left(1-\frac{v_{13}}{1-\sqrt{1-v_{11}}}\right)+2{\rm Li}_{2}\left(1-\frac{v_{13}}{1+\sqrt{1-v_{11}}}\right)\right.\right.\\ &\left.\left.-2{\rm Li}_{2}\left(1-\frac{v_{23}}{1-\sqrt{1-v_{22}}}\right)-2{\rm Li}_{2}\left(1-\frac{v_{23}}{1+\sqrt{1-v_{22}}}\right)\right)+\mathcal{O}(\varepsilon^{2})\right\}\;,\\ I_{1,1}^{(1)}(v_{12},v_{11})={}&\frac{\pi}{v_{12}}\left\{-\frac{1}{\varepsilon}-\log\frac{v_{11}}{v_{12}^{2}}-\varepsilon\left(\frac{1}{2}\log^{2}\frac{v_{11}}{v_{12}^{2}}\right.\right.\\ &\left.\left.+2{\rm Li}_{2}\left(1-\frac{v_{12}}{1-\sqrt{1-v_{11}}}\right)+2{\rm Li}_{2}\left(1-\frac{v_{12}}{1+\sqrt{1-v_{11}}}\right)\right)+\mathcal{O}(\varepsilon^{2})\right\}\;,\\ I_{2}^{(1)}(v)={}&\frac{2\pi}{v}\left\{1+\frac{\varepsilon}{\sqrt{1-v}}\log\frac{1+\sqrt{1-v}}{1-\sqrt{1-v}}+\mathcal{O}(\varepsilon^{2})\right\}\;.\end{split}$$ (51)

where the velocities are defined in the notation of Ref. Lyubovitskij et al. (2021), and where

$$\begin{split}v_{13}=&\;(1-\lambda)v_{11}+\lambda v_{12}\;,\qquad v_{23}=(1-\lambda)v_{12}+\lambda v_{22}\;,\qquad\lambda=\frac{v_{11}-v_{12}-\sqrt{v_{12}^{2}-v_{11}v_{22}}}{v_{11}+v_{22}-2v_{12}}\;.\end{split}$$ (52)

Note that spurious collinear poles appear on the right-hand side of Eq. (50), which cancel between $I_{1,1}^{(1)}(l_{i}l_{ik}/2,l_{ik}^{2}/4)/2$ and $I_{1,1}^{(1)}(l_{i}l_{k},l_{k}^{2})$. In addition, in the fully massless case $I_{1,1}^{(1)}$ is related to the massless two-denominator integral and gives the simple result Herren et al. (2022)

$$\begin{split}I_{1,1}^{(1)}(v_{12},v_{12})=&\;-\frac{\pi}{v_{12}^{1+\varepsilon}}\left\{\frac{1}{\varepsilon}+\varepsilon\,{\rm Li}_{2}\left(1-v_{12}\right)+\mathcal{O}(\varepsilon^{2})\right\}\;.\end{split}$$ (53)

V.2 Soft energy integrals

In this section, we introduce the basic techniques to solve the energy integrals in Eq. (40). After performing the angular integrals, we are left with the additional $z$-dependence induced by the energy denominator in Eq. (4). We focus on the cases relevant for QCD and QED soft radiation, where $\mu_{ij}=\mu_{i}$ and $\mu_{j}=0$. The differential emission probability per dipole is then given by

$$\begin{split}&{};\ldots,\tilde{p}_{ij},\ldots;p_{i},p_{j})\,\frac{n^{2}}{(p_{j}n)^{2}}\,\bar{W}_{ik,j}^{i}\\ &\quad=\frac{(4\pi)^{\varepsilon}}{16\pi^{2}}\,\frac{\Gamma(1-\varepsilon)}{\Gamma(1-2\varepsilon)}\,(2\tilde{p}_{ij}\tilde{K})^{-\varepsilon}\,\frac{\big{(}1-2\mu_{i}^{2}(2(\kappa-\mu_{i}^{2})-z+\sigma_{i,ij}(1+2\mu_{i}^{2}))\big{)}v_{p_{i},n}^{1-2\varepsilon}}{\big{[}\big{(}1+2\mu_{i}^{2}(1-\sigma_{i,ij})\big{)}^{2}-4\mu_{i}^{2}(1-z+\kappa)\big{]}^{3/2-\varepsilon}}\\ &\qquad\times\,\frac{z^{1-2\varepsilon}}{(1-z+\kappa)^{-\varepsilon}}\,\frac{{\rm d}z}{(1-z)^{1+2\varepsilon}}\,\frac{2}{\pi}\frac{{\rm d}\Omega_{j,n}^{2-2\varepsilon}}{\Omega(1-2\varepsilon)}\bar{W}_{ik,j}^{i}\;.\end{split}$$ (54)

In order to carry out the integration over the energy fraction, $z$, we expand the integrand into a Laurent series. The differential soft subtraction counterterm, summed over all dipoles, is given by

$$\begin{split}{\rm d}\sigma^{S}_{S}=&\;-8\pi\alpha_{s}\mu^{2\varepsilon}\sum_{\widetilde{\imath\jmath},k}{\rm d}\Phi_{+1}(\tilde{p}_{a},\tilde{p}_{b};\ldots,\tilde{p}_{ij},\ldots;p_{i},p_{j})\,{\bf T}_{\widetilde{\imath\jmath}}{\bf T}_{k}\frac{n^{2}}{(p_{j}n)^{2}}\,\bar{W}_{ik,j}^{i}\\ =&\;-\frac{\alpha_{s}}{2\pi}\frac{1}{\Gamma(1-\varepsilon)}\sum_{\widetilde{\imath\jmath},k}\frac{{\bf T}_{\widetilde{\imath\jmath}}{\bf T}_{k}}{{\bf T}_{\widetilde{\imath\jmath}}^{2}}\bigg{(}\frac{4\pi\mu^{2}}{2p_{i}p_{k}}\bigg{)}^{\varepsilon}\left\{-\frac{\delta(1-z)}{\varepsilon}(1-z_{-})^{-\varepsilon}+\frac{2z}{\left[1-z\right]_{+}}-4\varepsilon\,z\left[\frac{\log(1-z)}{1-z}\right]_{+}\right\}\,{\rm d}z\\ &\;\qquad\times\,\mathcal{J}_{ik}(z,\mu_{i}^{2},\kappa)\,\bigg{(}\frac{p_{k}n}{p_{i}n}\bigg{)}^{\varepsilon}\bigg{(}\frac{l_{ik}^{2}-l_{i}^{2}-l_{k}^{2}}{4}\bigg{)}^{\varepsilon}\,\frac{C_{\widetilde{\imath\jmath}}}{\pi}\frac{\Gamma(1-\varepsilon)^{2}}{\Gamma(1-2\varepsilon)}\,\bigg{[}\frac{{\rm d}\Omega_{j,n}^{2-2\varepsilon}}{\Omega(1-2\varepsilon)}\,\bar{W}_{ik,j}^{i}\bigg{]}\;.\end{split}$$ (55)

where $z_{-}=2\mu_{i}^{2}\big{(}\sqrt{1+(1+\kappa)/\mu_{i}^{2}}-1\big{)}$, and where the mass correction factor reads

$$\mathcal{J}_{ik}(z,\mu_{i}^{2},\kappa)=\frac{\big{(}1-2\mu_{i}^{2}(2(\kappa-\mu_{i}^{2})-z+\sigma_{i,ij}(1+2\mu_{i}^{2}))\big{)}v_{p_{i},n}^{1-2\varepsilon}}{\big{[}\big{(}1+2\mu_{i}^{2}(1-\sigma_{i,ij})\big{)}^{2}-4\mu_{i}^{2}(1-z+\kappa)\big{]}^{3/2-\varepsilon}}\,$$ (56)

and where the massless limit is given by $\mathcal{J}_{ik}(z,0,\kappa)=1$. All $1/\varepsilon$ poles have been extracted, and we can (after expanding the $\varepsilon$-dependent prefactors) compute the final result as an integral over the delta functions and plus distributions in $z$. In general, some of the terms must be computed numerically, as $n$ implicitly depends on $z$, see Eq. (14). In order to apply Eq. (56) to processes with resolved hadrons, we can make use of the formalism derived in Catani and Seymour (1997); Höche et al. (2019); Liebschner (2020). The $1/\varepsilon$ pole proportional to $2z/[1-z]_{+}$ can then be combined with the soft enhanced part of the collinear mass factorization counterterms. The extension to initial-state radiation requires a repetition of the derivation in Sec. IV.1 for initial-state kinematics. We will discuss this in a future publication.

V.3 Collinear integrals

To compute the purely collinear counterterms, we use the splitting kinematics and make use of the results in App. A.2. Note that we also use the corresponding definitions of the scaled masses. We define the purely collinear anomalous dimension in terms of the collinear remainders in Eq. (12)

$$\bar{\gamma}_{ab}(\varepsilon)=\frac{\Omega(2-2\varepsilon)}{\Omega(3-2\varepsilon)}\Big{(}\frac{z_{+}-z_{-}}{2}\Big{)}^{-1+2\varepsilon}\int_{z_{-}}^{z_{+}}{\rm d}z\,\big{(}(z-z_{-})(z_{+}-z)\big{)}^{-\varepsilon}\,C_{ab}(z,\varepsilon)\;,$$ (57)

where $z_{\pm}$ are given by Eq. (77). This leads to

$$\begin{split}\bar{\gamma}_{qq}(\varepsilon)=&\;C_{F}(1-\varepsilon)\left(1-\frac{z_{+}+z_{-}}{2}\right)\;,\\ \bar{\gamma}_{gg}(\varepsilon)=&\;C_{A}\left(\frac{z_{+}+z_{-}}{2}-\frac{(1-\varepsilon/2)(z_{+}+z_{-})^{2}-z_{+}z_{-}}{3-2\varepsilon}\right)\;,\\ \bar{\gamma}_{gq}(\varepsilon)=&\;T_{R}\left(1+\frac{z_{+}^{2}+z_{-}^{2}-(z_{+}+z_{-})}{1-\varepsilon}-\frac{(z_{+}-z_{-})^{2}}{3-2\varepsilon}\right)\;.\end{split}$$ (58)

The complete counterterm can now be obtained by Laurent series expansion, using Eq. (82)

$$\begin{split}&\int{\rm d}\Phi_{+1}(q;\tilde{p}_{ij},\tilde{K};p_{i},p_{j})\frac{8\pi\alpha_{s}\mu^{2\varepsilon}C_{ab}}{(p_{i}+p_{j})^{2}-m_{ij}^{2}}\\ &\qquad=\;\frac{\alpha_{s}}{2\pi}\,\bigg{(}\frac{\mu^{2}}{Q^{2}}\bigg{)}^{\varepsilon}\frac{\Omega(3-2\varepsilon)}{(4\pi)^{1-2\varepsilon}}\,(1-\hat{\mu}_{i}^{2}-\hat{\mu}_{j}^{2}-\hat{\kappa})^{2-2\varepsilon}\lambda(1,\hat{\mu}_{ij}^{2},\hat{\kappa})^{-1/2+\varepsilon}\\ &\qquad\qquad\times\int{\rm d}y\,(1-y)^{1-2\varepsilon}\,\frac{\big{[}\,y(1-\hat{\mu}_{i}^{2}-\hat{\mu}_{j}^{2}-\hat{\kappa})+\hat{\mu}_{i}^{2}+\hat{\mu}_{j}^{2}\,\big{]}^{-\varepsilon}}{y(1-\hat{\mu}_{i}^{2}-\hat{\mu}_{j}^{2}-\hat{\kappa})+\hat{\mu}_{i}^{2}+\hat{\mu}_{j}^{2}-\hat{\mu}_{ij}^{2}}\,\big{(}z_{+}-z_{-}\big{)}^{1-2\varepsilon}\,\bar{\gamma}_{ab}(\varepsilon)\;.\end{split}$$ (59)

There are three variants of Eq. (59) relevant for QCD. The first describes the collinear splitting of a massive quark into a quark and a gluon, and is characterized by $\hat{\mu}_{i}=\hat{\mu}_{ij}>0$ and $\hat{\mu}_{j}=0$. We obtain

$$\begin{split}&\int{\rm d}\Phi_{+1}(q;\tilde{p}_{ij},\tilde{K};p_{i},p_{j})\frac{8\pi\alpha_{s}\mu^{2\varepsilon}C_{qq}}{(p_{i}+p_{j})^{2}-m_{ij}^{2}}=\frac{\alpha_{s}}{2\pi}\,\frac{(1-\hat{\mu}_{i}^{2}-\hat{\kappa})}{\sqrt{\lambda(1,\hat{\mu}_{i}^{2},\hat{\kappa})}}\\ &\qquad\times\int\frac{{\rm d}y}{y(1-\hat{\mu}_{i}^{2}-\hat{\kappa})+\hat{\mu}_{i}^{2}}\,\sqrt{\big{[}(1-y)(1-\hat{\mu}_{i}^{2}-\hat{\kappa})+2\hat{\kappa}\big{]}^{2}-4\hat{\kappa}}\;\bar{\gamma}_{qq}(0)\;,\end{split}$$ (60)

Note that the result is infrared finite. This is expected because the only poles that occur in the radiation off massive quarks have their origin in soft gluon emission, which is fully accounted for by the eikonal integral in Eq. (54). This also shows that in our subtraction scheme, there is a hierarchy between the soft and the collinear enhancements, as required for a proper classification of leading and sub-leading logarithms.

The second case is the splitting of a gluon into two massive quarks, which is characterized by $\hat{\mu}_{i}=\hat{\mu}_{j}>0$ and $\hat{\mu}_{ij}=0$. The result is finite and reads

$$\begin{split}&\int{\rm d}\Phi_{+1}(q;\tilde{p}_{ij},\tilde{K};p_{i},p_{j})\frac{8\pi\alpha_{s}\mu^{2\varepsilon}C_{gq}}{(p_{i}+p_{j})^{2}}=\;\frac{\alpha_{s}}{2\pi}\,\bigg{(}1-\frac{2\hat{\mu}_{i}^{2}}{1-\hat{\kappa}}\bigg{)}\\ &\qquad\times\int{\rm d}y\,\frac{\sqrt{y^{2}(1-2\hat{\mu}_{i}^{2}-\hat{\kappa})^{2}-4\hat{\mu}_{i}^{4}}}{(y(1-2\hat{\mu}_{i}^{2}-\hat{\kappa})+2\hat{\mu}_{i}^{2})^{2}}\sqrt{\big{[}(1-y)(1-2\hat{\mu}_{i}^{2}-\hat{\kappa})+2\hat{\kappa}\big{]}^{2}-4\hat{\kappa}}\;\bar{\gamma}_{gq}(0)\;.\end{split}$$ (61)

The final case describes the collinear decay of a massless parton into two massless partons, and is characterized by $\hat{\mu}_{i}=\hat{\mu}_{j}=\hat{\mu}_{ij}=0$. This term is collinearly divergent, and we obtain

$$\begin{split}&\int{\rm d}\Phi_{+1}(q;\tilde{p}_{ij},\tilde{K};p_{i},p_{j})\frac{8\pi\alpha_{s}\mu^{2\varepsilon}C_{ab}}{(p_{i}+p_{j})^{2}}=\frac{\alpha_{s}}{2\pi}\bigg{(}\frac{\mu^{2}}{Q^{2}}\bigg{)}^{\varepsilon}\frac{(4\pi)^{\varepsilon}}{\Gamma(1-\varepsilon)}\\ &\qquad\times\int{\rm d}y\,\left\{-\frac{\delta(y)}{\varepsilon}\bigg{(}\frac{1-\sqrt{\hat{\kappa}}}{1+\sqrt{\hat{\kappa}}}\bigg{)}^{-\varepsilon}+\frac{1}{[y]_{+}}\right\}\frac{\big{(}\big{[}(1-y)(1-\hat{\kappa})+2\hat{\kappa}\big{]}^{2}-4\hat{\kappa}\big{)}^{1/2-\varepsilon}}{(1-\hat{\kappa})^{1-\varepsilon}}\,\frac{\Gamma(1-\varepsilon)^{2}}{\Gamma(2-2\varepsilon)}\,\bar{\gamma}_{ab}(\varepsilon)\;.\end{split}$$ (62)

The treatment of initial-state singularities will be discussed in a future publication.

VI Kinematical effects on sub-leading logarithmic corrections

We now demonstrate that our kinematics mapping satisfies the criteria for NLL accuracy laid out in  Dasgupta et al. (2018, 2020), extending the discussion of the massless case in Herren et al. (2022). We refer the reader to Herren et al. (2022) for more details on the general method of the proof. The analysis is based on the technique for general final-state resummation introduced in Banfi et al. (2005). Following the notation of Sec. 2 of Banfi et al. (2005), the observable to be resummed is denoted by $v$, while the hard partons and soft emission have momenta $p_{1},\ldots,p_{n}$ and $k$, respectively. The observable is a function of these momenta, $v=V(\{p\},\{k\})$, and for any radiating color dipole formed by the hard momenta, $p_{i}$ and $p_{j}$, the momentum of a single emission can be parametrized as

$$k=z_{i,j}p_{i}+z_{j,i}p_{j}+k_{T,ij}\;,\qquad\text{where}\qquad k_{T,ij}^{2}=2p_{i}p_{j}\,z_{i,j}\,z_{j,i}\;.$$ (63)

with rapidity $\eta_{ij}=1/2\ln(z_{i,j}/z_{j,i})$. An observable can be expressed as

$$V(k)=d_{l}\left(\frac{k_{T,l}}{Q}\right)^{a}e^{-b_{l}\eta_{l}}g_{l}(\phi_{l})\;,$$ (64)

with $k_{T,l}=k_{T,lj}$ and $\eta_{l}=\eta_{lj}$ for $j\nparallel l$ in the collinear limit. A natural extension of Eq. (63) to the case of massive emitters would be

$$k=(z_{i,j}-\bar{\mu}_{j}^{2}z_{j,i})\,p_{i}+(z_{j,i}-\bar{\mu}_{i}^{2}z_{i,j})\,p_{j}+k_{T,ij}\;,\qquad\text{where}\qquad k_{T,ij}^{2}=2p_{i}p_{j}\frac{2v_{p_{i},p_{j}}}{1+v_{p_{i},p_{j}}}\,z_{i,j}\,z_{j,i}\;.$$ (65)

In the quasi-collinear limit, this Sudakov decomposition agrees with the one given in Catani et al. (2001); Cacciari and Catani (2001) if the auxiliary (light-like) vector defining the anti-collinear direction is chosen as the direction of the color partner of the QCD dipole. In particular, for constant $z_{i,j}$, $z_{j,i}$ and small $k_{T}$, the gluon momentum behaves as

$$k\overset{k_{T}^{2}\ll p_{i}p_{j}}{\underset{m_{i/j}^{2}\varpropto k_{T}^{2}}{\longrightarrow}}z_{i,j}p_{i}+z_{j,i}p_{j}+k_{T,ij}+\mathcal{O}(k_{T}^{2})\;,$$ (66)

in complete agreement with Eq. (63). In the quasi-collinear limit, the value of the observable $V(k)$ will therefore be unchanged from the case of massless evolution. However, both the Sudakov radiator and the $\mathcal{F}$ function will change, due to the modified splitting functions, Eq. (10), and the effects of masses on the integration boundaries.

The Lund plane for gluon emission off a dipole containing a massive quark with mass $m$ and energy $E$ will have a smooth upper rapidity bound at $\eta=\ln(E/m)$, consistent with the well known dead-cone effect Marchesini and Webber (1990); Dokshitzer et al. (1991); Ghira et al. (2023). When the similarity transformation introduced in Sec. 2.2.3 of Banfi et al. (2005) is generalized to the quasi-collinear limit, and applied to a process with massive quarks, the location of this boundary in relative rapidity is unchanged, because the quasi-collinear limit requires $m\varpropto k_{T}$. The sub-leading logarithmic terms in the resummed cross section can then be extracted similarly to the massless case, but a change to the CMW scheme is required. We will not discuss the complete structure of the result, but address only the effects related to momentum mapping, which have been found to spoil NLL precision Dasgupta et al. (2018, 2020).

In order to prove that the momentum mapping of Sec. III.1 satisfies the criteria laid out in Dasgupta et al. (2018, 2020), we need to show that it has the same topological features as in the massless case Herren et al. (2022). This amounts to showing that hard, (quasi-)collinear emissions, i.e. highly energetic emissions in the direction of the hard partons, do not generate Lorentz transformations that change existing momenta by more than $\mathcal{O}((k_{T}^{2}/K^{2})^{\tilde{\rho}})$, where $\tilde{\rho}$ is a positive exponent. To show this, we analyze the behavior of Eq. (23) in the (quasi-)collinear region. We can split $K^{\mu}$ into its components along the original recoil momentum, $\tilde{K}^{\mu}$, the emitter, $\tilde{p}_{i}^{\mu}$, and the emission, ${p}_{j}^{\mu}$.

$$K^{\mu}=\tilde{K}^{\mu}-X^{\mu},\qquad{\rm where}\qquad X^{\mu}=p^{\mu}_{j}-(\tilde{p}_{ij}-p_{i})\;.$$ (67)

For the effect of the Lorentz transformation that determines the event topology after radiation, only the parametric form of $X^{\mu}$ is of interest. Using Eqs. (15) and (18), the small momentum $X^{\mu}$ for gluon radiation can be written as

$$\begin{split}X^{\mu}=&\;\bigg{(}\frac{1-z}{v_{\tilde{p}_{ij}\tilde{K}}\zeta}-(1-\bar{z})\bigg{)}\,\tilde{p}_{ij}^{\mu}-\frac{\mu_{ij}^{2}}{v_{\tilde{p}_{ij}\tilde{K}}}\bigg{(}\frac{1-z}{1+v_{\tilde{p}_{ij}\tilde{K}}}\,\tilde{K}^{\mu}+\frac{1-\bar{z}^{2}}{\bar{z}}\left(\tilde{K}^{\mu}-\bar{\kappa}\,\tilde{p}_{ij}^{\mu}\right)\bigg{)}\\ &\;+\bar{v}\,\frac{1+v_{\tilde{p}_{ij}\tilde{K}}}{2v_{\tilde{p}_{ij}\tilde{K}}}\left[\left(\tilde{K}^{\mu}-\bar{\kappa}\tilde{p}_{ij}^{\mu}\right)-\frac{1-\bar{z}+\bar{\kappa}}{\zeta}\left(\tilde{p}_{ij}^{\mu}-\bar{\mu}_{ij}^{2}\tilde{K}^{\mu}\right)\right]+k_{\perp}^{\mu}\;.\end{split}$$ (68)

In the quasi-collinear limit, this scales as $\mathcal{O}(k_{\perp})$, which is sufficient to make the effect of the Lorentz transformation vanish even for hard quasi-collinear splittings.

Finally, we note that the precise treatment of transverse recoil in the kinematics mapping of Sec. III.2 does not affect the resummed result at NLL precision, because the mapping is applied solely to the purely collinear parts of the splitting functions, Eq. (12). This can be understood with the help of Eq. (2.46) in Ref. Banfi et al. (2005), which will be structurally similar in the case of massive partons. The sub-leading logarithmic terms, which lead to a violation of NLL accuracy in many dipole showers, arise from accidental correlations between multiple soft-enhanced emissions. This means in particular, that the parton-shower equivalent of the integrals in Eq. (2.46) does not factorize in the strongly ordered limit. The differential probability for the emissions is given by the derivative of the radiator function in Eq. (2.21) of Ref. Banfi et al. (2005) and does not receive a contribution from the purely collinear remainders in Eq. (12), such that a change in the kinematics mapping cannot generate this particular type of NLL violation.

VII Comparison to experimental data

Figure 3: Alaric and Dire predictions in comparison to LEP data from Pfeifenschneider et al. (2000).

Figure 4: Alaric and Dire predictions in comparison to LEP data from Heister et al. (2004).

Figure 5: Alaric and Dire predictions in comparison to LEP data from Heister et al. (2001); Abdallah et al. (2011); Abbiendi et al. (2003); Abe et al. (2002).

In this section we present first numerical results obtained with the extension of the Alaric final-state parton shower to massive parton evolution. The implementation is part of the event generation framework Sherpa Gleisberg et al. (2004, 2009); Bothmann et al. (2019). We do not perform NLO matching or multi-jet merging, and we set $C_{F}=(N_{c}^{2}-1)/(2N_{c})=4/3$ and $C_{A}=3$. The quark masses are set to $m_{c}=1.42$ GeV and $m_{b}=4.75$ GeV and the same flavor thresholds are used for the evaluation of the strong coupling. The running coupling is evaluated at two loop accuracy, and we set $\alpha_{s}(m_{z})=0.118$. Following standard practice to improve the logarithmic accuracy of the parton shower, we employ the CMW scheme Catani et al. (1991). In this scheme, the soft eikonal contribution to the flavor conserving splitting functions is rescaled by $1+\alpha_{s}(t)/(2\pi)K$, where $K=(67/18-\pi^{2}/6)\,C_{A}-10/9\,T_{R}\,n_{f}$, and where $n_{f}$ is scale dependent with the same flavor thresholds as listed above. Velocity dependent corrections to $K$ should be included in principle, but the massless result provides an acceptable approximation at very large velocities Korchemsky and Radyushkin (1987); Kidonakis (2009), and can therefore be used for $b$-quark production at LEP, where $v\approx 0.99$. Our results include the simulation of hadronization using the Lund string fragmentation implemented in Pythia 6.4 Sjöstrand et al. (2006). We use the default hadronization parameters, apart from the following: PARJ(21)=0.3, PARJ(41)=0.4, PARJ(42)=0.45, PARJ(46)=0.5 and compare our predictions with those from the Dire parton shower Höche and Prestel (2015). All analyses are performed with Rivet Buckley et al. (2013).

Figure 3 displays predictions from the Alaric parton shower for differential jet rates in the Durham scheme compared to experimental results from the JADE and OPAL collaborations Pfeifenschneider et al. (2000). The $b$-quark mass corresponds to $y\sim 2.8\cdot 10^{-3}$, and hadronization effects dominate the predictions below $\sim 10^{-4}$. We observe good agreement with the experimental data. Overall, the quality of the results is similar to Ref. Herren et al. (2022), where heavy quark effects were modeled by thresholds. Figure 4 shows a comparison for event shapes measured by the ALEPH collaboration Heister et al. (2004). It can be expected that the numerical predictions will improve upon including matrix-element corrections, or when merging the parton shower with higher-multiplicity calculations. Again, the overall quality of the prediction is similar to Ref. Herren et al. (2022).

Finally, we show predictions for the $b$-quark fragmentation function as measured by the ALEPH Heister et al. (2001), DELPHI Abdallah et al. (2011), OPAL Abbiendi et al. (2003) and SLD Abe et al. (2002) collaborations. Figure 5 shows a fair agreement of both the Alaric and Dire predictions with experimental data. We note that both parton shower implementations use the same hadronization tune.

VIII Conclusions

We have introduced an extension of the recently proposed Alaric parton-shower model to the case of massive QCD evolution. An essential aspect of the new algorithm is the form of the collinearly matched massive eikonal, which is obtained by partial fractioning the angular component of the eikonal of a complete dipole. This technique preserves the positivity of the splitting function, thus leading to an excellent efficiency of the Monte-Carlo simulation. Inspired by the symmetry of the partonic final state in purely collinear splittings, we also introduced a dedicated kinematics mapping for this scenario and showed that it preserves the NLL precision of the overall simulation. We computed the infrared counterterms needed for the matching to fixed-order calculations at NLO accuracy, and discussed the logarithmic structure of the resummation in the case of heavy-quark evolution.

Several improvements of this algorithm are needed before it can be considered on par with the parton shower simulations used by past and current experiments. Clearly, spin correlations and dominant sub-leading color effects should be included. This can be achieved with the help of the techniques from Collins (1988); Knowles (1988a, b, 1990); Gustafson (1993); Dulat et al. (2018); Hamilton et al. (2021). An extension to initial-state evolution is needed for LHC phenomenology. It will need to account for the non-cancellation of certain types of singularities in processes with two massive initial states Caola et al. (2021). Finally, the algorithm should be extended to higher orders based on the techniques developed in Höche and Prestel (2017); Dulat et al. (2018); Gellersen et al. (2022).

In this context, we note that the all-orders (in $\varepsilon$) expressions from Sec. V, in conjunction with higher-order expressions for the angular integrals in Sec. V.1 that can be obtained from Lyubovitskij et al. (2021), can be used to compute the factorizable integrals at NNLO, thus providing a significant part of the components needed for an MC@NNLO matching Campbell et al. (2023). The computation of the remaining non-factorizable integrals is a further development needed in order to reach the precision targets of the high-luminosity LHC.

Appendix A Phase-space factorization in comparison to other infrared subtraction schemes

In this appendix, we compare the phase-space factorization in our method to the existing dipole subtraction schemes of Refs. Catani and Seymour (1997) and Catani et al. (2002). We will show that our generic form of the factorized phase space, derived from Eq. (35) can be used to obtain the relevant formulae, for pure final-state evolution.

A.1 Massless radiation kinematics

First, we show that Eq. (5.189) of Ref. Catani and Seymour (1997) can be derived from our generic expression, Eq. (35), using radiation kinematics. We start with the massless limit of Eq. (40)

$${\rm d}\Phi_{+1}(\tilde{p}_{a},\tilde{p}_{b};\ldots,\tilde{p}_{ij},\ldots;p_{i},p_{j})=\bigg{(}\frac{2\tilde{p}_{ij}\tilde{K}}{16\pi^{2}}\bigg{)}^{1-\varepsilon}\,\frac{\big{(}z(1-z)\big{)}^{1-2\varepsilon}}{(1-z+\kappa)^{1-\varepsilon}}\,{\rm d}z\,\frac{{\rm d}\Omega_{j,n}^{1-2\varepsilon}}{4\pi}\;.$$ (69)

Expressing the polar angle in the frame of $n$ in terms of $v$ and $z$ (see also Eq. (32) of Herren et al. (2022))

$$\cos\theta_{j,n}=1-2v\frac{1-z+\kappa}{1-z}=1-\frac{n^{2}\,v\,z}{(1-z)\,p_{a}n}\;.$$ (70)

we can perform a change of integration variables ${\rm d}\cos\theta_{j,n}\to{\rm d}v$, leading to Eq. (5.189) of Ref. Catani and Seymour (1997).

$$\begin{split}{\rm d}\Phi_{+1}(\tilde{p}_{a},\tilde{p}_{b};\ldots,\tilde{p}_{ij},\ldots;p_{i},p_{j})=&\;\frac{(2\tilde{p}_{ij}\tilde{K})^{1-\varepsilon}}{16\pi^{2}}\,\frac{\big{(}z(1-z)\big{)}^{1-2\varepsilon}}{(1-z+\kappa)^{-\varepsilon}}\,\frac{(\sin^{2}\theta_{j,n})^{-\varepsilon}}{1-z}\,{\rm d}z\,{\rm d}v\,\frac{{\rm d}\Omega_{j,n}^{1-2\varepsilon}}{(4\pi)^{1-2\varepsilon}}\\ =&\;\frac{(2p_{i}n)^{1-\varepsilon}}{16\pi^{2}}\,(1-z)^{-2\varepsilon}\,\bigg{(}\frac{vz}{1-z}\bigg{(}1-\frac{n^{2}vz}{(1-z)\,2p_{i}n}\bigg{)}\bigg{)}^{-\varepsilon}{\rm d}z\,{\rm d}v\,\frac{{\rm d}\Omega_{j,n}^{1-2\varepsilon}}{(2\pi)^{1-2\varepsilon}}\;.\end{split}$$ (71)

A.2 Massive splitting kinematics

In this appendix, we derive the phase-space factorization formula, Eq. (5.11) in Catani et al. (2002) from our generic expression, Eq. (35). We use the definitions of Eq. (24) Catani et al. (2002), and we set $Q^{2}=(\tilde{p}_{ij}+\tilde{K})^{2}$. In addition, we define the scaled masses

$$\hat{\mu}_{i}^{2}=\frac{m_{i}^{2}}{Q^{2}}\;,\qquad\hat{\mu}_{j}^{2}=\frac{m_{j}^{2}}{Q^{2}}\;,\qquad\hat{\mu}_{ij}^{2}=\frac{m_{ij}^{2}}{Q^{2}}\;,\qquad\hat{\kappa}=\frac{K^{2}}{Q^{2}}\;.$$ (72)

The single-emission phase space element is given by

$$\begin{split}{\rm d}\Phi_{+1}(q;\tilde{p}_{ij},\tilde{K};p_{i},p_{j})=\frac{{\rm d}\Phi_{2}(q;p_{ij},K)}{{\rm d}\Phi_{2}(q;\tilde{p}_{ij},\tilde{K})}\frac{{\rm d}p_{ij}^{2}}{2\pi}\,{\rm d}\Phi_{2}(p_{ij};p_{i},p_{j})\;.\end{split}$$ (73)

The decay of $p_{ij}$ is simplest to compute in its rest frame. In this frame, we can write

$$z=\frac{E_{i}^{(ij)}E_{K}^{(ij)}}{p_{ij}K}\Big{(}1-v_{ij,i}v_{ij,k}\cos\theta_{i,ij}\Big{)}=\frac{p_{i}p_{ij}}{p_{ij}^{2}}\Big{(}1-v_{ij,i}v_{ij,k}\cos\theta_{i,ij}\Big{)}\;,$$ (74)

where the velocities are given by

$$v_{ij,i}=\frac{\sqrt{y^{2}(1-\hat{\mu}_{i}^{2}-\hat{\mu}_{j}^{2}-\hat{\kappa})^{2}-4\hat{\mu}_{i}^{2}\hat{\mu}_{j}^{2}}}{y(1-\hat{\mu}_{i}^{2}-\hat{\mu}_{j}^{2}-\hat{\kappa})+2\hat{\mu}_{i}^{2}}\;,\qquad v_{ij,k}=\frac{\sqrt{\big{[}(1-y)(1-\hat{\mu}_{i}^{2}-\hat{\mu}_{j}^{2}-\hat{\kappa})+2\hat{\kappa}\big{]}^{2}-4\hat{\kappa}}}{(1-y)(1-\hat{\mu}_{i}^{2}-\hat{\mu}_{j}^{2}-\hat{\kappa})}$$ (75)

The decay phase space, written in the frame of $p_{ij}$, then reads

$$\begin{split}{\rm d}{\Phi}_{2}(p_{ij};p_{i},p_{j})=&\;\frac{(4\pi^{2})^{\varepsilon}}{16\pi^{2}}\,\frac{(E_{i}^{(ij)}v_{ij,i})^{1-2\varepsilon}}{(p_{ij}^{2})^{1/2}}\,{\rm d}\Omega_{i,ij}^{2-2\varepsilon}=\frac{(4\pi^{2})^{\varepsilon}}{16\pi^{2}}\,\frac{\big{(}p_{i,\perp}^{(ij)}\big{)}^{-2\varepsilon}}{v_{ij,k}}\,{\rm d}z\,{\rm d}\Omega_{i,ij}^{1-2\varepsilon}\;.\end{split}$$ (76)

We can use the factorization Ansatz $p_{i,\perp}^{(ij)\,2}=X(z-z_{-})(z_{+}-z)$, where the physical boundary condition gives the roots of the quadratic $|\cos\theta_{i,ij}|=1$, leading to

$$z_{\pm}=\frac{p_{i}p_{ij}}{p_{ij}^{2}}\Big{(}1\pm v_{ij,i}v_{ij,k}\Big{)}=\frac{y(1-\hat{\mu}_{i}^{2}-\hat{\mu}_{j}^{2}-\hat{\kappa})+2\hat{\mu}_{i}^{2}}{2\big{[}y(1-\hat{\mu}_{i}^{2}-\hat{\mu}_{j}^{2}-\hat{\kappa})+\hat{\mu}_{i}^{2}+\hat{\mu}_{j}^{2}\big{]}}\Big{(}1\pm v_{ij,i}v_{ij,k}\Big{)}\;.$$ (77)

The factor $X$ is determined by equating the transverse momentum at the extremal point $z_{i,\rm max}=p_{i}p_{ij}/p_{ij}^{2}$ to the total three-momentum of $p_{i}$ in the frame of $p_{ij}$. This leads to

$$p_{i,\perp}^{(ij)\,2}=\frac{p_{ij}^{2}}{v_{ij,k}^{2}}(z-z_{-})(z_{+}-z)\;,$$ (78)

and finally

$$\begin{split}{\rm d}{\Phi}_{2}(p_{ij};p_{i},p_{j})=&\;\frac{(4\pi^{2}/Q^{2})^{\varepsilon}}{16\pi^{2}\,v_{ij,k}^{1-2\varepsilon}}\,[y(1-\hat{\mu}_{i}^{2}-\hat{\mu}_{j}^{2}-\hat{\kappa})+\hat{\mu}_{i}^{2}+\hat{\mu}_{j}^{2}]^{-\varepsilon}\big{(}(z-z_{-})(z_{+}-z)\big{)}^{-\varepsilon}\,{\rm d}z\,{\rm d}\Omega_{i,ij}^{1-2\varepsilon}\;.\end{split}$$ (79)

We also have

$${\rm d}p_{ij}^{2}=Q^{2}(1-\hat{\mu}_{i}^{2}-\hat{\mu}_{j}^{2}-\hat{\kappa})\,{\rm d}y\;.$$ (80)

The last missing component of the emission phase space is the ratio of the production to the underlying Born phase-space element. It is most easily computed in the frame of $q=\tilde{p}_{ij}+\tilde{K}$ and results in

$$\frac{{\rm d}\Phi_{2}(q;p_{ij},K)}{{\rm d}\Phi_{2}(q;\tilde{p}_{ij},\tilde{K})}=\bigg{(}\frac{(2p_{ij}K)^{2}\,v_{ij,k}^{2}}{Q^{4}\lambda(1,\hat{\mu}_{ij}^{2},\hat{\kappa})}\bigg{)}^{1/2-\varepsilon}=\bigg{(}\frac{(1-y)^{2}(1-\hat{\mu}_{i}^{2}-\hat{\mu}_{j}^{2}-\hat{\kappa})^{2}\,v_{ij,k}^{2}}{\lambda(1,\hat{\mu}_{ij}^{2},\hat{\kappa})}\bigg{)}^{1/2-\varepsilon}\;.$$ (81)

where the Källen function is given by $\lambda(a,b,c)=(a-b-c)^{2}-4bc$. Combining Eqs. (79)-(81) gives the final result

$$\begin{split}{\rm d}\Phi_{+1}(q;\tilde{p}_{ij},\tilde{K};p_{i},p_{j})=&\;\frac{(Q^{2})^{1-\varepsilon}}{16\pi^{2}}(1-\hat{\mu}_{i}^{2}-\hat{\mu}_{j}^{2}-\hat{\kappa})^{2-2\varepsilon}\lambda(1,\hat{\mu}_{ij}^{2},\hat{\kappa})^{-1/2+\varepsilon}\\ &\qquad\times{\rm d}y\,(1-y)^{1-2\varepsilon}\,\big{[}\,y(1-\hat{\mu}_{i}^{2}-\hat{\mu}_{j}^{2}-\hat{\kappa})+\hat{\mu}_{i}^{2}+\hat{\mu}_{j}^{2}\,\big{]}^{-\varepsilon}\\ &\qquad\times{\rm d}z\,\big{(}(z-z_{-})(z_{+}-z)\big{)}^{-\varepsilon}\;\frac{{\rm d}\Omega_{i,ij}^{1-2\varepsilon}}{(2\pi)^{1-2\varepsilon}}\;.\end{split}$$ (82)