Fermi-gas correlators of ADHM theory and triality symmetry

We analytically study the Fermi-gas formulation of sphere correlation functions of the Coulomb branch operators for 3d $\mathcal{N}=4$ ADHM theory with a gauge group $U(N)$, an adjoint hypermultiplet and $l$ hypermultiplets which can describe a stack of $N$ M2-branes at $A_{l-1}$ singularities. We find that the leading coefficients of the perturbative grand canonical correlation functions are invariant under a hidden triality symmetry conjectured from the twisted M-theory. The triality symmetry also helps us to fix the next-to-leading corrections analytically.

1 Introduction The technique of supersymmetric localization allows for exact computations of observables in quantum field theories with a certain minimum amount of supersymmetry as pioneered by Pestun [1] (see [2] for reviews). It was shown in [3] that the partition function on a three-sphere S 3 of a threedimensional gauge theory with N ≥ 3 supersymmetry reduces to a matrix model via the localization.
For three-dimensional theories with N = 4 supersymmetry the sphere partition functions can be decorated by either of two types of half-BPS local operators; the Coulomb (resp. Higgs) branch operators whose expectation values define the Coulomb (resp. Higgs) branch. As shown in [4][5][6], the localization also allows for the evaluation of these protected correlators as generalized matrix integrals in such a way that a collection of the Coulomb or Higgs branch operators localize along a specific great circle S 1 in the S 3 .
The Ω deformation is very useful in the study of the supersymmetric gauge theories (See e.g. [7][8][9][10][11][12]). The protected correlation functions of the Higgs branch operators or the Coulomb branch operators in 3d N = 4 supersymmetric field theories are encoded by one-dimensional topological quantum mechanical models [13,14]. The associated topological quantum mechanical models arise from certain Ω deformations of the parent three-dimensional N = 4 supersymmetric field theories [15][16][17]. The topological quantum mechanics can be viewed as a non-trivial space of solution to OPE Ward identities equipped with the quantized Coulomb (resp. Higgs) branch algebra [15,18] which results from the quantization of the chiral ring of the Coulomb (resp. Higgs) branch operators.
It is shown in [19] that these sphere correlation functions of the Coulomb or Higgs branch operators as well as the sphere partition functions can be algebraically presented from the quantized Coulomb and Higgs branch algebras in terms of the twisted traces over the Verma modules without relying on the UV data.
The 3d N = 4 superconformal field theory (SCFT) appearing at low energy on a stack of N M2-branes on R × C at an A l−1 singularity probing the space C 2 × (C 2 /Z l ) has a UV description as a 3d N = 4 ADHM gauge theory with a gauge group U (N ), an adjoint hypermultiplet and l fundamental hypermultiplets [20,21]. In the near horizon limit of the M2-branes, one obtains the holographic dual AdS 4 × (S 7 /Z l ) background of M-theory, which provides us with an attractive example of the AdS/CFT correspondence [22][23][24]. For l = 1, the ADHM theory is self-mirror [25,26] and equivalent to the ABJM theory with k = 1 in the IR [27].
The large N limit of the partition function in the ADHM theory was studied well in [28,29]. In spite of the explicit expressions of the partition function, it is still tricky to evaluate it analytically in the large N regime. One of the analytic approach for the partition function is the Fermi-gas formalism [30] where the partition function is rewritten as the partition function of an ideal Fermi gas of N non-interacting particles. In the large N expansion of the free energy F = − log Z S 3 for the SCFT of the N M2-branes where Z S 3 is the sphere partition function, the leading coefficient can be evaluated from the two-derivative supergravity [31,32] and the next-to-leading coefficient is expected to be reproduced from higher derivative corrections in the supergravity. 1 Recently it has been proposed in [34] that the topological quantum mechanics that encodes the protected correlation functions of the world-volume theory of M2-branes are holographically dual to a certain protected sector of M-theory, that is the Ω-deformed or topologically twisted M-theory as an example of "twisted holography" [35].
Topologically twisted M-theory on an Ω-deformed background obeying the Calabi-Yau condition 1 + 2 + 3 = 0 can lead to a 5d theory on R × C 2 /Z l [36,37]. It is called the "twisted M-theory". The twisted M-theory is locally trivial and it is topological in R and holomorphic in the remaining four directions. It depends on the ratio 2 / 1 and it has a perturbative description as a non-commutative Chern-Simons theory at least for l = 1 and in some range of parameters. One interesting feature of this theory is a triality symmetry that permutes the Ω-deformation parameters [38] 1 → 2 , 2 → 3 , 3 → 1 . (1. 2) The twisted M-theory can contain M2-branes and M5-branes as line operators and surface operators in the 5d theory. 2 In the Ω-deformed background (1.1), when a stack of N M2-branes are placed on R × C 1 , the ADHM theory would acquire the mass parameter m for the adjoint hypermultiplet given by (1. 3) and it can be effectively described at low energy as topological quantum mechanics on R equipped with certain spherical part of the cyclotomic rational Cherednik algebras [42] with 1 being the quantization parameter as the quantized Coulomb branch algebra. The perturbative part of the protected correlation functions in the ADHM theory is expected to be holographically dual to a perturbative twisted M-theory background. It is shown from numerical and algebraic calculations in [43] that the perturbative part of protected correlation functions of the ADHM theory with l = 1 on a three-sphere enjoys the triality symmetry in the large N limit. In this paper we analytically evaluate the sphere correlators of the Coulomb branch operators for the ADHM theory in the Fermi-gas formulation. As argued in [44], the large N behavior of the correlation functions can be evaluated from averages of manybody operators in the Fermi-gas. The triality symmetry (1.2) that is associated to the adjoint mass (1.3) is a key in our analysis. This symmetry should be manifest in the perturbative part of the correlation functions in the large N limit [43]. In fact we show that the leading terms of the perturbative grand canonical correlation functions are generally invariant under the triality symmetry (1.2). Moreover, the triality symmetry constrains the forms of the perturbative part. In the Fermi-gas formulation, it is technically hard to compute the subleading terms of the perturbative part analytically. We are able to obtain it partially. The triality symmetry enables us to reconstruct the remaining missing pieces. In this way, we obtain consistent triality invariant subleading terms for higher-point functions.
The organization of this paper is as follows. In the next section, we review the Fermi-gas formulation by following the original argument in [30]. The technique there is extended to the correlators of the Coulomb branch operators, as shown in Section 3. We will derive the large N behavior for the multi-point correlators. Finally, we will give some remarks on related topics in Section 4.

Fermi-gas formulation
We start with a review of the Fermi-gas formulation [30] to analyze partition functions in 3d supersymmetric field theories on S 3 . In the Fermi-gas formulation, it is more convenient to go to the grand canonical ensemble rather than the canonical one. We show how to derive the grand potential in the large chemical potential limit.
Supersymmetric localization reduces path integrals to matrix models [3]. The partition function of the ADHM theory on S 3 takes the form where m is the mass of the adjoint hypermultiplet, and ζ is the Fayet-Iliopoulos (FI) parameter. For l = 1, this is equal to the matrix model of the ABJM theory with k = 1 [27]. Making use of the Cauchy identity, one can identify the matrix integral (2.1) with the canonical partition function of a non-interacting, one-dimensional Fermi-gas with N particles [30]: is the one-particle density matrix in the position representation. This leads to a systematic analysis of the large N limit of the partition function on S 3 as a thermodynamic limit of an ideal Fermi-gas. The analysis of the Fermi-gas of the partition function below has appeared in [45]. We generalize it to the computation of some correlation functions, and find new features on a manifestation of the hidden triality symmetry, expected in [38,43]. The thermodynamic limit of an ideal Fermi-gas can be obtained by considering the one-particle problem in the semi-classical approximation and the 1/N corrections to the thermodynamic limit can be obtained by evaluating the quantum corrections to the semi-classical limit. In the following discussion, we consider the grand canonical ensemble, in which the grand potential is introduced by Our goal in this section is to derive the large µ limit of J(µ) by using the Fermi-gas formulation. Letσ andp be canonically conjugate operators obeying [σ,p] = i (2.5) where = 1 2π . Then we can write the density matrix operator aŝ so that the kernel (2.3) can be realized as the matrix element of the operator (2.6) in the position space. From the density matrix operator (2.6), we can define a one-body HamiltonianĤ of a system by e −Ĥ :=ρ . (2.9) In the classical limit, this Hamiltonian reduces to where U (σ) is a potential term and T (p) is a kinetic term. In the following we analyze the Fermi-gas system in the case where the FI parameter ζ and the adjoint mass parameter m are pure imaginary so that the Hamiltonian is real. However, we expect that the expression with arbitrary ζ and m can be reached from our result by analytic continuation. Given a density matrix (2.3), one can quantize the Fermi-gas system by following the phase space formulation that is distinguished from the canonical quantization and the path integral formulation. The phase space quantization is based on the Wigner-Weyl transforms and the Weyl correspondence between c-number functions in the phase space and quantum mechanical operators in the Hilbert space so that quantum mechanical composition of functions relies on the star-product.
The Wigner transform of an operatorÂ with its matrix elements σ|Â|σ = A(σ, σ ) in the position space is the function [46][47][48][49] in the phase space. This maps a quantum mechanical operatorÂ in the Hilbert space to a function in the phase space. The inverse operation is the Weyl transform which relates a function B W (σ, p, ) in the phase space to a quantum operatorB in the Hilbert space with matrix elements . (2.12) When we deal with more than one particle, we need to include the effects of quantum statistics in the Wigner transform. For the Fermi-gas, the Wigner transform of the sbody operator O (s) is obtained by taking the anti-symmetrized operators P A O (s) where P A is the projection operator which anti-symmetrize the states [44]. The Wigner transform of a product of operatorsÂ andB is given by [46][47][48][49] (ÂB) W = A W B W . (2.14) Here is the star operation where the derivatives act on the left or on the right according to the directions indicated by the arrows. One can express all operators of quantum mechanics in terms of the Wigner transforms of operators in such a way that the semi-classical expansion of an operatorÂ is given by where A 0 is the classical limit ofÂ. From the Baker-Campbell-Hausdorff formula, (2.9) and (2.14) we find the Wigner transform of the Hamiltonian Furthermore, the semi-classical expansion of arbitrary function f (Ĥ) of the Hamiltonian operatorĤ is given by where f (r) (a) is the r-th derivative of f (x) evaluated at x = a and is the universal coefficients in the expansion around H W . It follows that G r is an even function of such that for the largest integer n < r 3 . For example, we have [50, 51] G 0 = 1, (2.21) In particular, the Wigner transform of a distribution operator at zero temperature is given by [50] θ W 2πµ where θ(x) is the Heaviside step function. By taking the trace of the distribution operator, we get the function n W (µ) that counts the number of eigenstates whose energy is less than 2πµ 1 . In the thermodynamic limit N → ∞, the behavior of the system is semi-classsical and the trace can be evaluated as an integral over the phase space. Therefore we obtain (2.23) Here the first term is the area of the quantum corrected Fermi surface defined by the equation and the second term is the quantum corrections arising from the semi-classical expansion of the distribution operator. Then the density of energy eigenstates is given by Let us first evaluate the area of the quantum corrected Fermi surface in the limit µ → ∞:  the potential term U (σ) and its derivative have the asymptotics and the kinetic term T (p) and its derivative have the asymptotics Let (σ + * , p + * ), (σ + * , p − * ), (σ − * , p + * ) and (σ − * , p − * ) be points in the quantum curve where .

(2.30)
It then follows from (2.24), (2.28) and (2.29) that where the exponentially small corrections in µ are power series in 2 . We divide the quantum corrected Fermi surface into four domains: O(e −σ ) for III (2.33) and the Wigner transform of the Hamiltonian Therefore we can solve σ along the quantum curve On the quantum curve in the regions II and IV, the exponential terms in p are larger than those in µ so that the kinetic term and its derivatives become for IV (2.36) and the Wigner transform of the Hamiltonian reduces to Thus we can solve for p along the quantum curve Putting the area of the quantum Fermi surface of the region I, II, III and IV together, we finally obtain the quantum corrected area Vol(µ) of Fermi surface 3 (2.39) The quantum corrections that correspond to the second term in (2.22), which are associated with the semi-classical expansion of a function of the Hamiltonian turn out to yield only the non-perturbative corrections of order e −µ [30]. Therefore we obtain where n np (µ) = O(µe −µ ) denotes the non-perturbative terms. Since our Hamiltonian is positive, it follows that n W (0) = 0 and n np (0) = −n 0 after resumming all the nonperturbative corrections.
To obtain the leading and next-to-leadingg coefficients that show up in the free energy, we observe that the grand canonical potential can be expressed in terms of the density (2.25) of eigenstates: According to the asymptotics of the trilogarithm .
We see that the leading coefficient (2.44) and the next-to-leading coefficient (2.45) are actually invariant under the triality transformation (1.2). The overall factor 1 1 2 3 can be interpreted as the equivariant volume of the Ω-deformed planes C 1 × C 2 × C 3 in the background (1.1) of the twisted M-theory.
In the next section, we extend this computation to correlation functions for Coulomb branch operators.

Coulomb branch correlators
The 3d N = 4 supersymmetric gauge theory generically contains two types of half-BPS local operators, i.e. the Coulomb and Higgs branch operators, which parametrize two branches of supersymmetric vacua, the Coulomb and Higgs branches respectively. The Coulomb branch operators can be built out of the monopole operators v n * 4 dressed by the vector multiplet scalar fields ϕ. They can be expanded as a sum over the monopole operators where R n * (ϕ, m C ) are polynomials in ϕ. The sphere correlation function of the Coulomb branch operators for the ADHM theory takes the form 5 The sphere correlation functions of the Coulomb branch operators can be universally expressed in an algebraic way in terms of the twisted traces over the Verma modules of the quantized Coulomb branch algebra [19]. The factor R 0 (−iσ, −im) inserted in the correlation function is pulled back from generators in the quantized Coulomb branch algebra. Since the non-trivial twisted traces involving the monopole operators or equivalently shift operators can appear only when they simply shift the vector multiplet scalar fields, only some insertion of non-periodic part without the shift, i.e. R 0 (−iσ, −im) in the integrand will lead to distinct Coulomb branch correlation functions with some change of residues.

Quantized Coulomb branch algebra
There exist two types of Ω-deformations for the 3d N = 4 supersymmetric gauge theory, in which two kinds of non-commutative algebras of the topological Coulomb and Higgs branch operators emerge. They are called the quantized Coulomb and Higgs branch algebras [15,16]. The quantized Coulomb branch algebra A C N,l; 1 , 2 of the ADHM theory is isomorphic to the spherical part SH cyc N,l of the cyclotomic rational Cherednik algebra 4 The monopole operators are labeled by the GNO charge ±A a , a = 1, · · · , r where A is the cocharacter that specifies an embedding of a U (1) monopole singularity into the gauge group G and r is the rank of G. For G = U (N ) we have A = (A 1 , · · · , A r ) ∈ Z r and we denote the integer-valued charge by n * .
5 See [5,6] for the result of supersymmetric localization. [42]. The algebra A C N,l; 1 , 2 can be also identified with the shifted Yangian Y l (m i ) of gl(1) which is obtained by deforming the subalgebra of an affine Yangian Y ( gl(1)) [52].
Let us introduce coordinates w a and shift operators v a , v −1 a , a = 1, · · · , N which obey The algebra A C N,l; 1 , 2 is generated by the operator where B n (x) is the Bernoulli polynomial as well as raising and lowering operators which take the forms: Here the operator h n can be determined by the relation

d n operators
There is an alternative presentation of the algebra in such a way that all generators can take the form of 1 times a triality-invariant expression [38]. We can introduce the Hamiltonian operator associated to any polynomial f (w) in w.
When we choose polynomials where B n (x) is the Bernoulli polynomial, which satisfy the recursion relation 6 we obtain the operator It has a generating function where ψ(z) is the digamma function. For example, we have Making use of the operator d n given by (3.20), one can also build the other generators in the quantized Coulomb branch algebra A C N,l; 1 , 2 which take triality-invariant fashion up to the overall 1 factor, as discussed in [38]. It manifests the symmetry of the algebra under the triality symmetry (1.2).

Fermi-gas formulation
In terms of the Fermi-gas formulation, we can also evaluate correlation functions of the Coulomb branch operators. The treatment is very similar to the previous work [53] for Wilson loop correlators in ABJM theory. We can rewrite the sphere correlation function (3.2) as First, we consider the sphere one-point function of a positive power function: As we will see later, this is an important building block to compute the thermodynamic limit N → ∞ of the one-point function d n . To study the large N limit of the onepoint function (3.29), we evaluate an average of the one-point function by integrating over the phase space in terms of the Wigner transform of the distribution operator (2.22): where the second line involves the corrections from the quantum Fermi surface and the third is associated to the corrections from the semi-classical expansion of the distribution operator.
The corrections from the quantum Fermi surface can be evaluated in the same manner in the previous section. Collecting the pieces (B.6), (B.7) (B.11) and (B.12), we obtain the leading and next-to-leading terms of the average over the quantum Fermi surface in the first line of (3.30): Next proceed to the G r corrections arising from the semi-classical expansion of the distribution operator. To order 2 corrections only come from G 2 and G 3 in (2.21): Consequently, only non-trivial corrections may come from G 3 . We find which involve µ n−2 and lower order terms. Hence the quantum corrections associated to the semi-classical expansion do not contribute to the leading and next-to-leading terms.
Putting all together, we finally arrive at For example, for l = 1 and ζ = 0 the expression reduces to a relatively simple form

Grand canonical one-point functions
In this and the next subsections, we would like to evaluate the large N limit of correlation functions of the operators d n . The k-point function of d n generically takes the form We go to the grand canonical ensemble, and study the large µ limit, as was done for the partition function. We can easily translate obtained results in the grand canonical ensemble into those in the canonical ensemble.
Let us consider the one-point function. According to the formula (3.18), the polynomial p n (σ) has the leading and next-to-leading terms: Thus the leading and next-to-leading terms of the one-point function in the thermodynamic limit can be obtained from (3.37) Making use of the Sommerfeld expansion one can express the one-point function of operator O in the grand canonical ensemble for an ideal Fermi-gas system as [54,55]

44)
By taking O = p n (σ) and β = 2π/ 1 in (3.43), we get from (3.41) the leading and next-to-leading terms of the grand canonical one-point function of the operator d n : = c n+2 µ n+2 + π 2 6 (n + 2)(n + 1) 2π We see that the leading and next-to-leading coefficients of the grand canonical one-point functions (3.45) are exactly invariant under the triality symmetry (1.2)! Indeed, it is simple to check that our analytic formula (3.45) of the grand canonical one-point function reproduces the numerical results in [43] when we specialize l = 1 and ζ = 0. We first encode the µ dependence into d 0 by replacing µ with τ 0 2π . Then we can obtain the perturbative correlation functions in [43] by taking the derivatives with respect to τ 0 and setting τ 0 to zero.
For example, the grand canonical one-point function of d 2 for l = 1 and ζ = 0 is given by The derivatives of (3.46) with respect to τ 0 lead to where

Grand canonical higher-point functions
Higher-point functions can be evaluated by taking the averages of many-body operators in the ideal Fermi-gas. The analysis is more involved than the one-point function.

Two-point functions
Let us see the grand canonical two-point functions of the operator d n . It has contributions from the following one-and two-body operators: From (2.11) we get the leading and next-to-leading terms of the Wigner transforms of the anti-symmetrized operators for (3.56) where the ellipsis indicates the terms at low orders in σ which do not contribute to the leading and next-to-leading coefficinets of the correlation functions. Plugging the Wigner transforms (3.57) and (3.58) into (3.55), we find the leading and next-to-leading terms of the two-point function In the second equality we have combined the average of one-body operator p n (σ) 2 with the average of the Wigner transform of the antisymmetrized two-body operator −2π δ(σ 1 − σ 2 )δ(p 1 − p 2 )p n (σ 1 ) 2 where = 1 2π so that they can be evaluated as the one-body integral involving ρ GC W (1 − ρ GC W ). The grand canonical connected two-point function of the operator d n can be obtained by subtracting the square of the grand canonical one-point functions. Thus we where we have used the relation in the second equality. We obtain from (3.60) and (3.37) the leading term in the connected two-point function: In fact, this is invariant under the triality symmetry (1.2)! The subleading terms also have contributions from the last term in (3.60). Unfortunately it seems difficult to evaluate it analytically because of the derivatives of ρ GC W in the integrands. However, we can guess a consistent next-to-leading term by requiring the triality invariance, We will check the validity of this guess below. It is straightforward to generalize the results (3.62) and (3.63) to the connected two-point functions for two distinct operators d n 1 and d n 2 with n 1 = n 2 by following the same argument. To make a result simpler, we introduce N 2 = n 1 + n 2 . The result is We should say again that the next-to-leading term in this expression is a guess based on the triality. Now we check that our analytic formula (3.64) of the grand canonical connected two-point function reproduces the numerical results in [43] as special cases with l = 1 and ζ = 0. For example, in a similar manner for the one-point function, we have by replacing µ with τ 0 2π where σ 2 and σ 3 are defined by (3.49) and (3.50). Then we get In fact, the leading coefficient (3.67) and the next-to-leading coefficients (3.68) precisely agree with the numerical results in [43]! 8 In order to verify our results further, we note that the grand canonical connected higher-point functions of the operator d 1 can be derived from the grand canonical potential (2.43) with non-zero FI parameter ζ = 0 since it can be viewed as a generating function of the correlation functions of the operator d 1 . By shifting the FI parameter ζ by τ 1 2π and expanding the grand canonical potential (2.43) in powers of τ 1 , we can extract the grand canonical connected higher-point functions of the operator d 1 From the expansion (3.71) we find that the leading and next-to-leading terms of (k +1)point functions of the operator d 1 are given by The result for k = 0 precisely agrees with (A.1) obtained from the formula (3.45) of the grand canonical one-point function. Also the result for k = 1 coincides with (A.5) obtained from the formula (3.64) of the grand canonical connected two-point function.
It is obvious to see that the grand canonical one-point functions (3.45) with nonzero FI parameter ζ = 0 can also be viewed as the generating function of the grand canonical connected correlation functions with an insertion of d n and an arbitrary number of d 1 . By replacing ζ with ζ + τ 1 2π in (3.45) and then expanding it in powers of τ 1 , we find Thus we obtain the grand canonical connected correlation functions of d n and an arbi-trary number of d 1 :   (3.64) and expanding it in powers of τ 1 , we obtain the leading and next-to-leading terms of higher-point functions:

Three-point functions
For the grand canonical three-point functions there are contributions from the one-, two-and three-body operators of the forms From the Wigner transforms of the antisymmetrized operators for (3.76) where the ellipsis stands for the terms which do not contribute to the leading term, we get the grand canonical three-point functions The grand canonical connected three-point function of the operator d n is obtained by subtracting the disconnected parts so that the leading term appears from the one-body integral. We obtain Here we have used the relation The next-to-leading term appearing from the one-body integral (3.81) is not still triality invariant since it involves further contributions from higher order terms in the Wigner transforms (3.77)-(3.79). Assuming the triality invariance, we find a consistent expression for the next-to-leading term i 3n 2 3n−5 (3n)(3n − 1) [(−1) n (l − 2iζ) 3n+1 + (l + 2iζ) 3n+1 ] 3 1 2 3 π 2 (l 2 + 4ζ 2 ) 3n+1 The same argument yields the connected three-point functions for generic three operators d n 1 , d n 2 and d n 3 where N 3 = n 1 +n 2 +n 3 . We find the triality invariant leading coefficient by analytically computing the one-body integral as in (3.81). Again the next-to-leading coefficient has additional contributions from the higher order Wigner transforms which we could not analytically evaluate. Instead of computing them explicitly, we obtain a consistent next-to-leading coefficient by restoring the triality invariance. As a consistency check of our expression (3.84) of the connected three-point function, note that when one of the three operators is taken as d 1 , say for n 3 = 1, it precisely coincides with the result obtained from (3.75) when k = 1.
We can also obtain the leading and next-to-leading coefficients of the grand canonical connected correlation functions with d n 1 , d n 2 , d n 3 and an arbitrary number of d 1 by shifting the FI parameter in (3.84) by τ 1 2π and expanding it in powers of τ 1 . We find

Four-point functions
The grand canonical four-point functions can be computed from the one-, two-, threeand four-body operators: The Wigner transforms of the antisymmetrized operators for (3.86) take the forms where the ellipsis indicates the terms which do not affect the leading term. Making use of (3.87)-(3.90), we can compute the grand canonical four-point func-tion of the operator d n (3.91) In particular, when the FI parameter is turned off, the one and three-point functions vanish so that the grand canonical four-point function (3.91) is simplified as The leading term that is proportional to µ 4n+2 appears from the disconnected term 3( d n d n GC /Ξ) 2 . By eliminating the disconnected terms from (3.91) and using the relation we obtain the grand canonical connected four-point function of the operator d n The leading term of the connected four-point function is proportional to µ 4n−1 . We see that the expression (3.94) is triality invariant! When n = 1, it coincides with the leading term in the previous result (3.72). The next-to-leading terms in the one-body integral (3.94) are not still triality invariant as there are additional contributions from the the higher order terms in the Wigner transforms (3.87)- (3.90). Assuming that the Wigner transforms complete the triality invariance so that it is consistent with (3.85), we find next-to-leading terms (3.95) More generally, we get the grand canonical connected four-point functions for d n 1 , d n 2 , d n 3 and d n 4 (3.96) where N 4 = n 1 +n 2 +n 3 +n 4 . For n 4 = 1 the expression (3.96) reproduces the connected four-point function obtained from (3.85) for k = 1 involving a single d 1 .
Again we can extract from (3.96) the leading and next-to-leading terms of the connected higher-point functions with additional insertion of d 1

k-point functions
We can compute more general higher-point functions by considering the many-body operators of d n . The grand canonical k-point function is obtained by summing over the averages of the Wigner transforms of the antisymmetrized l-body operators From the grand canonical k-point function (3.98) and the lower-point functions, the grand canonical connected k-point functions can be recursively computed as In particular, the leading terms of the grand canonical connected k-point functions can be easily evaluated by acting on the one-body integral with the differential operator 1 β k−1 ∂ k−1 µ . We find the leading term of the grand canonical connected k-point function of the operators d n i , i = 1, · · · , k: The leading term is proportional to µ N k −k+3 and its coefficient is proportional to 1 1 2 3 so that it is invariant under the triality symmetry (1.2). The next-to-leading terms appearing from (3.100) are not yet triality invariant. Provided that these altogether form the triality invariant expression, we get the consistent triality invariant next-to-leading terms of the grand canonical connected k-point function of the operators d n i , i = 1, · · · , k proportional to µ N k +1 : In fact, we have checked that the expressions (3.100) and (3.101) reproduce the leading and next-to-leading terms in all the previous results. It has been numerically found in [43] that the perturbative correlation function has a conjectural pattern in such a way that the non-vanishing terms have a power of σ 3 greater or equal to −1.
By setting µ = τ 0 2π in our results (3.98) and (3.101) and taking the derivatives with respective to τ 0 , we can write from them the perturbative correlation function in [43]: This is compatible with the conjectural pattern (3.102).

Concluding remarks
In this work, we evaluated the large N correlation functions of the Coulomb branch operators for 3d N = 4 ADHM theory in the Fermi-gas formulation. We confirmed that the leading perturbative part has the triality symmetry, expected from the dual twisted M-theory. Interestingly, the full analytic form of the next-to-leading order correction can be fixed by the requirement of this symmetry even though the Fermi-gas computation is technically hard at this order. This idea should be useful in higher order computations. We remark several related directions. The correlators have been studied in bootstrap program for 3d SCFTs which arise as the IR limit of the effective theory of multiple M2-branes [13,[58][59][60][61]. It would be nice to address the d n correlators via the bootstrap analysis to be compared with our results.
The twisted holography [35] (see also [34]) would relate our results to the perturbative calculations around a dominant semi-classical saddle point in the holographic dual five-dimensional holomorphic (symplectic)-topological theory on AdS 2 × S 3 . The k-point correlators that we have computed would correspond to an amplitude of the Feynman diagram in the 5d Chern-Simons theory with k-points insertion on the 1d defect, where the 1d topological quantum mechanics lives. Beyond the perturbative calculations, the triality symmetry may be broken due to the instanton corrections [43]. It would be nice to extend our Fermi-gas analysis by calculating the non-perturbative corrections to the correlators which capture the full geometry normal to AdS 2 × S 3 in AdS 4 × S 7 . Also our subleading terms could be used to test the holographic dual with higher derivative corrections as recently studied in [33,62].
The sphere correlators of the Coulomb and Higgs branch operators can be algebraically computed as a sum of the products of the twisted traces over the Verma modules [19]. For the ADHM theory, this sum is taken over a set of l Young diagrams with N k (k = 1, · · · , l) boxes obeying l k=1 N k = N . It would be interesting to analyze the large N behavior of the correlators in terms of the twisted traces.
In [43] it is conjectured from the numerical and algebraic computations that the generating function of connected correlation functions of the operator d n satisfies a recursion relation that leads to a quadratic constraint on the perturbative correlation functions. This is reminiscent of the "string equation" [63,64] in topological gravity and it may play a key role in the twisted holography. It is intriguing to give an analytical derivation or proof of this relation by extending our analysis.
The line operators in 3d N = 4 ADHM theory would be also realized in the twisted M-theory by introducing extra M2-branes intersecting with the original ones. The space of the local operators living at junctions of line operators is realized as Hom space which generalizes the bulk Coulomb and Higgs branch algebras [65]. It is interesting to figure out the triality symmetry in the presence of line operators by applying the Fermi-gas analysis as studied for the ABJM model in [53].