Two-loop Sunset Integrals at Finite Volume

We show how to compute the two-loop sunset integrals at finite volume, for non-degenerate masses and non-zero momentum. We present results for all integrals that appear in the Chiral Perturbation Therory ($\chi$PT) calculation of the pseudoscalar meson masses and decay constants at NNLO, including the case of Partially Quenched $\chi$PT. We also provide numerical implementations of the finite-volume sunset integrals, and review the results for one-loop integrals at finite volume.


Introduction
An analytical, ab initio description of Quantum Chromodynamics (QCD) in the hadronic low-energy regime remains elusive. One of the most promising alternatives involves numerical evaluation of the functional integral of QCD on a discretized space-time lattice. Known as Lattice QCD, this approach has long been restricted for computational reasons to large and unphysical values of the light quark masses. Recently, due to improvements in computing power and algorithmics, calculations with significantly smaller quark masses have become possible. A side effect of the lowered quark masses is an increase in the size of finite-volume corrections, and a detailed treatment of such effects is thus called for.
Fortunately, in many cases the finite-volume corrections can be evaluated analytically using Chiral Perturbation Theory (χPT) [1,2], which is the low-energy effective theory of QCD. The application of χPT at finite volume was first performed by Gasser and Leutwyler in Ref. [3], and a review of recent work in this area can be found in Ref. [4]. As many Lattice QCD simulations are performed with unequal valence and sea-quark masses [5], the properties of the light pseudoscalar mesons have also been calculated to next-to-nextto-leading order (NNLO) in Partially Quenched χPT (PQχPT) in Refs. [6][7][8]. Therefore, it is also of interest to extend the finite-volume description of the relevant loop integrals to account for the appearance of double poles in the PQχPT propagators. It should be noted that χPT is applicable at finite volume as soon as the typical momenta of a given process are sufficiently small. This imposes the restriction F π L > 1, where F π is the pion decay constant and the volume V ≡ L 3 . This study deals with the so-called p-regime, in which V is sufficiently large for zero-momentum fluctuations of the meson fields to be treated perturbatively, which introduces the additional requirement m 2 π F 2 π V ≫ 1, where m π is the pion mass. A multitude of finite-volume calculations exist at one-loop or next-to-leading order (NLO), and it should also be noted that some work at NNLO has recently appeared. This includes Ref. [9], where the finite-volume corrections to the quark condensate were calculated, and Ref. [10] which considered m π for the case of degenerate quark masses.
Our main objective is to show how the integrals needed in χPT calculations of pseudoscalar meson properties at NNLO and finite volume can be performed. As a starting point, the known results at one-loop order are reviewed, and we also show how these can be extended to higher order in d − 4. The methods for the one-loop integrals are then applied to the two-loop "sunset" integrals for arbitrary masses and momenta. We focus here on the integrals necessary for the calculation of form factors to NLO, and for the calculation of masses, decay constants and two-point functions to two loops (NNLO). This paper is structured as follows: Section 2 discusses a few preliminaries. In Section 3, the derivation of the one-loop integrals at finite volume is revisited, with emphasis on the treatment of PQχPT calculations at NLO. In Section 4, the two-loop sunset integrals are considered, and explicit expressions are given for the finite and divergent parts, for arbitrary values of the quark masses and with the propagator structure of PQχPT fully accounted for. Section 5 contains a numerical overview of the integrals presented in this study, along with a concluding discussion in Section 6. The appendices summarize the ingredients involving modified Bessel functions and theta functions, along with basic inte-grals in d dimensions and comments on the notational conventions in earlier work. Some preliminary results related to this study have been presented in Refs. [11,12].

Finite-volume sums
At finite volume in a cubic box, integrals over momenta should be replaced by sums over the allowed momenta. In one dimension of length L, with periodic boundary conditions, 1 this entails a summation over the allowed momenta p n ≡ 2πn/L, with n ∈ Z integer. The integrals over momenta should thus be replaced according to dp 2π where the latter notation will be used to indicate a finite-volume summation in the remainder of this paper. Infinities will be treated by dimensional regularization, using the convention d ≡ 4 − 2ε. The infinite-volume integrals have been treated extensively in the literature, see e.g. Ref. [14] including appendices and references therein. In practice, it is often desirable to study deviations from the infinite-volume limit, and we shall therefore use a framework in which the infinite-volume contribution can be easily identified. This can be achieved by application of the Poisson summation formula to Eq. (2.1), yielding 1 L n∈Z F (p n ) = lp dp 2π e ilpp F (p), (2.2) where the summation over l p spans a set of vectors of length nL such that n ∈ Z. The term with n = 0 then represents the infinite-volume result, while the sum of all the other terms is the finite-volume correction.
In the case of loop integrals over momenta in higher dimensions, Eq. (2.2) should be applied to all dimensions which have a finite extent. The four-vector l pµ then has components (0, n 1 L, n 2 L, n 3 L) when three of the dimensions have a finite extent L. The loop integrals in this paper are performed throughout in Euclidean space, with metric g µν = δ µν and signature (+, +, +, +). Throughout this paper, one of the dimensions (the "time" dimension) is assumed to be much larger in extent than the other three dimensions, which is the usual situation encountered in Lattice QCD.

Passarino-Veltman reduction
At infinite volume, a general method was developed by Passarino and Veltman [15] to obtain a minimal set of integrals by reduction of the tensor integrals H µν to a set of scalar integrals. This method relies on separation of the integrals into components that are scalars under Lorentz transformations and prefactors that contain δ µν and various momenta. Although Lorentz-invariance is explicitly broken by the introduction of a finite size, it is still possible, in the frame where p · l p = 0, to rewrite the integrals in scalar components, provided that a four-vector is introduced. The situation p · l p = 0 is referred to as the "center-of-mass" (cms) frame, which is a situation often realized in Lattice QCD. Because of the remaining symmetries, t µ is the only additional object required to rewrite the integrals in scalar components, but we also introduce the tensor as a convenient additional abbreviation.

One-loop integrals at finite volume
In general, the one-loop integrals in the NNLO expressions for the pseudoscalar meson masses and decay constants contain a maximum of two propagators with distinct masses. The simplest case with one propagator is denoted A, whereas the case with two distinct propagators is denoted B. In PQχPT, some three-propagator integrals denoted C also appear. These are due to the mixing of different lowest-order states in PQχPT, and they can always be re-expressed in terms of the B integrals. All of the integrals mentioned above have been extensively treated in the literature, see e.g. Refs. [3,[16][17][18]. However, it is instructive to review certain aspects of their derivation and numerical evaluation here, since they form building blocks in the calculation of the two-loop sunset integrals at finite volume.

One-propagator integrals
The basic one-loop, one-propagator integrals are where X = 1, r µ and r µ r ν . By application of the Poisson summation formula for the finite dimensions, Eq. (3.1) may be written as where the term with l r = 0 represents the infinite-volume contribution. In order to isolate the finite-volume part, Eq. (3.2) is decomposed according to where the first term represents the infinite-volume result and will not be considered further. The second term represents the finite-volume correction, and is free from divergences.
First, we consider the case of X = 1. We rewrite Eq. (3.1) using Eq. (A.1) as where the primed sum indicates that the term with l r = 0 is excluded. We next substitute r ≡r + il r /(2λ), and obtain where ther integral can be performed using Eq. (C.3) and by rescalingr ≡r/ √ λ, which gives The (triple) sum and integral can be evaluated in different ways. The technique used in Refs. [3,16] is to employ Eq. (A.2), which yields where the modified Bessel functions K ν are defined in App. A. The triple sum can be simplified by observing that l 2 r = kL 2 , with k integer. We further define the factor x(k), which indicates the number of times each value of k ≡ n 2 1 + n 2 2 + n 2 3 occurs in the triple sum. We then find x(k)f (k), (3.8) which reduces the triple sum to a single sum. The final result is where the arguments of K ν can be modified by rescaling λ before Eq. (A.2) is applied. Also, the sum over modified Bessel functions is found to converge fairly slowly. The second method considered here involves performing the summation, and leaving the integral to be evaluated numerically, see Ref. [18]. We observe that using the relation l 2 r = (l 2 1 + l 2 2 + l 2 3 )L 2 . The cubic power accounts for the summations over l 1 , l 2 and l 3 . The remaining sum in Eq. (3.10) can be performed in terms of the theta function θ 30 , which is defined in App. B. This gives where, as a final step, we rescale λ to obtain which is also valid for mL ∼ 1. Integrals with factors of r µ in the numerator are also required. Up to NNLO, these are ⌊r µ ⌋ and ⌊r µ r ν ⌋. Proceeding as above, we obtain where we note that integrals odd inr vanish, and that (3.14) The summations over the components of l r include both positive and negative contributions, and are symmetric under interchange of spatial directions. The sums which are odd in the components of l r then vanish, and lr l rµ l rν f (l 2 r ) = Thus, the final results for the ⌊r µ ⌋ and ⌊r µ r ν ⌋ integrals are where the remaining integration can again be performed in terms of the modified Bessel functions, giving (3.17) For the second method which involves the theta functions, we rewrite the sum using the identity which is also valid for the primed sums, as the term with n = 0 does not contribute. After rescaling λ, this gives and following the same steps as before, we also find

Two-propagator integrals
The basic one-loop, two-propagator integrals are where X = 1, r µ , r µ r ν and r µ r ν r α . By application of the Poisson summation formula for the finite dimensions, Eq. (3.21) may be written as where the term with l r = 0 represents the infinite-volume contribution. We again decompose Eq. (3.22) into the infinite-volume part and the finite-volume correction using where the latter term is obtained from Eq. (3.22) by replacing the unprimed sum with the primed one, indicating that the term with l r = 0 is excluded. The methods of Sect. 3.1 also apply here. We begin by introducing Gaussian parameterizations for both propagators in Eq. (3.22) in terms of the integration variables λ 1 and λ 2 . In a second step, we switch to a new set of variables (λ, x) with λ 1 ≡ xλ and λ 2 ≡ (1 − x)λ. Alternatively, we may first combine the two propagators using the Feynman parameterization 1 where y = 1 − x, and then treat the denominator according to Eq. (A.1). In both cases, the result is which is equivalent to Eq. (3.4). We now shift the integration variable to r ≡r+il/(2λ)+yp and obtain for the simplest case which differs from Eq. (3.6) by the integration over x and the factor e iylr·p . Due to the summation over components of l rµ with alternating signs, this factor always produces realvalued results. For the remaining integrals, we obtain

Center-of-mass frame
In the cms frame, p = (p, 0, 0, 0) such that p · l r = 0 for all l r . The integrals in the cms frame can be computed similarly to the one-propagator integrals, giving where the subscripts of the ⌊X⌋ V indicate the value of n in the one-propagator integrals given in Sect. 3.1. Also, the one-propagator integrals in the above expressions are functions ofm 2 rather than m 2 . We may then compute the integral over λ in ⌊X⌋ V and obtain a sum over modified Bessel functions. We are finally left with a single summation and an integral over x, to be performed numerically. The method introduced in Sect. 3.1 where the summations are performed in terms of theta functions is also applicable here, and yields a double integral over λ and x. In that case, the integral over x can be performed analytically. By setting the resulting integral with no additional powers of z is related to Dawson's integral or the error function (erf), depending on the sign of p 2 . The other cases are related to the (complex-valued) incomplete Gamma function by the substitution z 2 = u. However, a straightforward numerical evaluation of the double integral converges sufficiently fast for practical purposes.

Moving frame
In a general "moving frame", p can have non-zero components in the dimensions of finite length. In this case, the sums with odd powers of components of l r no longer vanish. In general, the finite-volume corrections can depend on all components of p, and no simple way of writing the result in terms of scalar functions of p 2 exists, as only a discrete subgroup of the three-dimensional rotation group remains as a symmetry in a finite cubic volume. Nevertheless, the relevant expressions can be evaluated numerically, albeit with some additional complications. For the formulation in terms of modified Bessel functions, the summation is no longer exclusively dependent on l 2 r , and thus the reduction of the triple sums using Eq. (3.8) is no longer possible. For the formulation in terms of theta functions, the summation over l r can still be performed separately for each dimension, provided that the factors of θ 3 30 are replaced by the product θ 3 (u 1 , q) θ 3 (u 2 , q)θ 3 (u 3 , q), where u i ≡ yp i L/(2π) and q ≡ e −1/λ . When factors of r µ appear in the integrands, derivatives w.r.t. u and q, as well as uncontracted factors of l rµ , also need to be accounted for.

Summary of one-loop results
Next, we discuss the relations between the various one-loop integrals and summarize the explicit expressions in a concise form. With the definition of Eq. (3.1) in mind, we introduce the more conventional notation where only the finite-volume correction has been retained. As discussed above, no simple rewriting in scalar components is possible for the momentum-dependent integrals, except in the cms frame with p = (p, 0, 0, 0). In that frame, we define which correspond to the usual definitions at infinite volume, except for the terms involving t µν , which appear only in the finite-volume contribution. The Passarino-Veltman construction [15] produces relations between the various integrals upon multiplication with p µ or δ µν . Using a number of relations can be obtained. These are and where we note that the relations in Eq. (3.36) are linearly dependent. Up to the order considered here, this leaves A V , A V 23 , B V and B V 23 as independent functions. We have checked the validity of the above relations numerically for n 1 , n 2 = 1, 2.
At NNLO in χPT, all one-loop integrals should be expanded around d = 4 up to and including terms of O(ε). This is necessary, since products of two one-loop integrals appear throughout the NNLO expressions, including the factorizable parts of the two-loop sunset integrals. We thus define with similar expansions for all functions A V i and B V i in Eqs. (3.32) and (3.33). The onepropagator integrals can then be written as using Eqs. (3.9), (3.12), (3.17) and (3.19). The integrands can be expressed either in terms of modified Bessel functions or theta functions, and are in each case given bŷ The expansion in ε = (4 − d)/2 can be performed using where the functionsK m are related to the modified Bessel functions and are defined in App. A. For all quantities in Eq. (3.39), the above results lead tô where K m →K m indicates that the functions K m should be replaced by the corresponding expressions forK m . For the one-loop two-propagator integrals, we find similar results, given bȳ The explicit expressions for the integrands arê where each case has again been given in terms of modified Bessel functions or theta functions. The functionsB V ε can be obtained from the above expressions using the equivalent of Eq. (3.40), along with corresponding changes in Eq. (3.41). However, the functionsĀ V ε andB V ε are expected to cancel completely in a full calculation within the M S scheme. This cancellation has already been demonstrated at NNLO for the scalar condensate in Ref. [9], and for m π in two-flavour ChPT in Ref. [10].

Two-loop sunset integrals at finite volume
First, we recall that some NNLO work at finite volume already exists. In Ref. [9], the finitevolume corrections were calculated for the quark condensate, and in Ref. [10] for m π . The former only involved products of one-loop integrals, while the latter only required consideration of the sunset integrals with degenerate masses. In this section, we provide completely general expressions for the sunset integrals, for arbitrary, non-degenerate masses. At finite volume, we define the basic sunset integral as where the required operators X are 1, r µ , s µ , r µ r ν , r µ s ν and s µ s ν . In Eq. (4.1), the n i are always non-zero and positive. If one of the n i is zero or negative, the integral becomes separable into a product of two one-loop integrals, which we have already dealt with in Section 3.
Application of the Poisson summation formula for all momenta in a finite dimension yields where X (1, 2, 3) will be used as a short-hand notation indicating which of the arguments (n i ,m 2 i ) are associated with the first, second and third propagators in Eq. (4.2), respectively. The vectors l r , l s are of the form (0, k 1 L, k 2 L, k 3 L) with k i ∈ Z. Eq. (4.2) can then be decomposed according to where X ∞ denotes the infinite-volume result with l r = l s = 0. The sunset integrals at infinite volume have been evaluated in several different ways (see e.g. Refs. [19][20][21][22]) and will not be considered further here. The second term in Eq. (4.3) represents the finitevolume correction. The present approach to the finite-volume correction is along the lines of Refs. [19,20], combined with an extension of the methods for the one-loop integrals in Section 3. We further decompose X V into terms where one of the possible loop momenta is not quantized and a contribution where both are quantized, according to where a "singly primed" sum indicates that the term with l = 0 has been excluded. For the "doubly primed" sums, all contributions with l r = 0, l s = 0 or l r = l s have been removed, i.e. the retained terms satisfy l r = 0, l s = 0 and l r = l s . The sum of all the terms in Eq. (4.5) reproduces the full sum in Eq. (4.2). Here, it should be taken into account that p is also quantized in the finite dimensions, such that the spatial momentum components satisfy We note that X rs is always finite, whereas X r , X s and X t may contain a non-local divergence, depending on the operator X and the values of the n i . If these integrals should be finite, they can be included in X rs by summation over all values of l r and l s (except of course l r = l s = 0).

Simplest sunset integral
We first restrict ourselves to the simplest case of 1 with n 1 = n 2 = n 3 = 1, which allows us to outline our procedure in a straightforward way. We will then proceed to give the expressions for the general case using the formalism established here.
From Eqs. (4.1), (4.2) and (4.5), and keeping in mind Eq. (4.6), we find that the sunset integrals exhibit a high degree of symmetry with respect to interchanges of r, s and t = p − r − s, together with l r , l s and l t . Substituting (r, s) → (s, r) and (r, t) → (t, r), including the respective l i , leads to the relations where we recall that the notation (1, 2, 3) refers to the propagators, as exhibited in Eq. (4.2). From the last relation in Eq. (4.7), we find that the evaluation of 1 r and 1 rs suffices to obtain the full result.

Simplest sunset integral with one quantized loop momentum
First, we calculate 1 r . We begin by combining two of the propagators with a Feynman parameter x, giving where we have shifted the integration variable according to s µ ≡s µ − x(r − p) µ , and defined The integration overs may then be performed in terms of standard d-dimensional integrals in Euclidean space, given in App. C.1. This gives where the expansion to O(ε) may be performed using where λ 0 ≡ 1/ε + log(4π) + 1 − γ. The term proportional to λ 0 involves the one-loop integral A V , which has been treated in Sect. 3. This also contains the nonlocal divergence, and contributes to 1 r . For clarity, we have added the arguments n 1 = 1 and m 2 1 to the notation for the one-loop integral. The remaining terms in Eq. (4.11) contribute where we can set d = 4 directly. In order to deal with the dependence of m 2 or r, we perform a partial integration in x to obtain (4.14) Here, the first two terms once more contain a one-loop integral, and we refer to this part as 1 r,G , with the remainder labeled 1 r,H . Further, we introduce the Gaussian parameters λ 1 and λ 4 according to Eq. (A.1) for the denominators (r 2 + m 2 1 ) and m 2 , respectively. This gives 1 r,F ≡ 1 r,G + 1 r,H , where we may complete the square in the exponential factor by substituting Ther integral can then be performed using Eq. (C.3), which gives Here, a more symmetric form can be obtained by substituting λ 2 ≡ (1 − x)λ 4 and λ 3 ≡ xλ 4 as integration variables, giving with which can be evaluated numerically with the methods discussed in Sect. 4.1.3.

Simplest sunset integral with two quantized loop momenta
Second, we calculate 1 rs . We introduce Gaussian parameterizations for all three propagators using Eq. (A.1) and set d = 4, giving after which we perform the redefinition and shift s by where we have again made use ofλ ≡ λ 1 λ 2 + λ 2 λ 3 + λ 3 λ 1 . We note that an analogous transformation results by first redefining s and then shifting r. The result is We note that the arguments of the exponential functions in Eqs. (4.18) and (4.23) coincide when l s = 0.

Numerical evaluation
Next, we discuss the numerical evaluation of Eq. (4.23). For this purpose, it is convenient to switch to the variables x, y, z and λ, where σ ≡ xy + yz + zx and x + y + z = 1. We also introduce the quantities As for the one-loop integrals, we may either perform the summations in terms of theta functions, or the λ integration in terms of modified Bessel functions. In terms of the latter, the result is where we note that in the cms frame where S rs = 0, we may write similarly to Eq. (3.8). Here, the factor x(k r , k s , k n ) denotes the number of times a given triplet of squares appears when the components of l r and l s are varied over all positive and negative integer values. In terms of theta functions, we find in the cms frame where the contributions with l 2 r , l 2 s or l 2 n equal to zero have been subtracted. The Jacobi and Riemann theta functions are defined in App. B, see also Eq. (4.91) and the accompanying discussion.
The expression for 1 r,H in Eq. (4.18) is clearly similar and can be treated along the same lines. In terms of modified Bessel functions, the terms with a single sum over l r may be treated similarly to the one-loop integrals using Eq. (3.8). Alternatively, the summation can be performed in terms of theta functions. The relevant expressions will be given when we summarize the full results for the sunset integrals.

Permutation properties
The finite-volume sunset integrals satisfy a number of relations which simplify the calculations, and provide useful checks on the numerics. These are the more general versions of Eq. (4.7). When applied to the full sunset integrals X , the variable interchanges (s, r), where the latter one follows from the identity from which it also follows that all parts of r µ s ν that are symmetric in µ and ν can be rewritten in terms of other integrals. In particular, at infinite volume r µ s ν can be expressed in terms of r µ r ν using various permutations of the m 2 i and n i . This also holds for the case of m 1 = m 2 and n 1 = n 2 . The relations (4.31) and (4.32) are also separately valid for X ∞ , X V and X rs , but not for the other components of Eq. (4.4).
From the above considerations, we can deduce what integrals should be calculated in order to obtain a complete description. As X r , X s and X t are closely related, we can obtain the required cases of X s using 1 s (1, 2, 3) = 1 r (2, 1, 3; l r → l s ), r µ s (1, 2, 3) = s µ r (2, 1, 3; l r → l s ), r µ r ν s (1, 2, 3) = s µ s ν r (2, 1, 3; l r → l s ), r µ s ν s (1, 2, 3) = s µ r ν r (2, 1, 3; l r → l s ), (4.34) and for the X t we find 2 1 t (1, 2, 3) = 1 r (3, 2, 1; l r → −l t ), from which we conclude that a complete description entails the calculation of X rs for X = 1, r µ , r µ r ν and r µ s ν , and of X r for X = 1, r µ , s µ , r µ r ν , r µ s ν and s µ s ν . We also note that the X r are symmetric under the interchange (m 2 , n 2 ) ↔ (m 3 , n 3 ) for X = 1, r µ and r µ r ν . For conciseness, we now introduce a set of functions to be used in the remainder of the text. In an arbitrary frame, we define and H V 28µν contain instances of the vectors l r or l s with uncontracted Lorentz indices. In the cms frame, such contributions with one Lorentz index vanish, and the bilinear ones become proportional to t µν . In the cms frame, we therefore have a simplified set of functions Because of this structure, r µ s ν is symmetric in µ, ν and can be obtained using Eq. (4.32). Still, we include r µ s ν as a useful check on our numerics, and because it appears in the expressions for the sunset integrals with one quantized loop momentum. Our numbering scheme for the sunset integrals has been chosen to be consistent with Ref. [19]. We also refer to the components of the functions H i by appending the indices (r, G), (r, H) etc., which were introduced in the detailed treatment of the simplest sunset integral.

Sunset integrals with one quantized loop momentum
Here, we follow along the lines of Sect. 4.1.1 and account for all needed cases of X r with X = 1, r µ , s µ , r µ r ν , r µ s ν and s µ s ν . Again, the first step is to combine the last two propagators with a Feynman parameter x and shift the integration variable by s µ ≡ s − x(r − p) µ . The integral overs can then be performed using Eq. (C.1). Using the notation f (r α ) for additional factors of r µ , r ν , this gives where the remaining integral over r is always finite because of the factor e ilr·r . It is then sufficient to expand thes integral in ε, while keeping only the singular and O(1) terms as in Eq. (4.11). We rewrite the singular terms using λ 0 ≡ 1/ε + ln(4π) + 1 − γ, and define the components of the sunset integrals proportional to λ 0 with the subscript A as in Eq. (4.12).
In terms of the one-loop integrals defined in Sect. 3, we find for the non-zero cases with n 2 , n 3 = 1, 2 the expressions where the superscripts denote the n i in the sunset integrals. Also, the one-loop integrals now show explicitly the m 2 i and n i of the denominator they involve. As the above integrals contain a non-local divergence, they should always cancel in physical results.
We now proceed to treat the terms containing log(m 2 ). As before, we first perform a partial integration in x, giving after which we denote the terms with negative powers of m 2 as X r,H , and the others as X r,G , as defined in Eq. (4.15) for the case of the simplest sunset integral. We note that the X r,G can again be expressed in terms of one-loop integrals. For n 2 , n 3 = 1, 2, the non-zero cases are We note that the decomposition of the parts of the sunset integrals which do not depend on λ 0 into X r,G and X r,H is clearly not unique, as it depends on the choice of Feynman parameterization. For example, had we chosen y = 1 − x instead of x as the Feynman parameter, we would have obtained terms containing log(m 2 2 ) in the X r,G . Also, the decomposition does not commute with derivatives w.r.t. masses, note e.g. that 1 n 1 12 r,G = 0 = −(∂/∂m 2 3 ) 1 n 1 11 r,G . The remaining part X r,H is algebraically the most complicated, but again follows exactly the procedure for the simplest sunset integral. First, we introduce Gaussian parameterizations for the negative powers of m 2 and (r 2 + m 2 1 ) using Eq. (A.1) with parameters λ 4 and λ 1 , respectively. While the expressions corresponding to 1 r,H in Eq. (4.15) are relatively lengthy, they all share the same basic structure. In particular, they all contain the same exponential factor, for which we may complete the square using the substitutions of Eq. (4.16). The resulting integrals can then be performed by means of Eq. (C.3). Finally, we define λ 2 ≡ (1 − x)λ 4 and λ 3 ≡ xλ 4 and perform the substitutions of Eq. (4.25) to obtain an integral in terms of x, y, z and λ.
Before we give explicit expressions for X r,H , we briefly discuss the methods used to obtain them. Due to the complexity of the required analytical manipulations, we have found it convenient to use FORM [23] according to the procedure outlined above. Alternatively, as described in Ref. [11], a number of tricks can be used to considerably simplify the task. For example, powers of r µ can be introduced into the numerators of the sunset integrals by taking derivatives w.r.t. l r , giving It is also noteworthy that integrals such as s µ r are very similar to the case of 1 r , differing only in an additional factor of x(r − p) µ . This leads to relations such as where the factor of x is understood to be included in the respective integrals. Due to the length and complexity of the resulting expressions for X r,H , we make use of the auxiliary quantities and and we also introduce the notation X n 1 n 2 n 3 r,H for X = r µ , s µ . With {a, b} µν = a µ b ν + a ν b µ , we find for the bilinear operators Given these expressions for X r,H , we may proceed as for the one-loop integrals and choose between performing the summations in terms of theta functions, or evaluating the λ integral in terms of modified Bessel functions. The results quoted in Eqs. (4.49)-(4.54) make no assumptions on the momentum p. Below, we restrict ourselves to the cms frame where p · l r = 0 or p = (p, 0, 0, 0). This case is the most commonly encountered, and the expressions for a moving frame can be obtained along similar lines.

Center-of-mass frame: Bessel functions
Here, we have performed the integration over λ in terms of the functions K ν (Y, Z) defined in App. (A.2). We note that the summation only depends on l 2 r , such that Eq. (3.8) is applicable. We have suppressed the arguments (Y, Z) in order to keep the expressions short and concise. The expressions always contain the abbreviated part

Center-of-mass frame: Theta functions
Next, instead of computing the integrals over x, y and λ, we have performed the summation in terms of the theta functions, previously encountered for the one-loop and simplest sunset integrals. In the cms frame, we make use of Eqs. (3.10), (3.18), and where we note that Eq. (4.68) can immediately be used for the primed sums by setting l r = nL, as the term with n = 0 does not contribute. We rescale λ such that the argument of all theta functions is e −1/λ , which we suppress for brevity. Further, we introduce the abbreviation H r,H;n 1 12

Sunset integrals with two quantized loop momenta
Here, we follow the treatment of Sect. 4.1.2, and generalize to all integrals X rs with X = 1, r µ , s µ , r µ r ν , r µ s ν and s µ s ν . All of these are not needed for completeness, but the redundant ones enable a check on our results by means of the relations given in Sect. 4.2.
We again introduce Gaussian parameterizations for the propagators using Eq. (A.1), and then shift the momenta using Eqs. (4.21) and (4.22). This leads to X rs = 1 Γ(n 1 )Γ(n 2 )Γ(n 3 )(4π) d Table 1. Numerical results for the one-propagator "tadpole" integrals, for m = 0.1395 GeV, which corresponds to mL ≈ 2.12 (L = 3 fm) and mL ≈ 2.83 (L = 4 fm). The corresponding continuum integrals are shown in the column labeled L = ∞. The continuum results employ the MS subtraction scheme with µ = 0.77 GeV. Note that the "23" case has no continuum counterpart. All results are given in units of the appropriate powers of GeV, and the pole configurations n of the propagators are given in App. D.

E Translation to Minkowski conventions
While we have used the Euclidean formalism throughout, it is also of interest to recall how the expressions for the one-loop and sunset integrals can be translated to Minkowski conventions. The required substitutions are where t µν corresponds to the spatial part of the metric.