A matrix model for the latitude Wilson loop in ABJM theory

In ABJ(M) theory, we propose a matrix model for the exact evaluation of BPS Wilson loops on a latitude circular contour, so providing a new weak-strong interpolation tool. Intriguingly, the matrix model turns out to be a particular case of that computing torus knot invariants in $U(N_1|N_2)$ Chern-Simons theory. At weak coupling we check our proposal against a three-loop computation, performed for generic framing, winding number and representation. The matrix model is amenable of a Fermi gas formulation, which we use to systematically compute the strong coupling and genus expansions. For the fermionic Wilson loop the leading planar behavior agrees with a previous string theory prediction. For the bosonic operator our result provides a clue for finding the corresponding string dual configuration. Our matrix model is consistent with recent proposals for computing Bremsstrahlung functions exactly in terms of latitude Wilson loops. As a by-product, we extend the conjecture for the exact $B^{\theta}_{1/6}$ Bremsstrahlung function to generic representations and test it with a four-loop perturbative computation. Finally, we propose an exact prediction for $B_{1/2}$ at unequal gauge group ranks.

Abstract: In ABJ(M) theory, we propose a matrix model for the exact evaluation of BPS Wilson loops on a latitude circular contour, so providing a new weak-strong interpolation tool. Intriguingly, the matrix model turns out to be a particular case of that computing torus knot invariants in U (N 1 |N 2 ) Chern-Simons theory. At weak coupling we check our proposal against a three-loop computation, performed for generic framing, winding number and representation. The matrix model is amenable of a Fermi gas formulation, which we use to systematically compute the strong coupling and genus expansions. For the fermionic Wilson loop the leading planar behavior agrees with a previous string theory prediction. For the bosonic operator our result provides a clue for finding the corresponding string dual configuration. Our matrix model is consistent with recent proposals for computing Bremsstrahlung functions exactly in terms of latitude Wilson loops. As a by-product, we extend the conjecture for the exact B θ 1/6 Bremsstrahlung function to generic representations and test it with a fourloop perturbative computation. Finally, we propose an exact prediction for B 1/2 at unequal gauge group ranks. Wilson loop operators are fundamental observables in any gauge theory. While in the case of pure Yang-Mills in two dimensions a complete solution for their vacuum expectations value exists [1,2], in higher dimensions only few examples of non-perturbative calculations can be found in literature. A notable exception is provided by ordinary Wilson loops in pure three-dimensional Chern-Simons theory [3], where the computation is made possible thanks to the topological nature of the model. A more significant class of examples is represented by the so-called supersymmetric/BPS Wilson loops, appearing quite ubiquitously in gauge theories with extended supersymmetry. Their main property consists in preserving a fraction of the supersymmetry charges, depending on the shape of the contour and on the couplings to the different fields appearing in the Lagrangian. Despite their BPS nature, the dependence of their vacuum expectation values on the coupling constant is generically non-trivial and interpolates between the weak and strong coupling regimes, thus providing a natural playground where to test the AdS/CFT correspondence and other non-perturbative methods.
In recent years the technique of supersymmetric localization (for recent reviews, see [4] and references therein) has allowed for the calculation of a large variety of BPS Wilson loops, in different theories and various dimensions [5][6][7][8][9][10][11][12]. Localization often reduces the computation of these observables to the average of suitable matrix operators, in terms of particular matrix integrals. These integrals can then be solved by applying the powerful machinery developed along the years, as for example large N expansion, orthogonal polynomials, loop equations, recursion relations.
In this paper we discuss a new example of such matrix model computations of Wilson loops. We focus on three-dimensional N = 6 superconformal U (N 1 ) × U (N 2 ) Chern-Simons theory with matter, also known as the ABJ(M) model [13,14], that can be viewed as an extended supersymmetric generalization of the familiar topological Chern-Simons theory. We propose a matrix model for calculating the exact quantum expectation value of certain 1/12-and 1/6-BPS Wilson loop operators, briefly referred to as bosonic and fermionic latitudes and parameterized by a real number ν ∈ [0, 1] [15]. For ν = 1 we recover the supersymmetric Wilson loops on the great circle [9,[16][17][18][19], for which a precise localization procedure has been derived in [8], leading to the ABJ(M) matrix model. Various aspects of such circular Wilson loops have been thoroughly studied in [20,21].
For generic ν we do not possess a direct derivation of our matrix model from supersymmetric localization, since these loop operators preserve supercharges that are different from the ones used in [8], and cannot be embedded in a natural way in the N = 2 superspace language employed there for the localization procedure. Rather, we formulate an ansatz for the matrix integral from symmetry considerations and show its non-trivial consistency with a three-loop perturbative calculation of the same observables. At equal gauge group ranks, we also prove that, without the operator insertion, the matrix model reproduces the ABJM partition function. This fact supports a possible interpretation of our matrix model as the result of localizing with the ν-dependent supercharge preserved by latitude operators. Encouraged by this fact, we perform a strong coupling analysis using the Fermi gas approach. In the fermionic case we obtain an important agreement with the existing string theory computation [22]. Instead, in the case of unequal ranks, N 1 = N 2 , we find a residual ν-dependence in the matrix model partition function. This would in principle prevent us form directly interpreting it as the result of localizing the ABJ theory on the sphere. Still, such a dependence is confined to a simple phase factor, which might stem from a framing anomaly. Therefore, the latter case requires a deeper analysis.
Latitude Wilson loops are important in order to study non-BPS observables, the Bremsstrahlung functions [15,23]. These are functions of the coupling constant, controlling the small angle limit of the anomalous dimension of a generalized cusp, constructed from supersymmetric Wilson lines. In particular, depending on the degree of supersymmetry, the various Bremsstrahlung functions can be obtained taking a suitable derivative of the latitude with respect to the ν parameter, at ν = 1 [15,22,23]. The fermionic and bosonic latitudes satisfy a cohomological equivalence [15]. Using the latter and under a suitable assumption (automatically satisfied by our matrix integral), remarkable relations can be derived for the different Bremsstrahlung functions at equal ranks N 1 = N 2 . Eventually, they are all related among themselves and expressible only in terms of the phase of the undeformed 1/6-BPS Wilson loop.
Besides the derivation and the strong coupling solution of our matrix model, the paper provides details of the perturbative three-loop computation of the bosonic latitude Wilson loop, for generic representation and winding number. This general result allows for comparisons in various limits. We also present a careful discussion of the framing dependence [24] at perturbative level. In particular, the perturbative result at generic framing is crucial for finding consistency with the prediction from our matrix model. The remarkable output is that the agreement between the perturbative and the matrix model calculations works upon the identification of framing with the effective parameter ν.

Summary of the results
The main focus of the paper is on the proposal of a matrix model computing the latitude expectation values at all orders. Based on symmetry considerations and the consistency with the perturbative results, we conjecture that the exact expectation value of the multiply wound bosonic latitude Wilson loop in the fundamental representation is given by where · · · stands for the normalized expectation value performed using the following matrix integral The expression for the fermionic 1/6-BPS latitude Wilson loop is obtained by computing a suitable linear combination of bosonic latitudes, as suggested by the cohomological equivalence discussed in [15] and reviewed in Section 2.
We also conjecture that the expectation values of latitude Wilson loops in higher dimensional representations can be obtained from this matrix model by the same generalization of the operator insertions as in the undeformed (ν = 1) case.
A crucial point for the consistency of the above proposal is that for N 1 = N 2 the partition function Z(ν) defined in (1.2) is independent of the ν parameter and coincides with the usual partition function of ABJM on S 3 . This important property is proved explicitly in Section 4.2. The ABJ case (i.e. N 1 = N 2 ) is subtler but still the full dependence on ν is confined to a simple phase factor (see Section 4.2). This points towards the possibility that a localization procedure performed with one of the ν-dependent supersymmetric charges preserved by the latitude could produce the above matrix integral.
Supported by this first evidence on the correctness of our proposal, in Section 5 we perform a careful analysis of the large N , strong coupling limit, of the expectation values of the bosonic and fermionic latitudes, as defined through our conjectured matrix model. We perform a Fermi gas analysis along the lines of [25] and obtain an explicit result re-summing the genus expansions both in the 1/12-BPS and 1/6-BPS cases. In particular, for the fermionic loop we obtain From this expression we can read the leading genus-zero term that consistently coincides with the semi-classical string computation of the 1/6-BPS loop performed in [22]. Recovering the string result is non-trivial and depends on a delicate cancellation of various contributions appearing in the intermediate steps of the calculation.
The matrix model proposal must also be consistent with explicit weak-coupling calculations. Therefore, we have performed a perturbative three-loop computation of the bosonic latitude, at generic framing f , representation R and winding number m. This Feynman diagram tour de force agrees with the third-order expansion of (1.1), provided the framing number is formally identified with ν. The agreement holds for generic winding and, in addition, we have explicitly verified the consistency of our matrix model proposal with the perturbative result for a few higher dimensional representations of the gauge group. The requirement of a specific choice of framing does not come as a surprise, in fact localization for Wilson loops on the great circle in ABJ(M) theory produces results at non-trivial framing (see [8]). Hence, it is reasonable to expect that our matrix model computation also yields results at a fixed framing number. The particular choice f = ν is natural, as it corresponds to the value at which the bosonic and fermionic latitude expectation values are related by the cohomological equivalence at the quantum level, as pointed out in [15]. This requires an analytic continuation of the framing parameter from an integer to a real number, which is perfectly legitimate at the matrix model level.
We discuss the different Bremsstrahlung functions that can be defined in ABJ(M) theories, and in particular we review how they are connected to suitable derivatives of the modulus of the latitude expectation value. The introduction of the modulus slightly modifies the prescription considered previously [15,22,23], allowing to eliminate some unexpected imaginary contributions to the Bremsstrahlung functions, appearing at three loops from the perturbative expansion of the latitudes at framing zero. The origin of these imaginary terms is quite peculiar and is rooted into an anomalous behavior of some correlation function of scalar composite operators, as we discuss carefully in subsection 3.4. Once we take into account this effect, it is straightforward to understand why the correct prescription for the Bremsstrahlung functions must contain the modulus, analogously to what was originally argued in [26]. We test the matrix model based expression for the Bremsstrahlung function associated to the internal angle of the bosonic cusp, B θ 1/6 , against its four-loop perturbative computation for a generic representation, which is presented in (2.30).
In the case of the ABJM theories, the following relation, was originally conjectured in [15] and is at the core of the relation between Wilson loops and Bremsstrahlung functions proposed there. This identity holds for the bosonic latitudes We explicitly check that (1.4) is obeyed by our three-loop perturbative result and, remarkably, it is an exact consequence of the proposed matrix model. In fact we prove that more generally our matrix model satisfies this identity not only at ν = 1 but for any value of ν. We also reproduce the three-loop result of [27] for the fermionic cusp anomalous dimension in the near-BPS regime. This chain of relations among the different Bremsstrahlung functions is also discussed in [28].
Finally, we briefly discuss a possible generalization of our approach for computing the fermionic Bremsstrahlung function from the 1/6-BPS fermionic latitude in the N 1 = N 2 case. We put forward an exact prediction for B 1/2 , whose expansion up to five loops is provided explicitly in (6.2). This result calls for future perturbative confirmation.
The paper is structured in the following way. We start in Section 2 by reviewing the construction of two general families of circular Wilson loops, whose contours are the latitudes on a sphere S 2 . The two families, bosonic and fermionic latitudes, differ by the nature of the connection appearing in the holonomy and the number of preserved supercharges. We discuss the framing dependence of Wilson loops expectation values, the cohomological relation between bosonic and fermionic latitudes and the connections among the different Bremsstrahlung functions, associated to fermionic and bosonic cusps. Subsequently, Section 3 illustrates the three-loop perturbative calculation of the bosonic latitude in the general situation described before. In particular, we elucidate the emergence of an imaginary term at framing zero at three loops and the requirement of considering the modulus of the latitude expectation value in computing the Bremsstrahlung function. Section 4 and Section 5 are the heart of the paper: here we propose our matrix model and study its main properties by using the Fermi gas technique. The strong coupling expansion is performed and successfully compared, in the 1/6-BPS fermionic case, with the semi-classical string computation. In Section 6 we briefly present and discuss a conjectured form of the fermionic Bremsstrahlung function for N 1 = N 2 , providing a prediction up to five loops. Six appendices complete the paper with conventions, details of the computation and further checks of our results.

BPS Wilson loops, Cusps and Bremsstrahlung functions: A review
We begin with a general review of the most fundamental properties of circular BPS Wilson loops in ABJ(M) theories and their non-trivial connections with cusped Wilson lines and the corresponding Bremsstrahlung functions. In U (N 1 ) k × U (N 2 ) −k ABJ(M) theory we consider the general class of Wilson loops that preserve a certain fraction of the original N = 6 supersymmetry 1 . Such operators can be constructed by generalizing the ordinary gauge holonomy with the addition of either scalar matter bilinears ("bosonic" Wilson loops) [16][17][18][19] or scalar bilinears and fermions ("fermionic" Wilson loops) [9]. For the case of straight line and maximal circular contours a general classification of such BPS operators based on the amount of preserved supercharges can be found in [10,11,29].
Latitude Wilson loops. We are primarily interested in the general class of bosonic and fermionic Wilson operators introduced in [30] (latitude Wilson loops). They feature a parametric dependence on a α-angle 2 that governs the couplings to matter in the internal R-symmetry space and a geometric angle θ 0 ∈ [− π 2 , π 2 ] that fixes the contour to be a latitude circle on the unit sphere Γ m : x µ = (sin θ 0 , cos θ 0 cos τ, cos θ 0 sin τ ) τ ∈ [0, 2mπ) for winding m (2.1) Note that here we are generalizing the definitions of [30] to Wilson loops with generic winding. As discussed in [15], these operators can be constructed in such a way that they depend uniquely on the effective "latitude parameter" The m-winding bosonic latitude Wilson loops corresponding to the two gauge groups are explicitly given by where the matrix describing the coupling to the (C I ,C I ) scalars reads The traces in (2.3) are taken over generic representations R,R of U (N 1 ) and U (N 2 ), respectively. The overall constants have been purposely chosen in order to normalize the tree level expectation values W m B (0) and Ŵ m B (0) to one. Similarly, the m-winding fermionic latitude Wilson loop for a generic representation R of the superalgebra U (N 1 |N 2 ) is defined as and The generalized prescription (2.5) that requires taking the supertrace of the superholonomy times a constant matrix assures invariance under super gauge transformations [30]. The overall constant in (2.5) can be chosen so as to normalize the expectation value to 1, if possible. In the rest of the paper we consider the fermionic operator only in the fundamental representation, for which We note that for N 1 = N 2 , if ν = 1 (ν = 0) and m is even (odd) this normalization becomes meaningless. In those cases one can simply compute the unnormalized expectation value, choosing R = 1. Whenever no confusion arises we will use W B,F as a shorthand for the single winding operators W 1 B,F . Moreover, if no explicit dependence on the representation is displayed, the Wilson loop is understood to be in the fundamental representation.
At classical level the fermionic latitude Wilson loop (2.5) is cohomologically equivalent to the following linear combination of bosonic latitudes where for simplicity we have restricted to Wilson loops in the fundamental representation. In the above formula Q(ν) is the linear combination of superpoincaré and superconformal charges [15] preserved by both bosonic and fermionic Wilson loops. If this equivalence survives at quantum level it allows to compute the vacuum expectation value W F (ν) of the fermionic operator as a combination of the bosonic ones. However, in three dimensions the problem of understanding how the classical cohomological equivalence gets implemented at quantum level is strictly interconnected with the problem of understanding framing, as we review below.
Matrix models for BPS Wilson loops. Using the procedure of supersymmetric localization the ABJ(M) partition function on the three-sphere can be reduced to the matrix model integral [8] where we are being cavalier on the precise normalization, which is unimportant for the computation of Wilson loops. At ν = 1 the expectation values of the 1/6-and 1/2-BPS Wilson loops can be computed as matrix model averages. In particular, the m-winding bosonic 1/6-BPS Wilson loop in the fundamental representation 3 is given by where the right-hand-side brackets stand for the integration using the matrix model measure defined in (2.12), normalized by the partition function. In this language, the 1/2-BPS Wilson loop can be computed as the average of a supermatrix operator, or equivalently using (2.10) [9]. For generic ν, no matrix model prescription has been found so far. In fact, even though the latitude Wilson loops are BPS operators, in this case the standard localization arguments of [8] cannot be directly applied [15]. We aim at filling this gap in Section 4, where we conjecture a matrix model for the latitude Wilson loops that turns out to be compatible with all the available data points at weak and strong coupling.
Framing. In three-dimensional Chern-Simons theories the computation of Wilson loop expectation values is affected by finite regularization ambiguities associated to singularities arising when two fields running on the same closed contour clash. In perturbation theory, this phenomenon is ascribable to the use of point-splitting regularization to define propagators at coincident points. For ordinary Wilson loops, following e.g. the prescription of [31], one allows one endpoint of the gluon propagator to run on the original closed path Γ on which the Wilson loop is evaluated, and the other to run on a framing contour Γ f . This is infinitesimally displaced from Γ and defined by the choice of a vector field on it. Then the one-loop Chern-Simons contribution is proportional to the Gauss linking integral which evaluates to an integer f (the framing number). This is a topological invariant that counts the number of times the additional closed contour Γ f introduced by the framing procedure winds around the original one Γ. This phenomenon has been first discovered and extensively discussed for pure U (N ) Chern-Simons theory [3], in connection to knot theory. In this case, the total effect of framing amounts to a phase, exponentiating the one-loop result (from now on we use the subscript on the expectation value to indicate a certain choice of framing) where λ is the 't Hooft Chern-Simons coupling shifted by the quadratic Casimir of the gauge group.
More recently, the same kind of framing dependence has been discussed also for non-topological Chern-Simons theories coupled to matter, in particular ABJ(M) theories [24,32]. In order to review framing effects in this context it is convenient to split the Wilson loop expectation value into its phase and its modulus. For the most general case of m-winding operators with a non-trivial latitude we set and similarly for W m F (ν). In the ν = 1 case, the bosonic 1/6-BPS operators have been computed up to three loops in the large N 1 , N 2 limit, for generic framing number and winding m [24,33]. As an effect of framing, their expectation values acquire imaginary contributions at odd orders, as well as additional real corrections at even orders. We stress that in this case the imaginary contributions are entirely due to framing.
For single winding the framing contributions can be still captured by a phase, precisely Φ B (f, 1, 1),Φ B (f, 1, 1) in (2.16), while the expectation values at framing zero are real quantities and coincide with the modulus. However, the phases are no longer a one-loop effect as in the pure Chern-Simons theory, but display non-trivial quantum corrections starting at three loops [24] Φ 17) For multiple windings the effect of framing is more complicated and ceases to be encapsulated into a phase [33]. This does not come as a surprise since the same pattern occurs also in the pure Chern-Simons theory.
For latitude Wilson loops the phases Φ B ,Φ B will depend in general on the framing and winding numbers, and the latitude parameter ν as well. According to the discussion above, they will be non-trivial functions of the couplings that reduce to the expansions (2.17) for ν = 1 and single winding. Having in mind the most general scenario, we may expect them not to necessarily account for all the framing effects (as in the multiple winding situation). Hence, in general the modulus in (2.16) does not coincide with the expectation value at framing zero. Moreover, for latitude operators the phase might not even be entirely produced by framing, as further framing independent imaginary contributions could arise. Checking if this is the case and better understanding the framing origin of Φ B ,Φ B is one of the goals of this paper.
We conclude this short review on framing by discussing its role in localization. It was argued in [8] that the matrix model (2.12) derived from localization computes the 1/6-and 1/2-BPS Wilson loops at framing one. This is because a point-splitting regularization, implied in the derivation, is compatible with supersymmetry only if the circular path and the framing contour are two Hopf fibers in the S 1 fibration of the three-sphere. These in turn have linking number one, which explains the particular framing number arising in this computation. The aforementioned studies on supersymmetric Wilson loops have provided a perturbative test of such an argument.
For the purposes of this paper we therefore stress that a matrix model computing the expectation value of latitude Wilson loops is expected to imply a particular choice of framing.
Cohomological equivalence and non-integer framing. As stated above, at classical level the fermionic and bosonic latitude Wilson loop operators are related by the cohomological equivalence (2.10). If relation (2.10) survives at the quantum level, then W F (ν) can be obtained as a linear combination of W B (ν) and Ŵ B (ν) . In particular, we expect that to be the case (namely no anomalies arise) in the localization approach, if the functional integral computing the Wilson loop expectation value is localized using its invariance under the same supercharge Q in (2.10).
For ν = 1 this reduces to the cohomological equivalence first discovered in [9]. In this case the localization computation is performed at framing 1, as recalled above, hence the equivalence is expected to hold for this particular choice of framing (in this section we restrict to the fundamental representation and set m = 1 for simplicity) but could be modified if another choice of framing is taken. In fact, this has been explicitly verified in ordinary perturbation theory at framing zero up to two loops, where the equivalent of (2.18) with f = 0 fails. On the other hand, from equation (2.17) it follows that up to two loops framing zero and one expectation values of bosonic 1/6-BPS Wilson loops are related as (2.19) Using this, and further defining identity (2.18) has been confirmed to hold, perturbatively [34][35][36].
In the latitude case, the analogous two-loop calculation [15] shows that the cohomological equivalence survives at quantum level in the form that is if we formally identify the framing number f with the latitude parameter ν. Therefore, in the general case, we allow the latitude to be evaluated at non-integer framing ν. Moreover, we expect that a matrix model computation of the latitude Wilson loop, respecting the cohomological equivalence at the quantum level, would imply framing ν.
Bremsstrahlung functions. The Bremsstrahlung function B is the physical quantity that measures the energy lost by a heavy quark slowly moving (|v| 1) in a gauge background. Generalizing the well-known law of electrodynamics, it is defined as [37] ∆E = 2πB dt(v) 2 (2.23) In a conformal field theory it also appears as the coefficient of the first non-trivial order in the small angle expansion of the cusp anomalous dimension, Γ cusp (ϕ) ∼ −Bϕ 2 , which governs the short distance divergences of a cusped Wilson loop.
In ABJ(M) theories, since we have bosonic and fermionic Wilson loops we can define different types of Bremsstrahlung functions [26,38]. Computing the divergent part of a fermionic, locally 1/2-BPS Wilson loop W ∠ F along a generalized cusped contour (two straight lines meeting at a point with a relative ϕ angle) one finds where Λ and are IR and UV cutoffs, respectively. Here θ is the internal angle that describes possible relative rotations of the matter couplings between the Wilson loops defined on the two semi-infinite lines. B 1/2 appears as a common factor in the small angles expansion as a consequence of the fact that for θ = ϕ the fermionic cusped Wilson loop is BPS and divergences no longer appear. Analogously, using a bosonic 1/6-BPS Wilson loop on a cusp we can define and similar relations for Ŵ ∠ B (ϕ, θ) that give rise toB θ 1/6 ,B ϕ 1/6 associated to the second gauge group. In the bosonic case we have in principle two different Bremsstrahlung functions since there are no BPS conditions for cusped bosonic Wilson loops.
A crucial problem consists in relating B to other physical quantities that in principle can be computed exactly using localization techniques, like for instance circular BPS Wilson loops. For the ABJM theory 4 , this problem was originally addressed in [26], where an exact prescription was given to compute B ϕ 1/6 in terms of a m-winding Wilson loop A similar prescription has been later derived for B 1/2 and B θ 1/6 in ABJM [15,22], in terms of single winding, latitude fermionic (2.5) and bosonic (2.3) Wilson loops, respectively These formulae were proven in [23] and [37], respectively. We note that in order to enforce the reality of the result we take the modulus of the expectation values 5 . We derive this prescription in Section 3.4. According to the previous discussion this is also supposed to remove framing ambiguities, henceforth the expectation values in (2.27) can be computed at any convenient framing. These prescriptions have already passed several tests at weak and/or strong coupling. At weak coupling, the lowest order term of B ϕ 1/6 computed using (2.26) agrees with the result obtained from a genuine two-loop calculation of Γ 1/6 cusp [38]. B θ 1/6 obtained from prescription (2.27) has been tested at weak coupling up to two loops [15] for Wilson loops in the fundamental representation of the gauge group. B 1/2 as computed via (2.27) has been tested at weak coupling up to two loops [15,38]. Moreover, the leading term at strong coupling is successfully reproduced by the string dual configuration of W F (ν) found in [22].
A direct perturbative calculation of B θ 1/6 at four loops has been performed in [39,40], and compared with the perturbative result for B ϕ 1/6 as obtained from prescription (2.26). Interestingly, it has been found that the simple relation is valid up to this order and has been conjectured to be true exactly. For the ABJM case (N 1 = N 2 ) this has been proved [28]. For the more general case, taking into account prescriptions (2.26) and (2.27) it amounts to conjecturing that We stress that this applies to generic N 1 and N 2 and no planar limit is assumed. In order to provide further support, in Appendix D we generalize the four-loop test to the case of Wilson loops in generic representations. The result for the Bremsstrahlung function reads and is in agreement with (2.27) upon using the circular Wilson loops in the appropriate representation. We observe that the four-loop contribution exhibits an explicit dependence on higher order Casimir invariants, thereby violating quadratic Casimir scaling, as recently observed in related four-dimensional contexts [41,42]. Concerning B 1/2 , a further point is worth mentioning separately. As discussed in [15], in the N 1 = N 2 situation one can derive from the first equation in (2.27) an exact expression in terms of 1/6-BPS winding Wilson loops This requires making use of the cohomological equivalence in (2.21), assuming a certain relation between W B (ν) and the undeformed m-winding Wilson loop and finally assuming the validity of the following identity As long as we do not know the exact expression for W B (ν) ν , Ŵ B (ν) ν this identity cannot be rigorously proved. However, it has been indirectly verified up to three loops by testing the prediction for B 1/2 from (2.31) against an explicit computation of Γ 1/2 cusp at this order [27]. In analogy with the m-winding case, it is likely to hold at any perturbative order and, in particular, at strong coupling. In fact, as an indirect check, in this regime (2.31) agrees with the explicit string theory computation of the Bremsstrahlung function performed up to the first subleading term [43,44]. Remarkably, as discussed in Section 4 our conjectured matrix model that computes ) not only at ν = 1 but as a functional identity. Assuming (2.32) to be true has far-reaching consequences, the main one being that B θ 1/6 and B 1/2 can be entirely expressed in terms of the phase Φ B introduced in (2.16) for the latitude Wilson loop at framing ν and single winding (we use the shorthand In particular, it follows that a genuine perturbative computation of B θ 1/6 directly from Γ 1/6 cusp allows us to make a prediction for the Φ(ν) function. In fact, exploiting the four-loop result for B θ 1/6 given in [39,40] the following prediction can be made for N 1 = N 2 = N and in the planar limit [28] Note that for f = ν = 1 it consistently reproduces (2.17). Now, merging (2.34) with the two-loop result for W B (ν) 0 [15], from (2.16) we obtain a three-loop prediction for W B (ν) ν in the case of ABJM theory In the next section we are going to test this prediction against a perturbative three-loop calculation done at framing ν. This turns out to be an indirect check of the validity of assumption (2.32) and, therefore, of identities (2.33) relating the Bremsstrahlung functions to the phases of bosonic latitude Wilson loops.

Perturbative result for the latitude Wilson loop
In this section we compute the expectation value of the bosonic latitude Wilson loop W B (ν) at weak coupling in perturbation theory. The evaluation ofŴ B (ν) easily follows by exchanging N 1 ↔ N 2 and sending k → −k.
We consider the most general case where a non-trivial framing number f is taken into account and the contour winds m times around the latitude circle. We also allow for the trace in definition (2.3) to be taken in a generic representation R of the U (N 1 ) gauge group. The U (N 1 ) color factors are expressed in terms of the Casimir invariants, as defined in Appendix A.1. We work at finite N 1 and N 2 , i. e. no planar limit is assumed.
The multiple windings and higher dimensional representations are not independent generalizations, as one can re-express the multiply wound Wilson loop as a linear combination of an alternative basis of operators in different representations [25,45]. Still, in perturbation theory we can treat these two properties independently and use the aforementioned relation as a consistency check of the computation.
In dealing with diagrams contributing to framing, we make use of general properties of the pure Chern-Simons perturbation theory. In particular, we apply the Alvarez-Labastida argument [46], stating that only diagrams with collapsible propagators can contribute to framing, to rule out their non-planar realizations. We argue that up to three loops the whole framing dependence of the Wilson loop can be effectively ascribed to and computed from the Gauss linking integral (2.14). We remark that although the linking number f for two closed curves is naturally an integer number, we will consider its continuation to real numbers.
In Section 3.1 we give details of the calculation for the single winding Wilson loop, whereas the generalization to multiple windings is discussed in Section 3.2. For readers who want to skip technical details we summarize our results in Section 3.3. Finally, in Section 3.4 we revisit the proof of identities (2.26) and (2.27) in light of the appearance of novel three-loop imaginary contributions not related to framing.

The computation
Bosonic and fermionic latitude Wilson loops in the fundamental representation and for single winding have been computed in [15], up to two loops and at framing zero. At non-trivial framing perturbative calculations up to three loops have been carried out in [24] only for W B (1), in the fundamental representation and in the planar limit.
Generalizing those results to W B (ν) in a generic representation, at one loop the only diagram contributing is a gluon exchange, which at non-trivial framing is proportional to the Gauss integral (2.14) . In order to draw higher loop diagrams in a more concise way we find convenient to define a "double-line exchange" given by the combination of a bi-scalar exchange and a one-loop corrected gluon exchange evaluated in dimensional regularization 2) The dependence on the latitude parameter comes from Tr(M 1 M 2 ), where M i stands for the coupling matrix evaluated at point τ i on the contour (see identity (B.10)). It follows that at two loops the following diagrams contribute In (3.5) we sum over all possible planar and non-planar permutations in order to factorize the two-loop diagram as half the squared one-loop graph (3.1). This does not contradict the Alvarez-Labastida argument, since the non-planar crossed configuration is identically vanishing. The non-trivial Feynman diagrams contributing at three loops are depicted in Figure 1, where for diagrams with multiple insertions a sum over all planar configurations arising from permutations of contour points has to be understood.
The triangle graph 1i is a new feature of the latitude Wilson loop stemming from the fact that even though Tr M vanishes, Tr M 3 does not, thus allowing for a non-trivial contribution for ν = 1 (see identity (B.11)).
Except for this graph, all the other diagrams in Figure 1 are structurally the same contributing to the expectation value of W B (1), which were already analyzed in some detail in [24]. We recall that all these diagrams vanish identically at framing zero [19], because of the antisymmetry of the ε tensors appearing ubiquitously in Chern-Simons perturbation theory, but give a non-vanishing contribution at framing f = 0. Other diagrams vanish identically even at non-trivial framing thanks to the tracelessness of the scalar coupling matrix, a property which is true in the undeformed case and remains true in the latitude case, as well. These diagrams have not been included in Figure 1. Diagrams with one-loop corrections to the bi-scalar correlator have also been neglected, as they have been argued to vanish identically [24], independently of framing.
In order to evaluate the diagrams in Figure 1 we can exploit several partial results from [24] to which we refer for more details on the computation. Those results are here generalized to include the latitude deformation, generic representations of the U (N 1 ) gauge group, non-planar contributions and generic framing. In particular, the latitude deformation does affect only diagrams which contain bi-scalar insertions, whereas all the others evaluate exactly as in the undeformed case.
Working with N 1 , N 2 generically different, we can group the diagrams on the basis of their color structures. In particular we find convenient to classify them according to their leading power in N 2 .
Diagrams with no N 2 powers. We start considering the subset of diagrams with no contributions from the U (N 2 ) sector, namely pure U (N 1 ) Chern-Simons contributions. These correspond to the class of diagrams 1a-1c, with all possible planar permutations in the sequence of insertion points, according to the Alvarez-Labastida argument [46]. Having no coupling to the bi-scalar fields, these graphs do not depend on the latitude parameter ν. Therefore they can be evaluated by observing that the combination of all permutations provide a factorization of the diagrams into elementary pieces involving the Gauss integral (2.14), which triggers the framing dependence. The result reads Diagrams with N 2 2 powers. A rather simple class of diagrams is the one with leading N 2 2 behavior, emerging from graphs 1d featuring gauge boson corrections at two loops [24,51]. Collecting the results of [24] for the gauge two-point function insertion diagram, and extending them to the most general color structure (no new topologies arise for this case and the relevant color factors can be found in Appendix B.1) we find Diagrams linear in N 2 . The most complicated contribution to the three-loop expectation value comes from diagrams with leading linear N 2 behavior. These include the factorized diagrams 1e, the interaction diagrams 1f-1h and the triangle graph 1i. The most efficient way to handle factorized diagrams is to sum over all possible planar and non-planar configurations, recalling that only their contractible configurations contribute to framing, whereas the remaining ones vanish identically. Referring to the diagram in Figure 1e, this leads to the following factorization where iπ 3 f is the value of the corresponding integral, whereas the rest comes from color and combinatorics. A latitude dependent part arises from the double-line exchange (3.2). The result can be correctly interpreted as emerging from the interference of the one-loop framing phase and the two-loop perturbative result from diagram (3.3), reproducing the expected exponentiation of framing. Interaction diagrams 1f-1h are the most complicated and we do not possess an exact expression for each individual graph. However, we can indirectly argue the value of their sum as follows. In [24] it was explained that for the 1/6-BPS Wilson loop W B (1) in the fundamental representation, consistency with the localization result requires the sum of 1f-1h to cancel a suitable piece of the gauge two-point function contribution (3.8). As a first consistency check we have verified that the same reasoning holds also for a generic representation of the gauge group, as each interaction diagram has precisely the same color factor as the gauge two-point function contribution. The total sum reads (3.10) Turning on the latitude deformation seemingly spoils such an argument, since diagram 1f acquires a ν-dependent factor from the trace of two M matrices, equation (B.10), that would sabotage the balance required for the aforementioned cancellation. However, the (ν 2 −1) term there, which would be absent in the undeformed case, is proportional to sin 2 τ 1 −τ 2 2 and therefore vanishes for colliding insertion points. Consequently, it protects the integrand from developing the singularity that might cause a potential dependence on framing. In fact, the corresponding integral can be shown to be framing independent and when evaluated at framing zero it vanishes. The remaining ν-independent term in (B.10) obviously yields the same contribution as in the undeformed case. Altogether, the latitude deformation plays no role in the analysis of the interaction diagrams and we are led to postulate that the whole contribution in the latitude case is precisely the same as in the ν = 1 case, that is diagrams 1f-1h cancel the same part of the two-loop gauge propagator (3.8).
Under this assumption, in the latitude case the only extra contribution is the triangle diagram with three bi-scalar insertions, Figure 1i. By direct inspection of the integrand, it is manifest that no framing dependence arises, as there are no singularities for coincident points. Therefore we can evaluate the diagram at framing zero and obtain where we have used the trace of three M matrices, equation (B.11), and the integral can be evaluated immediately observing that the integrand is actually constant. We remark that such a contribution is purely imaginary and therefore it mixes with the other imaginary three-loop corrections that are due to framing, though this part is framing independent. At this order this is the only imaginary contribution not arising from framing. In Section 3.4 we provide an additional interpretation of this term.

Multiple windings
Multiple windings introduce an overall m-dependent factor for each diagram, which only depends on the combinatorics of the insertions on the Wilson loops. Such factors can be computed recursively in an algorithmic manner as shown in [33]. The strategy involves simplifying iteratively the integration contours until landing on integrals which can be immediately computed in terms of single winding ones. This translates into a system of recursion relations, which supplied with an initial condition, that is the value of the integrals at m = 1, can be solved exactly. This procedure can be applied at any loop order and increases in complexity with the number of insertions on the Wilson contour. In some cases a computer implementation becomes necessary. For one-and two-loop diagrams we do not report the explicit computation, rather we state the final result in Section 3.3. Instead, we present some details of such a procedure for the three-loop diagrams in Figure 1, in particular stressing what is new compared to the single winding case.
The main difference that arises at multiple winding is the following. As stated in the previous section, for single winding only planar corrections contribute since the non-planar configurations, being non-contractible, vanish identically for the Alvarez-Labastida theorem. For multiple winding, instead, it is no longer guaranteed that non-contractible multiple winding diagrams do not contribute, since contractible single winding integrals can appear when resolving the multiply wound contours according to the procedure described above. Therefore, along the calculation we have to take into account all possible planar and non-planar configurations.
To better illustrate this point, we begin by considering the case of the factorized diagram in Figure 1e. As in the single winding case, it is convenient to complete the sum of planar configurations with the corresponding non-planar ones in order to reduce the three-loop structure to the product of lower-order ones. This implies adding and subtracting the integrals corresponding to the crossed contributions multiplied by the same color factors of the planar ones, in such a way to symmetrize the integral. For the symmetrized part we then obtain (the double line contour stands for multiple winding) For the non-planar crossed contributions, using the algorithm of [33], a recursive relation yields pictorially where the graphs on the right-hand-side represent integrals not diagrams, since the corresponding color and combinatorial factors have been already stripped out. We note that the term proportional to C 2 2 (R) is the relics of the symmetrization procedure in equation (3.13). Now, the combination of integrals in (3.14) is exactly the one appearing in (3.9) and evaluates iπ 3 f . Therefore, we find 15) The evaluation of diagrams 1d and 1f-1i for multiple winding is straightforward. In fact, diagrams with n insertions of fields on the contour have in general a dependence on winding through a polynomial in m 2 of degree n 2 . In particular, for diagrams with two and three insertions this boils down to a trivial m 2 factor. Consequently, at multiple winding the gauge two-point function diagram evaluates The same occurs for the interaction diagrams of Figures 1f-1h that simply acquire an overall m 2 factor compared to the single winding cousin. This is important since it allows to conclude that the addition of winding does not jeopardize the argument for the cancellation of these interaction diagrams against part of (3.16). Finally, the same m 2 overall factor arises for the triangle diagram as well. Diagrams with five and six insertions of fields along the contour, see Figures 1a-1c, are the most complicated. Their planar configurations were computed in [33] and, extended to generic representations, they read where we have summed over all possible planar permutations with the same topology, and hence the same color factor. The corresponding non-planar configurations give a non-trivial contribution since contractible planar configurations appear when decomposing the multiply wound contours. For five insertions we obtain +4 perms = 1 18 We note that in all the cases considered above the result of the recursive procedures organizes neatly in such a way that the singly wound integrals can be symmetrized and summed straightforwardly. This is not the case for the non-planar contributions with six insertions, where the partial results for the individual topologies read In all these formulae the pictures on the right-hand-side stand for integrals over the contour and not diagrams, since the color factors have been already extracted. Dealing with these integrals individually would be hard. However the non-planar sextic diagrams above combine in such a way that a symmetrized sum of integrals is reconstructed and simply evaluated. Amusingly, the final sum reads

The general three-loop result
Summing all the diagrams computed in the previous sections we obtain the three-loop expectation value for the bosonic latitude Wilson loop with parameter ν, framing f , winding number m and for a generic representation R where the Casimir invariants C 1 (R) and C 2 (R) for various representations are reported in (A.8) and (A.9). As already mentioned, the multiple windings and higher dimensional representations are not independent generalizations, rather they provide two different alternative bases of operators. Using the general expression in (3.24), we have verified explicitly that our result is in agreement with this expectation. In fact, by considering the first few windings of the Wilson loop in the fundamental representation it is easy to check that its expectation value can be obtained as a combination of single-winding operators in hook representations, according to the formula It is important to stress that this holds for any generic framing number.
In the undeformed case (ν = 1) this provides a generalization of the three-loop result of [24] to generic representations and with the inclusion of all non-planar contributions. In the case of totally symmetric and antisymmetric representations we have tested our result evaluated at f = 1 against the weak coupling expansion of the matrix models [52] (see also [53]) where the expectation values are defined in terms of the measure in (2.12). In Appendix C we supply the four-loop expansion of these matrix models up to rank-3 representations. We find perfect agreement with the perturbative result (3.24) at ν = 1, thus providing a strong mutual check of the correctness of (3.24) and of the localization prediction. Specifying result (3.24) to the fundamental representation of the gauge group U (N 1 ) we obtain and for single winding it boils down to Finally in the ABJM theory (N 1 = N 2 ≡ N ) this reads The crucial observation is that in the planar limit, setting f = ν this expression coincides with prediction (2.35). Moreover, it verifies relation (2.32) that was required for the consistency of expression (2.31) for the 1/2-BPS Bremsstrahlung function. However, we stress that in the more general case N 1 = N 2 identity (2.32) is no longer valid, as one can easily check from (3.28). Consequently, all its implications discussed in Section 2 stop working.
We conclude observing that the expectation value of the fermionic latitude can in principle be obtained from (3.28) by using the cohomological equivalence (2.21) conjectured in [24]. This applies only for the putative framing f = ν, while a generic framing would require a full-fledged computation of the fermionic diagrams at three loops, which is currently not known. At this special value of framing, using (2.21) we obtain the following prediction (with single winding, for simplicity) where the normalization factor R has been defined in (2.8). The first two terms of this expansion reproduce the perturbative expansion of [15] 6 .

The imaginary term at framing zero and the Bremsstrahlung functions
Imaginary contributions to the expectation values of Wilson loops in ABJM are usually associated to framing. Hence, the appearance of an imaginary term at three loops in the expectation value of the latitude Wilson loop at framing zero is a bit surprising, but not inconceivable. In fact, the operator in (2.3) that we are considering does not possess a definite hermiticity property, which would enforce the reality of its expectation value. Still, its appearance poses a question concerning the relation between the latitude and the Bremsstrahlung function associated to the internal angle of the generalized cusp. In the original proposal [15,22], B θ 1/6 was prescribed to be equivalent to where the latitude Wilson loop at framing zero on the right-hand-side was understood to be real. According to our findings it is not. This would induce an imaginary contribution in the Bremsstrahlung function that cannot be there, as explicitly checked by the computation in [39,40].
In order to resolve this tension we go back to the derivation of (3.31) in [22] and point out where this subtlety kicks in.
Using definition (2.25), the Bremsstrahlung function associated to θ is first determined by explicitly taking the derivatives with respect to θ of the generalized cusp Wilson loop. This is identified with the integral of two-point functions of operators constructed with the C 1 , C 2 fields [22] where the double bracket denotes the (normalized) correlation function of local operators at positions t 1 and t 2 on the 1/6-BPS Wilson straight line. Such a non-local operator defines a one-dimensional superconformal defect and the first two-point function above is fixed by conformal symmetry to possess the form where the coefficient γ encapsulates the quantum corrections and is ultimately proportional to the Bremsstrahlung function. The line configuration can be mapped to a circle via a conformal transformation, where the two-point function takes the form (3.34) in terms of the angles τ 1 and τ 2 on the circle. Conformal invariance of the theory assures that the γ factors in (3.33) and (3.34) are the same. The second correlation function in (3.32) is simply obtained from the first by exchanging t 1 and t 2 on the line or, equivalently, τ 1 and τ 2 on the circle. Since the operators obey bosonic statistics, we would be led to conclude that the two correlation functions are identical. However, expressions (3.33) and (3.34) are correct only at non-coincident points, whereas close to the singularity at t 1 = t 2 (or τ 1 = τ 2 ) a suitable regularization is required (for instance by the addition of contact terms) that might introduce parity-odd corrections. Since in (3.32) we are integrating over t 1 , t 2 this regularization can have sizable effects and, especially, lead to different results for the two integrals. Therefore, in the following we treat the two correlation functions as different objects.
For the latitude Wilson loop at framing zero an identity analogous to (3.32) reads where the exponentials in the integrand arise from the e iτ factors in the latitude operator (2.3). The correlation functions on the right-hand-side are on the 1/6-BPS circle, as we have set ν = 1 after taking the derivative. We see that keeping the correlation functions in (3.35) distinct an imaginary part arises, proportional to the antisymmetric combination of the two. For the purpose of computing the Bremsstrahlung function in terms of the latitude Wilson loop, we ascertain from (3.32) that only the symmetric combination of the correlation functions is relevant. This is equivalent to taking the real part in (3.35), which on the left-hand-side amounts to enforcing the modulus of the latitude expectation value. The precise coefficient of the relation between the Bremsstrahlung function and the latitude is computed by performing the integral in (3.35), after plugging the generic form of the correlation function (3.34) on the circle and equating the parameters γ appearing in (3.34) and (3.33). The steps are the same as in the original derivation of [22]. The final result reads We can thus put such a prediction which was anticipated in (2.27) on firmer grounds, from a first principles derivation 7 . Finally, we can explicitly verify the emergence of an imaginary contribution at three loops in the latitude (which in that context comes from the triangle diagram 1i), as arising from the imaginary part of the right-hand-side of (3.35), of the form + 2 path ordered perms (3.38) and performing exchanges in the integration variables, one can prove that this contribution precisely reconstructs the integrand of (3.12). This automatically verifies the equality in (3.35) and explains the imaginary contribution of the latitude from the two-point function perspective. This result hints at the fact that the two-point functions in (3.37) are actually distinct quantum mechanically due to the necessity of regularizing them at coincident points, as already discussed. In order to better clarify this point, we alternatively compute the one-loop two-point functions first and then plug them into (3.37). Calculating such a contribution on the circle is not straightforward, while it is immediate on the line. This computation reveals that indeed at one loop the integrands of the path-ordered correlation functions in (3.37) are opposite, due to R-symmetry index algebra and the properties of the matrix M = diag (−1, 1, −1, 1). However, after integrating over the insertion of the bi-scalar field along the line, the result can be shown to vanish. On the one hand this finding is in line with the symmetry expectations on the two-point functions, on the other hand it is seemingly in contradiction with the non-vanishing result obtained above in (3.38). This puzzle is explained observing that the integral in (3.37) is actually divergent and therefore the insertion of a vanishing integrand should be handled with care. One possibility consists in computing the one-loop two-point functions on the line in dimensional regularization where consistently we obtain an order correction. The other correlation function yields the opposite result, as stated above. We can plug this into (3.37), though this requires performing a conformal map that is not justified in non-integer dimensions. Ignoring this objection and pushing the computation ahead we ascertain that the resulting integral in (3.37) develops a pole in the regulator , thus exposing a finite contribution out of (3.39). We ultimately find which precisely agrees with the derivative of the latitude expectation value at three loops and framing zero, that is basically the triangle diagram (3.11). In conclusion, we have detected the interesting phenomenon that from the 1/6-BPS defect CFT perspective the triangle diagram in the latitude deformation arises from an anomalous behavior of the relevant two-point functions on the defect.

The matrix model
As reviewed in Section 2, a matrix model prescription for computing Wilson loops in ABJ(M) theory exists for the ν = 1 case [8], while it is still lacking for more general latitude operators. In this section we make a first attempt to fill this gap by proposing a matrix model to compute latitude Wilson loops. We then discuss some consistency checks to support our proposal. The bosonic latitude Wilson loop is partially supersymmetric and its preserved supercharges are given in equation (2.9) for ω 1 = ω 4 = 0. Consequently, and in parallel to the analogous situation in N = 4 SYM, it might be possible to compute its expectation value exactly using localization techniques. However, while in the fourdimensional case the latitude expectation value is obtained from the undeformed one by a simple rescaling of the coupling constant [7,[54][55][56][57][58][59][60][61], this is no longer true for ABJ(M), as already appears from the perturbative results of the previous section. Therefore, we expect a significant modification of the matrix model to take place as a result of the localization process. This program would involve localizing the ABJ(M) theory on the three-sphere using any of the supercharges in (2.9) preserved by the bosonic latitude. In particular, since we expect the localization procedure to be consistent with the cohomological equivalence (2.10), we should use the linear combination of supercharges (2.11). As already noticed in [15], this supercharge is non-chiral and differs in nature from those considered in the original analysis of [8]. Therefore, the generalization of the localization procedure of [8] to the latitude Wilson loop is not straightforward.
We are not going to pursue this direction here, rather we conjecture directly a matrix model that computes the latitude expectation value exactly, consistently with the perturbative results already available.
The idea is to start from the matrix model average (2.12) computing the expectation value of 1/6-BPS Wilson loops and try to deform it by introducing a suitable dependence on the ν parameter. As a route guidance we use the proposal (2.28) on the θ-Bremsstrahlung function, which in turn requires the ν-derivative of the matrix model to satisfy identity (2.29).
A natural way to satisfy this condition consists in requiring that the latitude Wilson loop is computed by inserting the operator Tr e 2π √ νλ into a matrix model which is symmetric under the inversion ν ↔ 1/ν. In fact, in this way taking the derivative with respect to ν evaluated at ν = 1, the ν dependence of the matrix model measure plays no role and the only non-trivial contribution comes from the derivative acting on the operator insertion. The result is a matrix model where the integrand, being evaluated at ν = 1, corresponds to the well-known ABJ(M) matrix model, except for the operator insertion ∂ ν Tr e 2π √ νλ | ν=1 which can be formally identified with the m-derivative of a multiply wound 1/6-BPS Wilson loop evaluated at m = 1, with m ≡ √ ν. In particular, trading ∂ ν with ∂ m provides the correct 1/2 factor appearing in (2.29). Such an argument is reminiscent of the one proposed in [26], but somehow with a reverse logic. In that case, for the ABJM theory a supersymmetric Wilson loop on a squashed sphere was considered, whose matrix model is invariant under the inversion of the squashing parameter [62]. This was used to argue that the derivative of this Wilson loop expectation value with respect to the squashing parameter b, evaluated at b = 1, could be traded with the derivative of the multiply wound 1/6-BPS Wilson loop with respect to the winding number m.
Driven by this discussion we are led to propose the following matrix model average for the expectation value of a multiply wound latitude Wilson loop where the average is evaluated and normalized using the matrix model partition function Similarly, Ŵ m B (ν) corresponds to the insertion of 1 1≤i≤N 2 e 2π m √ ν µ i . According to the discussion above, such a matrix model should arise from a suitable localization of the ABJ(M) theory. This is the simplest non-trivial deformation of the matrix model (2.12) that lands back on the usual expression at ν = 1, and whose kernel is symmetric under ν ↔ 1/ν. The precise dependence on ν in the hyperbolic functions and in the operator insertion is then fixed via comparison with the perturbative results. Indeed, we can evaluate this expression by expanding it at weak coupling. The main result of this analysis is that we recover precisely the expectation value (3.27) for the multiply wound latitude in the fundamental representation evaluated at framing ν. We stress that the agreement with the perturbative result holds separately for all different color structures at generically unequal and finite N 1 and N 2 and for generic winding number m. Moreover, in the N 1 N 2 approximation, the matrix model reconstructs the pure Chern-Simons result at framing ν, and in the specular N 2 N 1 limit it reproduces the expected behavior of [51].
In addition, we have considered the extension of the matrix model average defined in (4.1) to higher dimensional representations. We have explicitly ascertained that applying the prescriptions (3.26) (with an extra √ ν factor in the exponents) for the first few totally symmetric and antisymmetric representations, we do reproduce the corresponding perturbative results (3.24) for bosonic latitude Wilson loops. This lead us to conjecture, that the matrix model (4.1) computes also the expectation value of latitude Wilson loops in higher dimensional representations, upon applying the same prescriptions (3.26), as in the undeformed case. It is remarkable that the agreement works only at the specific choice of framing f = ν. On the one hand this does not come as a surprise. In fact this occurs already for the 1/6-BPS Wilson loop where localization implies non-trivial framing f = 1, which is the scheme compatible with the cohomological equivalence between bosonic and fermionic Wilson loops. For the latitude Wilson loop it is then highly suggestive that the agreement manifests at f = ν, since this is precisely the value at which the conjectured cohomological equivalence with the fermionic Wilson loop holds (see equation (2.18)). Since we expect a putative matrix model average to be able to compute both the bosonic and the fermionic operators, it is then natural that the matrix model indeed provides the result at framing ν.
Although it might sound a bit weird to consider non-integer framing, we note that at the level of the matrix model continuing the framing from an integer value to a generically real number is perfectly legitimate and is also a common occurrence (despite usually only for rational values), for instance when computing torus knot invariants in Chern-Simons theory [45].
We can finally draw a parallel with the four-dimensional N = 4 SYM theory. In that case, it is easy to realize that applying an analogous deformation procedure on the Gaussian matrix model computing the expectation value of Wilson loops in that theory, we reproduce the latitude operators. In fact, the ν dependent deformations in the matrix model measure cancel out completely and one is left with a modified operator insertion exhibiting an additional ν factor at the exponent. This eventually provides the coupling constant rescaling that characterizes the expectation value of latitude Wilson loops in N = 4 SYM.

Properties and relations with other matrix models
Supported by this first evidence on the correctness of our proposal, we devote the rest of this section to a discussion of its main properties and its possible interpretation.
First, we note that a striking similarity exists between expression (4.2) and the kind of matrix model emerging as the result of the so-called symplectic or SL(2, Z) transformation of the sphere partition function of pure Chern-Simons theory [45]. This depends on two coprime integer parameters (P, Q) and has been argued to compute torus knot invariants [45] (see also [63][64][65] for more references on torus knot invariants), using the celebrated relation between knots and Chern-Simons theory [3]. In our case, the symplectic transformation is rather performed on the supermatrix model (or equivalently on the lens space S 3 /Z 2 partition function [66,67]) and the result reads Consequently, the Wilson loop averages computed with this matrix model where s p R is the supersymmetric Schur function for the partition p R associated to the representation R, should arguably yield torus knot invariants of U (N 1 |N 2 ) Chern-Simons theory [68] (they are a generalization of colored HOMFLY polynomials [69] to U (N 1 |N 2 )/lens space). This interpretation in terms of knot invariants works only for coprime P, Q integers, this condition ensuring that the contour closes on the two-torus. The case we are considering is a (perfectly sensible) generalization of the torus knot matrix model to non-coprime integer values of the parameters, being P = √ ν and Q = 1/ √ ν. Hence, in general the latitude Wilson loop is not really computing knot invariants, as the corresponding contour does not close, rather it wraps the two-torus densely. An exception arises in the degenerate situations in which the torus reduces to one of the two cycles of length 2πP or 2πQ, respectively. In these two cases the factor in the Wilson loop operator of (4.1) is precisely of the form required for the closure of the contour (with P = √ ν being the length of the circle).
The matrix model in (4.2) can also be obtained by localizing the N = 2 U (N 1 +N 2 ) Chern-Simons theory on the squashed lens space S 3 √ ν /Z 2 with squashing parameter √ ν [62,70,71]. Selecting the particular vacuum that breaks the gauge group to U (N 1 ) × U (N 2 ), and then continuing it to the supermatrix version as done for the ABJ(M) matrix model [72][73][74], we land on (4.2). In fact, considering Chern-Simons theory with supergroup U (N 1 |N 2 ) is (at least perturbatively) equivalent to the lens space interpretation by rewriting its matrix model in the two-cut form and taking the analytic continuation N 2 → −N 2 .
After integrating out the N = 2 auxiliary fields one expects to recover the pure Chern-Simons observables and hence possibly compute knot invariants. In [75] this procedure was indeed applied to N = 2 Chern-Simons theory on a squashed sphere, which was observed to yield torus knot invariants of pure Chern-Simons on the threesphere. Our matrix model is formally a particular case of N = 2 Chern-Simons on a squashed lens space. The analysis of [75] then suggests that this matrix model should compute torus knot invariants on the lens space RP 3 . This is indeed consistent with the identification of the latitude matrix model and that for torus knot invariants, as described above.
In cauda venenum, we conclude with few critical arguments on our proposal. In particular we discuss whether the matrix model (4.1) can be interpreted as the result of localizing the ABJ(M) theory on S 3 , with the insertion of a ν-dependent operator corresponding to the latitude Wilson loop.
In principle, if this were the case one might expect the matrix model average to be computed with the ordinary ν-independent measure appearing in the ABJ(M) partition function. Instead, in our proposal (4.1) the kernel depends explicitly on ν. However, as already stated, if we require the localization procedure to be compatible with the cohomological equivalence, the path integral should be localized with the supercharge (2.11). Since the latter exhibits an explicit dependence on ν it might reasonably lead to the conjectured ν-dependence of the matrix model in (4.2).
If this interpretation is correct, then (4.2) itself should, after integration, give rise to the usual, ν-independent, result for the ABJ(M) partition function. In fact, in the N 1 = N 2 case, expression (4.2) can be rearranged in such a way that the ν-dependence disappears completely and it ends up coinciding with the ABJM partition function on S 3 . More details on this computation can be found in the next subsection.
Instead, for the most general N 1 = N 2 case a non-trivial ν dependence survives in the phase of the partition function (see equation (4.7)), whose appearance has been ascribed to a Chern-Simons framing anomaly in literature [76,77]. This leads to the conclusion that the deformation affects the partition function only in its somewhat unphysical part, whereas its modulus is ν independent. However, this is not a totally satisfactory explanation yet.
A more general and cautious attitude could be to assume that expression (4.1) is only a possible convenient way of rewriting the matrix model obtained properly via localization, which could arise as the result of applying some identities at the level of the matrix model average for the latitude Wilson loop. However, while we have verified that it is possible to perform a transformation that brings us back to the ordinary ABJ(M) partition function, the form of the resulting operator insertion is not very enlightening.
In conclusion, we do not have a definite clear explanation of whether and how the matrix model (4.1) could arise by performing an honest derivation of the latitude expectation value via localization, although there are some reassuring indications at least for the ABJM case. Nevertheless, supported by the striking agreement with the perturbative computation at weak coupling, we assume as a working hypothesis that the proposal in (4.1) correctly reproduces the latitude expectation value. Equipped with this tool, in Section 5 we perform a study of the matrix model average at strong coupling, where very little information is known on latitude Wilson loops [22].

Reformulation as a Fermi gas
The ABJ(M) matrix model can be reformulated in terms of a Fermi gas (see [78] for the original derivation in ABJM theory and [79][80][81][82] for its generalization to the ABJ model). This perspective provides a powerful tool for expanding systematically the partition function and Wilson loop observables in powers of 1/N at strong coupling, by using statistical mechanics technology. In this section we point out that the proposed matrix model for the latitude Wilson loop can also be given such an interpretation, paving the way for its study in the type IIA string and M-theory regimes. For simplicity we restrict the analysis to the ABJM slice, The partition function The crucial property that streamlines the Fermi gas reformulation is the Cauchy identity, which we present in a form that is suitable for our purposes Here the final sum is over all the permutations of the N eigenvalues and r is an arbitrary parameter.
We split the integrand of (4.2) into two combinations of hyperbolic functions with arguments containing the factors √ ν and 1/ √ ν respectively, and apply the Cauchy identity separately by choosing r in (4.5) appropriately. This procedure is similar for instance to the one used in [78] for N = 3 quiver models. After few algebraic steps we end up with 8 We thus see that the dependence on ν drops completely and the partition function lands on the same expression as for ABJM theory. We stress that these last steps are valid only for N 1 = N 2 . For different ranks of the gauge groups the starting point (4.5) must be replaced by a generalization of the Cauchy determinant identity discussed in [79][80][81]. Then if we repeat the steps leading from (4.5) to (4.6), we can isolate and evaluate the ν−dependent part of the partition function. We find The only effect of the deformed measure consists in altering the original phase of the ABJ partition function obtained in [21] by a trivial ν−dependent multiplicative factor. Going back to the ABJM case, it was observed in [78] that (4.6) can be formally interpreted as the canonical partition function of an ideal Fermi gas of N particles ρ(y a , y σ(a) ) ≡ Tr ρ (4.8) with a density matrix ρ completely factorized into non-trivial one-particle density matrices ρ(y a , y σ(a) ) = 1 We can then define the corresponding one-particle quantum Hamiltonian ρ = e −Ĥ such that y 1 |ρ|y 2 ≡ ρ(y 1 , y 2 ) (4.10) Using the explicit expression for ρ in (4.9) the Hamiltonian can be written as e −Ĥ = e −U (q) e −T (p) , U (q) = log 2 coshq 2 , T (p) = log 2 coshp 2 (4.11) in terms of a non-standard kinetic term T (p) and a potential U (q), whereq andp are canonically conjugate operators satisfying [q,p] = i , = 2πk (4.12) Accordingly, the quantum average of a single-particle operatorÔ ≡ O(q,p) is given by Introducing Wilson loops When a Wilson loop operator of the form N i=1 e 2πm √ νλ i , relevant for computing the bosonic latitude, is inserted into the matrix model average we can still use the Cauchy identity and perform the same steps as those leading to (4.6). In the case of a single winding operator (m = 1) the result we obtain reads As expected, the ν dependence does not drop any longer. In particular, it appears not only in the overall phase factor, but also in an interacting piece. This can be interpreted as the statement that for this Wilson loop the framing factor is in general non-trivial and related to the ν parameter. This is a generalization of what happens already for the undeformed operator (ν = 1), for which the matrix model computes the average at framing one. Within the Fermi gas approach reviewed above, the Wilson loop expectation value (4.14) maps to the quantum average (4.13) of the ν-dependent one-body operator with the ABJM density operator defined in (4.10), (4.11). We observe that for ν = 1 this operator reduces to that corresponding to the undeformed Wilson loop [25]. In general, the presence of the ν factor unbalances the (q,p) symmetry, which manifests in the ν = 1 case. We also point out that the insertion of operator (4.15) implies a normalization for the Wilson loop expectation value that differs by an overall N factor from the one used at weak coupling in Section 3.3. Here we find convenient to use this normalization for a better comparison with the formulae of [25].
Equations (4.14) and (4.15) can be generalized to the case of a multiply wound Wilson loop. In particular, its expectation value corresponds to the average of the single-body operatorÔ m (ν) = e m k (q + νp) (4.16)

A peculiarity of the matrix model
We conclude this section by highlighting a peculiar property of the matrix model average (4.1), valid for the ABJM theory (N 1 = N 2 ≡ N ). Precisely, we claim that its real part is independent of ν This property can be proven as follows. We consider the Wilson loop expectation value as in (4.14), after the application of the Cauchy identity. We take the real part of this expression and its derivative with respect to ν at the level of the integrand, obtaining the following expression Next we work directly on the integrand and prove that (4.17) holds before any integration.
The permutations in the symmetric group S N can be divided into those which are idempotent {σ ∈ S N |σ 2 = 1} and those which are not. The former are those constructed as products of cycles of maximal length 2, whereas the latter contain at least one cycle of length greater than 2. For permutations belonging to the first group we find that the sum vanishes identically since the summand is antisymmetric with respect to the exchange y c ↔ y σ(c) . In fact, due to this property a given term c =c in the sum either vanishes if σ(c) =c or, if σ(c) =c , it will be canceled by an opposite contribution for c =c . Next we consider permutations which are not idempotent, so that they do not coincide with their inverse. This set can be divided into pairs (σ, σ −1 ), which have the same signature and span the whole set. They give rise to the same factor N a=1 cosh in (4.18). Hence, restricting the sum in (4.18) to such a pair of permutations, we can focus on the two contributions since, by construction, σ −1 (c ) =c and consequently the latter term is equal and opposite to the former. Such a cancellation extends pairwise to all the terms in the sum (4.20), and therefore the terms in (4.18) associated to a given permutation and its inverse cancel completely. This argument in turn extends to all permutations which are not idempotent. Consequently, the whole expression (4.18) evaluates to zero, as claimed.
As an important corollary of identities (2.29) and (4.17) we find that This proves a property of both the multiply wound 1/6-BPS average and the latitude Wilson loops, which in particular allowed the steps leading to (2.33).

The matrix model at strong coupling
For the ABJM slice N 1 = N 2 = N , in this section we provide the expansion of the matrix model average (4.14) that computes the expectation value of the bosonic latitude Wilson loop at strong coupling. Automatically, this also provides the average of the fermionic operator at strong coupling, via the cohomological equivalence (2.21).
For simplicity, we confine the treatment to single winding operators (m = 1), and point out how to extend the calculation to the more general case of multiply wound Wilson loops, if need be. We work at large N but we do not restrict to the planar limit. We adopt the Fermi gas approach reviewed in the previous section, as this method has the virtue of granting a systematic control on both quantum corrections in the coupling and on the genus expansion around the large N limit 9 .
As recalled above, the matrix model (4.14) can be reformulated in terms of a onedimensional ideal (non-interacting) quantum gas of particles with Fermi statistics. In this setting the Chern-Simons level k plays the role of the Planck constant, equation (4.12), and the number of colors N corresponds to the number of particles. Therefore, large N is equivalent to the thermodynamic limit of the gas, whereas the expansion encoding quantum corrections corresponds to an expansion at small k. Consequently, the Fermi gas approach is suitable for studying the latitude expectation value in the M-theory regime, N → ∞ and k fixed. We will limit our discussion to perturbative corrections, neglecting exponentially small contributions, stemming from world-sheet and membrane instantons.

Thermodynamic limit and quantum corrections
When the number of particles becomes large the canonical partition function (4.8) is hard to deal with. We then resort to the grand-canonical ensemble with the grandcanonical partition function defined as where z = e µ is the fugacity and µ the chemical potential. Accordingly, the canonical average for the one-body operator (4.15) is substituted by the grand-canonical average with Fermi statistics The canonical average is then retrieved by an inverse transform in µ.
The large N expansion of the ABJM model can be argued to translate into a large chemical potential and energy expansion [25,78]. Resumming the power series from expanding (5.2) in this limit, one obtains [25] with Θ being the Heaviside function that defines the Fermi surface.
In order to keep track of quantum corrections for the operator average we use a semi-classical approach and look for a convenient way to expand (5.4) in powers of . As discussed in [25], this is easily accomplished using the phase space formalism of Quantum Mechanics. This amounts to trading operators with their Wigner transform, according to and operator products with -products In this formalism the trace is rewritten as an integral over the (p, q) coordinates TrÂ = dp dq 2π A W (q, p) (5.7) Applying the Wigner transform to the one-particle HamiltonianĤ defined in (4.10), we can expand any operator f (Ĥ) in powers of (Ĥ − H W (q, p)). Its semiclassical expansion is then obtained by taking the Wigner transform where the so-called Wigner-Kirkwood functions G r have an expansion of the form Applying this formalism to the distribution in (5.4), its semi-classical expansion reads The value H W (q, p) = µ defines the Fermi surface. In our case, from the explicit expression (4.11) of the one-body Hamiltonian we realize that H W is explicitly given by e −H W = e −U (q) e −T (p) , U (q) = log 2 cosh q 2 , T (p) = log 2 cosh p 2 (5.11) whereas for a generic operator of the formÔ a,b ≡ e (aq+bp)/k we find Equation (5.10) with entries (5.11) and (5.12) with a = 1, b = ν, are the key ingredients to obtain the expectation value W B (ν) ν at strong coupling from prescription (5.3). We devote the rest of this section to its explicit evaluation.

Expansion of Fermi gas at strong coupling
In order to evaluate (5.10) we closely follow the treatment of [25], generalizing the form of the operator's Wigner transform as in (5.12). The first step requires deriving the expression for the grand-canonical partition function of the Fermi gas. Since by construction this coincides with that of the undeformed case, we can read it from [25] Next, we concentrate on evaluating the occupation number distribution 14) The integrals are over the Fermi region and surface, whose picture is given in Figure  2a.
As suggested in [25], it is convenient to divide the Fermi surface into four sectors. The points where the separation occurs are where O(e −µ ) stands for exponentially suppressed terms. The logic behind such a separation is that along the curve bounding the regions, we have that alternatively |p| > µ or |q| > µ (5.16) This in turn means that exponentially small corrections in p and q are bounded by exponentially small corrections in µ and hence they can be neglected. Precisely, in regions I,III: e −|p| < e −µ T (p) ∼ p 2 in regions II,IV: e −|q| < e −µ U (q) ∼ q 2 (5.17) where the approximations are correct up to exponentially small terms. Due to the invariance of the Hamiltonian under p ↔ q, the different regions can be obtained from one another by exchanging the canonical coordinates. Upon the insertion of the O a,b operator this is equivalent to computing the contribution with the insertion of O b,a . Consequently, the idea is to explicitly compute the integration of O a,b along the curve in regions I and adding to it the contributions form the other regions obtained by changing signs and permuting labels. However, this procedure would overcount the contribution from the square (|q| < q * , |p| < p * ), which then needs to be subtracted. Naming this extra contribution n For the other contributions, we recall that for a = b, e.g. for the 1/6-BPS Wilson loop, the domain of integration where the operator does not get exponentially suppressed coincides with regions I and II only and the computation can be limited to those. For a = b this is not necessarily true. In particular, choosing a = 1, b = ν and assuming without loss of generality that 0 < ν < 1, we expect contributions from region III (see Figure 2b) not to be entirely exponentially suppressed. Working in regions I and III we can use the approximated expression for T (p) in (5.17). Following the steps of [25], the first term in (5.14) can be integrated in p to get an integral along the Fermi surface This expression is indeed p-independent as a consequence of (5.17). Solving for the δ function in the second term of (5.14), and adding it to (5.20) we obtain the resummed expression where we have used the definition of the Wigner-Kirkwood corrections (5.8). This object has been computed in [25] using the reduced Hamiltonian for regions I, III Plugging (5.23) into (5.22) the complete expression for n O a,b in region I reads The analogous expression for region III is obtained from (5.24) by sending b → −b.
One can now perform the change of variables u = exp (q/k) to put the integral in the form where the lower limit of integration has been lowered down to 0 at the affordable price of introducing spurious exponentially small corrections that we are neglecting anyway. Finally, the integralĨ a,b can be evaluated in full generality in terms of a hypergeometric functioñ Leading exponential asymptotics We now perform the large µ asymptotics of (5.27), discarding exponentially subleading contributions.
To this end we first invert the arguments of the hypergeometric function above using the following general identity and then expand each term in a power series. Since u * is exponentially large, the first term in the expansion is sufficient for retaining only the exponentially leading contributions. We will restrict here to the cases of interest, which areĨ 1,ν ,Ĩ ν,1 andĨ 1,−ν with 0 < ν < 1, although the same analysis could be carried out for generic a and b. In particular, this would allow to take into account the case of multiply wound latitude Wilson loops, for which we should set a = m and b = mν.
Restricting to the single winding operator, the large µ asymptotics in the various regions reads We note that the contribution from region III can be neglected, even though it is not exponentially suppressed (in the sense that it does not vanish exponentially for large µ). In fact, compared to the contributions from the other regions it is subdominant and bounded from above by the subleading exponential e µ/k in the whole 0 < ν < 1 range. Moreover, it possesses a different behavior in the ν → 1 limit. Summing up the contributions from regions I and II and subtracting the asymptotic expansion of n This expression is well-defined in the whole physical region 0 ≤ ν ≤ 1. In fact, for ν → 0 we explicitly obtain after discarding a constant term in µ arising from the first piece in (5.31). Although the two terms in (5.31) are separately singular for ν → 1, the singularities cancel, leaving the finite remainder This expression coincides with the result of [25] for the singly wound undeformed 1/6-BPS operator.
Genus expansion Having computed the asymptotic expansion of the n O 1,ν distribution, using prescription (5.3) we can finally evaluate the expansion at strong coupling of W B (ν) at framing ν through the grand-canonical average (setting = 2πk) k csc 2π ν k Γ(ν + 1) (5.35) with Ξ given in (5.13). The result can be expressed in terms of Airy functions as in [25] 10 As in the case of the m-winding 1/6 BPS Wilson loop computed in [25], this result is missing 1/k corrections. In fact, an alternative derivation of this equation [84] reveals in the second line should be modified and would contain a non-trivial dependence on k 11 .
Introducing the string coupling g s = 2πi k (5.37) 10 We recall that here the normalization of the operator has been chosen as in [25] and differs by a factor N from that used at weak coupling in Section 3. 11 We are grateful to Kazumi Okuyama for sharing with us his results.
we can now expand (5.36) at strong coupling and in the genus series While g s > 0 terms will be not reliable due to the subtlety mentioned above, we can safely compute the genus-zero term. To this end, it is convenient to define the new variable κ throught the identity In terms of this variable the genus-zero contribution reads To conclude, we mention that by taking the complex conjugate of these expression one obtains the genus-zero contribution for the bosonic Wilson loop Ŵ B (ν) ν corresponding to the second gauge group of the ABJM theory.
The fermionic operator As already discussed, once we know the expectation values for the bosonic latitude Wilson loops we can recover that of the fermionic operator thanks to the identity (2.21). At strong coupling, the cohomological equivalence can be implemented already at the stage of the occupation number distribution. Its fermionic version then reads This is the only surviving ν-exponential behavior, due to an unforeseen cancellation. Following the steps described above, we obtain the expression of the expectation value in terms of Airy functions and its genus expansion. The first term reads while higher genus contributions would be still affected by the lack of 1/k corrections in (5.42), as inherited from the bosonic result.
As already mentioned in Section 2, in this case the ν → 0 limit is ill-defined, due to the normalization factor R = 1/(e − iπν 2 − e iπν 2 ) that drives the limit to infinity, much alike what happens for the 1/2-BPS operator at even winding numbers. A sensible result at ν = 0 can be obtained by removing the R factor and replacing it with 1.
As a last comment, we note that taking the ν → 0 limit at the level of the Airy functions allows to derive the following curious relation valid at strong coupling between the singly wound bosonic latitude Wilson loop at ν = 0 and framing zero, and the fermionic 1/2-BPS Wilson loop at framing one.

Comparison with the string prediction
Classical string configurations that are dual to the latitude operators have been discussed in [22]. The fermionic operator maps to a type IIA string configuration in the AdS 4 × CP 3 background, whose endpoints are not fixed in the internal space, but rather move along a circle in CP 3 . This accounts for the non-trivial profile of the matter couplings (2.7) arising in the field theoretical definition of the Wilson loop. The semi-classical analysis of [22] reveals that the leading exponential behavior of such a configuration scales according to Expansion (5.43) for the matrix model at strong coupling remarkably agrees with this string prediction, thus providing a further non-trivial test of the correctness of our proposal. Beyond that, result (5.43) predicts the precise normalization of this exponential as well as its quantum corrections, which call for further string theory checks.
For the bosonic latitude operator no precise dual string configuration has been determined yet and therefore our findings constitute a brand new prediction, begging for a string theory confirmation. We remark that in the undeformed case the ratio between the bosonic and fermionic operators is simply proportional to √ λ, which has been interpreted as the volume of CP 1 inside CP 3 [21], in agreement with the proposed interpretation of the 1/6-BPS Wilson loop as a string smeared over that cycle [17]. For the latitude operator, instead, the bosonic expectation value displays a more complicated structure with two exponential behaviors (see equation (5.40)), potentially suggesting that the smearing over CP 1 interpretation does not carry through this case. This is in line with the comments in [22], which seem to rule out the possibility to describe the bosonic latitude through a simple geometric smearing in the internal space. 6 A conjecture for B 1/2 in the ABJ theory In this section we discuss a possible generalization of the Bremsstrahlung functions in (2.24) and (2.25) to the case of U (N 1 ) k × U (N 2 ) −k ABJ theory.
In general, for N 1 = N 2 less is known about the B-functions compared to the N 1 = N 2 case. Perturbative results for all the Bremsstrahlung functions exist from a direct evaluation of the corresponding cusped Wilson loop. Based on the two-loop result, in [38] it was argued that in the ABJ case the cusped Wilson loop has a doubleexponentiation structure. A different exponentiation structure still involving two terms has been further derived in [85] by resumming ladder diagrams to all orders. In particular, the result of [38] seems to point towards a non-trivial B 1/2 at first order, while having a two-loop vanishing contribution from both the exponents. It turns out instead, that the B ϕ 1/6 and B θ 1/6 expansions start at order two in the couplings. In [39,40] the evaluation of B θ 1/6 has been pushed to four loops by computing the corresponding cusped bosonic Wilson loop at that order. The remarkable output is that the result coincides with 1/2 a putative result for B ϕ 1/6 obtained by generalizing prescription (2.26) to the m-winding Wilson loop in ABJ whose expansion up to eighth order can be found in [33]. This property, although tested only at few perturbative orders, points towards the validity of identity (2.28) at any loop order also for N 1 = N 2 , and supports the conjecture that the general prescriptions (2.26, 2.27) hold also in the ABJ theory.
In fact, if we trust our matrix model as the correct prescription for computing latitude Wilson loops, we are guaranteed by construction that identity (2.29) is valid also for the ABJ theory. This identity, together with (2.28), conspires to sustain the conjecture that prescriptions (2.26, 2.27) have an obvious generalization to the ABJ case.
Therefore, supported by these reassuring facts, we use the prescription in (2.27) to make a prediction for B 1/2 in the N 1 = N 2 case. Inserting there the explicit matrix model expansion at weak coupling (here also including the fifth order term) we find as a sensible proposal for the fermionic Bremsstrahlung function in the ABJ theory.
There are of course some non-trivial aspects in this proposal that should be better understood. First of all, while the vanishing of the two-loop contribution is consistent with the result of [38], the one-loop term does not coincide and seems to suggest a different exponentiation structure. Moreover, the prescription of taking the modulus seems to be crucial to recover the vanishing of the second order term, a fact that needs a deeper understanding.
Further perturbative data, obtained from generalizing the three-loop computation of [27] to the ABJ case, would surely give more insights on the exponentiation procedure of the cusped Wilson loop, leading possibly to a check of conjecture (6.2). We leave this open problem for the future.

A ABJ(M) theory
In this section we summarize basic notions on the quiver U(N 1 ) k × U(N 2 ) −k ABJ(M) theory. Its field content includes two gauge fields (A µ ) i j and (Â µ )ˆiĵ belonging respectively to the adjoint of U(N 1 ) and U(N 2 ), and matter scalar fields (C I ) iĵ and (C I )ˆi j plus their fermionic superpartners (ψ I ) iĵ and (ψ I )ˆi j . The fields (C I ,ψ I ) transform in the (N 1 ,N 2 ) of the gauge group (small latin indices) while the pair (C I , ψ I ) belongs to the representation (N 1 , N 2 ). Index I = 1, 2, 3, 4 labels the fundamental representation of the SU (4) R-symmetry group.
After gauge fixing, the action reads with (c, c) and (c,ĉ) being the ghosts. We use spinor and group conventions of [27]. The invariant SU (4) -tensors are defined as 1234 = 1234 = 1 and the covariant derivatives are given by

A.1 Color conventions
The The structure constant are then defined by We perform computations associating to every generator T A in the given representation R of U (N 1 ) a matrix R j i with indices i, j = 1 . . . N 1 with commutation relation in such a way that the most generic Casimir invariant reads For the rank n totally symmetric and totally antisymmetric representations they evaluate For Hook representations with a total of m boxes and m − s boxes in the first row the quadratic Casimir invariants read

B Feynman rules
We use the Fourier transform definition In euclidean space we define the functional generator as Z ∼ e −S , with action (A.1). This gives rise to the following Feynman rules • Vector propagators in Landau gauge • Scalar propagator The one loop gauge propagators are given by The one-loop fermion propagator reads Being proportional to the difference (N 1 − N 2 ), it vanishes in the ABJM limit. color

B.1 Gauge two-point function at two loops
In Table 1 we list the non-vanishing diagrams contributing to the gluon two-loop self-energy and their properties. The color factors are relative to the diagram already inserted into the Wilson loop. The values up to order 0 are for the two-loop diagram only and a factor e 2γ E is understood.

B.2 Couplings to scalars
Up to three loop order the computation of the latitude expectation value involves the following traces of M matrices TrM (τ ) = 0 (B.9) Tr (M (τ 1 )M (τ 2 )) = 4 1 + (ν 2 − 1) sin 2 τ 1 − τ 2 2 (B.10) C Weak coupling expansion of the un-deformed matrix model In this section we expand the ABJ(M) matrix model at weak coupling and compute the expectation value of 1/6-BPS Wilson loops with m windings up to the fourth order, for generic N 1 and N 2 . This is performed by observing that every matrix model correlator with a total power of 2n eigenvalues scales as k −n . Therefore one can expand in power series the hyperbolic functions in the integrand and the exponential accounting for the operator insertion. This boils down to computing correlators in a Gaussian matrix model and for the purpose of the present expansion those listed in the appendices of [39] are sufficient. We evaluate the Wilson loop in higher rank totally symmetric and antisymmetric representations of U (N 1 ) using prescriptions (3.26). For the three lowest rank representations, up to four loops, we find for instance We can compare these results with the general three-loop expression derived from (3.24) by setting ν = 1, and framing f = 1 as required by comparison with localization predictions Selecting R = S 1 , S 2 , S 3 , A 2 and A 3 and using relations (A.8) we find perfect agreement with (C.2)-(C.6).

D B θ 1/6 at four loops for generic representations
In this section we provide the details of the computation of the θ-Bremsstrahlung function up to four loops, for generic representations of the U (N 1 ) gauge group. Most of the computation has already been addressed in [39,40], to which we refer for a more complete discussion. There the cusp has been evaluated for the fundamental representation of the gauge group. Here we extend it to general representations. This is aimed at verifying that the conjecture for the exact value of B θ 1/6 (cfr. equation (2.29)) put forward in [39,40] for operators in the fundamental representation, actually holds for any representation of the gauge group. We compute the θ-Bremsstrahlung function from a direct evaluation of the cusped 1/6-BPS Wilson loop using the definition (2.25), setting ϕ = 0. This amounts to computing the operator along a straight line, but at a non-trivial internal θ angle.
We focus only on graphs with an explicit dependence on θ. These are depicted in Figures 3 and 4. After evaluating the algebra of these diagrams using the Feynman rules in Appendix A and expressing color factors in terms of Casimir invariants using (A.7), we perform an integration-by-parts reduction of the corresponding Feynman integrals. The relevant master integrals have been evaluated in [39,40] up to the required order in .

D.1 Results for the four-loop diagrams
Here we report the results of the evaluation of the various diagrams. A common factor e −4 γ E k(4π) d/2 4 is understood. The planar topologies of Figure 3 yield  (d) =N 2 2 C 2 2 (R) 4 (π 2 (−9 + 7π 2 + 12 log 2(4 log 2 − 1)) C 2 θ ) 3 Diagrams (l)-(o) cancel pairwise and we have not shown their explicit result. The non-planar diagrams of Figure 4 read The corresponding contributions to diagram (p 1 ) in Figure 5 are obtained by multiplying these by 8 B(1 + 2 , 1) I(2, 1/2 + 3 ), where a factor of 2 stems from the two scalar propagators, a factor 4 comes from the normalization of HQET integrals and the indices of the bubble integrals are fixed by dimensional analysis.
The effect of applying the ν-derivative to the integrals in (5.26) (where for convenience we factor out (−1) −ν in the integral for region I) is to produce the following two new integralsĨ Note that for any given s these expansions stop after a finite number of terms. Plugging this expansion into (F.5) we finally obtain ∂ ν n O 1,ν ν=1 = ie 2µ k ((24 + π 2 ) k 2 − 48kµ + 48µ 2 ) 96π 2 k (F.7) If we now take the same expression (5.3) for the m-wound 1/6-BPS Wilson loop, compute its derivative with respect to m and set m = 1 we indeed reproduce (F.7) correctly. This constitutes a successful test of our procedure.
The n-th derivative In principle further checks can be performed applying an arbitrarily large number of derivatives, provided there is enough computational power. In particular, taking the n th derivative of the integrals in (5.26), evaluating their large µ asymptotic behavior and plugging the result in (5.24) we check that we obtain the same expression as from applying the n th derivative directly on the asymptotic expression