Circular Wilson loop in N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N} $$\end{document} = 2* super Yang-Mills theory at two loops and localization

We present a two-loop calculation of the supersymmetric circular Wilson loop in the N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N} $$\end{document} = 2* super Yang-Mills theory on the four-sphere. We develop an efficient framework for computing contributing Feynman graphs that relies on using the embedding coordinates combined with the Mellin-Barnes techniques for propagator-like integrals on the sphere. Our results exactly match predictions of supersymmetric localization providing a nontrivial consistency check for the latter in non-conformal settings.


Introduction
The supersymmetric localization offers a very powerful tool to study supersymmetric gauge theories both in perturbative and nonperturbative regimes, see, e.g., the review [1]. It allows one to compute various observables in theories defined on the four-sphere in terms of interacting matrix models. One notable example of this is the N = 2 * theory, which JHEP04(2021)089 arises as a massive deformation of the maximally supersymmetric Yang-Mills theory in four dimensions preserving N = 2 supersymmetry. The Lagrangian of the latter in conformallyflat spaces is uniquely fixed by its extended supersymmetry [2,3].
To date, the localization has been applied to calculations of partition functions, circular Wilson loops and certain correlators in various supersymmetric theories [1]. The corresponding matrix model predictions have been successfully checked at weak coupling against explicit perturbative analyses in superconformal [4][5][6][7][8][9][10][11][12][13][14] and in massless non-conformal supersymmetric gauge theories [15]. In all of these circumstances, perturbation theory was developed in flat four-dimensional Euclidean space and observables were shown to be the same (at least to the order considered) as the ones on the four-sphere.
The use of flat-space analysis is invalidated, however, when the theory contains an additional scale, like a deformation mass-scale in the N = 2 * theory. In this case, conformal symmetry is explicitly broken and the theory possesses different properties at short and long distances. In particular, observables exhibit dependence on a dimensionless parameter determined by the product of the mass and the radius of the four-sphere. To check the localization formulas on the sphere, we have to therefore generalize conventional Feynman integral technique to theories on curved backgrounds. This is what we are set to do in the present work. Namely, we develop an appropriate formalism for evaluating Feynman diagrams on the sphere and apply it to computing the circular Wilson loop in the N = 2 * theory to two-loop order.
A convenient way to organize the calculation in the N = 2 * theory is to decompose all contributing graphs into those present in the N = 4 theory and remaining ones, as it was done in a related but different context in refs. [4,13]. The contribution of the former diagrams is known explicitly whereas the latter define the difference between the circular Wilson loops in the two theories. To two-loop perturbative order, the difference only arises from graphs with gluon and scalar propagators dressed by one-loop corrections due to massive particles circulating in the virtual loop. This allows us to focus on propagator-like integrals containing massless external lines and massive loops.
Using the formalism of embedding coordinates [16][17][18][19], it is straightforward to obtain a representation for two-loop corrections to the circular Wilson loop in terms of scalar integrals built from propagators of massless and massive particles. However, due to rather involved form of the latter, evaluating them by brute force becomes prohibitively complicated. We demonstrate in this paper, that all calculations can significantly be simplified by employing the Mellin-Barnes transformation for propagators involved. This yields a very compact representation for contributing graphs as triple Mellin-Barnes integrals. Expanding them at small and large values of the mass, we obtain expressions for the circular Wilson loop to two loops. We verify below that it agrees with the prediction from the localization and, at the same time, it differs from the analogous expression for the circular Wilson loop in flat space.
The rest of the paper is organized as follows. In section 2, we present the N = 2 * super Yang-Mills theory on the four-sphere. We then introduce the embedding coordinates, which prove to be of extreme value in our calculations, and provide expressions for free propagators of fields belonging to the vector and hypermultiplets. In section 3, we JHEP04(2021)089 define the circular Wilson loop in N = 2 * theory. We compute it to two loops in flat space and compare the result with the analogous expression on the sphere predicted by the supersymetric localization. In section 4, we use the embedding formalism to express relevant two-loop corrections in terms of scalar integrals on the sphere. The calculation of these integrals is discussed in section 5. We outline our results and check them against the localization prediction both numerically, for finite mass, and analytically, in the regime of small and large mass. In both cases we observe a perfect agreement. Concluding remarks are presented in section 6. Calculational details, which are too bulky for the main text, are summarized in six appendices.

N = * super Yang-Mills theory on the sphere
In this section, we define the N = 2 * super Yang-Mills theory on the sphere. It arises as a mass deformation of the maximally supersymmetric N = 4 super Yang-Mills theory in four-dimensional space. The deformation breaks the conformal symmetry of the latter but preserves the residual N = 2 Poincaré supersymmetry. The Lagrangian of this theory on the four-sphere has been constructed in refs. [2,3] and we review it below.

Mass deformation
The maximally supersymmetric Yang-Mills theory can be obtained by compactifying the N = 1 super Yang-Mills in ten space-time dimensions down to four spatial dimensions. It describes gauge fields A µ , six real scalars X I (with I = 1, . . . , 6) and four gauginos λ A andλ A (with A = 1, . . . , 4) all belonging to the adjoint representation of the SU(N ) gauge group. 1 The action of the N = 4 on the four-sphere S 4 looks as where L 0 contains minimally-coupled kinetic terms of the fields and L int describes Yukawa interaction and quartic self-interaction of scalars, Here the scalar fields are grouped into skew-symmetric matrices φ AB = Σ I,AB X I / √ 2 and φ AB ≡ ABCD φ CD /2 (with A, B = 1, . . . , 4) with the expansion coefficients Σ I ,AB being the 1 The gauge fields and scalars arise as components of the gauge field AM = (Aµ, XI ) in D = 10 dimensional space-time with Minkowski signature. Rewriting the action of N = 4 theory in terms of these fields, one finds that the last component of AM has a different sign in front of the kinetic term. This implies that the action is no longer real for conventionally defined involution of fields. This requires either compexification of the path integral or introduction of a different involution such as symplectic Majorana condition [20] or just giving up hermiticity altogether.

JHEP04(2021)089
chiral blocks of Dirac matrices in six dimensions. The gauginos λ A α andλ α A transform under the two subgroups of the rotation group in four dimensions Spin(4) = SU(2) × SU (2). The covariant derivative for gauginos looks as D αβ λ β , is the spinor covariant derivative on the sphere (see eq. (C.2) below). Similar definition holds for derivatives acting on other fields, see, e.g., eq. (B.2). Finally, the last term in the expression for L 0 describes the conformal coupling of the scalars to the scalar curvature of the four-sphere with radius R. Upon the mass deformation, the massless N = 4 multiplet is decomposed into the N = 2 vector multiplet (A, λ, Φ) consisting of and the massive hypermultiplet (ψ, Z) consisting of where i = 1, 2. It is straightforward to rewrite the action (2.1) in terms of these fields. To save space, we refrain from doing it. The circular Wilson loop defined in eq. (3.1) below only depends on the scalar field Φ 1 . As a consequence, the scalar field Φ 2 does not appear in our calculation and we do not display it in all formulas that follow. The mass deformation of the N = 4 theory is driven by the following term whereZ i andψ α i are conjugated fields to their unbarred counterparts and Φ = Φ 1 . 2 In flat space, i.e., for R → ∞, it is generated by the superpotential [21], It corresponds to giving a vacuum expectation value to Z's via the mechanism of spontaneous symmetry breaking. Then, upon the shift by Z = O(m), the four-scalar interaction in (2.2) yields the first term in (2.5) and the Yukawa coupling induces a mass term for the gauginos. The last term in the right-hand side of (2.5) is required by supersymmetry [2,3]. It breaks, however, the reflection positivity [22]. Observables computed in such non-unitary setup do not, in general, possess standard reality properties. As we show below, this is not the case for the circular Wilson loop. The resulting action of the N = 2 * model on the sphere takes the form

JHEP04(2021)089
where the first two terms inside the brackets are given by eq. (2.2) with the fields redefined according to eqs. (2.3) and (2.4).
To reveal the particle mass spectrum in this theory, we examine the action (2.7) at zero value of the gauge coupling. We get where the first and the following two lines encode the field content of the vector and hypermultiplets, respectively. The kinetic terms for fields involve the covariant derivatives on the sphere defined in (B.2) and (C.2). Here we introduced two real (pseudo)scalar fields According to (2.4), they originate from the progenitor N = 4 fields (X 1 , X 4 ) and (X 2 , X 5 ), respectively. It follows from (2.8) that A i and B i acquire complex "masses', respectively. Different, though real, masses for (pseudo)scalars of the Wess-Zumino multiplet in the anti-de Sitter space were known to be enforced by supersymmetry [23]. By construction, the N = 2 * theory coincides with the N = 4 one for m → 0. In the opposite limit, i.e., m → ∞, the hypermultiplet becomes infinitely heavy and we recover the N = 2 super Yang-Mills.
To simplify the analysis that follows, we set the radius of the four-sphere to be R = 1. When needed, the dependence on R can be easily restored by rescaling the mass and distances.

Embedding coordinates
The action of the N = 2 * theory (2.7) is defined on the four-dimensional sphere. For the reason that will become clear in the next section, we have to rather generalize the definition (2.8) to the D-dimensional sphere.
Instead of using local coordinates on the sphere (see appendix A), it is advantageous to embed S D in a flat Euclidean space R D+1 . Introducing the embedding coordinates X M (with M = 0, . . . , D), the D-dimensional sphere of radius R = 1 is defined as

JHEP04(2021)089
Employing a stereographic projection we can map the sphere to a hyperplane X M = (0, x µ ) with µ = 1, . . . , D and x 2 = x 2 µ . An important property of this transformation is that it is conformal in R D+1 . Therefore, had the N = 2 * theory (2.8) possessed conformal symmetry, we could apply (2.12) to obtain the same theory in the flat space-time R D .
Since the conformal symmetry of the N = 2 * theory (2.7) is broken for m = 0, we are bound to study it on the sphere. Notice that at zero coupling, the Lagrangian of the vector multiplet, given by the first line in (2.8), does in fact enjoy the conformal symmetry. We exploit this fact in the next subsection to obtain a compact representation for free propagators of massless fields in the N = 2 * model.
A lot of activity has been devoted in the past to study field theories in the de Sitter space. This space can be obtained from (2.11) by a Wick rotation, X 0 → iX 0 , where X M are the embedding coordinates in R 1,D with the Minkowskian signature η M N = diag(+, −, . . . , −). A convenient choice of local coordinates in the de Sitter space, widely used in applications to cosmology, looks as where i = 1, . . . , D − 1 and z µ (with µ = 0, 1, . . . , D − 1) are the Minkowskian coordinates, The inverse transformation takes the form (2.15) We show in appendices B and C that calculations on the sphere can sometimes be simplified by Wick rotating to de Sitter space and using the coordinates (2.14) instead.

Propagators of vector multiplet
As a first step toward computing the circular Wilson loop in the N = 2 * theory, we derive in the next two subsections expressions for free propagators of fields entering (2.8) We start with the massless sector of the theory. As was mentioned in the previous subsection, free propagators of massless fields on the sphere S D can be obtained from analogous expressions in flat space by a conformal transformation (2.12). For the scalar field Φ belonging to the vector multiplet we have [18], (2.16) with the notation X 12 = X 1 − X 2 used throughout the paper.

JHEP04(2021)089
Let us turn to the gauge propagator. The gauge field on the sphere is related by a coordinate transformation to its analogue in R D+1 (2.17) where x µ are the local coordinates on the sphere S D and X M are the embedding coordinates in R D+1 . In the Feynman gauge, the propagator D M N A (X 12 ) = A M (X 1 )A N (X 2 ) obeys the same equation as the scalar propagator and, therefore, it reads [16,19] where η M N = diag(+, . . . , +) is the metric tensor in R D+1 .
To determine the massless gaugino propagators, it is convenient to combine two Weyl fermions into a Majorana one, Λ = (λ α ,λ α ), and to work with four-components spinors and Dirac matrices γ M (with M = 0, . . . , D) in R D+1 . As before, the fermion propagator on the sphere is related to that on R D+1 by [16,19] where the rotation matrices U i = U (X i ) satisfy the following defining relations with the usual notation for X ≡ γ M X M , x µ being the local coordinates on S D and the corresponding vielbeins e a µ (x) (with a = 0, . . . , D − 1). The explicit form of the matrix U (X) can be found in the appendix C, although is not important for our purposes.

Propagators of massive hypermultiplet
Moving on to the massive hypermultiplet, we can not rely on the conformal symmetry anymore. The details of calculations of the gaugino and scalar propagators entering the last two lines of (2.8) can be found in the appendices B and C, respectively.
We start with the scalars A i and B i (with i = 1, 2), which differ in their complex "masses" (2.10). Their propagators take the following form in the embedding coordinates where a new notation is introduced for to emphasize the relation between the two via the reflection transformation m → −m. We verify that for m → 0 the functions S ± (X 2 12 ) coincide with the massless propagator (2.16). As in the massles case, it is convenient to introduce a Majorana fermion Ψ i = (ψ i α ,ψ α ,i ). Referring to the appendix C for a thorough exposition of the calculations JHEP04(2021)089 involved, we only quote here the result for the free propagator of the Majorana fermion of the massive hypermultiplet It involves the same rotation matrices as in the massless case (2.19) but contains an additional Dirac matrix structure, X 12 X 2 with X 12 ≡ X 1 + X 2 . The coefficient functions A(X 2 12 ) and B(X 2 12 ) are given in terms of S ± introduced earlier, We verify that for m → 0 the function A(X 2 12 ) vanishes and the expression for D Ψ (X 1 , X 2 ) coincides with the massless propagator (2.19). The properties of the functions A(X 2 12 ) and B(X 2 12 ) are discussed in the appendix D. One can check that the scalar and gaugino propagators, eqs. (2.21) and (2.23), verify supersymmetric Ward identities that follow from δ Ψ i (X 1 )Z j (X 2 ) = 0 with δ being the supersymmetric variation which leaves the Lagrangian (2.8) invariant. This explains why the same functions S ± (X 2 12 ) enter (2.21) and (2.23).

Circular Wilson loop
In this work, we study the circular supersymmetric Wilson loop in N = 2 * theory on the four-sphere. It is defined as [24] where the gauge field A µ (x) and the real scalar Φ ≡ Φ 1 from the vector multiplet (2.3) are integrated along a great circle C on the four-sphere of radius R = 1. Here |ẋ(s)| = (g µνẋ µ (s)ẋ ν (s)) 1/2 withẋ µ (s) = ∂ s x µ (s) and g µν (x) being the metric tensor on the sphere. Using (2.17), we obtain an equivalent representation for the circular Wilson loop in the embedding coordinates where |Ẋ(s)| = (Ẋ 2 M (s)) 1/2 and the great circle on the sphere C can conveniently be parametrized as with an affine parameter s ranging in the interval 0 ≤ s < 2π.

JHEP04(2021)089
To lowest order in the coupling constant, we apply (2.16) and (2.18) for D = 4 to get from (3.2) is the quadratic Casimir in the fundamental representation of the SU(N ). Here the first term in the numerator of the integrand comes from the gauge field and the second one from the scalar. The one-loop correction to (3.4) is not sensitive to the hypermultiplet and, therefore, it is the same both in the N = 2 * and N = 4 theories. The massive hypermultiplet (2.4) contributes to (3.2) starting at order O(g 4 ) and, as a consequence, the circular Wilson loops in the two theories differ already at two loops. In the N = 4 theory, the circular Wilson loop can be found exactly at arbitrary value of the gauge coupling g 2 and finite N . At weak coupling, it was demonstrated by an explicit calculation in [25] (see also [26][27][28]) that W receives a nonzero contribution from ladder diagrams only and it can be expressed as a matrix integral. In the planar limit, the planar circular Wilson loop reads where I 1 is the modified Bessel function and λ = g 2 N is the 't Hooft coupling. The relation (3.5) was proved in ref. [2] using supersymmetric localization. It was also generalized there to the case of the massless superconformal N = 2 Yang-Mills theory with fundamental hypermultiplets and the N = 2 * theory on the sphere. The former result was tested in refs. [4,5] against explicit three-loop calculation of the Wilson loop in flat space. However, up to now not a single study was dedicated to testing the localization predictions in the N = 2 * theory. The reason for this is quite simple: in virtue of conformal symmetry, all of perturbative checks without exceptions were performed in flat space. In N = 2 * theory this simplification does not hold and one has to perform fullfledged calculations on the sphere with massive propagators. As we discussed in the previous section, the latter are quite involved on their own.

Prediction from localization
The circular Wilson loop in the N = 2 * theory is a finite function of the gauge coupling and the dimensionless parameter mR. Although we put R = 1 for simplicity, the dependence on R can be easily restored by replacing m → mR. The conformal symmetry of N = 2 * theory is explicitly broken by the mass deformation. As a consequence, in distinction to the N = 4 theory, one expects that the circular Wilson loop should admit different forms in the N = 2 * case on the sphere and in the flat space. We will show this explicitly below.
As was alluded to above, the circular Wilson loops are identical in the N = 4 and N = 2 * theories at one loop. The difference emerges only at two loops due to massive hypermultiplets circulating in virtual loops. Following refs. [4,13], we find it convenient to JHEP04(2021)089 consider the difference of the Wilson loops in the two models, W N =2 − W N =4 . This allows us to avoid computation of diagrams with the massless N = 2 vector multiplet in internal subgraphs and, thus, reduce significantly the number of contributing Feynman integrals.
The supersymmetric localization predicts the following result for the difference of the circular Wilson loops at weak coupling [2] where we inserted the subscript S 4 to indicate that the gauge theories are defined on the four-sphere. Writing down the color factor in the right-hand side of (3.6) in terms of the quadratic Casimir C F = (N 2 −1)/(2N ), we made use of the known property of non-Abelian exponentiation for Wilson loops in gauge theories [29,30].
The function f (m) carries the dependence on the mass parameter. 3 Since it is independent of the rank of the gauge group N , it can be derived from the localization formulas for the SU(2) gauge group. This gives where ψ(x) = (log Γ(x)) is the digamma function and γ is the Euler-Mascheroni constant.
For future reference, we present its expansion in two different asymptotic regions: The logarithm in the last relation is related to the beta-function in pure N = 2 super Yang-Mills theory. Indeed, for m → ∞, the hypermultiplet becomes infinitely heavy and W N =2 should match the Wilson loop in the pure N = 2 super Yang-Mills theory endowed with an ultraviolet cut-off m. The coupling constant in the latter theory depends on m and this dependence is driven by the known one-loop-exact beta-function. The logarithmically enhanced term in the expression for has to be proportional to the beta function in order for W N =2 to remain finite as m → ∞.
Our goal in this work is to verify (3.6) by computing ∆W S 4 to order O(g 4 ). As we just explained, the difference comes from graphs with the massive hypermultiplet circulating in loops. To two loop order, the only graphs we have to account for are those describing one-loop correction to propagators of the gauge field and massless scalar belonging to the vector multiplet, D   .6) is Here the last term in the right-hand side subtracts the contribution at m = 0, so that this difference vanishes for m → 0.
Each individual graph shown in figure 1 depends on the mass parameter m and develops an ultraviolet divergence. The divergences cancel, however, in the right-hand side of eq. (3.11), after we subtract the contribution of the same graphs with m = 0. Still, to define the contribution of individual diagrams, we have to introduce a regularization. For this purpose we employ conventional dimensional regularization and perform all calculations on the sphere S D with D = 4 − 2 . This could have been avoided had we repacked individual contributions differently. However, we found our way of performing calculation more advantageous.

Circular Wilson loop in flat space
In this subsection, we compute the difference W N =2 − W N =4 in flat space and compare it with the analogous expression on the sphere predicted by the localization (3.6). As was explained at the beginning of this section, we anticipate the two expressions to be different due to conformal symmetry breaking in the N = 2 * theory.
To lowest order in gauge coupling, the difference ∆W R 4 is given by the same relation (3.11) with X M = (x µ , 0) and x µ being local coordinates in R 4 . The most concise and efficient way to perform calculations in flat space is to take advantage of translational invariance and pass to the reciprocal momentum space via the Fourier transform

JHEP04(2021)089
where we factored out the dependence on the gauge coupling and the color factor. A similar relation holds for the scalar propagator D (1) m (x 12 ). Thus, there are two polarization operators we have to evaluate, one for the gauge field, Π µν G (p), and another for the massless scalar, Π Φ (p). Let us consider them in turn.
The Feynman diagrams contributing to the gluon polarization operator at one loop are depicted in the top row of figure 1. In dimensional regularization with D = 4 − 2 , the contribution from the graphs in G 1 , G 2 and G 3 reads, respectively, in the momentum representation (3.13) The relative multiplicative factors are as follows: (−1) in front of the first term comes from the Dirac-Fermi statistics, the factor of 2 in front of the second term is from two pairs of real (pseudo)scalars A i /B i in the hypermultiplet, and, finally, 4 is from the same two pairs of scalars inside the tadpole. Due to the gauge invariance, the polarization operator is transverse, (3.14) Contracting both sides with g µν , we find the integral representation for Π m (p 2 ), . (3.15) Integration over the large loop momentum yields a pole 1/(D − 4) but it cancels in the difference π m (p 2 ) − π 0 (p 2 ). The polarization operator for the massless scalar receives contributions from the Feynman diagrams in the second row of figure 1. It reads for the Yukawa coupling (S 1 ), three-scalar interaction (S 2 ) and tadpole (S 3 ), respectively.
Here the multiplicities of the hypermultiplet fields get altered by the strengths of couplings as governed by the Lagrangian (2.7) in the limit R → ∞. Adding up all contributions, we find with π m (p 2 ) given in eq. (3.15), in agreement with expectations based on supersymmetry. Combining together the above relations, we obtain from (3.11) JHEP04(2021)089 where the notation was introduced for Here in the second relation, we replaced π m (p 2 ) with its integral representation (3.15) and performed the loop momentum integration using Feynman parametrization. To proceed further, we use the Mellin-Barnes representation for the first term inside the brackets in (3.19) where the integration contour runs parallel to the imaginary axis and separates poles generated by the two Γ−functions in the integrand. Substituting (3.19) into (3.18) and performing sequentially the integrations over the momentum and parameters s 1 , s 2 , we finally obtain where Re j = −δ and 0 < δ < 1. We can use this relation to find asymptotic behavior of ∆W R 4 at small and large values of m. At large m, we move the integration contour to the left and pick up residues at poles j = −1, −3/2, −5/2, . . .
wherem = m e γ /2. At small m, we move the integration contour to the right and observe that the integrand has third order poles at j = 0, 1, . . . . Evaluating the residue at these poles, we obtain

JHEP04(2021)089
It is interesting to note that subleading, i.e., power suppressed, corrections to the circular Wilson loop have a form different on the sphere than in the flat space: for ∆W S 4 and ∆W R 4 they run, respectively, in even and odd powers of 1/m. The difference also persists at small m, though expansions of both ∆W S 4 and ∆W R 4 run in powers of m 2 , in the latter case the accompanying coefficients are enhanced by powers of log m.
The underlying reason why the circular Wilson loop is different in the flat space and on the sphere is due to conformal symmetry breaking in N = 2 * SYM for m = 0. Naively one may expect that ∆W S 4 and ∆W R 4 should coincide in the limit when the radius of the foursphere becomes large, or equivalently mR 1. We just demonstrated that this is not the case. To see why this happens, consider the massive scalar propagator on the sphere (2.22). Restoring the dependence on the radius R, we observe that the propagator is a function of two dimensionless parameters mR and X 2 12 /R 2 , with X 2 12 being a chordal distance on the sphere. Going to the limit R → ∞ with m held fixed, we have to distinguish two cases: X 2 12 /R 2 → 0 and X 2 12 /R 2 = O(R 0 ) corresponding to the short and large distances on the sphere. The scalar propagator coincides with its flat-space counterpart in the former case only. The latter situation is realized for the circular Wilson loop when the scalar propagator is stretched across the great circle of the sphere.

Circular Wilson loop at two loops
In this section, we apply the relation (3.11) and work out the representation for ∆W S 4 to order O(g 4 ) in terms of Feynman integrals on the sphere. The technique for computing these integrals is addressed in the following sections.
Applying (3.11), we have to compute one-loop corrections to propagators of the gluon and massless scalar. The corresponding Feynman diagrams are shown in figure 1. In close analogy with the flat-space analysis of section 3.2, their contribution is given by the product of propagators integrated over the position of the interaction vertices on the sphere. Performing calculations, we will use the embedding coordinates, which prove to be very efficient.
In what follows we use the notation X N (s) (with 0 ≤ s < 2π) for points on the great circle of the sphere and Z i,N for the coordinates of the integration vertices σ i . The integration measure on S D isˆd where where x µ are some local coordinates on S D and g µν (x) = η N M (∂Z M /∂x µ )(∂Z N /∂x ν ) is the corresponding metric tensor. We will also need the following tensor

Corrections to gluon propagator
Using the Feynman rules previously presented in section 2, we deduce the following representation for the contribution of diagrams displayed in the first line of figure 1 where we replaced the gauge propagator with (2.18) and introduced a notation for the one-loop polarization operator contracted with the tangent vectorsẊ 1 andẊ 2 . It can be split into the sum of three terms according to the graphs shown in the top row of figure 1, respectively. We start with the contribution of the fermion loop. It takes the form where the first two factors come from the contraction of the gauge field (2.17) with the tangent vectors and e a,µ (x) are the vielbeins. Replacing the fermion propagator D Ψ with its explicit expression (2.23), we encounter the following expression Here, in the first line, we applied (2.20), while in the second line, we used (C.6) and (4.3). Taking into account the last relation, we get from (4.6) where D Ψ (Z 1 , Z 2 ) is the fermion propagator (2.23) stripped down from the rotation matrices U and U −1 The relation (4.8) has a form that is similar to the fermion loop correction to the gauge propagator in a flat space. One important difference, however, is that all vector indices are contracted using the tensor Q N M . This ensures gauge invariance of Π G 1 [16,19]. Going through calculations of (4.8), we obtain 4 where the functions A(Z 2 ) and B(Z 2 ) are given by (2.24). 4 Arriving at this relation, we neglected terms containing total derivatives with respect to s1 and s2.
They yield vanishing contributions in the circular Wilson loop upon the integration over the contour in eq. (3.11).

JHEP04(2021)089
The calculation of the scalar loop goes along the same lines . (4.11) As in the previous case, this expression is similar to its analogues in the flat space. Here the two terms in the sum describe the contribution of the (pseudo)scalars A i and B i , respectively. Their propagators are defined in (2.21) and (2.22). The functions S ± satisfy nontrivial relations (see eqs. (D.1) and (D.9)) and they can be expressed in terms of the functions A(Z 2 ) and B(Z 2 ) introduced in eq. (2.24). Using these results, we can cast Π G 2 into the form Combining together (4.10) and (4.12), we find that the sum Π G 1 + Π G 2 takes a remarkably simple form where Z 2 12 = (Z 1 + Z 2 ) 2 = 4 − Z 2 12 and Z 12 = Z 1 − Z 2 . Finally, the tadpole contribution to the gauge propagator is localized at Z 1 = Z 2 and it is proportional to propagators of the massive scalars at zero separation where the delta function on the sphere is defined as´dσ 1 δ(Z 12 )f (Z 1 ) = f (Z 2 ) for an arbitrary test function.

Corrections to scalar propagator
Analogously, for the one-loop correction to the massless scalar propagator, we get where the scalar polarization operator is again a sum of three contributions for the Yukawa, triple scalar and tadpole contributions shown in the bottom row of figure 1, respectively. They read The first expression should be compared with the analogous equation in (4.13).

General expression
Combining together eqs. (4.4) and (4.15), we conclude that the one-loop correction to the linear combination of gauge and scalar propagators in (3.11) is given bẏ (4.18) where D loop and D tadpole describe contributions from diagrams in figure 1 containing loops and tadpoles, respectively. We take into account (4.14) and (4.17) to find Here in the first relation, we used δ(Z 12 ) to integrate over Z 2 and, in the second relation, we introduced notations for where Z 2 12 = 4 − Z 2 12 . These equations involve functions defined in (2.22) and (2.24). In the next section, we describe a technique for computing integrals entering (4.19).

Propagators at zero separation
Before we turn to computing all integrals, let us examine the contribution from the tadpoles (4.19). It depends on the scalar propagators S ± (0) at coincident points. To find them, it is convenient to use the relation (2.24) for X 2 12 → 0 (see also (D.1)) , (4.21) so that the sum of the two propagators becomes (4.22) To find A(0), it is instructive to use the Mellin-Barnes representation of A(X 2 12 ) derived in the appendix D. As a result, we find (see eq. (D.10)) . (4.23) For D = 4 − 2 , we expand (4.22) in powers of to get where H(x) = ψ(x + 1) + γ are the harmonic numbers. In the above equation, the pole 1/ has a UV origin. In the two-loop expression for the circular Wilson loop, it cancels against contributions of other diagrams in figure 1.

Comparison with localization
In this section, we compute the integrals (4.19) on the sphere and, then, use them to evaluate the circular Wilson loop (3.11). Echoeing (4.18), we can split ∆W S 4 into the sum of two terms coming from loop and tadpoles (5.1)

Contribution of tadpoles
We start with the contribution of tadpole, which reads Replacing D tadpole with its explicit expression (4.19), we obtain where S 0 is determined by eq. (4.24). Here we took into account that S 0 vanishes for m = 0 and introduced the notation for the integral where the integration contour separates increasing and decreasing poles stemming from Γ(−z + . . . ) and Γ(z + . . . ), respectively. The details of the calculation can be found in the appendix E. We recall that X i = X(s i ) parameterizes points on the great circle of the sphere (see eq. (3.3)). Substituting (5.4) into (5.3), we encounter the following contour integralŝ 2π 0 ds 1 ds 2 X 2 12 /4 Notice that these two only differ by a factor of (−z/(z + 1)) in the right-hand side.
Combining the above relations together, we arrive at .
The integral is finite for = 0 and it can be expanded at small using the MB Tools package [33]. In this way, we find

JHEP04(2021)089
This relation defines the contribution of the sum of the diagrams G 3 and S 3 in figure 1,

Contribution of loops
In a close analogy with (5.2), the contribution of graphs containing massive loops is where D loop is given by (4.19). Due to quite an intricate form of the functions involved, a direct closed-form evaluation of the integral in eq. (4.19) is way too complicated. As we show below, these difficulties can be alleviated by employing a Mellin-Barnes representation of the polarization operators (4.20) In this manner, replacing the scalar propagator in (4.19) with its expression (2.16), we find that the integrals over the sphere in (4.19) can be expressed in terms of a 'simple' integral Here, the amplitude I(z, j) is a meromorphic function of z and j Its derivation can be found in appendix E, (see eq. (E.12)). Combining together the above relations, we get from (5.8) the representation of W loop as a double Mellin-Barnes integral where the amplitudes Π 1 (j) and Π 2 (j) can be obtained from (5.9) and (4.20) by replacing the functions A, B and S ± by their Mellin-Barnes images (see eqs. (D.2) and (D.5)). The contour integrals in the second line of (5.12) can be evaluated using (5.5). This leads to

JHEP04(2021)089
where the notations were introduced for Here the functions M A , M B and M S describe contributions of various terms in the righthand side of eq. (4.20). Their explicit expressions in terms of the Mellin-Barnes amplitudes of the functions A, B and S ± are We recall that the sum of graphs (G 1 ) and (G 2 ) yields Π 1 in the first equation of (4.20) and, therefore, it produces the contribution to (5.13) proportional to M 1 . The diagram (S 1 ) generates the first two terms in the expression for Π 2 in (4.20) and its contribution to (5.13) is described by the sum of two terms M A +M B in the expression for M 2 in (5.14).
The remaining M S term in (5.14) comes from the graph (S 2 ). The calculation of the contribution to (5.13) from M B term in the expression for M 2 is presented in appendix F. The contribution of the remaining terms can be found in a similar manner. The resulting expressions for the individual graphs in (5.16) take the form of (complicated) triple Mellin-Barnes integrals. In the next subsection, we present their expansion in the limit of small and large masses. For finite mass, we used the MB Tools package [33] to compute them numerically.

Two loops versus localization
In this subsection, we summarize the obtained results and present the two-loop expressions for the circular Wilson loop where α = {G 1 + G 2 , S 1 , S 2 , G 3 + S 3 } enumerates the graphs shown in figure 1.
In the small mass limit, the individual contributions of the diagrams admit the following general form

JHEP04(2021)089
Explicit expressions for the leading coefficients w The subleading coefficients are given by w (1) where w (1) G 3 +S 3 was obtained by expanding the tadpole contribution (5.7) at small m. Substituting (5.18) into (5.17), we find in the limit → 0 In agreement with our expectations, UV poles cancel in the sum of all graphs. Moreover, finite rational terms and terms proportional to log 2 also cancel against each other as well.
The resulting expression for ∆W S 4 coincides with the localization prediction (3.6) and (3.8).
In the opposite large-mass limit, the contribution of the graphs in figure 1 looks like with the explicit expressions for the coefficients being w (∞) As in the previous case, we observe perfect agreement with the localization prediction, eqs. (3.6) and (3.9).
For finite values of the mass, we computed (5.17) numerically for various values of m and reproduced the expected result (3.6) and (3.7) to high precision.

Conclusions
In this paper, we developed a framework that allowed us to compute the circular Wilson loop in the N = 2 * super Yang-Mills on the four-sphere at two loop order. We verified that it perfectly agrees with the prediction of the supersymmetric localization but differs from analogous perturbative calculations in the flat space. In the latter case, the difference arises due to the presence of a mass scale in the theory. The reason being that the mass deformation explicitly breaks the conformal symmetry of the model and, as a consequence, a coordinate transformation from the sphere to the flat space becomes anomalous. This question deserves further investigation.
The main findings of our work are rather technical but they have a potential to be of value in other circumstances as well. We demonstrated that the calculation on the sphere can be simplified by employing the embedding coordinates instead of local coordinates (like spherical angles or stereographic coordinates). We derived, in particular, a new, very concise representation for a massive fermion propagator which proves to be very convenient in evaluating Feynman integrals on the sphere.
Having worked in the difference theory, we encountered two-loop graphs of propagator type only, that is, one-loop graphs with massless propagators modified by corrections from massive fields circulating in virtual loops. As a first step in computing these graphs, we proposed to use the Mellin-Barnes representation for the massive propagators. The main advantage of this representation is that it allows us to disentangle the dependence on the mass and the chordal distance on the sphere and, as a consequence, to reduce our calculations merely to the evaluation of integrals on the sphere containing products of powers of chordal distances. With all integrations done, we arrived at a compact representation for the two-loop circular Wilson loop as a sum of three-fold Mellin-Barnes integrals involving ratios of Euler Γ−functions. These integrals are performed in the complex plane along the contours separating increasing and decreasing sequences of poles generated by the latter.

JHEP04(2021)089
Using the Mellin-Barnes representation as a starting point, we can apply conventional techniques, reviewed in ref. [34], and compute emerging integrals by residues. We did not succeed in getting their closed-form expressions for finite value of the mass since a brute force use of Cauchy theorem is unwieldy and results in nested infinite series representation that we failed to resum. This is the reason why we limited ourselves only to their small and large mass expansions. For finite value of the mass, however, we computed contributing Mellin-Barnes integrals numerically to high precision. In all cases, we found that the two-loop result for the circular Wilson loop is in an exact agreement with the localization. Our consideration can thus be regarded as a first test of the latter in massive non-conformal settings.
The eventual two-loop result for the circular Wilson loop arises after massive cancellations between contributing two-loop graphs. This observation, as well as the simplicity of the localization formulae, hints that the Mellin-Barnes integrals mentioned above can be computed exactly for finite value of the mass as well. It would be interesting to observe this explicitly.
The presented techniques has a virtue of potential application even in massless setups and may resolve some existing puzzles. As was discovered in ref. [35], the bremsstrahlung function in N = 2 superconformal quiver theories, which is not immediately attainable by the localization proper, can be related to an elliptical Wilson loop, which is. While the comparison with perturbative calculations in flat space indicated complete agreement to the lowest three orders in the coupling, it was observed [36,37] that there is a deviation starting from four loops. It was argued in ref. [37] that the discrepancy can be accommodated by the effective mass term of the conformally coupled scalar in curved space-time. We hope that our method can be applied in these circumstances and shed light on the issue.
The formalism developed in this work can be used in computation of other observables in massive super Yang-Mills theories on the four-sphere, e.g., correlation functions of local operators or Wilson loops with insertions of these. Also, since de Sitter space can be obtained in the embedding coordinates by a Wick rotation, all ingredients of our perturbative construction can be used in cosmological applications as well.

A Local coordinates on the sphere
The spherical coordinates on the sphere X 2 = R 2 are defined as X 0 = R cos θ 1 , . . .
The volume element then reads where dΩ (D+1) is the element of the solid angle in the above spherical coordinates. The total volume of S D is

B Scalar propagator
In this appendix, we derive the massive scalar propagator on the sphere (2.21) and (2.22). As was mentioned in section 2.2, the calculation can be simplified by Wick rotating to the de Sitter space and using the coordinates (2.14). The de Sitter space is conformally flat and its metric can be cast into the form with η µν = diag (+, −, . . . , −) and z 0 being the conformal time (2.15). The Laplace-Beltrami operator admits the form The propagator of a conformally coupled real scalar of mass µ obeys the equation where the subscript 1 refers to the point z 1 . The function G µ (z 1 , z 2 ) depends on the coordinates only through the chordal distance, where a notation was introduced for M = 1 − 4µ 2 . In the vanishing mass limit, i.e., for µ → 0, we obtain After the analytic continuation to the sphere, for X 0 i → −iX 0 i , this relation yields the massless scalar propagator on the sphere (2.16). The the N = 2 * theory contains two pairs of scalars with masses µ ± , see eq. (2.10). Substitution of µ = µ ± into (B.6) yields the two propagators (with M = 1 ± 2im) Upon the analytic continuation to the sphere, these provide the propagators for A and B scalars introduced in (2.21) with S ± (X 2 12 ) = G µ ± (s) for s = X 2 12 = (X 0 1 − X 0 2 ) 2 + i (X i 1 − X i 2 ) 2 .

C Fermion propagator
In this appendix, we derive the massive fermion propagator (2.23). 5 We do it a step-wise manner, first, by finding the propagator in the de Sitter coordinates (2.15), and then, lifting it to the embedding space. It is the latter form that we use in the main body. The propagator of the fermion with the mass m obeys the following equation

JHEP04(2021)089
We start with the functions S ± (X 2 ) given by eq. (2.22). Using the properties of the hypergeometric functions (see eq. (E.8)), we can immediately obtain its Mellin-Barnes representation where its image is given by 3) The integration contour in (D.2) separates increasing and decreasing poles generated by Γ-functions of the form Γ(−z + . . . ) and Γ(z + . . . ), respectively. The main advantage of this form is that the integrand has a power-like dependence on the chordal distance X 2 and is better suited for performing integration on the sphere. Substituting (2.22) into (2.24), we can re-express the functions A(X 2 ) and B(X 2 ) in terms of the hypergeometric functions

JHEP04(2021)089
Large mass limit. Finally, in the large-mass limit, we replace sinh(πm) = exp(πm)/2 + . . . and substitute the product of the Γ−functions in the second line of (F.4) with its asymptotic behavior as m → ∞, where x = 2 − + z. The resulting integral can again be evaluated using the MB Tools package [33] leading to (F.12) The integral develops a simple UV pole and has the following unusual feature -it receives O(m) correction. It is ultimately related to poles generated by Γ(z 3 + 1 2 ) which arise from the integration over the circle, see eq. (5.5).
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.