Three-loop contributions to the $\rho$ parameter and iterated integrals of modular forms

We compute fully analytic results for the three-loop diagrams involving two different massive quark flavours contributing to the $\rho$ parameter in the Standard Model. We find that the results involve exactly the same class of functions that appears in the well-known sunrise and banana graphs, namely elliptic polylogarithms and iterated integrals of modular forms. Using recent developments in the understanding of these functions, we analytically continue all the iterated integrals of modular forms to all regions of the parameter space, and in each region we obtain manifestly real and fast-converging series expansions for these functions.


Introduction
Multiloop Feynman integrals are a cornerstone of perturbative Quantum Field Theory, as they are the main mathematical building blocks for the computation of higher-order corrections to physical observables. For this reason, a lot of effort has been put in understanding their mathematical structure and the class of special functions to which they evaluate. For instance, unitarity dictates that Feynman integrals must be multivalued functions with logarithmic branch cuts, and this must be reflected in the corresponding class of functions.
A lot of progress was made over the last decade in understanding the simplest class of special functions that appear in multiloop computations. It was realised that many Feynman integrals can be evaluated in terms of multiple polylogarithms [1][2][3]. These functions are by now well understood, including their analytic continuation to arbitrary kinematic regions and their efficient numerical evaluation for arbitrary complex arguments [4][5][6][7].
While multiple polylogarithms are sufficient to express all one-loop integrals in four space-time dimensions, it has been known for a long time that new classes of functions can appear starting from two-loop order. This was noted for the first time in the calculation of the two-loop corrections to the electron propagator in QED with massive electrons [8]. The simplest two-loop integral which cannot be expressed in terms of multiple polylogarithms is the so-called two-loop sunrise graph with three massive propagators, which has been extensively studied over the last decades [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27]. Recent interest in the sunrise graph was fuelled to a large extend by the fact that the same type of functions appears in many higherorder calculations for the Large Hadron Collider (LHC) at CERN where the mass of the top quark cannot be neglected, see for example refs [28][29][30][31][32][33][34][35][36]. By now, several representations of the sunrise graph in terms of different classes of special functions are known. They all have in common that they involve functions related to elliptic curves, most prominently elliptic multiple polylogarithms and iterated integrals of modular forms.
Since elliptic functions seem to be a feature of many two-loop diagrams with massive propagators, it is natural to expect that these functions also prominently show up when performing calculations in the electroweak sector of the Standard Model (SM), where the gauge bosons and the fermions acquire a mass through the Higgs mechanism. In this paper we consider one of the precision observables of the electroweak SM, namely the corrections to the ρ parameter, defined as the difference between the vacuum polarisations of the W and Z bosons. The ρ parameter is known through three-loop order in the SM in the limit where all quarks but the top quark are massless [37][38][39][40][41][42][43]. In ref. [44], corrections from three-loop diagrams were considered in the scenario where also the bottom quark has a non-vanishing mass. The results of ref. [44] where presented as an expansion in the ratio of the two quark masses. However, no closed analytic formula was presented, because the corresponding loop integrals were observed to involve functions of elliptic type, and at the time the theory of these functions was still largely underdeveloped.
While the series expansion of ref. [44] is sufficient to obtain reliable phenomenological predictions, it is interesting to revisit the computation of ref. [44] in the light of the recent developments in the understanding of elliptic Feynman integrals. In this way one can obtain fully analytic results at high-loop order for one of the most fundamental precision observables of the SM. First steps in this direction were taken in refs. [45,46], where the elliptic Feynman integrals of ref. [44] were computed in terms of a class of new transcendental functions, called iterative non-iterative integrals. The relationship between these functions and the functions that appear in the different known representations of the sunrise graph is not fully clear. It is an interesting question if completely new classes of elliptic-type functions show up in the computation of the ρ parameter which are not covered by the existing literature on the sunrise graph and on elliptic polylogarithms and iterated integrals of modular forms.
The purpose of this paper is to show that all the Feynman integrals that contribute to the three-loop corrections to the ρ parameter with two massive quark flavours can be expressed in terms of exactly the same class of functions that appear in the sunrise graph.
We show how all the elliptic Feynman integrals of ref. [44] can be performed in terms of iterated integrals of modular forms for the same congruence subgroup as for the sunrise integral. These iterated integrals can be analytically continued to the whole parameter space in such a way that they admit fast-converging series representations for all values of the quark masses.
Our paper is organised as follows: In Section 2 we review the results of ref. [44] and we identify the elliptic Feynman integrals that need to be computed. In Section 3 we introduce the mathematical background on elliptic curves and elliptic polylogarithms needed throughout the paper. In Section 4 we present our first main result, and we show how to evaluate the Feynman parameter representation for the simplest elliptic Feynman integral in terms of elliptic polylogarithms. In Section 5 we review the connection between elliptic polylogarithms and iterated integrals of modular forms, and in Section 6 we use this relationship to obtain analytic results for all elliptic Feynman integrals that contribute to the three-loop ρ parameter in the region where the ratio of the quark masses is small. In Section 7 we discuss the analytic continuation of these integrals to the whole parameter space. In Section 8 we present our final analytic result for the ρ parameter at three loops, and in Section 9 we draw our conclusions. We include several appendices where we collect formulas omitted throughout the main text.

Notations and computational setting
In this section we review the background to the three-loop QCD corrections to the electroweak ρ parameter with two massive quark flavours. We closely follow the presentation of ref. [44]. The ρ parameter can be written as: where the higher-order corrections are given by Σ Z (0) and Σ W (0) are, respectively, the transverse parts of the Z and W boson propagators at zero momentum. They are defined as where Π µν Z/W are the correlator functions for the Z and W bosons, and d is the space-time dimension.
We will consider the three-loop QCD corrections to the ρ parameter in n f -flavour QCD with n f − 2 massless quarks and two massive ones, whose masses we denote by m 1 and m 2 . In the SM, this corresponds to n f = 6, and m 1 and m 2 denote the masses of the top and bottom quarks respectively. It is possible to write the higher-order corrections to the ρ parameter as an expansion in the strong coupling constant α s : where G F is the Fermi constant and m t is the top-quark MS-mass. We will choose the renormalisation scale as µ 2 = m 2 t . The term δ (j) involves graphs with j + 1 loops, and this paper is devoted to the calculation of certain three-loop masters appearing in δ (2) .
Using Integration-By-Parts (IBP) identities [47,48], it is possible to reduce all threeloop integrals that contribute to δ (2) to the computation of a small set of master integrals [44]. All but six of the master integrals were evaluated analytically in terms of multiple polylogarithms. The remaining integrals were shown to involve elliptic functions. They can be embedded in two six-propagator topologies which we call A and B, and which only differ by exchanging the roles of m 1 and m 2 (see fig. 1). Topology A can be written as where the propagators are defined as (2.6) We work in dimensional regularisation with d = 4 − 2 , and choose the measure Topology A has ten master integrals, whose graphs we list in fig. 1. It is convenient to introduce a variable t defined as 1 For concreteness, we will choose m 2 2 < m 2 1 , i.e., 0 < t < 1. The case t > 1 can be obtained by analytic continuation, and we will discuss this at the end of the paper. Note that integrals from Topology B can be identified with integrals from Topology A with t > 1. Hence, it is sufficient to compute Topology A and to understand the analytic continuation to all positive values of t. We therefore only focus on Topology A with 0 < t < 1 in the first sections of this paper, and we only comment on Topology B when we discuss the analytic continuation of Topology A to t > 1. For the physical ρ parameter in the SM, we will be interested in topologies A and B evaluated at t phys = m 2 b /m 2 t , with It is convenient to work with the following basis of master integrals: (2.10) We have normalised all master integrals to be dimensionless, i.e., the functions f i only depend on t.
An efficient way to compute the master integrals is to use differential equations [3,[49][50][51][52]. The master integrals f 1 through f 7 satisfy a system of differential equations in so-called canonical form [53]:

11)
A 0 and A 1 are matrices of integer numbers, which we give explicitly in Appendix A. The canonical form of the differential equation in eq. (2.11) makes it manifest that the functions (f 1 , . . . , f 7 ) can be expressed in terms of a well-studied class of special functions called multiple polylogarithms (MPLs) [1,54]: G(a 1 , · · · , a n ; x) = x 0 du u − a 1 G(a 2 , · · · , a n ; u) , G(; x) = 1, a n = 0 . (2.12) In the case where a n = 0, the naive recursive definition in eq. (2.12) is divergent, and we define instead, Since eq. (2.11) only has singularities at t ∈ {0, 1, ∞}, we only need to consider the case where a i ∈ {0, 1} in eq. (2.12), which defines a class of functions known as harmonic polylogarithms [2]. Solving eq. (2.11) in terms of MPLs is standard, and the analytic results can be found in ref. [44]. We will therefore not discuss these integrals any further, and we only mention that we find complete agreement with the results of ref. [44].
Let us now turn to the remaining master integrals f 8 , f 9 and f 10 . It was already observed in ref. [44] that these integrals are not expressible in terms of MPLs. This can be seen from the perspective of differential equations. Before we discuss the form of the differential equations, it is convenient to change the basis of master integrals. More precisely, we consider the integrals f 8 and f 9 in d = 2 − 2 dimensions where they are finite. We denote the corresponding master integrals in d = 2 dimensions by f (where we set = 0) to obtain their four-dimensional analogues, which leads to simplifications in our computations.
Let us now have a closer look at the differential equations satisfied by f 8 and f (2) 9 . The differential equations remain coupled in the limit = 0. The corresponding homogeneous 2 × 2 system has the form 9,h . (2.14) It is convenient to transform this system into a second-order equation satisfied by f 8,h , with Remarkably, this differential operator is the same one as in the case of the sunrise integral with a massive external leg and three massive propagators of equal mass [13]. As a consequence, any element in the kernel of D 2 t can be written as a linear combination of the following two functions [13]: where λ(t) is given by 18) and K is the complete elliptic integral of the first kind, The functions in eq. (2.17) are respectively real and purely imaginary for 0 < t < 1. They can be analytically continued outside of this range [24,57]. We will return to the analytic continuation in Section 7. Equation (2.17) is sufficient to construct the general solution to the homogeneous system in eq. (2.14). The Wronskian matrix of the system is (2.20) The general solution to the homogeneous equation then takes the form While the solution of the homogeneous system in eq. (2.14) is well known, solving the inhomogeneous system satisfied by f and f 10 is not trivial. Here, we will compute f  is determined by the derivative of f  by differentiating f (2) 8 . Finally, we can compute f 10 from the differential equation that it satisfies, which reads where Li n (t) denote the classical polylogarithms, We see that the homogeneous part of the differential equation in eq. (2.22) is in canonical form. Nevertheless, f 10 cannot be expressed in terms of MPLs only, because it contains f in the inhomogeneous term. We can thus determine f 10 once f is known analytically. Note that f (2) 9 does not contribute to the differential equation for f 10 through O( 4 ), which is the order through which f 10 enters the three-loop corrections to the ρ parameter.
Before we discuss the computation of f and f 10 in the next section, let us mention that the system of differential equations satisfied by f (2) 8 and f (2) 9 was also analysed in refs. [45,46]. There the corresponding second-order differential operator was not directly in a form which matches the one for the sunrise graph in eq. (2.16). As a consequence, the homogeneous solutions of refs. [45,46] (i.e., the equivalent of our f (2) 8,h and f (2) 9,h ) take a different form, which cannot easily be recognised as being related to the sunrise graph.
3 eMPLs on the torus and on the elliptic curve As mentioned in the previous section, our strategy for the computation of f and f 10 relies on obtaining an analytic solution for f (2) 8 from its Feynman parametric representation. Before we discuss this in detail, we review in this section some of the mathematical background needed to perform all the integrals analytically. We keep the review to a strict minimum, and we refer to the literature for a more detailed discussion (see, e.g., refs. [58][59][60] and references therein).

Elliptic curves and elliptic functions
Let P n (x) be a polynomial of degree n. An elliptic curve can be defined (loosely) by the equation y 2 = P n (x) for n = 3, 4. In this paper we are naturally led to an elliptic curve defined by a cubic polynomial, so we will focus our discussion on the polynomial equation We call a 1 , a 2 and a 3 the branch points of the elliptic curve. Seen as points in CP 2 and using homogeneous coordinates, the elliptic curve is given by the points [x, y, 1] that satisfy eq. (3.1), together with the point [0, 1, 0] on the infinity line. It is important to establish our conventions for the branches of the square root. Throughout this paper we follow ref. [58]: if the branch points are real and ordered as a 1 < a 2 < a 3 , then Some ubiquitous quantities that appear in the study of elliptic curves are the periods ω i , with K(λ) defined in eq. (2.19), and the quasi-periods η i , where (3.5) E(λ) denotes the complete elliptic integral of the second kind, The periods and quasi-periods are related by the Legendre relation: Every elliptic curve is isomorphic to a complex torus. More precisely, if we define then the points in CP 2 that satisfy eq. (3.1) are isomorphic to the quotient C/Λ τ where the two-dimensional lattice Λ τ is defined as Note that τ can always be chosen to lie in the complex upper half-plane H = {τ ∈ C : Im τ > 0}. We can construct a map from the torus to the elliptic curve with a function is not relevant for this paper and we refer the reader to ref. [58] for a more explicit definition. We can also define a map which assigns to a point [x 0 , y 0 , 1] on the elliptic curve defined by

Elliptic multiple polylogarithms
A natural class of functions to consider when working with Feynman integrals are elliptic multiple polylogarithms (eMPLs). Loosely speaking, eMPLs can be thought of as a class of iterated integrals that generalises the complete elliptic integrals of eqs. (2.19) and (3.6), in the same way that the MPLs in eq. (2.12) generalise the logarithm.
Since we have two ways of describing an elliptic curve-as a torus C/Λ τ or as a set of points in CP 2 -there are also two equivalent ways of defining eMPLs. We start by defining them in terms of iterated integrals along a path on the torus as [61,62] The integer k is called the length, n i ∈ N * and i n i is the weight. The integration kernels g (n i ) (z, τ ) are defined through the Eisenstein-Kronecker series, where θ 1 is the Jacobi theta function and θ 1 is its derivative with respect to the first argument.
Unlike what happens for the MPLs in eq. (2.12), there is an infinite number of kernels one must consider for eMPLs. While all these kernels are necessary for the study of this class of functions, in any particular calculation of a Feynman integral only a limited set will be relevant. As we will see below, only g (1) (z, τ ) is relevant in our case. From now on we focus on this kernel and refer the reader to ref. [61] for a more complete exposition. First, we note that g (1) (z, τ ) has a simple pole with unit residue at each point of Λ τ . Iterated integrals over this kernel will thus only have logarithmic singularities. In particular, some of the integrals in eq. (3.12) may be divergent and require regularisation. We refer the reader to ref. [62] for details on how this procedure can be implemented. Second, is not periodic under translations in both directions of the lattice Λ τ : Finally, it is possible to write a closed form for the total differential of an eMPL. The total differential forΓ A; z, τ , where , is given by [60] dΓ . The hat means that the corresponding argumentÂ p is removed, and we use the conventions (z 0 , z k+1 ) = (z, 0), (n 0 , n k+1 ) = (0, 0). The ω (n) ij are differential one-forms given by Since we can equivalently represent an elliptic curve as the zero set of some polynomial equation in CP 2 , see eq. (3.1), there is an alternative definition of eMPLs which uses directly the coordinates (x, y) instead of the coordinate z on the torus. This representation is important in the context of our current discussion because, in applications to Feynman integral calculations, elliptic curves often arise via explicit polynomial equations. The change of variables from z to x is given by eq. (3.10). Inserting it into the definition of eMPLs in eq. (3.12) we arrive at the following class of iterated integrals [58,59] where the integration kernels are related to the kernels in eq. (3.12) via Here, z x is the image of [x, y, 1] on the torus, see eq. (3.11). The length and the weight of the eMPL in eq. (3.17) are respectively defined as k and k i=1 |n i |. Restricting ourselves to the kernels that will be relevant in this paper, we have We note that some of the integrals in eq. (3.17) may diverge and require regularisation. We refer to ref. [59] for a detailed discussion. Let us make some comments about the iterated integrals defined in this section. First, the functions defined in eqs. (3.12) and (3.17) satisfy the usual properties of iterated integrals. In particular, they form a shuffle algebra, i.e., any product of these functions evaluated at the same value of the upper integration limit can be written as a linear combination of the same class of functions. Second, using eq. (3.19) we can easily translate between the eMPLs defined in eqs. (3.12) and (3.17). Every linear combination of iterated integralsΓ can be written as a linear combination of E 3 functions, and vice-versa. Finally, eMPLs contain MPLs as a subspace. Indeed, the kernel ψ 1 (c, x, a) is precisely the kernel that appears in the definition of MPLs (c.f. eq. (2.12)) and thus in terms of eMPLs directly from its Feynman parameter representation. The basic idea is to perform all integrations in terms of MPLs, except for the last one which can be performed in terms of eMPLs, cf. refs. [63,64].
The Feynman parameter representation for f where Σ can be any non-empty subset of {1, 2, 3, 4}, and We recall that we can put = 0 because the integral is finite in d = 2 dimensions.
We find it convenient to choose Σ = {2}, which amounts to setting x 2 = 1 in F. The integrals over x 1 and x 3 can then be easily performed in terms of MPLs. After changing variables tox 4 where we used a compact notation which we now explain. First, we define y as (4.5) We note that in the region 0 < t < 1 the a i are real and a 1 < a 2 < a 3 . Then we define χ ± as Finally, we use the shorthand While we were able to integrate over three Feynman parameters without leaving the space of MPLs to reach eq. (4.3), the square-root y appearing in thex 4 integration means that the last integral will leave this space and should instead be carried in terms of eMPLs. The elliptic curve is defined by eq. (4.4), and our first step is to recast the integrand of eq. (4.3) in terms of eMPLs. We write the integrand of eq. (4.3) as The MPLs appearing in this expression can be written as eMPLs using eq. (3.21). Furthermore, since the kernels in eq. (4.9) are all of the form of the kernels in eq. (3.19), we can write Ω(x; t) in terms of eMPLs. This requires determining the boundary contribution of Ω(x; t) at x = 0, which is given by Starting from this representation, we can then compute the integral over the last Feynman parameterx 4 in terms of eMPLs simply using eq. (3.17), and we find where the overall factor Ψ 1 (t) is defined in eq. (2.17). Equation (4.11) is one of the main results of this paper and expresses the integral f 8 in terms of eMPLs. Let us make some comments about this result. First, we observe that f (2) 8 is proportional to the homogeneous solution Ψ 1 (t), which is a period of the elliptic curve defined by the polynomial equation in eq. (4.4). The period is multiplied by a linear combination of eMPLs of uniform weight two. If we assign weight one to the period Ψ 1 (t) [59], then f (2) 8 has uniform weight three. This is consistent with the fact that f can be interpreted as a banana graph in two dimensions with two distinct masses evaluated at zero external momentum. The banana graph in two dimensions is known analytically in the case of three equal masses, and it indeed evaluates to a function of uniform weight three [65]. Second, we can also express f (2) 8 in terms of the eMPLsΓ. We start by noting that the eMPLs in eq. (4.11) are evaluated at x ∈ {0, 1, t t−1 }. Under eq. (3.11) these points are mapped to where τ denotes the ratio of the two periods of the elliptic curve (see eq. (2.17)), Since there is never any ambiguity, we will in the following not explicitly write the dependence of the z i and τ on t. We can then use this to rewrite the integration kernels (3.19) in terms of those found in the definition of theΓ, and we find (4.14) Using this change of variables in eq. (3.17), we can express all the E 3 functions in eq. (4.11) in terms ofΓ functions. The procedure is straightforward, but the result is lengthy and not particularly illuminating, so we do not show it here explicitly. We only mention that since all points in eq. (4.12) are rational points of the form z i = r 6 + s τ 6 , with r and s integers, theΓ functions will have a very special form, namely they will all have arguments that are rational points.
Our next goal is to compute the master integrals f 9 and f 10 . As explained at the end of Section 2, this can be done by differentiating or integrating f (2) 8 with respect to t. These operations, however, are not straightforward to carry out on the expressions in eq. (4.11), because the eMPLs depend on t in a highly non-trivial way. It would be desirably to have a representation of f (2) 8 in terms of iterated integrals with a simple dependence on the kinematic variable. While in general such a form may not be easily obtained, the special rational form of the points in eq. (4.12) allows one to find such a representation in this case. This will be reviewed in the next section.

eMPLs and iterated Eisenstein integrals
In this section we review how eMPLs evaluated at rational points can be expressed in terms of another class of iterated integrals, namely the so-called iterated Eisenstein integrals. Eisenstein series are a special case of modular forms. In the first part of this section we review modular forms and Eisenstein series, and in a second part we review the relationship between iterated Eisenstein integrals and eMPLs.

Modular forms: a brief introduction
In Section 3.2 we have seen that that every elliptic curve is isomorphic to a torus C/Λ τ , for some value τ ∈ H. Different values of τ , however, do not necessarily describe different tori. Indeed, we can replace the basis of periods (ω 1 , ω 2 ) by an integer linear combination of them without changing the lattice they generate. More precisely, let τ ∈ H be obtained from τ via a modular transformation, defined as where Then τ and τ define the same lattice, Λ τ = Λ τ , and so they also define the same elliptic curve. In other words, τ and τ define the same elliptic curve if and only if they are related by a modular transformation. For this reason, modular transformations play a central role in the study of elliptic curves.
In applications it is often specific subgroups of SL(2, Z) that are of interest. Particularly important subgroups are the so-called congruence subgroups of level N , defined by . The equivalence class that contains i∞ is called the cusp at infinity. As an example, for γ ∈ Γ(1) = SL(2, Z) we find that γ · (i∞) = a/c and thus there is a single cusp for Γ(1), the cusp at infinity. For any N , the number of cusps of the congruence subgroups of eq. (5.3) is finite.
Our goal is to construct functions that transform nicely under some congruence subgroup Γ. A modular function for Γ is a meromorphic function f : H → C that is invariant under Γ. One can show that every modular function has at least one pole. If we want to consider holomorphic functions, i.e., functions without poles, then we need to consider more general transformations. A modular form of weight n for Γ is a function f : H → C that is holomorphic on H and at the cusps of Γ such that Here S n (Γ) denotes the space of cusp forms of weight n, i.e., the modular forms of weight n that vanish at all cusps of Γ. Its complement E n (Γ) is the space of Eisenstein series. The Eisenstein series are of particular interest when working with eMPLs, and so we discuss them in detail in the remainder of this section.
Of special importance here will be the functions [66] (0 ≤ r, s < N ), One can show that these functions are always Eisenstein series of weight n for Γ(N ) [60,66]. Moreover, they form a spanning set of E n (Γ(N )), i.e., every element of E n (Γ(N )) can be written as a linear combination of the functions in eq. (5.7). Note that they do not form a basis. The relations between them are however understood [60,66]. We also mention that these functions are real whenever τ is purely imaginary [66]. Finally, the Eisenstein series in eq. (5.7) are closely related to the coefficients g (n) (z, τ ) in eq. (3.13). More precisely, we have the relation (0 ≤ r, s < N ) [60,66] The previous equation shows that there is a connection between Eisenstein series for Γ(N ) and eMPLs evaluated at rational points of the form z = r N + s N τ . We know from eq. (4.12) that the arguments of the eMPLs that appear in f  to be closely connected to Eisenstein series. We review this connection in the next section.

Iterated Eisenstein integrals
Consider a set of modular forms f j (τ ) of weight n j for some congruence subgroup of level N . We define their iterated integral as [67,68] with I(; τ ) ≡ 1. A precise definition of these integrals requires a careful regularisation of divergences that can appear at the cusp at infinity τ = i∞, and we refer to ref. [68] for a detailed discussion. We define the length of I(f 1 , . . . , f k ; τ ) to be k and the weight is −k + k j=1 n j . In the case where all the modular forms f j (τ ) are Eisenstein series, we refer to the integral in eq. (5.9) as an iterated Eisenstein integral.
In ref. [60] it was shown that whenever an eMPL of weight n is evaluated at rational points of the form z = r N + s N τ , then this eMPL can be expressed as a linear combination of uniform weight n of iterated Eisenstein integrals of level N . More precisely, we see from eq. (3.15) that the total differential of a length-k eMPL is given by eMPLs of length k − 1 and one-forms ω ij . If all the arguments of the eMPL are rational, then we can write ω (n) ij in terms of Eisenstein series using eq. (5.8). Hence, if we assume recursively that the claim is true for eMPLs up to length k − 1, we see that the total differential of an eMPL of length k evaluated at rational points only involves iterated Eisenstein integrals and Eisenstein series. We can then integrate back in τ to obtain the desired representation at length k.
This recursive consideration also provides an efficient algorithm to express eMPLs evaluated at rational points in terms of iterated Eisenstein integrals. The starting point is to write an eMPL evaluated at rational points as the integral of its derivative with respect to τ ,Γ (5.10) Since the differential lowers the length by one, we can recursively express the integrand in terms of iterated Eisenstein integrals and integrate back using eq. (5.9). The integration constant is obtained by studying the behaviour at the cusp at infinity. For the cases of interest here, it is possible to compute Cusp Γ (A 1 · · · A k ; z, τ ) by performing a series expansion of the integrands inΓ (A 1 · · · A k ; z, τ ), and integrating only the leading terms (see the example below). This algorithm is iterative in the length of the eMPLs, and the starting point is the total differential of an eMPL of length one, where the differential one-forms ω (n 1 −r) ij can be expressed in terms of Eisenstein series via eq. (5.8) whenever z 1 and z are rational points.
We finish this section by noting that the procedure that we just presented to relate eMPLs and iterated Eisenstein integrals can be reformulated in terms of the coaction on eMPLs [60,69]. Indeed, the algorithm we presented is summarised by the relatioñ where ∆ Γ A; z, τ is the coaction on the eMPLs and we defined m[a ⊗ b] ≡ ab.
6 The master integrals for Topology A in the region 0 < t < 1 In this section we present our final result for the master integrals f 9 and f 10 in terms of iterated Eisenstein integrals. We focus in this section on the region 0 < t < 1, and we explore the analytic continuation to other regions in the next section.
We start by discussing f (2) 8 . We can follow the steps outlined in the previous section and express all the eMPLs that appear in the analytic expression for f (2) 8 in eq. (4.11) in terms of iterated Eisenstein integrals. We find 6,1,0 , a 6,1,0 ; τ + 16 I 1, a 6,1,0 , a 6,1,3 ; τ + 16 I 1, a 6,1,0 , a 6,2,0 ; τ 6,1,3 , a 6,1,0 ; τ + 16 I 1, a 6,1,3 , a 6,1,3 ; τ + 16 I 1, a 6,1,0 , a 6,2,0 ; τ 6,2,3 , a 6,1,0 ; τ − 40 I 1, a 6,2,3 , a 6,1,3 ; τ − 40 I 1, a 6,2,3 , a 6,2,0 ; τ +4 log 3 5 I 1, a 6,2,3 ; τ − 2 I 1, a 6,1,0 ; τ − 2 I 1, a where Cl 2 (x) denotes the Clausen function, The variable τ (t) is defined in eq. (4.13) and is we recalled here for convenience, Equation (6.4) can be inverted, and we find [17,70,71] Next, let us discuss f 9 . Using its differential equation, we obtain f (2) where Φ 1 (t) = ∂ t Ψ 1 (t) was defined in eq. (2.20) and J (t) denotes the Jacobian of the change of variables from t to τ , . (6.8) The derivative of f (2) 8,U with respect to τ can easily be carried out, as the iterated integrals in eq. (6.2) only depend on τ through the upper integration limit. We find f (2) 9,U (τ ) = 2πi ∂ τ f . (6.9) Finally, we discuss the calculation of f 10 . Aside from classical polylogarithms, the differential equation for f 10 in eq. (2.22) contains f (2) 8 as an inhomogeneous term. In order to solve the differential equation we follow the strategy of ref. [65]. We start by noting that we have expressed f (2) 8 in eq. (6.1) in terms of iterated Eisenstein integrals for the congruence subgroup Γ(6), for which a spanning set is given by the functions in eq. (5.7). It turns out that there is a smaller set of modular forms that is sufficient to express the result for f (2) 8 , namely Eisenstein series for Γ 1 (6). In ref. [72] it was shown that a basis for M n (Γ 1 (6)) is f n,p (τ ) = t(τ ) p Ψ 1 (t(τ )) n , 0 ≤ p ≤ n, f 0,0 (τ ) = 1 . (6.10) It is possible to write all polylogarithms that appear in eq. (2.22) in terms of iterated Eisenstein integrals for Γ 1 (6). Here we only sketch the argument, and we refer to refs. [65,71,72] for details. If we express the iterated integrals in eq. (6.2) in terms of the kernels in eq. (6.10), and we change variables from τ to t using eq. (6.4), we can write eq. (6.2) in terms of iterated integrals in t, with integration kernels of the form This class of iterated integrals contains at the same time MPLs (for n = 2) and iterated Eisenstein integrals for Γ 1 (6) (through eq. (6.10)). In other words, all contributions in eq. (2.22) can be expressed in terms of a unique class of iterated integrals, the iterated Eisenstein integrals for Γ 1 (6), and thus also in terms of Eisenstein series for Γ(6). Let us consider as an example G(0; t). Using eq. (6.8) it is possible to verify that the following identity holds where c is a boundary term associated with the regularisation of G(0; t). From the qexpansion of the iterated integrals on the right-hand side of the above equation [65], it is possible to write it in terms of iterated Eisenstein integrals : with c = log 9. After this step, the inhomogeneous term of eq. (2.22) only involves Eisenstein series and iterated Eisenstein integrals, and the differential equation for f 10 can easily be solved. We find f 10 (t) = 384 I a 6,2,0 ; τ + 384 I a 6,2,3 ; τ 6,1,0 ; τ + 2 I a 6,1,3 ; τ + I a 6,1,0 , 1; τ + 2 I a 6,1,3 , 1; τ + I a 6,2,3 , 1; τ (6.14) Equations (6.2), (6.9) and (6.14) are among the main results of this paper. They express all elliptic master integrals of Topology A in terms of a class of special functions that is well studied in both the mathematics and physics literature, namely eMPLs and iterated Eisenstein integrals for the congruence subgroup Γ(6). The discussion from the previous paragraph shows that we can even restrict the analysis to the larger congruence subgroup Γ 1 (6), at least through finite terms in the Laurent expansion in . This shows that all master integrals for Topology A can be expressed in terms of exactly the same class of functions as the well-known sunrise, kite and banana integrals with three or four equal masses. This extends the observation of Section 2 that f satisfy the same homogeneous differential equation as the sunrise graph. In particular, we find that there is no need to introduce new classes of transcendental functions beyond those already encountered for the sunrise and kite graphs. This is at variance with the analytic results for Topology A of refs. [45,46], where new classes of functions were introduced.
Since the integrals considered here and the sunrise graph seem to be so closely related, let us comment on the elliptic curves that appear in the computation of the two integrals. As already mentioned in Section 2, the second order differential operator in eq. (2.16) that describes the homogeneous solution is identical for the two sets of integrals. The solutions of eq. (2.16) are the two periods Ψ 1 (t) and Ψ 2 (t) of the family of elliptic curves parametrised by t. While in the case of the integrals considered here this curve is most naturaly defined by a cubic polynomial (cf. eq. (4.4)), the curve obtained from the Feynman parameter integral for the sunrise is defined by a quartic polynomial, cf. e.g. refs. [16,17,26]. There is no contradiction: the same elliptic curve may be represented as the zero set of different polynomial equations. An invariant that uniquely distinguishes different elliptic curves is its j-invariant. The j-invariant of the family of elliptic curves in eq. (4.4) is It is easy to check that eq. (6.15) agrees with the j-invariant for the family of elliptic curves obtained from the Feynman parametrisation of the equal-mass sunrise graph, with t = p 2 m 2 [17]. This shows that indeed the elliptic curves obtained from the Feynman parameter integrals of the sunrise integrals and the integrals considered here are identical. Finally, let us make a comment about the analytic structure of our results. We see that we can cast our results in the form, 9,U (t) and f 10,U (t) = f 10 (t) are defined in eqs. (6.2), (6.9) and (6.14). The matrix S is given by This form matches precisely the structure of elliptic Feynman integrals conjectured in ref. [59]. In particular, we see that the functions f 9,U (t) and f 10,U (t) are pure functions in the sense of ref. [59], and they have uniform transcendental weight two, three and four respectively.

Analytic continuation
As mentioned in Section 2, for the calculation of the ρ parameter we do not only need the integrals from Topology A (see fig. 1), but also those from Topology B, obtained by exchanging m 1 and m 2 . These integrals can equivalently be obtained by analytically continuing Topology A to the region t > 1. The analytic continuation of the non-elliptic integrals can be done using standard techniques. In this section we discuss the analytic continuation of the elliptic integrals f 9 and f 10 . The analytic continuation will be done following the steps described in ref. [66]. We start by discussing the analytic continuation of the homogeneous solutions in eq. (2.17). Given the singularities in the differential equation in eq. (2.16), there are four kinematic regions to consider: t < 0, 0 < t < 1, 1 < t < 0, t > 9. The homogeneous solution in eq. (2.17) is well behaved in the second region, by which we mean that they are local solutions to the differential equation that are respectively real and imaginary for t ∈ [0, 1]. The first step in the analytic continuation procedure is to obtain similarly well behaved solutions in the other three regions. Since the differential equation in eq. (2.16) is the same as that of the sunrise integral, we can simply reuse the results of refs. [24,57]: • Region 1, 0 < t < 1: • Region 2, 1 < t < 9: • Region 3, t > 9: • Region 4, t < 0: where λ(t) was defined in eq. (2.18) while λ 9 (t) and λ 0 (t) are given by: The functionsΨ (a,b) j are local solutions to the differential equation in eq. (2.16), but they do not extend individually to global solutions that define analytic functions with at most logarithmic singularities at the regular singular points of eq. (2.16). We can, however, construct a set of solutions with the desired analytic properties by patching together the local solutions in the right way. This leads to the correct analytic continuation of the functions in eq. (2.17) to all values of t, given by the piecewise definition Having obtained the correct analytic continuation of the homogeneous solutions to all values of t, we can also extend the definition of τ (t) in eq. (4.13) to arbitrary t, by replacing Ψ j (t) in eq. (4.13) by their piecewise definition in eq. (7.6). For the convenience of the reader, we reproduce here the definition of τ (t) as a function of the homogeneous solutions: With these definitions we can immediately evaluate the integrals f 9 and f 10 for all real values of t, for both Topology A and B. Indeed, we can insert the global definitions of Ψ j and τ from eqs. (7.6) and (7.7) into eq. (6.16). Each iterated Eisenstein integral admits a q-expansion (see eq. (5.5)), which is guaranteed to converge because Im τ (t) > 0 for all values of t. Therefore, we can obtain numerical results for all the integrals in eq. (6.16) for arbitrary t. The convergence of the q-expansion, however, might be very slow: it is controlled by the size of the imaginary part of τ which, depending on the value of t, might be very small. In the next section we address this shortcoming.

Numerical evaluation of iterated integrals of modular forms
In this section we discuss how we can speed up the numerical convergence of the q expansion of the iterated Eisenstein integrals. Indeed, for physical applications it is desirable to have representations for the master integrals that can be evaluated efficiently. This may however not be the case for our results, which are simply obtained by inserting the analytic continuation of the homogenous solutions into eqs. (7.7) and (6.16), without care for the convergence properties of the expression we obtain.
To make our point more concrete, let us consider the case where t takes its SM value t phys (cf. eq. (2.9)). For Topology A, we then find τ (t phys ) i 3.07. This corresponds to a value of q 6 = exp(πiτ (t)/3) 0.04 in the q-expansion of the iterated Eisenstein integrals (see eq. (5.5)) and so all integrals admit a fast converging q-expansion. For Topology B, however, the integrals are evaluated at τ (1/t phys ) 0.94 + i 0.13, which gives |q 6 | 0.88. The iterated Eisenstein integrals computed with τ (t phys ) will then converge one order of magnitude faster than those computed with τ (1/t phys ).
The convergence of the q-expansion can be substantially accelerated, for instance following the procedure of ref. [66] which we now briefly summarise. In a nutshell, the idea is to find a transformation γ ∈ SL(2, Z) that maximises the imaginary part of τ and thus the speed of the convergence of the q-expansion. For example, we can map τ to the so-called fundamental domain of SL(2, Z), defined as and |τ | > 1 {τ ∈ H : Re(τ ) ≤ 0 and |τ | = 1} .
(7.8) All τ ∈ F have Im(τ ) ≥ √ 3/2 and for every value of t, i.e., for every τ (t) ∈ H, there is a γ t ∈ SL(2, Z) such that γ −1 t · τ (t) ∈ F. While we might naively expect that the map to the fundamental domain in each of the four different regions for t in eq. (7.6) is given by a single transformation, or in other words In accordance with the previous discussion, τ A (t) admits a piecewise definition: to each of the eight regions bounded by the points t j in eq. (7.10) corresponds a different matrix γ t . We write γ t = γ (j,j+1) for t ∈ [t j , t j+1 ] , (7.13) where the indices are understood to be cyclically defined and (7.14) While we gave an explicit definition for τ A (t) for all values of t, doing the same for the different quantities in eq. (7.11) would lead us towards a very lengthy and repetitive discussion. Instead, we give these expression as ancillary files and focus our discussion in the remainder of this section on the two regions that are relevant for the calculation of the ρ parameter: t ∈ [t 2 , t 3 ] and t ∈ [t 7 , t 8 ].
We finish with a comment on Topology B. As already noted, it is obtained from Topology A by replacing t → 1/t. In eq. (7.11) we defined the vector of functions f A (t) with fast-converging representations of the master integrals for all real values of t. We adopt the same conventions for Topology B, i.e. all quantities with an index B admit a piecewise definition. They are simply obtained from those of Topology A using: Similarly, we have τ B (t) = τ A 1 t , and S B (t) = S A 1 t . (7.23) In particular, with these relations we can obtain Topology B at t phys ∈ [t 2 , t 3 ] from the expressions given above for Topology A for 1/t phys ∈ [t 7 , t 8 ].
As an example, f 6,0,1 , a 6,0,1 ; τ B (t phys ) − 40 I 1, b 6,0,1 , a 6,0,2 ; τ B (t phys ) 6,3,1 , a 6,0,1 ; τ B (t phys ) − 16 I 1, b 6,3,1 , a 6,0,2 ; τ B (t phys ) 6,3,2 , a 6,0,1 ; τ B (t phys ) + 16 I 1, b 6,3,2 , a 6,0,2 ; τ B (t phys ) − iπ 2 12 τ B (t) + ζ 3 π . (7.24) 8 Three-loop contributions to the ρ parameter We are now ready to give an expression for the three-loop contributions to the ρ parameter in terms of iterated integrals of modular forms. We start from the expression given in the ancillary files of ref. [44] renormalised in the MS scheme and set the regularisation scale µ 2 = m 2 t . The result depends on the SU (N c ) colour factors C A = N c and C F = (N 2 c − 1)/(2N c ), and on the number of massless quarks n l = n f − 2. 2 In their expressions, the authors of ref. [44] leave the elliptic master integrals associated with diagrams 8, 9 and 10 of fig. 1 (and their Topology B counterparts) unevaluated, making it particularly convenient to adapt their expression to our convention. For concreteness, we highlight the two main changes we make. First, we recall that we write our results in terms of t = m 2 2 /m 2 1 , whereas the expression of ref. [44] is written for x = m 2 /m 1 . Second, we found it more convenient to compute the masters associated with diagrams 8, 9 of fig. 1 in d = 2. Starting from their expression, we thus change variables to t and use the dimensionshifting relations given in Appendix B to rewrite them in terms of the masters we have computed. Then, we decompose the three-loop corrections as

Conclusion
In this paper we have presented for the first time fully analytic results in terms of eMPLs and iterated Eisenstein integrals for the three-loop corrections to the ρ parameter in the SM with two massive quark flavours. This computation was originally considered as an expansion in the ratio of the quark masses in ref. [44]. An important ingredient in our calculation is the realisation that the homogeneous second-order differential operator appearing in this computation is identical to the differential operator that appears in the computation of the well-known sunrise graph. As a consequence, all integrals can be expressed in terms of the exact same class of functions as the sunrise graph, which are also well-studied functions in pure mathematics. We can draw upon the knowledge of these functions to analytically continue them using tools and algorithms developed for the sunrise graph. This distinguishes our computation from the results of refs. [45,46], where a novel class of special functions was introduced for the same class of integrals.
Besides ref. [73], our computation is only the second time that iterated integrals of modular forms have been used to obtain fully analytic results for a complete physical observable. We believe that the techniques that we have used in our computation can have an impact also on the computation of other physical observables, and that they pave the way for obtaining more results involving this class of special functions.

A Differential equations for the non-elliptic master integrals
As written in eq. (2.11), the differential equations for the master integrals which are expressible in terms of MPLs, i.e. f 1 through f 7 of eq. (2.10), is given by The matrices A 0 and A 1 are given by