Generalised power series expansions for the elliptic planar families of Higgs + jet production at two loops

TL;DR

This paper develops a method using generalized power series expansions to efficiently compute two-loop master integrals for Higgs + jet production in QCD, accounting for heavy-quark mass effects across all kinematic regions.

Contribution

The authors introduce a novel differential equation approach with power series solutions for two-loop integrals, enabling high-precision calculations across physical thresholds.

Findings

High-precision integral evaluations with errors around 10^{-32}

Rapid computation time of about 1 second per integral

Method applicable to a broad class of Feynman integrals

Abstract

We obtain generalised power series expansions for a family of planar two-loop master integrals relevant for the QCD corrections to Higgs + jet production, with physical heavy-quark mass dependence. This is achieved by defining differential equations along contours connecting two fixed points, and by solving them in terms of one-dimensional generalised power series. The procedure is efficient and can be repeated in order to reach any point of the kinematic regions. The analytic continuation of the series is straightforward and we present new results below and above the physical thresholds. The method we use allows to compute the integrals in all kinematic regions with high precision. Performing a series expansion on a typical contour above the physical threshold takes on average $\mathcal{O}(1 \text{ second})$ per integral with worst relative error of $\mathcal{O}(10^{-32})$, on a single CPU core. After the series is found the numerical evaluation of the integrals in any point of the contour is virtually instant. Our approach is general and can be applied to Feynman integrals provided that a set of differential equations is available.

Tables

Table 1. Table 1: Timings for the computation of all the 73 master integrals in the point above threshold s 12 = 5 , s 13 = − 1 , p 4 2 = 13 15 , m 2 = 1 formulae-sequence subscript 𝑠 12 5 formulae-sequence subscript 𝑠 13 1 formulae-sequence superscript subscript 𝑝 4 2 13 15 superscript 𝑚 2 1 s_{12}=5,s_{13}=-1,p_{4}^{2}=\frac{13}{15},m^{2}=1 up to and including order ϵ 4 superscript italic-ϵ 4 \epsilon^{4} , on 1 CPU core. The point is reached by expanding along the contour γ thr subscript 𝛾 thr \gamma_{\text{thr}} defined in Section 4.1 . We use FIESTA 4.1 with the Vegas Hahn:2004fe integrator while pySecDec 1.4.3 with the quasi-Monte Carlo integrator. The ’1 segment’ lines correspond to the average time needed to perform the expansion along one of the four line segments of the contour γ thr subscript 𝛾 thr \gamma_{\text{thr}} . This timing is more representative of the typical time needed to reach a given point when several points have been already computed. Indeed, once several points have been computed (for example when performing a phase space integration) one can use as a boundary point the closest point to the next target point, which will likely be ’close’ to an already computed point.

	Truncation	Relative error	Time (73 MIs)	$\frac{\text{Time}}{\text{integral}}$
Expansion (${\gamma }_{\text{thr}}$)	$𝒪\left(t^{85}\right)$	$\le {10}^{-24}$	79 sec	1.11 sec
Expansion (${\gamma }_{\text{thr}}$)	$𝒪\left(t^{125}\right)$	$\le {10}^{-32}$	162 sec	2.21 sec
Expansion (${\gamma }_{\text{thr}}$, 1 segment)	$𝒪\left(t^{85}\right)$	$\le {10}^{-24}$	20 sec	0.27 sec
Expansion (${\gamma }_{\text{thr}}$, 1 segment)	$𝒪\left(t^{125}\right)$	$\le {10}^{-32}$	40 sec	0.55 sec
FIESTA 4.1 (Vegas)		$\le {10}^{-2}$	$60000$ sec	$821$ sec
pySecDec 1.4.3 (QMC)		$\le {10}^{-2}$	$25000$ sec	$342$ sec

Table 2. Table 2: Numerical results for the elliptic sectors in the point s 12 = 5 , s 13 = − 1 , p 4 2 = 13 15 , m 2 = 1 formulae-sequence subscript 𝑠 12 5 formulae-sequence subscript 𝑠 13 1 formulae-sequence superscript subscript 𝑝 4 2 13 15 superscript 𝑚 2 1 s_{12}=5,s_{13}=-1,p_{4}^{2}=\frac{13}{15},m^{2}=1 at order ϵ 4 superscript italic-ϵ 4 \epsilon^{4} .

#	Real part	Imaginary part	Rel. err.
66	$-0.3229462567669706555224248824$	$+0.0117911783550457146213317093$	$={10}^{-32}$
67	$-0.3039273344500928782645744918$	$-0.0091549103610263866197445332$	$={10}^{-32}$
68	$+0.4989288717476986372753776549$	$+0.3357385827462031355699714492$	$={10}^{-32}$
69	$-0.1216791054426520965987437711$	$-0.0043091895703585235915244984$	$={10}^{-32}$
70	$+2.0936824603121477765883288075$	$+0.0934305088779206383805166362$	$={10}^{-32}$
71	$+0.4453410820313442706096983608$	$-0.0041620964011457311182442998$	$={10}^{-32}$
72	$-0.2737950895536402954579381308$	$+0.0027361587278654015689589630$	$={10}^{-32}$
73	$-0.6420218523272138407244813424$	$+0.0011438044096963059034512574$	$={10}^{-32}$

Equations

$$\frac{\partial }{\partial x_i}\vec{f}(\vec{x},ϵ)={\mathbf{A}}_{x_i}(\vec{x},ϵ)\vec{f}(\vec{x},ϵ),$$

$$\frac{\partial {\mathbf{A}}_{x_j}}{\partial x_i}-\frac{\partial {\mathbf{A}}_{x_i}}{\partial x_j}+[{\mathbf{A}}_{x_i},{\mathbf{A}}_{x_j}]=0,$$

$$d\vec{f}(\vec{x},ϵ)=ϵd\tilde{\mathbf{A}}(\vec{x})\vec{f}(\vec{x},ϵ).$$

$$\vec{f}(\vec{x},ϵ)=\sum _{k=0}^{\infty }{ϵ}^k{\vec{f}}^{(k)}(\vec{x}),$$

$$\vec{f}(\vec{x},ϵ)=\vec{f}({\vec{x}}_0,ϵ)+\sum _{k\ge 1}{ϵ}^k\sum _{j=1}^k{\int }_0^1{\gamma }^{\ast }(d\tilde{\mathbf{A}})(t_1){\int }_0^{t_1}{\gamma }^{\ast }(d\tilde{\mathbf{A}})(t_2)\dots {\int }_0^{t_{j-1}}{\gamma }^{\ast }(d\tilde{\mathbf{A}})(t_j){\vec{f}}^{(k-j)}({\vec{x}}_0),$$

$$G(a_1,a_1\dots a_w,t)={\int }_0^{t}dt\frac{1}{t-a_1}G(a_2,\dots ,a_w,t),G({\vec{0}}_w,t)\equiv \frac{\log (t{)}^w}{w!},$$

$$\tilde{\mathbf{A}}(\vec{x})=\sum _{i=1}^{d_{\alpha }}{\mathbf{C}}_i\log ({\alpha }_i(\vec{x})),$$

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

$$\gamma (t):t\mapsto \{x_1(t),\dots ,x_m(t)\},\vec{x}(t_a)=\vec{a},\vec{x}(t_b)=\vec{b},$$

$$\frac{\partial }{\partial t}\vec{f}(t,ϵ)={\mathbf{A}}_t(t,ϵ)\vec{f}(t,ϵ),$$

$${\mathbf{A}}_t(t,ϵ)=\sum _{i=1}^m{\mathbf{A}}_{x_i}(t,ϵ)\frac{\partial x_i(t)}{\partial t}.$$

$${\vec{f}}_{\text{sing}}^{(i)}(t)=\sum _{j_1\in S_i}\sum _{j_2=0}^{\infty }\sum _{j_3=0}^{N_i}\vec{c}{}^{(i,j_1,j_2,j_3)}(t-\tau {)}^{w_{j_1}+j_2}\log {(t-\tau )}^{j_3},$$

$${\vec{f}}_{\text{reg}}^{(i)}(t)=\sum _{j=0}^{\infty }\vec{c}{}^{(i,j)}(t-\tau {)}^j.$$

$$(t-\tau {)}^w\log (t-\tau {)}^k\rho (t),w\in \mathbb{Q},k\in \mathbb{N},$$

$$\frac{\partial }{\partial t}\vec{f}(t,ϵ)=ϵ{\mathbf{A}}_t(t)\vec{f}(t,ϵ),$$

$$\vec{f}(t,ϵ)=\vec{f}(t_0,ϵ)+\sum _{k=1}^{\infty }{ϵ}^k\sum _{j=1}^k{\int }_{t_0}^{t}d{t}_1{\mathbf{A}}_t(t_1){\int }_{t_0}^{t_1}d{t}_2{\mathbf{A}}_t(t_2)\dots {\int }_{t_0}^{t_{j-1}}d{t}_j{\mathbf{A}}_t(t_j){\vec{f}}^{(k-j)}(t_0).$$

$${\mathbf{A}}_t(t)=\sum _{i=0}^{\infty }{\mathbf{A}}_t^{(i)}(t-\tau {)}^{w_i},w_i\in \mathbb{Q},$$

$$\vec{f}(t,ϵ)=\vec{f}(t_0,ϵ)$$

$$+\sum _{k=1}^{\infty }{ϵ}^k\sum _{j=1}^k\sum _{i_1,\dots ,i_j=0}^{\infty }{\mathbf{A}}_t^{(i_1)}\dots {\mathbf{A}}_t^{(i_j)}{\int }_{t_0}^{t}d{t}_1(t_1-\tau {)}^{w_{i_1}}\dots {\int }_{t_0}^{t_{j-1}}d{t}_j(t_j-\tau {)}^{w_{i_j}}{\vec{f}}^{(k-j)}(t_0).$$

$$\int (t-\tau {)}^w\log (t-\tau {)}^{k}dt,w\in \mathbb{Q},k\in \mathbb{N}.$$

$$\frac{{\partial }^k{g}^{(i)}(t)}{\partial t^k}+a_1(t)\frac{{\partial }^{k-1}{g}^{(i)}(t)}{\partial t^{k-1}}+\cdots +a_k(t)g^{(i)}(t)={\beta }^{(i)}(t),$$

$$\frac{{\partial }^{k}h(t)}{\partial t^k}+a_1(t)\frac{{\partial }^{k-1}h(t)}{\partial t^{k-1}}+\cdots +a_k(t)h(t)=0,$$

$$\hat{\mathbf{H}}(t)=\left(\begin{array}{ccccc}{h}_1(t)& & \cdots & & h_k(t)\\ \frac{\partial }{\partial t}{h}_1(t)& & \cdots & & \frac{\partial }{\partial t}{h}_k(t)\\ ⋮& & \ddots & & ⋮\\ \frac{{\partial }^{k-1}}{\partial t^k}{h}_1(t)& & \cdots & & \frac{{\partial }^{k-1}}{\partial t^{k-1}}{h}_k(t)\end{array}\right),$$

$$g^{(i)}(t)=\sum _{j=1}^k\left(h_j(t){\chi }_j^{(i)}+h_j(t){\int }_{t_0}^t\frac{W_j(\hat{\mathbf{H}}(s))}{W(\hat{\mathbf{H}}(s))}{\beta }^{(i)}(s)ds\right),$$

$$\frac{{\partial }^k{g}^{(i)}(t)}{\partial t^k}+\frac{b_1(t)}{t-\tau }\frac{{\partial }^{k-1}{g}^{(i)}(t)}{\partial t^{k-1}}+\cdots +\frac{b_k(t)}{(t-\tau {)}^k}{g}^{(i)}(t)={\beta }^{(i)}(t),$$

$$h_i(t)=\sum _{j_1\in S_i}\sum _{j_2=0}^{\infty }\sum _{j_3=0}^{K_i}{c}_i^{(j_1,j_2,j_3)}(t-\tau {)}^{{\lambda }_{j_1}+j_2}\log {(t-\tau )}^{j_3},$$

$$\lambda (\lambda -1)\cdots (\lambda -k+1)+b_1(\tau )\lambda (\lambda -1)\cdots (\lambda -k+2)+b_k(\tau )=0.$$

$$t^2\frac{{\partial }^{2}h(t)}{\partial t^k}+t{b}_1(t)\frac{\partial h(t)}{\partial t}+b_2(t)h(t)=0.$$

$$h_1(t)=t^{{\lambda }_1}\sum _{i=0}^{\infty }{c}_{1,i}{t}^i,h_2(t)=t^{{\lambda }_2}\sum _{i=0}^{\infty }{c}_{2,i}{t}^i,(c_{1,0},c_{2,0}=1)$$

$$h_2(t)=h_1(t)\int \frac{1}{h_1^2(t)}\exp \left(-{\int }^t\frac{b_1(t^{\prime })}{t^{\prime }}d{t}^{\prime }\right)dt,$$

$${\vec{f}}_{[\gamma ]}(t,ϵ)=\sum _{i\in S_e}{\vec{f}}_{[\gamma ]}(t,ϵ{)}_{[t_i-r_i,t_i,t_i+r_i^{\prime }]}+\sum _{i\in S_r}{\vec{f}}_{[\gamma ]}(t,ϵ{)}_{[{\kappa }_i-{\rho }_i,{\kappa }_i,{\kappa }_i+{\rho }_i^{\prime }]},$$

Full text

1 Introduction

The computation of Feynman integrals is a central ingredient for the prediction of collider experiments. In the past decades, we have seen an enormous progress in our capabilities to efficiently compute Feynman integrals in closed analytic form or in a purely numerical way. From the analytic side, several techniques are available. Some of the most effective techniques are the differential equations method  Kotikov:1990kg ; Kotikov:1991pm ; Bern:1993kr ; Remiddi:1997ny ; Gehrmann:1999as and the direct integration methods Brown:2009ta ; Panzer:2015ida . In dimensional regularisation, one is able to reduce a given (generally large) set of scalar Feynman integrals to a minimal set of linearly independent integrals, called master integrals (MIs), by using integration-by-parts identities (IBP) Tkachov:1981wb ; Chetyrkin:1981qh ; Laporta:1996mq ; Laporta:2001dd . Once a basis is identified, it is possible to define a closed system of first order linear differential equations, that can be solved in terms of iterated integrals Chen-iterated . Our understanding of differential equations for Feynman integrals has been further refined by the identification of canonical bases of integrals Henn:2013pwa , which make the solution of the equations in terms of iterated integrals completely algorithmic. In some cases Feynman integrals can be computed in terms of special functions known as multiple polylogarithms (MPLs) Goncharov:1998kja or, more recently, in terms of their elliptic generalisation, elliptic multiple polylogarithms (eMPLs) BrownLevin ; Broedel:2014vla ; Adams:2015ydq ; Broedel:2017kkb (for analytic results involving functions of elliptic type 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 ). Even though having a representation in terms of known functions is important from the conceptual and practical side, in recent years this approach has become challenging. For state-of-the-art computations one usually encounters multi-loop integrals depending on several mass scales. In this case the differential equations exhibit complicated analytic structures, and their solution in terms of known special functions is not well understood yet. Similar challenges are encountered when solving Feynman integrals by direct integration, e.g., in Feynman parameter space. From the purely numerical side, several methods are available to compute Feynman integrals by using Monte Carlo integration techniques Borowka:2015mxa ; Smirnov:2015mct . As opposed to the analytic approach, these methods are fully algorithmic. Nonetheless multi-loop multi-scale integrals generally present numerical instabilities that make their numerical integration challenging. A different numerical approach has been used in Boughezal:2007ny ; Czakon:2007qi ; Czakon:2008zk ; Mandal:2018cdj , where the solution of differential equations for Feynman integrals is obtained by using Runge-Kutta algorithms.

A third route of exploration has been methods based on series expansions. When it is difficult to obtain a closed form solution for a given integral, it is usually possible to obtain a (generalised) power series expansion of the solution. Series representations have a number of useful features. It is usually possible to compute several orders of the expansion and obtain an arbitrarily good approximation of the full solution. Moreover their numerical evaluation is virtually instant, since each term of the expansion is an elementary or a relatively simple function. All these features make series expansions a natural candidate to solve large classes of complicated Feynman integrals. Series expansion methods have been mostly applied to single scale problems in e.g. Pozzorini:2005ff ; Aglietti:2007as ; Mueller:2015lrx ; Lee:2017qql ; Lee:2018ojn ; Bonciani:2018uvv ; Mistlberger:2018etf . On the other hand, for integrals depending on several scales, series expansions have been performed with respect to one variable, parametrising special kinematic configurations, while keeping the dependence on the remaining variables exact (see e.g. Caffo:2002ch ; Caffo:2008aw ; Melnikov:2016qoc ; Melnikov:2017pgf ; Bonciani:2018omm ; Bruser:2018jnc ; Davies:2018ood ; Davies:2018qvx ), or to transport analytic boundary conditions in various regions Heller:2019gkq . Therefore, in the multivariate case, it is desirable to study a systematic approach to obtain results in all points of the kinematic regions.

In this paper we reduce the computation of a set of multivariate Feynman integrals Bonciani:2016qxi to a single scale problem, by defining differential equations along contours connecting two generic points of the kinematic regions. We then find generalised power series solutions by solving the (single-scale) differential equations with respect to the contour parameter (while replacing all the other variables with numbers). In this way the solution can be transported from a base point, where the integrals are assumed to be known, to a generic target point. We show that this approach is efficient, and can be repeated to compute the integrals in any point of the kinematic regions, with high numerical precision. More specifically, we apply this method to a family of planar (elliptic) Feynman integrals relevant for the two-loop QCD corrections to Higgs + jet production, below and above the heavy-quark threshold. Previously Bonciani:2016qxi these integrals were computed in the Euclidean region by using integral representations. Our results are new, and provide at the same time the analytic continuation of these integrals to the physical region, and an efficient method for their numerical evaluation. Further applications of these methods to the non-planar integral topologies are presented in a companion paper non-planars .

The paper is organised as follows. In Section 2 we review general properties of the differential equations for dimensionally regulated scalar Feynman integrals, and their solution in terms of iterated integrals. In Section 3 we describe the series expansion strategy used in this paper. We show that after defining the (multi-scale) differential equations along a (one-dimensional) contour, the series solution can be obtained by series expanding the differential equations and iteratively integrating them up to the desired order of the dimensional regulator. In Section 4 we apply this strategy to a family of planar integrals relevant for Higgs + jet production in the physical region. We show how, once a series expansion is found, the analytic continuation is performed in a straightforward manner. We finally show high precision numerical results and timings, with comparisons to sector decomposition programs. In Section 5 we draw our conclusions.

2 Differential equations for dimensionally regulated Feynman integrals

By using standard IBP reduction techniques, it is possible to identify a basis $\vec{f}(\vec{x},\epsilon)$ for a set of dimensionally regulated scalar Feynman integrals, where $\vec{f}(\vec{x},\epsilon)=\left\{f_{1}(\vec{x},\epsilon),\dots,f_{n}(\vec{x},\epsilon)\right\}$ and $\vec{x}=\left\{x_{1},\dots,x_{m}\right\}$ is the set of kinematic invariants. Given a basis, one is also able to define a system of first order linear differential equations satisfied by the basis, that in full generality takes the form

[TABLE]

where $\mathbf{A}_{x_{i}}$ is an $n$-by-$n$ matrix that depends rationally on its variables. If the set of basis integrals is minimal the differential equations satisfy the integrability condition,

[TABLE]

where the last term is a commutator. Nonetheless the applicability of our method does not rely on this condition, and it can be applied also to an overcomplete set of integrals.

The basis $\vec{f}$ is not unique. In Henn:2013pwa it was conjectured that, with a basis change, it is possible to cast differential equations for Feynman integrals in a canonical form, where the dependence on the dimensional regulator is factorised. In differential form the canonical equations have the following form

[TABLE]

In dimensional regularisation we are interested in a solution around $\epsilon=0$. By series expanding $\vec{f}(\vec{x},\epsilon)$,

[TABLE]

the solution of Eq. (3) can be written in terms of iterated integrals Chen:1977oja :

[TABLE]

where $\gamma$ is a path with domain $[0,1]$ in the space of external invariants and $\vec{f}(\vec{x}_{0},\epsilon)$ is a vector of boundary conditions. If $\vec{f}^{(i)}(x)$ admits a representation in terms of multiple polylogarithms Goncharov:1998kja ,

[TABLE]

with $G(,t)\equiv 1$, the transformation matrix to the canonical basis is algebraic, and the matrix $\tilde{\mathbf{A}}(x)$ is a $\mathbb{Q}$-linear combination of logarithms

[TABLE]

where $d_{\alpha}$ is the number of linearly independent logarithms, $\mathbf{C}_{i}$ are constant rational matrices, and $\alpha_{i}(\vec{x}),i\in\{1,\dots,d_{\alpha}\}$ are called letters of the differential equations. The set of letters is referred to as the alphabet. If the letters admit a representation in terms of rational functions, then the iterated integrals of (5) can be directly expressed in terms of multiple polylogarithms. On the other hand, if a rational representation is not found, it is often possible to find a polylogarithmic representation by using the knowledge of the symbol of iterated integrals Goncharov:1998kja ; Goncharov:2010jf . The symbol of the $i$-th basis element of (4) at order $\epsilon^{k}$ is obtained by the following recursive formula

[TABLE]

Once the symbol of the solution is known, it is possible to find a corresponding polylogarithmic expression by using an ansatz for the set of polylogarithmic functions, and by imposing boundary conditions (see e.g, Duhr:2011zq ; Bonciani:2016qxi ).

When considering integral sectors that do not admit a polylogarithmic representation, the properties of the transformation matrix to the canonical form are not yet well understood (but see e.g. Adams:2018yfj for recent progress in this direction). For this reason, when a canonical basis is not available, we will only assume that the differential equations are non-singular in the $\epsilon\rightarrow 0$ limit (note that it is always possible to find a basis of Feynman integrals satisfying such differential equations, see e.g. Henn:2014qga ).

3 Series expansion along a contour

Given a base point where the integrals are known, we show how the integrals are computed in a new point by series expanding the integrals along a contour connecting these points. We first define the differential equations along the contour, and then we show how a series expansion of the integrals is found by solving the differential equations around a set of points of the contour. The procedure can be repeated to reach any point of the kinematic regions. In this section we describe the general framework needed to perform series expansions along contours. In appendix B we discuss in detail a simple one-loop example.

Provided that a set of differential equations with respect to a complete set of kinematic invariants is available, we can define the differential equations along a contour $\gamma(t)$ connecting two fixed points $\vec{a}=\{a_{1},\dots,a_{m}\},\;\vec{b}=\{b_{1},\dots,b_{m}\}$. This is achieved by parametrising the contour with a parameter $t$:

[TABLE]

and by considering the differential equations with respect to $t$:

[TABLE]

where the new differential equations matrix can be readily obtained by using the chain rule,

[TABLE]

It is known that a set of master integrals, at a given order of the dimensional regulator, admit a solution in the vicinity of a singular point $\tau$ of the form (see e.g. Smirnov:1990rz )

[TABLE]

where $\vec{c}\;^{(i,j_{1},j_{2},j_{3})}$ are vectors of dimension $n$ (the number of master integrals), $S_{i}$ is a finite set of integers, $w_{k}$ is a complex constant (typically a rational number that accounts for the algebraic dependence of the matrix elements), and $N_{i}$ is the maximal power of the logarithms at order $\epsilon^{i}$. On the other hand, in the vicinity of a regular point $\tau$, the integrals admit a standard Taylor series representation

[TABLE]

More simply, each Feynman integral, in the vicinity of a singular point $\tau$, is expressed as a finite combination of terms of the form

[TABLE]

where $\rho(t)$ is a Taylor series. On the other hand, in the vicinity of a regular point, Feynman integrals admit a Taylor series representation.

In the remainder of this section we show how the series solutions (12) and (13) can be found by series expanding the differential equations matrices, and by iteratively integrating them until the desired order of $\epsilon$.

We know that, for many phenomenologically relevant processes, most integrals admit a polylogarithmic canonical basis. In the next subsection we discuss how we find a power series solution of a canonical set of differential equations. This form of the equations is usually much more compact then the differential equations for a generic basis of master integrals, and working with simplified differential equations renders their series expansion more efficient (for a discussion about timings see Section 4.4). In section 3.2 we discuss how we obtain series solutions for coupled sectors. We remark that the applicability of our method does not rely on (when it exists) a canonical basis of integrals, and a series solution for an arbitrary basis of integrals can be found by using the methods of section 3.2.

3.1 Canonical differential equations

When dealing with single scale problems, generalised power series solutions are usually obtained by defining generic power series, and by fixing the corresponding free coefficients by solving recurrence relations Wasow ; Lee:2017qql ; Mueller:2015lrx . Here we proceed in a more direct way, which is particularly suited for differential equations in canonical form. Specifically, we consider a canonical system of differential equations of the form

[TABLE]

and we assume that the solution is known (analytically or numerically with very high precision, e.g. from a previous expansion) for some $t=t_{0}$. From Eq. (5) it is easy to see that the solution is

[TABLE]

Since we are interested in the solution in the vicinity of a point $\tau$, we series expand the differential equations matrix around $\tau$,

[TABLE]

where $\mathbf{A}^{(i)}_{t}$ are constant matrices 111The series expansion for the matrix $\mathbf{A}^{(i)}_{t}$ can be obtained for example by using the built-in Mathematica function Series.. For a canonical basis, $w_{i},i\in\mathbb{N}$ is expected to be the set of all half integer numbers with $w_{i}\geq-1$. However the following discussion holds for generic complex numbers $w_{i}$. By plugging the previous expansion into (16), we get

[TABLE]

At each order $\epsilon^{k}$, we have to compute integrals of the form,

[TABLE]

These integrals can be computed analytically in terms of integer powers of $\log(t-\tau)$ and (rational) powers of $(t-\tau)$, by using recursively integration-by-parts identities222In practice, these integrals can be computed by using the built-in Mathematica function Integrate.

3.2 Coupled sectors

The problem of finding a canonical basis for integrals that do not admit a polylogarithmic representation is still poorly studied in the literature. In practice, we often deal with differential equations where only a subset of the equations is in canonical form, while the other sectors admit a generic rational dependence on the dimensional regulator and an algebraic dependence on the kinematic invariants. In this case, by iteratively taking derivatives $\partial/\partial_{t}$ of equation (10), it is possible to obtain a $k-$th order differential equation for a single integral, say $f_{j}^{(i)}(t)$, of the form,

[TABLE]

where, for ease of notation, we denoted the generic master integral at order $\epsilon^{i}$ as $g^{(i)}\equiv f_{j}^{(i)}(t)$. The inhomogeneous term $\beta^{(i)}(t)$ is a linear combination of $f_{\neq j}^{(\leq i)}(t)$, and it is known at every iteration $i$. The homogeneous equation associated to (20) is

[TABLE]

and it admits $k$ linearly independent solutions that we denote as $h_{1}(t),\dots,h_{k}(t)$. We note that the homogeneous equation does not depend on the order $\epsilon$ under consideration.

By defining the following fundamental matrix,

[TABLE]

the particular solution of (20) is

[TABLE]

where $\chi^{(i)}_{j},\;j\in\{1,\dots,k\}$ is a set of boundary constants, the Wronskian is defined by $W(\mathbf{\hat{H}}({}s{}))=\text{det}(\mathbf{\hat{H}}(t))$, while $W_{j}(\mathbf{\hat{H}}({}s{}))$ is the determinant obtained by replacing the $j$-th column of $\mathbf{\hat{H}}(t)$ by $(0,\dots,0,1)$. It can be shown that, for functions admitting a series solution of the form (12), the equation in the vicinity of a point $\tau$ can be written as,

[TABLE]

where the functions $b_{i}(t)$ are analytic in $t=\tau$. In this case a series solution in the vicinity of $t=\tau$ can be found by applying the Frobenius method, which is discussed in detail in Appendix A (see also e.g. Coddington ). More specifically, the Frobenius method allows to find a complete set of $k$ homogeneous series solutions in the vicinity of $\tau$. These solutions have the form

[TABLE]

where $c_{i}^{(j_{1},j_{2},j_{3})}$ are complex constants, $K_{i}$ is the maximal power of the logarithm, $S_{i}\subseteq\{1,\dots,k\}$ and $\lambda_{1},\dots,\lambda_{k}$ are the roots of the indicial equation, defined as

[TABLE]

From the form of the homogeneous series solutions, it is clear that the full solution (23) will require, at each order $i$, the computation of integrals of the form (3.1) which, as explained in the previous section, can be done analytically in terms of integer powers of logarithms and complex powers of $t$. This shows that, also in the case of coupled differential equations, we can explicitly find a series representation for Feynman integrals in the vicinity of regular or singular points of the form of Eqs. (13) and (12) respectively.

The Frobenius method is rather standard, and we review its general formulation for linear differential equations of generic order in Appendix A. Here we briefly show its application to second order equations, since this is the order encountered in this paper. Let us consider the case $k=2$ of (21) and, without loss of generality, let us assume that $t=0$ is a (regular) singular point,

[TABLE]

We have to distinguish two cases in order to proceed. Let us first assume that the two roots of the indicial equation $\lambda_{1},\lambda_{2}$ do not differ by an integer number, with $\lambda_{1}>\lambda_{2}$. In this case two linearly independent solutions are

[TABLE]

and the coefficients are fixed by requiring that the differential equation is satisfied order-by-order in $t$. If the two roots differ by an integer, the solution associated to $\lambda_{1}$ is obtained by (28) while the second solution is obtained by

[TABLE]

which can be expressed as a series by expanding all the integrands around $t=0$ and performing the integrations term-by-term.

3.3 Matching

Given a base point where the integrals are assumed to be known, we want to transport the solution to a new point, by using series expansions along a contour connecting these points.

In the previous sections we have seen that, given a singular or regular point of the differential equations, we can find a series solution valid in the vicinity of these points. By construction, these solutions converge only in a region that does not contain any singularity other then the expansion point. In the following we will consider truncated series that, within a given accuracy, provide a good approximation of the full series. When considering a generic contour $\gamma(t)$ for a range of values of $t$, we are interested in finding series expansions along the entire contour. The contour will in general contain multiple singular points of the differential equations, and it is necessary to find multiple series expansions and match them together in order to cover the entire contour. A good criterion for determining the domain of definition of a truncated series is that the series will converge fast enough only when considering it in a region such that the maximum distance from the expansion point is half the distance from the nearest singularity.

Let us assume that we have defined a contour $\gamma(t)$, for $t\in[0,1]$ and we want to find truncated series expansions that approximate the full solution within a given accuracy. $\gamma(0)$ is the known boundary point and $\gamma(1)$ is the target point we want to compute. There might be real and complex singularities. Let us denote the real singularities as $R=\{\tau_{i}|i\in N_{s}\}$ and the complex singularities by $C=\{\lambda_{i}=\lambda^{re}_{i}+i\lambda^{im}_{i}|i\in N_{c}\}$ where $N_{s},N_{c}\in\mathbb{N}$. Moreover we define the following set of real regular points $C_{r}=\bigcup_{i=1}^{N_{c}}\{\lambda^{re}_{i}-\lambda^{im}_{i},\lambda^{re}_{i},\lambda^{re}_{i}+\lambda^{im}_{i}\}$. We now consider a set of points $t_{i}\in R\bigcup C_{r}$ such that $t_{-1}<0<t_{1}<\dots<t_{N_{e}}<1<t_{N_{e}+1}<t_{N_{e}+2}$ with $N_{e}\in\mathbb{N}$. We then find truncated power series around $t_{0}=0$ and $t_{i}$ for $i\in\{1,\dots,N_{e}+1\}$. We denote these series as $\vec{f}_{[\gamma]}(t,\epsilon)_{[t_{i}-r_{i},t_{i},t_{i}+r_{i}]}$, where $t_{i}$ is the expansion point, and $r_{i}$ is the radius of the series defined as the distance between $t_{i}$ and the closest point $t_{j\neq i}$. Moreover we define $\vec{f}_{[\gamma]}(t,\epsilon)_{[t_{i}-r_{i},t_{i},t_{i}+r_{i}]}$ to be equal to the truncated series for $t_{i}-r_{i}\leq t\leq t_{i}+r_{i}$, while it is zero otherwise. Since the $t_{i}$ are in general not equally spaced it can happen that some segments of the contour are not cover by any series. In this case we iteratively introduce new (regular) expansion points $\kappa_{i},i\in\mathbb{N}$ , and corresponding series $\vec{f}_{[\gamma]}(t,\epsilon)_{[\kappa_{i}-\rho_{i},\kappa_{i},\kappa_{i}+\rho_{i}]}$, such that $\kappa_{i}$ is in the middle of an uncovered region and $\rho_{i}$ is the distance between $\kappa_{i}$ and the closest point among $\kappa_{j\neq i},t_{j}$. The procedure is repeated until the entire contour is covered.

There are two special cases that need some care. If there are no points $t_{i}$, we just expand around $t_{0}=0$ and the entire contour will be covered. If there are no points $t_{N_{e}+1},t_{N_{e}+2}$ we iteratively add new (regular) expansion points such that $0<k_{i}<t_{N_{e}}$. At the end of this procedure it can happen that the point $t=1$ in not covered. We then set $t_{N_{e}+1}=1$ and we add new (regular) expansion points as described above until the entire contour is covered.

We finally define the solution along the entire contour to be,

[TABLE]

$r^{\prime}_{i}\leq r_{i},\rho^{\prime}_{i}\leq\rho_{i}$ are such that there are no overlapping series, and $S_{e}\subseteq\{1,\dots,N_{e}\},S_{r}\subseteq\{1,\dots,N_{r}\}$ ($N_{r}$ is the total number of regular expansion points) are such that no series outside the unit interval $t\in[0,1]$ is considered in the sum.

In Sections 3.1 and 3.2 we have seen that the series solutions depend on a set of integration constants that have to be fixed by imposing boundary conditions. When considering a set of series along a contour, the integration constants of one series are fixed by knowing the boundary conditions at a given point (e.g. when the boundary conditions are known analytically by other means, or they are known because the contour under consideration intersects another contour along which the series expansion is already known), while the other series are fixed by imposing that two consecutive series have the same value in the contact point, i.e. the point where they are both defined.

4 A planar elliptic family for Higgs+jet production

The planar two-loop QCD corrections to Higgs+jet productions are mediated by the four integral families depicted in Fig. 1. These integrals were computed in Bonciani:2016qxi in the non-physical region, where all the Mandelstam invariants are negative reals. For the polylogarithmic sectors, the solution was expressed in terms of logarithms and dilogarithms up to weight two, while the weight-three and four solutions were expressed in terms of one-fold integrals over the lower weight solutions. Family A contains two elliptic sectors. The first elliptic sector is $I^{A}_{1,1,0,1,1,1,1,0,0}$; the homogeneous solutions are elliptic integrals, and the solution in dimensional regularisation requires iterated integrations over elliptic kernels. The second elliptic sector is $I^{A}_{1,1,1,1,1,1,1,0,0}$, and it admits a canonical form for the homogeneous differential equations, but it is coupled to the first elliptic sector via inhomogeneous terms. In Bonciani:2016qxi the solution for the elliptic sectors was expressed in terms of iterated integrals over elliptic kernels.

In this section we obtain generalised power series expansions for family A in the $p_{2}\leftrightarrow p_{3}$ channel, which exhibits the most complicated (spurious, see section 4.3) singularity structure of the channels needed for the scattering amplitude. All the other families are simpler, as they do not involve elliptic structures and can be solved with the same methods.

The family under consideration is defined by the following set of integrals

[TABLE]

with

[TABLE]

$d_{1}$-$d_{7}$ are propagators while $d_{8}$ and $d_{9}$ are numerators, therefore $a_{1},\dots,a_{7}$ are non-negative integers while $a_{9},a_{8}$ are non-positive integers. The kinematics is such that $p_{1}^{2}=p_{2}^{2}=p_{3}^{2}=0$ with

[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 propagators. The relevant physical region is defined by

[TABLE]

We solve the kinematic constraints for $s_{12},s_{13},p_{4}^{2}$ and we define the following scaleless variables

[TABLE]

The family contains 73 master integrals, and our choice for the integral basis is provided in the ancillary files of the arXiv submission. The first 65 master integrals satisfy a set of differential equations in canonical form,

[TABLE]

The alphabet of the differential equations consists of 42 letters, depending on the following set of 6 square roots333For ease of notation we joined products of square roots into one square root. However, as discussed in Section 4.2, in order to perform the analytic continuation, it is convenient to split them in square-root factors. The actual factorisation of the square roots used in this paper is the one of the integral basis provided in the ancillary files of the arXiv submission. We also note that, since the square roots appear as factors in the basis choice, one has the freedom to choose their factorisation provided that, once one choice is made, all the subsequent operations are carried out consistently with this choice.,

[TABLE]

We remark that our approach does not rely on the rational parametrisation of the set of square roots, and it works for general algebraic dependence of the differential equations. Integrals 66-73 are elliptic and satisfy coupled differential equations of the form

[TABLE]

where the vector $\vec{g}_{66-73}(\vec{x},\epsilon)$ depends on the polylogarithmic integrals $\vec{f}_{1-65}(\vec{x},\epsilon)$. Order-by-order in $\epsilon$ the polylogarithmic integrals $\vec{f}_{1-65}(\vec{x},\epsilon)$ satisfy completely decoupled differential equations. On the other hand, the elliptic sectors are coupled and, in general, one needs to solve higher order differential equations for single integrals. Specifically, the homogeneous matrix $\mathbf{B}^{(0)}_{x_{i}}$ has the following schematic form

[TABLE]

where the lines separate sector $I^{A}_{1,1,0,1,1,1,1,0,0}$ (first four rows) from sector $I^{A}_{1,1,1,1,1,1,1,0,0}$ (last four rows). We see that integrals 67 and 69 are coupled. For each of these integrals we can define a second order differential equation. Once integral 67 or 69 is known, all the other integrals are decoupled and can be solved by considering first order equations only.

4.1 Series solution of the differential equations

In order to solve the differential equations we need a set of boundary conditions. We use the point $x_{1}=x_{2}=x_{3}=0$, and the boundary conditions are Bonciani:2016qxi

[TABLE]

We are interested in a solution in the relevant physical region (34). Therefore we transport the boundary condition to the point $\left\{2,\;\frac{13}{25},\;-1\right\}$ by series expanding the integrals along the contour444In this paper we assume $p_{4}^{2}=13/25\approx(125/172)^{2}$, which is a good approximation of the Higgs-boson mass, when normalised to the top mass. However, the applicability of the method does not rely on the choice of this numerical value. For example, we could equally well consider values of the dimensionless ratios related to the bottom mass. More generally, we could apply the expansion method by considering the Higgs mass as an additional variable.,

[TABLE]

Here we discuss in detail how we expand along the contour $\gamma_{\text{thr}}(t)$ defined as

[TABLE]

This contour is interesting because it allows to analytically continue the integrals above the physical threshold $x_{1}=4$. In order to achieve the decomposition of the contour described in Section 3.3, we need to know where the singularities lie. When considering differential equations with algebraic dependence (square roots in the present case), the singularities are all the non-analytic points of the differential equations, i.e. points where the differential equations diverge but also the zeros of the square roots (branching points). Once the path is fixed, the problem is one dimensional and we can find the position of the singularities by solving for the zeros of the denominators of the differential equations and for the zeros of the square roots (see Section 4.2 for a detailed discussion about the different classes of singular points of the differential equations).

For the path $\gamma_{\text{thr}}(t)$, the relevant singularities are

[TABLE]

where $\tau_{1}$ is the physical threshold and $\tau_{2}$ is a spurious singularity outside the $[0,1]$ interval. We now partition the contour as described in Section 3.3. Specifically, we add the regular expansion points $\kappa_{1}=0,\;\kappa_{2}=1\;\kappa_{3}=\frac{1}{8}$ which results in the following partitioning of the contour,

[TABLE]

where we replaced rational numbers with (exact) real numbers. All the boundary constants are fixed by imposing the following chain of boundary conditions,

[TABLE]

4.2 Analytic continuation

In the previous sections we showed that in order to obtain converging power series expansions along a contour it is necessary to expand around the singular points of the differential equations (and regular points in order to ensure fast converging series, as described in Sec. 3.3). As already mentioned, by singular point we mean any non-analytic point of the differential equations and of the solution, i.e. power and logarithmic divergences, or branching points of the square roots. The singularities of the solution are a subset of the singularities of the differential equations. In this section we discuss the different classes of singularities encountered in the (series) solution of the differential equations and how to perform the analytic continuation across them. Moreover, for the sake of the following discussion, we remark that we solve the differential equations along contours entirely contained in the physical sheet.

The first class of singularities are the so-called physical singularities, which correspond to the singularities predicted by unitarity cuts. In correspondence of physical singularities the solution develops branch cuts, and in order to analytically continue the solution across them we use Feynman prescription, i.e. we assign a small imaginary part to the contour parameter in such a way that the invariant crossing the singularity acquires a small positive imaginary part. In order to show this, let us write down explicitly the first few terms of the expansion around the heavy-quark threshold $t=\frac{1}{2}$, discussed in the previous section, for, e.g., the (elliptic) integral 68 at order $\epsilon^{4}$ which reads555In what follows we show truncated numerical coefficients, while the full coefficients are computed with hundreds of digits. Keeping the coefficients numeric allows to drastically speed up the series expansion, especially at high truncation order.,

[TABLE]

When expanding along a path crossing a physical singularity, in this case $x_{1}=4$, the branch cuts are expressed by rational powers of the expansion parameter and logarithms. The analytic continuation is then performed by assigning a small imaginary part to the parameter $t$, consistently with Feynman prescription. In this case, since $x_{1}(t)=2+4t$, we have $t\rightarrow t+i\delta$, and,

[TABLE]

On the other hand, the expansion $\vec{f}_{[\gamma_{\text{thr}}]}(t,\epsilon)_{[0.81,\;1,\;1]}$ is a regular power series, as $t=1$ is a regular point, and the branch cut ambiguities of the solution in this region are fixed by imposing the boundary conditions from the analytically continued $\vec{f}_{[\gamma_{\text{thr}}]}(t)_{[0.19,\;0.5,\;0.81]}$, as in Eq. (4.1). In Fig. 2 we present the plots of integrals $\vec{f}_{66-73}$ along the contour $\gamma_{\text{thr}}$.

The second class of singularities are the so-called non-physical singularities, which are not physical (they are not predicted by unitarity) but they appear in the solution of a given integral basis. In order to explain the origin of these singularities let us consider a generic basis element of the integral basis666In this work we consider only integral basis with algebraic coefficients.,

[TABLE]

where $I_{a_{1ij},\dots,a_{nij}}$ is a scalar Feynman integral of the form,

[TABLE]

If one of the coefficients of Eq. (48), say $c_{ik}(\vec{x})$, is singular for $\vec{x}=\vec{x}_{s}$, where $\vec{x}_{s}$ is not a physical singularity and $I_{a_{1ik},\dots,a_{nik}}$ is non-zero in $\vec{x}_{s}$, $\vec{B}(\vec{x})$ is singular for $\vec{x}=\vec{x}_{s}$. This is what we call a non-physical singularity, since it originates only from the prefactors of the integral basis777In general, multiple coefficients $c_{ij}(\vec{x})$ might be singular in $\vec{x}_{s}$, or the corresponding $I_{a_{1ij},\dots,a_{nij}}$ might vanish in the same point. In this case it is not obvious whether or not $\vec{x}_{s}$ is a singular point for $\vec{B}(\vec{x})$, since the different singular terms might cancel out, and only after the series solution is obtained one can verify whether or not $\vec{x}_{s}$ is a singular point for $\vec{B}(\vec{x})$.. For integral bases with algebraic coefficients (as the ones considered in this paper) there are two kind of non-physical singularities: poles and branching points of square roots. Poles do not give rise to cuts, and no analytic continuation is needed across them. On the other hand, for square roots giving rise to non-physical branch cuts, we can choose an arbitrary analytic continuation across their branching points, since these cuts cancel at the level of the integrals $I_{a_{1ij},\dots,a_{nij}}$ and at the level of the scattering amplitude. More specifically, for the integrals of family A the only physical singularities are for $x_{1}=4$ and $x_{2}=4$. According to Feynman prescription we define $\sqrt{4-x_{1}}=-i\sqrt{x_{1}-4}$ and $\sqrt{4-x_{2}}=-i\sqrt{x_{2}-4}$ when $x_{1},x_{2}>4$. On the other hand, all the other square roots have no $x_{1}-4$, $x_{2}-4$ factors and we choose, for simplicity, to define them as $\sqrt{a}=i\sqrt{-a}$ when $a<0$, for any $a$ and in any region. In Appendix B we discuss a simple one-loop example and we show how non-physical singularities appear in the series solution.

Finally we have the so-called spurious singularities. These are singularities of the differential equations that do not correspond to singularities of the solution. These singularities are not physical and the coefficients of the integral basis are regular in these points. Therefore, since we consider contours entirely contained in the physical sheet, the solution is regular in these points, and the singular terms of the series solution corresponding to spurious singularities cancel within the truncation error, and no analytic continuation is needed.

4.3 Mapping the physical region to a finite region

In physical applications it is often necessary to consider integrals for very large or small values of the external invariants. When working with series expansions it is then convenient to map the relevant physical region to a finite region, so that it is possible to expand around limiting points in a straightforward manner. For two-to-two processes a convenient set of variables is the following,

[TABLE]

which, by keeping the Higgs mass fixed, maps the physical region (34) to the unit square in the $l,z$ space. The $l,z$ variables are also convenient when studying the singularities of the differential equations. The singularities are the zeros of the denominators of the differential equations and the zeros of the square roots. For our differential equations, the singular points are for $l=0,\;l=1,\;z=0,\;z=1$ and along the orbits defined by the following equations,888The singular orbits can be found, e.g., by using the built-in Mathematica function Solve.

[TABLE]

The orbit $z=0$ corresponds to points where $x_{1},x_{3}$ are infinity and these are singular points for our integral basis (see below). $l=0$ and $z=1$ correspond to $x_{3}=0$ and are not physical singularities but the series expansions are not analytic there because some of the square roots vanish. $l=1$ is a spurious singularity. The first singularity of Eq. (51) is the physical threshold. The second and last singular orbits of Eq. (51) are spurious. Indeed, the singular orbits do not correspond to physical singularities, and none of the square roots or the rational prefactors of the integral basis are singular along these orbits. The physical region, including the singular orbits, is represented in Fig. 3.

The $l,z$ variables are also the natural variables to perform expansions along contours reaching very large or small values of the invariants. Let us consider for example the contour,

[TABLE]

which corresponds to a contour from $\{x_{1}=\frac{24}{5},\;x_{2}=\frac{13}{25},\;x_{3}=-\frac{18083}{41040}\}$ to $\{x_{1}=\infty,\;x_{2}=\frac{13}{25},\;x_{3}=-\frac{169}{342}\}$. In this case the only singular point is for $t=1$, and it is sufficient to perform only one expansion around it to cover the full contour. $t=1$ corresponds to a (singular) point at infinity. Nonetheless, we never cross the points at infinity and possible branch cut ambiguities are fixed imposing boundary conditions at a finite point and by treating the square roots as explained in Section 4.2. The plots of the elliptic integral sectors along $\gamma_{\infty}$ are presented in Fig. 4.

4.4 Numerical results and timings

In this section we provide numerical results with accuracy and timings for the elliptic integrals $\vec{f}_{66-73}$, which are the most complicated ones. We tested different ways to estimate the accuracy of the numerical results obtained by truncating the relevant series. One estimate of the error is obtained by expanding along two different paths reaching the same end point, and then taking the difference of the results. Another way is to keep the path fixed but truncating the series at different orders. In all the cases we analysed these methods give the same estimate of the error. In Table 1 we provide a comparison of our expansion method against the c++ version of FIESTA 4.1 Smirnov:2015mct and pySecDec 1.4.3 Borowka:2017idc ; Borowka:2018goh . We remark that the timings include the time needed, starting from the set of differential equations (36) and (38), to define the differential equations along the contour, series expand the matrix elements and recursively integrate them up to the desired $\epsilon$ order, for all the expansion points along the contour. This is an accurate estimate of the timings for physical applications where, in principle, one has to repeat the procedure for each target point of interest. In Fig. 5 we report the time required to compute all the integrals up to and including the sectors defined in Eq. (4.4). The sectors are ordered according to their position in the integral basis under consideration and, at the level of the differential equations, they are coupled to one or more sectors with lower indices, and none of the sectors with higher indices. We note that most of the time is spent expanding the elliptic sectors, due to the lack of a canonical (hence simpler) basis for those integrals. In Table 2 we provide high-precision numerical results for integrals $\vec{f}_{66-73}$ at order $\epsilon^{4}$.

Finally, we remark that our method is easily parallelizable since, given one or more boundary points, the series expansion along different contours are completely independent operations.

[TABLE]

5 Conclusion

We showed that by defining the differential equations for a set of multi-scale Feynman integrals along contours connecting two generic points of the kinematic regions, it is possible to systematically obtain one-dimensional series expansions for the integrals along the contours. Specifically, we applied this method to obtain new results for a planar family of integrals relevant for the two-loop QCD corrections to Higgs + jet production. When the expansion is performed along a contour such that only one invariant changes along it, the analytic continuation above the physical thresholds becomes straightforward. We demonstrated that performing an expansion along a contour is fast, and makes it possible to repeat the procedure to compute the integrals over the entire kinematic regions, with arbitrary numerical precision. Our approach is algorithmic, and it seems plausible to implement it in a computer code that can be applied to solve complicated integrals of phenomenological integrals.

Acknowledgements

We thank Roberto Bonciani, Vittorio Del Duca, Hjalte Frellesvig, Johannes Henn, Martijn Hidding, Leila Maestri, Giulio Salvatori and Volodya Smirnov for useful discussions and collaboration on related work. We thank Roman Lee for insightful discussions about series expansion methods. We thank Armin Schweitzer for useful discussions about the analytic continuation of Feynman integrals of elliptic type, and for providing the (unpublished) canonical basis of the one-loop box example discussed in Appendix B. This research received funding from the European Research Council (ERC) under grant agreement No. 694712 (pertQCD) and by the Swiss National Science Foundation project No. 177632 (ElliptHiggs).

Appendix A General formulation of the Frobenius method

The Frobenius method has general validity, and in this appendix we consider a generic order-$k$ linear differential equation for an unknown function $f(t)$ (the discussion closely follows Coddington ). The Frobenius method can be used to find a complete set of homogeneous series solutions to the equation in the vicinity of a singular (or regular) point that we assume being located at $t=0$. The homogeneous equation can be written in general as:

[TABLE]

Since the $b_{i}(t)$ are analytic in $t=0$ they admit a Taylor series representation of the form,

[TABLE]

We look for a solution of the form,

[TABLE]

where lambda is a (complex in general) parameter to be fixed. By substituting the formal series solution into the equation we obtain,

[TABLE]

Where the $g(j)$ are linear in the $c^{(1)},c^{(2)},\cdots,c^{(j-1)}$ with polynomial coefficients in $\lambda$, while the polynomial $f(\lambda)$ is called the indicial polynomial and it reads,

[TABLE]

In order for $f(t)$ to be a solution, the coefficients of the right-hand side of (57) need to satisfy

[TABLE]

which is a recursion relation that can be solved for $c^{(j)}$ with $j\geq 1$ if $\sigma(\lambda+j)\neq 0$ for any $j\geq 1$. The $c^{(j)}$, $j\geq 1$ are then uniquely determined functions of $c^{(0)}$. In this case we get

[TABLE]

It is convenient to fix the value of $c^{(0)}$, and we conventionally choose $c^{(0)}=1$. If $\lambda=\lambda_{1}$ where $\lambda_{1}$ is a root of the indicial polynomial, $\sigma(\lambda_{1})=0$, and $\sigma(\lambda+j)\neq 0$ for any $j\geq 1$ then, $f_{\lambda_{1},1}(t)\equiv f(t)|_{\lambda=\lambda_{1}}$, is a solution of the equation $L[f(t)]=0$. If $\lambda_{1}$ has multiplicity 2, we need to find a second solution associated to it. If we differentiate Eq. (60) with respect to $\lambda$ we get,

[TABLE]

If $\lambda_{1}$ is a root of multiplicity 2 then $\sigma(\lambda_{1})=\left.\frac{\partial\sigma(t)}{\partial\lambda}\right|_{\lambda=\lambda_{1}}=0$, and $\frac{\partial f(t)}{\partial\lambda}$ is also a solution of Eq. (54). We have then that the second solution associated to the root $\lambda_{1}$ is

[TABLE]

If $\lambda_{1}$ has multiplicity $m_{1}$, all the $m_{1}$ solutions associated to $\lambda_{1}$ are obtained by taking $m_{1}-1$ derivatives.

Let us now suppose that $\lambda_{1}+k_{2}=\lambda_{2}$, for a positive integer $k_{2}$, i.e. the root $\lambda_{1}$ differs by an integer $k_{2}$ from another root $\lambda_{2}$. Moreover we assume that $\lambda_{2}$ is the only root that differs from $\lambda_{1}$ by an integer. In this case the solution associated to $\lambda_{1}$ cannot be determined as explained above, since the recursion (59) becomes ill defined when $j=k_{2}$. In order to find a solution we set $c^{(0)}=(\lambda-\lambda_{1})^{m_{2}}$, where $m_{2}$ is the multiplicity of the root $\lambda_{2}$,

[TABLE]

An explicit computation shows that Eq. (60) becomes

[TABLE]

and the recursion (59) can be written as,

[TABLE]

while

[TABLE]

and

[TABLE]

where the functions $G^{(i)}$ are non-zero for $\lambda=\lambda_{1}$. We note that the right hand side of (66) is now well defined when $\lambda$ approaches $\lambda_{1}$ since $\sigma(\lambda_{1}+k_{2})=\sigma(\lambda_{2})=0$ and the denominator has a factor $(\lambda-\lambda_{1})^{m_{2}}$ that cancels against the numerator. This shows that $f(t)|_{\lambda=\lambda_{1}}$ with $f(t)$ defined in Eq. (63) is a solution of Eq. (54). However we notice that $c^{(i)}=0$ for all $i\in\{1,\dots,k_{2}-1\}$ so that the solution has the form

[TABLE]

The leading term of the series is $t^{\lambda_{2}}$, therefore this solution is linearly dependent from the solution one would obtain by considering the root $\lambda_{2}$. In order to find a solution truly associated to $\lambda_{1}$ we take the $m_{2}$-th derivative of $f(t)$, which satisfies

[TABLE]

and $\frac{\partial^{m_{2}}f(t)}{\partial\lambda^{m_{2}}}|_{\lambda=\lambda_{1}}$ is indeed a solution. Moreover its leading term is $m!t^{\lambda_{1}}$ so that $\frac{\partial^{m_{2}}f(t)}{\partial\lambda^{m_{2}}}|_{\lambda=\lambda_{1}}$ is a genuine solution associated to $\lambda_{1}$. If $\lambda_{1}$ has multiplicity $m_{1}>1$ the other solutions are obtained by taking $m_{1}-1$ derivatives of $\frac{\partial^{m_{2}}f(t)}{\partial\lambda^{m_{2}}}$ with respect to $\lambda$ evaluated at $\lambda_{1}$.

The last case we need to consider is when there are more than one root differing from $\lambda_{1}$ by an integer. More precisely, we consider $r-1$ roots $\lambda_{2},\dots,\lambda_{r}$ with $Re(\lambda_{2})<\dots<Re(\lambda_{r})$ such that $\lambda_{1}+k_{i}=\lambda_{i}$ with $k_{i}$ a positive integer. Each root has multiplicity $m_{i}$ and we define $M=\sum_{i=2}^{r}m_{i}$. We set $c^{(0)}=(\lambda-\lambda_{1})^{M}$, so that

[TABLE]

By proceeding as in the case of only two roots differing by an integer, one obtains that a solution associated to $\lambda_{1}$ is given by

[TABLE]

while if $m_{1}>1$ the other solutions are obtained by taking $m_{1}-1$ derivatives.

Appendix B One-loop example

In this appendix we apply the methods of Sec. 3 to a one-loop integral family, and we discuss in detail the derivation of the series expansion along a given contour. Specifically, we consider the integral family relevant for the leading order corrections to Higgs + one jet production in QCD depicted in Fig. 6,

[TABLE]

with

[TABLE]

while the kinematics is $p_{1}^{2}=p_{2}^{2}=p_{3}^{2}=0$ with,

[TABLE]

where $p_{4}^{2}$ is the squared Higgs-mass, and $m^{2}$ is the squared mass of the propagators. We solve the kinematic constraints for $s_{12},s_{13},p_{4}^{2}$. A canonical basis for the integral family is given by,

[TABLE]

where we have used the following definitions for the square roots,

[TABLE]

The corresponding differential equations are,

[TABLE]

where the matrix $\tilde{A}$ is,

[TABLE]

and its matrix elements,

[TABLE]

The letters are defined as

[TABLE]

In order to discuss the solution we define, as in the main text, the following scaleless variables,

[TABLE]

and we set $m^{2}=1$. A convenient boundary point is $\vec{x}=\vec{0}$, where we have

[TABLE]

We consider the contour,

[TABLE]

The differential equations along the contour are,

[TABLE]

It is easy to see that the differential equations along the contour are singular for $t=0$ and $t=2/3$, the latter being the physical threshold $s_{12}=4$ while the former is a non-physical singularity (see discussion below). Following the prescription of Sec. 3.3 we have ,

[TABLE]

We start with the computation of $\vec{B}_{[\gamma_{1}]}(t,\epsilon)_{[0,0,1/3]}$. We series expand $\frac{\partial\tilde{A}(\vec{x}(t))}{\partial t}$ around $t=0$. By using the shorthand $a^{\prime}_{i,j}\equiv\frac{\partial a_{ij}}{\partial t}$ we have999As in the main text, we show truncated numerical coefficients, while the full coefficients are computed with hundreds of digits.,

[TABLE]

while the other matrix elements are identically zero. We are now able to find the series solution $\vec{B}_{[\gamma_{1}]}(t,\epsilon)_{[0,0,1/3]}$ by using the recursion,

[TABLE]

where $C^{(i)}$ is the boundary vector and the integral on the right-hand-side is an indefinite integral. We remark that by series expanding $\frac{\partial\tilde{A}(\vec{x}(t))}{\partial t}$ all the integrals in (87) are elementary and can be performed analytically. At order $\epsilon^{1}$ we have

[TABLE]

The order $\epsilon^{0}$ of the solution is given by the boundary conditions (82). By performing the integrals on the right-hand-side and by imposing the boundary condition $\vec{B}^{(1)}_{[\gamma_{1}]}(0)_{[0,0,1/3]}\equiv\vec{0}$ we get,

[TABLE]

Proceeding in the same way, at order $\epsilon^{2}$ we get,

[TABLE]

It is straightforward to iterate the procedure up to arbitrary order of $\epsilon$. We see that the solution is singular in $t=0$, since it develops a brunch cut for $t<0$. This is what we call, in Section 4.2, a non-physical singularity, since the only physical singularity along the contour is $s_{12}=4m^{2}$. The appearance of this singularity is due to the prefactors defining the canonical basis (B) and we see that, by inverting the basis, the integrals $I_{a_{1},a_{2},a_{3},a_{4}}$ admit a regular Taylor series expansion around $t=0$. We remark that, according to (B), the differential equations depend originally on $\sqrt{-t}$ and, as discussed in Section 4.2, we consider the standard brunch for $\sqrt{-t}$, i.e. $\sqrt{-t}>0$ for $t<0$ and $\sqrt{-t}=i\sqrt{t}$ for $t>0$.

The computation of $\vec{B}_{[\gamma_{1}]}(t,\epsilon)_{[1/3,2/3,1]}$ follows the same steps. One first expands $\frac{\partial\tilde{A}(\vec{x}(t))}{\partial t}$ around $t=2/3$. The series solution is then obtained by,

[TABLE]

and by using as boundary condition the value of the integrals in the contact point of the two series, in this case $t=1/3$,

[TABLE]

where the right-hand-side is known from the previous expansion. The first non-trivial $\epsilon$ orders of the solution are,

[TABLE]

and

[TABLE]

where $t^{\prime}=2/3-t$. We see that for $t>2/3$ the solution develops a brunch cut, and the cut ambiguity is resolved by Feynman prescription which, by (83), is $t\rightarrow t+i\delta$ with $\delta$ a small positive imaginary part.

Appendix C Plots

In this appendix we show plots of all the 73 master integrals at order $\epsilon^{4}$. The plots are obtained by series expanding along the contour $\gamma_{\text{thr}}$ defined in section 4.1 ($s_{13}=-1,\;p_{4}^{2}=\frac{13}{25},\;m^{2}=1$, $2\leq s_{12}\leq 6$). The solid dots represent numerical values computed with FIESTA 4.1. The real part of the integrals is in blue, the imaginary part is orange.