aainstitutetext: Instituto de Física Corpuscular, Universitat de València –
Consejo Superior de Investigaciones Científicas, Parc Científic, E-46980
Paterna, Valencia, Spain.bbinstitutetext: Departamento de Física Fundamental e
IUFFyM, Universidad de Salamanca, 37008 Salamanca, Spain.ccinstitutetext:
Escuela de Ciencias, Ingeniería y Diseño, Universidad Europea de Valencia,
Paseo de la Alameda 7, 46010 Valencia, Spain.

# Triple-collinear splittings with massive particles

Prasanna K. Dhani a Germán Rodrigo b,c and German F. R. Sborlini
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We analyze in detail the most singular behaviour of processes involving
triple-collinear splittings with massive particles in the quasi-collinear
limit, and present compact expressions for the splitting amplitudes and the
corresponding splitting kernels at the squared-amplitude level. Our
expressions fully agree with well-known triple-collinear splittings in the
massless limit, which are used as a guide to achieve the final expressions.
These results are important to quantify dominant mass effects in many
observables, and constitute an essential ingredient of current high-precision
computational frameworks for collider phenomenology.

††preprint: IFIC/23-46

##  1 Introduction

Owing to the non-Abelian nature of Quantum Chromodynamics (QCD) due to multi-
gluon interactions, a very relevant topic of study is the singular structure
of scattering amplitudes in various infrared (IR), soft and collinear, limits.
The soft limit is achieved when one or more gluons are emitted in a hard-
scattering process with vanishingly small energies and the collinear limit is
achieved when two or more partons, quarks and gluons, are emitted parallel to
each other i.e. with vanishingly small angles of separation. In these IR
limits, the scattering amplitude for the process becomes singular and it
exhibits a factorized structure Collins:1989gx ; Catani:2011st ;
Forshaw:2012bi . The singular behaviour is given by the QCD factorization
formulae Collins:1989gx  and it is captured by factors which have minimal
dependence on the hard-scattering process under consideration.

In quest for ever increasing precise theoretical predictions for observables
at high-energy colliders such as the CERN’s Large Hadron Collider (LHC), the
framework of perturbation theory plays a pivotal role. To this end, a higher
order theory prediction requires to combine contributions from both virtual
loop and real radiation Feynman diagrams for a meaningful comparison with the
experimental data. Virtual loop diagrams in the IR region lead to singular
configurations in four space-time dimensions. Similarly, real radiation of
soft and/or collinear partons produce kinematical singularities after the
phase space integration over the emitted partons. In the context of
dimensional regularization (DR) Bollini:1972ui ; tHooft:1972tcz ;
Cicuta:1972jf  where one regulates these singularities by performing the
computations in arbitrary space-time dimensions, the understanding of these IR
singular configurations and its isolation for the construction of efficient
methods to evaluate cross sections rely on the general factorization
properties of hard-scattering QCD amplitudes in both soft and collinear
limits. Based on this understanding, general algorithms for the computation of
IR finite jet cross sections at next-to-leading order (NLO) were proposed
through celebrated works in refs. Catani:1996vz ; Frixione:1995ms . Going
forward to the next perturbative orders i.e. next-to-NLO (NNLO) and next-to-
NNLO (N3LO), in the last twenty years or so, a gigantic effort has been put
forth and a variety of different techniques111See the review in refs.
TorresBobadilla:2020ekr ; Caola:2022ayt  and references therein. now exist
with the ability to make finite predictions. All of these methods rely on the
knowledge of underlying QCD factorization of hard-scattering amplitudes in IR
limits222For a review on IR structure of gauge theories, see ref.
Agarwal:2021ais  and references therein..

IR divergences emerging from phase space integration and virtual loop
integration cancel between themselves for a sufficiently inclusive observable
(there still remains collinear divergences in the presence of identified
partons and they are responsible for the evolution of non-perturbative parton
distribution functions (PDFs) in the case of initial states and fragmentation
functions (FFs) in the case of final states). However, close to the boundary
of phase space, real and virtual loop contributions are strongly imbalanced
due to differences in available phase space. As a result, the mechanism of IR
pole cancellation leaves residual effects in the form of large logarithmic
contributions for a wide class of hard-scattering observables. To obtain a
reliable perturbative prediction, these large logarithmic contributions have
to be unveiled to high powers in QCD strong coupling $\alpha_{\rm S}(\equiv
g_{\rm S}^{2}/(4\pi))$, and possibly resummed to all orders in perturbation
theory. Soft and collinear singular factors from the QCD factorization
formulae play a crucial role in achieving this goal through the computation of
observable dependent and process independent integrated quantities. For
instance, a novel methodology was recently developed in ref. Catani:2022sgr
to compute logarithmically-enhanced contributions of collinear origin for the
resummation in transverse momentum and N-jettiness distributions. Needless to
say, soft and collinear factorization formulae form the basis of parton shower
methods which are implemented in Monte Carlo event generators to describe the
exclusive structure of hadronic final states333Several examples of parton
shower generators can be found in the review ref. Hoche:2014rga , and
references therein..

Among IR factorization formulae, collinear factorization Bern:1993qk ;
Bern:1994zx ; Bern:1995ix ; Kosower:1999xi ; Bern:1999ry ; Catani:1999ss ;
Catani:2003vu ; Catani:2011st  is understood to be universal for theory
predictions at colliders. Precisely, the divergent factors called splitting
matrices at the amplitude level and splitting kernels at the squared amplitude
level, only depend on the momenta and quantum numbers of those external
partons involved in the collinear splitting process. In other words, splitting
matrices and splitting kernels do not depend on the specific hard-scattering
process under consideration and its production mechanism444Still, non-Abelian
interactions at higher orders might introduce correlations among collinear and
non-collinear states, breaking strict collinear factorization, as discussed in
refs. Catani:2011st ; Forshaw:2012bi ..

In massless QCD, i.e. QCD with $N_{f}=5$ massless flavors, tree-level double-
parton collinear limits for various splitting processes at
$\mathcal{O}(\alpha_{\rm S})$ are well known Altarelli:1977zs  and they are
proportional to the real radiation contribution to the Altarelli-Parisi
kernels for the LO evolution of the PDFs. At $\mathcal{O}(\alpha_{\rm
S}^{2})$, both triple-parton collinear limit at tree level Campbell:1997hg ;
Catani:1998nv ; Catani:1999ss  and double-parton collinear limit at one-loop
level Bern:1993qk ; Bern:1994zx ; Bern:1995ix ; Bern:1998sc ; Bern:1999ry ;
Kosower:1999rx ; Catani:2011st ; Sborlini:2013jba  are fully known. The
contributions to various splitting processes at $\mathcal{O}(\alpha_{\rm
S}^{3})$ (and beyond) are computed in refs. DelDuca:1999iql ;
Birthwright:2005ak ; Birthwright:2005vi ; Delto:2019asp ; DelDuca:2020vst ;
Catani:2003vu ; Catani:2011st ; Sborlini:2014mpa ; Sborlini:2014kla ;
Badger:2015cxa ; Czakon:2022fqi ; Bern:2004cz ; Badger:2004uk ; Dixon:2019lnw
; Vogt:2004mw ; Vogt:2005dw ; Moch:2004pa ; Hoche:2017iem ; Snigirev:2016uaq .

While there has been a lot of progress in understanding the collinear
factorization of hard-scattering matrix elements with massless partons, the
case with massive partons needs more dedicated efforts in view of modern day
precision colliders. The production of heavy flavors such as bottom and charm
quarks at high-energy colliders has been a topic of interest due to their
relevance in Standard Model (SM) and beyond the SM physics. For better theory
predictions, it is crucial to consider higher-order contributions beyond the
LO in perturbation theory. To this end, it is also important to have control
over explicit IR poles and possibly, large logarithmic contributions as a
result of scale hierarchies. QCD factorization with massive particles plays a
big role in achieving these goals.

It is well known that radiation from massive partons in the collinear limit is
suppressed, which constitutes the so-called dead-cone effect Dokshitzer:1991fd
, and the appearance of explicit collinear poles in perturbative computations
is screened due to the presence of mass. However, they can lead to IR finite
contributions proportional to ln${}^{n}\left(Q^{2}/m^{2}\right)$, where $Q$ is
the typical scale of the hard-scattering process and $m$ is the mass of the
heavy parton. In the kinematical region where $Q^{2}\gg m^{2}$, these large
logarithmic contributions can spoil both the numerical convergence and
reliability of theoretical prediction of a generic observable. The good news
is that these contributions can be computed in a process independent manner
and resummed to all orders in perturbation theory as they are related to
singular behaviour of matrix elements in the $m\rightarrow 0$ limit. The
singular behaviour is controlled by the QCD factorization formulae in the
quasi-collinear limit Catani:2000ef ; Keller:1998tf , in the same way that the
IR divergences are controlled by the soft and collinear factorization
formulae.

At $\mathcal{O}(\alpha_{\rm S})$, double-parton quasi-collinear limit for
various tree-level splitting processes are computed long back in refs.
Catani:2000ef ; Keller:1998tf . Consequently, general algorithms
Dittmaier:1999mb ; Roth:1999kk ; Phaf:2001gc ; Catani:2002hc  were proposed
analogously to the case of massless partons Catani:1996vz ; Frixione:1995ms
for the computation of jet observables. In this paper, we take the first step
in extending the QCD factorization in quasi-collinear limit to
$\mathcal{O}(\alpha_{\rm S}^{2})$ and we study simultaneous collinear limit of
three partons in sight of the importance for accurately predicting dominant
mass effects in phenomenological analysis. We compute the complete set of
splitting matrices at the amplitude level as well as splitting kernels at the
squared amplitude level in arbitrary space-time dimensions, $d=4-2\epsilon$,
and keep its dependence exact. At the squared amplitude level, we also fully
take into account azimuthal correlations, specific to the gluon channel, which
are necessary to construct some general algorithms to perform exact NNLO
computations of jet cross sections involving massive particles.

The outline of the paper is the following. We start in section 2 by reviewing
the basic concepts of collinear factorization and provide a formal definition
of the quasi-collinear limit. In particular, we provide relevant details about
the calculation of splitting amplitudes and polarized/unpolarized splitting
kernels with massive particles. Then, in section 3, we review the available
results for the double-collinear splitting processes with massive partons in
QCD. After that, we provide explicit compact expressions for the QCD triple-
collinear splitting amplitudes and splitting kernels with massive partons in
sections 4 and 5, respectively. We consider all the possible tree-level triple
splitting processes, thus having access to all the required splitting
amplitudes and splitting kernels for a complete double-real collinear
contribution to a NNLO cross-section calculation with massive partons. In
addition, we explain how to recover the corresponding QED results. Finally, in
section 6, we present the conclusions and discuss potential applications of
our results to high-precision QCD phenomenology.

##  2 Multiple collinear and quasi-collinear limits of scattering amplitudes

Let us briefly describe the multiple quasi-collinear limit of scattering
amplitudes at tree-level, and introduce various concepts and definitions. We
consider a generic hard-scattering process in QCD with massive quarks and
anti-quarks, that involves $n$ final-state partons along with any number of
non-coloured particles such as Higgs bosons, electroweak bosons, photons etc.
We denote as ${\cal M}^{(0)}\left(p_{1},\ldots,p_{m},\ldots,p_{n}\right)$ the
corresponding tree-level matrix element, with $p_{i}^{\mu}$ being the four-
momentum associated with $i^{\rm th}$-particle, while the dependence on
additional non-QCD particles through their momenta and quantum numbers is
always understood. We denote the flavor of parton $i$ with $a_{i}$ and we will
avoid indicating the flavor indices explicitly unless they are strictly
required.

The external partons $a_{i}$ are on-shell i.e. $p_{i}^{2}=m_{i}^{2}$. Also,
note that we always define the momenta $p_{i}$’s as outgoing, but the time
component $p_{i}^{0}$ of the four-momentum vector $p_{i}^{\mu}$ is not
constrained to be positive-definite. Different types of physical processes are
represented by applying crossing symmetry to the same matrix element ${\cal
M}^{(0)}$ (for brevity, we also avoid indicating the dependence on momenta
unless it is necessary). For instance, if $p_{i}$ has positive energy, ${\cal
M}^{(0)}$ describes a physical process that produces the parton $i$ in the
final state; if $p_{i}$ has negative energy, the matrix element ${\cal
M}^{(0)}$ describes a physical process produced by the collision of anti-
parton $\bar{i}$ in the initial state.

We are interested in analysing the behaviour of the hard-scattering matrix
element ${\cal M}^{(0)}$ in the kinematical configuration where $m$ out of $n$
external parton momenta become collinear to each other. Without any loss of
generality, we assume them to be $p_{1},\ldots,p_{m}$ and we define their sum
to be $p_{1\ldots m}=p_{1}+\ldots+p_{m}$. The mass of the parent decaying
particle with flavor $a$, which undergoes the splitting, is denoted by
$m_{1\ldots m}$. Since, the parent decaying particle is generally off-shell
($p_{1\dots m}^{2}\neq m_{1\ldots m}^{2}$), we define an on-shell vector

$\displaystyle\widetilde{P}^{\mu}=p_{1\dots m}^{\mu}-\frac{p_{1\ldots
m}^{2}-m_{1\ldots m}^{2}}{2\,n\cdot p_{1\ldots m}}n^{\mu}\,,$ (1)

such that it fulfils $\widetilde{P}^{2}=m_{1\ldots m}^{2}$. Then, we introduce
the Sudakov parametrization for the collinear momenta Catani:2002hc ;
Catani:2011st  as follows

$\displaystyle p_{i}^{\mu}=x_{i}\widetilde{P}^{\mu}+k_{\perp
i}^{\mu}-\frac{k_{\perp i}^{2}+x_{i}^{2}m_{1\ldots
m}^{2}-m_{i}^{2}}{x_{i}}\frac{n^{\mu}}{2\,n\cdot\widetilde{P}}\ ,\qquad i\in
C\equiv\\{1,\ldots,m\\},$ (2)

where the vector $\widetilde{P}$ denotes the collinear direction and the
introduction of an auxiliary light-like vector $n^{\mu}$ ($n^{2}=0$) is
necessary to specify how the collinear direction is approached or,
equivalently, to specify the transverse components $k_{\perp i}$ ($k_{\perp
i}\cdot\widetilde{P}=k_{\perp i}\cdot n=0$, with $k_{\perp i}^{2}<0$).

Since, the collinear limit is invariant under longitudinal boosts along the
direction of total momentum $p_{1\ldots m}$, we define following boost
invariant kinematic variables

$\displaystyle z_{i}$ $\displaystyle=\frac{x_{i}}{\sum_{j=1}^{m}x_{j}},$ (3)
$\displaystyle\tilde{k}_{i}^{\mu}$ $\displaystyle=k_{\perp
i}^{\mu}-z_{i}\sum_{j=1}^{m}k_{\perp j}^{\mu}~{}.$ (4)

Note that these boost invariant variables automatically satisfy the
constraints $\sum_{i=1}^{m}z_{i}=1$ and $\sum_{i=1}^{m}\tilde{k}_{\perp i}=0$
and hence only $2(m-1)$ of them are independent. In terms of the variables
introduced in eqs. (3) and (4), the sub-energies and the dot products are
expressed in the following way

$\displaystyle s_{ij}$
$\displaystyle\equiv(p_{i}+p_{j})^{2}=-z_{i}z_{j}\left(\frac{\tilde{k}_{j}}{z_{j}}-\frac{\tilde{k}_{i}}{z_{i}}\right)^{2}+(z_{i}+z_{j})\left(\frac{m_{i}^{2}}{z_{i}}+\frac{m_{j}^{2}}{z_{j}}\right)~{},$
(5) $\displaystyle\tilde{s}_{ij}$ $\displaystyle\equiv 2p_{i}\cdot
p_{j}=z_{i}z_{j}\left[-\left(\frac{\tilde{k}_{j}}{z_{j}}-\frac{\tilde{k}_{i}}{z_{i}}\right)^{2}+\frac{m_{i}^{2}}{z_{i}^{2}}+\frac{m_{j}^{2}}{z_{j}^{2}}\right]~{},$
(6)

and they are equivalent in the massless limit.

Quasi-collinear limit Catani:2000ef ; Keller:1998tf  is reached when the
transverse momenta $k_{\perp}$ is of $\mathcal{O}(m)$ or small. It turns out
that, by virtue of eqs. (5) and (6), the invariant mass of the splitting
system depends on $\\{k_{\perp i}^{2},m_{i}^{2}\\}$ and $m_{1\ldots m}^{2}$.
So, this allows us to formally define the quasi-collinear kinematical region
by uniformly rescaling these variables, i.e.

$m_{i}\rightarrow\rho\,m_{i}~{},\qquad k_{\perp i}\rightarrow\rho\,k_{\perp
i}~{},\qquad m_{1\ldots m}\rightarrow\rho\,m_{1\ldots m}~{},$ (7)

and taking the limit $\rho\rightarrow 0$. Keeping the most singular terms in
$\rho$ and neglecting sub-leading terms, we obtain the following amplitude
level factorization formula

$|{{\cal
M}^{(0)}\left(p_{1},\ldots,p_{m},\ldots,p_{n}\right)}\rangle\simeq{{\mbox{\boldmath$Sp$}}}^{(0)}_{a_{1}\dots
a_{m}}(p_{1},\ldots,p_{m};\widetilde{P})\,|{\overline{{\cal
M}}^{(0)}(\widetilde{P},p_{m+1},\ldots,p_{n})}\rangle\,,$ (8)

where the factor ${{\mbox{\boldmath$Sp$}}}^{(0)}_{a_{1}\dots
a_{m}}(p_{1},\ldots,p_{m};\widetilde{P})$ is a matrix in color + spin space,
called the splitting matrix that fully embodies the leading singular behaviour
when the parent particle $a$ with momentum $p_{1\ldots m}$ undergoes collinear
splitting into the particles $\\{a_{i}\\}_{i\in C}$, and it only depends on
information carried by the quasi-collinear particles. In eq. (8),
$\overline{{\cal M}}^{(0)}$ is called the reduced matrix element and it is
obtained from the original amplitude by replacing $m$ collinear partons by the
single parent parton of momentum $\widetilde{P}$. Taking the massless limit of
eq. (8), we obtain the well-known collinear factorization formula
Catani:2011st . It is useful to define the colour-projected splitting
matrices, and we have

$Sp_{a_{1}\ldots
a_{m}}^{(0)(c_{1},\ldots,c_{m};c)}\left(p_{1},\ldots,p_{m};\widetilde{P}\right)\equiv\langle
c_{1},\ldots,c_{m}|{{\mbox{\boldmath$Sp$}}}_{a_{1}\ldots
a_{m}}^{(0)}\left(p_{1},\ldots,p_{m};\widetilde{P}\right)|c\rangle\,,$ (9)

where $c$ is the colour of the parent parton and $c_{1},\ldots,c_{m}$ denote
the colour indices of the collinear partons.

An important observation related to eq. (8) is that, the parent parton in the
exact collinear or quasi-collinear limits is on-shell and it carries only
physical degrees of freedom. In this way, the original scattering amplitude
factorizes into two factors, namely the splitting matrix and the reduced
matrix element, both containing physical states as external particles. Notice
that this is valid both for massive and massless particles. The factorization
formula in eq. (8) is gauge independent, however, in order to unveil the
splitting amplitude directly we need to work in a physical gauge. Explicitly,
in physical gauges we have the following relations

$\displaystyle\frac{\imath\,d_{\mu\nu}(p,n)}{p^{2}+\imath 0}$
$\displaystyle\rightarrow$ $\displaystyle\frac{\imath}{p^{2}+\imath
0}\,\sum_{\lambda\ {\rm
phys.}}\,\varepsilon_{\mu}(\widetilde{P},\lambda)\varepsilon^{*}_{\nu}(\widetilde{P},\lambda)$
(10)

for gluons using the polarization tensor $d_{\mu\nu}(p,n)$, and

$\displaystyle\frac{\imath\,(\gamma^{\mu}p_{\mu}+m)}{p^{2}-m^{2}+\imath 0}$
$\displaystyle\rightarrow$ $\displaystyle\frac{\imath}{p^{2}-m^{2}+\imath
0}\,\sum_{\lambda\ {\rm
phys.}}\,u(\widetilde{P},\lambda)\bar{u}(\widetilde{P},\lambda)$ (11)

for quarks, where $\lambda$ denotes the polarization (or spin).

Besides the splitting matrices, we are also interested in the polarized and
unpolarized splitting kernels, which govern the collinear and quasi-collinear
limits of the squared matrix elements. In order to define the polarized
kernels, let us start from eq. (8) and consider a gluon-initiated splitting
process. Stripping the polarization vector from the splitting amplitude, we
have

$\displaystyle|{{\cal
M}^{(0)}\left(p_{1},\ldots,p_{m},\ldots,p_{n}\right)}\rangle\simeq$
$\displaystyle\sum_{\lambda\,{\rm
phys.}}\,{{\mbox{\boldmath$Sp$}}}^{(0),\mu}_{a_{1}\dots
a_{m}}(p_{1},\ldots,p_{m};\widetilde{P})\,\varepsilon_{\mu}(\widetilde{P},\lambda)\,$
(12) $\displaystyle\times|{{\cal
M}^{(0)}(\widetilde{P}^{-\lambda},p_{m+1},\ldots,p_{n})}\rangle\,,$

and computing the square, we obtain

$\displaystyle|{\cal
M}^{(0)}\left(p_{1},\ldots,p_{m},\ldots,p_{n}\right)|^{2}\simeq$
$\displaystyle\sum_{\lambda,\lambda^{\prime}\,{\rm phys.}}\,\langle{{\cal
M}^{(0)}(\widetilde{P}^{-\lambda^{\prime}},p_{m+1},\ldots,p_{n})}|\,\varepsilon_{\mu}(\widetilde{P},\lambda)\varepsilon^{*}_{\nu}(\widetilde{P},\lambda^{\prime})\,$
(13) $\displaystyle\times\left[{{\mbox{\boldmath$Sp$}}}^{(0),\nu}_{a_{1}\dots
a_{m}}(p_{1},\ldots,p_{m};\widetilde{P})\right]^{\dagger}{{\mbox{\boldmath$Sp$}}}^{(0),\mu}_{a_{1}\dots
a_{m}}(p_{1},\ldots,p_{m};\widetilde{P})\,\,$ $\displaystyle\times|{{\cal
M}^{(0)}(\widetilde{P}^{-\lambda},p_{m+1},\ldots,p_{n})}\rangle\,.$

The second line of eq. (13) contains the product of two amputated splitting
amplitudes, which leads us to define

$\hat{P}^{(0),\,\mu\nu}_{a_{1}\dots a_{m}}\equiv\left(\frac{s_{1\ldots
m}-m_{1\ldots m}^{2}}{8\pi\mu_{0}^{2\epsilon}\alpha_{\rm
S}}\right)^{m-1}\left[{{\mbox{\boldmath$Sp$}}}^{(0),\nu}_{a_{1}\dots
a_{m}}(p_{1},\ldots,p_{m};\widetilde{P})\right]^{\dagger}{{\mbox{\boldmath$Sp$}}}^{(0),\mu}_{a_{1}\dots
a_{m}}(p_{1},\ldots,p_{m};\widetilde{P})\,,$ (14)

that corresponds to the polarized splitting kernel, with the proper
normalization and $\mu_{0}$ denotes the DR scale. Note that analogous to the
splitting matrix ${{\mbox{\boldmath$Sp$}}}^{(0)}$ at the amplitude level, the
splitting kernel $\hat{P}^{(0)}$ is a matrix in color + spin space and it
encodes the leading singular behaviour of the squared amplitude in the quasi-
collinear limit. In eq. (14), we have introduced $s_{1\ldots m}\equiv
p_{1\ldots m}^{2}$ and summed over final-state colours and spins, averaged
over only colours of the parent decaying parton. By averaging over
spins/polarizations of the parent parton, we obtain the $m$-parton unpolarized
splitting kernel $\langle\hat{P}^{(0)}_{a_{1}\ldots a_{m}}\rangle$, which is a
generalization of the customary Altarelli-Parisi LO splitting function.
Explicitly, by fixing the normalization, the tree-level unpolarized splitting
kernels are given by

$\langle\hat{P}_{a_{1}\cdots a_{m}}^{(0)}\rangle=\left(\frac{s_{1\ldots
m}-m_{1\ldots m}^{2}}{8\pi\mu_{0}^{2\epsilon}\alpha_{\rm
S}}\right)^{m-1}\,{\overline{|{{\mbox{\boldmath$Sp$}}}_{a_{1}\dots
a_{m}}^{(0)}|^{2}}}~{}.$ (15)

In the quark-initiated case, we have the relation

$\langle\hat{P}_{a_{1}\cdots
a_{m}}^{(0)}\rangle=\frac{1}{2}\,\delta^{ss^{\prime}}\,\hat{P}^{(0),\,ss^{\prime}}_{a_{1}\dots
a_{m}}\,,$ (16)

because of helicity conservation. The spin indices $s$ and $s^{\prime}$ are of
the parent collinear quark (anti-quark) in the amplitude and conjugate
amplitude, respectively. In the gluon-initiated splitting process, the
unpolarized kernel is obtained by contracting the polarized one with the
polarization tensor of the parent parton, i.e. $d_{\mu\nu}(\widetilde{P},n)$,
and dividing by $d-2$ (number of initial-state polarizations in
$d=4-2\epsilon$ space-time dimensions). Explicitly, it is given by

$\displaystyle\langle\hat{P}^{(0)}_{a_{1}\ldots
a_{m}}\rangle=\frac{d_{\mu\nu}(\widetilde{P};n)}{d-2}\hat{P}^{(0),\mu\nu}_{a_{1}\ldots
a_{m}}.$ (17)

As we explained before, we shall choose a physical gauge to calculate the
splitting amplitudes. Therefore, we use the light-cone gauge (LCG) where

$d_{\mu\nu}(k,n)=-g_{\mu\nu}+\frac{k_{\mu}n_{\nu}+n_{\mu}k_{\nu}}{n\cdot
k}\,,$ (18)

represents the physical polarization tensor of a gauge vector boson (gluon or
photon) with momentum $k$. The auxiliary gauge vector $n$ in eq. (18) is taken
identical to the light-like vector $n$ introduced in eq. (1).

Finally, we compute the matrix elements associated with the evaluation of
massive splitting matrices and splitting kernels by using two different
frameworks. On the one side, we used QGRAF Nogueira:1991ex  to generate all
the necessary Feynman diagrams and an in-house code written in Form
Kuipers:2012rf  to perform Dirac and color algebra simplification. On the
other, we relied on the automatized packages FeynCalc Shtabovenko:2016sxi ;
Shtabovenko:2020gxv  and FeynArts Hahn:2000kx , exploiting their functions to
simplify Dirac’s chains. The results obtained are in total agreement within
both methods. For the extraction of the splitting amplitudes and kernels, we
follow the techniques described in refs. Catani:1999ss ; Catani:2002hc ;
Catani:2011st ; Czakon:2022fqi . In the following sections, we present
perturbative results for various tree-level splitting processes up to
$\mathcal{O}(\alpha_{\rm S}^{2})$ in the quasi-collinear limit described in
eq. (7).

##  3 Massive double-collinear splitting amplitudes and splitting kernels

We recall in this section the known results for the double-collinear splitting
matrices and kernels with massive quarks. There are two non-vanishing
splitting processes, i.e.

$\displaystyle Q\rightarrow Q_{1}g_{2}~{},$ (19) $\displaystyle g\rightarrow
Q_{1}\bar{Q}_{2}~{}.$ (20)

The remaining processes are related by charge conjugation invariance of QCD.
Their corresponding splitting matrices are

$\displaystyle{{\mbox{\boldmath$Sp$}}}_{Q_{1}g_{2}}^{(0)}$
$\displaystyle=g_{\rm
s}\mu_{0}^{\epsilon}\boldsymbol{t}^{c_{2}}_{c_{1}c}\frac{1}{\tilde{s}_{12}}\overline{u}(p_{1})\not{\varepsilon}(p_{2})u(\widetilde{P})~{},$
(21) $\displaystyle{{\mbox{\boldmath$Sp$}}}_{Q_{1}\bar{Q}_{2}}^{(0)}$
$\displaystyle=g_{\rm
s}\mu_{0}^{\epsilon}\boldsymbol{t}^{c}_{c_{1}c_{2}}\frac{1}{s_{12}}\overline{u}(p_{1})\not{\varepsilon}^{*}(\widetilde{P})v(p_{2})~{},$
(22)

where $\boldsymbol{t}^{c}_{c_{1}c_{2}}$ are the customary Gell-Mann colour
matrices in the fundamental representation and the Dirac spinors are
represented by $u(p)$ and $v(p)$. At the squared amplitude level, the
splitting kernels are given by Catani:2000ef ; Keller:1998tf

$\displaystyle\hat{P}^{(0),ss^{\prime}}_{Q_{1}g_{2}}$
$\displaystyle=\delta^{ss^{\prime}}C_{F}\left[\frac{1+z^{2}}{1-z}-\epsilon(1-z)-2\frac{m_{Q}^{2}}{\tilde{s}_{12}}\right],$
(23) $\displaystyle\hat{P}^{(0),\mu\nu}_{Q_{1}\bar{Q}_{2}}$
$\displaystyle=T_{R}\left[-g^{\mu\nu}-4\frac{k_{\perp}^{\mu}k_{\perp}^{\nu}}{s_{12}}\right],$
(24)

where $C_{F}=(N_{c}^{2}-1)/(2N_{c})$ is the Casimir in the fundamental
representation of SU($N_{c}$) QCD with $N_{c}$ colors. The corresponding
Casimir in the adjoint representation is given by $C_{A}=N_{c}$ and
$T_{R}=1/2$. The momentum fraction $z$ in eq. (23) is $z=z_{1}=1-z_{2}$ and
the transverse momentum vector $k_{\perp}^{\mu}$ in eq. (24) is
$k_{\perp}^{\mu}=k_{\perp 1}^{\mu}=-k_{\perp 2}^{\mu}$. As we anticipated in
section 2, the polarized splitting kernels are diagonal in the spin of the
parent quark (anti-quark) but are not diagonal in the polarization of the
parent gluon. To this end, we also obtain the corresponding unpolarized
splitting kernels by using eqs. (16-17) and they are as follows

$\displaystyle\langle\hat{P}^{(0)}_{Q_{1}g_{2}}\rangle$
$\displaystyle=C_{F}\left[\frac{1+z^{2}}{1-z}-\epsilon(1-z)-2\frac{m_{Q}^{2}}{\tilde{s}_{12}}\right],$
(25) $\displaystyle\langle\hat{P}^{(0)}_{Q_{1}\bar{Q}_{2}}\rangle$
$\displaystyle=T_{R}\left[1-\frac{2}{1-\epsilon}\left(z(1-z)-\frac{m_{Q}^{2}}{s_{12}}\right)\right].$
(26)

It is worth noticing that the mass dependence in
$\hat{P}^{(0),\mu\nu}_{Q_{1}\bar{Q}_{2}}$ appears implicitly inside the sub-
energy $s_{12}$; the functional structure is analogous to the massless case.
In the case of $\hat{P}^{(0),ss^{\prime}}_{Q_{1}g_{2}}$, besides a
contribution sharing the same functional dependence as in the massless case,
there is an additional term proportional to $m^{2}_{Q}$. We will see similar
features in the massive triple-collinear splitting kernels.

##  4 Massive triple-collinear splitting amplitudes

There are four non-vanishing splitting processes, which are given by

$\displaystyle Q$ $\displaystyle\rightarrow\bar{Q}_{1}Q_{2}Q_{3}~{},$ (27)
$\displaystyle Q$
$\displaystyle\rightarrow\bar{Q}_{1}^{\prime}Q_{2}^{\prime}Q_{3}~{},$ (28)
$\displaystyle Q$ $\displaystyle\rightarrow g_{1}g_{2}Q_{3}~{},$ (29)
$\displaystyle g$ $\displaystyle\rightarrow g_{1}Q_{2}\bar{Q}_{3}~{}.$ (30)

The quark-antiquark pair $\bar{Q}_{1}^{\prime}Q_{2}^{\prime}$ is of different
flavour than $Q$ and has a mass $m_{Q^{\prime}}$. All the other quarks and
antiquarks have mass $m_{Q}$. In the following sections, we present results
for the splitting matrices and splitting kernels corresponding to the
processes in eqs. (27-30).

###  4.1 Splitting amplitudes with four massive quarks

We consider first the splitting amplitude of four massive quarks of the same
flavour, eq. (27). In this case there are two colour structures and they are
related by the permutation $2\leftrightarrow 3$:

$\displaystyle\mathbf{C}_{1}:T_{R}\left(\delta_{c_{1}c_{3}}\delta_{cc_{2}}-\frac{1}{N_{c}}\delta_{cc_{3}}\delta_{c_{1}c_{2}}\right)~{},$
(31)
$\displaystyle\mathbf{C}_{2}:T_{R}\left(\delta_{c_{1}c_{2}}\delta_{cc_{3}}-\frac{1}{N_{c}}\delta_{cc_{2}}\delta_{c_{1}c_{3}}\right)~{}.$
(32)

The l.h.s. of these expressions refer to operators, whilst the r.h.s.
correspond to their projection according eq. (9). There are six distinct spin
structures out of which three are related to the other three by the
permutation $2\leftrightarrow 3$:

$\displaystyle\mathbf{S}_{1}:\bar{u}(p_{3})\gamma^{\mu}u(\widetilde{P})\bar{u}(p_{2})\gamma_{\mu}v(p_{1})~{},\quad\,\,\mathbf{S}_{2}:\bar{u}(p_{3})\not{p}_{2}u(\widetilde{P})\bar{u}(p_{2})\not{n}v(p_{1})~{},$
(33)
$\displaystyle\mathbf{S}_{3}:\bar{u}(p_{3})\not{p}_{1}u(\widetilde{P})\bar{u}(p_{2})\not{n}v(p_{1})~{},\qquad\mathbf{S}_{4}:\bar{u}(p_{2})\gamma^{\mu}u(\widetilde{P})\bar{u}(p_{3})\gamma_{\mu}v(p_{1})~{},$
(34)
$\displaystyle\mathbf{S}_{5}:\bar{u}(p_{2})\not{p}_{3}u(\widetilde{P})\bar{u}(p_{3})\not{n}v(p_{1})~{},\qquad\mathbf{S}_{6}:\bar{u}(p_{2})\not{p}_{1}u(\widetilde{P})\bar{u}(p_{3})\not{n}v(p_{1})~{}.$
(35)

The splitting amplitude for the process is given by

$\displaystyle{{\mbox{\boldmath$Sp$}}}_{\bar{Q}_{1}Q_{2}Q_{3}}^{(0)}=\frac{g_{\rm
s}^{2}\mu_{0}^{2\epsilon}}{s_{123}-m_{Q}^{2}}\left[\frac{\mathbf{C}_{1}}{s_{12}}\left(-\mathbf{S}_{1}+\frac{\mathbf{S}_{2}+\mathbf{S}_{3}}{n\cdot
p_{12}}\right)-\frac{\mathbf{C}_{2}}{s_{13}}\left(-\mathbf{S}_{4}+\frac{\mathbf{S}_{5}+\mathbf{S}_{6}}{n\cdot
p_{13}}\right)\right]~{},$ (36)

where $s_{123}\equiv(p_{1}+p_{2}+p_{3})^{2}$ and $p_{ij}\equiv p_{i}+p_{j}$.
If, instead, we consider the splitting process where the quark-antiquark pair
$\bar{Q}_{1}^{\prime}Q_{2}^{\prime}$ is of different flavour and mass than the
parent quark $Q$, eq. (28), the splitting amplitude is simplified to

$\displaystyle{{\mbox{\boldmath$Sp$}}}_{\bar{Q}_{1}^{\prime}Q_{2}^{\prime}Q_{3}}^{(0)}=\frac{g_{\rm
s}^{2}\mu_{0}^{2\epsilon}}{s_{123}-m_{Q}^{2}}\frac{\mathbf{C}_{1}}{s_{12}}\left(-\mathbf{S}_{1}+\frac{\mathbf{S}_{2}+\mathbf{S}_{3}}{n\cdot
p_{12}}\right)~{}.$ (37)

The explicit dependence on the mass $m_{Q^{\prime}}$ is hiding in $s_{123}$
and $s_{12}$, and the implicit dependence is through the momenta $p_{1}$ and
$p_{2}$.

###  4.2 $Q\rightarrow g_{1}g_{2}Q_{3}$

In this case there are two colour structures and they are related by the
permutation $1\leftrightarrow 2$:

$\displaystyle\mathbf{C}_{1}:\left(\boldsymbol{t}^{c_{2}}\boldsymbol{t}^{c_{1}}\right)_{c_{3}c}~{},\quad\mathbf{C}_{2}:\left(\boldsymbol{t}^{c_{1}}\boldsymbol{t}^{c_{2}}\right)_{c_{3}c}~{},$
(38)

and eight distinct spin structures

$\displaystyle\mathbf{S}_{1}:\bar{u}(p_{3})\not{p}_{1}\not{\varepsilon}(p_{1})\not{\varepsilon}(p_{2})u(\widetilde{P})~{},\qquad\;\;\;\mathbf{S}_{2}:\bar{u}(p_{3})\not{p}_{2}\not{\varepsilon}(p_{2})\not{\varepsilon}(p_{1})u(\widetilde{P})~{},$
$\displaystyle\mathbf{S}_{3}:\bar{u}(p_{3})\not{\varepsilon}(p_{2})u(\widetilde{P})\,p_{3}\cdot\varepsilon(p_{1})~{},\qquad\mathbf{S}_{4}:\bar{u}(p_{3})\not{\varepsilon}(p_{2})u(\widetilde{P})\,p_{2}\cdot\varepsilon(p_{1})~{},$
$\displaystyle\mathbf{S}_{5}:\bar{u}(p_{3})\not{\varepsilon}(p_{1})u(\widetilde{P})\,p_{3}\cdot\varepsilon(p_{2})~{},\qquad\mathbf{S}_{6}:\bar{u}(p_{3})\not{\varepsilon}(p_{1})u(\widetilde{P})\,p_{1}\cdot\varepsilon(p_{2})~{},$
$\displaystyle\mathbf{S}_{7}:\bar{u}(p_{3})\not{p}_{2}u(\widetilde{P})\,\varepsilon(p_{1})\cdot\varepsilon(p_{2})~{},\qquad\mathbf{S}_{8}:\bar{u}(p_{3})\not{p}_{1}u(\widetilde{P})\,\varepsilon(p_{1})\cdot\varepsilon(p_{2})~{}.$
(39)

Notice that these spin structures fulfil
$\mathbf{S}_{1}\leftrightarrow\mathbf{S}_{2}$,
$\mathbf{S}_{3}\leftrightarrow\mathbf{S}_{5}$,
$\mathbf{S}_{4}\leftrightarrow\mathbf{S}_{6}$ and
$\mathbf{S}_{7}\leftrightarrow\mathbf{S}_{8}$ under the transformation
$1\leftrightarrow 2$. The splitting amplitude for the process is given by

$\displaystyle{{\mbox{\boldmath$Sp$}}}^{(0)}_{g_{1}g_{2}Q_{3}}$
$\displaystyle=\frac{g_{s}^{2}\mu_{0}^{2\epsilon}}{\tilde{s}_{123}}\Bigg{[}\mathbf{C}_{1}\left(\frac{\mathbf{S}_{2}-2\mathbf{S}_{5}}{\tilde{s}_{23}}\right)+\mathbf{C}_{2}\left(\frac{\mathbf{S}_{1}-2\mathbf{S}_{3}}{\tilde{s}_{13}}\right)$
$\displaystyle+(\mathbf{C}_{2}-\mathbf{C}_{1})\frac{2}{\tilde{s}_{12}}\left(\mathbf{S}_{4}+\frac{n\cdot
p_{2}}{n\cdot p_{12}}\mathbf{S}_{8}-\mathbf{S}_{6}-\frac{n\cdot p_{1}}{n\cdot
p_{12}}\mathbf{S}_{7}\right)\Bigg{]}~{},$ (40)

being the full expression symmetric under the exchange $1\leftrightarrow 2$,
as expected. Here, we define
$\tilde{s}_{123}=\tilde{s}_{12}+\tilde{s}_{13}+\tilde{s}_{23}$, following the
notation introduced in eq. (6).

###  4.3 $g\rightarrow g_{1}Q_{2}\bar{Q}_{3}$

In this case there are two colour structures given by:

$\displaystyle\mathbf{C}_{1}:\left(\boldsymbol{t}^{c}\boldsymbol{t}^{c_{1}}\right)_{c_{2}c_{3}}~{},\qquad\mathbf{C}_{2}:\left(\boldsymbol{t}^{c_{1}}\boldsymbol{t}^{c}\right)_{c_{2}c_{3}}~{},$
(41)

and six distinct spin structures

$\displaystyle\mathbf{S}_{1}:\bar{u}(p_{2})\not{p}_{1}\not{\varepsilon}(p_{1})\not{\varepsilon}^{*}(\widetilde{P})v(p_{3})~{},\qquad\quad\mathbf{S}_{2}:\bar{u}(p_{2})\not{\varepsilon}^{*}(\widetilde{P})v(p_{3})\,p_{2}\cdot\varepsilon(p_{1})~{},$
$\displaystyle\mathbf{S}_{3}:\bar{u}(p_{2})\not{\varepsilon}^{*}(\widetilde{P})v(p_{3})\,p_{3}\cdot\varepsilon(p_{1})~{},\quad\quad\,\mathbf{S}_{4}:\bar{u}(p_{2})\not{\varepsilon}(p_{1})v(p_{3})\,p_{1}\cdot\varepsilon^{*}(\widetilde{P})~{},$
$\displaystyle\mathbf{S}_{5}:\bar{u}(p_{2})\not{p}_{1}(p_{1})v(p_{3})\,\varepsilon^{*}(\widetilde{P})\cdot\varepsilon(p_{1})~{},\,\,\,\mathbf{S}_{6}:\bar{u}(p_{2})\not{n}v(p_{3})\,\varepsilon^{*}(\widetilde{P})\cdot\varepsilon(p_{1})~{}.$
(42)

However, these spin structures can be grouped according to

$\mathbf{S}^{\prime}_{1}=\mathbf{S}_{1}-\mathbf{S}_{2}\,,\quad\mathbf{S}^{\prime}_{2}=2\,(\mathbf{S}_{2}+\mathbf{S}_{3}+\mathbf{S}_{4}-\mathbf{S}_{5})\,,$
(43)

in such a way that the splitting amplitude for the process is written as

$\displaystyle{{\mbox{\boldmath$Sp$}}}^{(0)}_{g_{1}Q_{2}\bar{Q}_{3}}$
$\displaystyle=\frac{g_{s}^{2}\mu_{0}^{2\epsilon}}{s_{123}}\Bigg{[}\mathbf{C}_{1}\left(\frac{\mathbf{S}^{\prime}_{1}+\mathbf{S}^{\prime}_{2}}{\tilde{s}_{13}}\right)+\mathbf{C}_{2}\,\frac{\mathbf{S}^{\prime}_{1}}{\tilde{s}_{12}}+(\mathbf{C}_{1}-\mathbf{C}_{2})\left(\frac{\mathbf{S}^{\prime}_{2}}{s_{23}}+\frac{s_{123}}{s_{23}}\frac{\mathbf{S}_{6}}{n\cdot
p_{23}}\right)\Bigg{]}.$ (44)

It is worth appreciating that, in the massless limit, we recover the same
expressions presented in ref. Czakon:2022fqi : in fact, we explicitly defined
the spinor bases to ease the comparison and check our results. Besides, we
appreciate that both the massive and massless splitting amplitudes share the
same spinor and colour structures.

##  5 Massive triple-collinear splitting kernels

In order to obtain compact expressions, we exploit the results provided in
ref. Catani:1999ss  for the massless triple-collinear kernels, including the
definition of the kinematical variable

$t_{ij,k}=2\frac{z_{i}\tilde{s}_{jk}-z_{j}\tilde{s}_{ik}}{z_{i}+z_{j}}+\frac{z_{i}-z_{j}}{z_{i}+z_{j}}\tilde{s}_{ij}\,.$
(45)

The unpolarized splitting kernel involving two quark-antiquark pairs of
different flavours is given by

$\displaystyle\langle\hat{P}^{(0)}_{\bar{Q}_{1}^{\prime}Q_{2}^{\prime}Q_{3}}\rangle$
$\displaystyle=C_{F}T_{R}\Bigg{\\{}\frac{\tilde{s}_{12}\tilde{s}_{123}}{2s_{12}^{2}}\left[-\frac{t_{12,3}^{2}}{\tilde{s}_{12}\tilde{s}_{123}}+\frac{4z_{3}+(z_{1}-z_{2})^{2}}{1-z_{3}}+(1-2\epsilon)\left(z_{1}+z_{2}-\frac{\tilde{s}_{12}}{\tilde{s}_{123}}\right)\right]$
$\displaystyle+\frac{2m_{Q^{\prime}}^{2}}{s_{12}^{2}}\bigg{[}\frac{z_{3}\tilde{s}_{123}}{(1-z_{3})^{2}}(1+2z_{3}-3z_{3}^{2}+4z_{1}z_{2})-\frac{\tilde{s}_{23}}{1-z_{3}}(2-3z_{1}-5z_{2}+z_{1}^{2}+z_{2}^{2})$
$\displaystyle-\frac{\tilde{s}_{13}}{1-z_{3}}(2-5z_{1}-3z_{2}+z_{1}^{2}+z_{2}^{2})-\epsilon\Big{(}\tilde{s}_{123}(1-z_{3})-\tilde{s}_{12}(1+z_{3})\Big{)}\bigg{]}-2\frac{m_{Q}^{2}\tilde{s}_{12}}{s_{12}^{2}}$
$\displaystyle+\frac{4m_{Q^{\prime}}^{4}}{s_{12}^{2}}z_{3}\left[\epsilon+\frac{2z_{1}z_{2}}{(1-z_{3})^{2}}+\frac{2z_{3}}{1-z_{3}}\right]-4\frac{m_{Q}^{2}m_{Q^{\prime}}^{2}}{s_{12}^{2}}\Bigg{\\}}~{}.$
(46)

The splitting kernel involving two quark-antiquark pairs of the same flavour
is decomposed as

$\displaystyle\langle\hat{P}^{(0)}_{\bar{Q}_{1}Q_{2}Q_{3}}\rangle=\left[\langle\hat{P}^{(0)}_{\bar{Q}_{1}^{\prime}Q_{2}^{\prime}Q_{3}}\rangle\Big{|}_{m_{Q^{\prime}}=m_{Q}}+(2\leftrightarrow
3)\right]+\langle\hat{P}_{\bar{Q}_{1}Q_{2}Q_{3}}^{\rm{(0,id)}}\rangle~{},$
(47)

where the interference contribution is given by

$\displaystyle\langle\hat{P}_{\bar{Q}_{1}Q_{2}Q_{3}}^{\rm{(0,id)}}\rangle$
$\displaystyle=C_{F}\left(C_{F}-\frac{C_{A}}{2}\right)\Bigg{\\{}\frac{\tilde{s}_{12}\tilde{s}_{13}}{s_{12}s_{13}}\Bigg{[}(1-\epsilon)\left(\frac{2\tilde{s}_{23}}{\tilde{s}_{12}}-\epsilon\right)+\frac{\tilde{s}_{123}}{\tilde{s}_{12}}\bigg{\\{}\frac{1+z_{1}^{2}}{1-z_{2}}-\frac{2z_{2}}{1-z_{3}}$
$\displaystyle-\epsilon\left(\frac{(1-z_{3})^{2}}{1-z_{2}}+1+z_{1}-\frac{2z_{2}}{1-z_{3}}\right)-\epsilon^{2}(1-z_{3})\bigg{\\}}-\frac{\tilde{s}_{123}^{2}}{\tilde{s}_{12}\tilde{s}_{13}}\frac{z_{1}}{2}\bigg{\\{}\frac{1+z_{1}^{2}}{(1-z_{2})(1-z_{3})}$
$\displaystyle-\epsilon\left(1+2\frac{1-z_{2}}{1-z_{3}}\right)-\epsilon^{2}\bigg{\\}}\Bigg{]}+\frac{m_{Q}^{2}}{s_{12}s_{13}}\Bigg{[}2\tilde{s}_{23}\frac{3-7z_{2}+z_{2}^{2}+z_{2}^{3}+4z_{2}^{2}z_{3}}{(1-z_{2})(1-z_{3})}$
$\displaystyle+2\tilde{s}_{12}\frac{(1-2z_{2})^{2}+(1-2z_{3})^{2}-z_{3}(1-z_{3}^{2})-z_{2}(1+2z_{3})+z_{2}z_{3}(2z_{2}+3z_{3})}{(1-z_{2})(1-z_{3})}$
$\displaystyle-2\epsilon\tilde{s}_{12}\frac{2(1-z_{2})^{2}+(1-z_{3})^{2}-z_{3}(1-z_{3}^{2})-2z_{2}z_{3}(1+z_{1})}{(1-z_{2})(1-z_{3})}$
$\displaystyle-2\epsilon\tilde{s}_{23}\frac{1-z_{2}^{2}+z_{2}(1-3z_{3})}{1-z_{3}}+\epsilon^{2}\Big{(}2\tilde{s}_{12}(1-z_{3})-\tilde{s}_{23}(1-z_{1})\Big{)}\Bigg{]}$
$\displaystyle+\frac{2m_{Q}^{4}}{s_{12}s_{13}}\Bigg{[}2\frac{(1-2z_{2})^{2}-z_{2}(1-z_{2}z_{3})}{(1-z_{2})(1-z_{3})}-\epsilon\frac{(1-2z_{2})(3-2z_{2}-z_{2}z_{3})}{(1-z_{2})(1-z_{3})}+\epsilon^{2}\Bigg{]}\Bigg{\\}}$
$\displaystyle+(2\leftrightarrow 3)~{}.$ (48)

The splitting kernel involving two gluons, and a massive quark as parent
parton, admits the following decomposition

$\displaystyle\langle\hat{P}^{(0)}_{g_{1}g_{2}Q_{3}}\rangle=C_{F}^{2}\langle\hat{P}_{g_{1}g_{2}Q_{3}}^{\rm{(0,ab)}}\rangle+C_{A}C_{F}\langle\hat{P}_{g_{1}g_{2}Q_{3}}^{\rm{(0,nab)}}\rangle~{},$
(49)

where

$\displaystyle\langle\hat{P}_{g_{1}g_{2}Q_{3}}^{\rm{(0,ab)}}\rangle$
$\displaystyle=\Bigg{\\{}\frac{\tilde{s}_{123}^{2}z_{3}}{2\tilde{s}_{13}\tilde{s}_{23}}\left[\frac{1+z_{3}^{2}}{z_{1}z_{2}}-\epsilon\frac{z_{1}^{2}+z_{2}^{2}}{z_{1}z_{2}}-\epsilon(1+\epsilon)\right]+\frac{\tilde{s}_{123}}{\tilde{s}_{13}}\bigg{[}\frac{z_{3}(1-z_{1})+(1-z_{2})^{3}}{z_{1}z_{2}}$
$\displaystyle-\epsilon(z_{1}^{2}+z_{1}z_{2}+z_{2}^{2})\frac{1-z_{2}}{z_{1}z_{2}}+\epsilon^{2}(1+z_{3})\bigg{]}+(1-\epsilon)\left[\epsilon-(1-\epsilon)\frac{\tilde{s}_{23}}{\tilde{s}_{13}}\right]$
$\displaystyle+\frac{m_{Q}^{2}}{\tilde{s}_{23}}\bigg{[}\frac{2(z_{2}-2z_{3}+z_{2}^{2}+2z_{3}^{2}+3z_{2}z_{3})}{z_{1}z_{2}}-\frac{2(2z_{2}+4z_{3}-z_{1}z_{3})}{z_{1}}$
$\displaystyle-\frac{2\tilde{s}_{123}}{\tilde{s}_{23}}\frac{1+z_{2}^{2}+z_{3}^{2}+2z_{2}z_{3}}{z_{1}}-\frac{\tilde{s}_{12}}{\tilde{s}_{13}}\frac{z_{2}(1-z_{3})+4z_{3}}{z_{2}}+2\epsilon
z_{1}\frac{\tilde{s}_{123}}{\tilde{s}_{23}}$
$\displaystyle+\epsilon(1-z_{3})\frac{\tilde{s}_{123}}{\tilde{s}_{13}}\bigg{]}+\frac{4m_{Q}^{4}}{\tilde{s}_{23}^{2}}+\frac{4m_{Q}^{4}}{\tilde{s}_{13}\tilde{s}_{23}}\Bigg{\\}}+(1\leftrightarrow
2)~{},$ (50)

and

$\displaystyle\langle\hat{P}_{g_{1}g_{2}Q_{3}}^{\rm{(0,nab)}}\rangle$
$\displaystyle=\Bigg{\\{}(1-\epsilon)\left(\frac{t_{12,3}^{2}}{4\tilde{s}_{12}^{2}}+\frac{1}{4}-\frac{\epsilon}{2}\right)+\frac{\tilde{s}_{23}^{2}}{2\tilde{s}_{12}\tilde{s}_{13}}\bigg{[}\frac{(1-z_{3})^{2}(1-\epsilon)+2z_{3}}{z_{2}}$
$\displaystyle+\frac{z_{2}^{2}(1-\epsilon)+2(1-z_{2})}{1-z_{3}}\bigg{]}-\frac{\tilde{s}_{123}^{2}}{4\tilde{s}_{13}\tilde{s}_{23}}z_{3}\left[\frac{(1-z_{3})^{2}(1-\epsilon)+2z_{3}}{z_{1}z_{2}}+\epsilon(1-\epsilon)\right]$
$\displaystyle+\frac{\tilde{s}_{123}}{2\tilde{s}_{12}}\left[(1-\epsilon)\frac{z_{1}(2-2z_{1}+z_{1}^{2})-z_{2}(6-6z_{2}+z_{2}^{2})}{z_{2}(1-z_{3})}+2\epsilon\frac{z_{3}(z_{1}-2z_{2})-z_{2}}{z_{2}(1-z_{3})}\right]$
$\displaystyle+\frac{\tilde{s}_{123}}{2\tilde{s}_{13}}\bigg{[}(1-\epsilon)\frac{(1-z_{2})^{2}+z_{3}^{2}-z_{2}}{z_{2}(1-z_{3})}-\epsilon\left(\frac{2(1-z_{2})(z_{2}-z_{3})}{z_{2}(1-z_{3})}-z_{1}+z_{2}\right)$
$\displaystyle-\frac{z_{3}(1-z_{1})+(1-z_{2})^{3}}{z_{1}z_{2}}+\epsilon(1-z_{2})\left(\frac{z_{1}^{2}+z_{2}^{2}}{z_{1}z_{2}}-\epsilon\right)\bigg{]}-\frac{2m_{Q}^{2}}{\tilde{s}_{12}}\frac{z_{1}^{2}+z_{2}^{2}}{z_{1}z_{2}}$
$\displaystyle+\frac{m_{Q}^{2}}{\tilde{s}_{13}\tilde{s}_{23}}\left[\frac{\tilde{s}_{12}}{2}\frac{z_{2}+4z_{3}-z_{2}z_{3}}{z_{2}}+\tilde{s}_{13}\frac{2z_{2}+4z_{3}-z_{1}z_{3}}{z_{1}}-\frac{\epsilon}{2}\tilde{s}_{123}(1-z_{3})\right]$
$\displaystyle-\frac{m_{Q}^{2}}{\tilde{s}_{12}\tilde{s}_{23}}\bigg{[}\tilde{s}_{12}\frac{3z_{1}z_{2}^{2}+z_{1}^{2}z_{2}-2z_{1}^{2}z_{3}+2z_{2}^{2}z_{3}+4z_{2}^{3}}{z_{1}z_{2}(1-z_{3})}$
$\displaystyle+2(\tilde{s}_{13}z_{2}-\tilde{s}_{23}z_{1})\frac{z_{1}^{2}+z_{1}z_{2}+z_{2}^{2}}{z_{1}z_{2}(1-z_{3})}\bigg{]}-\frac{m_{Q}^{4}}{\tilde{s}_{13}\tilde{s}_{23}}\Bigg{\\}}+(1\leftrightarrow
2)~{}.$ (51)

Finally, the splitting kernel with a gluon as a parent parton is decomposed as

$\displaystyle\hat{P}^{(0),\mu\nu}_{g_{1}Q_{2}\bar{Q}_{3}}=C_{F}T_{R}\hat{P}^{\rm{(0,ab)},\mu\nu}_{g_{1}Q_{2}\bar{Q}_{3}}+C_{A}T_{R}\hat{P}^{\rm{(0,nab)},\mu\nu}_{g_{1}Q_{2}\bar{Q}_{3}}~{},$
(52)

where

$\displaystyle\hat{P}^{\rm{(0,ab)},\mu\nu}_{g_{1}Q_{2}\bar{Q}_{3}}$
$\displaystyle=g^{\mu\nu}\bigg{[}2-\frac{2\tilde{s}_{123}\tilde{s}_{23}+(1-\epsilon)(\tilde{s}_{123}-\tilde{s}_{23})^{2}}{\tilde{s}_{12}\tilde{s}_{13}}+2m_{Q}^{2}\tilde{s}_{123}\frac{(\tilde{s}_{12}-\tilde{s}_{13})^{2}}{\tilde{s}_{12}^{2}\tilde{s}_{13}^{2}}$
$\displaystyle+4m_{Q}^{4}\left(\frac{1}{\tilde{s}_{12}^{2}}+\frac{1}{\tilde{s}_{13}^{2}}\right)\bigg{]}+\frac{4\tilde{s}_{123}}{\tilde{s}_{12}\tilde{s}_{13}}\left[\tilde{k}_{2}^{\mu}\tilde{k}_{3}^{\nu}+\tilde{k}_{3}^{\mu}\tilde{k}_{2}^{\nu}-(1-\epsilon)\tilde{k}_{1}^{\mu}\tilde{k}_{1}^{\nu}\right]$
$\displaystyle+8m_{Q}^{2}\left[\frac{\tilde{k}_{2}^{\mu}\tilde{k}_{2}^{\nu}}{\tilde{s}_{13}^{2}}+\frac{\tilde{k}_{3}^{\mu}\tilde{k}_{3}^{\nu}}{\tilde{s}_{12}^{2}}\right]-\frac{8m_{Q}^{2}}{\tilde{s}_{12}\tilde{s}_{13}}(1-\epsilon)\tilde{k}_{1}^{\mu}\tilde{k}_{1}^{\nu}~{},$
(53)

and

$\displaystyle\hat{P}^{\rm{(0,nab)},\mu\nu}_{g_{1}Q_{2}\bar{Q}_{3}}$
$\displaystyle=\Bigg{\\{}-\frac{2g^{\mu\nu}}{s_{23}^{2}}\bigg{[}\frac{z_{3}\tilde{s}_{12}(z_{2}\tilde{s}_{13}-z_{3}\tilde{s}_{12})}{(1-z_{1})^{2}}+\frac{z_{2}(1-2z_{2})(2-2z_{2}-5z_{3})}{2z_{1}(1-z_{1})^{2}}\tilde{s}_{23}^{2}$
$\displaystyle+\frac{\tilde{s}_{13}\tilde{s}_{23}}{z_{1}(1-z_{1})^{2}}\Big{(}z_{2}+z_{3}-2z_{2}^{2}-3z_{3}^{2}-3z_{2}z_{3}+2z_{2}^{3}+3z_{3}^{3}+5z_{2}^{2}z_{3}+6z_{2}z_{3}^{2}\Big{)}\bigg{]}$
$\displaystyle-\frac{2}{s_{23}^{2}}\bigg{[}\left(\frac{4z_{3}^{2}\tilde{s}_{123}}{z_{1}(1-z_{1})}-\frac{2(1-z_{2})\tilde{s}_{23}}{z_{1}}\right)\tilde{k}_{2}^{\mu}\tilde{k}_{2}^{\nu}-\left(\frac{4z_{2}z_{3}\tilde{s}_{123}}{z_{1}(1-z_{1})}-\frac{\tilde{s}_{23}(1-z_{1})}{z_{1}}\right)\tilde{k}_{2}^{\mu}\tilde{k}_{3}^{\nu}$
$\displaystyle+\epsilon\tilde{s}_{23}\tilde{k}_{1}^{\mu}\tilde{k}_{1}^{\nu}\bigg{]}-4g^{\mu\nu}\frac{m_{Q}^{2}}{s_{23}^{2}}\bigg{[}\frac{2z_{2}-4z_{2}^{2}-3z_{2}z_{3}}{(1-z_{1})^{2}}\tilde{s}_{23}+\tilde{s}_{12}\bigg{(}\frac{z_{3}(z_{2}-z_{3})}{(1-z_{1})^{2}}+\frac{z_{1}}{1-z_{1}}$
$\displaystyle+\frac{1-z_{1}}{z_{1}}\bigg{)}\bigg{]}+\frac{4m_{Q}^{2}}{s_{23}^{2}}\bigg{[}\left(2+\frac{z_{3}(z_{2}-z_{3})}{z_{1}(1-z_{1})}\right)2\tilde{k}_{2}^{\mu}\tilde{k}_{2}^{\nu}+\left(2-\frac{(z_{2}-z_{3})^{2}}{z_{1}(1-z_{1})}\right)\tilde{k}_{2}^{\mu}\tilde{k}_{3}^{\nu}-\epsilon\tilde{k}_{1}^{\mu}\tilde{k}_{1}^{\nu}\bigg{]}$
$\displaystyle+4g^{\mu\nu}\frac{m_{Q}^{4}}{s_{23}^{2}}\left[1-\frac{z_{1}}{1-z_{1}}-\frac{z_{2}z_{3}}{(1-z_{1})^{2}}\right]+\frac{g^{\mu\nu}}{2s_{23}}\bigg{[}\frac{2-3z_{1}}{1-z_{1}}\frac{\tilde{s}_{23}^{2}}{\tilde{s}_{13}}-\frac{z_{3}(1-2z_{1})}{z_{1}(1-z_{1})}\frac{\tilde{s}_{12}^{2}}{\tilde{s}_{13}}$
$\displaystyle+\frac{z_{2}\tilde{s}_{13}}{z_{1}(1-z_{1})}-\frac{1-2z_{2}}{z_{1}}\tilde{s}_{12}-\frac{1-3z_{2}-4z_{3}+2z_{2}^{2}+5z_{3}^{2}+7z_{2}z_{3}}{z_{1}(1-z_{1})}\frac{\tilde{s}_{12}\tilde{s}_{23}}{\tilde{s}_{13}}$
$\displaystyle-\frac{1-2z_{2}-z_{3}+3z_{2}^{2}+4z_{3}^{2}+7z_{2}z_{3}}{z_{1}(1-z_{1})}\tilde{s}_{23}\bigg{]}+\frac{1}{s_{23}}\bigg{[}\left(\frac{\tilde{s}_{123}}{\tilde{s}_{13}}-\frac{1-z_{3}}{z_{1}}\right)2\tilde{k}_{3}^{\mu}\tilde{k}_{3}^{\nu}$
$\displaystyle+\left(\frac{z_{2}-z_{3}-z_{2}^{2}+z_{3}^{2}+2z_{2}z_{3}}{z_{1}(1-z_{1})}\frac{\tilde{s}_{123}}{\tilde{s}_{13}}-\frac{z_{3}}{z_{1}}\frac{\tilde{s}_{23}}{\tilde{s}_{13}}+\frac{1-z_{1}}{z_{1}}\right)\left(\tilde{k}_{2}^{\mu}\tilde{k}_{3}^{\nu}+\tilde{k}_{3}^{\mu}\tilde{k}_{2}^{\nu}\right)$
$\displaystyle+\left(\frac{1-z_{2}}{z_{1}}\frac{\tilde{s}_{23}-\tilde{s}_{13}}{\tilde{s}_{13}}-\frac{2z_{3}^{2}}{z_{1}(1-z_{1})}\frac{\tilde{s}_{123}}{\tilde{s}_{13}}\right)2\tilde{k}_{2}^{\mu}\tilde{k}_{2}^{\nu}+\epsilon\bigg{(}\frac{\tilde{s}_{123}}{\tilde{s}_{13}}\left(\tilde{k}_{1}^{\mu}\tilde{k}_{3}^{\nu}+\tilde{k}_{3}^{\mu}\tilde{k}_{1}^{\nu}\right)$
$\displaystyle+\frac{\tilde{s}_{23}}{\tilde{s}_{13}}\left(\tilde{k}_{1}^{\mu}\tilde{k}_{2}^{\nu}+\tilde{k}_{2}^{\mu}\tilde{k}_{1}^{\nu}\right)+2\tilde{k}_{1}^{\mu}\tilde{k}_{1}^{\nu}\bigg{)}\bigg{]}+\frac{2m_{Q}^{2}}{s_{23}\tilde{s}_{23}}\bigg{[}\left(3+\frac{z_{3}(z_{2}-z_{3})}{z_{1}(1-z_{1})}\right)2\tilde{k}_{2}^{\mu}\tilde{k}_{2}^{\nu}+2\tilde{k}_{3}^{\mu}\tilde{k}_{3}^{\nu}$
$\displaystyle+\left(\frac{z_{3}(z_{2}-z_{3})}{z_{1}(1-z_{1})}+\frac{3z_{2}+z_{3}}{1-z_{1}}\right)\left(\tilde{k}_{2}^{\mu}\tilde{k}_{3}^{\nu}+\tilde{k}_{3}^{\mu}\tilde{k}_{2}^{\nu}\right)-2\epsilon\tilde{k}_{1}^{\mu}\tilde{k}_{1}^{\nu}\bigg{]}+2g^{\mu\nu}\frac{m_{Q}^{4}}{s_{23}\tilde{s}_{13}}\frac{2-3z_{1}}{1-z_{1}}$
$\displaystyle-g^{\mu\nu}\frac{m_{Q}^{2}}{s_{23}}\bigg{[}\frac{z_{1}^{2}+2z_{3}-z_{1}z_{2}-4z_{1}z_{3}}{z_{1}(1-z_{1})}\frac{\tilde{s}_{12}}{\tilde{s}_{13}}+\frac{z_{1}^{2}+2z_{2}^{2}+4z_{3}^{2}+z_{1}z_{3}+6z_{2}z_{3}}{z_{1}(1-z_{1})}$
$\displaystyle-\frac{2(2-3z_{1})}{1-z_{1}}\frac{\tilde{s}_{23}}{\tilde{s}_{13}}\bigg{]}+g^{\mu\nu}\bigg{[}\frac{\tilde{s}_{23}^{2}}{2\tilde{s}_{12}\tilde{s}_{13}}+\frac{z_{1}-2z_{2}}{4z_{1}}\frac{\tilde{s}_{23}}{\tilde{s}_{12}}-\frac{z_{2}}{4z_{1}}\frac{\tilde{s}_{13}}{\tilde{s}_{12}}-\frac{z_{3}}{4z_{1}}\frac{\tilde{s}_{12}}{\tilde{s}_{13}}+\frac{1-z_{1}}{4z_{1}}$
$\displaystyle+\frac{z_{1}-2z_{3}}{4z_{1}}\frac{\tilde{s}_{23}}{\tilde{s}_{13}}-\frac{\epsilon}{2}\bigg{]}+\left[\frac{1-z_{2}}{z_{1}}\frac{\tilde{s}_{23}}{\tilde{s}_{12}\tilde{s}_{13}}-\frac{1-z_{2}}{z_{1}\tilde{s}_{12}}-\frac{z_{3}}{z_{1}}\frac{\tilde{s}_{123}}{\tilde{s}_{12}\tilde{s}_{13}}\right]\tilde{k}_{2}^{\mu}\tilde{k}_{2}^{\nu}$
$\displaystyle+\left[\frac{\tilde{s}_{123}}{\tilde{s}_{12}\tilde{s}_{13}}-\frac{1-z_{3}}{z_{1}\tilde{s}_{12}}\right]\tilde{k}_{3}^{\mu}\tilde{k}_{3}^{\nu}+\left[\frac{1-2z_{3}}{2z_{1}}\frac{\tilde{s}_{23}}{\tilde{s}_{12}\tilde{s}_{13}}+\frac{1-2z_{2}}{2z_{1}}\frac{\tilde{s}_{123}}{\tilde{s}_{12}\tilde{s}_{13}}+\frac{z_{2}-z_{3}}{z_{1}\tilde{s}_{12}}\right]\tilde{k}_{2}^{\mu}\tilde{k}_{3}^{\nu}$
$\displaystyle+\left[\frac{1}{2z_{1}}\frac{\tilde{s}_{12}+\tilde{s}_{13}}{\tilde{s}_{12}\tilde{s}_{13}}-\frac{\tilde{s}_{123}}{\tilde{s}_{12}\tilde{s}_{13}}\right]\tilde{k}_{3}^{\mu}\tilde{k}_{2}^{\nu}+\epsilon\left[\frac{\tilde{k}_{3}^{\mu}\tilde{k}_{1}^{\nu}}{\tilde{s}_{13}}+\frac{\tilde{k}_{1}^{\mu}\tilde{k}_{2}^{\nu}}{\tilde{s}_{12}}-\frac{\tilde{s}_{23}}{\tilde{s}_{12}\tilde{s}_{13}}\tilde{k}_{1}^{\mu}\tilde{k}_{1}^{\nu}\right]$
$\displaystyle+g^{\mu\nu}m_{Q}^{2}\left[\frac{\tilde{s}_{23}}{\tilde{s}_{12}\tilde{s}_{13}}-\frac{z_{1}+2z_{2}}{2z_{1}\tilde{s}_{12}}-\frac{z_{1}+2z_{3}}{2z_{1}\tilde{s}_{13}}\right]+\frac{2m_{Q}^{2}}{\tilde{s}_{12}\tilde{s}_{13}}\bigg{[}\tilde{k}_{2}^{\mu}\tilde{k}_{3}^{\nu}-\tilde{k}_{3}^{\mu}\tilde{k}_{2}^{\nu}-\tilde{k}_{1}^{\mu}\tilde{k}_{3}^{\nu}$
$\displaystyle-\tilde{k}_{2}^{\mu}\tilde{k}_{1}^{\nu}-\epsilon\tilde{k}_{1}^{\mu}\tilde{k}_{1}^{\nu}\bigg{]}\Bigg{\\}}+(2\leftrightarrow
3).$ (54)

By following the description given in eq. (17), which is equivalent to
performing an azimuthal average, we obtain the unpolarized gluon-initiated
splitting kernel. It is given by

$\displaystyle\langle\hat{P}^{(0)}_{g_{1}Q_{2}\bar{Q}_{3}}\rangle=C_{F}T_{R}\,\langle\hat{P}^{\rm{(0,ab)}}_{g_{1}Q_{2}\bar{Q}_{3}}\rangle+C_{A}T_{R}\,\langle\hat{P}^{\rm{(0,nab)}}_{g_{1}Q_{2}\bar{Q}_{3}}\rangle~{},$
(55)

with

$\displaystyle\langle\hat{P}^{\rm{(0,ab)}}_{g_{1}Q_{2}\bar{Q}_{3}}\rangle$
$\displaystyle=$
$\displaystyle\Bigg{\\{}\left[\frac{s_{123}^{2}\left(z_{1}^{2}(1-\epsilon)+2(1-z_{2})z_{3}-\epsilon\right)}{\tilde{s}_{12}\tilde{s}_{13}(1-\epsilon)}\right.$
(56) $\displaystyle+$
$\displaystyle\left.\frac{2s_{123}((z_{1}+1)\epsilon+z_{2}-1)}{\tilde{s}_{12}(1-\epsilon)}+\frac{\tilde{s}_{13}(1-\epsilon)}{\tilde{s}_{12}}-\epsilon\right]$
$\displaystyle+$
$\displaystyle\left[\frac{m_{Q}^{2}}{\tilde{s}_{12}}\,\left(\frac{2s_{123}(2(1-z_{3})z_{3}+\epsilon-1)}{\tilde{s}_{12}}+\frac{2s_{123}(z_{1}+2z_{2}z_{3}+\epsilon)}{\tilde{s}_{13}}-4\right)\right.$
$\displaystyle-$
$\displaystyle\left.\,\frac{4\,m_{Q}^{4}}{\tilde{s}_{13}}\left(\frac{1}{\tilde{s}_{12}}+\frac{1}{\tilde{s}_{13}}\right)\right]\frac{1}{1-\epsilon}\Bigg{\\}}\,+(2\leftrightarrow
3)\,,$

and

$\displaystyle\langle\hat{P}^{\rm{(0,nab)}}_{g_{1}Q_{2}\bar{Q}_{3}}\rangle$
$\displaystyle=$
$\displaystyle\Bigg{\\{}\left[\frac{s_{123}^{2}z_{3}}{2s_{23}\tilde{s}_{13}}\left(\frac{(1-z_{1})^{3}-z_{1}^{3}}{z_{1}(1-z_{1})}-\frac{2z_{3}(-2z_{1}z_{2}-z_{3}+1)}{z_{1}(1-z_{1})(1-\epsilon)}\right)\right.$
(57) $\displaystyle-$
$\displaystyle\left.\frac{s_{123}^{2}}{2\tilde{s}_{12}\tilde{s}_{13}}\left(z_{1}^{2}-\frac{z_{1}+2z_{2}z_{3}}{1-\epsilon}+1\right)-\frac{(\tilde{t}_{23,1})^{2}}{4s_{23}^{2}}+\frac{\epsilon}{2}-\frac{1}{4}\right.$
$\displaystyle+$
$\displaystyle\left.\frac{s_{123}}{2s_{23}}\left(\frac{z_{1}^{3}+1}{z_{1}(1-z_{1})}+\frac{z_{1}(z_{3}-z_{2})^{2}-2(z_{1}+1)z_{2}z_{3}}{z_{1}(1-z_{1})(1-\epsilon)}\right)\right.$
$\displaystyle+$
$\displaystyle\left.\frac{s_{123}(1-z_{2})}{2\tilde{s}_{13}}\left(-\frac{2(1-z_{2})z_{2}}{z_{1}(1-z_{1})(1-\epsilon)}+\frac{1}{z_{1}(1-z_{1})}+1\right)\right]$
$\displaystyle+$
$\displaystyle\left[\frac{m_{Q}^{2}}{s^{2}_{23}}\,\left(\frac{2s_{123}^{2}(1-z_{2})}{\tilde{s}_{13}(1-z_{1})z_{1}}-\frac{s_{123}^{3}(z_{1}+2z_{2}z_{3}+\epsilon)}{\tilde{s}_{12}\tilde{s}_{13}}-\frac{2s_{123}^{2}z_{1}^{2}(1-2z_{2})}{\tilde{s}_{13}(1-z_{1})}\right.\right.$
$\displaystyle+$
$\displaystyle\left.\left.\frac{2s_{123}^{2}(4(1-z_{2})z_{2}+z_{2}+2\epsilon-2)}{\tilde{s}_{13}}-\frac{2s_{123}\tilde{s}_{12}z_{2}}{\tilde{s}_{13}(1-z_{1})}-\frac{2s_{123}\tilde{s}_{12}z_{3}}{\tilde{s}_{13}z_{1}}\right.\right.$
$\displaystyle-$
$\displaystyle\left.\left.s_{123}(z_{1}(1-4z_{2})+4(1-z_{2})z_{2}+2\epsilon+3)+\frac{s_{123}}{z_{1}}+4\tilde{s}_{12}\right.\right.$
$\displaystyle-$
$\displaystyle\left.\left.\frac{2s_{123}\tilde{s}_{12}\left(z_{2}-2z_{2}^{2}+\epsilon\right)}{\tilde{s}_{13}}\right)\,+\,\frac{2\,m_{Q}^{4}}{\tilde{s}_{12}\tilde{s}_{13}}\right]\frac{1}{1-\epsilon}\Bigg{\\}}\,+(2\leftrightarrow
3)\,.$

As an independent and highly non-trivial check on our results presented in
eqs. (5-57), in the massless limit, we confirm the expressions for the triple-
collinear splitting kernels given in ref. Catani:1999ss .

Finally, let us comment on the extension of the present results to other
theories, in particular to QED. All the expressions presented here were
computed within QCD, but we carefully stripped the color factors. Besides, in
$g\rightarrow gQ\bar{Q}$ and $Q\rightarrow ggQ$, we explicitly separated the
splitting kernels into Abelian and non-Abelian contributions. Thus, we can
compute the corresponding QED splitting kernels by replacing gluons by
photons, taking the Abelian part and replacing the color factors by electric
charges. For the four-quark splittings, the QCD results are already Abelian
(i.e. gluon self-interactions are absent), so we have

$\displaystyle\langle\hat{P}^{\rm{(0,QED)}}_{\bar{Q}_{1}^{\prime}Q_{2}^{\prime}Q_{3}}\rangle$
$\displaystyle=$ $\displaystyle
e_{Q}^{2}e_{Q^{\prime}}^{2}\left[\langle\hat{P}^{\rm{(0)}}_{\bar{Q}_{1}^{\prime}Q_{2}^{\prime}Q_{3}}\rangle\right]_{C_{A}\rightarrow
0,\,C_{F}\rightarrow 1,\,T_{R}\rightarrow 1}\,,$ (58)
$\displaystyle\langle\hat{P}^{\rm{(0,QED)}}_{\bar{Q}_{1}Q_{2}Q_{3}}\rangle$
$\displaystyle=$ $\displaystyle
e_{Q}^{4}\left[\langle\hat{P}^{\rm{(0)}}_{\bar{Q}_{1}Q_{2}Q_{3}}\rangle\right]_{C_{A}\rightarrow
0,\,C_{F}\rightarrow 1,\,T_{R}\rightarrow 1}\,,$ (59)

whilst for the remaining cases, we write

$\displaystyle\langle\hat{P}^{\rm{(0,QED)}}_{\gamma_{1}\gamma_{2}Q_{3}}\rangle$
$\displaystyle=$ $\displaystyle
e_{Q}^{4}\langle\hat{P}^{\rm{(0,ab)}}_{g_{1}g_{2}Q_{3}}\rangle\,,$ (60)
$\displaystyle\hat{P}^{\rm{(0,QED)}}_{\gamma_{1}Q_{2}\bar{Q}_{3}}$
$\displaystyle=$ $\displaystyle
e_{Q}^{4}\hat{P}^{\rm{(0,ab)}}_{g_{1}Q_{2}\bar{Q}_{3}}\,.$ (61)

These expressions are in agreement with the relations presented in eq. (71)
from ref. Catani:1999ss  for the massless case. The procedure described before
is known in general as Abelianization deFlorian:2015ujt ; deFlorian:2016gvk ;
deFlorian:2018wcj ; Ajjath:2019ixh ; Ajjath:2019vmf , and it is independent on
the mass of the particles.

##  6 Conclusions and outlook

In this article, we have presented the first calculation of triple-collinear
splitting amplitudes and (un)polarized splitting kernels with massive partons
at tree-level in QCD. We have carefully compared the kinematics of the
collinear and quasi-collinear limits, and have provided very compact explicit
expressions for all the possible processes involving QCD partons.

Although collinear singularities are strictly absent for splitting processes
with massive particles, mass-dependent logarithmically enhanced terms emerge
in scattering amplitudes from the quasi-collinear region. Since these
contributions depend on ratios of masses to the typical energy scale of the
process, they could lead to large corrections in the high-energy limit. So, in
order to improve theoretical predictions it is necessary to use subtraction
counter-terms involving massive splitting kernels. At NLO, the massive dipole
formalism Catani:2002hc  required double-collinear splittings to build the
corresponding counter-terms. The massive triple-collinear splittings provided
in this article allow to extend the collinear subtraction formalism to reach
NNLO accuracy, taking into account parton mass effects. This might play a
crucial role in collider phenomenology involving heavy quarks. For instance,
in forward charm-quark production Maciula:2022lzk , there are no complete
calculations that keeps track of the mass of the charm quark, which is
expected to have an important phenomenological impact.

Also, having access to the full set of triple-collinear splitting kernels with
massive partons paves the road for more precise parton shower generators. Up
to now, their improvement relied mostly on the inclusion of higher-logarithmic
corrections, since mass effects were always consider as sub-dominant. However,
as we are reaching next-to-next-to-leading logarithmic accuracy or even
higher-orders, mass effects start to compete and can no longer be neglected.

Finally, the knowledge of triple-collinear splitting processes with massive
partons could play an important role in the consistent treatment of mass
effects in parton densities evolution. For this purpose, higher-order
corrections to the Altarelli-Parisi kernels of splitting functions, keeping
the exact mass dependence, and generalising the DGLAP equations to the quasi-
collinear limit are required. The results provided in this work constitute a
very first step towards this direction, and might be relevant in sight of the
current precision high-energy collider phenomenology program.

## Acknowledgements

We would like to thank the Galileo Galilei Institute (Firenze), where this
work was completed, and S. Catani, for hospitality and discussing during the
Workshop Theory Challenges in the Precision Era of the Large Hadron Collider.
We would also like to thank A. Szczurek and M. Wiesemann for their fruitful
comments when presenting this work at EPS-HEP 2023. This work is supported by
the Spanish Government (Agencia Estatal de Investigación MCIN
/AEI/10.13039/501100011033) Grants No. PID2020-114473GB-I00,
PID2022-141910NB-I00 and Generalitat Valenciana Grant No. PROMETEO/2021/071.
The work of PKD is supported by European Commission MSCA Action COLLINEAR-
FRACTURE, Grant Agreement No. 101108573. The work of GS is partially supported
by EU Horizon 2020 research and innovation programme STRONG-2020 project under
Grant Agreement No. 824093 and H2020-MSCA-COFUND-2020 USAL4EXCELLENCE-
PROOPI-391 project under Grant Agreement No 101034371.

## References

  * (1) J. C. Collins, D. E. Soper and G. F. Sterman, _Factorization of Hard Processes in QCD_ , _Adv. Ser. Direct. High Energy Phys._ 5 (1989) 1–91, [hep-ph/0409313].
  * (2) S. Catani, D. de Florian and G. Rodrigo, _Space-like (versus time-like) collinear limits in QCD: Is factorization violated?_ , _JHEP_ 07 (2012) 026, [1112.4405].
  * (3) J. R. Forshaw, M. H. Seymour and A. Siodmok, _On the Breaking of Collinear Factorization in QCD_ , _JHEP_ 11 (2012) 066, [1206.6363].
  * (4) C. G. Bollini and J. J. Giambiagi, _Dimensional Renormalization: The Number of Dimensions as a Regularizing Parameter_ , _Nuovo Cim. B_ 12 (1972) 20–26.
  * (5) G. ’t Hooft and M. J. G. Veltman, _Regularization and Renormalization of Gauge Fields_ , _Nucl. Phys. B_ 44 (1972) 189–213.
  * (6) G. M. Cicuta and E. Montaldi, _Analytic renormalization via continuous space dimension_ , _Lett. Nuovo Cim._ 4 (1972) 329–332.
  * (7) S. Catani and M. H. Seymour, _A General algorithm for calculating jet cross-sections in NLO QCD_ , _Nucl. Phys. B_ 485 (1997) 291–419, [hep-ph/9605323].
  * (8) S. Frixione, Z. Kunszt and A. Signer, _Three jet cross-sections to next-to-leading order_ , _Nucl. Phys. B_ 467 (1996) 399–442, [hep-ph/9512328].
  * (9) W. J. Torres Bobadilla et al., _May the four be with you: Novel IR-subtraction methods to tackle NNLO calculations_ , _Eur. Phys. J. C_ 81 (2021) 250, [2012.02567].
  * (10) F. Caola, W. Chen, C. Duhr, X. Liu, B. Mistlberger, F. Petriello et al., _The Path forward to N 3LO_, in _Snowmass 2021_ , 3, 2022.  2203.06730.
  * (11) N. Agarwal, L. Magnea, C. Signorile-Signorile and A. Tripathi, _The infrared structure of perturbative gauge theories_ , _Phys. Rept._ 994 (2023) 1–120, [2112.07099].
  * (12) S. Catani and P. K. Dhani, _Collinear functions for QCD resummations_ , _JHEP_ 03 (2023) 200, [2208.05840].
  * (13) S. Höche, _Introduction to parton-shower event generators_ , in _Theoretical Advanced Study Institute in Elementary Particle Physics: Journeys Through the Precision Frontier: Amplitudes for Colliders_ , pp. 235–295, 2015.  1411.4085.  DOI.
  * (14) Z. Bern, G. Chalmers, L. J. Dixon and D. A. Kosower, _One loop N gluon amplitudes with maximal helicity violation via collinear limits_ , _Phys. Rev. Lett._ 72 (1994) 2134–2137, [hep-ph/9312333].
  * (15) Z. Bern, L. J. Dixon, D. C. Dunbar and D. A. Kosower, _One loop n point gauge theory amplitudes, unitarity and collinear limits_ , _Nucl. Phys. B_ 425 (1994) 217–260, [hep-ph/9403226].
  * (16) Z. Bern and G. Chalmers, _Factorization in one loop gauge theory_ , _Nucl. Phys. B_ 447 (1995) 465–518, [hep-ph/9503236].
  * (17) D. A. Kosower, _All order collinear behavior in gauge theories_ , _Nucl. Phys. B_ 552 (1999) 319–336, [hep-ph/9901201].
  * (18) Z. Bern, V. Del Duca, W. B. Kilgore and C. R. Schmidt, _The infrared behavior of one loop QCD amplitudes at next-to-next-to leading order_ , _Phys. Rev. D_ 60 (1999) 116001, [hep-ph/9903516].
  * (19) S. Catani and M. Grazzini, _Infrared factorization of tree level QCD amplitudes at the next-to-next-to-leading order and beyond_ , _Nucl. Phys. B_ 570 (2000) 287–325, [hep-ph/9908523].
  * (20) S. Catani, D. de Florian and G. Rodrigo, _The Triple collinear limit of one loop QCD amplitudes_ , _Phys. Lett. B_ 586 (2004) 323–331, [hep-ph/0312067].
  * (21) G. Altarelli and G. Parisi, _Asymptotic Freedom in Parton Language_ , _Nucl. Phys. B_ 126 (1977) 298–318.
  * (22) J. M. Campbell and E. W. N. Glover, _Double unresolved approximations to multiparton scattering amplitudes_ , _Nucl. Phys. B_ 527 (1998) 264–288, [hep-ph/9710255].
  * (23) S. Catani and M. Grazzini, _Collinear factorization and splitting functions for next-to-next-to-leading order QCD calculations_ , _Phys. Lett. B_ 446 (1999) 143–152, [hep-ph/9810389].
  * (24) Z. Bern, V. Del Duca and C. R. Schmidt, _The Infrared behavior of one loop gluon amplitudes at next-to-next-to-leading order_ , _Phys. Lett. B_ 445 (1998) 168–177, [hep-ph/9810409].
  * (25) D. A. Kosower and P. Uwer, _One loop splitting amplitudes in gauge theory_ , _Nucl. Phys. B_ 563 (1999) 477–505, [hep-ph/9903515].
  * (26) G. F. R. Sborlini, D. de Florian and G. Rodrigo, _Double collinear splitting amplitudes at next-to-leading order_ , _JHEP_ 01 (2014) 018, [1310.6841].
  * (27) V. Del Duca, A. Frizzo and F. Maltoni, _Factorization of tree QCD amplitudes in the high-energy limit and in the collinear limit_ , _Nucl. Phys. B_ 568 (2000) 211–262, [hep-ph/9909464].
  * (28) T. G. Birthwright, E. W. N. Glover, V. V. Khoze and P. Marquard, _Multi-gluon collinear limits from MHV diagrams_ , _JHEP_ 05 (2005) 013, [hep-ph/0503063].
  * (29) T. G. Birthwright, E. W. N. Glover, V. V. Khoze and P. Marquard, _Collinear limits in QCD from MHV rules_ , _JHEP_ 07 (2005) 068, [hep-ph/0505219].
  * (30) M. Delto and K. Melnikov, _Integrated triple-collinear counter-terms for the nested soft-collinear subtraction scheme_ , _JHEP_ 05 (2019) 148, [1901.05213].
  * (31) V. Del Duca, C. Duhr, R. Haindl, A. Lazopoulos and M. Michel, _Tree-level splitting amplitudes for a gluon into four collinear partons_ , _JHEP_ 10 (2020) 093, [2007.05345].
  * (32) G. F. R. Sborlini, D. de Florian and G. Rodrigo, _Triple collinear splitting functions at NLO for scattering processes with photons_ , _JHEP_ 10 (2014) 161, [1408.4821].
  * (33) G. F. R. Sborlini, D. de Florian and G. Rodrigo, _Polarized triple-collinear splitting functions at NLO for processes with photons_ , _JHEP_ 03 (2015) 021, [1409.6137].
  * (34) S. Badger, F. Buciuni and T. Peraro, _One-loop triple collinear splitting amplitudes in QCD_ , _JHEP_ 09 (2015) 188, [1507.05070].
  * (35) M. Czakon and S. Sapeta, _Complete collection of one-loop triple-collinear splitting operators for dimensionally-regulated QCD_ , _JHEP_ 07 (2022) 052, [2204.11801].
  * (36) Z. Bern, L. J. Dixon and D. A. Kosower, _Two-loop g — $>$ gg splitting amplitudes in QCD_, _JHEP_ 08 (2004) 012, [hep-ph/0404293].
  * (37) S. D. Badger and E. W. N. Glover, _Two loop splitting functions in QCD_ , _JHEP_ 07 (2004) 040, [hep-ph/0405236].
  * (38) L. J. Dixon, E. Herrmann, K. Yan and H. X. Zhu, _Soft gluon emission at two loops in full color_ , _JHEP_ 05 (2020) 135, [1912.09370].
  * (39) A. Vogt, S. Moch and J. A. M. Vermaseren, _The Three-loop splitting functions in QCD: The Singlet case_ , _Nucl. Phys. B_ 691 (2004) 129–181, [hep-ph/0404111].
  * (40) A. Vogt, S. Moch and J. Vermaseren, _Photon-parton splitting functions at the next-to-next-to-leading order of QCD_ , _Acta Phys. Polon. B_ 37 (2006) 683–688, [hep-ph/0511112].
  * (41) S. Moch, J. A. M. Vermaseren and A. Vogt, _The Three loop splitting functions in QCD: The Nonsinglet case_ , _Nucl. Phys. B_ 688 (2004) 101–134, [hep-ph/0403192].
  * (42) S. Höche and S. Prestel, _Triple collinear emissions in parton showers_ , _Phys. Rev. D_ 96 (2017) 074017, [1705.00742].
  * (43) A. M. Snigirev, _Triple parton scattering in collinear approximation of perturbative QCD_ , _Phys. Rev. D_ 94 (2016) 034026, [1603.08187].
  * (44) Y. L. Dokshitzer, V. A. Khoze and S. I. Troian, _On specific QCD properties of heavy quark fragmentation (’dead cone’)_ , _J. Phys. G_ 17 (1991) 1602–1604.
  * (45) S. Catani, S. Dittmaier and Z. Trocsanyi, _One loop singular behavior of QCD and SUSY QCD amplitudes with massive partons_ , _Phys. Lett. B_ 500 (2001) 149–160, [hep-ph/0011222].
  * (46) S. Keller and E. Laenen, _Next-to-leading order cross-sections for tagged reactions_ , _Phys. Rev. D_ 59 (1999) 114004, [hep-ph/9812415].
  * (47) S. Dittmaier, _A General approach to photon radiation off fermions_ , _Nucl. Phys. B_ 565 (2000) 69–122, [hep-ph/9904440].
  * (48) M. Roth, _Precise predictions for four fermion production in electron positron annihilation_ , other thesis, 1999.
  * (49) L. Phaf and S. Weinzierl, _Dipole formalism with heavy fermions_ , _JHEP_ 04 (2001) 006, [hep-ph/0102207].
  * (50) S. Catani, S. Dittmaier, M. H. Seymour and Z. Trocsanyi, _The Dipole formalism for next-to-leading order QCD calculations with massive partons_ , _Nucl. Phys. B_ 627 (2002) 189–265, [hep-ph/0201036].
  * (51) P. Nogueira, _Automatic Feynman graph generation_ , _J. Comput. Phys._ 105 (1993) 279–289.
  * (52) J. Kuipers, T. Ueda, J. A. M. Vermaseren and J. Vollinga, _FORM version 4.0_ , _Comput. Phys. Commun._ 184 (2013) 1453–1467, [1203.6543].
  * (53) V. Shtabovenko, R. Mertig and F. Orellana, _New Developments in FeynCalc 9.0_ , _Comput. Phys. Commun._ 207 (2016) 432–444, [1601.01167].
  * (54) V. Shtabovenko, R. Mertig and F. Orellana, _FeynCalc 9.3: New features and improvements_ , _Comput. Phys. Commun._ 256 (2020) 107478, [2001.04407].
  * (55) T. Hahn, _Generating Feynman diagrams and amplitudes with FeynArts 3_ , _Comput. Phys. Commun._ 140 (2001) 418–431, [hep-ph/0012260].
  * (56) D. de Florian, G. F. R. Sborlini and G. Rodrigo, _QED corrections to the Altarelli–Parisi splitting functions_ , _Eur. Phys. J. C_ 76 (2016) 282, [1512.00612].
  * (57) D. de Florian, G. F. R. Sborlini and G. Rodrigo, _Two-loop QED corrections to the Altarelli-Parisi splitting functions_ , _JHEP_ 10 (2016) 056, [1606.02887].
  * (58) D. de Florian, M. Der and I. Fabre, _QCD $\oplus$QED NNLO corrections to Drell Yan production_, _Phys. Rev. D_ 98 (2018) 094008, [1805.12214].
  * (59) A. H. Ajjath, P. Banerjee, A. Chakraborty, P. K. Dhani, P. Mukherjee, N. Rana et al., _NNLO QCD $\oplus$QED corrections to Higgs production in bottom quark annihilation_, _Phys. Rev. D_ 100 (2019) 114016, [1906.09028].
  * (60) A. H. Ajjath, P. Mukherjee and V. Ravindran, _Infrared structure of $SU(N)\times U(1)$ gauge theory to three loops_, _JHEP_ 08 (2020) 156, [1912.13386].
  * (61) R. Maciula and A. Szczurek, _Far-forward production of charm mesons and neutrinos at forward physics facilities at the LHC and the intrinsic charm in the proton_ , _Phys. Rev. D_ 107 (2023) 034002, [2210.08890].


<EMAIL_ADDRESS>[url]https://www.asc.tuwien.ac.at/arnold/
<EMAIL_ADDRESS>
aff1]Inst. f. Analysis u. Scientific Computing, Technische Universität Wien,
Wiedner Hauptstr. 8, A-1040 Wien, Austria. [cor1]Corresponding author:
<EMAIL_ADDRESS>
# Hypocoercivity for Linear ODEs and Strong Stability for Runge–Kutta Methods

Franz Achleitner  Anton Arnold  Ansgar Jüngel  [

###### Abstract

In this note, we connect two different topics from linear algebra and
numerical analysis: hypocoercivity of semi-dissipative matrices and strong
stability for explicit Runge–Kutta schemes. Linear autonomous ODE systems with
a non-coercive matrix are called hypocoercive if they still exhibit uniform
exponential decay towards the steady state. Strong stability is a property of
time-integration schemes for ODEs that preserve the temporal monotonicity of
the discrete solutions. It is proved that explicit Runge–Kutta schemes are
strongly stable with respect to semi-dissipative, asymptotically stable
matrices if the hypocoercivity index is sufficiently small compared to the
order of the scheme. Otherwise, the Runge–Kutta schemes are in general not
strongly stable. As a corollary, explicit Runge–Kutta schemes of order $p\in
4{\mathbb{N}}$ with $s=p$ stages turn out to be _not_ strongly stable. This
result was proved in [4], filling a gap left open in [7]. Here, we present an
alternative, direct proof.

##  1 Hypocoercive ODEs

For linear autonomous ordinary differential equations (ODEs),

$\frac{\,{\rm d}u}{\,\,{\rm d}t}=\mathbf{L}u,\quad t>0,\quad
u(0)=u^{0}\in{\mathbb{C}}^{n},$ (1)

with Lyapunov stable matrices $\mathbf{L}\in{\mathbb{C}}^{n\times n}$ (i.e.,
all eigenvalues of $\mathbf{L}$ have nonpositive real part and the purely
imaginary eigenvalues are non-defective), we are concerned with characterizing
their short-time decay behavior. To this end, we review first hypocoercivity
properties of such systems [1, 9]:

###### Definition 1

  1. (a)

The matrix $\mathbf{L}\in{\mathbb{C}}^{n}$ is called _dissipative_ (resp.
_semi-dissipative_) if its Hermitian part,
$\mathbf{L}_{H}:=(\mathbf{L}+\mathbf{L}^{*})/2$ is negative definite (resp.
negative semi-definite).

  2. (b)

$-\mathbf{L}$ is called _hypocoercive_ (or _positive stable_) if there are
constants $\lambda>0$ and $c\geq 1$ such that the matrix exponential satisfies

$\|e^{\mathbf{L}t}\|_{2}\leq ce^{-\lambda t},\quad t\geq 0.$

  3. (c)

Let $\mathbf{L}\in{\mathbb{C}}^{n}$ be semi-dissipative. Its _hypocoercivity
index (HC-index)_ ${m_{HC}}$ is defined as the smallest integer
$m\in{\mathbb{N}}_{0}$ (if it exists) such that

$T_{m}:=\sum_{j=0}^{m}\mathbf{L}_{S}^{j}\mathbf{L}_{H}(\mathbf{L}_{S}^{*})^{j}<0,$

where $\mathbf{L}_{S}:=(\mathbf{L}-\mathbf{L}^{*})/2$ denotes the skew-
Hermitian part of $\mathbf{L}$.

The HC-index of $\mathbf{L}$ characterizes the structural complexity of the
interplay between $\mathbf{L}_{H}$ and $\mathbf{L}_{S}$, and it is bounded by

$0\leq\frac{n-\rm rank\,\mathbf{L}_{H}}{\rm
rank\,\mathbf{L}_{H}}\leq{m_{HC}}(\mathbf{L})\leq n-{\rm
rank\,\mathbf{L}_{H}}\leq n-1;$

see [1, 3, 5]. Hence the matrix $\mathbf{L}$ is dissipative if and only if
${m_{HC}}(\mathbf{L})=0$. For later use, we consider the example

$\mathbf{L}=\left(\begin{array}[]{rrrrr}0&-1&&&\\\ 1&\ddots&\ddots&&\\\
&\ddots&\ddots&-1&\\\ &&1&0&-1\\\
&&&1&-1\end{array}\right)\in{\mathbb{R}}^{N\times N}.$ (2)

Using the Kalman rank condition [1, Proposition 1], one easily verifies that
$\mathbf{L}$ has (maximal) HC-index $N-1$. Moreover, the matrix $\mathbf{L}$
is asymptotically stable, i.e., it has only eigenvalues with negative real
part.

We recall that the HC-index gives a precise characterization of the short-time
behavior of the solutions to (1):

######  Proposition 1 ([2])

Let the matrix $\mathbf{L}\in{\mathbb{C}}^{n\times n}$ be semi-dissipative.
Then $\mathbf{L}$ is asymptotically stable (with HC-index
${m_{HC}}\in{\mathbb{N}}_{0}$) if and only if

$\|e^{t\mathbf{L}}\|_{2}=1-ct^{a}+\mathcal{O}(t^{a+1})\quad\mbox{ for
}t\in[0,\varepsilon),$ (3)

for some $a,c,\varepsilon>0$. In this case, necessarily $a=2m_{HC}+1$.

The sharp multiplicative factor $c$ in (3) has been determined explicitly in
[2, Theorem 2.7(b)].

##  2 Strong stability for Runge–Kutta methods

It is known that a matrix $\mathbf{L}\in{\mathbb{C}}^{n\times n}$ is _Lyapunov
stable_ if and only if there exists a positive definite Hermitian matrix
$\mathbf{P}$ such that

$\mathbf{L}^{*}\mathbf{P}+\mathbf{P}\mathbf{L}\leq 0.$ (4)

In this case, the solution $u(t)$ of (1) is nonincreasing in the norm
$\|\cdot\|_{\mathbf{P}}:=\sqrt{\langle\cdot,\mathbf{P}\cdot\rangle}$ :

$\frac{\,{\rm d}}{\,\,{\rm d}t}\|u(t)\|_{\mathbf{P}}^{2}=\langle
u,(\mathbf{L}^{*}\mathbf{P}+\mathbf{P}\mathbf{L})u\rangle\leq 0.$ (5)

It is often desirable that a numerical scheme for (1) reproduces this decay
behavior on the discrete level.

Following [4], we consider here explicit Runge–Kutta schemes, where $u^{k}$ is
an approximation for $u(k\tau),\,k\in{\mathbb{N}}_{0}$, and $\tau$ is the
uniform time step:

$u^{k}=u^{k-1}+\tau\sum_{i=1}^{s}b_{i}K_{i}^{k},\quad
K_{i}^{k}=\mathbf{L}\bigg{(}u^{k-1}+\tau\sum_{j=1}^{i-1}a_{ij}K_{j}^{k}\bigg{)},\quad
i=1,\ldots,s,$ (6)

where $b_{i}\in{\mathbb{C}}$ are the weights, $a_{ij}\in{\mathbb{C}}$ are the
coefficients of the Runge–Kutta matrix, and $s\in{\mathbb{N}}$ is the number
of stages. This scheme can be rewritten in compact form as
$u^{k}=R(\tau\mathbf{L})u^{k-1}$, using its _stability function_ $R(z)$. For
the scheme (6) to have order $p$, one needs at least $s\geq p$ stages. In this
case, its stability function takes the form

$R(z)=\sum_{j=0}^{p}\frac{z^{j}}{j!}+\sum_{j=p+1}^{s}c_{j}\frac{z^{j}}{j!},\quad
z\in{\mathbb{C}},\quad c_{p+1}\neq 1,$ (7)

with some constants $c_{p+1},\ldots,c_{s}\in{\mathbb{C}}$.

For the discrete analog of the monotonicity estimate (5), we shall use the
following notions:

###### Definition 2

  * (a)

The Runge–Kutta scheme (6) is strongly stable if for all matrix dimensions
$n\in{\mathbb{N}}$, for all Lyapunov stable matrices
$\mathbf{L}\in{\mathbb{C}}^{n\times n}$, and for all Hermitian matrices
$\mathbf{P}>0$ such that (4) holds, the numerical solution to (1) satisfies
$\|u^{1}\|_{\mathbf{P}}\leq\|u^{0}\|_{\mathbf{P}}$ for all initial data
$u^{0}\in{\mathbb{C}}^{n}$ and sufficiently small time steps.

  * (b)

The Runge–Kutta scheme (6) is strongly stable w.r.t. a subset
$\mathcal{L}_{0}$ of Lyapunov stable matrices (of any dimension $n$), if the
condition from (a) holds for all $\mathbf{L}\in\mathcal{L}_{0}$.

Strong stability of explicit Runge–Kutta schemes was studied in, e.g., [8, §4]
and [7]. We shall focus on explicit Runge–Kutta schemes with $s=p$ stages.
Their stability function is given by the first sum in (7), and their stability
behavior was analyzed in [7]:

$\mbox{They are }\quad\left\\{\begin{array}[]{ll}\mbox{strongly
stable}&\quad\mbox{if }p\in 4{\mathbb{N}}_{0}+3\ ,\\\ \mbox{not strongly
stable}&\quad\mbox{if }p\in 4{\mathbb{N}}_{0}+1\mbox{ or }p\in
4{\mathbb{N}}_{0}+2\ ,\\\ \end{array}\right.$

but the following instability result for the case $p\in 4{\mathbb{N}}$ was
only found and proved recently in [4]:

######  Theorem 2 ([4])

Explicit Runge–Kutta schemes of order $p\in 4{\mathbb{N}}$ with $s=p$ stages
are _not_ strongly stable.

While the proof of this result in [4] was based on the quite technical result
in Proposition 1, we shall give here an independent direct proof. In order to
motivate our subsequent proof, we first cite another result from [4]. To this
end, we define the following subset of asymptotically stable (and thus
Lyapunov stable) matrices:

$\mathcal{L}_{AS}^{m}:=\\{\mathbf{L}\mbox{ is semi-dissipative and
asymptotically stable}:{m_{HC}}(\mathbf{L})\leq m\\},\quad
m\in{\mathbb{N}}_{0}.$

######  Proposition 3 ([4])

All explicit Runge–Kutta schemes of order $p\in{\mathbb{N}}$ (with $s\geq p$
stages) are strongly stable w.r.t. $\mathcal{L}_{AS}^{m}$, if
$m\in{\mathbb{N}}_{0}$ satisfies $2m+1\leq p$.

_Proof (of Theorem 2)._ For each fixed $p\in 4{\mathbb{N}}$, we shall analyze
an asymptotically stable matrix $\mathbf{L}$ as a counterexample. Due to
Proposition 3, the HC-index of such $\mathbf{L}$ must be “large enough”. More
precisely, we choose $\mathbf{L}$ of the form (2) with $N:=1+p/2$. It
satisfies $\mathbf{L}\in\mathcal{L}_{AS}^{m}$ with $m=p/2$, which violates the
index condition in Proposition 3, and may hence serve as a counterexample.

When choosing $\mathbf{P}=\mathbf{I}$, inequality (4) is satisfied, and it
remains to show that

$\|R(\tau\mathbf{L})\|_{2}:=\sup_{\|u_{0}\|=1}\|R(\tau\mathbf{L})u^{0}\|\leq
1,\quad\mbox{where}\ R(z)=\sum_{j=0}^{p}\frac{z^{j}}{j!},$ (8)

does _not_ hold on any interval $\tau\in[0,\varepsilon)$, with $\varepsilon>0$
arbitrarily small, see Definition 2(a). Equivalently, we shall show that the
matrix function
$\mathbf{M}(\tau):=\mathbf{I}-R(\tau\mathbf{L})^{*}R(\tau\mathbf{L})$ has a
negative determinant on $(0,\varepsilon)$ for sufficiently small
$\varepsilon>0$. Hence, $R(\tau\mathbf{L})^{*}R(\tau\mathbf{L})$ has at least
one eigenvalue larger than one, and $\|R(\tau\mathbf{L})\|_{2}>1$ follows. We
shall obtain the inequality $\det\mathbf{M}(\tau)<0$, for $\tau$ small enough,
from the subsequent lemma, thus closing this proof. $\hfill\square$

###### Lemma 4

Let $p\in 4{\mathbb{N}}$, and $\mathbf{L}$ of the form (2) with $N:=1+p/2$.
Then, $\mathbf{M}(\tau):=\mathbf{I}-R(\tau\mathbf{L})^{*}R(\tau\mathbf{L})$
satisfies

$\det\mathbf{M}(\tau)=c\tau^{N^{2}}+\mathcal{O}(\tau^{N^{2}+1})\quad\mbox{as
}\tau\to 0\,,$ (9)

with some $c<0$ and $R(z)$ defined in (8).

_Proof._ We only give a sketch of the proof here, the full details are given
in the Appendix.

First, we consider the Runge–Kutta method with $p=s=4$, and hence $N=3$. In
this case, we compute the matrix $\mathbf{M}(\tau)$ explicitly:

$\mathbf{M}(\tau)=\left(\begin{array}[]{rrr}\tau^{5}/12&\tau^{4}/4&\tau^{3}/3\\\
\tau^{4}/4&2\tau^{3}/3&\tau^{2}\\\
\tau^{3}/3&\tau^{2}&2\tau\end{array}\right)+\left(\begin{array}[]{rrr}\mathcal{O}(\tau^{6})&\mathcal{O}(\tau^{5})&\mathcal{O}(\tau^{4})\\\
\mathcal{O}(\tau^{5})&\mathcal{O}(\tau^{4})&\mathcal{O}(\tau^{3})\\\
\mathcal{O}(\tau^{4})&\mathcal{O}(\tau^{3})&\mathcal{O}(\tau^{2})\\\
\end{array}\right)\quad\mbox{as }\tau\to 0\,.$

The determinant of the first matrix yields the leading order coefficient:

$\det\mathbf{M}(\tau)=\det\left(\begin{array}[]{rrr}\tau^{5}/12&\tau^{4}/4&\tau^{3}/3\\\
\tau^{4}/4&2\tau^{3}/3&\tau^{2}\\\
\tau^{3}/3&\tau^{2}&2\tau\end{array}\right)+\mathcal{O}(\tau^{10})=-\frac{\tau^{9}}{216}+\mathcal{O}(\tau^{10})\quad\mbox{as
}\tau\to 0.$

Consequently, $\det\mathbf{M}(\tau)<0$ for sufficiently small $\tau>0$, which
disproves (8) for $p=4$.

For the general case $p\geq 8$, we insert the stability matrix (7),

$\mathbf{M}(\tau)=\mathbf{I}-\bigg{(}\sum_{j=0}^{p}\frac{\tau^{j}}{j!}(\mathbf{L}^{*})^{j}\bigg{)}\,\bigg{(}\sum_{\ell=0}^{p}\frac{\tau^{\ell}}{\ell!}\mathbf{L}^{\ell}\bigg{)},$

expand the matrix coefficients $M_{ij}$ of $\mathbf{M}(\tau)$ in powers of
$\tau$ and use properties of the matrices $\mathbf{L}$ and $\mathbf{L}_{H}$,
leading to

$M_{ij}=\bar{m}_{ij}\tau^{p+3-(i+j)}-\bar{n}_{ij}\tau^{p+1}+\mathcal{O}(\tau^{p+4-(i+j)}),$

where $\bar{m}_{ij}$ and $\bar{n}_{ij}$ are numbers depending on $p$ and
$\mathbf{L}_{H}$. The determinant of $\mathbf{M}(\tau)$ can be computed by
using the Leibniz formula for determinants and definition $N=p/2+1$:

$\det\mathbf{M}(\tau)=\det(\bar{m}_{ij}-\bar{n}_{ij})_{1\leq i,j\leq
N}\tau^{N^{2}}+\mathcal{O}(\tau^{N^{2}+1}).$

The remaining determinant can be calculated by taking into account the
explicit formula of the Hankel determinant:

$c:=\det(\bar{m}_{ij}-\bar{n}_{ij})_{1\leq i,j\leq
N}=2^{N}\bigg{(}\prod_{i=1}^{N}\frac{1}{(p/2-i+1)!}\bigg{)}^{2}\frac{\prod_{1\leq
i<j\leq
N}(i-j)^{2}}{\prod_{i,j=1}^{N}(i+(j-1))}\bigg{(}1-\binom{p}{p/2}\bigg{)}.$

We deduce from $\binom{p}{p/2}>1$ that $c<0$ as claimed in Lemma 4. In
particular, for sufficiently small $\tau>0$, $\det\mathbf{M}(\tau)$ is
negative, which also finishes the proof of Theorem 2. $\hfill\square$

Finally, we note that $\|R(\tau\mathbf{L})\|_{2}>1$, the condition to violate
strong stability, can also be tested numerically. But only for $p=4$ the
computation can be carried out in the standard double precision of Matlab: It
shows that $\varepsilon=0.304$ can be used, yielding
$\max_{[0,\varepsilon]}\|R(\tau\mathbf{L})\|_{2}-1\approx 1.3$E-6. For larger
values of $p$, this computation is numerically so sensitive that we used
octuple precision: For $p=8$ one can use $\varepsilon=0.027$, yielding
$\max_{[0,\varepsilon]}\|R(\tau\mathbf{L})\|_{2}-1\approx 9.0$E-22, and for
$p=12$, $\varepsilon=0.027$ with
$\max_{[0,\varepsilon]}\|R(\tau\mathbf{L})\|_{2}-1\approx 7.2$E-46.

##  3 ACKNOWLEDGMENTS

The authors were supported by the Austrian Science Fund (FWF) via the FWF-
funded SFB # F65. The third author acknowledges support from the FWF, grant
P33010. This work received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme, ERC
Advanced Grant NEUROMORPH, no. 101018153.

## References

  * [1] F. Achleitner, A. Arnold, and E. Carlen.  On multi-dimensional hypocoercive BGK models. Kinet. Relat. Models 11 (2018), 953–1009.
  * [2] F. Achleitner, A. Arnold, and E. A. Carlen.  The hypocoercivity index for the short time behavior of linear time-invariant ODE systems. J. Differ. Equ. 371 (2023), 83–115.
  * [3] F. Achleitner, A. Arnold, and E. Carlen.  The hypocoercivity index for the intermediate and large time behavior of ODEs. In preparation (2023).
  * [4] F. Achleitner, A. Arnold, A. Jüngel.  Necessary and sufficient conditions for strong stability of explicit Runge–Kutta methods. Submitted 2023. https://arxiv.org/abs/2308.05689
  * [5] F. Achleitner, A. Arnold, and V. Mehrmann.  Hypocoercivity and hypocontractivity concepts for linear dynamical systems. Electron. J. Linear Algebra 39 (2023) 33–61.
  * [6] C. Krattenthaler.  Advanced determinant calculus. A complement. Lin. Alg. Appl. 41 (2005), 68–166.
  * [7] Z. Sun and C.-W. Shu.  Strong stability of explicit Runge–Kutta time discretizations. SIAM J. Numer. Anal. 57 (2019), 1158–1182.
  * [8] E. Tadmor.  From semidiscrete to fully discrete: Stability of Runge–Kutta schemes by the energy method. II. In: D. Estep and S. Tavener (eds.), Collected Lectures on the Preservation of Stability under Discretization. Proc. Appl. Math. 109, pp. 25–49. SIAM, Philadelphia, 2002.
  * [9] C. Villani.  Hypocoercivity.  Memoirs Amer. Math. Soc., vol. 202. Amer. Math. Soc., Providence, 2009.

##  4 Appendix: Proof of Lemma 4

Let $p\in 4{\mathbb{N}}$ and let $e_{j}$ denote the $j$-th Euclidean basis
vector. By exploiting the special staircase form of the matrix $\mathbf{L}$ in
(2), we find that

$\mathbf{L}e_{1}=e_{2},\quad\mathbf{L}e_{j}=-e_{j-1}+e_{j+1}\quad\mbox{for
}j=2,\ldots,N-1,\quad\mathbf{L}e_{N}=-e_{N-1}-e_{N}.$

A computation shows that, for $j=1,\ldots,N-1$ (with $m=N-1$),

$e_{j}\in\ker\mathbf{L}_{H},\quad\sqrt{-\mathbf{L}_{H}}e_{j}=\sqrt{-\mathbf{L}_{H}}\mathbf{L}e_{j}=\cdots=\sqrt{-\mathbf{L}_{H}}\mathbf{L}^{m-j}e_{j}=0,\quad\sqrt{-\mathbf{L}_{H}}\mathbf{L}^{m-j+1}e_{j}\neq
0$ (10)

and $e_{N}\not\in\ker\mathbf{L}_{H}$. The coefficients of the matrix
$\mathbf{M}:=I-R(\tau\mathbf{L})^{*}R(\tau\mathbf{L})$ are given by
$M_{ij}=\langle e_{i},\mathbf{M}e_{j}\rangle$ for $i,j=1,\ldots,N$. By
definition (7) of the stability function and the condition $s=p$, we have

$\displaystyle\mathbf{M}$
$\displaystyle=\mathbf{I}-\bigg{(}\sum_{j=0}^{p}\frac{\tau^{j}}{j!}(\mathbf{L}^{*})^{j}\bigg{)}\,\bigg{(}\sum_{\ell=0}^{p}\frac{\tau^{\ell}}{\ell!}\mathbf{L}^{\ell}\bigg{)}$
$\displaystyle=-\sum_{\ell=1}^{p+1}\frac{\tau^{\ell}}{\ell!}\sum_{k=0}^{\ell}\binom{\ell}{k}(\mathbf{L}^{*})^{k}\mathbf{L}^{\ell-k}+\frac{\tau^{p+1}}{(p+1)!}\big{(}(\mathbf{L}^{*})^{p+1}+\mathbf{L}^{p+1}\big{)}+\mathcal{O}(\tau^{p+2})\quad\mbox{as
}\tau\to 0.$

The first term on the right-hand side can be rewritten by using [2, (A.6)]: as

$\displaystyle\mathbf{M}$
$\displaystyle=-2\sum_{\ell=1}^{p+1}\frac{\tau^{\ell}}{\ell!}\sum_{k=0}^{\ell-1}\binom{\ell-1}{k}(\mathbf{L}^{*})^{k}\mathbf{L}_{H}\mathbf{L}^{\ell-1-k}+\frac{\tau^{p+1}}{(p+1)!}\big{(}(\mathbf{L}^{*})^{p+1}+\mathbf{L}^{p+1}\big{)}+\mathcal{O}(\tau^{p+2}).$

Therefore, for $i,j=1,\ldots,N$,

$\displaystyle M_{ij}$
$\displaystyle=2\sum_{\ell=1}^{p+1}\frac{\tau^{\ell}}{\ell!}\sum_{k=0}^{\ell-1}\binom{\ell-1}{k}\Big{\langle}\sqrt{-\mathbf{L}_{H}}\mathbf{L}^{k}e_{i},\sqrt{-\mathbf{L}_{H}}\mathbf{L}^{\ell-1-k}e_{j}\Big{\rangle}$
$\displaystyle\phantom{xx}{}+\frac{\tau^{p+1}}{(p+1)!}\big{\langle}e_{i},\big{(}(\mathbf{L}^{*})^{p+1}+\mathbf{L}^{p+1}\big{)}e_{j}\big{\rangle}+\mathcal{O}(\tau^{p+2})=:I_{1}+I_{2}+\mathcal{O}(\tau^{p+2}).$

For $I_{1}$, due to (10), the terms in the sum with $k\leq m-i$ and
$\ell-1-k\leq m-j$ vanish, but those with $k\geq m-i+1$ and $\ell-1-k\geq
m-j+1$ may not vanish (recall that $m=p/2$). Since the $\tau$-order of the
leading nonvanishing term equals $\ell=m+2+k-j=2m+3-(i+j)$, we find that

$\displaystyle I_{1}$
$\displaystyle=-\frac{2\tau^{2m+3-(i+j)}}{(2m+3-(i+j))!}\binom{2m+2-(i+j)}{m-i+1}\big{\langle}e_{i},(\mathbf{L}^{*})^{m-i+1}\mathbf{L}_{H}\mathbf{L}^{m-j+1}e_{j}\big{\rangle}+\mathcal{O}(\tau^{2m+4-(i+j)}).$

Because of $p+1=2m+1$, the term $I_{2}$ is of leading order
$\tau^{2m+3-(i+j)}$ if and only if $i=j=1$. Therefore,

$\displaystyle M_{ij}$
$\displaystyle=-\frac{2\tau^{2m+3-(i+j)}}{(2m+3-(i+j))!}\binom{2m+2-(i+j)}{m-i+1}\big{\langle}e_{i},(\mathbf{L}^{*})^{m-i+1}\mathbf{L}_{H}\mathbf{L}^{m-j+1}e_{j}\big{\rangle}$
$\displaystyle\phantom{xx}{}+\delta_{1i}\delta_{j1}\frac{\tau^{p+1}}{(p+1)!}\big{\langle}e_{1},\big{(}(\mathbf{L}^{*})^{p+1}+\mathbf{L}^{p+1}\big{)}e_{1}\big{\rangle}+\mathcal{O}(\tau^{2m+4-(i+j)}),$

where $\delta_{i1}$ is the Kronecker delta.

We claim that the inner products in the first term of $M_{ij}$ are all equal
to $-1$. The inclusion
$\ker\sqrt{-\mathbf{L}_{H}}\supseteq\operatorname{span}\\{e_{1},\ldots,e_{N-1}\\}$
implies that

$\displaystyle\big{\langle}e_{i}$
$\displaystyle,(\mathbf{L}^{*})^{m-i+1}\mathbf{L}_{H}\mathbf{L}^{m-j+1}e_{j}\big{\rangle}=-\Big{\langle}\sqrt{-\mathbf{L}_{H}}\mathbf{L}^{m-i+1}e_{i},\sqrt{-\mathbf{L}_{H}}\mathbf{L}^{m-j+1}e_{j}\Big{\rangle}$
$\displaystyle=-\Big{\langle}\sqrt{-\mathbf{L}_{H}}\mathbf{L}^{N-i}e_{i},\sqrt{-\mathbf{L}_{H}}\mathbf{L}^{N-j}e_{j}\Big{\rangle}=-\Big{\langle}\sqrt{-\mathbf{L}_{H}}e_{N},\sqrt{-\mathbf{L}_{H}}e_{N}\Big{\rangle}=\langle
e_{N},\mathbf{L}_{H}e_{N}\rangle=-1\,.$ (11)

For the second term of $M_{ij}$ we compute similarly:

$\langle e_{1},\mathbf{L}^{p+1}e_{1}\rangle=\langle
e_{1},\mathbf{L}^{2m+1}e_{1}\rangle=\langle(\mathbf{L}^{*})^{m}e_{1},\mathbf{L}\mathbf{L}^{m}e_{1}\rangle=\langle(\mathbf{L}^{*})^{m}e_{1},\mathbf{L}_{H}\mathbf{L}^{m}e_{1}\rangle=-\langle\sqrt{-\mathbf{L}_{H}}(\mathbf{L}^{*})^{m}e_{1},\sqrt{-\mathbf{L}_{H}}\mathbf{L}^{m}e_{1}\rangle\,,$
(12)

where we used $\langle y,\mathbf{L}_{S}y\rangle=0$ for all
$y\in{\mathbb{C}}^{n}$. We recall from (4) that
$\sqrt{-\mathbf{L}_{H}}\mathbf{L}^{m}e_{1}=e_{N}$, and obtain analogously that

$\sqrt{-\mathbf{L}_{H}}(\mathbf{L}^{*})^{m}e_{1}=(-1)^{m}e_{N}=e_{N}\,,$

where we used the property $m=p/2\in 2{\mathbb{N}}$. Thus (12) equals $-1$. So
we conclude

$\begin{split}M_{ij}&=\frac{2\tau^{2m+3-(i+j)}}{(2m+3-(i+j))!}\binom{2m+2-(i+j)}{m-i+1}-2\delta_{1i}\delta_{j1}\frac{\tau^{2m+1}}{(2m+1)!}+\mathcal{O}(\tau^{2m+4-(i+j)})\\\
&=\bigg{(}m_{ij}\tau^{2m+3-(i+j)}-2\delta_{1i}\delta_{j1}\frac{\tau^{2m+1}}{(2m+1)!}\bigg{)}+\mathcal{O}(\tau^{2m+4-(i+j)})\,,\end{split}$

where we defined

$m_{ij}:=\frac{2}{(2m+3-(i+j))!}\binom{2m+2-(i+j)}{m-i+1}=\frac{2}{(2m+3-(i+j))(m-i+1)!(m-j+1)!}\,.$
(13)

We proceed with the computation of the determinant of $\mathbf{M}$. We use
$M_{ij}=\mathcal{O}(\tau^{2m+3-(i+j)})$ and Leibniz’ formula for determinants:

$\det\mathbf{M}=\sum_{\sigma\in
S_{N}}\operatorname{sgn}(\sigma)\prod_{i=1}^{N}M_{i,\sigma(i)}=\sum_{\sigma\in
S_{N}}\operatorname{sgn}(\sigma)\,\mathcal{O}(\tau^{\sum_{i=1}^{N}(2m+3-(i+\sigma(i)))})=\mathcal{O}(\tau^{(m+1)^{2}}),$

where $\operatorname{sgn}$ is the sign function of permutations in the
permutation group $S_{N}$ and

$\displaystyle\sum_{i=1}^{N}(2m+3-(i+\sigma(i)))=(m+1)(2m+3)-2\sum_{i=1}^{m+1}i=(m+1)(2m+3)-(m+1)(m+2)=(m+1)^{2}.$

This shows that

$\det\mathbf{M}=\tau^{(m+1)^{2}}\det\bigg{(}m_{ij}-\delta_{1i}\delta_{1j}\frac{2}{(2m+1)!}\bigg{)}_{1\leq
i,j\leq N}+\mathcal{O}(\tau^{(m+1)^{2}+1})\,,$ (14)

matching the claim (9). So it remains to determine the sign of the leading
term.

The Laplace expansion with respect to the first column yields for the
determinant on the right-hand side:

$\det\bigg{(}m_{ij}-\delta_{1i}\delta_{1j}\frac{2}{(2m+1)!}\bigg{)}_{1\leq
i,j\leq N}=\det(m_{ij})_{1\leq i,j\leq N}-\frac{2}{(2m+1)!}\det(m_{ij})_{2\leq
i,j\leq N}\,.$ (15)

The first determinant on the right-hand side becomes, by Leibniz’ formula for
determinants and definition (13) of $m_{ij}$,

$\displaystyle\det(m_{ij})_{1\leq i,j\leq N}$ $\displaystyle=\sum_{\sigma\in
S_{N}}\operatorname{sgn}(\sigma)\prod_{i=1}^{N}\frac{2}{(2m+3-(i+\sigma(i)))(m-i+1)!(m-\sigma(i)+1)!}$
$\displaystyle=2^{N}\bigg{(}\prod_{i=1}^{N}\frac{1}{(m-i+1)!}\bigg{)}^{2}\sum_{\sigma\in
S_{N}}\operatorname{sgn}(\sigma)\prod_{i=1}^{N}\frac{1}{2m+3-(i+\sigma(i))}$
$\displaystyle=2^{N}\bigg{(}\prod_{i=1}^{N}\frac{1}{(m-i+1)!}\bigg{)}^{2}\det\bigg{(}\frac{1}{2m+3-(i+j)}\bigg{)}_{1\leq
i,j\leq N}.$

The determinant of the Hankel matrix can be determined explicitly from
Cauchy’s double alternant (see, e.g., [6, Theorem 12]):

$\Delta_{1}:=\det\bigg{(}\frac{1}{2m+3-(i+j)}\bigg{)}_{1\leq i,j\leq
N}=\det\bigg{(}\frac{1}{i+(j-1)}\bigg{)}_{1\leq i,j\leq N}=\frac{\prod_{1\leq
i<j\leq N}(i-j)^{2}}{\prod_{i,j=1}^{N}(i+(j-1))}>0\,.$

The second determinant on the right-hand side of (15) is treated similarly:

$\det(m_{ij})_{2\leq i,j\leq
N}=2^{N-1}\bigg{(}\prod_{i=2}^{N}\frac{1}{(m-i+1)!}\bigg{)}^{2}\Delta_{2},$

where

$\Delta_{2}:=\det\bigg{(}\frac{1}{2m+3-(i+j)}\bigg{)}_{2\leq i,j\leq
N}=\det\bigg{(}\frac{1}{i+(j-1)}\bigg{)}_{1\leq i,j\leq N-1}\\\
=\frac{\prod_{1\leq i<j\leq N-1}(i-j)^{2}}{\prod_{i,j=1}^{N-1}(i+(j-1))}.$

In fact, we can express $\Delta_{1}$ in terms of $\Delta_{2}$:

$\Delta_{1}=\frac{\prod_{1\leq i<j\leq
N-1}(i-j)^{2}}{\prod_{i,j=1}^{N-1}(i+(j-1))}\
\frac{(2N-1)\prod_{i=1}^{N-1}(i-N)^{2}}{\prod_{i=1}^{N}(i+N-1)^{2}}=\frac{\Delta_{2}}{\Delta_{3}},$

where, because of $N=m+1$,

$\displaystyle\Delta_{3}$
$\displaystyle:=\frac{\prod_{i=1}^{N}(i+N-1)^{2}}{(2N-1)\prod_{i=1}^{N-1}(i-N)^{2}}=\frac{\big{(}(2N-1)!/(N-1)!\big{)}^{2}}{(2N-1)(N-1)!^{2}}$
$\displaystyle=\frac{(2N-1)!}{(N-1)!^{2}}\frac{(2N-2)!}{(N-1)!^{2}}=\frac{(2m+1)!}{m!^{2}}\frac{(2m)!}{m!^{2}}=\frac{(2m+1)!}{m!^{2}}\binom{2m}{m}.$

Therefore, (15) becomes

$\displaystyle
c:=\det\bigg{(}m_{ij}-\delta_{1i}\delta_{1j}\frac{2}{(2m+1)!}\bigg{)}_{1\leq
i,j\leq N}$
$\displaystyle=2^{N}\bigg{(}\prod_{i=1}^{N}\frac{1}{(m-i+1)!}\bigg{)}^{2}\Delta_{1}-\frac{2^{N}}{(2m+1)!}\bigg{(}\prod_{i=2}^{N}\frac{1}{(m-i+1)!}\bigg{)}^{2}\Delta_{2}$
$\displaystyle=2^{N}\bigg{(}\prod_{i=1}^{N}\frac{1}{(m-i+1)!}\bigg{)}^{2}\Delta_{1}\bigg{(}1-\frac{(m!)^{2}}{(2m+1)!}\Delta_{3}\bigg{)}$
$\displaystyle=2^{N}\bigg{(}\prod_{i=1}^{N}\frac{1}{(m-i+1)!}\bigg{)}^{2}\Delta_{1}\bigg{(}1-\binom{2m}{m}\bigg{)}<0\,.$

Thus, (14) shows that, for sufficiently small $\tau>0$, the sign of
$\det\mathbf{M}$ is

$\operatorname{sgn}(\det\mathbf{M})=\operatorname{sgn}\bigg{(}1-\binom{2m}{m}\bigg{)}<0,$

finishing the proof. $\hfill\square$