## 1 Introduction

Many single scale Feynman integrals arising in massless and massive multi-loop calculations in Quantum Chromodynamics (QCD) [QCD] have been found to be expressible in terms of harmonic polylogarithms (HPLs) [Remiddi:1999ew], generalized harmonic polylogarithms [Moch:2001zr, Ablinger:2013cf], cyclotomic harmonic polylogarithms [Ablinger:2011te], square-root valued iterated integrals [Ablinger:2014bra], as well as more general functions, entering the corresponding alphabet in integral iteration. After taking a Mellin transform

(1.1) |

they can be equivalently expressed in terms of harmonic sums [Blumlein:1998if, Vermaseren:1998uu] in the simpler examples and finite sums of different kinds in the other cases [Moch:2001zr, Ablinger:2013cf, Ablinger:2011te, Ablinger:2014bra], supplemented by special numbers like the multiple zeta values [Blumlein:2009cf] and others appearing in the limit of the nested sums, or the value at of the iterated integrals in Refs. [Blumlein:1998if, Vermaseren:1998uu, Remiddi:1999ew, Moch:2001zr, Ablinger:2013cf, Ablinger:2011te, Ablinger:2014bra].

In many higher order calculations a considerable reduction of the number of integrals to be calculated is obtained using integration by parts identities (IBPs) [IBP], which allow to express all required integrals in terms of a much smaller set of so called master integrals. Differential equations satisfied by these master integrals [DEQ, Caffo:1998du] can then be obtained by taking their derivatives with respect to the parameters of the problem and inserting the IBPs in the result. What remains is to solve these differential equations, given initial or boundary conditions, if possible analytically. One way of doing this is to derive an associated system of difference equations [DIFFEQ, Ablinger:2010ty, Ablinger:2014nga, Ablinger:2015tua] after applying a mapping through a formal Taylor series or a Mellin transform. If these equations factorize to first order equations, we can use the algorithm presented in Ref. [Ablinger:2015tua] for general bases to solve these systems analytically and to find the corresponding alphabets over which the iterated integrals or nested sums are built. The final solution in and space is found by using the packages Sigma [SIG1, SIG2], EvaluateMultiSums and SumProduction [EMSSP].

However, there are physical cases where full first order factorizations cannot
be obtained for either the differential equations in or the difference equations in .^{1}^{1}1There
are cases in which factorization fails in either - or -space, but not in
both, cf. [Ablinger:2013eba]. This opportunity has to be always checked. The next level of complexity
is given by non-factorizable differential or difference equations of second order. Examples of this are the massive
sunrise and kite integrals [SABRY, TLSR1, Caffo:1998du, Caffo:2002ch, Laporta:2004rb, TLSR2, TLSR2a, Bailey:2008ib, Broadhurst:2008mx, TLSR3, BLOCH2, TLSR4, Adams:2014vja, Adams:2015gva, Adams:2015ydq, TLSR5, Remiddi:2016gno, Adams:2016xah, TSLR6].

In the present paper, we will address the analytic solution of typical cases of this kind, related to a series of
master integrals appearing in the 3-loop corrections of the -parameter in [Grigo:2012ji]. It turns
out that these integrals are more general than those appearing in the sunrise and kite diagrams, due to the
appearance of also the elliptic integral of the second kind, , which cannot be transformed away.
The corresponding second order differential equations have more than three singularities, as in the case
of the Heun equation [HEUN]. For the sake of generality, we will seek solutions of the second order
homogeneous differential equations which are given in terms of Gauß’ functions [GAUSS] within
the class of globally bounded solutions [VANH1], cf. also [BOSTAN]. Here the parameters of the
function are rational numbers and the argument is a rational function of . The complete elliptic integrals
and [TRICOMI, ELLINT, ABEL, JAC1]^{2}^{2}2As a convention, the modulus is
chosen in this paper, also used within the framework of Mathematica. are special cases of this class.

The hypergeometric function obeys different relations like the Euler- and Pfaff-transformations [HYP, SLATER], the 24 Kummer solutions [KUMMER, RIEMANN] and the 15 Gauß’ contiguous relations [HYP, SLATER]. There are more special transformations for higher than first order in the argument [KUMMER, GOURSAT, VANHVID, VID2]. In the present case, equivalent representations are obtained by applying arithmetic triangle groups [TAKEUCHI]. The corresponding algorithm has been described in Ref. [IVH] in its present most far reaching form. The relations of this type may be useful to transform a found solution into another one, which might be particularly convenient. In the case a function space of more solutions is considered, these relations have to be exploited to check the independence of the basis elements.

The main idea of the approach presented here is to obtain the factorization of a high order scalar difference or
differential equation, after uncoupling [Zuercher:94, ORESYS, NewUncouplingMethod] the corresponding linear
systems, to all first order parts and its second order contributions. While the first order parts have been
algorithmically
solved in Ref. [Ablinger:2015tua], the treatment of second order differential equations shall be
automated.^{3}^{3}3For more involved physical problems also irreducible higher order differential equations may
occur. The class of solutions has an algorithmic automation to a wide extent [IVH] and it seems that
this class constitutes the next one following the iterated square-root valued letters in massive single-scale
3-loop integrals. Applying this method, we obtain the corresponding functions with (partly) fractional
parameters and rational argument, and ir(rational) pre-factors, forming the new letters of the otherwise
iterated integrals. These letters contain a definite integral by virtue of the integral representation of
the function, which cannot be fully transformed into an integral depending on the follow-up integration
variable only through its integration boundaries. In general, we have therefore to iterate new letters of this kind.
Through this we obtain a complete algorithmic automation of the solution also when second order
differential operators contribute, having solutions.

As it will be shown, in a series of cases the reduction of the functions to complete elliptic integrals and is possible. Therefore we also study special representations in terms of -series, which have been obtained in the case of the sunrise graph, cf. [BLOCH2, Adams:2014vja, Adams:2015gva, Adams:2015ydq, Adams:2016xah], before. More general representations are needed for the integrals considered in the present paper and we describe the necessary extension.

In performing a higher loop calculation, in intermediary steps usually more complicated nested
integrals and sums occur than in the final result^{4}^{4}4For a simple earlier case, see e.g.
[Vermaseren:2005qc, Ablinger:2010ty].. The various necessary decompositions of the problem that have
to be performed, such as the integration by parts reduction and others, account for this in part. It appears
therefore necessary to have full control on these occurring structures first, which finally may simplify
in the result. Moreover, experience tells that in more general situations, more and more of these structures
survive, cf. [Ablinger:2014nga] in comparison to [Vermaseren:2005qc]. If the mathematical properties
of the quantities occurring are known in detail, various future calculations in the field will be more easily
performed.

The paper is organized as follows. In Section 2 we present the linear systems of first order differential equations for master integrals in Ref. [Grigo:2012ji] which cannot be solved in terms of iterated integrals. We first perform a decoupling into a scalar second order equation and an associated equation for each system. Using the algorithm of Ref. [Ablinger:2015tua], the non-iterative solution both in - and Mellin--space is uniquely established. In Section 3 we first determine the homogeneous solutions of the second order equations, which turn out to be solutions [VANH1] and obey representations in terms of weighted complete elliptic integrals of first and second kind at rational argument. In Section 4 we derive the solutions in the inhomogeneous case, which are given by iterated integrals in which some letters are given by a higher transcendental function defined by a non-iterative, i.e. definite, integral in part. We present numerical representations for deriving overlapping expansions around and . The methods presented apply to a much wider class of functions than the ones being discussed here specifically. These need neither to have a representation in terms of elliptic integrals, nor of just a function. The respective letter can be given by any multiple definite integral.

Owing to the fact that we have elliptic solutions in the present cases we may also try to cast the solution in terms of series in the nome

(1.2) |

where

(1.3) |

denotes the ratio of two complete elliptic integrals of first kind and is a rational function associated to the elliptic curve of the problem. It is now interesting to see which closed form solutions the corresponding series in obey. All contributing quantities can be expressed in terms ratios of the Dedekind function [DEDEKINDeta], cf. Eq. (6.14). However, various building blocks are only modular forms [KF, KOECHER, MILNE, RADEMACHER, DIAMOND, SCHOENENBERG, APOSTOL, KOEHLER, ONO, MIYAKE, SERRE, MART1] up to an additional factor of

(1.4) |

We seek in particular modular forms which have a representation in terms of Lambert–Eisenstein series [LAMBERT, EISENSTEIN] and can thus be represented by elliptic polylogarithms [ELLPOL]. However, the -factor (1.4) in general remains. Thus the occurring -integrands are modulated by a -factorial [SLATER, ERDELYI1] denominator.

Structures of the kind for are frequent even in the early literature. A prominent case is given by the invariant , see e.g. [HURWITZ],

(1.5) |

with an Eisenstein series, cf. Eq.(LABEL:eq:EISEN1), and the discriminant

(1.6) |

In the more special case considered in [BLOCH2, Adams:2014vja, Adams:2015gva, Adams:2015ydq, Adams:2016xah] terms of this kind are not present.

For the present solutions, we develop the formalism in Section 5. We discuss possible extensions of integral classes to the present case in Section 6 and of elliptic polylogarithms [ELLPOL], as has been done previously in the calculation of the two-loop sunrise and kite-diagrams [BLOCH2, Adams:2014vja, Adams:2015gva, Adams:2015ydq, Adams:2016xah]. Here the usual variable is mapped to the nome , expressing all contributing functions in the new variable. This can be done for all the individual building blocks, the product of which forms the desired solution. Section 7 contains the conclusions.

In Appendix LABEL:sec:8

we briefly describe the algorithm finding for second order ordinary differential equations

solutions with a rational function argument. In Appendix LABEL:sec:10 we present for convenience details for the necessary steps to arrive at the elliptic polylogarithmic representation in the examples of the sunrise and kite integrals [BLOCH2, Adams:2014vja, Adams:2015gva, Adams:2015ydq, Adams:2016xah]. Here we compare some results given in Refs. [BLOCH2] and [Adams:2014vja]. In Appendix LABEL:sec:11 we list a series of new sums, which simplify the recent results on the sunrise diagram of Ref. [Ananthanarayan:2016pos].In the present paper we present the results together with all necessary technical details and we try to refer to the related mathematical literature as widely as possible, to allow a wide community of readers to apply the methods presented here to other problems.

## 2 The Differential Equations

The master integrals considered in this paper satisfy linear differential equations of second order

(2.1) |

with rational functions , which may be decomposed into

(2.2) |

The homogeneous equation is solved by the functions , which are linearly independent, i.e. their Wronskian obeys

(2.3) |

The homogeneous Eq. (2.1) determines the well-known differential equation for

(2.4) |

which, by virtue of (2.2), has the solution

(2.5) |

normalizing the functions accordingly. A particular solution of the inhomogeneous equation (2.1) is then obtained by Euler-Lagrange variation of constants [EULLAG]

(2.6) |

with

(2.7) |

and two constants to be determined by special physical requirements. We will consider indefinite integrals for the solution (2.6), which allows for more singular integrands. For the class of differential equations under consideration, can be expressed by harmonic polylogarithms and rational functions, is a polynomial, and the functions turn out to be higher transcendental functions, which are even expressible by complete elliptic integrals in the cases considered here. Therefore Eq. (2.6) constitutes a nested integral of known functions [Remiddi:1999ew, Moch:2001zr, Ablinger:2013cf, Ablinger:2011te, Ablinger:2014bra] and elliptic integrals at rational argument.

We consider the systems of differential equations [Grigo:2012ji] for the terms of the master integrals

(2.8) |

and

(2.9) |

with

(2.10) | |||||

(2.11) | |||||

(2.12) | |||||

(2.13) | |||||

By applying decoupling algorithms [Zuercher:94, ORESYS, NewUncouplingMethod] one obtains the following scalar differential equation

(2.16) | |||||

(2.17) | |||||

The above differential equations of second order contain more than three singularities. We seek solutions in terms of Gauß’ hypergeometric functions with rational arguments, following the algorithm described in Appendix LABEL:sec:8. It turns out that these differential equations have solutions.

Two more master integrals are obtained as integrals over the previous solutions. They obey the differential equations

(2.18) | |||||

(2.19) | |||||

with the values of the Riemann -function at integer argument and the harmonic polylogarithms are defined by [Remiddi:1999ew]

(2.20) | |||||

Subsequently, we will use the shorthand notation . The harmonic polylogarithms occurring in the inhomogeneities of Eqs. (2.18, 2.19) can be rewritten as polynomials of

(2.21) |

cf. [Blumlein:2003gb].

## 3 Solution of the homogeneous equation

In the following we will derive the solution of the homogeneous part of Eqs. (2, 2.16) as examples in detail, using the algorithm outlined in Ref. [IVH], see also Appendix LABEL:sec:8.

The homogeneous solutions of Eq. (2) read

(3.1) | |||||

(3.2) |

with

(3.3) |

The solutions (3.1, 3.2) are close integer series [VANH1] obeying

(3.4) |

with . The Wronskian for this system is

(3.5) |

The solutions are shown in Figure 1.

Equivalent solutions are found by applying relations due to triangle groups [TAKEUCHI], see Appendix LABEL:sec:8,

(3.6) | |||||

(3.7) | |||||

where

(3.8) |

These solutions have the Wronskian (3.5) up to a sign^{5}^{5}5The sign can be adjusted by . but differ from those in (3.1, 3.2).
The ratios of the homogeneous solutions are given by

(3.9) | |||||

(3.10) |

The hypergeometric functions appearing in (3.6, 3.7) are given in terms of complete elliptic integrals [TRICOMI]

(3.11) | |||||

(3.12) |

Both solutions obey (3.4) with .

We also used the relation [HYP1]

(3.13) |

noting that it is always possible to map a function with into complete elliptic integrals. Their integral representations in Legendre’s normal form [LEGENDRE] read

(3.14) | |||||

(3.15) |

In going from (3.1, 3.2) to (3.6, 3.7) also a contiguous relation had to be applied, leading to a linear combination of two hypergeometric functions. The solutions are shown in Figure 2.

The ratio exhibits the interesting form

(3.16) |

with

(3.17) |

Whether solutions emerging in single scale Feynman integral calculations as solutions of differential equations for master integrals are always of the class to be reducible to complete elliptic integrals a priori is not known. However, one may use the algorithm given in Appendix LABEL:sec:8 to map a solution to one represented by elliptic integrals, if the parameters of the respective solution match the required pattern.

The homogeneous solutions of (2.16) read

(3.18) | |||||

(3.19) |

with

(3.20) |

The argument appeared already in complete elliptic integrals by A. Sabry in Ref. [SABRY], Eq. (68), with , calculating the so-called kite-integral at 2 loops, 55 years ago; see also Ref. [Laporta:2004rb], Eq. (A.11), for the sunrise-diagram and [Remiddi:2016gno], Eq. (D.18), with for the kite-diagram. The latter aspect also shows the close relation between the elliptic structures appearing for both topologies, which has been mentioned in Ref. [Adams:2016xah].

Using the Legendre identity [LEGENDRE]

(3.21) |

one obtains the Wronskian of the system (3.18, 3.19)

(3.22) |

cf. (2.5).

The homogeneous solutions (3.18, 3.19), which are complex for , are shown in Figure 3. The real part of has a discontinuity at moving from to , while Re vanishes for . Likewise, Im vanishes for .

We finally consider the ratio

(3.23) |

where

(3.24) |

This structure is the same as in (3.16, 3.17) up to the pre-factor.

The solution of the inhomogeneous equations (2, 2.16) are obtained form (2.6) specifying the constants by physical requirements. The previous calculation of the corresponding master integrals in [Grigo:2012ji] used expansions of the propagators [Q2EXP, Steinhauser:2000ry], obtaining series representations around and . The first expansion coefficients of these will be used to determine and . The inhomogeneous solutions are given by

(3.25) |

with

(3.26) |

Eq. (3.25) is an integral which cannot be represented within the class of iterative integrals. It therefore requires a generalization. We present this in Section 4. Efficient numerical representations using series expansions are given in Section 5.

## 4 Iterated Integrals over Definite Integrals

The elliptic integrals (3.14, 3.15) cannot be rewritten as integrals in which their argument
only appears in one of their integral boundaries.^{6}^{6}6Iterative non-iterative integrals have been
introduced by the 2nd author in a talk on the 5th International Congress on Mathematical Software, held at
FU Berlin, July 11-14, 2016, with a series of colleagues present, cf. [ICMS16].
Therefore, the integrals of the type of Eq. (3.25)
do not belong to the iterative integrals of the type given in
Refs. [Remiddi:1999ew, Ablinger:2011te, Ablinger:2013cf, Ablinger:2014bra] and generalizations thereof to
general alphabets, which have the form

(4.1) |

For a given difference equation, associated to a corresponding differential equation, the algorithms of [SIG1, SIG2] based on [SIGBASE] allow to decide whether or not the recurrence is first order factorizable. In the first case the corresponding nested sum-product structure is returned. In the case the problem is not first order factorizable, integrals will be introduced whose integrands depend on variables that cannot be moved to the integration boundaries and over which one will integrate by later integrals. This is the case if the corresponding quantity obeys a differential equation of order , not being reducible to lower orders. Examples of this kind are irreducible Gauß’ functions, to which also the complete elliptic integrals and belong.

The new iterative integrals are given by

(4.2) | |||||

and cases in which more than one definite integral appears. Here the are the usual letters of the different classes considered in [Remiddi:1999ew, Ablinger:2011te, Ablinger:2013cf, Ablinger:2014bra] multiplied by hyperexponential pre-factors

(4.3) |

and is given by

(4.4) |

such that the -dependence cannot be transformed into one of the integration boundaries completely. We have chosen here as a rational function because of concrete examples in this paper, which, however, is not necessary. Specifically we have

(4.5) | |||||

The new iterated integral (4.2) is not limited to the emergence of the functions (4.5). Multiple definite integrals are allowed as well. They emerge e.g. in the case of Appell-functions [HYP, SLATER] and even more involved higher functions. These integrals also obey relations of the shuffle type w.r.t. their letters , cf. e.g. [SHUF, Blumlein:2003gb].

Within the analyticity region of the problem one may derive series expansions of the corresponding solutions around special values, e.g. and other values to map out the function for its whole argument range. In many cases, one will even find convergent, widely overlapping representations, which are highly accurate and provide a numerical solution in terms of a finite number of analytic expansion coefficients. We apply this method to the solution of the differential equations in Section 2 in the following section and return to the construction of a closed form analytic representation using -series and Dedekind functions in Section 6.

## 5 The Solution of the Inhomogeneous Equation by Series Expansion

The inhomogeneous solutions of type (3.25) can be expanded into series around
and analytically using computer algebra packages like Mathematica or maple.
One either obtains Taylor series or superpositions of Taylor series times a factor . For all
solutions both expansions have a wide overlap^{7}^{7}7This technique has also been used in Ref. [TLSR2a].
and one may obtain in this way a highly accurate representation of all solutions in the complete region .

In the following we present the first terms of the series expansion for the functions around and .

For we obtain

(5.1) | |||||