3d $\mathcal{N}=4$ OPE Coefficients from Fermi Gas

The partition function of a 3d $\mathcal{N}=4$ gauge theory with rank $N$ can be computed using supersymmetric localization in terms of a matrix model, which often can be formulated as an ideal Fermi gas with a non-trivial one-particle Hamiltonian. We show how OPE coefficients of protected operators correspond in this formalism to averages of $n$-body operators in the Fermi gas, which can be computed to all orders in $1/N$ using the WKB expansion. We use this formalism to compute OPE coefficients in the $U(N)_k\times U(N)_{-k}$ ABJM theory as well as the $U(N)$ theory with one adjoint and $N_f$ fundamental hypermultiplets, both of which have weakly coupled M-theory duals in the large $N$ and finite $k$ or $N_f$ regimes. For ABJM we reproduce known results, while for the $N_f$ theory we compute the all orders in $1/N$ dependence at finite $N_f$ for the coefficient $c_T$ of the stress tensor two-point function.

Wilson loops in terms of (N − 1)-dimensional matrix model integrals [1]. Localization was then applied by [2] to 3d N = 2 Chern-Simons-matter theories, such as the N = 6 ABJM theory with gauge group U (N ) k × U (N ) −k and Chern-Simons level k [3], to compute matrix model observables as N -dimensional integrals. For low N , these integrals can be computed by hand, but this becomes infeasible for large N . The large N regime is of particular interest for conformal field theories (CFTs) with holographic duals, such as 4d N = 4 SYM dual to Type IIB string theory on AdS 5 × S 5 , and 3d ABJM dual to Type IIA string theory on AdS 4 ×CP 3 for large k and M-theory on AdS 4 ×S 7 /Z k at small k. When N is large, the bulk is described by weakly coupled supergravity, so if the localization integrals can be computed at large N , then they could be compared to bulk calculations. For N = 4 SYM, the matrix model is Gaussian, so observables in the matrix model can be computed to all orders in 1/N or even finite N using standard matrix model methods like orthogonal polynomials and topological recursion and then successfully matched to bulk calculations [4][5][6][7][8]. For ABJM the matrix model is much more complicated, but nevertheless topological recursion was successfully applied in [9] to compute matrix model observables in the 't Hooft limit of large N, k and fixed N/k. This limit could not be used to study the M-theory regime, however, which is only weakly coupled in the bulk for large N and finite k.
To study the M-theory limit of ABJM, [10] showed that the ABJM matrix model could be written in the form which describes a gas of non-interacting fermions with a certain 1-particle density matrix ρ and ∼ k. The partition function could then be computed in a small ∼ k WKB expansion using the Wigner approach to quantum statistical mechanics, and the all orders in 1/N result could be written in terms of an Airy function. At leading order in large N , the resulting sphere free energy scaled as N 3 2 as predicted from M-theory. Wilson loops could also be computed as thermal averages of 1-body operators in the Fermi gas [11], and the all orders in 1/N result also took the form of Airy functions. This Fermi gas method was soon generalized to a wide variety of 3d N = 4 theories (for a review with extensive references, see [12]). For instance, the N = 4 U (N ) gauge theory coupled to one hypermultiplet in the adjoint representation and N f hypermultiplets in the fundamental is holographically dual to M-theory on AdS 4 × S 7 /Z N f for finite N f [13,14]. In [15,16], the sphere partition function of this N f matrix model was written in the form of (1.1), and was then computed to all orders in 1/N in terms of Airy functions.
Applications of the Fermi gas method have so far been limited to the sphere partition function and the expectation value of Wilson loops. Recently, it has become clear that a much wider class of physical observables can be computed using localization. One such class are integrated correlators of conserved currents. Since localization does not require conformal symmetry, the sphere partition function can be computed for a theory deformed by real masses m, which couple to flavor multiplets, as well as on the squashed sphere with squashing parameter b [17][18][19], which couples to the stress tensor multiplet. Taking n derivatives of the free energy in terms of m (or b) and then setting m = 0 (or b = 1) then computes the correlator of n flavor currents (or stress tensors) integrated over the sphere.
For instance, consider the coefficient c T of the two-point function of canonically normalized stress-tensors: which is related by the conformal Ward identity to the OPE coefficient λ OOT for any scalar operator O as c T ∝ 1/λ 2 OOT [20]. We can compute c T by taking two derivatives of b as which expresses c T in terms of the expectation value of a certain 2-body operator in the matrix model. This method was used in [21][22][23] to compute c T in a wide variety of 3d N = 2 theories whose matrix models involved only a small number of integrals that could be performed by hand. For ABJM, the three real masses m i couple to the R-symmetry currents that are in the same multiplet as the stress tensor, so c T can be computed from The mass deformed partition function Z(m i ) was computed to all orders in 1/N using the Fermi gas method in [24], which was used in [25] to compute c T to all orders in 1/N . The integrated stress tensor multiplet 4-point function was similarly computed from [26] and used to constrain the correlator in the large N expansion. Another class of observables that can be computed using matrix model expectation values are OPE coefficients of half-BPS scalar operators in 3d N = 4 CFTs. As shown in [27,28], all 3d N = 4 CFTs contain a 1d topological sector that describes half-BPS operators in the theory. For certain theories, such as ABJM with k = 1, a Lagrangian was derived for this 1d sector that could be used to compute OPE coefficients of half-BPS operators as expectation values of n-body operators in the matrix model [29]. For small N , the integrals in these expectation values were evaluated explicitly in [30,31]. A certain half-BPS operator that appears in the stress-tensor correlator was related to ∂ 4 m i log Z m=0 using the 1d theory in [25], so this OPE coefficient could be computed to all orders in 1/N for any k using the Fermi gas result for Z(m i ).
In this paper, we will show how to use the Fermi gas method to compute these n-body operator expectation values to all orders in 1/N in terms of Airy functions. We will first apply these methods to ∂ 2 m i log Z m=0 and ∂ 4 m i log Z m=0 for ABJM theory, where we consider these quantities as expectation values of n-body operators with respect to the m i = 0 matrix model. Our results match those previously computed using the Fermi gas expression for Z(m i ), which serves as a non-trivial check of our methods. We will then consider c T as computed from ∂ 2 b log Z b=1 in the N f matrix model described above. We will compute the all orders in 1/N result for this quantity up to an additive N -independent constant, which in particular allows us to extract the first two leading orders at large N and finite N f . This is a new result that could be compared to bulk calculations in the future.
The rest of this paper is organized as follows. In Section 2, we discuss the calculation of n-body operator expectation values using the Fermi gas technique, and explain why the answers take the form of Airy functions. In Section 3 we then apply these methods to ABJM theory and recover the known results for ∂ 2 m i log Z m=0 and ∂ 4 m i log Z m=0 . In Section 4 we apply our methods to the N f matrix model to compute c T . We end with a discussion of our results and future directions in Section 5. Details of our calculations are given in various Appendices.

Fermi gas for n-body operators
We begin by reviewing the calculation of quantum-mechanical averages of n-body operators in a free Fermi gas, using the Wigner formalism. This is a classic subject that has been recently reviewed in [11], so we will be brief. For the n-body operators that we consider for 3d N = 4 matrix models, we show that these averages in the grand canonical ensemble are polynomial in the chemical potential µ, so that the N -dependence of the canonical ensemble averages are captured by Airy functions, as was previously shown for the partition function.
Let us first consider a system of N distinguishable particles, whose Hilbert space is spanned by the basis of eigenstates |x 1 , . . . , x N . An n-body operator O n is defined as an operator that is invariant under any permutation of the particles and acts on the Hilbert space as For a system of N identical fermions, we must anti-symmetrize the states using the projection operator As discussed in the Introduction, we want to compute the thermal averages of these operators in the canonical ensemble of a free Fermi gas with temperature T = 1. The canonical density matrixρ gas in a free gas factorizes into single particle density matricesρ whereĤ is the 1-particle Hamiltonian. For the thermal average of an n-body operator, we can integrate out (N − n) particles to write the average in terms of the n-particle density matrix ρ n as where O n is now restricted to the n-particle Hilbert space. These thermal averages are unnormalized, so that the canonical partition function Z N in (1.2) is by definition the average of the zero-body unit operator Z N = 1 .
The canonical ensemble is difficult to use in practice due to the restriction to a definite number N of particles. Instead, we can eliminate N using a Legendre transform to define the grand canonical ensemble with chemical potential µ, so that averages in each ensemble are related as For instance, the grand partition function Ξ corresponds to the average of the zero-body unit operator, and so takes the form As in the canonical ensemble, the density matrix of a free Fermi gas factorizes and can be written in terms of single particle Fermi distributions aŝ The grand canonical average of an n-body operator in a free Fermi gas can then be written in the n-particle Hilbert space as This average can be computed in a semiclassical small expansion using the Wigner formalism. The Wigner transform of an n-body operator O n is defined as where p i is the canonical conjugate of x i with the standard commutation relations The Wigner transform satisfies the product rule The utility of the Wigner transform is that it computes quantum statistical averages as phase-space integrals We can apply it to the grand canonical average in (2.8) to get where we integrated by parts n times to remove the in the products of O n andρ GC , and we have chosen to have the antisymmetrizer act on the states in the Wigner transform (2.9) of O n , since the density matrix is usually a more complicated operator.
The Wigner transform ofρ GC (Ĥ) involves two distinct small expansions. The first arises from taking the Wigner transform of a nontrivial function ofĤ. This can be computed using the so-called Wigner-Kirkwood expansion, which expresses the Wigner transform of any function f (Ĥ) as (2.14) The lowest couple Wigner-Kirkwood coefficients G r are trivially while G r for r ≥ 2 can be computed using the product rule (2.11) and have an expansion For instance, the lowest couple are and we show higher orders in Appendix C. The second expansion involved in the computation of (ρ GC (Ĥ)) W comes from computing the Wigner transform ofĤ itself. Let us assume that thep andx dependence ofĤ factorizes, as it does in all the matrix models we consider, so that the canonical density matrixρ iŝ The Wigner transform ofĤ is then computed as where log gives an expansion in 2 , and we show higher orders in Appendix C.
To sum up, the calculation of the canonical average O n of an n-body operator requires the following steps: 1. Compute the Wigner transform (P n O n ) W of the anti-symmetrized operator using (2.2) and (2.9).
2. Compute the Wigner transform (ρ GC (Ĥ)) W of the single particle grand canonical density matrix in a small expansion using the Wigner-Kirkwood (2.14) and quantum 3. Compute the 2n phase space integrals in (2.13) and multiply by the grand potential Ξ to get the grand canonical average O n GC (µ).
4. Recover the canonical average O n by performing the contour integral over µ in (2.5).
The grand potential Ξ for the 3d N = 4 theories that are amenable to the Fermi gas technique generically take the form where c i are various theory-dependent constants. One can perform the contour integral in (2.5) and use the identity to get the famous Airy function behavior of the partition function Z. The n-body operators that we consider will all have grand canonical averages that are polynomial in µ: we see that the N -dependence of the expectation values of O n will be completely captured by Ai(x) and Ai (x). We will demonstrate this for specific models in the following sections.

ABJM matrix model
We will start by discussing n-body operators in the N = 6 U (N ) k × U (N ) −k ABJM theory.
These theories describe N M2 branes probing a C 4 /Z k singularity and for finite k the holographic dual is described by M-theory on AdS 4 × S 7 /Z k . If we parametrize S 7 by complex coordinates on C 4 subject to 4 i=1 |z i | 2 = 1, then the Z k acts as In the large N and large k limit with fixed N/k the bulk description is Type IIA string theory on AdS 4 × CP 3 [3].

Matrix model
ABJM theory is an N = 6 SCFT, so it has an SO(6) symmetry as well as a U (1) flavor symmetry. In N = 2 language, ABJM theory consists of vector multiplets for each U (N ) k × U (N ) −k gauge group, as well as four chiral multiplets Z a , W a for a = 1, 2 which transform under the gauge groups and the SU (2) R ×SU (2) R ×U (1) flavor symmetry as shown in Table 1. The partition function on a squashed sphere with squashing parameter b and deformed by masses m i for the chiral fields can then be computed by assembling the standard N = 2 ingredients, as reviewed in [32], to get where Q = b + 1 b , the λ i , µ i correspond to the Cartans of the two gauge fields, and each chiral field contributes a factor with masses m i determined by the charge assignments in Table 1, so that m 1 corresponds to U (1) and m 2 + m 3 , m 2 − m 3 correspond to the Cartans of each factor in SU (2) × SU (2), respectively. The functions s b (x) are reviewed in Appendix A. If we restrict to the round sphere with b = 1, and set m 3 (or m 2 ) zero, then we can simplify the partition function using identities in A and write it in terms of m ± = m 2 ± m 1 (or As shown in [33], the partition function can furthermore be simplified using the Cauchy determinant formula to take the form where we integrated out µ i and rescaled λ i → x i /(2πk). For k = 1, this is the matrix model for a U (N ) gauge theory with an adjoint hypermultiplet and a fundamental hypermultiplet, deformed by an adjoint mass m + and an FI parameter m − 2 , which is expected from mirror symmetry [33]. 1 For k > 1, this mathematical reformulation has no physical explanation.
Lastly, as shown in [24], one can use the Cauchy determinant formula to further write this partition function (up to an unimportant overall numerical factor) in the form of (1.1) with single particle density matrixρ where the canonical variables are normalized so that The grand partition function Ξ was then computed in [24] using the standard Fermi gas methods to all orders in to get where the constant map function A is given by [35] A and in the second line we wrote A in the small k expansion, where B n+1 is a Bernoulli number.
The mass deformed partition function Z(m + , m − ) can then be obtained by computing the µ contour integral in (2.6) to get where the non-perturbative corrections come from both the e −µ terms that were left undetermined in (3.7), which translate into e − √ N k terms, as well as e − √ N/k terms that are invisible in the small k expansion of the Fermi gas method.
As discussed in the Introduction, taking n derivatives of the free energy ∂ n m ± log Z(m + , m − ) m=0 for ABJM gives n-point functions of stress tensor multiplet operators integrated on the sphere S 3 . For instance, taking two derivatives we find that c T in the normalization of [25] is (3.10) Taking four derivatives, we get the integrated stress tensor multiplet four point function. As discussed in [25,36], this quantity can be considered in the 1d protected subsector, where only a finite number of protected operators appear. It can then be used to compute the OPE coefficient of a certain protected operator in the 84 of SO(6) R as which from (3.9) can also be written to all orders in 1/N as a complicated sum of Airy functions and their derivative, similar to c T , whose explicit form we will not use.
Instead of computing these mass derivatives using the all orders in 1/N Airy function expression for the mass deformed partition function in (3.9), one could write them as n-body operators in the m ± = 0 matrix model with Hamiltonian (2.18) by taking derivatives of the exact mass deformed partition function (3.4). In particular, we find that ∂ 2 m ± Z m ± =0 can be expressed as (n ≤ 2)-body operators (3.13) Similarly, we can express ∂ 4 m − Z m ± =0 as (n ≤ 4)-body operators: as well as a similar more complicated operator from taking four derivatives of m + . From for the 2-body operators and for the 4-body operator, where the small k expansion of A(k) is given in (3.8). In the next few subsections we will use our n-body formalism to recover all the polynomial in µ terms shown here as well as the constant term to the lowest few orders in k. This will test all aspects of our formalism, which we will then use in the next section to compute completely new physical data.

∂
We start by taking the Wigner transform (2.9) of the antisymmetrized one-and two-body operators for O − in (3.13) in their respective n-particle Hilbert spaces to get where x ij ≡ x i − x j and p ij ≡ p i − p j . We then compute the phase space integral (2.13) for both operators to get where we denote ρ W ≡ (ρ GC (µ)) W (x, p), the leading term in (P 2 O − 2 ) W vanished since it is odd in each variable, and the delta functions allowed us to simplify the 2-body integrals. We will be able to compute the expansion of the first term A analytically, while we will need to use numerics for the second term B. The reason A can be computed analytically is that the products of ρ GC (µ) = 1 e H−µ +1 can be simplified as where in the second equality we used the Sommerfeld expansion as shown in [11] to write ρ in terms of derivatives of the Heaviside function θ. We can then apply the Wigner-Kirkwood expansion to the simpler function θ(µ −Ĥ) to express A as dxdp 2π 20) where recall that the coefficients G r have an expansion given in (2.16), and the quantum Hamiltonian H W also has an expansion (C.2). The resulting phase space integrals are very similar to those performed in the original Fermi gas paper [10], and as shown in Appendix We can then expand the Wigner-Kirkwood coefficients and the quantum Hamiltonian in to get two-dimensional integrals at each order in that we write explicitly in Appendix B. These integrals involve derivatives of ρ GC that we were unable to evaluate analytically.
Instead, we computed them numerically for many values of µ to high precision to get Finally, we combine A and B to find which matches the expectation (3.15) to the order shown.

∂ 2 m + Z m=0
Next, we consider the other two body operator O + , which is expected to have the same expectation value as O − . The integrals that appear in this calculation are more challenging than those of O − , so we will work just to leading order in . The operator O + consists of a zero body operator O + 0 whose expectation value is trivial, as well as the 2-body operator O + 2 whose anti-symmetrized Wigner transform is The phase space integral of O + 2 is then where in the second line we expanded sech 2 (x) for large x and replaced ρ W by its classical value to leading order in . We then evaluate this integral numerically and set = 2πk to (3.27) − π 2 4 N term, but N is strictly speaking not defined in the grand canonical ensemble. The physical quantity is the canonical ensemble average, which we get by taking the µ integral in (2.5) of (3.27) and adding − π 2 4 N to get where the second line followed from applying integration by parts using the definition of Ξ in (3.7) for m ± = 0. The term in the square brackets matches the expected grand canonical average in (3.15) to leading order in k.

∂ 4 m − Z m=0
Finally, we consider the expectation value of the (n ≤ 4)-body operator O 4 defined in (3.14), whose grand canonical expectation value is expected to take the form (3.16). As with the previous subsection, we will consider this calculation to leading order in . From taking the anti-symmetrized Wigner transforms of the operators in (3.14), we find that the leading order in contribution comes from the terms

(3.29)
We use these to compute the leading order phase space integral

N f matrix model
We will now discuss n-body operators in the so-called N f matrix model introduced in [15], which is an N = 4 theory of a U (N ) vector multiplets coupled to one hypermultiplet in the adjoint representation and N f hypermultiplets in the fundamental. These theories describe N M2 branes probing a C 2 × C 2 /Z N f singularity and for finite N f the holographic dual is described by M-theory on AdS 4 × S 7 /Z N f . If we parametrize S 7 by complex coordinates on In other words, we have M2s probing a KK monopole of charge N f . In the large N and large N f limit with fixed N/N f the bulk description is Type IIA string theory, where the KK monopole is realized by N f D6 branes wrapping the fixed loci of Z N f , which is AdS 4 × S 3 [13,14].
We will first review the matrix model of this theory, and discuss how c T can be computed from the matrix model partition function deformed by a squashing parameter b, by taking two derivatives of b and setting b = 1. This defines a 2-body operator that we will compute in a small = 2π expansion to several orders, which will end up corresponding to a small N f expansion. The corresponding canonical average can then written as Airy functions that capture the full 1/N expansion to the lowest few orders in N f . Unlike the previously discussed ABJM cases, this is a new result.

Matrix model
The deformed by masses for the matter fields as well as a FI parameter for the gauge field can then be computed by assembling the standard N = 4 ingredients in e.g. [32] to get As shown in [15], one can use the Cauchy determinant formula to further write this partition function (up to an unimportant overall numerical factor) in the form of (1.1) with single particle density matrixρ where the canonical variables are normalized so that = 2π . The Hamiltonianρ = e −Ĥ can then be written aŝ where these are the same U (x) and T (p) that appeared in the ABJM Hamiltonian (3.12).
Unlike ABJM, for the N f model does not depend on any physical parameters, so the small expansion would seem to be unphysical. In practice, however, one finds that the polynomial in µ terms in grand canonical averages only receive a finite number of corrections, just as in the ABJM case, so the small expansion still serves as a useful way of deriving the Airy function behavior of matrix model observables, which captures the all orders in 1/N expansion. Furthermore, the corrections tend to correlate with N f corrections. For instance, the grand partition function Ξ was computed in [15,16] using the standard Fermi gas methods to all orders in to get where this is the same constant map function A that was defined for ABJM in (3.8). For To compute c T , however, one must take two derivatives of the squashing parameter b as in (1.3), and the Fermi gas formalism cannot be applied to (4.2) for finite b. 2 Instead, we will apply our formalism to the (n ≤ 2)-body operator we derive from (4.2): (4.8) In the next subsection we will use our n-body formalism to compute all the corrections to the polynomial in µ terms in the grand canonical averages of these expressions, as well as the constant term to the lowest few orders in , which will correspond to the lowest few The operator O b consists of a zero body operator O b 0 whose expectation value is trivial, as well as the 1-body O b 1 and 2-body O b 2 operators whose anti-symmetrized Wigner transforms in their respective n-body Hilbert spaces are The phase space integrals for csch (x 12 ) 2x 12 − 2x 2 12 + π 2 coth (x 12 ) + π 2 csch (x 12 ) + δ (x 12 ) π 2 p 12 csch p 12 2 + 1 cosh p 12 2 + 1 , where we applied the Wigner-Kirkwood expansion to the densities. We can then expand the Wigner-Kirkwood coefficients and the quantum Hamiltonian in for the n-body operators to get 2n-dimensional integrals at each order in that we write explicitly in Appendix B.
We computed these numerically for many values of µ to high precision to get 3 where we set = 2π and we were unable to guess analytic formulae for the constant in µ terms. Note that just as in the ABJM case, the polynomial in µ terms only received a finite number of corrections, while the constant in µ term receives corrections to all orders.
We can now get the canonical average c T = O b by taking the µ integral in (2.5) and adding the zero-body term − π 2 4 N to get where we performed the µ integral as in (2.21) using the definition of Ξ in (4.7), and the constant term c was determined to the first couple orders in N f : If we set N f = 1 in (4.12), then we get the ABJM value of c T in (3.10) as expected, up to the constant term c that is not yet completely fixed. We can also expand (4.12) at large N to get where note that these first two terms do not depend on c. The leading term in fact matches the leading large N term for the ABJM c T if we identify N f = k, which is expected since to leading order in Newton's constant supergravity on

Conclusion
In this work we showed how OPE coefficients of protected operators in 3d N = 4 gauge theories with rank N can be computed as statistical averages of n-body operators in a free Fermi gas, which can be computed to all orders in 1/N in terms of Airy functions using the WKB expansion. This generalizes the Fermi gas formalism originally derived in [10]  shown in (4.14). It would be nice to compare this to a bulk calculation in the future, or to the 3d N = 4 numerical bootstrap results of [37]. It would also be useful to generalize the lowest few orders in N f numerical result for c shown in (4.13) to an all orders analytic result, perhaps in terms of the constant map function (3.8) as was shown in the related context of the sphere free energy [16].
The Fermi gas calculations in this work resemble those of the previous work [11] on expectation values of Wilson loops in the fundamental representation, which could be written as averages of one-body operators in the Fermi gas, but are much simpler. While the grand canonical averages of Wilson loops were exponential in the chemical potential µ, so that all orders in the WKB expansion were required to derive the Airy function behavior, the grand canonical averages of the n-body operators that appear in our formalism always take the form of polynomials in µ of finite degree, so only a finite number of terms need be computed to derive the all orders in 1/N Airy function behavior. Furthermore, in many cases the coefficients of the polynomial in µ terms receive only a few perturbative corrections in the WKB expansion, as was previously observed in Fermi gas calculations of the sphere free energy [10]. On the other hand, a technical difficulty that occurs in the case of n-body operators for n > 1 is that there is no Sommerfeld-like expansion of products of Fermi distributions in a large µ expansion that can be used to simplify the phase space integrals, as was shown for one-body operators in [11]. As a result, we were forced to compute many of the phase space integrals numerically. In most cases, we found that the numerical results were consistent with simple analytic expressions to high precision, which suggests that an analytic formula for the large µ expansion should exist.
There are many future applications of the formalism in this work. As discussed in the introduction, all 3d N = 4 CFTs contain a 1d topological sector that can be used to compute protected OPE coefficients in terms of the n-body matrix model averages that we consider.
For instance, in [31] the protected OPE coefficients that appear in the correlators of low lying half-BPS four-point functions in U (N ) 1 × U (N ) −1 ABJM theory were written explicitly as (n ≥ 4)-body operators. If these quantities could be computed to all orders in 1/N using our formalism, then the result could be compared to the numerical bootstrap result of [31], as was done for the correlator of stress tensor multiplet operators in [25]. One could also use our formalism to compute the integrated correlators of n conserved currents in 3d N = 4 CFTs, which as discussed in the introduction are related to n-body operator expectation values in the CFT's matrix model. These integrated correlators could then be used to constrain the 3d N = 4 numerical bootstrap [37], as was done in a related 2d context in [38].
The Fermi gas expansion is ideal for studying the M-theory regime where N is large and another physical parameter, such as k for ABJM or N f for the N f theory, is finite. In the large N and large k or N f regime with finite k/N or N f /N , which is called the 't Hooft limit, the Fermi gas results no longer apply. Instead, other methods such as topological recursion [39,40] have been used to compute quantities such as the sphere partition function and the expectation value of Wilson loops [9]. It would nice to apply these methods to the n-body operators considered in this work, which would complement the Fermi gas analysis and provide a window on a larger range of physical parameters. More ambitiously, if one could compute the non-perturbative corrections in both methods, then one might able to find exact expressions for these OPE coefficients, as was previously achieved for the ABJM sphere partition function [9,12,[41][42][43][44][45].

A The double sine function
The double sine function s b (x) (for reviews see for instance [34,46]) is defined as where the integration contour evades the the pole at t = 0 by going into the upper half-plane.
This function obeys several identities: The last identity is important when simplifying the partition function on the round sphere b = 1.

B Phase space integrals
In this appendix we give details for the phase space calculations in the main text. The analytic calculations are similar to those performed in [10]. We first discuss the phase space integrals for the ABJM matrix model, and then for the N f matrix model.
dxdp 2π We perform the calculation in an expansion in . The integrals which contribute order-byorder in are: 1. The 1 term in A is given by the integral 2. The term is given by the integral 3. Finally, the 3 term comes from the integral The expressions for H W and G r to order 4 are given in Appendix C.
These integrals can be done using the methods described in [10]. It turns out that all Wigner-Kirkwood contributions above vanish at this order, so we can set G i = 0. Thus the only non-vanishing phase-space integral we need to compute to order 4 is Since we only need derivatives of this integral, we will only calculate the µ-dependent terms.
Following [10], we can write this integral as n = 4(I 1 + I 2 ) , (B.6) Where we have split the phase-space integral into two regions: x(µ, p) 3 dp , (B.7) where p(µ, x) is the solution to H W (x, p) = µ in terms of p, and similarly for x(µ, p). In addition, x * (µ) = p * (µ) = µ up to exponentially small corrections. We briefly explain how this is done. First, due to symmetry under x → −x, p → −p, we can restrict ourselves to x, p ≥ 0. We are looking for solutions to the equation with H W defined in Appendix C. This equation defines a curve in the space (x, p). Consider the point (x * , p * ) along the curve such that p * = µ. Then from the form of H W we see that x * = µ as well, up to exponentially small corrections in µ. We plot the integration region in Figure 1. The phase space integral n can now be calculated as a sum of the integral over the regions X 1 , X 2 , X 3 . Instead, we would like to perform the integral separately in the region where p ≤ µ (which is X 2 ∪ X 3 ) and the region where x ≤ µ (which is X 1 ∪ X 2 ). The integral over X 2 ∪ X 3 is precisely what we called I 1 . As for the integral over X 1 ∪ X 2 , we would like to call this I 2 , but it is clear that we are overcounting since we are integrating over X 2 twice.
To remove its contribution, we simply subtract the integral over X 2 , which is simply the square integral Subtracting this from the integral over X 1 ∪ X 2 , we are left precisely with I 2 . Thus the sum 4(I 1 + I 2 ) is exactly the phase-space integral n.
We can find p(µ, x), x(µ, p) using an expansion in (using the expression for the Hamiltonian in Appendix C). We start by calculating I 1 . Up to exponentially small corrections, we can write We can now perform the integral I 1 . At leading order in we need to perform the following integrals: x * (µ) 0 (2µ − x) 3 dx, which can be calculated explicitly and produces 5µ 4 4 .
Similar considerations allow us to calculate the rest of the integrals at subleading order in . We find and we can use similar considerations to calculate I 2 . This calculation is much simpler, since the only µ-dependence appears at leading order in . We find I 2 = µ 4 12 up to µ-independent terms.
Summing I 1 and I 2 , we find 2. At order 3 the integrals that contribute are:
The expressions for the Hamiltonian H W and the coefficients G r appear in Appendix C. We computed these integrals numerically for many values of µ to get the result in the main text.