The scalar sunset diagram at finite temperature with imaginary square masses

We evaluate the finite temperature scalar sunset diagram with imaginary square masses, that appears in the Gribov-Zwanziger approach to Yang-Mills (YM) theory beyond one-loop order. Since YM theory at finite temperature is governed by center-symmetry and the Polyakov loop, we also include the possibility of a constant temporal background gauge field in the form of color-dependent imaginary chemical potentials.


I. INTRODUCTION
In recent years, much valuable progress has been made towards the understanding of non-abelian gauge theories at finite temperature using background field gauge (BFG) methods [1,2] in the Landau-DeWitt gauge, in combination with several functional methods [3][4][5][6][7][8][9][10][11][12][13][14]. On the one hand, BFG methods provide an efficient way to describe the confinement/deconfinement order parameter (the Polyakov loop or any of its proxies [3]) because the related center symmetry is explicit at the quantum level and is easily maintained in approximation schemes [15][16][17]. On the other hand, functional methods provide a method of choice when investigating infrared, non-perturbative properties of non-abelian theories [18].
However, most functional approaches take as a starting point the usual Faddeev-Popov version of the gauge fixing which is known to be a valid description of non-abelian gauge arXiv:2006.12184v1 [hep-th] 22 Jun 2020 theories at high energies but which is also expected to be modified in the infrared due to the influence of Gribov copies [19]. It is then an interesting question whether a complete gaugefixing procedure in the infrared (IR) regime could capture some genuine non-perturbative effects, beyond those that are captured by the infinite hierarchies of equations considered in functional methods. Even more, it has been suggested that a resolution of the IR gaugefixing may open the way to a new perturbative perspective on certain aspects of the infrared dynamics of non-abelian gauge fields [20,21].
Several models have been put forward in order to implement the BFG formalism in the Landau-deWitt gauge while restricting the number of Gribov copies. In Refs. [16,[22][23][24][25], the formalism was used within the Curci Ferrari (CF) model [26] to compute the background potential and Polyakov loop up to two-loop order, both in pure Yang-Mills theories and in heavy-quark QCD. It was argued that this model could be part of a complete gauge-fixing in the Landau gauge, since a CF gluon mass term may arise after the Gribov copies have been accounted for via an averaging procedure [27], see also Ref. [28] for a related discussion in a different gauge. One salient feature of the results obtained within the CF model is that, not only various aspects of the phase structure are already accounted for at leading one-loop order, but the two-loop corrections turn out to be small and tend to improve the results, supporting the existence of a "perturbative way" lurking behind the gauge-fixing problem.
A more explicit way to account for Gribov copies is the Gribov-Zwanziger (GZ) method [29][30][31]. At the cost of introducing some new fields, the functional integral is restricted to a region that contains less Gribov copies, the so-called Gribov region. In Ref. [32], a GZ type action for the Landau-DeWitt gauge was proposed, but it was later established in Ref. [33] that this model is not invariant under background gauge transformations and an alternative proposal was made where both background gauge invariance and BRST symmetry are manifest. It is however not clear how to extend this proposal at finite temperature while maintaining the background gauge invariance. Here, we will follow the framework of Ref. [34], where BRST symmetry is sacrificed (just as in the CF model) to establish a background gauge invariant GZ type model that is easy to implement at finite temperature.
In Ref. [34], the one-loop background potential and the Polyakov loop up to first order were determined within this model and in Ref. [35] these calculations were extended to the case of QCD with heavy quarks, leading to the the best agreement to date with the available lattice data regarding the description of the upper boundary line in the so-called Columbia plot. A natural question is whether these promising results at leading order resist the inclusion of higher order corrections, which would support similar results within the CF approach.
The present work is a modest contribution towards this goal: we address the calculation of the scalar sunset diagram and the mass derivatives thereof that appear in the two-loop background potential in the model at finite temperature. This potential puts forward new challenges because imaginary square masses appear with the introduction of the auxiliary fields needed to localize the GZ action. Indeed, the tree-level gluon propagator in the GZ model reads with γ the Gribov parameter. Though the existence of imaginary masses in the GZ model is a well-known fact, to our knowledge there is no literature on the proper handling of imaginary masses in higher-order loop calculations at finite temperature. The full calculation of the two-loop potential in the GZ model as well as the Polyakov loop will be treated in a different work [36]. For related work at zero temperature, see [37][38][39].
A convenient tool to make sense of the finite temperature contributions to the potential is thermal splitting which is commonly used in calculations that involve Matsubara sums [40,41]. By decomposing sum-integrals according to the number of thermal factors, UV divergences become much easier to handle. Moreover, we can separate a vacuum piece, which will equal the zero-temperature contribution. The vacuum two-loop sunset amplitude for real masses was calculated in [42] and the finite temperature contributions have been known for a long time [43], with a recent generalization in the presence of the Polyakov loop [16,23,44]. Part of this work will therefore be an extension of these results to the case of imaginary square masses. We do not aim at a full generalization, however, instead limiting ourselves to the cases that appear in the two-loop calculation in the GZ model [36].
This work is organized as follows. In Sec. II, we look at the scalar tadpole sum-integral as a pedagogical introduction to the techniques that will be used to deal with the sunset sum-integral. In particular, we introduce the spectral representation and give a first trivial example of thermal splitting. In Sec. III, we look at the scalar sunset sum-integral. In Sec. IV, we investigate the relevant mass derivatives of the sunset sum-integrals and their respective thermal splittings that are also needed for the evaluation of the GZ potential at two-loop order. More technical details are gathered in the Appendices.

II. THE SCALAR TADPOLE AS A SIMPLE EXAMPLE
In what follows, we denote Euclidean momenta by capital letters Q, K, L, . . . Each of these comprises a bosonic Matsubara frequency ω n ≡ 2πT n, with n ∈ Z, and a spatial momentum q , with q ≡ |q|. Integration over Euclidean momenta is encoded in sum-integrals, which we keep denoting, however, as standard integrals for simplicity: We work in dimensional regularization, meaning that the integral over spatial momenta with d = 4 − 2 .
In the context of Yang-Mills theory at finite temperature, it is crucial to take into account the order parameter for the confinement/deconfinement transition also known as the Polyakov loop , or, equivalently, the corresponding constant, temporal and diagonal gluonic backgroundĀ 0 such that ∝ tr exp{igβĀ 0 }, with β ≡ 1/T the inverse temperature.
In this situation, the Matsubara frequencies are shifted by a color-dependent imaginary chemical potential, ω n → ω κ n ≡ ω n +r j κ j , where ther j denote the components ofĀ 0 along the diagonal part {t j } of the su(N ) algebra,Ā 0 =r j t j , while the κ j denote the weights of the adjoint representation, that arise as one diagonalizes the adjoint action of all the t j : [t j , t κ ] = κ j t κ . For the present paper, we do not need to know more about the precise way the Polyakov loop appears in explicit calculations. In what follows, we denote by Q κ = (ω κ n ; q) the shifted Euclidean momentum and we also introduce the notationr · κ ≡r j κ j .
In this first section, as a pedagogical example, we treat the scalar tadpole sum-integral with assuming that the square mass α is purely imaginary. The procedure that follows might seem unnecessarily complicated for the evaluation of such a simple sum-integral. However, it introduces the basic ingredients that make the corresponding evaluation of the scalar sunset sum-integral in the next section much simpler.

A. Spectral representation
The first step is to evaluate the Matsubara sum in Eq. (3). To this purpose, we decompose the propagator as with ε q,α ≡ q 2 + α. It proves useful to rewrite the previous identity in the form of a "spectral representation" where q 0 ≡ dq 0 /(2π) and We mention that, in the presence of imaginary square masses, the notations q 0 and δ(q 0 ∓ ε q,α ) are understood as mere bookkeeping devices allowing to select the two complex energies ±ε q,α . Similarly, sign(q 0 ) selects the corresponding sign in front of ±ε q,α , and should therefore be understood as the sign of the real part of q 0 . 1 Using the spectral representation (6) in Eq. (3), we find with a simple Matsubara sum. We stress that, even though q 0 takes complex values, it does not interfere with the Matsubara frequencies because its real part never vanishes. Using 1 The spectral representation can be given a rigorous meaning by defining the Dirac and sign distributions along the appropriate contour. We will not need these technicalities here though.
standard techniques for the evaluation of Matsubara sums, we then arrive at 2 with n x ≡ 1/(e x/T − 1) the Bose-Einstein distribution function.

B. Thermal splitting
One problem with the expression above is that it involves thermal factors with energies whose real parts can be as negative as possible. In particular, this does not facilitate the extraction of UV divergences. To remedy this situation, we write where sign(q 0 ) is to be understood as the sign of the real part of q 0 , see above, θ(q 0 ) is equal to 1 if the real part of q 0 is positive and zero otherwise, and |q 0 | = q 0 sign(q 0 ). Plugging Eq. (11) into Eq. (10) and then back into Eq. (8), we arrive at the "thermal splitting" of the tadpole sum-integral: J κ α = J α (0n) + J κ α (1n). Here, J α (0n) denotes the pure vacuum contribution (no thermal factor), depending neither on the temperature nor on the background, while where σ κ α (q 0 ; q) ≡ ρ α (q 0 ; q) sign(q 0 ) n |q 0 |−i sign(q 0 )r·κ . In the contribution with one thermal factor, one can perform the frequency integral.
Moreover, because this contribution is UV finite, one can take the limit d → 4 and evaluate the angular integral analytically. One obtains where σ α ∈ {−1, +1}. On the other hand, the vacuum contribution is conveniently computed by rewriting it as a standard d-dimensional Euclidean integral Seen as a function of a complex α, this integral is analytic, with a branch cut for α ∈ Re − . Therefore, its value for α imaginary can be obtained by analytic continuation of the known expression for α ∈ Re + . We simply find whereμ 2 ≡ 4πµ 2 e −γ and γ is the Euler-constant. Using similar techniques, we now would like to evaluate the (0-leg) sunset sum-integral where momentum and color conservation imply respectively Q+K +L = 0 and κ+λ+τ = 0, see Fig. 1. 3 We consider the case where the square masses α, β and γ are purely imaginary.
In fact, we restrict to those cases that are relevant for the GZ framework, where the square masses are either 0 or ±im 2 . More precisely, it can be shown that the relevant scalar sunset sum-integrals that appear in the GZ framework are S κλτ ααα , S κλτ αα(−α) , S κλτ α00 , with α = ±im 2 , together of course with the corresponding permutations of masses [36].
Using the spectral representation (6) in Eq. (16), we find with a double Matsubara sum. Standard techniques for the evaluation of Matsubara sums together with color conservation κ + λ + τ = 0, lead then to where, in going from the first to the second line, we have used the well known identity n x n y = (1 + n x + n y )n x+y . Finally, by making use of Eq. (11), we arrive at the thermal splitting of the scalar sunset sum-integral: S κλτ αβγ = S αβγ (0n) + S κλτ αβγ (1n) + S κλτ αβγ (2n). As in the previous example, S αβγ (0n) denotes the pure vacuum contribution (no thermal factor), depending neither on the temperature nor on the background, while and S κλτ αβγ (2n) = cyclic q 0 ;q σ κ α (q 0 ; q) where cyclic stands for the cyclic permutations of the pairs (α, κ), (β, λ) and (γ, τ ) and σ κ α (q 0 ; q) was defined in the previous section. We note that the thermal splitting considered here assumes that the denominator l 0 + k 0 + q 0 never vanishes. We show in Appendix A that this is indeed so in those cases that are relevant for the GZ framework. More generic cases may require a regularization of the denominator but we shall not consider them here.

B. Contribution with one thermal factor
Integration over the frequencies leads this time to where we note that the dependence on σ α has dropped in the last line, which explains a posteriori why we did not include it in our notation forĨ βγ (ε q,α ; q). We show in Appendix B that this quantity does not depend on q either. It follows that S κλτ αβγ (1n) = cyclic J κ α (1n)I α βγ (0n), with The case I 0 βγ (0n) can be evaluated immediately as As for the general case I α βγ (0n), we show in Appendix B that it can be obtained from the analytic continuation of the vacuum Euclidean integral We find with α − ≡ α − 0 + and More details on the continuation are given in Appendix B.

C. Vacuum contribution (no thermal factor)
The vacuum contribution can be written as a standard d-dimensional Euclidean integral This integral is known analytically for positive square masses α, β and γ [42]. One strategy to obtain the corresponding integral for imaginary square masses is to perform an analytic continuation of the result of [42] with respect to the masses. This is an efficient strategy in the case where the non-vanishing masses are all equal.
Let us take for instance S α00 (0n). Seen as function of a complex α, it is analytic with a branch cut for α ∈ Re − . Because the expression for α ∈ Re + , has branch cuts only for α ∈ Re − , it can be immediately used to represent S α00 (0n) in the case where α is purely imaginary. Using a similar argument, we obtain The evaluation of S αα(−α) (0n) is trickier because it involves the analytic continuation of a function of two complex variables. Although this can be done in principle, we here chose a more direct evaluation by adapting the technique in Ref. [42] to the case of imaginary square masses. This technique is based on the derivation of a differential equation satisfied by S αβγ (0n). Although the derivation of the differential equation is not affected by the presence of imaginary square masses, we reproduce it here for completeness and because we shall use it for other purposes later.
Let us write the scalar sunset vacuum integral symbolically as with L = −Q − K. Expanding the identities we find Then, writing Q 2 = Q 2 + α − α, as well as we arrive at Using the first relation in order to replace {G 2 α G β G γ } in the second, we find Finally, using that In the case where γ = 0, we can use to arrive at the differential equation In Appendix C, we use this equation to obtain S ααγ (0n) in some appropriate range of values for γ. We follow the approach of Ref. [42] by carefully adapting it to the case of imaginary square masses. We find in particular We mention that the vacuum integrals S +++ (0n), S +00 (0n) and S ++− (0n) have also been used in [37,39] but the individual results are not quoted, only their final combination in the vacuum GZ horizon condition at two-loop order. For similar integrals involving both real and imaginary square masses, see [38].

IV. MASS DERIVATIVES
The various sunset diagrams that appear in the GZ framework lead also to mass derivatives of the scalar sunset, in the limit where the corresponding mass is taken to zero. In fact, what appear are the limits [36] ∆S κλτ where, as already introduced above, the squaring of the mass indices corresponds to the doubling of the associated propagators, or more generally to taking minus the derivative with respect to the associated square mass. The relevant cases for the GZ framework are ∆S κλτ 0 2 γ(−γ) , ∆S κλτ 0 2 0γ , ∆S κλτ 0 2 0 2 γ , with γ = ±im 2 , together with the corresponding permutations of the masses [36]. We also mention that ∆S 0 2 ,αβ is nothing but the functionT defined in [45], for zero external momentum but generalized to the case of finite temperature.
It is easily seen that, even though S κλτ α 2 βγ , S κλτ α 2 β 2 γ , J κ α 2 and J λ β 2 are singular in the limits α → 0 or β → 0, the combinations of sum-integrals in Eqs. (53) and (54) admit regular limits. To verify this, let us first note that the potential singularities originate either from the vacuum pieces which do not depend on the color weights (κ, λ, τ ), or from the thermal pieces in the case where the weights are equal to zero. Therefore, one can safely ignore the weights in order to check the regularity of the above limits. For instance, in the limit α → 0, the sum-integral S α 2 βγ is dominated by the Q → 0 region and behaves consequently as The divergent behavior in the RHS is precisely what is subtracted in Eq. (53) to ensure that the limit is regular. 4 To understand the singular structure of S α 2 β 2 γ , we first write it identically as The term with K · Q in the second line leads to regular contributions in the limit α → 0 and β → 0, while the term with K 2 (resp. Q 2 ) leads to singular contributions as α → 0 (resp. β → 0), controlled by the Q → 0 (resp. K → 0) region of the integral. We find These are precisely the terms that are subtracted in Eq. (54) to obtain a regular limit.
The regular limits ∆S κλτ 0 2 βγ and ∆S κλτ 0 2 0 2 γ admit thermal splittings that one derives from the corresponding splitting (51) of the scalar sunset, and which we now discuss.
A. Thermal splitting From Eqs. (51) and (53), we find the thermal splitting Similarly, from Eqs. (51) and (54), we find the thermal splitting with where we have used that J 0 (0n) = 0. In what follows, we shall also make extensive use of the property J 0 2 (0n) = 0, valid in dimensional regularization []. This property might be more difficult to grasp than the previous one because J α 2 (0n) diverges in the limit α → 0. However this just means that the function J α 2 (0n), although defined for α = 0, is not continuous at α = 0. Then, we shall always make sure that when the property J 0 2 (0n) = 0 is used, it corresponds to J α 2 (0n) being evaluated for α = 0 and not to a limit being taken. We mention finally that the results to be presented below can be obtained without ever using J 0 2 (0n) = 0 although the calculations are lengthier.
C. Contributions with one thermal factor ∆S κλτ 0 2 βγ (1n): The contribution with one thermal factor to ∆S κλτ 0 2 βγ can be rewritten as Owing to Eq. (29), it seems that the contribution in the first round bracket can be neglected in the limit α → 0. This turns out to be true, although one needs to pay a little bit of attention, since, at finite temperature and in the case where κ = 0, J κ α 2 (1n) diverges in the same limit. Fortunately, the divergence goes as α −1/2 , which is not enough to compensate the vanishing of the round bracket ∼ α. We find eventually that with where the contribution within brackets is regular in the limit α → 0 and we have used J 0 2 (0n) = 0 in the last step.
The first quantity can be computed using similar tricks as for I 0 βγ (0n) in the previous section, namely where we have used that together with T =0 K 1 = 0 in dimensional regularization.
As for ∆I β 0 2 γ (0n), we can proceed in many different ways, either by acting with −∂/∂α on the previously determined expression for I β αγ (0n), followed by the α → 0 limit after appropriate subtraction of the α → 0 singular part, or by computing the appropriate subtracted Euclidean integral and analytically continuing it from K 2 > 0 to K 2 = −β imaginary. Here we proceed with this second strategy but instead of continuing the explicit expression of the integral, we continue the corresponding differential equation, with the advantage that ∆I β 0 2 γ (0n) will be expressed in terms of already computed quantities.
To derive the differential equation, we basically consider the same equations as in Eqs. (37) but with the propagator G β (K) and the integral T =0 K missing. It is easily seen that one can follow the steps below Eqs. (37) by removing the factors G β (K) (those terms that did not have such a factor need to be discarded) and to replace the explicit occurrences of β by −K 2 . It follows that This identity is valid for α = 0, in which case, we obtain After continuation, we find eventually ∆S κλτ 0 2 0 2 γ (1n): Similarly, the contribution with one thermal factor to ∆S κλτ 0 2 0 2 γ can be rewritten as Using Eq. (29), this rewrites We note that Using Eqs. (74) and (79) and after subtracting the β → 0 and γ → 0 singular parts, according to Eqs. (84) and (85) respectively, we find eventually D. Contributions with two thermal factors ∆S κλτ 0 2 βγ (2n): The contribution with two thermal factors to ∆S κλτ 0 2 βγ can be rewritten as The terms with the α-derivative acting on the thermal factor can be treated using an integration by parts after noticing that df (ε q,α )/dα = df (ε q,α )/dq 2 = (df (ε q,α )/dq)/(2q). The boundary term vanishes both for q → ∞ (due to the thermal factor). The boundary at (note that α → 0 is taken only after the limit q → 0). We obtain where ∂/∂q denotes the partial derivative at ε q,α fixed. Using ∂ ∂q ln αβ; γ = 2(k + q) γ − α − β − 2σ α σ β ε q,α ε k,β + 2qk as well as we find eventually Note that the first two integrals remain safe when κ = 0 provided we first sum over σ α .

V. CONCLUSIONS
In this work, we have evaluated the scalar sunset diagrams with imaginary square masses that appear in the two-loop background potential in the GZ type model with a background gauge invariance from Ref. [34]. This also involves some mass derivatives of the scalar sunset in the limit where the corresponding mass is taken to zero. In fact, what appear are not mass derivatives by themselves, but specific combinations with tadpole integrals and their mass derivatives that admit regular limits in the zero mass limit. The evaluated cases include three scalar sunsets and three mass derivative combinations. In each case the square masses are either 0 or ±im 2 .
Through thermal splitting we have decomposed the sum-integrals into contributions with 0, 1 and 2 thermal factors. For the terms with 0 thermal factors, the vacuum contribution, we obtained the integral by analytic continuation from the results for real masses from Ref. [42], for the cases where the non-vanishing masses were equal. In the other cases, i.e. the cases with square masses of opposite sign, we have made a direct evaluation adapting the technique of Ref. [42] to imaginary square masses.
Instead of considering only the scalar sunset sum-integrals that appear in the GZ framework, one could make a broader study of all scalar sunset diagrams with purely imaginary masses. However, when using thermal splitting this requires some regularization of the denominator for the contributions with 1 or 2 thermal factors, in order to avoid singularities.
This problem is particular for imaginary masses: in the case of real masses one simply adds a regulator to the denominator in the form of an infinitesimal imaginary number. In some cases, like the cases considered here, the imaginary masses themselves work as a regulator, making it impossible for the denominator to vanish, but this is not true in general. Even when we limit ourselves to cases where the square masses are either 0 or ±im 2 , there are examples where the denominator can vanish, e.g. when one square mass is 0 and the other two square masses are im 2 . In principle, it is possible to find a consistent regularization for each case but one should investigate how this affects the subsequent steps of the calculation.
Since this lies beyond the scope of the GZ application that we are pursuing, we leave this question for a future study.
The sunset diagrams that have been calculated in this work make up a substantial part of the calculation of the two-loop background potential in the GZ framework. We are currently evaluating the full two-loop potential [36] in the presence of a temporal background in order to study the deconfinement transition in YM theory using this framework, as well as the interplay between the Polyakov loop and the Gribov parameter. It will be interesting to compare these results with the one-loop results in the same model from Ref. [34], as well as with the two-loop studies in the CF model at finite temperature from Refs. [16,[22][23][24][25]. Théorique at Ecole Polytechnique, where this work was developed.
For real masses, one needs to add an imaginary regulator i0 + to the denominator to avoid divergences. For imaginary square masses, the discussion is more intricate. As we now argue, however, in all cases of interest for the GZ framework, a regulator is not necessary.
Since m > 0, we can solve the second equation as qk cos θ = −(q 2 + k 2 ) and plug it back into the first condition to arrive at which has again no solution if m > 0.
Appendix B: Evaluation ofĨ βγ (ε q,α ; q) Let us consider the vacuum Euclidean integral with imaginary square masses β and γ. In order to make contact withĨ βγ (ε q,α ; q), let us evaluate the frequency integral in (B1) using the residue theorem. To this purpose, we write it as (B5) In the last step, we have written q 2 4 as −ε 2 q,−Q 2 to emphasize the fact this last integral is similar to the one definingĨ βγ (ε q,α ; q) in Eq. (27), with α replaced by −Q 2 . More precisely, if we introduce the function seen now as a function of a complex Q 2 for a fixed q, we have both F q (Q 2 ≥ 0) = I βγ (0n)(Q 2 ≥ 0) and F q (Q 2 = −α ∈ iR) =Ĩ βγ (ε q,α ; q).

Analytic continuation
To turn this observation into a practical way to determineĨ βγ (ε q,α ; q), we note first that ifĨ βγ (ε q,α ; q) makes sense for a given value of α, it makes sense for any other value close to it. In particular, we can writẽ Second, it is easily seen that F q (Q 2 ) is analytic in the semi-plane Re Q 2 > 0. Indeed, the potential singularities are restricted to the region defined by the condition (ε k,β + ε l,γ ) 2 − ε 2 q,−Q 2 = 0, which corresponds to and whose real part obeys From these considerations, it follows that, if we know explicitly a function G(Q 2 ) which is analytic over an open connected subset Ω of Re Q 2 > 0 containing both the Q 2 = −α + 0 + and the Q 2 > 0 axis, and which agrees with I βγ (0n)(Q 2 ≥ 0) along this axis, then G(Q 2 ) = F q (Q 2 ) over Ω. In particular,Ĩ βγ (ε q,α ; q) can be obtained as G(Q 2 = −α + 0 + ).

Explicit expression
It remains to construct explicit examples of G(Q 2 ). This is easily done by evaluating I βγ (0n)(Q 2 ≥ 0) using the Feynman trick and by extending the result to complex values of Q 2 . Of course, as long as Q 2 ≥ 0, there are many equivalent forms of G(Q 2 ) that one can write. In order to determineĨ βγ (ε q,α ; q), we should only use those forms that obey the above mentioned analytic properties.
The Feynman trick allows us to write To perform the integral over x, it is convenient to write If we want to split the logarithm of the product of the two factors in Eq. (B11), we need to check that the sum of the arguments of the factors lies between −π and π. Since, the product of the factors never crosses the branch cut, it is enough to work at x = 0 and Q 2 = 0. One finds which lies between −π/2 and π/2. We can then split the logarithm and compute the xintegral to obtain, after some trivial simplifications, It will be convenient to rewrite this Euclidean expression as where we have changed the scale under the logarithms to the price of changing lnμ 2 /Q 2 by Let us now analyse the singularities of G (1) (Q 2 ). We mention that we are not after the precise determination of the singularities. Rather we want to check that they comply with the above mentioned requirements, allowing one to extract I βγ (ε q,α ; q) as G (1) (Q 2 = −α + 0 + ).
First, there is a branch cut along the negative real axis, originating from lnμ 2 /Q 2 . A second branch cut originates from R(−Q 2 , β, γ) which contains a square root. More precisely, this branch cut corresponds to The determination of theV (j) 's proceeds recursively: one first determinesV (−2) by integrating the corresponding differential equation with the explicit expression (C3) forḡ (−2) .
KnowingV (−2) , one can then determineḡ (−1) from (C4) and repeat the procedure, until all theV (j) 's have been determined. We mention that the integration of each differential equation gives eachV (j) (m 2 , c 2 ), in terms of a boundary valueV (j) (m 2 , c 2 0 ). There seems to be a circular reasoning a priori. We see below how this problem is avoided.
Following [42], we note that ∂ ∂c 2 It follows that the differential equation can be rewritten as The benefit of this rewriting is two-fold. First, it can be integrated to provide an expression forV (j) (m 2 , c 2 ) in terms of an integral involvingḡ (j) (m 2 , c 2 ) and a boundary valuē which provides a one-dimensional integral representation forV (j) (m 2 , c 2 ).