# Leading-Color Two-Loop QCD Corrections for Three-Photon Production at Hadron Colliders

S. Abreu, B. Page, E. Pascual, V. Sotnikov

► To cite this version:

S. Abreu, B. Page, E. Pascual, V. Sotnikov. Leading-Color Two-Loop QCD Corrections for Three-Photon Production at Hadron Colliders. Journal of High Energy Physics, 2021, 2021 (01), pp.078. 10.1007/JHEP01(2021)078. hal-03011024

HAL Id: hal-03011024

https://hal.science/hal-03011024v1

Submitted on 10 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Leading-color two-loop QCD corrections for three-photon production at hadron colliders

S. Abreu,$^{a,b}$ B. Page,$^c$ E. Pascual$^d$ and V. Sotnikov$^e$

$^a$Theoretical Physics Department, CERN,
1211 Geneva 23, Switzerland

$^b$Mani L. Bhaumik Institute for Theoretical Physics,
Department of Physics and Astronomy, UCLA, Los Angeles, CA 90095, U.S.A.

$^c$Laboratoire de physique de l'Ecole normale supérieure, ENS,
Université PSL, CNRS, Sorbonne Université, Université Paris-Diderot,
Sorbonne Paris Cité, 24 rue Lhomond, 75005 Paris, France

$^d$Physikalisches Institut, Albert-Ludwigs-Universität Freiburg,
Hermann-Herder-Str. 3, D-79104 Freiburg, Germany

$^e$Max Planck Institute for Physics (Werner Heisenberg Institute),
D-80805 Munich, Germany

E-mail: samuel.abreu@cern.ch, bpage@ipht.fr,
evgenij.pascual@saturn.uni-freiburg.de, sotnikov@mpp.mpg.de

**ABSTRACT:** We compute the two-loop helicity amplitudes for the production of three photons at hadron colliders in QCD at leading-color. Using the two-loop numerical unitarity method coupled with analytic reconstruction techniques, we obtain the decomposition of the two-loop amplitudes in terms of master integrals in analytic form. These expressions are valid to all orders in the dimensional regulator. We use them to compute the two-loop finite remainders, which are given in a form that can be efficiently evaluated across the whole physical phase space. We further package these results in a public code which assembles the helicity-summed squared two-loop remainders, whose numerical stability across phase-space is demonstrated. This is the first time that a five-point two-loop process is publicly available for immediate phenomenological applications.

**KEYWORDS:** NLO Computations, QCD Phenomenology

**ARXIV EPRINT:** 2010.15834
# Contents

<table><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Notation and conventions</td><td>3</td></tr><tr><td>2.1 Helicity amplitudes</td><td>3</td></tr><tr><td>2.2 Finite remainders</td><td>6</td></tr><tr><td>3 Calculation of helicity amplitudes</td><td>7</td></tr><tr><td>3.1 Two-loop numerical unitarity</td><td>7</td></tr><tr><td>3.2 Pentagon functions</td><td>8</td></tr><tr><td>3.3 Functional reconstruction of master integral coefficients</td><td>9</td></tr><tr><td>3.4 Two-loop finite remainders</td><td>11</td></tr><tr><td>3.5 Reference values</td><td>12</td></tr><tr><td>3.6 Validation</td><td>14</td></tr><tr><td>4 Numerical evaluation</td><td>14</td></tr><tr><td>5 Conclusion</td><td>16</td></tr><tr><td>A Master integral basis</td><td>17</td></tr><tr><td>B Collinear checks</td><td>20</td></tr><tr><td>B.1 One-loop photon splitting amplitudes</td><td>23</td></tr></table>

## 1 Introduction

Precise theoretical predictions for scattering experiments at particle colliders crucially rely on the availability of higher-order scattering amplitudes. Over the next few years, the large increase in the amount of data collected by the experiments at the Large Hadron Collider (LHC) at CERN will translate into measurements made at unprecedented accuracy. To make the most of the physics program of the LHC it is fundamental that theoretical predictions reach comparable levels of precision. A large number of $2 \to 2$ scattering processes are now available including next-to-next-to-leading (NNLO) order QCD corrections (see e.g. [1] for a recent review). However, reaching the same level of precision for higher multiplicity processes still remains a formidable challenge, with the best available predictions only including next-to-leading order (NLO) QCD and electro-weak corrections.

In this work, we focus on the double-virtual NNLO QCD corrections to the production of three photons at hadron colliders. This process is important for studying various beyond-the-Standard-Model (BSM) phenomena. In particular, it can be used to constrain
anomalous quartic gauge [2] and Higgs couplings [3–5]. Furthermore, it is a major contribution to the irreducible background in searches for associated production of BSM particles and a photon [6–8]. Similar to the related two-photon production process [9–14], three-photon production exhibits very slow perturbative convergence [15, 16] and non-reliable estimates of uncertainty from missing higher orders. In particular, a significant tension between LHC Run 1 data and NLO QCD predictions has been reported [17]. The computation of NNLO QCD corrections is thus crucial for this process. While several NNLO infrared subtraction approaches are capable of tackling triphoton production [15, 16], obtaining the two-loop five-point amplitudes required for the double-virtual contributions in a representation that allows for an efficient on-the-fly computation remains a major obstacle. Our goal is to provide such expressions, and by doing so making these amplitudes widely available for any future phenomenological studies involving this process.

In the last few years, there has been remarkable progress in the calculation of two-loop five-point amplitudes. A basis of master integrals relevant for the scattering of five massless particles has been known for a number of years in the planar limit [18, 19], and more recently including non-planar topologies [20, 21]. These master integrals evaluate to a class of multi-valued special functions with logarithmic branch cuts. Given the large number of scales in five-point kinematics, it is essential to find a representation of the master integrals in terms of special functions that manifests the physical properties of amplitudes. This was achieved by defining a set of so-called pentagon functions, first for the planar integral topologies [22], and more recently for a complete set of five-point massless master integrals [23]. The latter work also provides an efficient code for their numerical evaluation across all of phase space. At the same time, there has been substantial progress in the reduction of five-point two-loop scattering amplitudes to master integrals [24–36], often building on the use of finite-field arithmetic [26, 37], and functional reconstruction techniques [26, 31, 38–40]. As a result, all five-parton two-loop planar QCD amplitudes were computed numerically [41–45] and in analytic form [46–49]. Despite the remarkable progress of the last few years, five-point amplitudes are not yet broadly available in a form suitable for direct phenomenological applications. Indeed, the only phenomenological study involving five-point two-loop amplitudes is the work of [15], where the authors evaluated planar two-loop $q\bar{q} \to \gamma\gamma\gamma$ amplitudes on a sufficient set of phase-space points to construct interpolating functions.

In this paper, we calculate a complete set of independent planar two-loop helicity amplitudes required for the double-virtual NNLO QCD corrections to three-photon production at hadron colliders. We present them for the first time in a form which is suitable for direct and flexible phenomenological applications. Indeed, our results have already been employed in ref. [16] to implement the first on-the-fly computation of NNLO QCD corrections to triphoton production within the MATRIX framework [50]. We obtain the two-loop helicity amplitudes by following a similar approach to that used in refs. [47, 49], based on the two-loop numerical unitarity approach [43, 51] as implemented in the recently released C++ library CARAVEL [36]. Unlike in refs. [47, 49], where only the finite remainders were obtained, in this paper we obtain analytic expressions for the decomposition of the amplitudes into master integrals, marking the first time that such results are
available for five-point two-loop amplitudes. This decomposition is valid to all orders in the dimensional regulator. They are valuable for studying the analytic complexity of the master-integral decomposition and for future computations of higher-order corrections. We also obtain analytic expressions for the two-loop finite remainders, decomposed into the pentagon functions of ref. [23]. We confirm that the latter are more suitable for efficient and stable numerical evaluations. All the analytic results we obtain are made available in the supplementary material, and we provide a public C++ library [52] for the efficient numerical evaluation of the finite remainders.

The paper is organized as follows. In section 2 we establish our conventions and define the objects we will be computing. In section 3, we present our approach to the calculation of the two-loop amplitudes and their finite remainders. We discuss some properties of the analytic results we obtained, give reference evaluation values and discuss the checks that were performed. In section 4, we present a public implementation of the numerical evaluation of the analytic results for finite remainders and demonstrate its efficiency and numerical stability.

# 2 Notation and conventions

## 2.1 Helicity amplitudes

We consider a complete set of helicity amplitudes required for the computation of the double-virtual next-to-next-to leading order (NNLO) QCD corrections to the production of three photons at hadron colliders. More precisely, we consider the parton-level scattering process

$$q(p_1, h_1) + \bar{q}(p_2, h_2) \rightarrow \gamma(p_3, h_3) + \gamma(p_4, h_4) + \gamma(p_5, h_5). \qquad (2.1)$$

This is the only (sub-)process required for these corrections, as the loop-induced process $gg \rightarrow \gamma\gamma\gamma$ vanishes to all orders in the coupling constants due to the charge conjugation symmetry of QCD$\otimes$QED. We denote the momenta and the helicity states of the particles as $p_i$ and $h_i$ respectively. All particles are massless, and the kinematics of the process are thus fully specified by the five independent Mandelstam invariants

$$\begin{align}
s_{12} &= (p_1 + p_2)^2, & s_{23} &= (p_2 + p_3)^2, & s_{34} &= (p_3 + p_4)^2, \nonumber \\
s_{45} &= (p_4 + p_5)^2, & s_{15} &= (p_1 + p_5)^2, & & \tag{2.2}
\end{align}$$

and the parity-odd contraction of four momenta

$$\mathrm{tr}_5 := 4i\epsilon(p_1, p_2, p_3, p_4), \qquad (2.3)$$

where $\epsilon(\cdot, \cdot, \cdot, \cdot)$ is the fully anti-symmetric Levi-Civita symbol. The physical region associated with the process in eq. (2.1), where momenta $p_1$ and $p_2$ are incoming, is characterized by

$$s_{12}, s_{34}, s_{45} > 0, \quad s_{23}, s_{15} < 0, \quad \mathrm{tr}_5^2 < 0. \qquad (2.4)$$

The latter condition is equivalent to the negativity of the five-point Gram determinant and it is trivially satisfied by real-valued momenta in eq. (2.3).
The amplitudes for this process, denoted as $M(1_q^{h_1}, 2_{\bar{q}}^{h_2}, 3_{\gamma}^{h_3}, 4_{\gamma}^{h_4}, 5_{\gamma}^{h_5})$, can be decomposed into a color factor, a helicity-dependent spinor weight, and a Lorentz invariant kinematic factor. That is, we can write

$$M(1_q^{h_1}, 2_{\bar{q}}^{h_2}, 3_{\gamma}^{h_3}, 4_{\gamma}^{h_4}, 5_{\gamma}^{h_5}) := e_q^3 \delta_{i_1 i_2} \Phi(1_q^{h_1}, 2_{\bar{q}}^{h_2}, 3_{\gamma}^{h_3}, 4_{\gamma}^{h_4}, 5_{\gamma}^{h_5}) A(1_q^{h_1}, 2_{\bar{q}}^{h_2}, 3_{\gamma}^{h_3}, 4_{\gamma}^{h_4}, 5_{\gamma}^{h_5}), \quad (2.5)$$

where $i_1$ and $i_2$ are color indices of the quark and the antiquark, $e_q$ is the electric charge of the external quark in the process, and $\Phi$ denotes the spinor-weight factor. In the following, we will call $A(1_q^{h_1}, 2_{\bar{q}}^{h_2}, 3_{\gamma}^{h_3}, 4_{\gamma}^{h_4}, 5_{\gamma}^{h_5})$ the helicity amplitudes for the process in eq. (2.1). For simplicity, we will often suppress the arguments of $A$. We employ the 't Hooft-Veltman scheme of dimensional regularization with $D = 4-2\epsilon$ space-time dimensions to regularize infrared and ultraviolet divergences. We define the dimensionally regularized helicity amplitudes with external quarks as in [45].

The bare helicity amplitudes have a perturbative expansion in powers of the bare strong coupling $\alpha_s^0 = (g_s^0)^2/(4\pi)$, which we write as

$$A = A^{(0)} + \frac{\alpha_s^0}{2\pi} A^{(1)} + \left(\frac{\alpha_s^0}{2\pi}\right)^2 A^{(2)} + \mathcal{O}\left((\alpha_s^0)^3\right). \quad (2.6)$$

The renormalized coupling $\alpha_s$ is related to the bare $\alpha_s^0$ through

$$\alpha_s^0 \mu_0^{2\epsilon} S_{\epsilon} = \alpha_s \mu^{2\epsilon} \left(1 - \frac{\beta_0}{\epsilon} \frac{\alpha_s}{2\pi} + \mathcal{O}\left(\alpha_s^2\right)\right), \quad S_{\epsilon} = (4\pi)^{\epsilon} e^{-\epsilon \gamma_E}, \quad (2.7)$$

where $\gamma_E$ is the Euler-Mascheroni constant, and $\mu_0$ and $\mu$ are the dimensional regularization and renormalization scales, which from now on we assume to be equal. $\beta_0$ is the first coefficient of the QCD $\beta$-function,

$$\beta_0 = \frac{11C_A - 4T_F N_f}{6}, \quad (2.8)$$

where $C_A = N_c$ is the quadratic Casimir of the adjoint representation of the SU($N_c$) group, and $T_F = 1/2$ is the normalization of fundamental representation generators. Below we will also need the quadratic Casimir of the fundamental representation, $C_F = \frac{N_c^2-1}{2N_c}$. We define the perturbative expansion of the renormalized amplitudes as

$$A_R = A_R^{(0)} + \frac{\alpha_s}{2\pi} A_R^{(1)} + \left(\frac{\alpha_s}{2\pi}\right)^2 A_R^{(2)} + \mathcal{O}(\alpha_s^3). \quad (2.9)$$

The coefficients $A_R^{(k)}$ are then related to their bare counterparts as

$$A_R^{(0)} = A^{(0)}, \quad A_R^{(1)} = S_{\epsilon}^{-1} A^{(1)}, \quad A_R^{(2)} = S_{\epsilon}^{-2} A^{(2)} - \frac{\beta_0}{\epsilon} S_{\epsilon}^{-1} A^{(1)}. \quad (2.10)$$

There are 16 different helicity configurations to consider. However, it can be easily shown that only 2 of them are independent. In this paper we choose the independent configurations to be

$$\begin{aligned} A_{++}(1, 2, 3, 4, 5) &:= A(1_q^+, 2_{\bar{q}}^-, 3_{\gamma}^+, 4_{\gamma}^+, 5_{\gamma}^+), \\ A_{--}(1, 2, 3, 4, 5) &:= A(1_q^-, 2_{\bar{q}}^-, 3_{\gamma}^-, 4_{\gamma}^+, 5_{\gamma}^+), \end{aligned} \quad (2.11)$$
<table><thead><tr><th>Helicity</th><th>Expression</th></tr></thead><tbody><tr><td>++</td><td>A<sub>+++</sub>(1, 2, 3, 4, 5)</td></tr><tr><td>+-</td><td>A<sub>-++</sub>(1, 2, 3, 4, 5)</td></tr><tr><td>+-</td><td>A<sub>-++</sub>(1, 2, 4, 3, 5)</td></tr><tr><td>+-</td><td>A<sub>-++</sub>(1, 2, 5, 3, 4)</td></tr><tr><td>+-</td><td>P A<sub>-++</sub>(2, 1, 5, 3, 4)</td></tr><tr><td>+-</td><td>P A<sub>-++</sub>(2, 1, 3, 4, 5)</td></tr><tr><td>+-</td><td>P A<sub>-++</sub>(2, 1, 4, 3, 5)</td></tr><tr><td>+-</td><td>P A<sub>+++</sub>(2, 1, 3, 4, 5)</td></tr></tbody></table>

<table><caption>**Table 1.** Relation between all helicity configurations and the two independent basis elements in eq. (2.11). **P** denotes a parity transformation.</caption><thead><tr><th>Helicity</th><th>Expression</th></tr></thead><tbody><tr><td>++</td><td>A<sub>+++</sub>(2, 1, 3, 4, 5)</td></tr><tr><td>+-</td><td>A<sub>-++</sub>(2, 1, 3, 4, 5)</td></tr><tr><td>+-</td><td>A<sub>-++</sub>(2, 1, 4, 3, 5)</td></tr><tr><td>+-</td><td>A<sub>-++</sub>(2, 1, 5, 3, 4)</td></tr><tr><td>+-</td><td>P A<sub>-++</sub>(1, 2, 5, 3, 4)</td></tr><tr><td>+-</td><td>P A<sub>-++</sub>(1, 2, 3, 4, 5)</td></tr><tr><td>+-</td><td>P A<sub>-++</sub>(1, 2, 4, 3, 5)</td></tr><tr><td>+-</td><td>P A<sub>+++</sub>(1, 2, 3, 4, 5)</td></tr></tbody></table>

where we indexed each independent amplitude by the photon helicities. In table 1, we show how permutations of momenta, parity and charge conjugation can be used to relate all helicity configurations to these two amplitudes. Note that all momentum permutations are within the scattering region defined in eq. (2.4). We choose $\Phi_{-++}$ such that $A_{-++}^{(0)} = 1$ and we choose

$$ \Phi_{+++} = \frac{[31]\langle 12 \rangle^3 \langle 13 \rangle}{\langle 14 \rangle^2 \langle 15 \rangle^2 \langle 23 \rangle^2}. \qquad (2.12) $$

Each helicity amplitude can be further decomposed into individually gauge-invariant contributions as

$$ \begin{gathered} \mathcal{A}^{(1)} = C_F A^{(1)}, \\ \mathcal{A}^{(2)} = C_F^2 B^{(2,0)} + C_F C_A B^{(2,1)} + C_F T_F N_f A^{(2,N_f)} + C_F T_F \left( \sum_{f=1}^{N_f} Q_f^2 \right) \tilde{A}^{(2,N_f)}, \quad (2.13) \end{gathered} $$

where $N_f$ is the number of light quarks, and $Q_f$ is the ratio of the electric charges of the quark with flavor $f$ and the quark in the initial state. Here we do not consider the contributions from heavy quark loops. In figure 1 we depict representative diagrams for each of these contributions.<fnref n="1" />

In this paper, we compute the two-loop amplitudes $\mathcal{A}^{(2)}$ in the leading-color approximation, where the number of colors $N_c$ is large with the ratio $N_f/N_c$ kept constant. In this approximation only planar topologies contribute. We also do not consider the gauge-invariant term $\tilde{A}^{(2,N_f)}$ which includes non-planar contributions. The validity of this approximation for phenomenology is discussed in ref. [16]. We write the two-loop amplitudes $\mathcal{A}^{(2)}$ as

$$ \mathcal{A}^{(2)} = \frac{N_c^2}{4} (A^{(2,0)} + \mathcal{O}(N_c^{-2})) + C_F T_F N_f A^{(2,N_f)}, \quad A^{(2,0)} = B^{(2,0)} + 2B^{(2,1)}. \quad (2.14) $$

In summary, the main goal of this paper consists in the calculation of four functions: $A_{+++}^{(2,0)}$, $A_{+++}^{(2,N_f)}$, $A_{-++}^{(2,0)}$ and $A_{-++}^{(2,N_f)}$. We note that since $A_{+++}^{(0)} = 0$, only the latter two are required to compute NNLO QCD corrections to the process in eq. (2.1).

<footnote n="1">¹Contributions similar to those in figure 1(d) but where either a single or three photons attach to the closed quark loop vanish due to charge-conjugation symmetry of QED$\otimes$QCD.</footnote>
<figure><img src="image_2.png" /></figure>

(a) Contribution proportional to $C_F^2$

<figure><img src="image_3.png" /></figure>

(b) Contribution proportional to $C_F C_A$

<figure><img src="image_4.png" /></figure>

(c) Contribution proportional to $N_f$

<figure><img src="image_5.png" /></figure>

(d) Contribution proportional to $\sum_f Q_f^2$

**Figure 1.** Representative diagrams for each of the contributions in eq. (2.13). Photons ($\gamma$) are denoted by wavy lines, gluon by curly lines, and quarks by straight lines.

## 2.2 Finite remainders

At NNLO, the two-loop amplitudes contribute to the cross section only through finite remainders (see e.g. [53]), which can be expanded similarly to eq. (2.9) as

$$ \mathcal{R} = \mathcal{R}^{(0)} + \frac{\alpha_s}{2\pi} \mathcal{R}^{(1)} + \left(\frac{\alpha_s}{2\pi}\right)^2 \mathcal{R}^{(2)} + \mathcal{O}(\alpha_s^3). \qquad (2.15) $$

The coefficients $\mathcal{R}^{(i)}$ are obtained from the expansion of the renormalized amplitudes $\mathcal{A}_R$ by subtracting the remaining infrared singularities [54]. More precisely, we define

$$ \begin{aligned} \mathcal{R}^{(0)} &= \mathcal{A}_R^{(0)}, \\ \mathcal{R}^{(1)} &= \mathcal{A}_R^{(1)} - \mathbf{I}^{(1)}\mathcal{A}_R^{(0)} + \mathcal{O}(\epsilon), \\ \mathcal{R}^{(2)} &= \mathcal{A}_R^{(2)} - \mathbf{I}^{(1)}\mathcal{A}_R^{(1)} - \mathbf{I}^{(2)}\mathcal{A}_R^{(0)} + \mathcal{O}(\epsilon), \end{aligned} \qquad (2.16) $$

where the functions $\mathbf{I}^{(1)}$ and $\mathbf{I}^{(2)}$ are defined as (see e.g. [9, 10])

$$ \begin{aligned} \mathbf{I}^{(1)}(\epsilon) &= -C_F \frac{e^{\gamma_E \epsilon}}{\Gamma(1-\epsilon)} \left(\frac{1}{\epsilon^2} + \frac{3}{2\epsilon}\right) \left(-\frac{s_{12}}{\mu^2} - i0\right)^{-\epsilon}, \\ \mathbf{I}^{(2)}(\epsilon) &= -\frac{1}{2} \mathbf{I}^{(1)}(\epsilon) \mathbf{I}^{(1)}(\epsilon) - \frac{\beta_0}{\epsilon} \mathbf{I}^{(1)}(\epsilon) + \frac{e^{-\gamma_E \epsilon} \Gamma(1-2\epsilon)}{\Gamma(1-\epsilon)} \left(\frac{\beta_0}{\epsilon} + K\right) \mathbf{I}^{(1)}(2\epsilon) + \mathbf{H}(\epsilon). \end{aligned} \qquad (2.17) $$

In $\mathbf{I}^{(2)}(\epsilon)$, we have introduced the functions

$$ \begin{aligned} K &= \left(\frac{67}{18} - \frac{\pi^2}{6}\right) C_A - \frac{10}{9} T_F N_F, & \mathbf{H}(\epsilon) &= \frac{e^{\gamma_E \epsilon}}{2\epsilon\Gamma(1-\epsilon)} H_q, \\ H_q &= \left(\frac{\pi^2}{2} - 6\zeta_3 - \frac{3}{8}\right) C_F^2 + \left(\frac{13}{2}\zeta_3 + \frac{245}{216} - \frac{23}{48}\pi^2\right) C_A C_F + \left(\frac{\pi^2}{12} - \frac{25}{54}\right) T_F C_F N_f. \end{aligned} \qquad (2.18) $$
The finite remainders can be decomposed similarly to eq. (2.14), and we write

$$ 
\begin{align} 
\mathcal{R}^{(1)} &= C_F R^{(1)}, \nonumber \\ 
\mathcal{R}^{(2)} &= \frac{N_c^2}{4} (R^{(2,0)} + \mathcal{O}(N_c^{-2})) + C_F T_F N_f R^{(2,N_f)}. \tag{2.19}
\end{align} 
$$

Analogously to eq. (2.11), whenever convenient we will also write $R_{++}^{(i,j)}$ and $R_{-++}^{(i,j)}$ to denote the contributions to each of the two independent helicity states.

## 3 Calculation of helicity amplitudes

To compute the functions $A_{++}^{(2,0)}$, $A_{++}^{(2,N_f)}$, $A_{-+}^{(2,0)}$ and $A_{-+}^{(2,N_f)}$ defined in eq. (2.14), we use the framework of two-loop numerical unitarity [25, 43, 51, 55] coupled with functional reconstruction techniques. The same approach was already used previously to compute the five-parton two-loop amplitudes [47, 49]. We build on the implementation of this framework in CARAVEL [36], which we modify to handle amplitudes with external photons. In this section, we summarize the main steps of our calculation.

Before delving into this, however, we note that our approach requires knowing the one-loop amplitudes to order $\epsilon^2$. We have computed them to all orders in $\epsilon$ using the same techniques as those used for the two-loop amplitudes. This is by now an easy calculation so we will not discuss it further, and simply include the results in the supplementary material.

### 3.1 Two-loop numerical unitarity

Our approach to the calculation of two-loop amplitudes is built on their numerical evaluation within the framework of two-loop numerical unitarity. We target independently each of the helicity amplitudes $A_h^{(2,j)}$ with $j = \{0, N_f\}$, $h = \{+++, -++\}$, see eq. (2.14). The approach starts from a parametrization of the integrand of the amplitude $A_h^{(2,j)}(\ell_l)$ in terms of master integrands and surface terms [25] ($\ell_l$ denotes the loop momenta). The master integrands are associated with master integrals, and the surface terms integrate to zero. This decomposition is naturally organized in terms of propagator structures. More precisely, we write

$$ A_h^{(2,j)}(\ell_l) = \sum_{\Gamma \in \Delta} \sum_{i \in M_{\Gamma} \cup S_{\Gamma}} c_{\Gamma, i} \frac{m_{\Gamma, i}(\ell_l)}{\prod_{j \in P_{\Gamma}} \rho_j}, \quad (3.1) $$

where $\Delta$ is the set of propagator structures $\Gamma$, $P_\Gamma$ is the multiset of inverse propagators $\rho_j$ in $\Gamma$, and $M_\Gamma$ and $S_\Gamma$ denote the sets of master integrands and surface terms. The coefficients $c_{\Gamma,i}$ are determined using the factorization properties of the integrand $A^{(2,j)}(\ell_l)$ in specific configurations $\ell_l^\Gamma$ of the loop momenta where the inverse propagators $\rho_j \in P_\Gamma$ are on-shell, that is $\rho_j(\ell_l^\Gamma) = 0$ iff $j \in P_\Gamma$. In this limit, the leading contribution to eq. (3.1) factorizes as

$$ \sum_{\text{states}} \prod_{i \in T_{\Gamma}} A_i^{(0)}(\ell_l^\Gamma) = \sum_{\Gamma' \ge \Gamma, i \in M_{\Gamma'} \cup S_{\Gamma'}} \frac{c_{\Gamma',i} m_{\Gamma',i}(\ell_l^\Gamma)}{\prod_{j \in (P_{\Gamma'} \setminus P_{\Gamma})} \rho_j(\ell_l^\Gamma)}. \quad (3.2) $$

The sum on the right-hand side is over the propagator structures $\Gamma'$ such that $P_\Gamma \subseteq P_{\Gamma'}$. On the left-hand side, $T_\Gamma$ denotes the set of tree amplitudes associated with the vertices in
the diagram corresponding to $\Gamma$ and the sum is over the scheme-dependent physical states of each internal line of $\Gamma$.

We refer to eq. (3.2) as cut equations, and to its left hand side as cuts. The coefficients $c_{\Gamma,i}$ are determined by sampling the cut equations over enough values of $\ell_l^\Gamma$. To construct the color-stripped products of tree amplitudes on the left-hand side of eq. (3.2), we use the unitarity-based color decomposition approach of [56, 57] to include colorless particles. The dimensional-regulator dependence of cuts is determined with the approach of decomposition by particle content [49, 58, 59], based on dimensional reduction. We evaluate the color-stripped and dimensional-regulator-free tree amplitudes through Berends-Giele recursion [60]. The resulting system of equations is then solved numerically on a given phase-space point. All numerical operations are done using finite-field arithmetic. This allows us to obtain exact coefficients $c_{\Gamma,i}$ for rational phase-space points. The latter can be generated, for instance, by using momentum-twistor variables [61]. Once the $c_{\Gamma,i}$ have been determined, we obtain the decomposition of the amplitude in terms of master integrals,

$$A_h^{(2,j)} = \sum_{\Gamma \in \Delta} \sum_{i \in M_{\Gamma}} c_{\Gamma,i} m_{\Gamma,i}, \qquad (3.3)$$

where $m_{\Gamma,i}$ is the master integral associated with the numerator $m_{\Gamma,i}(\ell_l)$.

Compared to five-parton amplitudes, we note that the fact that photons cannot be ordered leads to a proliferation of topologies. To be more explicit, consider as an example the two-loop leading-color color-ordered five-gluon amplitudes. Due to the color ordering, there are only 5 different pentagon-box topologies<fnref n="2" /> to consider, which can be labeled e.g. by the vertex attached to the gluon of momentum $p_1$. For the photon amplitudes we are considering in this paper, we should multiply this number by 3!, corresponding to the possible orderings of the photons. We thus have to consider 30 different pentagon-box topologies.

## 3.2 Pentagon functions

A basis of master integrals for planar five-point massless amplitudes is known [18, 19]. Order by order in $\epsilon$, they can be expressed in terms of multiple polylogarithms. MPLs form a special class of functions with logarithmic singularities, and this class of functions can be equipped with algebraic structures which allow one to find relations between them [62–64]. As a consequence, we can define a set of functions, called pentagon functions, which form a basis for the MPLs that appear in the master integrals contributing to planar two-loop five-point massless amplitudes [22, 23].

After evaluating the coefficients $c_{\Gamma,i}$ in eq. (3.1), we directly obtain the decomposition of the helicity amplitudes in terms of master integrals, see eq. (3.3). Assuming that the decomposition of the master integrals in terms of pentagon functions is known, we in turn obtain a decomposition of the amplitude in terms of pentagon functions. If we denote the pentagon functions by $\{h_i\}_{i \in B}$, with $B$ the associated set of labels, we can then decompose

<footnote n="2">²We call pentagon-box the topology corresponding to the diagram shown in figure 1(a).</footnote>
the amplitude as

$$A_h^{(2,j)} = \sum_{i \in B} \sum_{k=-4}^{0} \epsilon^k d_{k,i} h_i + \mathcal{O}(\epsilon^4), \quad (3.4)$$

where we make explicit that two-loop amplitudes have at most poles of order $\epsilon^{-4}$. The decomposition of eq. (3.4) presents a major advantage compared to the decomposition of an amplitude in terms of master integrals: it allows us to write one-loop and two-loop amplitudes in terms of a common basis of functions. It then follows from eqs. (2.16) and (2.19) that remainders themselves can be decomposed in terms of pentagon functions,<fnref n="3" /> that is,

$$R_h^{(2,j)} = \sum_{i \in B} r_i h_i. \qquad (3.5)$$

In this paper, we will adopt the pentagon functions defined in ref. [23]. Aside from the fact that they can be efficiently evaluated across phase-space, they are also defined to be a basis for the whole orbit under the symmetry group of five-point kinematics in the {1, 2}-channel (see eq. (2.4)). This is crucial for our calculation. As already mentioned previously, the amplitudes receive contributions from all orderings of the photons in the final state. An important consequence of this observation is that there is no Euclidean region for the amplitude, that is no region where the amplitude is real.

### 3.3 Functional reconstruction of master integral coefficients

In refs. [47, 49], we argued that having two-loop corrections in mind one should reconstruct the analytic form of the coefficients $r_i$ of the decomposition of the remainders in eq. (3.5). Indeed, these are expected to be much simpler than the master integral coefficients in the decomposition of the amplitude in eq. (3.3). To better understand the simplicity of the remainder in comparison to the two-loop amplitudes for five-point QCD amplitudes, in this paper we reconstruct the coefficients $c_{\Gamma,i}$ in eq. (3.3). We thus obtain a decomposition for the amplitude that is valid to all orders in $\epsilon$. This is the first time that such a decomposition has been made available for a two-loop five-point QCD amplitude. While this expression contains more information than required for the computation of NNLO corrections to the process in eq. (2.1), it can be used for defining the three-loop remainder required for N³LO corrections. For NNLO applications, we can use it to extract two-loop remainders in the form given in eq. (3.5).

Let us briefly discuss how the master-integral coefficients are computed. We proceed following the same steps as in ref. [49], which we adapt to the reconstruction of master-integral coefficients. For simplicity of the expressions, we use a single label to index the master integrals and rewrite eq. (3.3) as

$$A_h^{(2,j)} = \sum_i c_i(\epsilon, \vec{s}, \text{tr}_5) m_i(\epsilon, \vec{s}, \text{tr}_5), \quad (3.6)$$

<footnote n="3">³It is clear that the coefficients of the Laurent expansion of $I^{(1)}$ and $I^{(2)}$ around $\epsilon = 0$ also belong to this space of functions.</footnote>
where $\vec{s}$ denotes the set of five independent $s_{ij}$ defined in eq. (2.2) and the Levi-Civita contraction $\text{tr}_5$ is defined in eq. (2.3). We discuss the details of our choice of master-integral basis in appendix A. We then make the ansatz that

$$c_i(\epsilon, \vec{s}, \text{tr}_5) = \frac{1}{P_i(\epsilon)} \sum_{k=0}^{\kappa_i} \epsilon^k c_{i,k}(\vec{s}, \text{tr}_5), \qquad (3.7)$$

where we used the fact that there are no poles in $\epsilon$ that are kinematic dependent, and $\kappa_i$ is the maximal power of $\epsilon$ in the numerator. The polynomials $P_i(\epsilon)$ are trivial to determine by sampling the coefficients at enough values of $\epsilon$ (for details see e.g. [36]). Since $\text{tr}_5$ can be written as the square root of a polynomial in the Mandelstam variables, it follows that the most generic coefficient $c_{i,k}$ can be written as

$$c_{i,k}(\vec{s}, \text{tr}_5) = c_{i,k}^{+}(\vec{s}) + \text{tr}_5 c_{i,k}^{-}(\vec{s}), \qquad (3.8)$$

where $c_{i,k}^{+}(\vec{s})$ and $c_{i,k}^{-}(\vec{s})$ are rational functions of the $s_{ij}$. We note that

$$\begin{aligned} c_{i,k}^{+}(\vec{s}) &= \frac{1}{2} (c_{i,k}(\vec{s}, \text{tr}_5) + c_{i,k}(\vec{s}, -\text{tr}_5)), \\ c_{i,k}^{-}(\vec{s}) &= \frac{1}{2 \text{tr}_5} (c_{i,k}(\vec{s}, \text{tr}_5) - c_{i,k}(\vec{s}, -\text{tr}_5)), \end{aligned} \qquad (3.9)$$

which means that we can access individually each of these rational functions of the $s_{ij}$ by evaluating the amplitudes at parity-conjugate phase-space points.

To determine the analytic form of the rational functions $c_{i,k}^{\pm}$, we numerically evaluate the master integral coefficients using the two-loop numerical unitarity method outlined in section 3.1. We recall that these numerical evaluations are performed using finite-field arithmetic which allows us to obtain exact values for the coefficients [26]. We use the variables defined in appendix C of [49] to obtain a rational parametrization of phase space. This parametrization presents the major advantage of having all but one twistor variable being Mandelstam invariants, which removes any ambiguity in mapping an expression from twistor variables to Mandelstam invariants. Finally, in all calculations we set $s_{12} = 1$ and work with dimensionless variables. The dependence on $s_{12}$ is reintroduced by dimensional analysis.

As in ref. [49], we find that the denominator of the $c_{i,k}^{\pm}(\vec{s})$ can be easily determined. Indeed, we find that

$$c_{i,k}^{\pm}(\vec{s}) = \frac{n_{i,k}^{\pm}(\vec{s})}{\prod_j W_j^{q_{i,k,j}}(\vec{s})}, \qquad (3.10)$$

where the $W(\vec{s})$ are a subset of the so-called letters of the symbol alphabet associated with the contributing master integrals, namely the subset that is polynomial in $\vec{s}$ and not $\text{tr}_5$.<fnref n="4" /> As the photons are not ordered, it is not sufficient to consider the planar alphabet which is only closed under cyclic permutations. The full non-planar alphabet [65], however, is

<footnote n="4">⁴Denominators of the form $\text{tr}_5^n$ with $n$ even are allowed since they are polynomials in $\vec{s}$, see also the discussion in section 3.4.</footnote>
<table><caption>**Table 2.** Characterizing data for analytic expressions of master integrals, see eqs. (3.6) and (3.7). ‘Max degree’ is the highest polynomial degree across the different numerators $c_{i,k}^{\pm}$. The column ‘# independent $c_{i,k}^{\pm}$’ gives the dimension of the space of rational functions required to write all the $c_{i,k}^{\pm}$.</caption><thead><tr><th>Helicity</th><th>Max degree</th><th># independent c<sub>i,k</sub><sup>&pm;</sup></th></tr></thead><tbody><tr><td>A<sub>-++</sub><sup>(2,0)</sup></td><td>32</td><td>1320</td></tr><tr><td>A<sub>-++</sub><sup>(2,N<sub>f</sub>)</sup></td><td>20</td><td>203</td></tr><tr><td>A<sub>+++</sub><sup>(2,0)</sup></td><td>27</td><td>1244</td></tr><tr><td>A<sub>+++</sub><sup>(2,N<sub>f</sub>)</sup></td><td>18</td><td>130</td></tr></tbody></table>

closed under all photon permutations and sufficient to determine all the denominators. The powers $q_{i,k,j}$ are determined by evaluating the amplitudes on a one-dimensional line in phase-space [47, 49].

The determination of the analytic form of the coefficients $c_{i,k}$ is then reduced to the determination of the polynomials $n_{i,k}^{\pm}(\vec{s})$, which depend on four variables (once we set $s_{12} = 1$). With our in-house implementation of the multivariate Newton method [36], we obtain their analytic form with coefficients in a finite field. Following the simplification procedure of the expressions described in ref. [49], we were able to reconstruct all rational coefficients from a small number of numerical evaluations in additional finite fields. We provide these results in the supplementary material.

In table 2 we compile information that characterizes the complexity of the master integral coefficients. The column labeled ‘Max degree’ is a measure of the complexity of the functional reconstruction step. The column ‘# independent $c_{i,k}^{\pm}$’ is a measure of the complexity of the final result. We find that the complexity of the reconstruction is comparable to the two-loop remainders of five-parton amplitudes reconstructed in ref. [49], but the complexity of the final result is much higher than in the five-parton remainder case (where the number of independent rational functions was at most $O(120)$).

## 3.4 Two-loop finite remainders

Once the coefficients $c_i$ in eq. (3.6) have been determined, and the expressions for the master integrals in terms of pentagon functions are known, we can obtain the decomposition of the two-loop remainder in terms of pentagon functions, see eq. (3.5). The coefficients $r_i$ have a form similar to the $c_{i,k}$, that is

$$r_i(\vec{s}, \text{tr}_5) = r_i^+(\vec{s}) + \text{tr}_5 r_i^-(\vec{s}), \qquad (3.11)$$

where the $r_i^+(\vec{s})$ and $r_i^-(\vec{s})$ are rational functions of the $s_{ij}$. As for the $c_{i,k}^{\pm}$, the denominator of the $r_i^{\pm}$ is given by products of the letters in the subset of the (non-planar) symbol alphabet that are polynomial in the $s_{ij}$. It is interesting to note that the denominator $\text{tr}_5^2$ is absent, even though it is polynomial in the Mandelstam invariants. This is not the case for the $c_{i,k}^{\pm}$, that is there are poles in the coefficients $c_i(\epsilon, \vec{s}, \text{tr}_5)$ that are absent in the $r_i(\vec{s}, \text{tr}_5)$, and the $r_i$ are thus expected to be numerically more stable than the $c_i$.
<table><caption>**Table 3.** Characterizing data for analytic expressions of two-loop remainders, see eq. (2.19) and the discussion below for the notation. We use the same labels as in table 2. ‘Max weight’ is the highest transcendental weight of the pentagon functions appearing in each remainder.</caption><thead><tr><th>Helicity</th><th>Max degree</th><th># independent r<sub>i</sub><sup>&pm;</sup></th><th>Max weight</th></tr></thead><tbody><tr><td>R<sub>-++</sub><sup>(2,0)</sup></td><td>30</td><td>171</td><td>4</td></tr><tr><td>R<sub>-++</sub><sup>(2,N<sub>f</sub>)</sup></td><td>13</td><td>57</td><td>3</td></tr><tr><td>R<sub>+++</sub><sup>(2,0)</sup></td><td>16</td><td>62</td><td>2</td></tr><tr><td>R<sub>+++</sub><sup>(2,N<sub>f</sub>)</sup></td><td>12</td><td>12</td><td>1</td></tr></tbody></table>

We also note that, owing to the complexity of the rational functions in $c_{i,k}^{\pm}$, it is a non-trivial task to obtain the expressions for the $r_i$ from them. Finally, it is important to use the simplification procedure described in ref. [49] to obtain relatively compact expressions suitable for efficient numerical evaluation. Our expressions for the remainders are provided in the supplementary material.

In table 3 we compile some data to characterize the complexity of the analytic expressions. Compared to table 2, we added the column ‘Max weight’ which gives a measure of the complexity of the contributing pentagon functions.<fnref n="5" /> The highest possible transcendental weight in a two-loop remainder is 4, and only one remainder saturates this bound. This number is also important for numerical evaluations, as the highest weights dominate the evaluation time of the pentagon functions. It is interesting to note that the maximal polynomial degree of the numerator we need to reconstruct is only slightly lower than in the case of master integral coefficients, see table 2, which means that the complexity of the reconstruction of remainders and of master integral coefficients is not substantially different. The complexity of the final result for master integrals is however much higher, since they depend on a larger number of independent rational functions. As expected, for NNLO applications it is thus much more efficient to work at the level of the remainders.

## 3.5 Reference values

In order to facilitate the comparison with our results and to explicitly demonstrate the pole structure of the amplitudes, we present a numerical evaluation of the remainders and loop amplitudes on a randomly-chosen phase-space point in the physical region. Given that we are computing the Lorentz-invariant quantities defined in eq. (2.14), it is sufficient to specify the five independent Mandelstam variables together with the value of tr$_5$,

$$ \begin{aligned} s_{12} &= 1.322500000, & s_{23} &= -0.994109498, & s_{34} &= 0.264471591, \\ s_{45} &= 0.267126049, & s_{15} &= -0.883795230, & tr_5 &= -0.113828361. \end{aligned} \quad (3.12) $$

The one- and two-loop amplitudes and remainders evaluated on this point are presented in tables 4 to 6. Expressions for the amplitudes in terms of master integrals for the one-

<footnote n="5">⁵Roughly speaking, the transcendental weight corresponds to the number of iterated integrations in the definition of the pentagon functions.</footnote>
<table><caption>**Table 4.** Reference evaluations of all independent bare one-loop amplitudes on the phase-space point of eq. (3.12).</caption><thead><tr><th></th><th>ε<sup>-2</sup></th><th>ε<sup>-1</sup></th><th>ε<sup>0</sup></th><th>ε<sup>1</sup></th><th>ε<sup>2</sup></th></tr></thead><tbody><tr><td>A<sub>-++</sub><sup>(1)</sup></td><td>-1.000000000</td><td>-3.174284697<br>-3.141592654i</td><td>-3.437681197<br>-16.69077768i</td><td>-4.542364174<br>-48.29215997i</td><td>-28.34154945<br>-104.73071151i</td></tr><tr><td>A<sub>+++</sub><sup>(1)</sup></td><td>0</td><td>0</td><td>-122.4876141<br>-218.2099911i</td><td>-613.1620024<br>-1772.249665i</td><td>-1264.781477<br>-6727.583766i</td></tr></tbody></table>

<table><caption>**Table 5.** Reference evaluations of all independent bare two-loop amplitudes on the phase-space point of eq. (3.12).</caption><thead><tr><th></th><th>ε<sup>-4</sup></th><th>ε<sup>-3</sup></th><th>ε<sup>-2</sup></th><th>ε<sup>-1</sup></th><th>ε<sup>0</sup></th></tr></thead><tbody><tr><td>A<sub>-++</sub><sup>(2,0)</sup></td><td>0.500000000</td><td>2.257618031<br>+3.141592654i</td><td>-3.317245357<br>+20.90350063i</td><td>-55.54942686<br>+44.34772277i</td><td>-248.7699347<br>-87.79211669i</td></tr><tr><td>A<sub>-++</sub><sup>(2,N<sub>f</sub>)</sup></td><td>0</td><td>0.1666666667</td><td>1.335872677<br>+1.047197551i</td><td>4.646264515<br>+12.872514370i</td><td>10.33373680<br>+83.15472523i</td></tr><tr><td>A<sub>+++</sub><sup>(2,0)</sup></td><td>0</td><td>0</td><td>122.4876141<br>+218.2099911i</td><td>-132.6755953<br>+2049.613188i</td><td>-9927.845724<br>+3575.607623i</td></tr><tr><td>A<sub>+++</sub><sup>(2,N<sub>f</sub>)</sup></td><td>0</td><td>0</td><td>0</td><td>81.65840942<br>+145.4733274i</td><td>895.9475013<br>+2327.538099i</td></tr></tbody></table>

<table><caption>**Table 6.** Reference evaluations of all independent one-loop and two-loop remainders on the phase-space point of eq. (3.12).</caption><tbody><tr><td>R<sub>-++</sub><sup>(1)</sup></td><td>-5.281908761 - 6.718468192i</td></tr><tr><td>R<sub>+++</sub><sup>(1)</sup></td><td>-122.4876141 - 218.2099911i</td></tr><tr><td>R<sub>-++</sub><sup>(2,0)</sup></td><td>-17.93042514 - 84.48074943i</td></tr><tr><td>R<sub>-++</sub><sup>(2,N<sub>f</sub>)</sup></td><td>8.536235118 + 25.51694192i</td></tr><tr><td>R<sub>+++</sub><sup>(2,0)</sup></td><td>-2043.581205 - 3461.464426i</td></tr><tr><td>R<sub>+++</sub><sup>(2,N<sub>f</sub>)</sup></td><td>327.6279319 + 861.8112864i</td></tr></tbody></table>

and two-loop amplitudes are presented in the supplementary material in the directories `anc/oneLoopAmplitudes/` and `anc/twoLoopAmplitudes/` respectively. Expressions for the finite remainders in terms of pentagon functions are presented in a series of files in the directory `anc/remainders/`. We show how to assemble these files into the full amplitudes and remainders in `anc/example_assembly.m`, where, using an included numerical evaluation of the master integrals and pentagon functions, the expressions are combined to compute the amplitudes and remainders at the reference phase-space point (3.12) and reproduce the numbers in tables 4 to 6.
## 3.6 Validation

In this section we discuss the checks we performed on our setup to compute the two-loop amplitudes for three-photon production at hadron colliders.

We first computed analytic four- and five-point one-loop amplitudes. The four-point amplitudes were checked against the results of ref. [10], and we found full agreement with the quoted remainders. For the five-point amplitudes, we reproduced the results obtained from OpenLoops [66] to order $\epsilon^0$ and numerically verified all the relations in table 1. Finally, we checked that five-point one-loop amplitudes have the correct collinear behavior up to order $\epsilon^2$, i.e., that in these limits they are either regular or factorize into products of splitting amplitudes and four-point amplitudes.

At the two-loop level, we first recomputed the four-point amplitudes, reproducing the two-loop remainders given in ref. [10]. For the five-point amplitudes, we verified that they have the correct pole structure and numerically checked that they satisfy the relations in table 1. To further check the results at order $\epsilon^0$, we also verified that the amplitudes have the correct collinear behavior as described in section B. Finally, since ref. [15] does not include explicit results for the two-loop amplitudes our expressions cannot be directly compared to theirs. Nevertheless, our results were shown to lead to consistent predictions for the production of three photons at the LHC [16].

# 4 Numerical evaluation

**Squared finite remainders.** The contribution of two-loop helicity amplitudes to physical cross sections is constructed from the finite remainders (see e.g. [53]) — specifically the squared finite remainders, summed over helicity and color states. This object, which we denote $H$, admits a perturbative expansion in powers of the renormalized coupling $\alpha_s$,

$$H = H^{(0)} + \frac{\alpha_s}{2\pi} H^{(1)} + \left(\frac{\alpha_s}{2\pi}\right)^2 H^{(2)} + \mathcal{O}(\alpha_s^3), \quad (4.1)$$

which we normalize such that $H^{(0)} = 1$. The $\mathcal{O}(\alpha_s)$ contribution $H^{(1)}$ is given by

$$H^{(1)} = \frac{1}{\sum_h |\mathcal{M}_h^{(0)}|^2} C_F \left( \sum_h |\mathcal{M}_h^{(0)}|^2 2 \operatorname{Re}[R_h^{(1)}] \right), \quad (4.2)$$

and the $\mathcal{O}(\alpha_s^2)$ contribution $H^{(2)}$ by

$$H^{(2)} = \frac{1}{\sum_h |\mathcal{M}_h^{(0)}|^2} \left( \sum_h \frac{N_c^2}{4} N_c |\Phi_h R_h^{(1)}|^2 + \sum_h |\mathcal{M}_h^{(0)}|^2 2 \operatorname{Re} \left[ \frac{N_c^2}{4} R_h^{(2,0)} + C_F T_F N_f R_h^{(2,N_f)} \right] \right). \quad (4.3)$$

In this expression, the first line corresponds to the one-loop squared contribution, and the second line to the interference of two-loop and tree-level amplitudes. In the one-loop squared contributions, we have expanded the $C_F^2$ factor for consistency with the limit in
which the two-loop contributions were computed. The normalization factor in eq. (4.2)–(4.3) is consistent with setting $H^{(0)} = 1$, and in our normalization

$$ |\mathcal{M}_h^{(0)}|^2 = e_q^6 N_c |\Phi_h|^2 |A_h^{(0)}|^2. \qquad (4.4) $$

The helicity sums in eq. (4.2)–(4.3) are performed using the relations from table 1. We note that the remainders $R_{++}^{(2,0)}$ and $R_{+++}^{(2,N_f)}$ do not contribute to $H^{(2)}$, in the same way that $R_{+++}^{(1)}$ does not contribute to $H^{(1)}$. Nevertheless, $R_{+++}^{(1)}$ does contribute to $H^{(2)}$ through the one-loop squared contributions.

**Numerical evaluation.** Having phenomenological applications in mind, we have implemented in a C++ library the numerical evaluation of the finite remainders $R_h^{(1,0)}$, $R_h^{(2,0)}$, $R_h^{(2,N_f)}$ (see eq. (2.19)), and of $H^{(1)}$ and $H^{(2)}$ as defined in eq. (4.2)–(4.3). This library can be obtained from a git repository [52]. The library relies on PentagonFunctions++ [23] for the numerical evaluation of pentagon functions. For installation and usage instructions we refer to the README.md file which can be found in the root directory of the repository.

We recall that the rational coefficients in the decomposition of the remainders in terms of pentagon functions are simplified using the multivariate partial fractioning procedure outlined in [49]. Furthermore, we optimize the evaluation of rational coefficients using FORM [67, 68]. As a result, the time spent on their numerical evaluation is negligible compared to the time spent on the evaluation of transcendental functions. On average, the evaluation time of $H^{(2)}$ is on the order of a few seconds per phase-space point in double precision, using the default settings of PentagonFunctions++.

Besides the evaluation speed, the calculation of the $H^{(i)}$ must be numerically stable. To demonstrate the stability of our results, we compare the numerical evaluation of $H^{(2)}$ in double precision, which we denote $H_{\text{double}}^{(2)}$, with the evaluation in quadruple precision, which we denote $H_{\text{quad}}^{(2)}$, on a sample of 90000 phase-space points from the distribution employed in [16] for the computation of predictions at the center-of-mass energy of 8 TeV. The phase-space cuts correspond to the ones used in the ATLAS 8 TeV measurement [17] (we refer to the table 1 of [16] for the explicit definitions). We set $N_c = 3$, $N_f = 5$, and the renormalization scale $\mu_R$ to the invariant mass of the three-photon system $m_{\gamma\gamma\gamma}$. Assuming $H_{\text{quad}}^{(2)}$ to be correct at least up to a relative error of $\sim 10^{-16}$, we define

$$ d = -\log_{10} \left| \frac{H_{\text{double}}^{(2)}}{H_{\text{quad}}^{(2)}} - 1 \right| \qquad (4.5) $$

as a measure of the number of correct decimal digits in $H_{\text{double}}^{(2)}$. In figure 2 we show a histogram of this quantity for the 90000 sampled points on a logarithmic scale. We observe that less than 0.1% of the sampled points have an accuracy of less than four digits. This level of accuracy is more than adequate for phenomenological applications. Indeed, our implementation was already used for a Monte Carlo phase-space integration in [16], which converged to an overall integration error below 1% in NNLO differential distributions.

It is interesting to note that a good understanding of the physical properties that govern the analytic structure of scattering amplitudes can be used to explain and improve the
<figure><img src="image_6.png" /><figcaption>**Figure 2.** The logarithmic distribution of correct decimal digits (as defined in eq. (4.5)) for 90000 double-precision evaluations of the $H^{(2)}$ function. The phase-space points are sampled from a distribution representative of typical phenomenological studies.</figcaption></figure>

behavior of numerical algorithms. In particular, we recall that our analytic representation of the remainders (and thus of the $H^{(i)}$) only has poles that are associated with a subset of the letters of the two-loop five-point massless alphabet [65], see the discussion in section 3.4. While the letters corresponding to unphysical or spurious singularities can vanish inside the physical phase space, the amplitude must stay regular As an example, $s_{12} + s_{23}$ is a letter which can vanish in the physical phase space, see eq. (2.4), but is not a physical threshold. This implies that large cancellations can potentially occur when the amplitude is evaluated on phase-space points that are close to the surfaces where those letters vanish. Conversely, away from these small neighborhoods, the numerical evaluations should be accurate. We observe precisely this behavior. We verified that all of the unstable phase-space points in figure 2 are close to surfaces associated with spurious singularities. This shows that it is important to organize the analytic structure of the rational coefficients in the expressions for amplitudes or remainders using physical considerations, as it allows us to obtain compact expressions whose numerical evaluation is under full control. Indeed, once such a form is found, we can setup a robust precision-rescue system based on the fact that all the problematic regions of phase-space are explicitly known. We leave this for future work.

# 5 Conclusion

In this paper we have computed the two-loop planar corrections to the production of three photons at hadron colliders. This was achieved within the framework of two-loop numerical unitarity, coupled with analytic reconstruction techniques. Our results include expressions for both one- and two-loop amplitudes valid to all orders in the dimensional regulator.
In both cases, the amplitudes are written in terms of a set of master integrals. To our knowledge, this is the first time that master-integral coefficients have been obtained in analytic form for physical two-loop five-point scattering amplitudes. All our results are presented in the supplementary material.

As is well known, for NNLO (two-loop) phenomenology only the finite remainders of two-loop amplitudes are needed. By writing one- and two-loop amplitudes in terms of a basis of pentagon functions up to the required order in $\epsilon$, we obtained a decomposition of the remainder in terms of these functions. After simplifying the coefficients in this decomposition using multivariate partial-fractioning techniques, we obtained compact expressions for the two-loop remainders of the two independent helicity amplitudes. We demonstrated that the remainders are simpler functions than the all-order amplitudes.

While over the last years there has been a number of new results for massless five-point amplitudes, when having in mind phenomenological applications it is also important that the expressions for the remainders are numerically stable and can be efficiently evaluated across the relevant physical phase space. This point has been a major obstacle to computing the NNLO corrections for those processes. The expressions we obtain in this paper are ready to be used for phenomenological studies. This was demonstrated by verifying the numerical stability of the remainders when combined in a code that computes the color-and helicity-summed squared remainders. We have made this code public in a format that can be interfaced with real-radiation programs and employed to compute complete NNLO theoretical predictions. This is the first time that a two-loop massless five-point process is available in this form. Our results have already been used in ref. [16], and we expect that they will be instrumental in any future studies.

## Acknowledgments

We thank Fernando Febres Cordero for many enlightening discussions and assistance in setting up the cluster runs, and we thank Harald Ita for many discussions and comments on this manuscript. The work of B.P is supported by the French Agence Nationale pour la Recherche, under grant ANR-17-CE31-0001-01. V.S. is supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme, *Novel structures in scattering amplitudes* (grant agreement No. 725110). The authors acknowledge support by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant no INST 39/963-1 FUGG (bwForCluster NEMO).

## A Master integral basis

Here we give details on the bases of master integrals employed in the supplementary material. We consider a master integral decomposition at both one and two loops of the form in eq. (3.3), where the master integrals within the basis are labelled $m_{\Gamma,i}$. Each propagator structure $\Gamma$ can be understood diagrammatically, and they can be grouped into sets where the elements are related by different permutations of the external legs. We call each such
<table><caption>**Table 7.** Notation used in specifying the numerators of master integrals.</caption><tr><td>l[i]</td><td>l<sub>i</sub></td></tr><tr><td>p[k]</td><td>p<sub>k</sub></td></tr><tr><td>sq[v]</td><td>v<sup>2</sup></td></tr><tr><td>sp[v[i], v[j]]</td><td>v<sub>i</sub> &middot; v<sub>j</sub></td></tr><tr><td>tr5[1, 2, 3, 4]</td><td>tr<sub>5</sub></td></tr><tr><td>mu[i, j]</td><td>&mu;<sub>ij</sub></td></tr><tr><td>D</td><td>D</td></tr></table>

grouping a topology. We recall that the index $i$ in $m_{\Gamma,i}$ labels the different master integrals associated with a given $\Gamma$. In the supplementary material, master integrals are labelled by a topology name, the index $i$, and the relevant scattering kinematics in terms of five-point Mandelstam invariants. This information appears in the format

MI[topologyName, i, mandelstam1, mandelstam2, ...].

In this representation, the permutations act on the Mandelstam arguments and hence permutations of master integrals are codified as different Mandelstam invariants appearing in the final arguments.

Each $m_{\Gamma,i}$ denotes a numerator. These are polynomial in the loop momenta, and rational in external kinematic variables and in the dimensional regulator $D$. The notation used in specifying the numerators in the supplementary material is defined in table 7. There, we make use of the scalar product of the components of the loop momenta beyond 4-dimensions

$$ \mu_{ij} = \ell_i^{D-4} \cdot \ell_j^{D-4}. \qquad (A.1) $$

We collect the definitions of the master integrals employed in the one- and two-loop amplitudes in the files `anc/oneLoopAmplitudes/MasterIntegralDefinitions.m` and `anc/twoLoopAmplitudes/MasterIntegralDefinitions.m`, in the format

MI[...] -> num,

where `num` is the expression of the numerator. All master integrals are “unitarity compatible”, in that they have unit propagator powers, and we specify the numerator $m_{\Gamma,i}(\ell)$ for a single representative of each topology. We present the external-momenta conventions and loop-momenta routing for a representative permutation of the one-loop integrals in table 8, and of the two-loop integrals in tables 9 to 11. In all cases, we take all external momenta as outgoing and we omit the loop momentum labels when they are not required to disambiguate the numerators.
<figure><img src="image_7.png" /><figcaption>**Table 8.** Permutation representatives for one-loop master integral topologies. For one-loop master integrals we use the basis of [69], which includes only a certain subset of IBP relations (those that are independent of $D$). The basis thus includes integrals that would be reducible if one were to use the full set of IBP relations, such as all "OneMassTriangle" masters. The advantage of using this basis is that the coefficients in (3.3) are $D$-independent.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Table 9.** Permutation representatives for five-point two-loop master integral topologies.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Table 10. Permutation representatives for four-point two-loop master integral topologies.</figcaption></figure>

## B Collinear checks

The behavior of scattering amplitudes with massless external legs in the limit where two external legs become collinear is well understood. In this appendix we discuss how this was used to check our results for the amplitudes for the production of three photons at hadron colliders.

Let us first review the universal behavior of scattering amplitudes when two massless external legs with four-momenta $p_i$ and $p_j$ become collinear. We denote this limit as $i||j$. As they are massless, in the collinear limit the momenta become proportional. We introduce the massless four-momentum $p_P$, such that when $i||j$

$$p_i \rightarrow z p_P \quad \text{and} \quad p_j \rightarrow (1-z) p_P. \qquad (B.1)$$

The variable $z$ can be understood as the fraction of energy of $p_P$ which is contributed by $p_i$. We note that $z \to 0$ or 1 corresponds, respectively, to the limit where particles $i$ or $j$ go soft.

The study of the collinear behavior of bare helicity amplitudes is better formulated at the level of color-stripped amplitudes. In our case, we define

$$\bar{\mathcal{A}}(1_{q^{h_1}}, 2_{q^{h_2}}, 3_{\gamma^{h_3}}, 4_{\gamma^{h_4}}, 5_{\gamma^{h_5}}) := \Phi(1_{q^{h_1}}, 2_{q^{h_2}}, 3_{\gamma^{h_3}}, 4_{\gamma^{h_4}}, 5_{\gamma^{h_5}}) \mathcal{A}(1_{q^{h_1}}, 2_{q^{h_2}}, 3_{\gamma^{h_3}}, 4_{\gamma^{h_4}}, 5_{\gamma^{h_5}}), \quad (B.2)$$
<table><caption>**Table 11.** Permutation representatives for two- and three-point two-loop master topologies.</caption><tr><td>MI["TwoMassSlashedTriangle", 1, 4, s23, s45]</td><td><img src="image_0.png" alt="Diagram for MI["TwoMassSlashedTriangle", 1, 4, s23, s45]" /></td></tr><tr><td>MI["OneMassTriangleBubble", i, s45]</td><td><img src="image_1.png" alt="Diagram for MI["OneMassTriangleBubble", i, s45]" /></td></tr><tr><td>MI["TwoMassTriangleBubble", i, s12, s45]</td><td><img src="image_2.png" alt="Diagram for MI["TwoMassTriangleBubble", i, s12, s45]" /></td></tr><tr><td>MI["ThreePointFactorizedBubbleBubble", i, s23, s45]</td><td><img src="image_3.png" alt="Diagram for MI["ThreePointFactorizedBubbleBubble", i, s23, s45]" /></td></tr><tr><td>MI["FactorizedBubbleBubble", i, s45]</td><td><img src="image_4.png" alt="Diagram for MI["FactorizedBubbleBubble", i, s45]" /></td></tr><tr><td>MI["Sunrise", i, s34]</td><td><img src="image_5.png" alt="Diagram for MI["Sunrise", i, s34]" /></td></tr></table>

which is trivially related to the $M$ and $A$ we have used elsewhere in the paper, see eq. (2.5). In particular, we will use the same abbreviated notational style whenever it is unambiguous. In the $i||j$ limit, the amplitudes can diverge. The divergent behavior is captured by the factorization formula

$$ \bar{\mathcal{A}}(i^{h_i}, j^{h_j}, \dots) \xrightarrow{i||j} \sum_h \mathrm{Split}_{Ph}(i^{h_i}, j^{h_j}) \bar{\mathcal{A}}(\bar{P}^{-h}, \dots). \qquad (B.3) $$

Here, the factors $\mathrm{Split}_{Ph}(i^{h_i}, j^{h_j})$ are known as the splitting amplitudes for the particle $P$ with heliciy $h$ to split into particles $i$ with heliciy $h_i$ and $j$ with helicity $h_j$. The splitting functions depend only on the momenta, helicity and color quantum numbers of particles $i$, $j$ and $P$, and the dependence on $p_P$ is solely in terms of the energy fraction $z$. If the amplitude on the left-hand side has $n$ external legs, the amplitude on the right-hand side is a $(n-1)$-point amplitude, where $\bar{P}$ denotes the conjugate particle to the particle $P$ in the subscript of the splitting function. Finally, we note that it might be that the right-hand side of eq. (B.3) is zero, in which case the amplitude on the left-hand side is regular in the $i||j$ limit.

The behavior in eq. (B.3) is valid to all orders in the coupling constant. Naturally, the splitting amplitudes have a perturbative expansion in powers of $\alpha_s^0$,

$$ \mathrm{Split}_{Ph}(i^{h_i}, j^{h_j}) = \mathrm{Split}_{Ph}^{(0)}(i^{h_i}, j^{h_j}) + \frac{\alpha_s^0}{2\pi} \mathrm{Split}_{Ph}^{(1)}(i^{h_i}, j^{h_j}) \\ + \left(\frac{\alpha_s^0}{2\pi}\right)^2 \mathrm{Split}_{Ph}^{(2)}(i^{h_i}, j^{h_j}) + \mathcal{O}\left((\alpha_s^0)^3\right). \qquad (B.4) $$
<figure><img src="image_10.png" /><figcaption>**Figure 3.** Graphical representation of equation (B.5). White blobs represent the splitting amplitudes and hatched blobs represent the four-point amplitudes. The number of contained circles represents the loop order of the object. The types of particles are omitted from this representation, the external and internal lines representing only momenta.</figcaption></figure>

We can then expand both sides of (B.3), and we are particularly interested in the order $(\alpha_s^0)^2$ contributions:

$$ \bar{\mathcal{A}}^{(2)}(i^{h_i}, j^{h_j}, \dots) \xrightarrow{i \parallel j} \sum_h \left( \begin{aligned} & \mathrm{Split}_{\mathrm{ph}}^{(2)}(i^{h_i}, j^{h_j}) \bar{\mathcal{A}}^{(0)}(\bar{P}^{-h}, \dots) + \mathrm{Split}_{\mathrm{ph}}^{(1)}(i^{h_i}, j^{h_j}) \bar{\mathcal{A}}^{(1)}(\bar{P}^{-h}, \dots) \\ & + \mathrm{Split}_{\mathrm{ph}}^{(0)}(i^{h_i}, j^{h_j}) \bar{\mathcal{A}}^{(2)}(\bar{P}^{-h}, \dots) \end{aligned} \right). \qquad (B.5) $$

It can be useful to understand this equation graphically, as presented in figure 3. As in eq. (B.3), if the right-hand side vanishes, then the left-hand side is regular in the $i \parallel j$ limit.

In the context of the three-photon production amplitudes at hand, many of the collinear limits are equivalent due to the Bose symmetry of the final state. As the process is unordered, the limits can be categorized according to the particles which become collinear — either $(\gamma, \gamma)$, $(q, \bar{q})$, $(q, \gamma)$ or $(\bar{q}, \gamma)$. These can then be further classified according to the helicity of the involved particles. Our approach to check that our results satisfy eq. (B.5) is to evaluate both sides of the equation numerically. Note that this also tests the pentagon functions in near-singular regions of phase-space. To check their numerical behaviour, we constructed differential equations for the pentagon functions from their analytic representation, which we then solved numerically to high precision using [70]. In the following, we discuss each of the different collinear configurations. Before doing so, we note that we have performed the same checks on the one-loop five-point amplitudes up to order $\epsilon^2$. This allows one to check the implementation of the one-loop splitting amplitudes up to $O(\epsilon^2)$ (see section B.1 below), and is also a non-trivial check of the $O(\epsilon)$ and $O(\epsilon^2)$ contributions to the one-loop amplitudes which we have computed.

($q, \bar{q}$) limit: this limit is not accessible in the physical region for three-photon production at hadron colliders, see eq. (2.4). For this reason equation eq. (B.5) cannot be checked in the $(q, \bar{q})$ limit.

($\gamma, \gamma$) limit: the only splitting amplitude that can contribute is $\gamma^{h_1} \rightarrow \gamma^{h_2} \gamma^{h_3}$ which vanishes to all orders for any helicity. It then follows that the amplitudes should be regular when any of the photons become collinear, independently of the helicity. We have verified that this is indeed the case for our amplitudes.

($q, \gamma$) limit: for the amplitude $\bar{\mathcal{A}}_{+++}$, given the fact that all the photons have the same helicity, there is a single limit to check corresponding to the splitting $q \rightarrow \gamma^+q$. We note that the four-point amplitudes on the right-hand side of eq. (B.5) would then have a vanishing tree, and the checks thus only involve the tree-level and one-loop splitting amplitudes. For
the amplitude $\bar{\mathcal{A}}_{-++}$, there are two different limits to check, $q \to \gamma^-q$ and $q \to \gamma^+q$. In the first case, only tree-level and one-loop splitting amplitudes are needed. In the latter case, two-loop splitting amplitudes also contribute. For the two cases that only require up to one-loop splitting amplitudes, we checked eq. (B.5) numerically, using the expressions given in section B.1 below. For the case that involves the two-loop splitting amplitudes, we performed a different type of check. We first note that loop-level splitting amplitudes, when normalized to the corresponding tree-level, depend only on the vanishing $s_{ij}$ and the energy fraction $z$ (see e.g. eqs. (B.10) and (B.9) below). We can then consider two different collinear phase-space points which have the same values of $s_{ij}$ and $z$ but with all other kinematic parameters different. By evaluating our amplitudes on these two points we can verify that the amplitudes are consistent with the expected dependence of the normalized two-loop splitting amplitude.

($\bar{q}, \gamma$) limit: the situation is exactly the same as for the $(q, \gamma)$ limit.

## B.1 One-loop photon splitting amplitudes

We write the splitting amplitudes as

$$ \text{Split}_{P_h}^{(l)}(i^{h_i}, j^{h_j}) = \rho^{(l)}(i^{h_i}, j^{h_j})\text{Split}^{(0)}(i^{h_i}, j^{h_j}). \qquad (B.6) $$

The $\rho^{(l)}(i^{h_i}, j^{h_j})$ are transcendental functions of $s_{ij}$ and $z$, see eq. (B.1). Furthermore, they depend on the scheme-defining parameter $D_s$, which is equal to $D$ in the HV scheme but is kept generic here. Most of our checks only depend on the one-loop splitting amplitudes, which can be obtained from the one-loop QCD splitting amplitudes [71] with color-algebra manipulations. We first define

$$ r_{\Gamma}(\epsilon) = \frac{1}{(4\pi)^{\epsilon}} \frac{\Gamma(1-\epsilon)^2 \Gamma(1+\epsilon)}{\Gamma(1-2\epsilon)}, \qquad (B.7) $$

and

$$ F(z, s_{ij}, \epsilon) = r_{\Gamma}(\epsilon) \left(\frac{-s_{ij}}{\mu^2}\right)^{-\epsilon} \frac{1}{\epsilon^2} \sum_{m=1}^{\infty} \epsilon^m \text{Li}_m \left(\frac{-z}{1-z}\right). \qquad (B.8) $$

We then find that the “distinct-helicity” splitting amplitudes are

$$ \rho_{q^+}^{(1)}(i_{\gamma}^{+}, j_{q}^{-}) = \rho_{q^{-}}^{(1)}(i_{\gamma}^{-}, j_{q}^{+}) = F(z, s_{ij}, \epsilon), \qquad (B.9) $$

and the “like-helicity” splitting amplitudes are

$$ \rho_{q^+}^{(1)}(i_{\gamma}^{-}, j_{q}^{-}) = \rho_{q^{-}}^{(1)}(i_{\gamma}^{+}, j_{q}^{+}) = F(z, s_{ij}, \epsilon) - z r_{\Gamma}(\epsilon) \left(\frac{-s_{ij}}{\mu^2}\right)^{-\epsilon} \frac{D_s - 2}{4(1-\epsilon)(1-2\epsilon)}. \qquad (B.10) $$

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.
References

[1] S. Amoroso et al., *Les Houches 2019: physics at TeV colliders: standard model working group report*, in the proceedings of the 11th Les Houches Workshop on Physics at TeV Colliders: PhysTeV Les Houches, June 10-28, Les Houches, France (2019) [arXiv:2003.01700] [INSPIRE].

[2] ATLAS collaboration, *Search for new phenomena in events with at least three photons collected in pp collisions at $\sqrt{s}$ = 8 TeV with the ATLAS detector*, Eur. Phys. J. C 76 (2016) 210 [arXiv:1509.05051] [INSPIRE].

[3] H. Denizli, K.Y. Oyulmaz and A. Senol, *Testing for observability of Higgs effective couplings in triphoton production at FCC-hh*, J. Phys. G 46 (2019) 105007 [arXiv:1901.04784] [INSPIRE].

[4] H. Denizli and A. Senol, *Sensitivity reach on anomalous Higgs couplings via triphoton production for the post-LHC circular high-energy hadron colliders*, arXiv:2005.08760 [INSPIRE].

[5] J.A. Aguilar-Saavedra, J.M. Cano and J.M. No, *More light on Higgs flavor at the LHC: Higgs couplings to light quarks through $h + \gamma$ production*, arXiv:2008.12538 [INSPIRE].

[6] A. Zerwekh, C. Dib and R. Rosenfeld, *Triple photon production at the Tevatron in technicolor models*, Phys. Lett. B 549 (2002) 154 [hep-ph/0207270] [INSPIRE].

[7] N. Toro and I. Yavin, *Multiphotons and photon jets from new heavy vector bosons*, Phys. Rev. D 86 (2012) 055005 [arXiv:1202.6377] [INSPIRE].

[8] G. Das and P. Mathews, *Neutral triple vector boson production in Randall-Sundrum model at the LHC*, Phys. Rev. D 92 (2015) 094034 [arXiv:1507.08857] [INSPIRE].

[9] C. Anastasiou, E.W. Glover and M.E. Tejeda-Yeomans, *Two loop QED and QCD corrections to massless fermion boson scattering*, Nucl. Phys. B 629 (2002) 255 [hep-ph/0201274] [INSPIRE].

[10] E.W. Glover and M.E. Tejeda-Yeomans, *Two loop QCD helicity amplitudes for massless quark massless gauge boson scattering*, JHEP 06 (2003) 033 [hep-ph/0304169] [INSPIRE].

[11] S. Catani, L. Cieri, D. de Florian, G. Ferrera and M. Grazzini, *Diphoton production at hadron colliders: a fully-differential QCD calculation at NNLO*, Phys. Rev. Lett. 108 (2012) 072001 [Erratum ibid. 117 (2016) 089901] [arXiv:1110.2375] [INSPIRE].

[12] J.M. Campbell, R. Ellis, Y. Li and C. Williams, *Predictions for diphoton production at the LHC through NNLO in QCD*, JHEP 07 (2016) 148 [arXiv:1603.02663] [INSPIRE].

[13] T. Gehrmann, N. Glover, A. Huss and J. Whitehead, *Scale and isolation sensitivity of diphoton distributions at the LHC*, arXiv:2009.11310 [INSPIRE].

[14] S. Alioli et al., *Precise predictions for photon pair production matched to parton showers in GENEVA*, arXiv:2010.10498 [INSPIRE].

[15] H.A. Chawdhry, M.L. Czakon, A. Mitov and R. Poncelet, *NNLO QCD corrections to three-photon production at the LHC*, JHEP 02 (2020) 057 [arXiv:1911.00479] [INSPIRE].

[16] S. Kallweit, V. Sotnikov and M. Wiesemann, *Triphoton production at hadron colliders in NNLO QCD*, Phys. Lett. B 812 (2021) 136013 [arXiv:2010.04681] [INSPIRE].
[17] ATLAS collaboration, Measurement of the production cross section of three isolated photons in pp collisions at $\sqrt{s} = 8$ TeV using the ATLAS detector, Phys. Lett. B 781 (2018) 55 [arXiv:1712.07291] [INSPIRE].

[18] C.G. Papadopoulos, D. Tommasini and C. Wever, The pentabox master integrals with the simplified differential equations approach, JHEP 04 (2016) 078 [arXiv:1511.09404] [INSPIRE].

[19] T. Gehrmann, J.M. Henn and N.A. Lo Presti, Analytic form of the two-loop planar five-gluon all-plus-helicity amplitude in QCD, Phys. Rev. Lett. 116 (2016) 062001 [Erratum ibid. 116 (2016) 189903] [arXiv:1511.05409] [INSPIRE].

[20] S. Abreu, L.J. Dixon, E. Herrmann, B. Page and M. Zeng, The two-loop five-point amplitude in $\mathcal{N} = 4$ super-Yang-Mills theory, Phys. Rev. Lett. 122 (2019) 121603 [arXiv:1812.08941] [INSPIRE].

[21] D. Chicherin, T. Gehrmann, J.M. Henn, P. Wasser, Y. Zhang and S. Zoia, All master integrals for three-jet production at next-to-next-to-leading order, Phys. Rev. Lett. 123 (2019) 041603 [arXiv:1812.11160] [INSPIRE].

[22] T. Gehrmann, J.M. Henn and N.A. Lo Presti, Pentagon functions for massless planar scattering amplitudes, JHEP 10 (2018) 103 [arXiv:1807.09812] [INSPIRE].

[23] D. Chicherin and V. Sotnikov, Pentagon functions for scattering of five massless particles, JHEP 20 (2020) 167 [arXiv:2009.07803] [INSPIRE].

[24] K.J. Larsen and Y. Zhang, Integration-by-parts reductions from unitarity cuts and algebraic geometry, Phys. Rev. D 93 (2016) 041701 [arXiv:1511.01071] [INSPIRE].

[25] H. Ita, Two-loop integrand decomposition into master integrals and surface terms, Phys. Rev. D 94 (2016) 116015 [arXiv:1510.05626] [INSPIRE].

[26] T. Peraro, Scattering amplitudes over finite fields and multivariate functional reconstruction, JHEP 12 (2016) 030 [arXiv:1608.01902] [INSPIRE].

[27] A. Georgoudis, K.J. Larsen and Y. Zhang, Azurite: an algebraic geometry based package for finding bases of loop integrals, Comput. Phys. Commun. 221 (2017) 203 [arXiv:1612.04252] [INSPIRE].

[28] A. Georgoudis, K.J. Larsen and Y. Zhang, Cristal and Azurite: new tools for integration-by-parts reductions, PoS(RADCOR2017)020 [arXiv:1712.07510] [INSPIRE].

[29] H.A. Chawdhry, M.A. Lim and A. Mitov, Two-loop five-point massless QCD amplitudes within the integration-by-parts approach, Phys. Rev. D 99 (2019) 076011 [arXiv:1805.09182] [INSPIRE].

[30] J. Böhm, A. Georgoudis, K.J. Larsen, H. Schönemann and Y. Zhang, Complete integration-by-parts reductions of the non-planar hexagon-box via module intersections, JHEP 09 (2018) 024 [arXiv:1805.01873] [INSPIRE].

[31] T. Peraro, FiniteFlow: multivariate functional reconstruction using finite fields and dataflow graphs, JHEP 07 (2019) 031 [arXiv:1905.08019] [INSPIRE].

[32] T. Peraro and L. Tancredi, Physical projectors for multi-leg helicity amplitudes, JHEP 07 (2019) 114 [arXiv:1906.03298] [INSPIRE].

[33] D. Bendle et al., Integration-by-parts reductions of Feynman integrals using Singular and GPI-space, JHEP 02 (2020) 079 [arXiv:1908.04301] [INSPIRE].
[34] Y. Wang, Z. Li and N. Ul Basat, *Direct reduction of multiloop multiscale scattering amplitudes*, Phys. Rev. D 101 (2020) 076023 [arXiv:1901.09390] [INSPIRE].

[35] X. Guan, X. Liu and Y.-Q. Ma, *Complete reduction of integrals in two-loop five-light-parton scattering amplitudes*, Chin. Phys. C 44 (2020) 093106 [arXiv:1912.09294] [INSPIRE].

[36] S. Abreu et al., *Caravel: a C++ framework for the computation of multi-loop amplitudes with numerical unitarity*, arXiv:2009.11957 [INSPIRE].

[37] A. von Manteuffel and R.M. Schabinger, *A novel approach to integration by parts reduction*, Phys. Lett. B 744 (2015) 101 [arXiv:1406.4513] [INSPIRE].

[38] J. Klappert and F. Lange, *Reconstructing rational functions with FireFly*, Comput. Phys. Commun. 247 (2020) 106951 [arXiv:1904.00009] [INSPIRE].

[39] J. Klappert, S.Y. Klein and F. Lange, *Interpolation of dense and sparse rational functions and other improvements in FireFly*, arXiv:2004.01463 [INSPIRE].

[40] J. Klappert, F. Lange, P. Maierhöfer and J. Usovitsch, *Integral reduction with Kira 2.0 and finite field methods*, arXiv:2008.06494 [INSPIRE].

[41] S. Badger, H. Frellesvig and Y. Zhang, *A two-loop five-gluon helicity amplitude in QCD*, JHEP 12 (2013) 045 [arXiv:1310.1051] [INSPIRE].

[42] S. Badger, C. Brønnum-Hansen, H.B. Hartanto and T. Peraro, *First look at two-loop five-gluon scattering in QCD*, Phys. Rev. Lett. 120 (2018) 092001 [arXiv:1712.02229] [INSPIRE].

[43] S. Abreu, F. Febres Cordero, H. Ita, B. Page and M. Zeng, *Planar two-loop five-gluon amplitudes from numerical unitarity*, Phys. Rev. D 97 (2018) 116014 [arXiv:1712.03946] [INSPIRE].

[44] S. Badger et al., *Applications of integrand reduction to two-loop five-point scattering amplitudes in QCD*, PoS LL2018 (2018) 006 [arXiv:1807.09709] [INSPIRE].

[45] S. Abreu, F. Febres Cordero, H. Ita, B. Page and V. Sotnikov, *Planar two-loop five-parton amplitudes from numerical unitarity*, JHEP 11 (2018) 116 [arXiv:1809.09067] [INSPIRE].

[46] S. Badger, C. Brønnum-Hansen, H.B. Hartanto and T. Peraro, *Analytic helicity amplitudes for two-loop five-gluon scattering: the single-minus case*, JHEP 01 (2019) 186 [arXiv:1811.11699] [INSPIRE].

[47] S. Abreu, J. Dormans, F. Febres Cordero, H. Ita and B. Page, *Analytic Form of Planar Two-Loop Five-Gluon Scattering Amplitudes in QCD*, Phys. Rev. Lett. 122 (2019) 082002 [arXiv:1812.04586] [INSPIRE].

[48] S. Badger et al., *Analytic form of the full two-loop five-gluon all-plus helicity amplitude*, Phys. Rev. Lett. 123 (2019) 071601 [arXiv:1905.03733] [INSPIRE].

[49] S. Abreu, J. Dormans, F. Febres Cordero, H. Ita, B. Page and V. Sotnikov, *Analytic form of the planar two-loop five-parton scattering amplitudes in QCD*, JHEP 05 (2019) 084 [arXiv:1904.00945] [INSPIRE].

[50] M. Grazzini, S. Kallweit and M. Wiesemann, *Fully differential NNLO computations with MATRIX*, Eur. Phys. J. C 78 (2018) 537 [arXiv:1711.06631] [INSPIRE].

[51] S. Abreu, F. Febres Cordero, H. Ita, M. Jaquier, B. Page and M. Zeng, *Two-loop four-gluon amplitudes from numerical unitarity*, Phys. Rev. Lett. 119 (2017) 142001 [arXiv:1703.05273] [INSPIRE].
[52] https://gitlab.com/five-point-amplitudes/FivePointAmplitudes-cpp.git

[53] S. Weinzierl, Does one need the $O(\epsilon)$- and $O(\epsilon^2)$-terms of one-loop amplitudes in an NNLO calculation?, Phys. Rev. D 84 (2011) 074007 [arXiv:1107.5131] [INSPIRE].

[54] S. Catani, The singular behavior of QCD amplitudes at two loop order, Phys. Lett. B 427 (1998) 161 [hep-ph/9802439] [INSPIRE].

[55] S. Abreu, F. Febres Cordero, H. Ita, M. Jaquier and B. Page, Subleading poles in the numerical unitarity method at two loops, Phys. Rev. D 95 (2017) 096011 [arXiv:1703.05255] [INSPIRE].

[56] A. Ochirov and B. Page, Full colour for loop amplitudes in Yang-mills theory, JHEP 02 (2017) 100 [arXiv:1612.04366] [INSPIRE].

[57] A. Ochirov and B. Page, Multi-quark colour decompositions from unitarity, JHEP 10 (2019) 058 [arXiv:1908.02695] [INSPIRE].

[58] F.R. Anger and V. Sotnikov, On the dimensional regularization of QCD helicity amplitudes with quarks, arXiv:1803.11127 [INSPIRE].

[59] V. Sotnikov, Scattering amplitudes with the multi-loop numerical unitarity method, Ph.D. thesis, Freiburg University, Freiburg, Germany (2019).

[60] F.A. Berends and W.T. Giele, Recursive calculations for processes with n gluons, Nucl. Phys. B 306 (1988) 759 [INSPIRE].

[61] A. Hodges, Eliminating spurious poles from gauge-theoretic amplitudes, JHEP 05 (2013) 135 [arXiv:0905.1473] [INSPIRE].

[62] A.B. Goncharov, M. Spradlin, C. Vergu and A. Volovich, Classical polylogarithms for amplitudes and Wilson loops, Phys. Rev. Lett. 105 (2010) 151605 [arXiv:1006.5703] [INSPIRE].

[63] C. Duhr, H. Gangl and J.R. Rhodes, From polygons and symbols to polylogarithmic functions, JHEP 10 (2012) 075 [arXiv:1110.0458] [INSPIRE].

[64] C. Duhr, Hopf algebras, coproducts and symbols: an application to Higgs boson amplitudes, JHEP 08 (2012) 043 [arXiv:1203.0454] [INSPIRE].

[65] D. Chicherin, J. Henn and V. Mitev, Bootstrapping pentagon functions, JHEP 05 (2018) 164 [arXiv:1712.09610] [INSPIRE].

[66] F. Buccioni et al., OpenLoops 2, Eur. Phys. J. C 79 (2019) 866 [arXiv:1907.13071] [INSPIRE].

[67] B. Ruijl, T. Ueda and J. Vermaseren, FORM version 4.2, arXiv:1707.06453 [INSPIRE].

[68] J. Kuipers, T. Ueda and J.A.M. Vermaseren, Code optimization in FORM, Comput. Phys. Commun. 189 (2015) 1 [arXiv:1310.7007] [INSPIRE].

[69] W.T. Giele, Z. Kunszt and K. Melnikov, Full one-loop amplitudes from tree amplitudes, JHEP 04 (2008) 049 [arXiv:0801.2237] [INSPIRE].

[70] M. Hidding, DiffExp, a Mathematica package for computing Feynman integrals in terms of one-dimensional series expansions, arXiv:2006.05510 [INSPIRE].

[71] 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] [INSPIRE].