Evaluating two-loop non-planar master integrals for Higgs + jet production with full heavy-quark mass dependence

TL;DR

This paper analytically computes non-planar two-loop master integrals with full heavy-quark mass dependence, crucial for precise Higgs + jet production predictions at NNLO and NLO levels in QCD.

Contribution

It provides a canonical form basis for differential equations of these integrals and offers solutions up to weight 4, including elliptic sectors, enabling accurate numerical evaluations.

Findings

Analytic solutions for integrals up to weight 4.

Elliptic sector solutions via power series.

High-precision numerical results across kinematic regions.

Abstract

We present the analytic computation of a family of non-planar master integrals which contribute to the two-loop scattering amplitudes for Higgs plus one jet production, with full heavy-quark mass dependence. These are relevant for the NNLO corrections to inclusive Higgs production and for the NLO corrections to Higgs production in association with a jet, in QCD. The computation of the integrals is performed with the method of differential equations. We provide a choice of basis for the polylogarithmic sectors, that puts the system of differential equations in canonical form. Solutions up to weight 2 are provided in terms of logarithms and dilogarithms, and 1-fold integral solutions are provided at weight 3 and 4. There are two elliptic sectors in the family, which are computed by solving their associated set of differential equations in terms of generalized power series. The resulting series may be truncated to obtain numerical results with high precision. The series solution renders the analytic continuation to the physical region straightforward. Moreover, we show how the series expansion method can be used to obtain accurate numerical results for all the master integrals of the family in all kinematic regions.

Tables

Equations

$$\displaystyle f^{ade}\,\text{tr}\big{(}T^{d}T^{b}T^{e}T^{c}\big{)}=0\,.$$

$$I_{a_1{a}_2{a}_3{a}_4{a}_5{a}_6{a}_7{a}_8{a}_9}$$

$$P_1$$

$$P_2$$

$$P_3$$

$$s$$

$$\frac{\partial }{\partial s_i}\vec{\text{I}}={\mathbf{M}}_{s_i}\vec{\text{I}},$$

$$\frac{\partial }{\partial s_i}\vec{B}=\left[\left({\partial }_{s_i}\mathbf{T}\right){\mathbf{T}}^{-1}+\mathbf{T}{\mathbf{M}}_{s_i}{\mathbf{T}}^{-1}\right]\vec{B}.$$

$$\left({\partial }_{s_i}\mathbf{T}\right){\mathbf{T}}^{-1}+\mathbf{T}{\mathbf{M}}_{s_i}{\mathbf{T}}^{-1}=ϵ{\mathbf{A}}_{s_i},$$

$$d\vec{B}=ϵd\tilde{\mathbf{A}}\vec{B},$$

$${\tilde{\mathbf{A}}}_1:=\int {\mathbf{A}}_{s_1}d{s}_1,$$

$$\displaystyle\tilde{\mathbf{A}}_{i}:=\int\bigg{(}\mathbf{A}_{s_{i}}-\partial_{s_{i}}\sum_{j=1}^{i-1}\tilde{\mathbf{A}}_{j}\bigg{)}ds_{i}\,,~{}~{}~{}~{}i=2,...,4\,.$$

$$\frac{a+\text{Alg}}{a-\text{Alg}}$$

$$\log \left(a\pm b\sqrt{c}\right)=\frac{1}{2}\left[\pm \log \left(\frac{a/b+\sqrt{c}}{a/b-\sqrt{c}}\right)+\log \left(a^2-b^{2}c\right)\right](\mathrm{mod}i\pi ).$$

$$\displaystyle\log\left(a+b\sqrt{c}+d\sqrt{e}\right)=\frac{1}{2}\bigg{[}$$

$$\displaystyle+\log\left(\frac{a+b\sqrt{c}+d\sqrt{e}}{a-b\sqrt{c}+d\sqrt{e}}\right)\bigg{]}\,\pmod{i\pi}\,.$$

$$\log \left(\frac{a+b\sqrt{c}+d\sqrt{e}}{a-b\sqrt{c}+d\sqrt{e}}\right)$$

$$\displaystyle-\log\Big{[}\left(a-b\sqrt{c}+d\sqrt{e}\right)\left(a+b\sqrt{c}-d\sqrt{e}\right)\Big{]}\,\pmod{i\pi}\,.$$

$$\left(a+b\sqrt{c}-d\sqrt{e}\right)\left(a+b\sqrt{c}+d\sqrt{e}\right)=a^2+b^{2}c-d^{2}e+2ab\sqrt{c}.$$

$$\displaystyle\log\bigg{(}a+b\sqrt{c}+d\sqrt{e}\bigg{)}$$

$$-\frac{1}{4}\log \left(\frac{a^2-b^{2}c-d^{2}e+2bd\sqrt{c}\sqrt{e}}{a^2-b^{2}c-d^{2}e-2bd\sqrt{c}\sqrt{e}}\right)$$

$$+\frac{1}{4}\log \left(\frac{a^2-b^{2}c+d^{2}e+2ad\sqrt{e}}{a^2-b^{2}c+d^{2}e-2ad\sqrt{e}}\right)$$

$$+\frac{1}{4}\log \left(a^4-2a^2\left(b^{2}c+d^{2}e\right)+{\left(b^{2}c-d^{2}e\right)}^2\right)(\mathrm{mod}i\pi ).$$

$$\mathcal{E}:m^2>0\&p_4^2<2m^2\&p_4^2-4m^2<t<0\&-2m^2+p_4^2-t<s<2m^2,$$

$${\tilde{\mathcal{R}}}_{+,+,\dots ,+}$$

$${\tilde{\mathcal{R}}}_{-,+,\dots ,+}$$

$$⋮$$

$${\tilde{\mathcal{R}}}_{-,-,\dots ,-}$$

$$l_{25}=\frac{1+\frac{\sqrt{4m^2-s}}{\sqrt{-s}}}{1-\frac{\sqrt{4m^2-s}}{\sqrt{-s}}}.$$

$$(4m^2-s\ge 0\&-s\ge 0)||(4m^2-s\le 0\&-s\le 0),$$

$$\mathcal{R}:t\le$$

$$\displaystyle\left(t<s\,\&\,\frac{4m^{2}(s+t)-st}{4m^{2}}\leq p_{4}^{2}\leq\frac{-4m^{2}t+s^{2}+st}{s}\right)\bigg{)}\,\&\,m^{2}\geq 0\,.$$

$$\mathcal{S}\left(B_i^{(k)}\right)=\sum _j\mathcal{S}\left(B_j^{(k-1)}\right)\otimes {\tilde{\mathbf{A}}}_{ij}.$$

$$\{1,2,3,4,5,6,7,8,9,10,11,12,16,20,23,26,33,35,38,40,49\}.$$

$$\displaystyle\begin{array}[]{rrrrrrrr}l_{1}>0\,,&l_{2}<0\,,&l_{3}<0\,,&l_{4}<0\,,&l_{5}<0\,,&l_{6}<0\,,&l_{7}<0\,,&l_{8}<0\,,\\ l_{9}>0\,,&l_{10}>0\,,&l_{11}>0\,,&l_{13}>0\,,&l_{25}<0\,,&l_{26}<0\,,&l_{27}<0\,,&l_{28}>0\,,\\ l_{29}<0\,,&l_{39}>0\,,&l_{40}>0\,,&l_{43}>0\,,&l_{44}>0\,,&l_{46}>0\,,&l_{48}>0\,,&l_{49}>0\,,\\ l_{53}>0\,,&l_{54}>0\,,&l_{55}>0\,,&l_{56}>0\,,&l_{60}>0\,,&l_{61}>0\,.\end{array}$$

$${\text{Li}}_2\left(\pm l_i{l}_j\right),{\text{Li}}_2\left(\pm \frac{l_i}{l_j}\right),{\text{Li}}_2\left(\pm \frac{1}{l_i{l}_j}\right)\text{for }{l}_i,l_j\in {\mathcal{A}}_2\cup \{l_{33},l_{38},l_{41}\},$$

Full text

††preprint: MPP-2019-157

  TCDMATH 19-16

1 Introduction

After the discovery of the Higgs boson Aad:2012tfa ; Chatrchyan:2012xdj at the Large Hadron Collider (LHC) at CERN, the LHC physics program has been centered around measuring the properties, couplings and quantum numbers of the Higgs boson, and looking for footprints of New Physics (NP) effects. A process which may be useful in this respect is the production of a Higgs boson in association with a jet. At the LHC, the main production mode of the Higgs boson is via gluon-gluon fusion. However, the Higgs boson does not couple directly to the gluons, the interaction being mediated by a heavy-quark loop. Thus, leading-order production requires the evaluation of one-loop amplitudes, the next-to-leading order (NLO) QCD corrections involve the evaluation of two-loop amplitudes, and so on.

In many NP models, the high $p_{T}$ tail of the Higgs $p_{T}$ distribution is sensitive to deviations of the Higgs-top coupling Harlander:2013oja ; Banfi:2013yoa ; Azatov:2013xha ; Grojean:2013nya ; Schlaffer:2014osa ; Buschmann:2014twa ; Dawson:2014ora ; Buschmann:2014sia ; Ghosh:2014wxa ; Dawson:2015gka ; Langenegger:2015lra ; Azatov:2016xik ; Grazzini:2016paz ; Grazzini:2018eyk ; Gorbahn:2019lwq from the Standard Model (SM) coupling. In the full theory, Higgs production in association with one jet Ellis:1987xu and the Higgs $p_{T}$ distribution Kauffman:1991jt are known only at leading order; at NLO they are known Jones:2018hbb when including only the dependence on the top-quark mass and neglecting the bottom-quark mass. The computations are made more tractable by using the Higgs Effective Field Theory (HEFT), i.e. by replacing the loop-mediated Higgs-gluon coupling by a tree-level effective coupling, which lowers by one loop the amplitudes to be computed. In the HEFT, Higgs production in association with one jet Boughezal:2015dra ; Boughezal:2015aha and the Higgs $p_{T}$ distribution Chen:2016zka are known at next-to-next-to-leading order (NNLO) accuracy. However, for Higgs production in association with one jet or for the Higgs $p_{T}$ distribution one can show that the HEFT can be applied when the Higgs mass is smaller than the heavy-quark mass, $m_{H}\lesssim m_{Q}$, and when the jet or Higgs transverse momenta are smaller than the heavy-quark mass, $p_{T}\lesssim m_{Q}$ Baur:1989cm ; deFlorian:1999zd , by using the leading-order results Ellis:1987xu ; Kauffman:1991jt as a benchmark. In the $p_{T}\lesssim m_{Q}$ region, the HEFT approximation may be improved by including the top-bottom interference at NLO, which is estimated by interfering a top-quark loop computed in HEFT with a bottom-quark loop computed as an expansion in a small bottom-quark mass Lindert:2017pky .

When the jet or Higgs transverse momenta are of the order or larger than the heavy-quark mass, $p_{T}\gtrsim m_{Q}$, the HEFT is not a viable approximation. In that region, approximate computations at NLO exist Lindert:2018iug ; Neumann:2018bsx , based on the two-loop amplitudes for Higgs plus three partons, computed in the limit $p_{T}\gg m_{Q}$ Kudashkin:2017skd , and, as outlined above, a numerical NLO computation in the full theory Jones:2018hbb , which includes only the dependence on the top-quark mass.

To this date, no NLO computation of the whole $p_{T}$ spectrum for Higgs production in association with one jet, where quark-mass effects are included for all flavours, is available. In Ref. Bonciani:2016qxi an analytic computation is presented of the planar master integrals which contribute to the two-loop scattering amplitudes for Higgs+jet production, with full heavy-quark mass dependence. In this paper, we compute analytically one of the two remaining families of non-planar master integrals.

The remainder of the paper is organized as follows. In Section 2 we review the integral families which are necessary to compute the two-loop Higgs+jet amplitudes in QCD with full heavy quark mass dependence, and we set up the notation for the non-planar family of integrals that is computed in this paper. In Section 3 we present the computation of the master integrals of the polylogarithmic sectors. We briefly review some properties of the canonical basis, and then we discuss how to obtain a minimal-size alphabet for the canonical master integrals. We then discuss how to obtain a region where the canonical basis and the alphabet is manifestly real-valued, and we provide results in this region in terms of logarithms and dilogarithms up to weight 2, and 1-fold integrals over polylogarithms up to weight 4. In Section 4 we describe how (generalized) power series expansions are obtained for the polylogarithmic and the elliptic sectors. Series expansions are applied to obtain results below and above the physical threshold. Lastly in Section 6 we conclude and summarize the results.

2 The Higgs + jet integral families

The evaluation of the two-loop QCD contribution to $gg\rightarrow gH$, for instance with the package FeynArts Hahn:2000kx , results in $\mathcal{O}(300)$ Feynman diagrams. After performing the Dirac algebra, these diagrams get expressed in terms of $\mathcal{O}(20000)$ scalar Feynman integrals. Using symmetries and IBP identities, the integrals can all be expressed in terms of a basis of master integrals which are members of one of the seven-propagator integral families depicted on Figure 1, that we name with consecutive letters from A to G. The quark channels $q\bar{q}\rightarrow gH$, $qg\rightarrow qH$, and $\bar{q}g\rightarrow\bar{q}H$ map into a subset of the same seven families.

Of these families, the planar families A to D were computed in Ref. Bonciani:2016qxi . The non-planar family E turns out to not contribute to the amplitude, as each diagram in that family gets multiplied by the vanishing color-structure,

[TABLE]

Thus, only non-planar families F and G need still to be computed. This paper is devoted to the integrals in family F, while family G is postponed to a future publication.

2.1 Definition of the integral family

The integral family we are considering is depicted in Figure 2, and is given by

[TABLE]

where $\gamma_{E}=-\Gamma^{\prime}(1)$ is the Euler-Mascheroni constant, and where

[TABLE]

In the above, $P_{1}$-$P_{7}$ are propagators while $P_{8}$ and $P_{9}$ are numerator factors, so $a_{8}$ and $a_{9}$ are restricted to non-positive integers. The kinematics is such that $p_{1}^{2}=p_{2}^{2}=p_{3}^{2}=0$, and additionally,

[TABLE]

where $p_{4}^{2}$ is the squared Higgs-mass, $p_{4}^{2}=m_{H}^{2}$, and $m^{2}$ is the squared mass of the massive quark coupling to the Higgs.

The family contains 73 master integrals, defined in Appendix B. It is shown in Appendix C that the maximal cuts associated to integral sectors 126 and 127 are elliptic integrals. We refer to these sectors as the elliptic sectors, while the remaining sectors will be referred to as the polylogarithmic sectors. The starting point for solving all master integrals will be the differential equation method. To obtain closed form systems of differential equations for the master integrals, we have to perform IBP-reductions, for which we used the programs FIRE Smirnov:2008iw ; Smirnov:2013dia ; Smirnov:2014hma ; Smirnov:2019qkx and Kira Maierhoefer:2017hyi ; Maierhofer:2018gpa .

3 The polylogarithmic sectors

This section contains the computation of the master integrals of the polylogarithmic sectors, which relies on expressing the integrals in terms of a canonical basis. The integration of the polylogarithmic sectors is performed at weight 2 in terms of manifestly real-valued logarithms and dilogarithms in a sub-region $\mathcal{R}$ of the Euclidean region, where the canonical basis integrals and the alphabet are real-valued. The weight 3 and 4 expressions are then obtained as 1-fold integrals over the weight 2 result.

3.1 Canonical basis

Suppose that $\vec{\text{I}}$ is a set of master integrals, each of which is given by a Feynman integral times an overall normalization factor that does not depend on the external scales. We may then write down a closed form system of first order linear differential equations,

[TABLE]

where $\{s_{i}\}$ denotes the set of external scales, and where the matrix $\mathbf{M}_{s_{i}}$ consists of rational functions of the external scales and the dimensional regulator $\epsilon$. We parametrize our kinematics by the external scales $\{s_{i}\}=\{s,t,m^{2},p_{4}^{2}\}$ defined in Section 2.1.

One may perform a change of basis $\vec{B}=\mathbf{T}\vec{\text{I}}$, to find a new system of differential equations,

[TABLE]

It has been conjectured in Ref. Henn:2013pwa that for some matrix $\mathbf{T}$, it is possible to put the system of differential equations in a canonical form. In this form all dependence on the dimensional regulator $\epsilon$ is in an overall prefactor,

[TABLE]

where the matrices $\mathbf{A}_{s_{i}}$ do not depend on $\epsilon$. In the polylogarithmic case it is conjectured that $\mathbf{T}$ is a matrix of rational or algebraic functions, and that it is possible to combine the partial derivative matrices into a $d\log$-form. The system of differential equations then becomes,

[TABLE]

where the entries of the matrix $\tilde{\mathbf{A}}$ are $\mathbb{Q}$-linear combinations of logarithms of rational or algebraic functions, and where in particular $\partial_{s_{i}}\tilde{\mathbf{A}}=\mathbf{A}_{s_{i}}$. The canonical basis of the polylogarithmic sectors consists of 65 integrals, which are given in Appendix B.

3.1.1 Deriving the $\tilde{\mathbf{A}}$-matrix

Starting from equations (6), (7), the matrix $\tilde{\mathbf{A}}$ may be found by defining the following matrices,

[TABLE]

The matrix $\tilde{\mathbf{A}}$ is then obtained as the sum of the matrices $\tilde{\mathbf{A}}_{i}$. The ordering of the variables can be chosen arbitrarily. Hence, trying this procedure for different orderings might lead to simpler expressions. We benefited from the fact that $\tilde{\mathbf{A}}_{i}$ should not depend on the variables $s_{j}$, with $j<i$ and used numerical evaluations on those variables for which the simplification was not explicit in each of the intermediate steps. Further simplifications were done by rationalizing different roots, changing the differential equation matrices in accordance with the corresponding chain rule of the substitution and replacing back the square roots in the end. In this way we obtained $\tilde{\mathbf{A}}$ as a $\mathbb{Q}$-linear combination of logarithms, for which the correctness of all the partial derivatives was checked analytically.

3.1.2 Simplifying the alphabet

The form of $\tilde{\mathbf{A}}$ that is found by repeatedly integrating and projecting out components, in the manner of Eq. (9), contains many linearly dependent combinations of logarithms. It is desirable to express these logarithms in terms of a linearly independent set. We will generally use the term ‘letter’  to refer to the arguments of the logarithms, and ‘alphabet’ to a set of such arguments. However, when we speak about linear (in-)dependence of a set of letters we refer to the linear (in-)dependence of the logarithms of the letters. For example, if $\mathcal{A}^{\text{example}}=\{a,b\}$, then we let $\text{Span}_{\mathbb{Q}}\mathcal{A}^{\text{example}}=\text{Span}_{\mathbb{Q}}\{\log a,\log b\}$.

We seek to express $\tilde{\mathbf{A}}$ in terms of a linearly independent alphabet, which is also manifestly symmetric under flipping signs of the square roots. As a first step we enumerate the irreducible factors of the arguments of the logarithms that appear in $\tilde{\mathbf{A}}$. An alphabet consisting of these letters spans the space of logarithms in $\tilde{\mathbf{A}}$, but is still overcomplete. Let us denote this overcomplete set by $\mathcal{A}^{\text{oc}}$. We may obtain a linearly independent basis of $\text{Span}_{\mathbb{Q}}\mathcal{A}^{\text{oc}}$ by starting with an empty set $\emptyset=\mathcal{A}^{\text{idp}}$, and iteratively adding to $\mathcal{A}^{\text{idp}}$ a letter from $\mathcal{A}^{\text{oc}}$ that is independent of the ones already contained in $\mathcal{A}^{\text{idp}}$. After adding a letter to $\mathcal{A}^{\text{idp}}$ we remove from $\mathcal{A}^{\text{oc}}$ the elements that are contained in $\text{Span}_{\mathbb{Q}}\mathcal{A}^{\text{idp}}$. For simplicity we let each choice of new letter be the one that is of smallest (notational) size in $\mathcal{A}^{\text{oc}}$.

After iterating this procedure we are left with a linearly independent alphabet $\mathcal{A}^{\text{idp}}$, such that all elements of $\tilde{\mathbf{A}}$ lie in $\text{Span}_{\mathbb{Q}}\mathcal{A}^{\text{idp}}$. Next we seek to find an alphabet whose letters are manifestly symmetric under flipping signs of their square roots. In particular, we seek to write the algebraic letters in the form,

[TABLE]

where $a$ is some rational function of the external scales, and Alg denotes a product of square roots, and inverses of square roots. Flipping a sign of a square root in Alg sends the letter to its reciprocal, and hence the logarithm of the letter to its negative. Note that there is the freedom to combine the square roots in Alg together, but we prefer to split up the roots in terms of irreducible factors. This can be done unambiguously if we work in a region where each irreducible radicand has a definite sign, which is described in more detail in Section 3.2.

To find an alphabet $\mathcal{A}^{\text{sym}}$ spanning the elements of $\tilde{\mathbf{A}}$, and in which all algebraic letters are symmetric of the form of Eq. (10), we encountered two cases. The first case consists of letters of the type $\log\left(a+b\,\text{Alg}\right)$, which for simplicity we denote by $\log\left(a+b\sqrt{c}\right)$, as the following arguments go through in the same way for a combination of square roots or a single square root. One may substitute letters of this type in terms of symmetric ones by considering the following relation, obtained from multiplying by the conjugate,

[TABLE]

If the term $a^{2}-b^{2}c$ contains irreducible factors that are linearly independent of the original alphabet, we add these to the alphabet as well.

The second case consists of letters of the form $\left(a+b\sqrt{c}+d\sqrt{e}\right)$, where $\sqrt{c}$ and $\sqrt{e}$ may again denote a combination of square roots. Multiplying by the conjugate with respect to the first square root, i.e. $\left(a-b\sqrt{c}+d\sqrt{e}\right)$, leads to the following relation,

[TABLE]

Furthermore, multiplying the fraction by the term $1=\left(a+b\sqrt{c}-d\sqrt{e}\right)/\left(a+b\sqrt{c}-d\sqrt{e}\right)$, composed of the conjugate with respect to the second square root, yields the relation,

[TABLE]

After expanding the products of conjugate terms, their algebraic dependence is captured in a single term of the sum. For example,

[TABLE]

Such terms can be dealt with in the same manner as Eq. (11). Putting everything together leads to the final relation,

[TABLE]

It can be verified that the letters on the right-hand side are in fact sufficient to rewrite every term of the form $\log\left(\pm a\pm b\sqrt{c}\pm d\sqrt{e}\right)$, where the plus and minus signs may differ from each other.

At this point the alphabet $\mathcal{A}^{\text{sym}}$ consists of linearly independent letters with manifest symmetry properties under flip of sign of any square root, and the alphabet covers the entries of $\tilde{\mathbf{A}}$, i.e. $\text{Span}_{\mathbb{Q}}\{\tilde{\mathbf{A}}_{ij}\}\subseteq\text{Span}_{\mathbb{Q}}\mathcal{A}^{\text{sym}}$. The alphabet is still larger than necessary when $\text{Span}_{\mathbb{Q}}\{\tilde{\mathbf{A}}_{ij}\}$ is a proper subspace. In this case some letters only appear in fixed combinations in $\{\tilde{\mathbf{A}}_{ij}\}$.

Indeed, we found at this stage an alphabet $\mathcal{A}^{\text{sym}}$ that contains 75 letters, while the rank of the vector space spanned by the entries of the canonical matrix $\{\tilde{\mathbf{A}}_{ij}\}$ is equal to 69. To reduce the alphabet to 69 independent letters, we sorted $\{\tilde{\mathbf{A}}_{ij}\}$ by (notational) complexity and picked out the first 69 independent entries. We could identify that 12 letters of $\mathcal{A}^{\text{sym}}$ only appear together in pairs of two in these entries. Combining these pairs into 6 letters yields the final alphabet $\mathcal{A}$ which is written out fully in Appendix A.

Letters $l_{63},l_{64},l_{65},l_{67},l_{68}$ and $l_{69}$ are the result of combining pairs of letters of $\mathcal{A}^{\text{sym}}$. Each pair contains the same square roots, so that the algebraic dependence of the ‘combined’  letters is still symmetric: changing the sign of a square root sends the letter to its reciprocal, and hence changes the overall sign of its logarithm.

3.2 A manifestly real region

It is convenient to work in a kinematic region in which the integrals are real-valued and free of branch cuts. Such a region can be found for the Feynman integrals in Eq. (2) by requiring that their second Symanzik polynomial is positive in the whole integration domain, i.e. $\mathcal{F}>0$. It is sufficient to consider the scalar integral with maximal number of propagators, which provides the region

[TABLE]

which we refer to as the Euclidean region. The canonical basis integrals may be complex-valued in the Euclidean region as they are algebraic combinations of Feynman integrals, and the square roots in the prefactors may be evaluated at negative argument. The alphabet also contains these square roots, and we seek to work in a region where the letters are manifestly real-valued as well.

In Appendix B we label 15 square roots whose radicands are irreducible polynomials. These roots appear in the canonical basis and the alphabet of the polylogarithmic sectors, and are given explicitly in Eq. (79). As the roots only appear in certain fixed combinations, they may be combined together into just 10 roots. We decided to keep them separated. The letters are real-valued when the ratios of roots which they contain are real-valued. To find a region where this is the case for all letters, let us first decompose the phase space into $2^{15}$ regions $\tilde{\mathcal{R}}_{\sigma_{1},\ldots,\sigma_{15}}$, depending on the signs of the radicands of the 15 roots in Eq. (79),

[TABLE]

Next, we filter out regions where (some of) the letters are complex valued. For example, consider the letter

[TABLE]

Looking at the ratio of two roots in the letter, we see that the letter is real-valued if and only if

[TABLE]

since the external scales are real-valued and with any other assignment of signs to the radicands, the ratio evaluates to an imaginary number. Note that this is also the case if the roots in $l_{25}$ had been combined together. Hence, we filter out all the regions with $(\sigma_{2},\sigma_{7})=(+-)$ or $(\sigma_{2},\sigma_{7})=(-+)$.

After selecting all regions where the alphabet is manifestly real-valued, we consider their intersection with the Euclidean region. Just one region has a non-empty intersection, which is the one where all the radicands are non-negative: $\mathcal{R}=\tilde{\mathcal{R}}_{+\ldots+}\cap\mathcal{E}$. In this region the canonical basis is real-valued as well. After simplifying, we find that the region is specified by the following constraints,

[TABLE]

3.3 Analytic integration at weight 2

The solutions of a canonical form system of differential equations may be written, order by order in $\epsilon$, in terms of Chen iterated integrals Chen . Because our family of integrals has multiple square roots that cannot be simultaneously rationalized through a variable change, it is not straightforward to rewrite these iterated integrals in terms of multiple polylogarithms, although such a representation is not necessarily excluded Heller:2019gkq . In the following section it is outlined how to perform the integration of the basis integrals at weight 2, in terms of a basis of manifestly real-valued logarithms and classical polylogarithms. The results are provided in the region $\mathcal{R}$ defined in the previous section, where the canonical basis integrals and the letters are real-valued.

It is useful to first consider the symbol Brown:2009qja ; Goncharov.A.B.:2009tja ; Goncharov:2010jf of the canonical basis integrals, and integrate the differential equations up to terms that lie in the kernel of the symbol, which will be referred to as ‘integrating the symbol’. Terms in the kernel of the symbols may be fixed afterwards by imposing that our solutions satisfy the system of differential equations, and by fixing overall transcendental constant from boundary conditions. Let us expand the basis integrals in $\epsilon$ as $\vec{B}=\sum_{k=0}^{\infty}B^{(k)}\epsilon^{k}$. The symbol of the $i$-th basis integral at order $k$ may then be obtained from the following recursive formula

[TABLE]

Note that the leading order of the canonical integrals $\vec{B}^{(0)}$ is constant and hence equal to the leading order of the boundary term given in Eq. (5).

Consider next the weight 2 integration of the symbol. At weight 2, i.e. order $\epsilon^{2}$, some canonical integrals are identically zero. Furthermore, the symbols of the remaining 40 nonzero integrals can be expressed in terms of the symbols of basis integrals

[TABLE]

By considering permutations of $p_{1},p_{2},p_{3}$ the number of independent symbols can be reduced even further, decreasing the number of integrals that need to be studied. We aim to integrate the symbol by writing a sufficiently general ansatz of logarithms and dilogarithms with undetermined prefactors, in the spirit of Ref. Duhr:2011zq . We then equate the symbol of the ansatz with the symbol of the individual master integrals and solve the resulting linear system. This can be done unambiguously since we express the symbol in terms of the linearly independent alphabet of Appendix A. We pick the basis of logarithms and dilogarithms in the ansatz such that they are manifestly real-valued in the region $\mathcal{R}$. Thus we require the arguments of the dilogarithms to lie in the range $(-\infty,1]$ for all of $\mathcal{R}$, while we require the arguments of the logarithms to be positive. This also guarantees that no branch-cuts will be crossed, such that the final expressions will be valid in at least the region $\mathcal{R}$. If one moves outside of $\mathcal{R}$ the basis functions may cross spurious branch cuts, which leads to incorrect results.

We denote the set of letters that appear at weight 2 by $\mathcal{A}_{2}$. In the region $\mathcal{R}$ the signs of these letters are completely fixed in the following way,

[TABLE]

For the logarithmic terms in the ansatz we therefore consider products of the type $\log(\pm l_{i})\log(\pm l_{j})$ with $l_{i},l_{j}\in\mathcal{A}_{2}$, where a minus may be included to fix a positive sign for the argument. Furthermore, we include dilogarithms with the following arguments in the ansatz,

[TABLE]

where we filter out dilogarithms whose argument does not lie between $(-\infty,1]$ in the region $\mathcal{R}$. We included the spurious letters $l_{33},l_{38}$ and $l_{41}$ in the ansatz, which do not appear in the symbol at weight 2, but are necessary for the ansatz to be sufficiently general. We identified these letters by using direct integration methods, outlined at the end of this section. Without knowledge of these spurious letters, we could have proceeded with an ansatz that includes all letters in $\mathcal{A}$.

After equating the symbol of the ansatz with the symbol of the canonical integrals, and solving the resulting system of equations, the following products of logarithms survive,

[TABLE]

Furthermore, the following dilogarithms are contained in the final result,

[TABLE]

It remains to fix terms that lie in the kernel of the symbol. At weight 2 these are transcendental constants and terms that have the form of a transcendental number times a logarithm, such as $i\pi\log(l_{i})$ or $\log(2)\log(l_{i})$. Since the canonical basis is real-valued in $\mathcal{R}$, it is already guaranteed there will be no contributions of the form $i\pi\log(\ldots)$. Furthermore, if terms of the type $\log(2)\log(l_{i})$ were missing in the solutions that we determined at the symbol level, they would not satisfy the system of differential equations. However, we find that our solutions already satisfy the differential equations without adding such terms, so that only additive transcendental constants are left undetermined. Note that to check if our weight 2 solutions satisfy the differential equations, we first need to derive the basis integrals at weight 1, which can be done in a similar manner as described here for the weight 2 case.

Next, we fix the overall transcendental constants from boundary conditions. It is not directly possible to use the heavy mass limit for this purpose, derived in Section 5, as it lies outside of region $\mathcal{R}$ where our solution is valid. Using the series solution strategy described in Section 4, we can transport the heavy mass limit up to high numerical precision to a regular point in $\mathcal{R}$, and use this to fix the remaining constants in our solution. We find in this way that some integrals carry additive constants proportional to $\pi^{2}$. The final result is provided in Appendix D.

Lastly, we note that the weight 2 integration of the symbol may also be performed by direct integration methods. Firstly, note that the most complicated Feynman integrals $I_{a_{1},a_{2},\ldots}$ appearing in Eq. (22) have 5 propagators but at weight 2 we only need to compute them at order $1/\epsilon$ because they carry a prefactor proportional to $\epsilon^{3}$ in the canonical basis. Setting up the Feynman parametrization for these integrals, and regularizing them using analytic regularization Panzer:2014gra , we find that they are linearly reducible Brown:2009ta ; Panzer:2015ida at order $1/\epsilon$ and can be computed algorithmically by direct integration, for example using HyperInt Panzer:2014caa . All other Feynman integrals appearing in the basis elements of Eq. (22) may also be computed using direct integration up to the required order in $\epsilon$. The resulting solutions are then given in terms of Goncharov multiple polylogarithms Goncharov:1998kja . One may use the Mathematica package of Ref. Frellesvig:2016ske , to convert these Goncharov polylogarithms to classical polylogarithms. The resulting expressions are generally larger than those obtained from a suitable ansatz, but they can be used to identify the spurious letters needed for the ansatz, such as the letters $l_{33},l_{38}$ and $l_{41}$ in Eq. (28).

3.4 One-fold integrals for weights 3 and 4

At weights 3 and 4, the polylogarithmic sectors have not been integrated out analytically. For those weights, we use the method of Ref. Caron-Huot:2014lda (see also Appendix E of Bonciani:2016qxi ) to represent the result in terms of one-fold integrals. The method works as follows. Firstly, it is clear that

[TABLE]

where $\gamma:[0,1]\rightarrow\mathbb{C}^{4}$ is some path in the phase space of $(s,t,m^{2},p_{4}^{2})$. Order by order in $\epsilon$ we may write

[TABLE]

Since the polylogarithmic sectors were integrated up to weight 2, we directly obtain the weight 3 expression as a one-fold integral over the weight 2 expression. By performing integration by parts, it is also possible to write the $i$-th order as a one-fold integral over the $(i-2)$-th order result,

[TABLE]

where $[\tilde{\mathbf{A}}]_{\gamma(0)}^{\gamma(1)}=\tilde{\mathbf{A}}(\gamma(1))-\tilde{\mathbf{A}}(\gamma(0))$. In this way one may obtain the weight 4 result as a one-fold integral, without needing an analytic expression of the basis integrals at weight 3.

For the base point of the integration we pick the point

[TABLE]

at which we may obtain high precision results for the canonical integrals using the series solution strategy of Section 4. We note that the diagonal entries of $\tilde{\mathbf{A}}$ at positions $38,40,42,45,$ and $49$ are divergent in the point $r^{\mathcal{R}}$. However, we may read off from the canonical basis that these integrals are identically zero in $r^{\mathcal{R}}$, so that the term $[\tilde{\mathbf{A}}]_{\gamma(0)}^{\gamma(1)}\vec{B}^{(i-1)}(\gamma(0))$ in Eq. (3.4) is well-defined.

Lastly, let us explicitly define a path $\gamma$ with basepoint $r^{\mathcal{R}}$, some endpoint in $\mathcal{R}$, and which lies fully inside region $\mathcal{R}$. A simple choice would be a straight line, but $\mathcal{R}$ is not convex. To ensure that $\gamma$ lies fully in $\mathcal{R}$, we consider a path that moves along a straight line in the $s,t$ and $m^{2}$ direction, but averages the $p_{4}^{2}$-coordinate between the upper and lower bounds in Eq. (3.2). We work this out in detail next.

Let us assume that $s\leq t$, so that the bounds on $p_{4}^{2}$ in region $\mathcal{R}$ are given by

[TABLE]

Defining a path for $s\geq t$ can be done in the same manner, if we change to the corresponding upper bound. Let us assume the endpoint of the path $\gamma$ is given by $(s^{\prime},t^{\prime},m^{\prime 2},p_{4}^{\prime 2})\in\mathcal{R}$. Next, write $p_{4}^{\prime 2}$ as a combination of the upper and lower bounds evaluated at $(s^{\prime},t^{\prime},m^{\prime 2})$,

[TABLE]

which yields

[TABLE]

Next, consider the following straight line path in the phase space of $(s,t,m^{2})$,

[TABLE]

We may then define the integration path by

[TABLE]

where $p_{\text{up}}^{2}(\gamma^{(3)}(\lambda))$ indicates the upper bound of $p_{4}^{2}$ evaluated at the point given by $\gamma^{(3)}(\lambda)$, and similarly for $p_{\text{down}}^{2}(\gamma^{(3)}(\lambda))$. By construction $\gamma$ has the endpoint $(s^{\prime},t^{\prime},m^{\prime 2},p_{4}^{\prime 2})$ at $\lambda=1$, and lies inside $\mathcal{R}$ for all $\lambda\in[0,1]$.

3.5 Series solutions

In Section 4.1 a strategy is outlined to obtain series solutions (truncated at high order) for Feynman integrals from their system of differential equations. Here we briefly point out that the method can be straightforwardly applied to compute the polylogarithmic sectors. Furthermore, the strategy can be used to obtain results in the physical region.

4 The elliptic sectors

There is a lot of recent progress in the theory of Feynman integrals evaluating to elliptic polylogarithms BrownLevin ; Broedel:2014vla ; Adams:2016xah ; Broedel:2017kkb (for applications of related classes of functions to high energy physics see e.g. Laporta:2004rb ; Kniehl:2005bc ; Adams:2013kgc ; Bloch:2013tra ; Bloch:2014qca ; Adams:2014vja ; Adams:2015gva ; Adams:2015ydq ; Remiddi:2016gno ; Primo:2016ebd ; Bonciani:2016qxi ; Adams:2016xah ; Passarino:2016zcd ; Harley:2017qut ; vonManteuffel:2017hms ; Ablinger:2017bjx ; Chen:2017pyi ; Hidding:2017jkk ; Bogner:2017vim ; Bourjaily:2017bsb ; Broedel:2017siw ; Laporta:2017okg ; Broedel:2018iwv ; Mistlberger:2018etf ; Lee:2018jsw ; Broedel:2018qkq ; Adams:2018bsn ; Adams:2018kez ; Broedel:2019hyg ; Bogner:2019lfa ; Kniehl:2019vwr ; Broedel:2019kmn ). However, most of these cases deal with relatively restricted kinematic dependence. More work would be needed to extend these methods to our situation with several independent square roots. Moreover, analytic continuation to physical regions is not yet very well understood. For this reason, we adopt here an approach based on series expansions methods.

4.1 Series solution of the differential equations

It is well known that when a closed form solution of a Feynman integral is not available, it is often possible to find a series expansion which provides a good approximation of the full solution. For an application of series expansion methods to single scale integrals see e.g. Pozzorini:2005ff ; Aglietti:2007as ; Mueller:2015lrx ; Lee:2017qql ; Lee:2018ojn ; Bonciani:2018uvv ; Mistlberger:2018etf . For integrals depending on several scales, multivariate series expansions have been used for special kinematic configurations (typically small or large energy limits, see e.g. Melnikov:2016qoc ; Melnikov:2017pgf ; Bonciani:2018omm ; Bruser:2018jnc ; Davies:2018ood ; Davies:2018qvx ). In this paper we reduce the multivariate problem to a single scale problem by defining a contour $\gamma(\lambda)$ that connects a boundary point to a fixed point of the phase space. We then find generalized power series solutions by solving the (single-scale) differential equations with respect to $\lambda$ (all the other variables are replaced by numbers).

More specifically, given a boundary point and a fixed target point, we define a contour $\gamma(\lambda)$ which connects the two points. It is then possible to define the differential equations with respect to $\lambda$. The problem is now one dimensional, and we can apply e.g. the methods of Wasow ; Lee:2018ojn , which provide a formula for the solution of the differential equations near a singular point, whose coefficients are fixed by solving recurrence relations. Here we proceed in a different way, by directly integrating the differential equations in terms of (generalized) power series. Note that in order to find series expansions for the elliptic sectors we also need to series expand the polylogarithmic sectors, which contain the majority of the integrals. For the polylogarithmic sectors we have, along a generic contour $\gamma(\lambda)$, differential equations of the form

[TABLE]

A series solution around a singular or regular point $\lambda_{0}$ is obtained in the following way. We first series expand the matrix elements of $\mathbf{A}_{\lambda}$ around $\lambda_{0}$, and then we recursively integrate the right hand side until the desired order of the dimensional regularization parameter is reached. Note that at each step one needs to perform only integrals of the form

[TABLE]

for $w\in\mathbb{Q}$ and $n\in\mathbb{Z}_{\geq 0}$, which can be solved again in terms of powers of $\lambda$ and of $\log(\lambda-\lambda_{0})$.

For the elliptic sectors the differential equations are not fully decoupled. However it is possible to take multiple derivatives of the differential equations and obtain higher order differential equations for single integrals. For an integral at order $\epsilon^{i}$ , say $g^{(i)}(\lambda)$, the differential equation has the form

[TABLE]

where by construction the inhomogeneous term is known at each iteration. In order to solve the equation we need to find a complete set of homogeneous (series) solutions. These can be found by using the Frobenius method Coddington . Once the homogeneous solutions $g_{h,j}(\lambda),\;j\in\{1,\dots,k\}$ are known, we obtain a series solution for the full equation by using the formula

[TABLE]

where $\chi^{(i)}_{j},\;j\in\{1,\dots,k\}$ is a set of boundary constants, $W({}\lambda{})$ is the Wronskian determinant, and $W_{j}({}\lambda{})$ is the determinant obtained by replacing the $j$-th column of the Wronskian matrix by $(0,\dots,0,1)$. Since the homogeneous solutions are a power series in the vicinity of the expansion point $\lambda_{0}$, order by order in the dimensional regularization parameter, we only need to consider integrals of the form (56), which are performed analytically in terms of powers of $\lambda$ and of $\log(\lambda-\lambda_{0})$ as explained above.

In principle the choice of the contour $\gamma(\lambda)$ is arbitrary, however we observed that choosing the contour to be a straight line renders the series expansion faster when compared to non-linear contours. In general the series expansion around a given point will converge only until the next singular point, therefore the choice of the expansion points requires some care. We start by considering real singular points, i.e. singular points corresponding to real values of $\lambda$. Given the set of real singular points along $\gamma(\lambda)$, we perform a series expansion around each singular point. Moreover, we define the domain of these expansions to be a segment centered at the expansion point, with radius equal to half the distance from the closest singular point (this ensures that the series converge sufficiently fast in their domain of definition). Some singularities might be complex, i.e. of the form $\lambda=\lambda_{r}+i\lambda_{i}$. For each complex singularity we expand around the three real points $\lambda_{r}-\lambda_{i},\lambda_{r},\lambda_{r}+\lambda_{i}$, and we define their domain as for the case of real singularities described before. At the end of this procedure, there might be segments along the contour not patched by any series. We then consider the uncovered regions and we expand around the middle point of each region, and we define the domain of these expansions to be, as before, a segment centered at the expansion point, with radius equal to half the distance from the closest singular point. We repeat this procedure until the entire contour is patched.

We remark that by singular point we mean any non-analytic point of the differential equations, i.e. points where the differential equations diverge, and branching points of possible square roots of the differential equations. We identify three classes of singularities. The first class of singularities are the so-called physical singularities, i.e. the singularities predicted by unitarity. When expanding around these singularities the solution will develop branch cuts across the expansion point, and the cut ambiguity is resolved by the Feynman prescription. Another class of singularities are the so-called non-physical singularities, i.e. singularities that are not physical, but which appear in the solution of the differential equations for a given basis of integrals (for example a canonical basis). These singularities correspond to singularities of the basis choice, but they cancel after rotating the basis back to a set of scalar Feynman integrals without prefactors. Specifically, non-physical singularities might originate from poles in the coefficients defining the integral basis, or from branching points of possible square roots appearing in the integral basis111In this work we consider only integral basis with algebraic coefficients.. Poles do not give rise to cuts, and no care is needed in their definition. On the other hand, for square roots giving rise to non-physical branch cuts, we can arbitrarily choose an analytic continuation across their branching points, since these cuts cancel at the level of physical quantities. Finally, there are the so-called spurious singularities. These singularities are not predicted by unitarity, and they do not appear in the solution of the differential equations on the physical sheet. Therefore we choose our contours to be entirely contained in the physical sheet (see also moriello_expansions for a detailed discussion about the singularities of differential equations in the context of our expansion strategy).

Finally, we remark that the procedure described above is entirely algorithmic, and can be used, given a boundary point, to transport the solution to any point of interest.

4.2 Differential equations for the elliptic sectors

Our basis for the elliptic sectors is provided in Appendix B. The corresponding differential equations have the form

[TABLE]

where the vector $\vec{G}_{66-73}(\vec{x},\epsilon)$ depends on the polylogarithmic integrals $\vec{B}_{1-65}(\vec{x},\epsilon)$, and the homogeneous matrix has the schematic form

[TABLE]

where the lines separate sectors $I_{0,1,1,1,1,1,1,0,0}$ and $I_{1,1,1,1,1,1,1,0,0}$. We see that the integrals $\vec{B}_{66-67}$ and $\vec{B}_{70-71}$ are coupled, and they satisfy second order differential equations. By defining a contour $\gamma(\lambda)$ we obtain a single set of differential equations with respect to $\lambda$,

[TABLE]

where the matrix $\mathbf{A}_{\lambda}(\lambda,\epsilon)$ is given by,

[TABLE]

Given the structure of the homogeneous matrix it is clear how to perform a series expansion along a contour. We start with the expansion around a given point of the contour. For a given $\epsilon$ order we series expand the polylogarithmic sectors as described above. We then use Eq. (58) to expand the coupled integrals. We iterate until the required $\epsilon$ order is solved ($\epsilon^{4}$ in our case). We repeat the procedure for a suitably chosen set of points until the entire contour is covered. In the next section we discuss in detail the expansion along two contours of interest.

4.3 Numerical results

In this section we provide a numerical analysis of our solution. We use $p_{4}^{2}=13/25\approx(125/172)^{2}$, which is a good approximation to the physical value of the Higgs mass, assuming a unit mass for the top quark. However, the suitability and accuracy of the method is not related in any way to the numerical values of the kinematical points chosen. We could just as well have chosen values of $p_{4}^{2}$, which are related to the ratio of the Higgs mass to the charm or $b$-quark mass, or to anything else.

We solve the basis integrals along a contour crossing a particle production threshold, in particular we chose a path from the kinematic point $\textrm{P}_{\textrm{below}}=(s=2,t=-1,p_{4}^{2}=13/25)$ to the kinematic point $\textrm{P}_{\textrm{above}}=(s=6,t=-1,p_{4}^{2}=13/25)$. Our starting point is the heavy mass limit $\textbf{0}=(0,0,0)$ where we can use the results of Section 5. We perform an expansion along a straight path $\gamma_{1}$ from 0 to the point $\textrm{P}_{\textrm{regular}}=(2.5,-1,13/25)$, which does not cross any singularity. However, by inspection of the differential equation matrices we see that they possess a singularity in $t=-\frac{52}{495}$, which spoils the convergence of the expansion in 0. To obtain a better precision we then patch ten expansions on the line from 0 to $\textrm{P}_{\textrm{regular}}$, by making sure that we never evaluate series at more than half the distance from the closest singularity (either physical or spurious).

We now have high precision values at $\textrm{P}_{\textrm{regular}}$ that allows us to start expanding from here. Therefore we write an expansion along a path $\gamma_{2}$ joining $\textrm{P}_{\textrm{regular}}$ and $\textrm{P}_{\textrm{singular}}=(4,-1,13/25)$ which is a singular point since it lies on the threshold $s=4m^{2}$. Thus, we use this expansion only up to the midpoint $\textrm{P}_{\textrm{mid}}=(3,-1,13/25)$ (Expansion 1). To go beyond this point, we perform an expansion (Expansion 2) from $\textrm{P}_{\textrm{singular}}$ to $\textrm{P}_{\textrm{mid}}$ using the values obtained there to fix the constants of integration. Finally, we perform an analytic continuation of the second expansion to obtain an expression valid above threshold up to $\textrm{P}_{\textrm{above}}$. The analytic continuation is straightforward since, when considering a series along a contour such that only one of the external scales varies ($s$ in our case), only terms of the form $\log{(s-4m^{2})}$ and $\sqrt{s-4m^{2}}$ have to be analytically continued. The situation is depicted in Figure 3. The expansion is carried up to order $\epsilon^{4}$ and the series are truncated at order $t^{50}$. In Figure 4 the integrals $B^{(4)}_{72-73}$ are plotted together with some values obtained with standard numerical software. The remaining plots can be found in Appendix E. To measure the numerical error of our method we can compare values obtained patching different contours. For example, we compared the values at $\textrm{P}_{\textrm{mid}}$ obtained expanding from 0 to $\textrm{P}_{\mathrm{regular}}$ and then from $\textrm{P}_{\textrm{regular}}$ to $\textrm{P}_{\textrm{mid}}$ with those obtained expanding directly from 0 to $\textrm{P}_{\textrm{mid}}$. The worst error is for $B_{72}^{(4)}$ where we get difference of $\mathcal{O}(10^{-16})$. For comparison, the numerical evaluation with FIESTA Smirnov:2015mct of all the 73 master integrals in one point above the $s=4m^{2}$ threshold requires $\mathcal{O}(10^{4})$ seconds on 48 cores. On the other hand, finding a series expansion above threshold requires $\mathcal{O}(100s)$ on an ordinary laptop on 1 core, which is $\mathcal{O}(1\text{ second})$ per integral, after which the plots are virtually instantly generated.

5 Boundary terms

In this section we describe how to obtain boundary conditions for the system of differential equations associated with our basis of integrals. Similar to the computation of the planar Higgs + jet families in Ref. Bonciani:2016qxi , we consider the heavy mass limit, where $s$, $t$ and $p_{4}^{2}$ tend to zero. For the planar families graph theoretical prescriptions Smirnov:1990rz ; Smirnov:1994tg ; Smirnov:2002pj were used to compute the limit. Although the same strategy may in principle be applied for the current non-planar family, we perform the computation of the boundary terms for the current family using the method of expansions by regions.

The method of expansions by regions was originally developed in Ref. Beneke:1997zp for the case of near threshold expansions, and has since been further clarified (for various limits) in various works, see for example Ref. Jantzen:2011nz , and also Ref. Semenova:2018cwy for a recent overview. In Ref. Smirnov:1999bza the strategy was considered in the Feynman parametric representation, and a systematic geometry based algorithm for finding the relevant regions was presented in Refs. Pak:2010pt ; Jantzen:2012mw . We will use the formulation in the Feynman parametrization, and the package asy2.1.m from Ref. Jantzen:2012mw to find the regions.

First we present the results. We parametrize the heavy mass limit by a straight line,

[TABLE]

with $x\rightarrow 0$. It turns out that only the first two integrals of the polylogarithmic sectors are nonzero in the limit. More specifically, we have

[TABLE]

In the elliptic sectors all but the last integral are identically zero in the heavy mass limit. In particular, we have

[TABLE]

Next, we work out the computation of the last boundary term $B_{73}$ as an illustrative example. This basis integral is given explicitly by

[TABLE]

All the Feynman integrals that appear in $B_{73}$ lie in the same sector, but with different numerators. Using asy2.1.m we find that the contributing regions as $x\rightarrow 0$ are

[TABLE]

where the $\alpha_{i}$’s denote the Feynman parameters associated to each propagator in the Feynman parametrization. The asymptotic limit of the integrals then takes the form,

[TABLE]

where the superscripts denote the contributions of both regions respectively. In other words, for the integrals with superscripts the kinematics has been rescaled in the manner of Eq. (63), the Feynman parameters are rescaled according to Eq. (5), and the overall dependence on $x$ has been explicitly factored out of the integrand.

Rescaling the prefactors in Eq. (5), and plugging in the results of Eq. (5) yields the asymptotic limit of $B_{73}$,

[TABLE]

There are 3 surviving contributions, which lie in region 2. It is sufficient to compute these at leading order in $x=0$, as the higher orders vanish in the limit. The integrations themselves can be done by picking a suitable integration order and integrating one Feynman parameter at a time, while keeping full dependence on the dimensional regulator $\epsilon$.222We found that some of the integrations converge in disjoint domains of $\epsilon$. In situations without a domain of $\epsilon$, where all parts of an integral are convergent, one usually considers each part in its own domain of convergence from which it can analytically be continued to a desired common final domain after an evaluation. Still to be on the safe side, we prefer to proceed with an auxiliary analytical regularization, and assign the third propagator the exponent $\nu_{3}$, which serves as an additional regulator. After performing all integrations we take the limit $\nu_{3}\rightarrow 1$. The end result is found to be the same. The results are given by

[TABLE]

Plugging these expressions into Eq. (69) yields the boundary term of $B_{73}$, stated in Eq. (5). We used the assumption $t<0$ during the integration to avoid branch cuts, and we may analytically continue the expression to the physical region by using the Feynman prescription, which tells us to interpret $t$ as having an infinitesimally small positive imaginary part.

The remaining boundary terms may be computed using the same analysis as was done here for the most complicated sector. Although we found that all basis integrals except for $B_{1}$, $B_{2}$, and $B_{73}$ are zero, showing this requires the computation of numerous integrals which cancel with each other at the end.

6 Conclusion

The Feynman diagrams of the $gg\rightarrow gH$ amplitudes with full heavy quark mass dependence fit into seven integral families. The computation of the planar families was presented in Ref. Bonciani:2016qxi , leaving two non-planar integral families to be evaluated. In this paper we presented the computation of one of the two families of non-planar master integrals. The family under consideration consists of 73 master integrals, of which 65 are polylogarithmic, while the remaining 8 integrals in the top sectors are elliptic, implying that their solution may not be presented in terms of multiple polylogarithms. Therefore the paper was split up into two parts, the first dealing with the computation of the polylogarithmic integrals, and the second with the computation of the elliptic sectors.

The polylogarithmic sectors may be computed by choosing a basis that puts their associated system of differential equations in canonical form, in which all dependence on the dimensional regularization parameter is in an overall prefactor, and the partial derivative matrices are combined into a matrix that is in $d\log$-form. Although the formal solution of such a system may be immediately given in terms of Chen’s iterated integrals, the presence of many square roots, that cannot be simultaneously rationalized, makes it unclear how to convert the iterated integrals to multiple polylogarithms which admit fast numerical evaluation. We showed that at weight 2 the canonical basis integrals may nonetheless be solved by starting from an ansatz of dilogarithms and logarithms. The weight 3 and 4 expressions could then be written as 1-fold integrals over polylogarithmic expressions.

Alternatively, we showed how the system of differential equations may be solved in terms of (generalized) series solutions. The approach relies on reducing the multidimensional system of differential equations to a one-dimensional system, by integrating along a line towards a given kinematic point of interest. The differential equations are then expanded in terms of a (generalized) series, and the integration is performed analytically up to very high order. This expansion based approach furthermore works to obtain high precision numerical results for both the polylogarithmic and the elliptic sectors, and in a sense trivializes the analytic continuation of the family of integrals to the physical region as well. At any stage it is only required to perform the analytic continuation of powers of logarithms across their branch cuts, and the direction in which to cross the branch cut is determined by the Feynman prescription.

To compute the Higgs + jet cross section it is necessary to perform the evaluation of another non-planar integral family, which we believe may be done using the same methods as have been presented for the current family. We leave the computation of the remaining family to a future publication.

Acknowledgements

HF, MH, LM and GS would like to thank ETH Zürich for the hospitality during the preparation of this work. The work of HF, MH, LM and GS was supported by STSM Grants from the COST Action CA16201 PARTICLEFACE. The work of GS was supported by ETH Zürich and by the Pauli Center for Theoretical Studies, Zürich. The work of FM was supported by the Swiss National Science Foundation project No. 177632 (ElliptHiggs). The work of HF is part of the HiProLoop project funded by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement 747178. The work of VS was supported by RFBR, grant 17-02-00175. This research received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme, grant agreements No 725110 (Amplitudes) and 647356 (CutLoops).

Appendix A Alphabet of the polylogarithmic sectors

The full alphabet is given by the following 69 letters,

[TABLE]

[TABLE]

[TABLE]

The $q_{i}$ are given by the following polynomials,

[TABLE]

The terms labelled by $\{r_{i}\}$ are the square roots. The first 15 roots appear in the polylogarithmic sectors, while $r_{16}$ appears in the elliptic sectors (see next section),

[TABLE]

The labelling has been chosen so that the radicands of the roots are irreducible polynomials. In the polylogarithmic basis elements $B_{1},\ldots,B_{65}$, the 15 roots only appear in 10 combinations,

[TABLE]

It may also be verified that these same 10 combinations are sufficient to express all ratios of roots appearing in the letters. Hence, in principle is it possible to combine these and work with just 10 independent square roots for the polylogarithmic sectors.

In the choice of basis for the elliptic sectors, roots $r_{2}$ and $r_{16}$ appear separately. Therefore, there are 12 independent combinations of roots in the full basis, given in Appendix B.

Appendix B Canonical basis and basis for elliptic sectors

The set of 73 master integrals is represented in Figure 5. The canonical basis of the polylogarithmic sectors is given by the following 65 integrals,

[TABLE]

In addition we consider the following choice of basis for the elliptic sectors,

[TABLE]

Appendix C Maximal cut of the elliptic sectors

Of the discussed integral families, only two are genuinely elliptic as we have seen. These are the six-propagator sector 126 consisting of integrals $I_{66}$ and $I_{67}$, and the seven-propagator sector 127 consisting of the six integrals $I_{68}$-$I_{73}$. Representing these sectors by their most basic members $I_{011111100}$ and $I_{111111100}$ we may derive their maximal cut using the loop-by-loop Baikov representation Frellesvig:2017aai , in both cases yielding a one-fold integral representation. In $d=4$ these cut integrals become

[TABLE]

of which the first has been analyzed in refs. Primo:2016ebd ; Hidding:2017jkk , and the second discussed in ref. Frellesvig:2019kgj . These two integrals both evaluate to elliptic integrals, and can be expressed as integrals over an elliptic curve, proving that genuinely elliptic structures have to be present in representations of the integrals in these sectors.

We see that this elliptic curve,

[TABLE]

is in common between the two sectors, implying that it might be possible to express all integrals in the integral family considered in this paper, in terms of polylogarithmic functions defined on that elliptic curve in the spirit of refs. BrownLevin ; Broedel:2017kkb .

Appendix D Polylogarithmic sectors up to weight 2

We provide below the results for the canonical basis integrals up to weight 2 in terms of logarithms and dilogarithms, provided in the region $\mathcal{R}$ defined in Eq. (3.2), where the result is manifestly real-valued. Let $\vec{B}=\sum_{k=0}^{\infty}\vec{B}^{(k)}\epsilon^{k}$. At weight 0 we have

[TABLE]

At weight 1 we have

[TABLE]

and $B_{i}^{(1)}=0$ for all other $i\leq 65$. Finally, at weight 2 we have

[TABLE]

Appendix E Plots of basis integrals

In this appendix we provide plots of all the integrals of our basis along the line going from the kinematic point $\textrm{P}_{\textrm{below}}=(s=2,t=-1,p_{4}^{2}=13/25)$ to $\textrm{P}_{\textrm{above}}=(6,-1,13/25)$, which are obtained using the series solution strategy. The blue and orange lines are the real and imaginary parts and the dots are values obtained with FIESTA for comparison. Note that the only noticeable mismatch is for $B_{68}^{(4)}$ and $B_{69}^{(4)}$, but it is well inside the uncertainty produced by FIESTA (not reported here).