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

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.

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 ¹¹1 Note that in contrast to standard DGLAP notation, we separate the gluon splitting function into two parts, associated with the soft singularities at $z\to 0$ and $z\to 1$.

$$\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.

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).

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}&{\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})=\\ &\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)