Matrix Model of Chern-Simons Matter Theories Beyond The Spherical Limit

A class of matrix models which arises as partition function in U(N) Chern-Simons matter theories on three sphere is investigated. Employing the standard technique of the 1/N expansion we solve the system beyond the planar limit. In particular we study a case where the matrix model potential has 1/N correction and give a general solution thereof up to the order of 1/N^2. We confirm that the general solution correctly reproduces the past exact result of the free energy up to the order in the case of pure Chern-Simons theory. We also apply to the matrix model of N=2 Chern-Simons theory with arbitrary numbers of fundamental chiral multiplets and anti-fundamental ones, which does not admit the Fermi gas analysis in general.


Introduction
Recent progress in supersymmetric Chern-Simons matter theories has been made on the basis of the exact result by means of the supersymmetric localization. This technique allows one to compute the partition function of supersymmetric theories exactly by the steepest descent method, which reduces the path integral calculation to a certain matrix model one [1]. This calculation was done generically on S 3 [2,3,4] (see also [5,6]) and on S 2 × S 1 called superconformal index [7,8] (see also [9,10,11]). These exact results were in precise agreement with the prediction from AdS 4 /CFT 3 duality [12,13] (see also [14,15]). See [16] for review and references therein.
On the other hand, progress in non-supersymmetric Chern-Simons matter theories has also been made not relying on the localization technique but on the 1/N expansion technique by restricting the class of matter fields to vector fields. It was conjectured that such system is exactly soluble in the 't Hooft large N limit [17,18]. The thermal partition function in Chern-Simons vector models on S 2 × S 1 was determined exactly in the leading of the 1/N expansion near the critical high temperature [19] (see also [20,21]), which was used to show the three dimensional bosonization duality [22,23,24].
In contrast the exact large N analysis of three sphere partition function for any nonsupersymmetric Chern-Simons matter theory is not performed due to technical difficulty. So far analysis of the three sphere partition function was done perturbatively near the weak coupling limit of the Chern-Simons coupling constant [25] to confirm that the system obeys the F-theorem [3]. Perturbative analysis is, however, not enough to provide evidence for duality and exact large N analysis is strongly awaited.
In this situation we change gears to study a class of matrix models which is to be obtained as the three sphere partition function of Chern-Simons matter theories in order to capture a generic feature of such a class of matrix models toward a big goal to show bosonization duality on the three sphere. We are also interested in analyzing such a class of matrix models beyond the planar limit because the bosonization duality is expected to hold at the subleading in the 1/N expansion [21]. (See [26,27,28] for recent arguments. ) To illustrate what kind of matrix models is to be studied, let us consider the partition function of generic U(N) k Chern-Simons matter theories on the three sphere with unit radius.
where Φ denotes all the matter fields collectively and S[Φ, A] is the action for the matter fields. This may be computed perturbatively as follows. (See [29] for details.) We first expand the gauge field by the vector spherical harmonics We take the Lorenz gauge ∇ µ A µ = 0, which kills the modes a s,s l,r with s > 0. The residual gauge can kill the mode a 0,0 0,0 except its Cartan part, which we denote by σ. We expand the matter fields in a similar manner. Taking into account the Faddeev-Popov determinant we integrate out all massive modes such as a s,s+1 l,r , a s+1,s l,r and all the modes coming from the matter fields, which are massive on the three sphere. Then the partition function reduces to the finite-dimensional integration of the effective action over σ.
where the ellipsis is some function of σ s generated by integrating out the massive modes. When the matter fields vanish, for example, the effective action can be exactly computed by the supersymmetric localization or the cohomological localization [30] by adding the auxiliary fields so as to complete the gauge field in the N = 2 vector multiplet. The result is [31] Z ∼ . (1.4) Then by denoting the correction of the matter fields to the effective action by V [σ], the partition function is given by a form such that (1.5) In this note we analyze this class of matrix models by restricting the form of potential to be consisting of single trace operators: V [σ] = N N s=1 W σ (σ s ). 1 The goal of this paper is to solve this class of matrix models incorporating the standard technique of 1/N expansion developed in the study of ordinary hermitian matrix models [32] beyond the spherical limit [33,34]. See [35,36] for review and further references.
The rest of this paper is written as follows. In §2 we perform some preliminary analysis of the matrix models. In §3 we derive the loop equation for this class of matrix models. In §4 we solve the loop equation by using the 1/N expansion. We give a general solution for the planar limit ( §4.1) and for the genus one ( §4.3). Then we apply this solution to a few examples in §5. We first apply to pure Chern-Simons theory and compare with the known exact result to test the validity of the presented framework ( §5.1). We then apply to N = 2 Chern-Simons theory with arbitrary numbers of fundamental and anti-fundamental chiral multiplets ( §5.2). §6 is devoted to discussion and future direction. In §A we give a brief review on the partition function of pure Chern-Simons theory on three sphere in order for this paper to be self-contained.

Matrix model of Chern-Simons matter theories
Throughout this paper we investigate a class of matrix models such that where N is a normalized constant and W σ (σ) is a non-singular function of σ, which can be determined at least perturbatively by integrating out the massive modes for the original Chern-Simons matter theory. When the original theory has N = 2 supersymmetry, the matrix model potential can be determined exactly by using the localization method [2]. For example, for N = 2 U(N) k Chern-Simons theory with N F fundamental chiral multiplets with the canonical R-charge, the partition function is of a form such that [2,3] , and ζ is the FI parameter, which is equivalent to the real mass parameter. Adding the same number of anti-fundamental chiral multiplets to this system, the partition function becomes See [37,38] for detailed analysis of this type of matrix models. For explicit computation we write the form of potential in (2.1) as where λ, t p are parameters. In this parametrization the partition function of pure Chern-Simons theory on S 3 , which we briefly review in §A, is given by See (A.4) for the normalization. The first example (2.2) is formally given by where p ≥ 1. 2 The second example (2.3) is The matrix model (2.1) can be recast in the same form as hermitian matrix models with positive eigenvalues. By changing integration variables so that φ s = e σs , the partition function becomes where R + represents the positive real axis and where δW (φ s ) is analytic on the positive real axis. Note that the matrix model potential has the logarithmic cut on the negative real axis. As a result the partition function can be written by using a positive definite hermitian matrix Φ as . This suggests that this class of matrix models can be analyzed by using the standard technique employed in ordinary hermitian matrix models. In what follows we show that the free energy and correlators of some sector can be determined in order in the 1/N expansion.

Loop equation
There is a well-known method to determine the free energy in matrix models by using the so-called resolvent, which is in the current situation defined by the the vacuum expectation value of the generating function of regular single trace operators: We remark that the resolvent is well-defined around the infinity and formally behaves as ω(z) ∼ 1 z in the vicinity of the infinity though this asymptotic behavior is not guaranteed due to the fact that the potential (2.9) has the logarithmic cut, which ends in the infinity. This suggests that the behavior of the resolvent around the infinity generally gets logarithmic correction. Still we can expect that there exists a limit approaching to the infinity such that the resolvent behaves as ω(z) ∼ 1 z on a certain patch. This request will be important to determine the resolvent later.
Once the resolvent is determined, the coupling dependence on {t p } of the free energy given by F = − log Z is determined by Consider a one-to-one transformation on R + denoted by φ s → φ ′ s . Then we obtain the identity such that Suppose the (infinitesimal) transformation φ ′ s = φ s + aδφ s . Expanding the right-hand side in terms of a, we find that the zero-th order term cancels the left-hand side and the equation at the linear order gives the Schwinger-Dyson equation The first two terms are computed as The third term is where we define the density function Hereafter we assume z is outside the support of the density function in order to exclude the case where (3.5) is divergent. The density function satisfies R + dxρ(x) = 1. The density function can be computed by evaluating the discontinuity of the resolvent across the real axis. Indeed, by using the formula where x is a real number, ǫ is an infinitely-small positive number and P denotes the principal value, the discontinuity of the resolvent between x ± iǫ is computed as By using this relation, (3.5) can be rewritten as where C R + denotes a circle encircling R + counterclockwise. The 4th term vanishes because Collecting these we obtain the loop equation Note that this is of the same form as that of ordinary matrix models except the integration region.
Once the resolvent is determined by the loop equation, so are the density function by (3.8) and the coupling dependence on {t p } of the free energy by (3.2), and similarly for other coupling dependence. For example, when δW (Φ) = p≥0 u p Φ p log Φ, the other coupling dependence on {u p } of the free energy is determined by ∂ ∂up and υ(z) is the vacuum expectation value of the generating function of singlet operators of the form Tr(Φ p log Φ): which is computed by using the resolvent ω(z) as (3.14) Similarly the λ dependence of the free energy is determined as Acting d dTz and d dUz on the free energy gives correlators of singlet operators such that with n > 1, where the subscript conn means the connected part of the correlator. It is also possible to compute correlators including some number of the operator Tr(log Φ) 2 by differentiating the free energy several times with respect to λ. We remark that it is not guaranteed and has to be confirmed that correlators computed this way agree with those computed from the original theory by using the path integral. 3 In other words, the potential of the matrix model generally depends on operators inserted in the path integral. This will be easily seen by considering a partition function of some supersymmetric theory computed by using the localization method. Since correlators of non-supersymmetric operators cannot be computed by the exact method, the potential of the matrix model cannot be reused to compute correlators of non-supersymmetric operators. In this context the matrix model potential is available when operators inserted in the path integral or parameters deforming the original theory maintain supersymmetry.

Solution of the loop equation
In this section we solve the loop equation (3.11) in the 1/N expansion. Analysis will depend on the large N behavior of the potential. Generically it can be written as For explicit calculation maintaining certain extent of generality we study the case where the potential is given by (2.9) with λ, t p of order one: In the examples of supersymmetric Chern-Simons theory in §2, this case corresponds to the 't Hooft limit with the number of flavors N F of order N. Accordingly the consistent 1/N expansion of the resolvent will be such that Plugging this into the loop equation and expanding with respect to 1/N the loop equation is decomposed as follows. For g = 0, which we call the genus zero and half loop equation respectively for convenience, and for g ≥ 1, where we defineK From these equations ωḡ(z) can be determined in order fromḡ = 0. Once the resolvent is determined at the orderḡ in the 1/N expansion, so is the density function from (3.8) as The coupling dependence of the free energy on t p is determined from (3.2) as

Planar solution
Let us solve the planar loop equation (4.4). First we show that the planar loop equation contains the saddle point equation of the starting matrix model in the large N limit. For this purpose we compute discontinuity of both sides in (4.4) between x − iǫ and x + iǫ.
The discontinuity of the left-hand side is where we used (4.8). That of the right-hand side is Therefore we obtain with x in the support of the density function in the leading of the 1/N expansion. This is the same as the saddle point equation derived from the starting matrix model (2.8) in the large N limit. Suppose that the support of the density function consists of s distinct connected in- Taking account of the fact that the loop planar equation (4.4) is quadratic, we make each interval correspond to a square root cut of the solution. Under this ansatz we solve (4.10). Let us consider a trial function H(z) which sees the deviation of the resolvent from the s-cut square root function As mentioned in the previous section, we solve the loop equation so that the resolvent behaves as ω 0 (z)∼ 1 z in a limit approaching the infinity. This suggests that the trial function behaves as H(z)∼ 1 z s+1 up to the signature and thus is analytic around the infinity.
Therefore using the Cauchy theorem we find 4 where z is a complex number outside supp(ρ 0 ) and C ∞ is an infinitely large circle. Assuming further that the trial function is analytic except the support of the leading density function, we can compute the left-hand side by deforming the contour to the non-analytic region: where we used (4.10) in advance and C supp(ρ 0 ) denotes a circle encircling the intervals supp(ρ 0 ) counterclockwise. Therefore the trial function is determined as so is the planar resolvent: .
Then the planar density function is computed from (4.8) as The endpoints of the cuts are determined in the following way. Assume the solution obtained above behaves asymptotically as ω 0 (z) = 1 z +· · · approaching to the infinity. This is satisfied if and only if where the signature is chosen suitable. These give s + 1 constraints for 2s endpoints of the cuts, which is not sufficient unless s = 1. For s ≥ 2 case, the residual conditions are provided by stability against the tunneling of eigenvalues between different cuts [39]. We 4 One may more generally conclude that, for example, C∞ dw 2πi are polynomials of z. The resolvent obtained from this form in the same way as described below at first looks different from (4.15) but reduce to the same form by using the boundary condition, which is the same as (4.17). We give a comment on this point below.
demonstrate the residual condition following [40]. First we write the total matrix model potential in terms of the density function in the large N limit.
where µ is Lagrange multiplier and ̺ 0 is the dynamical planar density function determined by the saddle point equation Differentiating this with respect to x leads to (4.10), which we solved as (4.15). This suggests that integrating (4.10) with respect to x does not get back to (4.19) because the density function is not analytic on the edges of the cuts so integration of (4.10) takes different values on every intervals in general. Requiring those values to be the same (as µ) gives a nontrivial condition. 5 To write down the condition we define a function µ on R + by where ρ 0 (y) is defined by the analytic continuation of ρ 0 (y) from an interval to the whole positive real axis. Then the function ρ 0 (x) takes pure imaginary values outside supp(ρ 0 ). Differentiating with respect to x gives µ(x) ′ = Re[−2πi ρ 0 (x)], which suggests that the derivative of µ(x) vanishes on each cut and thus µ(x) is constant on each cut, as expected.
There is a comment on the solution (4.15). By using the condition (4.17) the solution can be rewritten in a different form such as where k = 0, 1, · · · , s. On the other hand, this form of solution can be obtained directly by starting with a different equation from (4.12) as mentioned in the footnote 4. Then the cuts can be determined not only by the asymptotic condition ω 0 (z) = 1 z + · · · with z ∼ ∞, but also the fact that ω 0 (z) is non-singular around z ∼ 0. These two conditions give the condition (4.17).
The free energy in the leading of the 1/N expansion is easily determined as with (4.19). From this expression, we can reproduce the derivatives of the free energy with respect to the coupling constants obtained in the previous section. For example, acting d dTz on (4.23) yields where we used (4.19), R + dxρ 0 (x) = 1 and (4.2). This is nothing but the equation (4.9).

Hole correction
The hole correction of the resolvent is determined by the genus half loop equation (4.5). Let us derive the saddle point equation at this order from (4.5). For this purpose let us rewrite (4.5) as As the planar case we compute discontinuity of the left-hand side between x−iǫ and x+ iǫ.
The discontinuity of the first term is computed as That of the second term is Therefore we obtain with x ∈ supp(ρ 0 ). Combining this with the planar saddle point equation (4.4) we obtain . This is the same form as the planar saddle point equation (4.4) replacing W 0 (x) with W (x) and the previous argument to solve this equation holds without any modification. Therefore a solution of the usual saddle point equation is correct up to the order of hole correction! This is a nice simplification, while the caveat is the region which (4.27) holds. That is, the region where we need to solve (4.27) is on supp(ρ 0 ). However when we solve (4.27) as done in §4 the cut appears as the support of the density function including the hole correction, supp(ρ 0, 1   2 ). This small discrepancy may imply that the loop equation can be solved by assuming that the support of the planar density function matches the one including the hole correction. This assumption may be important to separate out the genus half one from ω 0, 1 2 (z). We expect that the discussion above will hold in more general matrix models such as two-matrix models. We leave the proof of this conjecture to future work.

Genus one correction
Let us determine the genus one correction of the resolvent. For simplicity we first study the case where W 1 = 0, so ω 1 2 = 0. In this case the genes one loop equation, (4.6) with g = 1, reduces toK This can be solved in the same manner as in the ordinary hermitian matrix model [33,34]. Let us first compute the right-hand side.
Then the second term is computed as ).
The third term is (4.29) Therefore we obtain In order to compute da i dTz , we act d dTz on the constraint equations of the edges of the cuts, (4.17), (4.21).
where we used 2πi ρ 0 (x) = W ′ 0 (x) − 2ω 0 (x). These are computed as The solution can be written as where α i,l are determined by plugging this back and setting the coefficients of polynomials with respect to z to zero. The determining equations of α i,l are Substituting (4.35) into (4.30) we find This can be rewritten as the image of the linear operatorK in such a way that is constructed inductively as follows. Start with the identity such that h(w) on both sides and computing the right-hand side results in (4.43) The case with n = 1 implies χ (z−a i ) n is valid for n ≤ n − 1, we can rewrite the right-hand side as i (z). Therefore we find thenKχ (n) i (z) = 1 (z−a i ) n is valid for n = n. By using this function we finally solve the genus one loop equation as 46) up to terms in the kernel of the operatorK such as z m h(z) with m = 0, · · · , s. Remark that ω 1 (z) behaves at most 1 z s+1 and thus the leading asymptotic behavior of the total resolvent ω(z) is unchanged, so is the cut. From (4.46) the coupling dependence of the genus one free energy can be determined. For example the dependence on t p is determined by (4.9) and that on λ is by (3.15). 7 Next we consider the case where the matrix model potential contains 1/N correction: W 1 (w) = 0. In this case, as derived in (4.6), the genus one loop equation is corrected by the genus half resolvent so thatKω 1 . This equation is more involved and there may be simplification in a way that the planar resolvent and the genus half one can be determined at the same time as shown in §4.2. To see this let consider the deviation of the resolvent from the solution ω(z) = ω 0, 1 2 (z) + δω(z) and substitute this into the original loop equation (3.11). We obtain . Therefore at the leading order of the 1/N expansion this reduces tô where we defineK byK This equation is of the same form as the one without the hole correction by replacing ω 0 , W 0 with ω 0, 1 2 , W , respectively. Since the above argument to solve this equation holds as it is by performing the replacement, a solution of (4.48) is given by the same form as (4.46) where a i and M (n) i are replaced with the ones including the hole correction. We emphasize that this simplification happens only at the genus one order and at higher order one may need to solve (4.6) in general.

Applications
In this section we apply the presented formulation developed in the previous section to a few examples. Firstly we apply to the three sphere partition function in U(N) k pure Chern-Simons theory in order to test the presented framework by comparing the exact result known for pure Chern-Simons theory as reviewed in §A. Secondly we apply to N = 2 U(N) k Chern-Simons theory with n F fundamental chiral multiplets andn F antifundamental ones. This system does not admit the Fermi gas analysis in general and there may be no systematic way to study the system beyond the spherical limit except our formulation at present.

Pure Chern-Simons theory
The matrix model potential for pure Chern-Simons theory is given by (2.8) with (2.5). For simplicity we first study a case such that there is no hole correction. Then W ′ (w) = W ′ 0 (w) = log w w λ + 1 w . Let us first determine the planar resolvent. For λ > 0, the potential has only one stable minimum so we have only to consider a solution with one cut: supp(ρ 0 ) = [a − , a + ] with 0 < a − < a + . In order to simplify the integration we start with the solution of the form (4.22) with k = 1: . Inflating the contour we can compute the right-hand side as where we computed the first term as The edges of the cut a − , a + are determined by the asymptotic behavior around the infinity and the regularity around the origin. ω 0 (z) can approach to 1 z when z → −∞, which is achieved if and only if log These can be solved as a ± = (e λ 2 ± e λ − 1) 2 . Note that a − a + = 1. Then the planar resolvent can be simplified as The planar density function is computed as  Figure 1: In Fig.1(a), the blue and yellow curves depict the matrix model potential and the planar density function, respectively, when λ = 0.5. Eigenvalues tend to clump around the potential minimum. In Fig.1(b), the differentiation of the planar free energy with respect to λ is plotted. The blue curve depicts the result obtained by the resolvent method and the yellow one from the past exact result. They almost coincide.
This solution matches the one given in [36], where the solution is expressed in the original coordinates. We plot the density function as well as the potential in Fig.1(a). The planar free energy is given by (4.23) and its λ derivative is (3.15): We could not perform the integration in the right-hand side analytically, so instead we evaluated it by numerics. The numerical result is in good agreement with the past exact result (A.20) as seen in Fig.1(b). Next we study the genus one correction. Now we consider the one cut solution, so the genus one correction of the resolvent is given by where χ ± are defined by (4.45). M (1) ± are computed by (4.29). As done in the computation of the resolvent, the integration can be simplified by using (4.17). .
By inflating the contour to the infinity, the first term is computed as and the second term vanishes. Therefore M (1) . In the same way M ± are computed as M The genus one correction of the free energy is computed by using (3.15): This time we could perform the integral analytically: This result is in precise agreement with the past exact result of (A.21) with ζ = 0.
Next we consider a case where the matrix model potential has a hole correction by the , which is still integrable as shown in §A. As discussed in §4.2 we solve the usual saddle point equation ω(x−iǫ) + ω(x+ iǫ) = W ′ (x), which is correct up to the hole order. To emphasize the difference from the previous computation we denote the cut with prime added so that supp(ρ 0, 1 . Then a solution of this saddle point equation, ω 0, 1 2 (x), is given by . This can be computed in the same way as previously and we obtain 9 The edges of the cut are also determined in the same way. The result is a ′ ± = ca ± , where c = e − λ ζ N and a ± are the same as previously. By using this the resolvent up to the hole order is simplified as ω 0, 1 2 (z) = c −1 ω 0 (ź), whereź = c −1 z. As argued in §4.3 the genus one resolvent is given bỹ The usage of this terminology can be justified by adding some auxiliary fields into pure Chern-Simons theory so that the theory has N = 2 supersymmetry. 9 The planar resolvent and the hole one are determined from this bý which should be done before the edges of the cuts are determined.

whereχ (n)
i (z) is given by (4.45) with a i replaced into a ′ i . By using a ′ ± = ca ± the genus one resolvent including the FI term can be written asω 1 (z) = c −2 ω 1 (ź). Finally we compute the differentiation of the total free energy with respect to λ.
where the ellipsis represents the terms of order N −3 . Then the 2nd term is computed as The 3rd term is As a result we obtain This is in perfect agreement with the past exact result of (A.21).

N = 2 Chern-Simons theory with arbitrary numbers of fundamental and anti-fundamental chiral multiplets
As another example we consider the matrix model of N = 2 Chern-Simons theory with n f fundamental chiral multiplets andn f anti-fundamental ones with the canonical R-charge. Let us set n Without losing generality, we can assume that n f ≥n f . The matrix model potential of this system is given by combining (2.2) and (2.3) Its derivative is where we used ℓ ′ (z) = −πz cot(πz). Thus it is clear that the matrix model potential takes values in the complex number for a general number of chiral multiplets.
In order to determine the number of the cuts in the resolvent by identifying that of the potential minimum, we regard the matrix model potential as an analytic function with all the parameters. When n (−) f is pure imaginary and λ, ζ, n (+) f are real, the potential becomes real. We fix the number of the cuts of the resolvent in this situation.
We consider the large k, N limit holding its ratio and other parameters ζ, n (±) f fixed. In this limit the potential has only one stable minimum as in the case of pure Chern-Simons theory so we have only to consider a solution with one cut: supp We compute the resolvent up to the hole correction by (4.22) with k = 1: , (5.19) where h(z) = (z − a − )(z − a + ). Inflating the contour we can compute the right-hand side by picking up the pole as , (5.20) where the first term is the contribution of the logarithmic branch cut (−∞, 0], the second one is the one of the pole at w = −1 and the third one is at w = z. We compute integrations in a manner that where Then the resolvent becomes The edges of the cut a − , a + are determined by the asymptotic behavior around the infinity and the regularity around the origin. ω 0, 1 2 (z) can approach to 1 z when z → −∞, which is achieved if and only if is regular at the origin if and only if (5.26) From these equations the edges of the cut are determined order by order in 1/N.
As argued in §4.2 the planar resolvent should be determined so as to have the same cut as that of ω 0, 1 2 (z). Since the leading part of the potential in the large N limit is unchanged, the form of the planar resolvent is unchanged except the edges of the cut: (5.2) with a i → a i . The genus half resolvent is determined before the edges of the cut are expanded in the 1/N power series and given by ) .
(5.27) Then as argued in §4.3, the genus one resolvent is given bỹ where χ (n) i (z) are given by (4.45) with a ± replaced into a ± . The λ derivative of the free energy up to the genus one order is given by It is known that this system has the dual description known as Seiberg-like duality [41]. The dual theory is U(N ′ ) −k Chern-Simons theory with n f fundamental andn f antifundamental chiral multiplets with n fnf mesonic operators as well as some monopole operators with a suitable superpotential, where N ′ depends generally on N, k, n f ,n f , which is still in the class investigated in this paper. It would be interesting to test the duality from our general solution. We hope to come back to this problem in a future publication.

Discussion
In this paper we have performed a general analysis on a class of matrix models describing Chern-Simons matter theories on three sphere incorporating the standard technique of the 1/N expansion developed in the study of ordinary hermitian matrix models. We have derived the loop equation for all orders in the 1/N expansion and presented its explicit solution up to the genus one order when the potential has 1/N correction. We have applied the formulation to pure Chern-Simons theory and confirmed that the presented solution reproduces the exact result known in the past. We have also applied the framework to N = 2 Chern-Simons theory with arbitrary numbers of fundamental and anti-fundamental chiral multiplets and obtained a formal expression of the solution up to the genus one order in the 1/N expansion.
This paper mainly focused on construction of the framework to solve a class of matrix models. We are very much interested in applying the formula obtained in this paper to a duality pair of Chern-Simons matter systems and testing the bosonization duality holds at the next leading order in the 1/N expansion. In particular it would be interesting to develop the presented large N technique in a class of unitary matrix models which arises as partition function of Chern-Simons matter theories on S 2 × S 1 . For a class of Chern-Simons vector models, the effective matrix model potential was determined exactly in the leading of the large N limit [19] and the three dimensional bosonization was confirmed at the order. We hope that the formulation developed in this paper is useful for future study in this direction.
In this paper in order to study beyond the planar limit we adopted the iterative procedure given in [33,34]. Another iterative approach has been proposed by using the Feynman graph of the trivalent vertexes [42,43]. It would be interesting to reformulate the presented formula in this note in terms of the different approach.
Another interesting question is whether this class of matrix models has the equivalent description of some two dimensional CFT as the ordinary hermitian matrix models [44,45,46]. (See also [47].) Naively this answer seems to be no due to the fact that the degrees of freedom in 3d system are much bigger than those of 2d one in a generic situation. However we have a suspicion that the answer could be yes for a certain matrix model of this kind, intuitively because vector models coupling to Chern-Simons theory appear as effective field theory of anyonic system [48,49] and wave function describing a quantum Hall state known as Laughlin wave function [50] is given by a correlator of certain 2d (rational) CFT [51]. In fact it was shown that this answer becomes yes for a similar class of matrix models to the one studied in this paper [52], where the corresponding CFT is identified with q-deformed one. Exploring this question is left to future work.
There is a straight-forward generalization of the presented formulation to a different gauge group [53] or two matrices. This generalization to two matrices is important for the application to higher supersymmetric Chern-Simons matter theories such as the ABJM theory [54]. The 1/N correction of the free energy in the ABJM theory was computed in [55,56,57]. In this development a new technique called Fermi gas approach was invented [56]. This approach is powerful to study non-perturbative aspects of the ABJM theory from the S 3 partition function [58,59]. It is an important problem to test whether the traditional techniques of matrix models can reproduce the results obtained by new ones in recent development beyond the spherical limit.
We hope to come back to these issues in the near future.

Acknowledgments
The author would like to thank S. Sugimoto and T. Takayanagi for valuable discussions and comments on the draft. The author would also like to thank Y. Imamura for a helpful comment on the first version of this paper.
A Partition function of pure Chern-Simons theory on S 3 In this appendix we briefly overview the three sphere partition function in U(N) k pure Chern-Simons theory and a derivation of its large N expansion used in the main text. The partition function is defined formally by a path integral over the gauge field on S 3 such that 10 Z The classic paper [60] demonstrated explicitly in the case of SU(2) that this can be exactly determined as a function of the Chern-Simons level without performing the path integral by clarifying its relation to a modular transformation matrix of the characters in the corresponding Affine Lie algebra. Generalization to arbitrary gauge group is straightforward. Since modular transformation matrices were already determined in general Affine Lie algebras [61], the exact result of the partition function for U(N) k pure Chern-Simons theory was given by After this exact result was studied in terms of the 1/N and 1/k expansions [62,63], it was insightfully observed that the Chern-Simons partition function (A.2) exactly matches that of the topological string theory on a Calabi-Yau three-fold background by identifying the string coupling constant with the pure imaginary Chern-Simons level [64]. This lead to the conjecture of gauge/geometry duality [65] between Chern-Simons/Topological string theories [66,67].
Subsequently it was pointed out that the partition function (A.1) reduces to a matrix model such that [31] This matrix model was extensively studied in relation to the topological string theory [68,69]. With the help of the Weyl denominator formula the matrix integral was explicitly performed by gaussian integration in perfect agreement with (A.2) [2]. This matrix model was also evaluated exactly by the orthogonal (or characteristic) polynomial method in accordance with (A.2) [70]. The orthogonal polynomials associated with this matrix model were found to be Stiltjes-Wigert polynomials. Let us compute the matrix model (A.4) in a Fermi gas like approach [56] including the "FI term": The inside of the determinant is computed by gaussian integration as follows. The 1/N expansion was done as follows [64,66]. The expansion coefficients are determined as functions of λ = N k . Here B n is the n-th Bernoulli number defined by 11 x e x − 1 = ∞ n=0 B n n! x n . (A.14) The summation in the last term in (A.13) can be done by using a formula such that 12 Our definition of the Bernoulli number is different from the one adopted in several past literatures such as [63,64,66]. The difference is B where [x] is the integer part of x.
Plugging these back, we obtain the free energy as As a result, the coefficients in the 1/N expansion of the form (A.11) are determined as ζ(2m) m λ 2m (2m + 1)(2m + 2) , (A.17) A few comments are in order. The leading term in the 1/N expansion, F 0 , can be obtained directly from (A.10) by taking the large N limit [62].
where Li s (z) is the polylogarithm defined by Li s (z) = ∞ n=1 z n n s . The next leading term F 1 can be simplified by using (A.12) as follows.