Wilson loops in 3d $\mathcal{N}=4$ SQCD from Fermi gas

We study 1/2 BPS Wilson loops in 3d $\mathcal{N}=4$ $U(N)$ Yang-Mills theory with one adjoint and $N_f$ fundamental hypermultiplets from the Fermi gas approach. By numerical fitting, we find the first few worldsheet instanton corrections to the Wilson loops with winding numbers 1, 2 and 3. We verify that our Fermi gas results are consistent with the matrix model results in the planar limit.


Introduction
Fermi gas approach [1] is a powerful technique to study large N behavior of the partition functions of various three dimensional theories. Of particular interest is the so called ABJ(M) theory [2,3] which describes the worldvolume theory on N M2-branes on the orbifold C 4 /Z k . It turns out that the partition function of ABJ(M) theory on a three sphere correctly reproduces the N 3/2 scaling of free energy expected from the holographically dual M-theory on AdS 4 × S 7 /Z k [4]. The Fermi gas formalism enables us to go further beyond this perturbative behavior, and now we have a complete understanding of the instanton corrections to the free energy coming from wrapped M2-branes on the bulk M-theory side (see [5,6] for reviews). Also, Fermi gas formalism can be applied to the computation of Wilson loops in ABJ(M) theory and many interesting properties were uncovered in the last few years [7][8][9][10][11][12]33]. However, for many theories with less supersymmetry than the ABJ(M) case, we are still lacking a detailed understanding of the large N behavior of free energy.
In [13][14][15], some progress has been made in the so called N f matrix model which describes the S 3 partition function of d = 3 N = 4 U (N ) Yang-Mills theory coupled to one adjoint and N f fundamental hypermultiplets. This theory naturally appears as the worldvolume theory of N D2-branes in the presence of N f D6-branes, which in turn is holographically dual to M-theory on AdS 4 × S 7 /Z N f in the large N limit. Here the Z N f action on S 7 is induced from the orbifold C 2 × (C 2 /Z N f ) with A N f −1 ALE singularity. By the 3d N = 4 mirror symmetry, this theory is dual to a A N f −1 quiver gauge theory [16,17]. As emphasized in [14], the N f matrix model admits not only the ordinary 't Hooft limit (type IIA limit on the bulk side) Here 1/N f plays the role of string coupling. An important consequence of the Fermi gas description is that these two limits are actually smoothly connected since the grand partition function of N f matrix model is a completely well defined quantity for any value of N f . In [15], the first few instanton corrections to the grand potential were determined in a closed form as a function of N f ; it is found that the structure of instanton corrections of N f matrix model is quite different from the ABJM case, and in particular there is no obvious connection with the topological string. However, the pole cancellation mechanism between worldsheet instantons and membrane instantons, originally found in the ABJM case [18], works also in the N f matrix model [15].
In this paper, we will consider the Wilson loops in the N f matrix model from the Fermi gas approach. We will focus on the 1/2 BPS Wilson loop on S 3 which wraps m times around the equator. By the numerical fitting, we find the first few worldsheet instanton corrections to the vacuum expectation value (VEV) of winding Wilson loops with winding number m = 1, 2, 3.
We find that our Fermi gas result is consistent with the planar VEV of winding Wilson loops in the 't Hooft limit (1.1) obtained from the resolvent of N f matrix model [14]. Our Fermi gas result suggests that there is no "pure" membrane instanton corrections to the Wilson loop VEVs except for the contributions from bound states of membrane instantons and worldsheet instantons. This is reminiscent of the instanton corrections to the 1/2 BPS Wilson loops in ABJ(M) theory [8,10].
This paper is organized as follows. In section 2, we consider VEVs of 1/2 BPS winding Wilson loops in N f matrix model from the Fermi gas approach and explain our numerical method to compute them. In section 3, we find the perturbative part of winding Wilson loops in a closed from, and in section 4 we determine the first few worldsheet instanton corrections for the Wilson loops with winding number m = 1, 2, 3. In section 5, we reproduce the perturbative part of Wilson loop in the fundamental representation from the WKB expansion of spectral traces. In section 6, we summarize the planar solution of N f matrix model obtained in [14]. We find that the planar resolvent in [14] can be vastly simplified. In section 7, we compare the Fermi gas results and the matrix model results of the Wilson loop VEVs in the planar limit and find a perfect agreement. Finally, we conclude in section 8. In appendix A we list some exact values of Wilson loop VEVs and in appendix B we mention some curious properties of Wilson loops for N f = 4.

Review of N f matrix model
First we review the S 3 partition function of d = 3 N = 4 U (N ) Yang-Mills theory with one adjoint and N f fundamental hypermultiplets. In d = 3, the gauge coupling has mass dimension 1/2 and it flows to infinity in the IR. Thus the gauge kinetic term is irrelevant in IR and the S 3 partition function of this theory is independent of the gauge coupling. This theory flows to a superconformal fixed point in the IR, which is conjectured to be holographically dual to M-theory on AdS 4 × S 7 /Z N f .
Using the supersymmetric localization, the S 3 partition function is reduced to a finite dimensional integral [19], which is dubbed the N f matrix model in [14] Using the Cauchy identity, (2.1) can be rewritten as a partition function of N fermions on a real line where the density matrix ρ is given by It is also useful to express ρ(x, y) as a matrix element of the quantum mechanical operator ρ ρ(x, y) = x| ρ|y , with x and p being the coordinate and momentum of a fermion obeying the canonical commutation relation As discussed in [1], it is more convenient to consider the grand canonical ensemble by summing over N with fugacity κ = e µ . It turns out that the grand partition function is written as a Fredholm determinant which in turn is physically interpreted as a grand canonical ensemble of ideal Fermi gas [1]. The large N behavior of canonical partition function Z(N, N f ) can be deduced from the large µ behavior of the modified grand potential J(µ, N f ), which is related to the grand partition function by [18] Then the canonical partition function is recovered from J(µ, N f ) by the integral transformation where C is a contour on the µ-plane running from e πi 3 ∞ to e − πi 3 ∞. The modified grand potential J(µ, N f ) can be decomposed into several pieces: The first term of (2.9) is the perturbative part where C, B, and A are given by [13][14][15] Here A c (N f ) is a certain resummation of the constant map contributions of topological string [15,20,21] A As shown in [15], A c (N f ) can be evaluated in a closed form for integer N f . In particular, for N f = 1 we find (2.9) denote the worldsheet instanton and membrane instanton corrections, respectively, while the last term J bound (µ, N f ) in (2.9) represents the contributions from the bound states of worldsheet instantons and membrane instantons. In [15], the first few instanton corrections are determined: the worldsheet instantons are given by and the membrane 1-instanton is given by At present, we do not know the exact form of the bound states in the N f matrix model. In the large N limit, the integral (2.8) can be evaluated by the saddle point approximation, where the saddle point value µ * of the chemical potential is given by Plugging this value µ * into the perturbative part of grand potential, we recover the N 3/2 behavior of free energy Also, the free energy receives instanton corrections of order where σ is one of the three scalar fields in the N = 4 vectormultiplet and x µ (s) parametrizes the equator of S 3 . The VEV of such BPS Wilson loops can be reduced to a finite dimensional integral by the supersymmetric localization [19]. Here we will focus on the Wilson loop wound around the equator of S 3 m times, which we will call the winding Wilson loop with winding number m. Now, let us consider the (un-normalized) VEV of winding Wilson loop 2 where the expectation value is defined by As in the case of ABJM theory [7,8], we can study the Wilson loop VEV of N f matrix model from the Fermi gas approach. As discussed in [8], using the relation (1 + εe mx i ) = Det 1 + κρ(1 + εe mx ) (2.25) and picking up the O(ε) term on both sides, we find that the grand canonical VEV of winding Wilson loop is written as In the following, we will mainly consider the grand canonical VEV of winding Wilson loops normalized by the grand partition function (2.27) 2 We define the Wilson loop VEV without dividing by the dimension N of representation Our definition will simplify the grand canonical VEV of Wm.

Computation of trace
To compute the grand canonical VEV of Wilson loop (2.27), we have to evaluate the trace Tr(ρ e mx ) This can be done systematically using the Tracy-Widom lemma [24]. Noticing that ρ is written as one can show the th power of ρ is given by [24] ρ where the functions ψ j (x) (j = 0, 1, · · · ) are determined recursively starting from ψ 0 (x) = 1 When N f is an integer, one can compute the exact values of VEV by closing the contour and picking up the residue of poles, as in the case of ABJM theory [18,[25][26][27]. In appendix A, we list some examples of the exact values of Wilson loop VEVs. As explained in [33], we can also compute the sequence of functions {ψ j } numerically with high precision by discrete Fourier transformations 3 . We can estimate the accuracy of our numerics by comparing with the exact values in appendix A. As an example, let us consider the Wilson loop VEV in the fundamental representation (or winding number m = 1) for N f = 6. We have computed the exact values of W (N, N f = 6) 4 up to N = 20, and in Fig. 1 we plotted the relative error On the other hand, the instanton factors (2.19) for N f = 6, N = 20, are given by From (2.33) and (2.34), one can see that our numerical computation has enough accuracy to study instanton corrections to the Wilson loop VEVs. Also, from (2.34) we observe that the membrane instanton correction is highly suppressed compared to the worldsheet instanton correction. This is a generic phenomenon in the convergence region (2.24). Thus it is difficult to study membrane instanton corrections to the Wilson loops numerically; in this paper we will mainly consider the worldsheet instanton corrections to the Wilson loop VEVs.

Odd winding number
It turns out that the grand canonical VEV in (2.27) for odd m can be rewritten in a simpler form. Let us first consider the fundamental representation. We notice that Wilson loop factor e mx for m = 1 reduces to M (x) in (2.29). Then the trace Tr(ρ M ) can be simplified using the formal operator relation Then we have Using the cyclicity of trace, we find Finally, we find that the grand canonical VEV of fundamental Wilson loop is written as This is reminiscent of the grand canonical VEV of 1/2 BPS Wilson loops in ABJM theory [8]. More generally, for odd m one can show that in a similar manner as (2.36). On the other hand, we could not find a similar expression for even winding numbers.

Perturbative part
In this section, we consider the large µ expansion of the grand canonical VEV (2.27) of winding Wilson loops. It is useful to introduce the modified version of the grand canonical As in the case of modified grand potential, W m (µ, N f ) is written as a sum of perturbative part and the exponentially suppressed instanton corrections. We find that the leading term (perturbative part) of W m (µ, N f ) in the large µ limit is given by The coefficient c m (N f ) in (3.2) can be determined as follows. As in the case of partition function (2.8), the canonical picture and the grand canonical picture of Wilson loop VEV are related by the integral transformation Then expanding the integrand of (3.3) the canonical VEV is written as a sum of Airy function and its derivatives In particular, the perturbative part of canonical VEV is given by By fitting the numerical value of W m (N, N f ) with the expression (3.6), we can determine the coefficient c m (N f ). In this way, we find that the coefficient in the perturbative part of for winding numbers m = 1, · · · , 4. For general winding number m, we conjecture We have checked this behavior numerically for m = 1, · · · , 8. In Fig. 2, we show the plot of Wilson loop VEV with m = 2 as an example. As we can see from

Instanton corrections
We can continue the numerical fitting beyond the perturbative part and fix the instanton coefficients in the expansion (3.4) and (3.5). As we will see below, we determine the first few worldsheet instanton corrections to W m (µ, N f ) for m = 1, 2, 3 in a closed form as a function of N f . We conjecture that there is no "pure" membrane instanton corrections to W m (µ, N f ) except for the contributions of bound states.

Fundamental representation
Let us first consider the fundamental representation. We find that the worldsheet instanton corrections are given by In Fig. 3 we plot the quantity is the perturbative part given by the Airy function (3.6) and µ * is the saddle point value of the chemical potential in (2.16), and W inst (N, N f ) is the instanton correction in the canonical picture up to worldsheet 3-instantons, obtained from the grand canonical picture (4.1) using (3.4) and (3.5). If we have subtracted worldsheet instantons correctly in (4.2), δ should decay exponentially as N becomes large. Indeed, in Fig. 3 we find that δ decays exponentially for N f = 7, 8, 11, as expected. We have also checked a similar behavior of δ for other values of N f . This confirms the correctness of the instanton corrections in (4.1). In (4.1), we observe that the worldsheet 1-instanton and 2-instanton have no poles in the convergence region N f > 2 (2.24). This suggests that there is no "pure" membrane instanton corrections as in the case of 1/2 BPS Wilson loops in ABJM theory [8], since there is no need for the membrane instantons to appear to cancel the poles. On the other hand, the worldsheet 3-instanton has a pole at N f = 4. We conjecture that this pole is canceled by the bound state of order e −2µ−4µ/N f . It would be interesting to determine the coefficient of this bound state contribution as a function of N f . Also, we observe that the grand canonical VEV (4.1) has two pieces which scale differently in the large µ limit: the constant term −N f /4 and the remaining part whose leading term gives rise to the perturbative part (3.2). One can translate this decomposition into the canonical picture where the second term behaves in the large N limit as In section 7, we will show that this decomposition (4.3) is consistent with the genus-zero result in [14].

Winding number m = 2, 3
Next we consider the worldsheet instanton corrections to the winding Wilson loops W m (µ, N f ) for winding numbers m = 2, 3.
Winding number m = 2 For m = 2 we find We observe that the last line of (4.5) is related to the derivative of the modified grand potential −∂ µ J(µ, N f ) + 2B, (4.6) where B is the coefficient in the perturbative part (2.11). As in the case of fundamental representation in the previous subsection, W 2 (µ, N f ) consists of two parts with different scaling behavior in the large µ limit, which implies that the Wilson loop VEV W 2 (N, N f ) in the canonical picture can be decomposed as where the second term in (4.7) behaves in the large N limit as Winding number m = 3 For m = 3 we find (4.9) In this case, we observe that W 3 (µ, N f ) consists of three parts with different scalings in the large µ limit. Noticing that the last term of (4.9) is proportional to W in (4.1), we conjecture that the canonical VEV W 3 (N, N f ) for winding number m = 3 is decomposed as where the last term scales as We have confirmed our result of instanton corrections (4.5) and (4.9) for m = 2, 3 by performing a similar test as in Fig. 3. We have checked a correct exponential decay of the quantity δ for various N f 's in the convergence region (2.24).

WKB expansion
In this section, we will consider the WKB expansion (small expansion) of spectral trace Tr(ρ s e x ) and try to reproduce the perturbative part of fundamental representation. Our starting point is the Mellin-Barnes representation of the grand canonical VEV [6,28] where c is a positive constant in the region 2/N f < c < 1. By picking up the poles at s = ∈ N we recover the small κ expansion (2.28). On the other hand, closing the contour in the direction Re(s) < c we can find the large µ expansion of W (µ, N f ).
In the quantum mechanical description of density matrix (2.4), the Planck constant is fixed to = 2π (2.5). However, one can formally introduce the parameter in the canonical commutation relation [ x, p] = i and perform the WKB expansion of the spectral trace Tr(ρ s e x ). Finally we set = 2π at the end of computation. This procedure was successfully applied to several examples [29,30].
At the zero-th order of WKB expansion, the spectral trace is given by the classical phase space integral and the higher order corrections can be systematically computed by using the Wigner transformation of operatorρ s e x [30][31][32][33]. In this way, we find the WKB expansion of spectral trace as where D(s) is a formal power series in We find that the n th order term D n (s) has the following structure: where p n (s) is a (3n − 1) th order polynomial of s. The first three terms are given by The coefficient c n of the highest order term N 2n f s 3n−1 in p n (s) (5.6) is found to be where B 2n (1/2) denotes the Bernoulli polynomial B 2n (z) evaluated at z = 1/2. We have computed D n (s) up to n max = 10. As mentioned above, the large µ expansion of grand canonical VEV can be found by closing the contour of (5.1) in the direction Re(s) < c. The perturbative part comes from the pole at s = 2/N f of Z 0 (s) in (5.2) In order to reproduce the coefficient c 1 (N f ) of the perturbative part of fundamental representation in (3.7), we need to show that Although we do not have an analytic proof of this relation, we can check this numerically by using the Padé approximation 2n D n (s) ≈ 1 + a 1 2 + · · · + a n max/2 nmax 1 + b 1 2 + · · · + b n max/2 nmax , (5.10) and set = 2π at the end. As we can see from Fig. 4a, the Padé approximation exhibits a nice agreement with the expected behavior in (5.9). One can also repeat the same analysis for the constant part −N f /4 in (4.1), which comes from the pole at s = 0 in (5.1) After setting = 2π, we find D(0) = 0 as expected. In order to reproduce the constant −N f /4, we need Again, we can check this relation numerically using the Padé approximation for D (0). We find a nice agreement with the right hand side of (5.13) (see Fig. 4b).
We expect that the worldsheet -instanton comes from the pole at s = (2 − 4 )/N f . In principle, we can find the instanton coefficients from the WKB analysis. However, it is difficult in practice since both the classical part Z 0 (s) and the quantum corrections D(s) should have poles at s = (2 − 4 )/N f in order to reproduce the Fermi gas result (4.1) 5 . This is different from the situation for the perturbative part and the constant term considered above, where the poles only come from the classical part. By the same reason, it is difficult to study the winding Wilson loop W m (µ, N f ) with m ≥ 2 from the WKB analysis.

Planar solution of N f matrix model
In this section, we will summarize the planar solution of the N f matrix model. This section is mostly a review of the result in [14], but we find that the the planar resolvent in [14] can be vastly simplified. Using this simplified expression of resolvent, we will directly show that the resolvent satisfies the planar loop equation without referring to the relation to the O(n) matrix model [34][35][36].
Let us consider the planar resolvent of N f matrix model in the 't Hooft limit (1.1) along the cut z ∈ [a, b]. Here V (z) is given by Also, from the definition (6.1) and the symmetry of the eigenvalue distribution we find Mapping to the torus C/(Z + τ Z) To find the resolvent, it is convenient to map z to the variable u by a Jacobi elliptic function z = a sn(u, k), (6.6) where the elliptic modulus k is given by In what follows we will suppress the dependence of k. We also need the derivative of z ∂ u z = a cn(u)dn(u) = a (1 − a 2 z 2 )(1 − a −2 z 2 ). (6.8) Furthermore, the variable u is related to the flat coordinate v on the torus C/(Z + τ Z) by v = u 2K , (6.9) where K = K(k 2 ) is the elliptic integral of the first kind and the complex structure of the torus is given by Then the coordinate z and the end-point of the cut a is written in terms of the Jacobi theta functions (6.11) We will use the variables z, u and v interchangeably. One can show that z(v) satisfies 12) and the various points on the z-plane are mapped to the points on the v-plane as Conditions obeyed by the resolvent ω(v) Here we will write down the conditions for the resolvent ω(v) in terms of the variable v. The loop equation becomes We also require that there are no cuts along In terms of the variable v these conditions become (6.14) Using the above conditions, ω(v) can be extended to a function on the torus obeying the following functional relations (6.17)

Planar resolvent
The resolvent ω(v) was obtained in [14] by taking a limit of the solution of O(n) matrix model [34][35][36]. After massaging the expression in [14], we find a very simple formula for the resolvent The function A(z) in (6.18) is given by

and A(u) denotes the function introduced in [14]
A(u) = πu 2K + K Z(u) + K ∂ u log z, (6.20) with Z(u) being the Jacobi zeta function In terms of the variable v, we find that A(v) in (6.18) is written as One can remove the linear term of v by performing the modular S-transformation As discussed in [37], we could use the S-dual variables (v , τ ) from the beginning. However, we will not do so and we will continue to use the original variables (v, τ ).
In the following, we will show that ω(z) in (6.18) indeed satisfies the necessary functional relations.
Symmetries of ω(v) First, let us consider the relation (6.5). In terms of v, (6.5) is written as This is satisfied since z(v + τ /2) = 1/z(v) and which follows from the identity of Jacobi theta functions Then, the τ /2-shift relation (6.25) implies that ω(v) is periodic (6.17) with period τ . One can also show the relation (6.16) by using (6.28) Next let us consider the normalization condition (6.4). So far λ in (6.18) is just a formal parameter in the ansatz of solution. The relation between λ and the 't Hooft coupling λ is fixed by the condition (6.4). Using A(0) = π/2, we find This shift by 1/8 is consistent with the Fermi gas result (2.11) at the leading order in the 't Hooft limit.
Absence of poles at z = ±1 We also want to show that the resolvent is regular at z = ±1. This requires One can show that (6.31) is satisfied by using The second equality in (6.32) is a consequence of the formula in the Appendix A of [14].
Absence of poles at z = ±i The regularity of resolvent at z = ±i determines λ as a function of a. One can show that A(u) in (6.19) is regular at z = ±i, and hence near z = ±i the resolvent behaves as From this behavior, the condition for the absence of pole at z = ±i is found to be where t is given by Using (6.8) we find that t is written as This reproduces the result in [14].
In [14], this was shown implicitly by taking a limit of the resolvent of O(n) matrix model. Here we will show the loop equation (6.15) directly.
To do this, it is convenient to introduce the operator T ± shifting v by ±1 Then the planar loop equation is written as (6.38) where L is given by 7 On the right hand side of (6.38), we have also used T · z = −z and V (v + 1) = V (−z).
For a rational function f (z) of z, the action of L reads Then it is natural to decompose the resolvent into the eigenfunctions G±(v) with eigenvalues T = −e ±iν . However, for n = 2 L has a double root at T = 1 which is a somewhat degenerate case. This is similar to solving the linear differential equation (d/dx − 1) 2 y = 0: it is well known that one of the solution y = xe x is not an eigenfunction of d dx . One can think of the function A(z) as an analogue of this solution y = xe x .
However, we should be careful about A(z) which transforms inhomogeneously under T ± T ± A(z) = A(z) ± π, (6.42) due to the linear term πv in (6.22). Then we find Now we are ready to prove (6.15). It is natural to decompose ω(z) in (6.18) into three parts where , One can easily show that the action of L on ω 1,2,3 is given by L ω 3 (z) = 0. (6.46) where we have used (6.43) in the first line. Adding these three equations (6.46), finally we arrive at the desired loop equation (6.38).

Planar free energy
In [14], the second derivative of the planar free energy is found explicitly as From (6.34), (6.36) and (6.47), we find that the derivative of λ and ∂F 0 ∂λ with respect to t have a simple form Note that the role of A-period and B-period is opposite from the standard definition. Using this relation (6.48), one can find the planar free energy F 0 as a function of 't Hooft coupling λ. In particular, in the small λ or large λ regime F 0 can be explicitly found as a power series. Let us first consider the large λ behavior of F 0 . In the large λ limit the size of the cut [−T, T ] in the original variable x becomes large, which implies a = e −T → 0. More precisely, we find that the large t expansion of a is given by Note that t and the shifted 't Hooft coupling λ = λ + 1/8 are related by (6.34) and the exponential correction e −t in (6.49) is identified with the worldsheet instanton factor (2.19) Then, integrating the relation (6.48) the planar free energy becomes where a 0 is a constant coming from the constant term A in the perturbative part (2.10) of grand potential [15] a 0 = lim On the other hand, in the small λ limit the size of cut [−T, T ] in the x-variable becomes small, which implies a → 1. From (6.34) and (6.36), we find that the small λ expansion of log a is given by and the free energy becomes . (6.54)

't Hooft limit of Wilson loops
In this section, we consider the 't Hooft expansion of normalized Wilson loop VEV Note that our normalization of VEV is different from [14] (see footnote 2.2 for our definition). As we will see below, we find a perfect agreement between the matrix model result and the Fermi gas result for the genus-zero part W (g=0) m of Wilson loop in (7.1).

Results of matrix model
We can read off the genus-zero VEV of Wilson loops from the small z expansion of the resolvent (6.18) 8 To write down the small z expansion of ω(z) in (6.18), let us first consider the small u expansion of A(z) in (6.19) where t is defined in (6.36) and the factor I is given by Here E = E(k 2 ) denotes the elliptic integral of the second kind. By inverting the relation z = a sn(u) in (6.6), we can write down the small z expansion of u Combining (7.4) and (7.6), we find that A(z) is expanded as We notice that the coefficients in this expansion (7.7) are some linear combinations of I and t. Plugging this expansion (7.7) into (6.18) and read off the coefficient of z m in (7.3), we can find the planar VEV of winding Wilson loops W (g=0) m up to arbitrary winding number m, in principle. For instance, the planar VEV of Wilson loop in the fundamental representation is given by .
(7.8) 8 Using the symmetry ω(z −1 ) = −ω(z), one can read off the Wilson loop VEVs from the large z expansion of resolvent as well This agrees with the result of [14]. For the higher winding numbers we find 9π 2 a 2 − a 2 + 1 7a 4 + 12a 2 + 7 tI 45π 2 a 4 + 8a 8 + 30a 6 + 97a 4 + 30a 2 + 8 I 2 180π 2 a 6 . (7.9) We observe that the planar VEV of winding Wilson loop is a linear polynomial of I for odd m, and quadratic in I for even m. This structure originates from the linear dependence of A(z) on I in (7.7).
Small λ expansion From the small λ expansion of log a in (6.53), one can easily find the small λ expansion of Wilson loop VEVs in (7.9). For general winding number m, we find that the Wilson loop VEVs are expanded as (7.10) For instance, the small λ expansion of the fundamental representation is given by (7.11) which reproduces the result in [14] 9 . We have checked that (7.10) correctly reproduces the expansion of W (g=0) m for m = 1, · · · , 6 in (7.9). We conjecture that (7.10) holds for any winding number m. It would be interesting to reproduce this expansion (7.10) from the perturbative calculation of matrix model along the lines of [14].
Large t expansion One can also study the large 't Hooft coupling, or large t behavior of Wilson loop VEVs (7.9). This large t regime is directly related to the Fermi gas result, which we will consider in the next subsection.
The large t behavior of W (g=0) m can be found from the small a expansion of I in (7.5) together with the large t expansion of a in (6.49). For the first three winding numbers m = 1, 2, 3, the large t expansion of W (g=0) m is given by ) . (7.13)

Comparison with Fermi gas
Let us compare the matrix model results (7.13) with the Fermi gas results in section 4. In the grand canonical picture, the 't Hooft limit is given by As discussed in [14], at the level of genus-zero the Wilson loop VEV in the canonical picture can be obtained by plugging the saddle point value µ * of the chemical potential into the grand canonical VEV W m (µ * , N f ) of Wilson loop. However, to study the instanton corrections we have to include the exponentially small corrections to the saddle point value µ * of the chemical potential, beyond the perturbative expression in (2.16). This is achieved by identifying the saddle point value µ * with the derivative of the planar free energy F 0 [14] µ ). (7.15) Plugging this expansion of µ * (7.15) into the grand canonical VEV W m (µ = µ * , N f ) in section 4 (eqs.(4.1), (4.5), (4.9) for m = 1, 2, 3, respectively), we have confirmed that the Fermi gas results perfectly match the matrix model results (7.13) in the planar limit. Let us take a closer look at the correspondence between the Fermi gas results and the matrix model results. The perturbative part (3.2) in the Fermi gas picture corresponds to the term I γ /a m in the matrix model result (7.9), where γ = 1, 2 for odd m and even m, respectively. We find that the coefficient of this term in W One can see that this matrix model result (7.17) is correctly reproduced from the perturbative part in the Fermi gas picture (3.2) in the 't Hooft limit (7.14) lim where µ should be identified with the saddle point value µ * in (7.15). Also, we observe that the matrix model result W (g=0) m (7.13) contains several pieces with different scalings in the large t limit, which naturally corresponds to the similar decomposition in the Fermi gas picture. For instance, the constant −1/4 in the fundamental representation (7.8) corresponds to the first term in the decomposition (4.3) observed in the Fermi gas picture. Similarly, the first two terms in the m = 2 VEV in (7.13) corresponds to the genuszero part of the first term in (4.7)

Conclusions
In this paper, we have studied the Wilson loops in the N f matrix model from the Fermi gas approach. We have determined the first few worldsheet instanton corrections to the winding Wilson loops for the winding number m = 1, 2, 3, and found that our Fermi gas result is consistent with the planar limit of matrix model result. We find that the Wilson loop VEVs can be decomposed into several pieces with different scaling behavior in the large N limit. Also, we conjecture that the grand canonical VEVs of winding Wilson loops do not receive "pure" membrane instanton corrections except for the bound state contributions. This is reminiscent of the instanton corrections to the 1/2 BPS Wilson loops in the ABJM theory [8,10]. There are many interesting open problems. To study the partition functions and Wilson loops in the N f matrix model further, it is very important to understand the structure of bound states. In the case of ABJM theory, the effect of bound states can be removed by introducing the "effective" chemical potential µ eff [38], which in turn is related to the quantum period of the quantized mirror curve of local P 1 × P 1 [39,40]. It would be interesting to see if one can define a similar "effective" chemical potential in the N f matrix model as well.