Higgs boson decay into gluons in a 4D regularization: IR cancellation without evanescent fields to NLO

Abstract

Higgs decay using an effective Higgs–Yang-Mills interaction in terms of a dimension five operator as well as usual QCD interactions is revisited in the context of Implicit Regularization (IReg) and compared with conventional dimensional regularization (CDR), four dimensional helicity (FDH) and dimensional reduction (DRED) schemes. The decay rate for $H \rightarrow gg(g)$ is calculated in this strictly four-dimensional set-up to $\alpha _s^3$ order in the strong coupling. Moreover we include joint processes that contribute at the same perturbative order in the real emission channels consisting of 3 gluons as well as gluon quark–antiquark final states with light (zero mass) quarks. Unambiguous identification and separation of UV from IR divergences is achieved putting at work the renormalization group scale relation inherent to the method. UV singularities are removed as usual by renormalization, the IR divergences are cancelled due to the method’s compliance with the Kinoshita–Lee–Nauenberg (KLN) theorem. Most importantly, we verify that no evanescent fields such as $\epsilon$-scalars need be introduced as required by some mixed regularizations that operate partially in the physical dimension.

1 Introduction

The Standard Model of Particle Physics (SM) is a good description of the physics at and below the electroweak scale, the discovery of the Higgs boson at the Large Hadron Collider (LHC) being the confirmation of this framework [1, 2]. It is also clear that the SM does not provide a complete description of particle interactions. Phenomena such as dark matter and dark energy, a consistent quantum theory of gravity, the stability of the electroweak vacuum up to the Planck scale and the hierarchy problem, to name a few, motivated the development of models beyond the SM (BSM).

The Future Circular Collider aims at reaching collision energies of around 100 TeV demanding higher precision in theoretical computations for Standard Model phenomena and electroweak pseudo-observables (EWPOs). New calculational methods have been developed through computer algebraic algorithms for analytical and numerical methods. At a given loop order, the complexity is measured by the number of virtual massive particles in the evaluation of Feynman integrals. On the other hand, a theoretical library for perturbation theory calculations of Feynman integrals in closed form beyond one loop does not exist [3]. While numerical integration methods seem to be the main tool to address those problems, analytical techniques play an important role to systematically and unambiguously remove infrared and ultraviolet divergences from physical observables.

Regularization frameworks that operate partially or entirely in the physical dimension can bring some advantages and simplifications in the evaluation of Feynman amplitudes, whereupon extensions such as the dimensional reduction (DRED), the Four Dimensional Helicity scheme (FDH), and the Implicit Regularization (IReg), among others [4] have been constructed:

• CDR (Conventional Dimensional Regularization): internal and external gluons are all treated as d-dimensional.

• HV (’t Hooft Veltman scheme): internal gluons are d-dimensional and external ones are strictly 4-dimensional.

• DRED (Dimensional Reduction): internal and external gluons are all treated as quasi-4-dimensional.

• FDH (Four-Dimensional Helicity): internal gluons are treated as quasi 4-dimensional and external ones are treated as strictly 4-dimensional.

• IReg (Implicit Regularization): all fields as well as internal momenta are defined in the physical dimension.

Therefore IReg acts directly on the dimension of the theory and can be systematically implemented to all orders in perturbation theory [5, 6], in consonance with Bogoliubov’s recursion formula. From the point of view of applications, IReg has been used to obtain the two-loop gauge coupling beta function of abelian and non-abelian theories [7, 8], also in the context of supersymmetry [9,10,11]. In all cases, compliance with gauge (and supersymmetry) was established, with an all-order proof for abelian gauge symmetry already developed [12, 13].

CDR relies on the assumption that the quantities of interest depend smoothly on the spacetime dimensionality. Such an assumption leads to well known problems for dimension specific theories, such as supersymmetric and chiral theories, as discussed for example in the review [14], covering issues related to the $\gamma _5$ matrix¹. The history of dimensional methods and schemes to handle the chiral anomaly and preserve gauge invariance and renormalizability in the Standard Model is vast and to date in frank development, see [18,19,20] for recent works. Mostly this is effected through the construction of appropriate counterterms (CT) imposed by Ward–Takahashi and Slavnov–Taylor identities order by order in perturbation theory, which can be carried out consistently in the Breitenlohner – Maison scheme [21]. A method that dispenses CT, as it preserves the BRST symmetry at all loop-orders has been advocated in [22, 23]

In mixed regularization schemes that operate partially in the physical dimension such as DRED or FDH, an auxiliary space which has the characteristics of a 4-dimensional space, Q4S, is introduced. Such a quasi-4-dimensional space is decomposed as $Q4S=QdS\oplus Q n_{\epsilon } S$, where QdS is formally d-dimensional and its complement $Q n_{\epsilon } S$ has dimension $n_{\epsilon }=4-d$ [24]. The metric tensor for the original 4-dimensional space 4S is denoted by $\bar{g}^\mu _\nu$ whereas the metric tensor of the spaces Q4S, QdS and $Q n_{\epsilon } S$ are respectively written as $g^\mu _\nu$, $\hat{g}^\mu _\nu$ and $\tilde{g}^\mu _\nu$. They satisfy

$$\begin{aligned} g^\mu _\nu= & {} \hat{g}^\mu _\nu +\tilde{g}^\mu _\nu , \,\, g^{\mu \nu }\tilde{g}^\rho _\nu =\tilde{g}^{\mu \rho }, \,\,\nonumber \\ g^{\mu \nu }\hat{g}^\rho _\nu= & {} \hat{g}^{\mu \rho }, \,\, \hat{g}^{\mu \rho }\bar{g}_\rho ^\nu = \bar{g}^{\mu \nu }, \,\, \hat{g}^{\mu \nu }\tilde{g}^\rho _\nu = 0 , \end{aligned}$$

(1.1)

with

$$\begin{aligned} g^\mu _\mu =4, \, \tilde{g}^\mu _\mu =n_{\epsilon }=2\epsilon \,\,\, \text {and} \,\,\, \hat{g}^\mu _\mu =d. \end{aligned}$$

(1.2)

Furthermore, mathematical consistency and d-dimensional gauge invariance require that $Q4S \supset QDS \supset 4S$ and forbid to identify $g^{\mu \nu }$ as $\bar{g} ^{\mu \nu }$ [25].

Due to the decomposition of Q4S, the gauge field is aligned as $A^a_\mu = \hat{A} ^a _\mu + \epsilon ^a _\mu$, where $\hat{A}^a_\mu \in QdS$ and $\epsilon ^a_\mu \in Q n_{\epsilon } S$. The $\epsilon$-dimensional field $\epsilon ^a_\mu$ is a scalar under d-dimensional Lorentz transformations and transforms in the adjoint representation. Due to this decomposition, the Lagrangian is modified to incorporate $\epsilon$-scalars, giving rise to extra Feynman diagrams with $\epsilon$-scalars, which contribute additional finite terms to divergent loop amplitudes [26].²

Ideally a fully mathematical consistent regularization scheme that prevents the emergence of symmetry breaking terms or spurious anomalies and that is valid to arbitrary higher order is necessary. IReg is advantageous in the sense that gauge symmetry breaking terms are linked to momentum routing violation in Feynman diagram loops. Such symmetry breaking terms accompany surface terms whose structure is built to arbitrary loop order, in a regularization independent fashion. Ultraviolet renormalization is constructed by subtracting loop integrals which need not be explicitly evaluated to define renormalization group functions. On the other hand an invariant regularization scheme must comply with infrared finiteness as stated by the Kinoshita–Lee–Nauenberg theorem [28, 29]. The KLN theorem, in a nutshell, states that in gauge theories the infrared divergences coming from loop integrals are canceled by the IR divergences coming from phase space integrals. As a result it is of theoretical interest to test the applicability of IReg in a practical calculation involving IR divergences that only cancel at the level of cross sections or decay rates, verifying the KLN theorem.

The $H \rightarrow gg$ decay described by an effective model in which the top quark is integrated out, provides a simple and reliable model to test this regularization scheme, as shown in [30]. Such a model has been used in [31] to test the four-dimensional regularization/renormalization (FDR) method, while in [32] it was studied in the context of the four-dimensional-unsubtraction method (FDU).

The main goal of this work is the computation of the total decay rate of this process in IReg to NLO to show that, contrarily to regularization schemes that work partially in the physical dimension, no modification of the Lagrangian such as through the inclusion of $\epsilon$-fields is required. This amounts to a significant simplification which we expect to prevail in calculations beyond NLO. We compute the one-loop virtual diagrams that arise from this effective model and apply IReg to extract the UV divergences which are absorbed in the process of renormalization, using the remaining UV finite amplitudes to obtain the regularized virtual decay rate. Then we apply the spinor-helicity formalism to compute the real diagrams of the process and obtain the real decay rate. It is expected that the IR divergences cancel after combining these decay rates, leading to a finite final result.

Our work is organized as follows: in Sect. 2 we briefly review the IReg method, which is applied in Sect. 3 to the decay $H \longrightarrow gg$ at NLO. Section 4 is devoted to a comparison of our results to dimensional schemes. Finally, we conclude in Sect. 5.

2 The IReg method: UV/IR identification and UV renormalization

IReg is a regularization method that operates on the momentum space and was shown to respect unitarity, locality and Lorentz invariance [6]. This procedure operates on the specific physical dimension of the theory, therefore we do not need to extend the space-time dimensions. IReg also does not require any changes to the Lagrangian and can be applicable to arbitrary n-loop calculations, making it an alternative to dimensional schemes. In this work we will be concerned with one-loop examples only, a recent review of the method applicable to n-loop amplitudes can be found in [33].

The main idea of IReg is to use an algebraic identity at integrand level recursively until the UV divergent behavior is only present in irreducible loop integrals that depend on internal momentum. In this way, the UV finite content of the amplitude (that may still be IR divergent) will contain denominators with dependence on physical parameters (external momenta and masses). To be concrete, we illustrate the method with the following massless toy integral in Minkowsky momentum space, with p denoting an arbitrary external 4-momentum

$$\begin{aligned} \int _k \dfrac{1}{k^2 (k-p)^2}, \quad \int _k \equiv \int \frac{d^{4}k}{(2\pi )^{4}}. \end{aligned}$$

(2.1)

We start³ by introducing an infrared regulator $\mu$ in the denominator like

$$\begin{aligned} \lim _{\mu \rightarrow 0}\int _k \dfrac{1}{(k^2-\mu ^2) [(k-p)^2 - \mu ^2]}. \end{aligned}$$

(2.2)

In the case of IR safe integrals, the regulator $\mu$ is needed to avoid spurious IR divergences in the course of the evaluation. It will cancel in the end result. In the case of IR divergent integrals, the $\mu$ will survive and parameterize the IR divergences.

By power counting, we notice that as $k \longrightarrow \infty$ this integral diverges, but there is a dependence both on internal and external momenta. We want to isolate the UV divergent content in an integral that is solely dependent on the internal momentum k. We notice that this is possible by rewriting the portion of the integrand that depends on external momentum as

$$\begin{aligned} \dfrac{1}{(k-p)^2 - \mu ^2} = \dfrac{1}{k^2-\mu ^2} + \dfrac{2k\cdot p - p^2}{(k^2-\mu ^2)((k-p)^2-\mu ^2)}, \end{aligned}$$

(2.3)

where on the right hand side the second term diminishes the divergence of the integral by one order and the first term leads to an integral that depends only on the internal momentum. The procedure exemplified above works in general, by repetitive usage of Eq. (2.3).

Following the separation of the divergences of the amplitude, the UV divergent content of the amplitude can be expressed by integrals with denominators that depend only on the internal momentum k. These integrals are classified as Basic Divergent Integrals (BDI’s) and they can take either logarithmic or quadratic forms which are respectively

$$\begin{aligned}{} & {} I_{log}^{\nu _1...\nu _{2r}}(\mu ^2) = \int _k \dfrac{k^{\nu _1}...k^{\nu _{2r}}}{(k^2 - \mu ^2)^{r+2}},\quad \text{ and }\quad \nonumber \\{} & {} I_{quad}^{\nu _1...\nu _{2r}}(\mu ^2) = \int _k \dfrac{k^{\nu _1}...k^{\nu _{2r}}}{(k^2 - \mu ^2)^{r+1}}. \end{aligned}$$

(2.4)

Here ${\nu _1}...\nu _{2r}$ represent Lorentz indices, there are as many momenta in the numerator as the number of integers comprised in the interval 1, ...2r, counted in unit steps for $r \ge 1$. The case $r=0$ corresponds to no momenta in the numerator. Any BDI with odd power of k in the numerator is automatically zero once the integral goes over the entire 4-momentum space and all the denominators have even powers of k. The BDI’s are written in terms of Lorenz indices and can be rewritten as scalar integrals, multiplying metric tensors, after setting surface terms (ST) to zero, as we will explain shortly. Scalar logarithmic and quadratic divergent integrals are finally given as

$$\begin{aligned} I_{log}(\mu ^2)= & {} \int _k \dfrac{1}{(k^2 - \mu ^2)^2},\quad \text{ and }\nonumber \\ I_{quad}(\mu ^2)= & {} \int _k \dfrac{1}{(k^2 - \mu ^2)}. \end{aligned}$$

(2.5)

The ultraviolet renormalization using IReg has been shown to comply with Lorentz and gauge symmetry by restraining local and momentum routing dependent ST to zero. IReg can be systematized order by order in the loop expansion in such a way that renormalization constants are given in terms of loop integrals defining a local renormalization scheme, [6, 12, 13]. In the method, symmetry breaking terms can be expressed as a well defined difference between divergent integrals with the same superficial degree of freedom. These are called ST and they are not originally fixed, which indicates that they are related to momentum routing invariance in Feynman diagrams (the possibility to perform a shift in the integration variables). As their value is associated with symmetry breaking terms, they play a critical role in IReg for the preservation of the symmetries of the system and we must carefully choose a value that allows the symmetries of the underlying theory to be preserved. Nonetheless, in a constrained version of IReg, it has been proven that these regularization dependent ST may be set to zero, complying with gauge invariance, [7, 34]. This will actually allow us to reduce the BDI’s with Lorenz indices $\nu _1...\nu _{2r}$ to linear combinations of scalar integrals with the same degree of divergence (multiplied by metric tensor combinations), plus well defined ST. Generally in the four dimensional Minkowskian space-time a ST of order i can be written as

$$\begin{aligned} \Gamma _i^{\nu _1 ... \nu _j} = \int _k \dfrac{\partial }{\partial k_{\nu _1}} \dfrac{k^{\nu _2}...k^{\nu _j}}{(k^2 - \mu ^2)^{(2+j-i)/2}}, \end{aligned}$$

(2.6)

with $\textit{j} \ge 2$ and $\mu$ the infrared regulator introduced earlier. The general formula allows for the computation of ST that appear in any order of perturbation theory. For instance, with the subscript $\textit{i}=0,1,2$ are designated the ST that involve differences of divergent integrals of logarithmic, linear or quadratic order respectively. As an example, one obtains

$$\begin{aligned} \Gamma _0^{\nu _1 \nu _2}= & {} \int _k \dfrac{\partial }{\partial k_{\nu _1}} \dfrac{k^{\nu _2}}{(k^2 - \mu ^2)^2} \nonumber \\= & {} 4 \Big ( \dfrac{g^{\nu _1 \nu _2}}{4} I_{log}(\mu ^2) - I_{log}^{\nu _1 \nu _2}(\mu ^2) \Big ) = 0 \end{aligned}$$

(2.7)

Finally, when considering IR-safe integrals, one must still take the limit in which the infrared regulator $\mu$ is set to zero. In this case, one rewrites the BDI’s in terms of a positive arbitrary constant $\lambda$ which will play the role of the renormalization group scale. It is achieved by using the equation below

$$\begin{aligned} I_{log}(\mu ^2) = I_{log}(\lambda ^2) + b \ln \Big ( \dfrac{\lambda ^2}{\mu ^2} \Big ), \quad b = \dfrac{i}{(4 \pi )^2}. \end{aligned}$$

(2.8)

It is worth noticing that a minimal subtraction renormalization scheme emerges naturally from this formalism, in which the infinite divergences that depend only on the internal momentum are subtracted from the theory. This means that the $I_{log}(\lambda ^2)$ will be subtracted via renormalization whereas the IR divergent part $\ln (\mu ^2)$ will cancel in the final amplitude for infrared safe processes and in the cross section/decay rate, which are IR-safe observables.

3 NLO corrections to $H\rightarrow gg$ in the large top mass limit

As discussed in the introduction, we will use the decay $H\rightarrow gg$ as a working example to test the method of IReg in the presence of both UV/IR divergences using an effective non-abelian field theory approach. Our objective is twofold: (a) the renormalization of an effective field theory is highly non-trivial, in particular, when considering alternatives to CDR [24, 30]. Thus, it is essential to understand how IReg can be applied in this context. (b) The presence of both UV/IR divergences requires a precise match between virtual and real contributions in order that a finite and regularization independent result occurs. Thus, it is a stringent test for any regularization scheme. Finally, the decay $H\rightarrow gg$ has served as a benchmark for other regularization methods, allowing for a clear comparison among methods [24, 30, 31].

As usual, since the decay we consider is mainly due to the top quark loop, it is reliable to consider the limit in which its mass is infinite. Thus, we add the following term to the massless QCD Lagrangian [35,36,37]

$$\begin{aligned} L_{eff} = \dfrac{1}{4}AH G_{\mu \nu }^{a} G^{a,\mu \nu }, \end{aligned}$$

(3.1)

where H represents the Higgs boson field, $G_{\mu \nu }^a$ is the field strength tensor of the SU(3) gluon field given by

$$\begin{aligned} G_{\mu \nu }^a = \partial _\mu A_\nu ^a - \partial _\nu A_\mu ^a + g_{s} f^{abc} A_\mu ^b A_\nu ^c \end{aligned}$$

(3.2)

and $f^{abc}$ are the anti-symmetric SU(3) structure constants. The effective coupling A can be obtained by performing the matching of the full theory to its effective version, so that [38,39,40],

$$\begin{aligned} A = \dfrac{\alpha _s}{3 \pi v}\Big ( 1 + \dfrac{11}{4} \dfrac{\alpha _s}{\pi }\Big ), \end{aligned}$$

(3.3)

with ${\alpha _s}=\frac{g_s^2}{4\pi }$ denoting the strong coupling and v the electroweak vacuum expectation value, $v^2=(G_f\sqrt{2})^{-1}$ with $G_f$ the Fermi constant. The Feynman rules can be straightforwardly obtained, [37]. For the diagrams involving only gluons the Feynman rules are given by the Yang Mills Lagrangian.

Once the model is defined, we present in the next subsections its UV renormalization at one-loop level, as well as the calculation of the virtual and real contributions for the decay $H\rightarrow gg$.

3.1 UV renormalization

As usual we adopt multiplicative renormalization, rewriting the effective Lagrangian as

$$\begin{aligned} (L_{eff})_{ren} = \dfrac{1}{4} Z_{\alpha _s} Z_A AH G_{\mu \nu } G^{\mu \nu }, \end{aligned}$$

(3.4)

where $Z_{A}$ and $Z_{\alpha _s}$ are the renormalization constants for the gluon-field and coupling constant respectively. Notice that we do not renormalize the Higgs field, since it can only appear as an external leg in the process we consider. The part of the Lagrangian corresponding to massless QCD is renormalized in the standard way, implying that $Z_{A}$ and $Z_{\alpha _s}$ are already known. In the framework of IReg, they are given in [41]. At first order in $\alpha _{s}$, $Z_{A}$ and $Z_{\alpha _s}$ are given by

$$\begin{aligned} Z_A&= 1\! + \!\alpha _s \dfrac{1}{(4\pi )} \dfrac{1}{b} I_{log}(\mu ^2)\Big [ \Big ( \dfrac{13}{6}\! -\! \dfrac{\zeta }{2} \Big )C_A - \dfrac{4}{3} T_F N_F \Big ] + O(\alpha _s^2) \end{aligned}$$

(3.5)

$$\begin{aligned} Z_{\alpha _{s}}&= 1 - \alpha _s \dfrac{1}{(4\pi )} \dfrac{1}{b} I_{log}(\lambda ^2) \Big [ \dfrac{11}{6} C_A - \dfrac{2}{3} T_F N_F \Big ] +O(\alpha _s^2) \end{aligned}$$

(3.6)

where $\zeta$ is the gauge parameter, $N_F$ is the number of light quarks flavours, $T_F=\frac{1}{2}$ results from the trace over colour matrices and $C_A=3$ is a color factor. Notice that, while for $Z_{\alpha _{s}}$ we have the UV divergence expressed as $I_{log}(\lambda ^2)$ as usual (minimal subtraction scheme), for $Z_{A}$ we have $I_{log}(\mu ^2)$. This happens because we are considering on-shell gluons, $\mu ^{2}$ playing the role of their fictitious mass.

Finally, by considering the terms just discussed, the counterterm to be added to our process is (we adopt Feynman Gauge, $\zeta =1$)

$$\begin{aligned} V_{count}= & {} \dfrac{\alpha _s}{b \pi } \left[ C_A \Big (\dfrac{5}{12} I_{log}(\mu ^2) - \dfrac{11}{12} I_{log}(\lambda ^2) \Big )\right. \nonumber \\{} & {} \left. - \dfrac{1}{3} T_F N_F \Big ( I_{log}(\lambda ^2) - I_{log}(\mu ^2) \Big ) \right] V_0\,, \end{aligned}$$

(3.7)

where $V_{0}$ corresponds to the tree-level amplitude for $H\rightarrow gg$, given below

$$\begin{aligned} V_0= i A \delta ^{ab} H^{\mu \nu }\,\quad H^{\mu \nu }(p1,p2)=-p_1^\nu p_2^\mu + g^{\mu \nu } p1\cdot p2\, . \end{aligned}$$

(3.8)

Notice that $V_{0}$ has an implicit dependence on color and Lorentz indexes, which we suppress for simplicity. The same argument holds for all amplitudes $V_{i}$ to be defined in the next section.

3.2 Virtual contributions

We can now compute the one-loop virtual diagrams. There are 5 diagrams that contribute to the one-loop order correction, which are represented in Fig. 1⁴

Fig. 1

Virtual diagrams contribution to the decay rate $H \longrightarrow gg(g)$. From left to right they are respectively $V_1$,$V_2$,$V_3$,$V_4$,$V_5$. The dashed line represents the Higgs field, the curly lines represent the gluon field

For all the diagrams we choose the external momenta of the two gluons to be $p_1$ and $p_2$, the momentum of the Higgs boson to be q, and the internal momentum of the loop to be k. All the external momenta are inwards, therefore we can write the equation of momentum–energy conservation as $p_1+p_2+q=0$. We apply the on-shell conditions by imposing $p_1^2=p_2^2=0$. The results for the integrals evaluated in this section can be found in Appendix A.

We begin with the diagram $V_1$ whose amplitude is given by

$$\begin{aligned} \begin{aligned} V_1&=- A g^2 C_A \delta ^{ab} \\&\quad \times \Bigg ( I_4 g^{\mu \nu } + 2I^{\alpha _1 \alpha _2} g^{\mu \nu } (p_1^{\alpha _1}p_1^{\alpha _2} + p_2^{\alpha _1}p_2^{\alpha _2}) + 11 I_2^{\mu \nu } \\&\quad - I_2^\nu ( 6p_1^\mu + p_2^\mu ) + I_2^\mu (p_1 + 6 p_2)^\nu \\&\quad - I (p_1 \cdot p_2) (4 p_2^\mu p_1^\nu + p_1^\mu p_2^\nu - 4(p_1 \cdot p_2)g^{\mu \nu })\\&\quad + I_2 (-9 (p_1 \cdot p_2) g^{\mu \nu } + 9 p_2^\mu p_1^\nu - 3 p_1^\mu p_2^\nu + p_2^\mu p_2^\nu ) \\&\quad + 10 I^{\alpha _1 \mu \nu } (p_1 + p_2)^{\alpha _1} + I^{\alpha _1 \mu } (-2 p_1^{\alpha _1} p_1^\nu \\&\quad - 6 p_1^{\alpha _1} p_2^\nu -6 p_2^{\alpha _1} p_1^\nu + p_2^{\alpha _1} p_2^\nu )\\&\quad + I^{\alpha _1 \nu } (-6 p_2^{\alpha _1} p_1^\mu - 2 p_2^{\alpha _1} p_2^\mu \\&\quad + p_1^{\alpha _1} p_1^\mu -6 p_1^{\alpha _1} p_2^\mu ) \\&\quad + I^{\alpha _1} (-6 p_1^{\alpha _1} p_2^\mu p_1^\nu + p_1^{\alpha _1} p_1^\mu p_2^\nu - 3 p_1^{\alpha _1} p_2^\mu p_2^\nu \\&\quad + 3 p_2^{\alpha _1} p_1^\mu p_1^\nu + 6 p_2^{\alpha _1} p_2^\mu p_1^\nu - p_2^{\alpha _1} p_1^\mu p_2^\nu \\&\quad + 6 p_1 \cdot p_2 g^{\mu \nu } (p_1^{\alpha _1}-p_2^{\alpha _1})) + 3 I_2^{\alpha _1} g^{\mu \nu } (-p_1^{\alpha _1} \\&\quad + p_2^{\alpha _1}) -(p_1 \cdot p_2) p_1^\mu I^{\nu } + I^{\mu } (p_1 \cdot p_2) p_2^\nu \Bigg ) \end{aligned} \end{aligned}$$

(3.9)

We have already applied the on-shell conditions, $p_i^2=0$, $i=1,2$, and we adopted the following convention for the integrals

$$\begin{aligned} I^{\alpha _{1}\cdots \alpha _{n}}&= \int _{k}\frac{k^{\alpha _{1}}\cdots k^{\alpha _{n}}}{k^2(k-p_1)^2(k+p_2)^2} \end{aligned}$$

(3.10)

$$\begin{aligned} I_{2}^{\alpha _{1}\cdots \alpha _{n}}&= \int _{k}\frac{k^{2} k^{\alpha _{1}}\cdots k^{\alpha _{n}}}{k^2(k-p_1)^2(k+p_2)^2} \end{aligned}$$

(3.11)

$$\begin{aligned} I_{4}^{\alpha _{1}\cdots \alpha _{n}}&= \int _{k}\frac{k^{4} k^{\alpha _{1}}\cdots k^{\alpha _{n}}}{k^2(k-p_1)^2(k+p_2)^2} \end{aligned}$$

(3.12)

After using the results collected in Appendix A, the final result reads

$$\begin{aligned} V_1&= A g^2 C_A \delta ^{ab} \Big [(-I_{quad}(\mu ^2) \bigg (\dfrac{13}{2} g^{\mu \nu }\bigg ) \nonumber \\&\quad -I_{log}(\mu ^2) \bigg (- \dfrac{43}{6} p_1 \cdot p_2 g^{\mu \nu } + \dfrac{1}{4} p_1^\mu p_1^\nu - \dfrac{1}{6} p_1^\mu p_2^\nu \nonumber \\&\quad + \dfrac{29}{6}p_1^\nu p_2^\mu + \dfrac{1}{4}p_2^\mu p_2^\nu \bigg )\nonumber \\&\quad + \dfrac{1}{12} \ln (-\mu _0)\bigg (-2 p_1^\nu p_2^\mu - \dfrac{13}{2} p_1 \cdot p_2 g^{\mu \nu }\bigg )\nonumber \\&\quad + \ln (-\mu _0)^2 (- p_1^\nu p_2^\mu + p_1 \cdot p_2 g^{\mu \nu })\nonumber \\&\quad - \dfrac{5}{18} ( p_1^\nu p_2^\mu + 4 p_1 \cdot p_2 g^{\mu \nu })\Big ] \end{aligned}$$

(3.13)

where $\mu _{0}=\mu ^{2}/m_{H}^{2}$.

For the diagram $V_2$ we obtain

$$\begin{aligned} \begin{aligned} V_2&= A g^2 C_A \delta ^{ab} \dfrac{1}{2} \int _k \dfrac{1}{k^2(k-p_1-p_2)^2} \\&\quad \times \Big ( 4 k^2 g^{\mu \nu } - 4k \cdot (p_1+p_2) g^{\mu \nu }\\&\quad + 2k^\mu k^\nu -k^\mu (p_1^\nu + p_2^\nu ) -k^\nu (p_1^\mu +p_2^\mu ) \Big ) \end{aligned} \end{aligned}$$

(3.14)

and the regularized amplitude is

$$\begin{aligned} V_2&= A g^2 C_A \delta ^{ab} \Big [ I_{quad}(\mu ^2) \bigg (\dfrac{5}{2} g^{\mu \nu }\bigg ) \nonumber \\&\quad + \dfrac{1}{2} I_{log}(\mu ^2) \bigg (- \dfrac{13}{3} p_1 \cdot p_2 g^{\mu \nu } - \dfrac{1}{3} p_1^\mu p_1^\nu \nonumber \\&\quad - \dfrac{1}{3}p_1^\mu p_2^\nu - \dfrac{1}{3}p_1^\nu p_2^\mu - \dfrac{1}{3}p_2^\mu p_2^\nu \bigg )\nonumber \\&\quad + \dfrac{1}{12} \ln (-\mu _0)\bigg (2 p_1^\nu p_2^\mu - \dfrac{13}{2} p_1 \cdot p_2 g^{\mu \nu } \bigg )\nonumber \\&\quad + \dfrac{5}{18} ( p_1^\nu p_2^\mu + 4 p_1 \cdot p_2 g^{\mu \nu })\Big ]. \end{aligned}$$

(3.15)

For the diagram $V_3$ we obtain

$$\begin{aligned} V_3= & {} \dfrac{1}{2} Ag^2C_A \delta ^{ab} \int _k \dfrac{1}{k^2(k-p_2)^2}\Big (2k^2 g^{\mu \nu } -2 k \cdot p_2 g^{\mu \nu } \nonumber \\{} & {} -3 p_1 \cdot p_2 g^{\mu \nu } + 2 p_2^2 g^{\mu \nu } + 10 k^\mu k^\nu -5 k^\nu p_{2}^\mu \nonumber \\{} & {} -5 k^\mu p_{2}^\nu + 3 p_{1}^\nu p_{2}^\mu + p_{2}^\mu p_{2}^\nu \Big ) \end{aligned}$$

(3.16)

and are evaluated in IReg as

$$\begin{aligned} \begin{aligned} V_3&= A g^2 C_A \delta ^{ab}\Bigg [ (I_{quad}(\mu ^2) \bigg (\dfrac{7}{2} g^{\mu \nu }\bigg ) \\&\quad + \dfrac{1}{2} I_{log}(\mu ^2) \bigg (3 p_1^\nu p_2^\mu - 3p_1 \cdot p_2 g^{\mu \nu } - \dfrac{2}{3} p_2^\mu p_2^\nu \bigg )\Bigg ]. \end{aligned} \end{aligned}$$

(3.17)

The diagram $V_4$ can be obtained from the result of $V_3$ by substituting $p_{1}\leftrightarrow p_{2}$, $\mu \leftrightarrow \nu$

$$\begin{aligned} \begin{aligned} V_4&= A g^2 C_A \delta ^{ab}\Bigg [( I_{quad}(\mu ^2) \bigg (\dfrac{7}{2} g^{\mu \nu }\bigg ) \\&\quad + \dfrac{1}{2} I_{log}(\mu ^2) \bigg (3 p_1^\nu p_2^\mu - 3p_1 \cdot p_2 g^{\mu \nu } - \dfrac{2}{3} p_1^\mu p_1^\nu \bigg )\Bigg ]. \end{aligned} \end{aligned}$$

(3.18)

Finally, the regularized amplitude of diagram $V_5$ is given by

$$\begin{aligned} V_{5}=A g^2 C_A \delta ^{ab} \int _k \dfrac{-3 g^{\mu \nu }}{k^2-\mu ^2} = -3 A g^2 C_A \delta ^{ab} I_{quad}(\mu ^2) g^{\mu \nu }. \end{aligned}$$

(3.19)

Once all amplitudes are regularized, we obtain an UV divergent part given by

$$\begin{aligned} V_{div} = \dfrac{\alpha _s}{\pi } C_A\left[ \dfrac{I_{log}(\mu ^2)}{2b}\right] V_0\;, \end{aligned}$$

(3.20)

in terms of the tree level amplitude, Eq. (3.8). We emphasize that only terms up to $\mathcal {O}(\alpha _s^{2})$ are retained, implying that only the first term of Eq. (3.3) is to be considered in the definition of $V_{0}$ above. It is worth noticing that quadratic divergent integrals cancel in the sum, being an illustration of consistency of the method. Finally, the UV finite part yields (as can be easily checked, only the diagrams $V_1$ and $V_2$ contribute)

$$\begin{aligned} V_{rest} = \dfrac{\alpha _s}{\pi }C_{A} \left[ -\dfrac{\ln (\mu _0)^2}{4} - \dfrac{i \pi \ln (\mu _0)}{2} + \dfrac{\pi ^2}{4}\right] V_0\,, \end{aligned}$$

(3.21)

where we used $\ln (-\mu _0)^2 = \ln (\mu _0)^2 + 2 i \pi \ln (\mu _0) - \pi ^2$.

At this point we can check if our result is UV finite after adding the countertem obtained in the last subsection (Eq. 3.7)

$$\begin{aligned} V_{ren}= & {} V_{div}+V_{count}\nonumber \\= & {} \dfrac{\alpha _s}{b \pi } \left[ \Big ( I_{log}(\lambda ^2) - I_{log}(\mu ^2) \Big ) \Big ( \dfrac{11}{12}C_A - \dfrac{1}{3} T_f N_F \Big ) \right] V_0. \end{aligned}$$

(3.22)

Using the scaling relation, Eq. 2.8, we obtain

$$\begin{aligned} V_{ren} = \dfrac{\alpha _s}{\pi } \left[ \Big ( \dfrac{11}{12}C_A - \dfrac{1}{3} T_f N_F \Big ) \ln \Big ( \dfrac{\lambda ^2}{\mu ^2} \Big ) \right] V_0 , \end{aligned}$$

(3.23)

rendering an UV finite result as expected.

Finally, the virtual decay rate can be obtained by considering the sum of the tree-level amplitude with the one-loop radiative correction

$$\begin{aligned} V = V_0\left( 1 + \dfrac{11}{4} \dfrac{\alpha _s}{\pi }\right) + V_{ren} + V_{rest}. \end{aligned}$$

(3.24)

In the above equation, the tree level amplitude encompasses also the correction due to the matching of the effective theory to the full SM as given by Eq.  3.3. By squaring this result one obtains, up to the order $\alpha _{s}^{3}$

$$\begin{aligned} |V|^2= & {} |V_0|^2 \left[ 1 + \dfrac{\alpha _s}{\pi }\left( \dfrac{11}{2} + \Big (\dfrac{11}{6}C_A - \dfrac{2}{3} T_f N_F\Big )\right. \right. \nonumber \\{} & {} \left. \left. \times \ln \Big (\dfrac{\lambda ^2}{\mu ^2} \Big ) +\dfrac{C_A}{2}\Big (-\ln (\mu _0)^2+ \pi ^2\Big ) \right) \right] \end{aligned}$$

(3.25)

and the virtual decay rate is given by

$$\begin{aligned} \Gamma _v= & {} \Gamma _0 \left[ 1 + \dfrac{\alpha _s}{\pi }\left( \dfrac{11}{2} -\Big (\dfrac{11}{6}C_A - \dfrac{1}{3} N_F\Big )\right. \right. \nonumber \\{} & {} \left. \left. \times \ln (\mu _0) +\dfrac{C_A}{2}\Big (-\ln (\mu _0)^2+ \pi ^2\Big ) \right) \right] \,, \end{aligned}$$

(3.26)

where we are choosing the renormalization scale at the Higgs mass ($\lambda ^{2}=m_{H}^{2}$) since $\mu _{0}=\mu ^{2}/m_{H}^{2}$.

3.3 Real decay rate

The diagrams that will contribute to the real decay are represented in Fig. 2 up to $\alpha _s \sqrt{\alpha _s}$ order

Fig. 2

Real diagrams contributing to the decay $H \longrightarrow ggg$ and $H \longrightarrow gq\bar{q}$. From left to right they are respectively $R_1$, $R_2$ and $R_3$. The dashed line represents the Higgs field and the curly lines represent the gluon field. The $\{p_i,p_j,p_k\}$ correspond to the three permutations of $p_i$, $p_j$ and $p_k$, so $R_1$ stands for 3 diagrams

We consider first the diagrams with only gluons as external legs, which are the more involved to be obtained. In order to simplify the calculations we adopt the spinor helicity formalism in this case. For the diagram with external (light) quarks, we use the standard procedure. We begin with the s-channel diagram $R_1$, whose amplitude is given by

$$\begin{aligned} iM_{R_1}= & {} \, g_{s}f^{bcd}V^{\nu \tau \delta }(p_2,p_3,-(p_2+p_3)) i \dfrac{g^{\delta \delta ^{'}} \delta ^{d d^{'}}}{(p_2+p_3)^2} \nonumber \\{} & {} iA \delta ^{ad^{'}} H^{\delta ^{'} \mu }(-(p_2+p_3),p_1) \epsilon _1^{\mu } \epsilon _2^{\nu } \epsilon _3^{\tau }\nonumber \\= & {} -Ag_{s}f^{bca} \dfrac{1}{(p_2+p_3)^2} V^{\nu \tau \delta } H^{\delta \mu } \epsilon _1^{\mu } \epsilon _2^{\nu } \epsilon _3^{\tau }\,, \end{aligned}$$

(3.27)

where the expanded tensors are given by

$$\begin{aligned} V^{\nu \tau \delta }(p_2,p_3,-(p_2+p_3)) = (p_2-p_3)^\delta g^{\nu \tau } - p_2^\nu g^{\tau \delta } + p_3^\tau g^{\delta \nu } \end{aligned}$$

(3.28)

and

$$\begin{aligned} H^{\delta \mu }(-(p_2+p_3),p_1) = - g^{\mu \delta }p_1 \cdot (p_2+p_3) + p_1^\delta (p_2+p_3)^\mu . \end{aligned}$$

(3.29)

Contracting the indices and using the Lorenz condition $\epsilon ^\mu p_\mu =0$, one obtains

$$\begin{aligned} \begin{aligned} iM_{R_1}&= \dfrac{-Ag_{s}f^{bca}}{s_{23}} (s_{12}+s_{13})(\epsilon _1 \cdot \epsilon _3 \, p_3 \cdot \epsilon _2 - p_2 \cdot \epsilon _3 \, \epsilon _1 \cdot \epsilon _2 )\\&\quad - s_{12} \, p_3 \cdot \epsilon _1 \, \epsilon _2 \cdot \epsilon _3 + s_{13} \, p_2 \cdot \epsilon _1 \, \epsilon _2 \cdot \epsilon _3\\&\quad + 2 (p_2 \cdot \epsilon _1 + p_3 \cdot \epsilon _1) (p_1 \cdot \epsilon _2 \, p_2 \cdot \epsilon _3 - p_1 \cdot \epsilon _3 \, p_3 \cdot \epsilon _2). \end{aligned} \end{aligned}$$

(3.30)

Now we proceed using the spinor helicity formalism. Firstly, we define three auxiliary momenta $r_i$, $i=1,2,3$, one for each of the massless gluons

$$\begin{aligned} r(\epsilon _1)&=p_3 \equiv 3,&r(\epsilon _2)&=p_1 \equiv 1,&r(\epsilon _3)&=p_2 \equiv 2. \end{aligned}$$

(3.31)

With this choice, the terms $p_2 \cdot \epsilon _3=p_1 \cdot \epsilon _2=p_3 \cdot \epsilon _1=0$ are automatically zero, which allows for a great deal of simplification of the amplitude,

$$\begin{aligned} iM_{R_1}= & {} \dfrac{-Ag_{s}f^{bca}}{s_{23}}\Big [ (s_{12}+s_{13})(\epsilon _1 \cdot \epsilon _3 \, p_3 \cdot \epsilon _2)\nonumber \\{} & {} + s_{13} \, p_2 \cdot \epsilon _1 \, \epsilon _2 \cdot \epsilon _3 - 2 p_2 \cdot \epsilon _1 \, p_1 \cdot \epsilon _3 \, p_3 \cdot \epsilon _2 \Big ].\nonumber \\ \end{aligned}$$

(3.32)

The t channel is obtained by making the replacements $1 \longleftrightarrow 3$, $2 \longleftrightarrow 1$ and $3 \longleftrightarrow 2$ while the u channel is obtained by making $1 \longleftrightarrow 2$, $2 \longleftrightarrow 3$ and $3 \longleftrightarrow 1$.

The amplitude for $R_2$ can be obtained in a similar fashion. The end result is

$$\begin{aligned} iM_{R_2} = -Ag_{s}f^{abc} \Big (p_1 \cdot \epsilon _3\, \epsilon _1 \cdot \epsilon _2 + p_2 \cdot \epsilon _1\, \epsilon _2 \cdot \epsilon _3 + p_3 \cdot \epsilon _2\, \epsilon _1 \cdot \epsilon _3 \Big ). \end{aligned}$$

(3.33)

Summing the s, t and u channels from $R_1$ with the one vertex diagram from $R_2$ we have

$$\begin{aligned} \begin{aligned} iM_g&= - Ag_s \Bigg ( \dfrac{f^{bca}}{s_{23}} \big [(s_{12}+s_{13})(\epsilon _1 \cdot \epsilon _3 \, p_3 \cdot \epsilon _2 )\\&\quad + s_{13} \, p_2 \cdot \epsilon _1 \, \epsilon _2 \cdot \epsilon _3 - 2 p_2 \cdot \epsilon _1 \, p_1 \cdot \epsilon _3 \, p_3 \cdot \epsilon _2 \big ] \\&\quad + \dfrac{f^{abc}}{s_{12}} \big [(s_{13}+s_{23}) (\epsilon _3 \cdot \epsilon _2 \, p_2 \cdot \epsilon _1 ) + s_{23} \, p_1 \cdot \epsilon _3 \, \epsilon _1 \cdot \epsilon _2 \\&\quad - 2 p_1 \cdot \epsilon _3 \, p_3 \cdot \epsilon _2 \, p_2 \cdot \epsilon _1 \big ]\\&\quad + \dfrac{f^{cab}}{s_{13}} \big [ (s_{23}+s_{12})(\epsilon _1 \cdot \epsilon _2 \, p_1 \cdot \epsilon _3) + s_{12} \, p_3 \cdot \epsilon _2 \, \epsilon _3 \cdot \epsilon _1 \\&\quad - 2 p_3 \cdot \epsilon _2 \, p_2 \cdot \epsilon _1 \, p_1 \cdot \epsilon _3 \big ]\\&\quad + f^{abc} ( p_1 \cdot \epsilon _3\, \epsilon _1 \cdot \epsilon _2 + p_2 \cdot \epsilon _1\, \epsilon _2 \cdot \epsilon _3 + p_3 \cdot \epsilon _2\, \epsilon _1 \cdot \epsilon _3 ) \Bigg ). \end{aligned} \end{aligned}$$

(3.34)

To obtain the unpolarized absolute squared value of the amplitude we need to sum over all possible colors and helicities,

$$\begin{aligned} |\overline{M_g}|^2= & {} \sum _{col, polr}|M_g|^2 =\sum _{col} 2(|M_g^{+++}|^2\nonumber \\{} & {} +|M_g^{+--}|^2+|M_g^{-+-}|^2+|M_g^{--+}|^2) \end{aligned}$$

(3.35)

and use the spinor helicity formalism to perform the spin sums (see e.g. [37]). A massless real valued four momentum vector $p^\mu$ in the representation

$$\begin{aligned} p^{\alpha \dot{\alpha }}= \sigma _\mu ^{\alpha \dot{\alpha }}p^\mu , \quad p_{ \dot{\alpha } \alpha } = \bar{\sigma }_{\dot{\alpha } \alpha }^\mu p_\mu \end{aligned}$$

(3.36)

where $\sigma _\mu = \big ( \mathbb {I}, \vec {\sigma } \big )$ and $\bar{\sigma }^\mu = \big ( \mathbb {I},- \vec {\sigma } \big )$ allows a bispinor decomposition,

$$\begin{aligned} p^{\alpha \dot{\alpha }} = p \rangle [p, \quad p_{ \dot{\alpha } \alpha } = p] \langle p \end{aligned}$$

(3.37)

with the properties

$$\begin{aligned} \langle p q \rangle = \sqrt{2 p \cdot q} e^{i \phi }, \quad [ q p ] = \sqrt{2 p \cdot q} e^{-i \phi } \end{aligned}$$

(3.38)

for an arbitrary $\phi$ phase and

$$\begin{aligned} p \cdot q = q^\mu p_\mu = \dfrac{1}{2} q_{\dot{\alpha } \alpha } p^{\alpha \dot{\alpha }} = \dfrac{1}{2} \langle q p \rangle [p q]. \end{aligned}$$

(3.39)

The relation between these objects and the usual Mandelstam variables which we use to compute the amplitudes is

$$\begin{aligned} \langle i j \rangle [j i] = 2p_i p_j = (p_i + p_j)^2 = s_{ij}\,. \end{aligned}$$

(3.40)

The polarization vectors are given in this notation by

$$\begin{aligned}{} & {} [\epsilon _p^-(r)]^{\alpha \dot{\alpha }} = \sqrt{2} \dfrac{p \rangle [r}{[pr]}, \end{aligned}$$

(3.41)

$$\begin{aligned}{} & {} [\epsilon _p^+(r)]^{\alpha \dot{\alpha }} = \sqrt{2} \dfrac{r \rangle [p}{\langle r p \rangle }, \end{aligned}$$

(3.42)

from which follow the inner products

$$\begin{aligned}{} & {} \epsilon ^-_{p}(r_1) \cdot \epsilon ^+_q(r_2)= \dfrac{\langle p r_2 \rangle [q r_1]}{[p r_1] \langle r_2 q \rangle },\nonumber \\{} & {} \epsilon ^-_{p}(r_1) \cdot \epsilon ^-_q(r_2) = \dfrac{\langle p q \rangle [r_2 r_1]}{[p r_1] [ q r_2]},\nonumber \\{} & {} \epsilon ^+_{p}(r_1) \cdot \epsilon ^+_q(r_2) = \dfrac{[ r_1 r_2 ] \langle q p \rangle }{\langle r_1 p \rangle \langle r_2 q \rangle },\end{aligned}$$

(3.43)

$$\begin{aligned}{} & {} \epsilon _p^-(r_1) \cdot q = \dfrac{1}{\sqrt{2}} \dfrac{\langle p q \rangle [q r_1]}{[p r_1]}, \quad \nonumber \\{} & {} \epsilon _p^+(r_1) \cdot q = \dfrac{1}{\sqrt{2}} \dfrac{\langle p q \rangle [q r_1]}{\langle r_1 p \rangle } \end{aligned}$$

(3.44)

with p, q being the physical momenta and $r_1$ and $r_2$ the reference momenta.

Using our previous choice of reference momenta, we have, for the ${+ - -}$ helicity configuration

$$\begin{aligned} M_g^{+ - -}= & {} -Ag_{s}f^{abc} \Bigg [ \Big (\dfrac{\langle 23\rangle }{[32]} \dfrac{\langle 23\rangle [12]}{\langle 31\rangle }\Big ) \Big ( \dfrac{s_{13}}{s_{23}} + \dfrac{s_{13}}{s_{12}} + \dfrac{s_{23}}{s_{12}} + 1\Big )\nonumber \\{} & {} -2\Big (\dfrac{1}{\sqrt{2}} \dfrac{\langle 23\rangle [12]}{\langle 31\rangle }\Big )\Big (\dfrac{1}{\sqrt{2}} \dfrac{\langle 13\rangle [21]}{[32]}\Big )\Big (\dfrac{1}{\sqrt{2}} \dfrac{\langle 32\rangle [13]}{[21]}\Big )\nonumber \\{} & {} \Big (\dfrac{1}{s_{23}} + \dfrac{1}{s_{12}} + \dfrac{1}{s_{13}}\Big ) \Bigg ]\,, \end{aligned}$$

(3.45)

and, for the ${+ + +}$ helicity case, one obtains

$$\begin{aligned} M_g^{+ + +}= & {} -Ag_{s}f^{abc} \nonumber \\{} & {} \times \Bigg [ \dfrac{[13]}{\langle 31\rangle } \dfrac{1}{\sqrt{2}} \dfrac{\langle 31\rangle [23]}{\langle 12\rangle } \Big (\dfrac{s_{12}+s_{13}}{s_{23}} + \dfrac{s_{12}}{s_{13}} + 1 \Big )\nonumber \\{} & {} + \dfrac{[23]}{\langle 32\rangle } \dfrac{1}{\sqrt{2}} \dfrac{\langle 23\rangle [12]}{\langle 31\rangle } \Big (\dfrac{s_{13}+s_{23}}{s_{12}} + \dfrac{s_{13}}{s_{23}} + 1 \Big ) \nonumber \\{} & {} + \dfrac{[12]}{\langle 21\rangle } \dfrac{1}{\sqrt{2}} \dfrac{\langle 12\rangle [31]}{\langle 23\rangle } \Big (\dfrac{s_{23}+s_{12}}{s_{13}} + \dfrac{s_{23}}{s_{12}} + 1 \Big )\nonumber \\{} & {} + \dfrac{1}{\sqrt{2}} \dfrac{\langle 23\rangle [12]}{\langle 31\rangle } \dfrac{1}{\sqrt{2}} \dfrac{\langle 12\rangle [31]}{\langle 23\rangle } \dfrac{1}{\sqrt{2}} \dfrac{\langle 31\rangle [23]}{\langle 12\rangle }\nonumber \\{} & {} \times \Big (\dfrac{1}{s_{12}} + \dfrac{1}{s_{13}} + \dfrac{1}{s_{23}}\Big ) \Bigg ]. \end{aligned}$$

(3.46)

Squaring the amplitudes, one gets

$$\begin{aligned} |M_g^{+ - -}|^2= & {} \dfrac{1}{2} A^2g_{s}^2 f^2 \dfrac{s_{23}^3}{s_{12}s_{13}}, \nonumber \\ |M_g^{+ + +}|^2= & {} \dfrac{1}{2} A^2g_{s}^2 f^2 \dfrac{m_H^8}{s_{12} s_{13} s_{23}}. \end{aligned}$$

(3.47)

where we used the momentum–energy conservation condition $m_H^2 = s_{12} + s_{13} + s_{23}$.

The remaining helicity configurations can be obtained by permutation of the momenta. The structure constants can be written as $f^2 = f_{abc} f^{abc} = 2 C_A^2 C_F$ with $C_F = \dfrac{N^2 - 1}{2N}$ and $C_A = N$ with $N=3$ for the SU(3) group. Finally, we obtain the unpolarized amplitude considering gluons only, $M_g$,

$$\begin{aligned} |\overline{M_g}|^2 = A^2 \pi \alpha _s 8 C_A^2 C_F \dfrac{1}{s_{12} s_{13} s_{23}} (s_{12}^4 + s_{13}^4 + s_{23}^4 + m_H^8), \end{aligned}$$

(3.48)

that can be written as

$$\begin{aligned} |\overline{M_g}|^2= & {} A^2 \pi \alpha _s 16 C_A^2 C_F \Bigg [ \dfrac{s_{12}^3}{s_{13} s_{23}} + \dfrac{s_{13}^3}{s_{12} s_{23}} + \dfrac{s_{23}^3}{s_{12} s_{13}}\nonumber \\{} & {} + 6(s_{12}+s_{13}+s_{23}) + \dfrac{2(s_{12}^2+s_{13}^2)+3s_{12}s_{13}}{s_{23}}\nonumber \\{} & {} +\dfrac{2(s_{13}^2+s_{23}^2)+3s_{13}s_{23}}{s_{12}}+ \dfrac{2(s_{12}^2+s_{23}^2)+3s_{12}s_{23}}{s_{13}} \Bigg ].\nonumber \\ \end{aligned}$$

(3.49)

For the case of light quarks in the external legs, diagram $R_3$, the calculation is easier. Its amplitude $M_q$ is given by

$$\begin{aligned} M_q= & {} i A \alpha _s t_b \epsilon _\nu (p_3) H^{\rho \nu }((-p_1-p_2),p_3) \bar{u}^s(p_1) \nonumber \\{} & {} \times \gamma _\rho v^{s'}(p_2) \dfrac{1}{(-p_1-p_2)^2 + i \epsilon }. \end{aligned}$$

(3.50)

By considering the unpolarized amplitude, one must sum over the color and spin degrees of freedom of the external particles. In this case, it is essential to recall that, in IReg, the IR divergences are parametrized as fictional masses for otherwise massless particles. Thus, when using the completeness relation for the sum over the spins of the massive light quarks in the limit of vanishing quark masses, one must retain the massive contribution until the integration over the phase space is effected, before taking the limit. By doing so, one obtains

$$\begin{aligned} |\overline{M_q}|^2 = A^2 \alpha _s 4 \pi \text {Tr}(t_b t^b) \Big ( \dfrac{(s_{13}^2 + s_{12}^2)}{s_{23}} + \dfrac{4 \mu ^2}{s_{23}^2} \dfrac{(s_{13} + s_{12})^2}{2} \Big ), \end{aligned}$$

(3.51)

where $\text {Tr}(t^b t_b) = \text {Tr} \Big ( C_F \mathbb {I} \Big ) = C_F C_A$.

The total unpolarized amplitude of the decay is finally given by

$$\begin{aligned} |\overline{M}|^{2} = |\overline{M_g}|^{2} + |\overline{M_q}|^{2} \end{aligned}$$

(3.52)

Our next task is to perform the phase space integral considering the external particles (gluons and light quarks) with mass $\mu$. As customary, one defines dimensionless variables which, in our case, are given by

$$\begin{aligned} s_{13}= & {} q^2(\chi _1 + \mu _0),\quad s_{23} = q^2(\chi _2 + \mu _0),\nonumber \\ s_{12}= & {} q^2(1-\chi _1-\chi _2 + \mu _0), \end{aligned}$$

(3.53)

where $\mu _0 = \dfrac{\mu ^2}{q^2}$ and $q^{2}=m_{H}^{2}$ (we are using the frame of reference in which the Higgs boson is at rest). In terms of these variables the decay rate is given by

$$\begin{aligned}&\Gamma _r (H \longrightarrow gg(g), gq \bar{q}) \nonumber \\&\quad = \Gamma _0 \dfrac{\alpha _s}{\pi } \int \Big [ C_A\Big (2 + 3 \chi _2 - \dfrac{4}{(\chi _2 + \mu _0)} + \dfrac{5\chi _1}{(\chi _2 + \mu _0)} \nonumber \\&\qquad - \dfrac{\chi _1^2}{(\chi _2 + \mu _0)} + \dfrac{1}{(\chi _1 + \mu _0) (\chi _2 + \mu _0)} \Big )\nonumber \\&\qquad + N_F \Big ( \dfrac{2\mu _0}{(\chi _2 + \mu _0)^2} + \dfrac{1}{(\chi _2 + \mu _0)} - \dfrac{2 \chi _1}{(\chi _2 + \mu _0)}\nonumber \\&\qquad + \dfrac{2\chi _1^2}{(\chi _2 + \mu _0)} -2 + 3\chi _2\Big ) \Big ] d\chi _1 d\chi _2, \end{aligned}$$

(3.54)

where an overall factor $\dfrac{C_A C_F}{4}$ was set to 1. The integration is over a massive phase space and we already used the energy momentum conservation condition $\chi _1+\chi _2+\chi _3=1$. The integrals are evaluated in Appendix B using results collected in [4, 31]. Notice that we have already extracted the tree-level decay rate, which is given by

$$\begin{aligned} \Gamma _0 = \dfrac{|M_{Hgg}|^2}{32 \pi m_H} = \dfrac{A^2 m_H^3}{8 \pi }. \end{aligned}$$

(3.55)

The final result is

$$\begin{aligned} \begin{aligned}&\Gamma _r (H \longrightarrow gg(g), gq \bar{q}) \\&= \Gamma _0 \dfrac{\alpha _s}{\pi }\left[ C_A \Big (\dfrac{73}{12} + \dfrac{11}{6} \ln (\mu _0 ) \right. \\&\left. + \dfrac{\ln ^2(\mu _0 )}{2} - \dfrac{\pi ^2}{2} \Big ) + N_{F} \Big ( \dfrac{-\ln {(\mu _0)}}{3} - \dfrac{7}{6} \Big )\right] . \end{aligned} \end{aligned}$$

(3.56)

By combining Eqs. 3.26 and 3.56, we get the final decay rate for the process

$$\begin{aligned} \Gamma _T((H \longrightarrow gg(g), gq \bar{q})) = \Gamma _0 \left[ 1 + \dfrac{\alpha _s}{\pi }\left( \dfrac{95}{4} - \dfrac{7}{6}N_{F} \right) \right] , \end{aligned}$$

(3.57)

which, after specializing to QCD, complies with the well-known result of the literature [40].

As a final comment, in IReg it is easy to obtain the decay rate at any other desired renormalization point. This is achieved by leaving the renormalization scale $\lambda ^2$ as a free parameter until the very end. In our case, one rewrites $\log \Big (\dfrac{\lambda ^2}{\mu ^2} \Big ) = \log \Big (\dfrac{\lambda ^2}{m_H^2} \Big ) - \log (\mu _0)$ in 3.25, which produces the result below,

$$\begin{aligned}{} & {} \Gamma _T((H \longrightarrow gg(g), gq \bar{q})) \nonumber \\{} & {} \quad = \Gamma _0 \left[ 1 + \dfrac{\alpha _s}{\pi }\left( \dfrac{95}{4} - \dfrac{7}{6}N_{F} +\dfrac{11 C_A - 2N_{F}}{6} \ln \left( \dfrac{\lambda ^{2}}{m_{H}^{2}} \right) \right) \right] .\nonumber \\ \end{aligned}$$

(3.58)

4 Comparison with dimensional schemes

In this section we intend to compare our results with the ones obtained in the context of dimensional schemes, in a similar way to the discussions presented in [4, 32]. For the results in dimensional schemes, we will mainly use the analysis performed in [24, 30]. In these references, the form factors for the decay were computed up to two-loop order. As discussed there, although one needs to consider a broader set of operators in the effective Lagrangian at two-loop order in all schemes, this issue is relevant for DRED already at one-loop. This is due to the presence of additional, fictitious particles, denoted $\epsilon$-scalars. For our particular example, when using DRED, the one-loop contribution due to light quarks can only be consistently obtained when additional operators are taken into account. As we have shown in the previous section, in the case of IReg the inclusion of $\epsilon$-scalars (or additional operators) is not necessary.

Following [24, 30], the form factors in all dimensional schemes already UV renormalized can be obtained. In CDR the form factor related to the process $H\rightarrow gg$ is given by⁵

$$\begin{aligned} F_{\scriptscriptstyle {\text {CDR}}} (\alpha _s) =\Big (\frac{\alpha _s}{4\pi }\Big ) \Bigg \{C_A \Bigg [ -\frac{2}{\epsilon ^2} -\frac{11}{3\epsilon } +\frac{\pi ^2}{6} \Bigg ] +\frac{2 N_{F}}{3 \epsilon }\Bigg \}+\mathcal {O}(\epsilon ). \end{aligned}$$

(4.1)

In the case of FDH, one needs also to consider the $\epsilon$-scalars when performing the UV renormalization, which amounts to

$$\begin{aligned} F_{\scriptscriptstyle {\text {FDH}}} (\alpha _s)&=F_{\scriptscriptstyle {\text {CDR}}} (\alpha _s) + \Big (\frac{\alpha _s}{4\pi }\Big ) n_\epsilon C_A \Bigg [ 1+\frac{1}{6\epsilon } \Bigg ] +\mathcal {O}(n_\epsilon \epsilon ). \end{aligned}$$

(4.2)

Finally, in DRED one must consider both processes $H\rightarrow gg$ and $H\rightarrow {\tilde{g}}{\tilde{g}}$ where ${\tilde{g}}$ stands for the $\epsilon$-scalars. In this case, one needs the form factor of FDH and also the form factor with $\epsilon$-scalars

$$\begin{aligned} F_{{\tilde{g}}{\tilde{g}}} (\alpha _s)= & {} \Big (\frac{\alpha _s}{4\pi }\Big ) \Bigg \{C_A \Bigg [ -\frac{2}{\epsilon ^2} -\frac{4}{\epsilon } +2+\frac{\pi ^2}{6}-2n_\epsilon \Bigg ] +\frac{N_{F}}{3 \epsilon }\Bigg \}\nonumber \\{} & {} +\mathcal {O}(n_\epsilon \epsilon ). \end{aligned}$$

(4.3)

One should notice that we are already setting the couplings related to $\epsilon$-scalars equal to their corresponding values in usual QCD. This is justified because the UV renormalization was already carried out. In all cases, the form factors are normalized to their tree-level values. Recovering them, the virtual contribution to the decay rate can be readily obtained for each scheme after setting $n_\epsilon =2\epsilon$

$$\begin{aligned} \Gamma _{v}^{\scriptscriptstyle {\text {CDR}}}=&\Gamma _{0}\left\{ 1+\frac{\alpha _s}{\pi }\left[ \dfrac{11}{2}+C_{A}\left( -\frac{1}{\epsilon ^{2}}-\frac{11}{6\epsilon }+\frac{\pi ^{2}}{12}\right) +\frac{N_{F}}{3\epsilon }\right] \right\} \nonumber \\&+\mathcal {O}(\epsilon ) \end{aligned}$$

(4.4)

$$\begin{aligned} \Gamma _{v}^{\scriptscriptstyle {\text {FDH}}}=&\Gamma _{0}\left\{ 1+\frac{\alpha _s}{\pi }\left[ \dfrac{11}{2}+C_{A}\left( -\frac{1}{\epsilon ^{2}}-\frac{11}{6\epsilon }+\frac{\pi ^{2}}{12}+\frac{1}{6}\right) +\frac{N_{F}}{3\epsilon }\right] \right\} \nonumber \\&+\mathcal {O}(\epsilon ) \end{aligned}$$

(4.5)

$$\begin{aligned} \Gamma _{v}^{\scriptscriptstyle {\text {DRED}}}=&\Gamma _{0}\left\{ 1+\frac{\alpha _s}{\pi }\left[ \dfrac{11}{2}+C_{A}\left( -\frac{1}{\epsilon ^{2}}-\frac{11}{6\epsilon }+\frac{\pi ^{2}}{12}\right) +\frac{N_{F}}{3\epsilon }+\frac{N_{F}}{6}\right] \right\} \nonumber \\&+\mathcal {O}(\epsilon ) \end{aligned}$$

(4.6)

We notice that the first term in the square brackets of each equation corresponds to the correction due to the matching of the effective theory to the full SM as given by Eq. 3.3. We can compare these results to Eq. 3.26, and see that the correspondence $\epsilon ^{-1}\rightarrow \log {\mu _{0}}$, $\epsilon ^{-2}\rightarrow \log ^{2}{\mu _{0}}/2$ is fulfilled as first noticed in [4]. Moreover, the result in CDR does not have any finite term (apart from factors of $\pi ^2$ that will be cancelled against the real contribution). The same statement holds true for IReg. For FDH and DRED, on the other hand, there is the appearance of finite terms proportional to $C_{A}$ and $N_{F}$ respectively.

Regarding the real contributions, to the best of our knowledge, the results are not available in the literature for all of the schemes. Nevertheless, the part proportional to $N_{F}$ can be readily obtained. For CDR, there is only one diagram (the diagram on the right of Fig. 2) which produces the following unpolarized amplitude

$$\begin{aligned} |\bar{M}_{q}^{\scriptscriptstyle {\text {CDR}}}|^{2}\propto \dfrac{(s_{13}^2 + s_{12}^2)}{s_{23}} + \dfrac{(d-4)}{2}\dfrac{(s_{13} + s_{12})^2}{s_{23}}. \end{aligned}$$

(4.7)

For FDH, we have the same diagram of CDR and obtain the same result. This happens because the diagram is not one particle irreducible, implying that the internal gluon is regular in the notation of [4]. Therefore, although the external gluon is split into a d-dimensional gluon and a $\epsilon$-scalar, there is no way to produce a diagram with only $\epsilon$-scalars. In the case of DRED this happens, since all vector bosons must be split. The unpolarized amplitude in DRED is given by

$$\begin{aligned}{} & {} |\bar{M}_{q}^{\scriptscriptstyle {\text {DRED}}}|^{2}\propto \dfrac{(s_{13}^2 + s_{12}^2)}{s_{23}} + \dfrac{(d-4)}{2}\dfrac{(s_{13} + s_{12})^2}{s_{23}}\nonumber \\{} & {} \quad +\dfrac{n_\epsilon }{2}\dfrac{(s_{13} + s_{12})^2}{s_{23}}. \end{aligned}$$

(4.8)

It should be noticed that by setting $n_\epsilon =2\epsilon$ one obtains the result of a strictly four-dimensional calculation. One can also compare these results with Eq. 3.51 obtained within IReg which we reproduce below

$$\begin{aligned} |\bar{M}_{q}^{\scriptscriptstyle {\text {IReg}}}|^2 \propto \dfrac{(s_{13}^2 + s_{12}^2)}{s_{23}} + 2 s \mu _{0}\dfrac{(s_{13} + s_{12})^2}{s_{23}^{2}}. \end{aligned}$$

(4.9)

As can be seen, the result in IReg is similar to CDR/FDH in the sense that there is an extra term. In the case of CDR/FDH it comes from extending the physical dimension to d, while in IReg it is encoded in the fictitious mass that we have added for the massless particles. Once the unpolarized amplitude is known in all schemes, one can obtain the part proportional to $N_{F}$ of the real contribution to the decay rate

$$\begin{aligned} \Gamma _{q,r}^{\scriptscriptstyle {\text {CDR/FDH}}}&=\Gamma _{0}\frac{\alpha _s}{\pi }\left[ -\frac{1}{3\epsilon }-\frac{7}{6}\right] N_{F}+\mathcal {O}(\epsilon ) \end{aligned}$$

(4.10)

$$\begin{aligned} \Gamma _{q,r}^{\scriptscriptstyle {\text {DRED}}}&=\Gamma _{0}\frac{\alpha _s}{\pi }\left[ -\frac{1}{3\epsilon }-\frac{4}{3}\right] N_{F}+\mathcal {O}(\epsilon ) \end{aligned}$$

(4.11)

Once again, the result in IReg can be mapped to the one of CDR/FDH after the identification $\epsilon ^{-1}\rightarrow \log {\mu _{0}}$.

5 Conclusion

In conclusion, the decay rate $\Gamma ((H \longrightarrow gg(g), gq \bar{q}))$ at $\alpha _s^3$ order in the strong coupling and large top quark mass limit has been computed using an effective interaction Lagrangian of Higgs to gluons added to the QCD Lagrangian in the framework of the fully quadri-dimensional regularization scheme IReg and compared to dimensional schemes CDR, FDH and DRED. The purpose was twofold. Firstly to achieve not only a full separation of BDI from the UV finite integrals (which IReg accomplishes to arbitrary loop order) but to single out the IR content as well. Secondly to verify whether the additional degrees of freedom associated to epsilon scalars in some of the dimensional schemes have a counterpart in the non-dimensional scheme IReg.

The present calculation provided a proof of concept example involving a non-abelian effective theory Lagrangian with sufficient complexity to allow to infer that crucial steps of the procedure are in compliance with fundamental requirements such as gauge invariance and the removal of IR singularities fulfilling the KLN theorem. In particular the use of a mass regulator in the propagators is adequately implemented in IReg, as well as the renormalization schemes adopted for an effective theory. The latter point in particular is non-trivial, as it involves the use of the method’s renormalization scale relations impacting on the overall cancellation of IR singularities. In addition, by comparing with different dimensional schemes, one concludes that IReg does not require the use of evanescent fields at one loop level.

Appendices

Appendix A: Integrals used in the evaluation of the virtual contribution to $H\rightarrow gg$

In this appendix we collect the integrals needed for the evaluation of the virtual contribution to the process $H\rightarrow gg$. We organize them in terms of the diagrams presented in Fig. 1. The UV finite part of the integrals were evaluated using Package-X [43].

We define $b=\dfrac{i}{(4 \pi )^2}$ and apply the on-shell limit for each integral. We also omit all non-contributing integrals where we already applied the limit $\mu _0 = \mu ^2/m_H^2 \rightarrow 0$.

1.1 A.1 Integrals of diagram $V_1$

The integrals are defined in Eqs. (3.10)–(3.12) and their results read

$$\begin{aligned} I= & {} b \dfrac{\ln ^2{(-\mu _0)}}{4 p_1\cdot p_2} \end{aligned}$$

(A.1)

$$\begin{aligned} I^{\mu }= & {} - b \dfrac{\left[ 2+ \ln {(-\mu _0)}\right] }{2 p_1 \cdot p_2}(p_{1}^{\mu } - p_{2}^{\mu }) \end{aligned}$$

(A.2)

$$\begin{aligned} I^{\mu \nu }= & {} \dfrac{I_{log}(\mu ^2)}{4} g^{\mu \nu } + \dfrac{b}{4 p_1 \cdot p_2} \Big [ \ln (-\mu _0) p_1.p_2 g^{\mu \nu }\nonumber \\{} & {} +3 p_1p_2\,g^{\mu \nu }- 2 p_1^\mu \nonumber \\{} & {} \left( (\ln (-\mu _0)+2) p_{1}^ {\nu }+p_{2}^ {\nu }\right) -2 \ln (-\mu _0) p_{2}^ {\mu }\, p_{2}^{\nu }\nonumber \\{} & {} -2 p_{1} ^{\nu } p_{2}^ {\mu }-4 p_{2}^ {\mu }\, p_{2}^ {\nu } \Big ] \end{aligned}$$

(A.3)

$$\begin{aligned} I^{\alpha \mu \nu }= & {} \dfrac{I_{log}(\mu ^2)}{12} ((p_1-p_2)^\nu g^{\alpha \mu }\nonumber \\{} & {} + (p_1-p_2)^\mu g^{\alpha \nu } + (p_1 - p_2)^\alpha g^{\mu \nu } )\nonumber \\{} & {} + b \Big [ \frac{1}{18 p_1 \cdot p_2} 4 (-p_1.p_2 p_{2}^ \mu g^{\alpha \nu }\nonumber \\{} & {} + p_{1}^ \nu p_1 \cdot p_2 g^{\alpha \mu }-p_1 \cdot p_2 p_{2 }^\nu g^{\alpha \mu } -p_1 \cdot p_2 p_{2 }^\alpha g^{\mu \nu })\nonumber \\{} & {} + 3 \ln (-\mu _0)(- p_1 \cdot p_2 p_{2 }^\mu g^{\alpha \nu }+ p_{1 } ^\nu p_1 \cdot p_2 g^{\alpha \mu }\nonumber \\{} & {} - p_1 \cdot p_2 p_{2 }^\nu g^{\alpha \mu }) \nonumber \\{} & {} + p_{1 }^\mu ((3 \ln (-\mu _0)+4) p_1 \cdot p_2 g^{\alpha \nu }-3 p_{1 }^\nu p_{2}^{\alpha }\nonumber \\{} & {} +3 p_{2}^{\alpha }p_{2}^{\nu })\nonumber \\{} & {} -3 \ln (-\mu _0) p_1 \cdot p_2 p_{2 }^\alpha g^{\mu \nu }+p_{1 }^\alpha (4 p_1 \cdot p_2 g^{\mu \nu }\nonumber \\{} & {} +3 \ln (-\mu _0) p_1 \cdot p_2 g^{\mu \nu }-3 p_{1 }^\nu p_{2 }^\mu \nonumber \\{} & {} -p_{1 }^\mu ((6 \ln (-\mu _0)+13) p_{1 }^\nu +3 p_{2}^{\nu })+3 p_{2 }^\mu p_{2 }^\nu )\nonumber \\{} & {} +3 p_{1 }^\nu p_{2 }^\alpha p_{2 }^\mu \nonumber \\{} & {} +13 p_{2 }^\alpha p_{2 }^\mu p_{2 }^\nu +6 \ln (-\mu _0) p_{2 }^\alpha p_{2}^{\mu } p_{2}^{\nu }) \Big ] \end{aligned}$$

(A.4)

$$\begin{aligned} I_2= & {} I_{log}(\mu ^2) + b( 2 + \ln {(-\mu _0)} ) \end{aligned}$$

(A.5)

$$\begin{aligned} I_2^{\nu }= & {} \dfrac{I_{log}(\mu ^2)}{2}(p_1^\nu - p_2^\nu ) + \dfrac{b}{2}(2 +\ln {(-\mu _0)})(p_{1 }^\nu - p_{2 }^\nu )\nonumber \\ \end{aligned}$$

(A.6)

$$\begin{aligned} I_2^{\mu \nu }= & {} \dfrac{I_{quad}(\mu ^2)}{2} g^{\mu \nu } + \dfrac{I_{log}(\mu ^2)}{3} (p_2^\mu p_2^\nu + p_1^\mu p_1^\nu ) \nonumber \\{} & {} - \dfrac{I_{log}(\mu ^2)}{6}(p_1 \cdot p_2 g^{\mu \nu } + p_1^\mu p_2^\nu + p_1^\nu p_2^\mu ) \nonumber \\{} & {} + \frac{b}{36} \Big [ \Big ( \ln (-\mu _0) (12(p_{1 }^\mu p_{1}^{\nu } + p_{2}^{\mu } p_{2}^ {\nu })\nonumber \\{} & {} -6 (p_{1}^ {\nu } p_{2}^ {\mu } + p_{2}^{\mu } p_{2}^ {\nu } + p_1 \cdot p_2 \,g^{\mu \nu }) \Big )\nonumber \\{} & {} +26( p_{1}^{ \mu } p_{1}^{\nu }+p_{2}^{\mu }\, p_{2}^{\nu })-10 (p_{1}^{\nu }\, p_{2}^{\mu }+p_{1}^{\mu } p_{2}^{\nu })\nonumber \\{} & {} \quad - 16 p_1 \cdot p_2 \,g^{\mu \nu } ) \Big ] \end{aligned}$$

(A.7)

$$\begin{aligned} I_4= & {} I_{quad}(\mu ^2) - I_{log}(\mu ^2) p_1 \cdot p_2 \nonumber \\{} & {} - b (p_1 \cdot p_2) (2 + \ln {(-\mu _0)}) \end{aligned}$$

(A.8)

1.2 A.2 Integrals of diagram $V_2$

$$\begin{aligned}{} & {} \int _k \dfrac{1}{(k-p_1-p_2)^2} \nonumber \\{} & {} \quad = \int _k \dfrac{1}{k^{2}-2k\cdot (p_1+p_2) + 2p_{1}\cdot p_{2}}\nonumber \\{} & {} \quad = I_{quad}(\mu ^2) + 2b p_1 \cdot p_2 \ln (-\mu _0) \end{aligned}$$

(A.9)

$$\begin{aligned}{} & {} \int _k \dfrac{k^\mu }{k^2 (k-p_1-p_2)^2} = \int _k \dfrac{k^\mu }{k^{2}(k^{2}-2k\cdot (p_1+p_2) + 2p_{1}\cdot p_{2})} \nonumber \\{} & {} \quad = \dfrac{I_{log}(\mu ^2)}{2} (p_1 + p_2)^\mu + \frac{b}{2} (\ln (-\mu _0)+2) \left( p_{1}^{\mu }+p_{2}^{\mu }\right) \end{aligned}$$

(A.10)

$$\begin{aligned}{} & {} \int _k \dfrac{k^\mu k^\nu }{k^2(k-p_1-p_2)^2} \nonumber \\{} & {} \quad = \int _k \dfrac{k^\mu k^\nu }{k^{2}(k^{2}-2k\cdot (p_1+p_2) + 2p_{1}\cdot p_{2})} \nonumber \\{} & {} \quad = \dfrac{I_{quad}(\mu ^2)}{2} g^{\mu \nu } - \dfrac{I_{log}(\mu ^2)}{2} p_1 \cdot p_2 g^{\mu \nu } \nonumber \\{} & {} \qquad + \dfrac{I_{log}(\mu ^2)}{3} (g^{\mu \nu } p_1 \cdot p_2 + (p_1 + p_2)^\mu (p_1 + p_2)^\nu )\nonumber \\{} & {} \qquad + b \Big [ \frac{1}{18} p_1\cdot p_2 g^{\mu \nu }+\frac{1}{3} \ln (-\mu _0)\big ( p_1 \cdot p_2 g^{\mu \nu }+ p_{1}^{\mu } p_{1}^{\nu }\nonumber \\{} & {} \qquad + p_{1}^{\nu } p_{2}^{\mu }+ p_{1}^{\mu } p_{2}^{\nu }+ p_{2}^{\mu } p_{2}^{\nu }\big )\nonumber \\{} & {} \qquad +\frac{2}{9} \big (p_{1}^{\mu } p_{1}^{\nu } + p_{1}^{\nu } p_{2}^{\mu } + p_{1 }^{\mu } p_{2}^{\nu } + p_{2}^{\mu } p_{2}^{\nu } \big ) \Big ] \end{aligned}$$

(A.11)

1.3 A.3 Integrals of diagrams $V_3$ and $V_4$

$$\begin{aligned}{} & {} \int _k \dfrac{1}{(k-p_2)^2} = \int _k \dfrac{1}{k^{2}-2k.p_2} = I_{quad}(\mu ^2) \end{aligned}$$

(A.12)

$$\begin{aligned}{} & {} \int _k \dfrac{1}{k^2(k-p_2)^2} = \int _k \dfrac{1}{k^{2}(k^{2}-2k.p_2)} = I_{log}(\mu ^2) \end{aligned}$$

(A.13)

$$\begin{aligned}{} & {} \int _k \dfrac{k^\mu }{k^2(k-p_2)^2} = \int _k \dfrac{k^\mu }{k^{2}(k^{2}-2k.p_2)} \nonumber \\{} & {} \qquad \qquad = \dfrac{I_{log}(\mu ^2)}{2} p_2^\mu \end{aligned}$$

(A.14)

$$\begin{aligned}{} & {} \int _k \dfrac{k^\mu k^\nu }{k^2(k-p_2)^2} = \int _k \dfrac{k^\mu k^\nu }{k^{2}(k^{2}-2k.p_2)}\nonumber \\{} & {} \qquad \qquad =\dfrac{I_{quad}(\mu ^2)}{2} g^{\mu \nu } + \dfrac{I_{log}(\mu ^2)}{3} p_2^\mu p_2^\nu \end{aligned}$$

(A.15)

Appendix B: Integrals used in the evaluation of the real contribution to $H\rightarrow gg$

The phase space integral is found to be

$$\begin{aligned} \rho = \dfrac{q_0^2}{(4\pi )^3} \int _{\chi _1^{-}}^{\chi _1^{+}} d\chi _1 \int _{\chi _2^{-}}^{\chi _2^{+}} d\chi _2 \end{aligned}$$

(B.1)

where the boundaries for $\chi _2^\pm$ are given by

$$\begin{aligned} \chi _2^\pm = \dfrac{1 - \chi _1}{2} \pm \sqrt{\dfrac{(\chi _1 - 3\mu _0)[(1-\chi _1)^2-4\mu _0]}{4(\chi _1 + \mu _0)}}, \end{aligned}$$

(B.2)

and the boundaries for $\chi _1^\pm$ are

$$\begin{aligned} \chi _1^{+} = 1 - 2\sqrt{\mu _0},\quad \text{ and }\quad \chi _1^{-} = 3\mu _0. \end{aligned}$$

(B.3)

Finally, the integrals were evaluated using the following equations, [4, 31]

$$\begin{aligned} I(s) = \int d \chi _1 d \chi _2 \dfrac{1}{(\chi _1+\mu _0)(\chi _2+\mu _0)} \end{aligned}$$

(B.4)

and

$$\begin{aligned} J_p(s) = \int d\chi _1 d\chi _2 \dfrac{\chi _1^p}{(\chi _2+\mu _0)} \end{aligned}$$

(B.5)

with $p\ge 0$. Using our limit of integrations, the integrals are evaluated to be

$$\begin{aligned} I(s) = \dfrac{\ln ^2(\mu _0)-\pi ^2}{2} \end{aligned}$$

(B.6)

and

$$\begin{aligned} J_p(s) =&- \dfrac{1}{p+1} \ln (\mu _0) + \int _0^1 d\chi _1 \chi _1^p [\ln (\chi _1) + 2 \ln (1-\chi _1)]\nonumber \\ =&- \dfrac{1}{p+1} \ln (\mu _0) - \dfrac{1}{p+1} \Big [\dfrac{1}{p+1} +2 \sum _{n=1}^{p+1} \dfrac{1}{n} \Big ]. \end{aligned}$$

(B.7)