Exact results in a N=2 superconformal gauge theory at strong coupling

We consider the $\mathcal{N}=2$ SYM theory with gauge group SU($N$) and a matter content consisting of one multiplet in the symmetric and one in the anti-symmetric representation. This conformal theory admits a large-$N$ 't Hooft expansion and is dual to a particular orientifold of $AdS_{5}\times S^{5}$. We analyze this gauge theory relying on the matrix model provided by localization a la Pestun. Even though this matrix model has very non-trivial interactions, by exploiting the full Lie algebra approach to the matrix integration, we show that a large class of observables can be expressed in a closed form in terms of an infinite matrix depending on the 't Hooft coupling $\lambda$. These exact expressions can be used to generate the perturbative expansions at high orders in a very efficient way, and also to study analytically the leading behavior at strong coupling. We successfully compare these predictions to a direct Monte Carlo numerical evaluation of the matrix integral and to the Pad\'e resummations derived from very long perturbative series, that turn out to be extremely stable beyond the convergence disk $|\lambda|<\pi^2$ of the latter.


Introduction and summary of results
Four-dimensional gauge theories with extended supersymmetry are intensively explored to obtain exact results in quantum field theory and improve our understanding of the strong-coupling regime. Important progress in these directions has been achieved, over the years, in the case of N = 4 SYM theories by exploiting integrability, localization techniques and the AdS/CFT correspondence.
The analysis of cases with reduced supersymmetry is of course more difficult, but for N = 2 theories significant steps have been taken. In particular, for these theories localization techniques can still be used and some instances of holographic correspondence are known. In particular, as shown in [1], by putting a N = 2 SYM theory on a compact space like a four-sphere, one can localize the infinite-dimensional path-integral on a finite-dimensional locus and reduce the calculation to an interacting matrix model. If conformal symmetry is present 1 , this matrix model also encodes information on the observables of the theory in flat space. In this way the partition function and the vacuum expectation value of the BPS Wilson loop have been computed [1,[3][4][5][6][7][8]. Many other observables of the gauge theory can be described by the matrix model, such as the extremal correlators of chiral/anti-chiral fields [9][10][11][12][13][14][15][16][17][18][19][20][21][22], the correlators of chiral operators with the BPS Wilson loop [16,[23][24][25] and the brehmsstrahlung function [26][27][28][29]. The same localization techniques can also be applied to study various properties of the massive N = 2 * SYM theory [30][31][32][33] and of quiver theories [34][35][36][37][38]. These matrix-model computations can be carried out following two different strategies, called the Cartan sub-algebra and the full Lie algebra approach in [39]. In the first approach, the matrix model is written as an integral over the matrix eigenvalues and its large-N limit is described by the asymptotic eigenvalue distribution determined within a saddle-point approximation. The second approach, instead, keeps the matrix integral over the full Lie algebra and exploits recursion relations [18] to evaluate correlation functions. This technique, originally developed because of its effectiveness at finite N , has been applied also to study the large-N and the large-charge sectors of gauge theories [40][41][42][43].
For the N = 4 SYM theory with group SU(N ) the matrix model arising from localization is Gaussian, which makes relatively easy to treat it both at finite N and in the large-N limit, and also to discuss its strong-coupling behavior. On the contrary, for N = 2 SYM theories the matrix model has a very complicated potential and its use for an analytical description of the strong-coupling regime is far from obvious.
Previous strong-coupling results have been obtained in the Cartan sub-algebra approach. In particular, in the seminal work [5] the expectation value of the BPS Wilson loop in N = 2 SYM with N F = 2N fundamental flavors (SQCD) has been computed in the planar approximation at leading order for large values of the 't Hooft coupling λ, by solving an associated Wiener-Hopf problem. The planar free energy of SQCD has been later evaluated in [7] and generalizations to N = 2 superconformal theories with gauge group SO(N ) or Sp(N ) have been obtained in [13]. An important outcome is that the planar expectation value of the Wilson loop scales proportionally to λ η with an exponent η determined by the large-N limit of the ratio N F /N , confirming planar equivalence with N = 4 SYM (where N F = 0) in all models where N F does not grow with N . The strong-coupling limit of the chiral two-point functions at planar level has been initially studied in [15] and further explored in [17]. More recently, the Wiener-Hopf method has been successfully exploited to determine the strong-coupling expectation value of the BPS Wilson loop in the SU(N )×SU(N ) quiver theory [35]. Although still at planar level, this calculation is at the next-to-leading order in the 't Hooft coupling, and represents an important extension of the leadingorder results of [4,5]. Despite attempts to simplify the mathematical structure of the Wiener-Hopf approach (see for example [6]), this method has not been extended beyond the planar level, the main technical obstruction being that it is not based on a controlled expansion in terms of a parametrically small quantity.
The study of the strong-coupling regime in the large-N limit for N = 2 superconformal theories is clearly an important ingredient also to achieve a full understanding of the AdS/CFT correspondence in systems with a reduced amount of supersymmetry, at least for those models that possess an holographic dual that can be accessed at strong coupling from the gauge side.
To address these issues it is convenient to consider N = 2 theories that are as close as possible to N = 4 SYM. Therefore, in this paper we focus on a particular N = 2 conformal model with gauge group SU(N ), called E theory in [22,25], which is very close to N = 4 SYM in many ways. In the 't Hooft limit, it possesses a large class of observables, comprising the free energy and the expectation value of the BPS Wilson loop, which agree with the corresponding ones in the N = 4 SYM at the planar level [22,39] and deviate from the N = 4 results only at order 1/N 2 . Moreover, it admits a holographic dual of the form AdS 5 × S 5 /Z 2 which is obtained by means of a Z 2 orbifold/orientifold procedure from the holographic dual of N = 4 SYM (see for instance [44]). Actually, the E theory itself can be realized by taking a suitable orientifold projection of a twonode quiver model, which in turn can be engineered with a system of fractional D3-branes in a Z 2 orbifold [45]. It is natural to distinguish the observables of the E theory in two classes, those corresponding to operators that belong to the untwisted sectors in the string construction, and those associated to operators that belong to the twisted sectors. For this reason, throughout this paper we will use the terminology of "untwisted" and "twisted" observables.
In the holographic dictionary an untwisted observable of the E theory, which we generically denote by U , is described by a supergravity excitation of AdS 5 × S 5 which is even under the orbifold/orientifold parity. Therefore, the same excitation exists also in the maximally supersymmetric case where it corresponds to an observable of the N = 4 SYM theory which we may call U 0 . One then expects that the E-theory result be identical to the one in N = 4 at the leading order in the large-N expansion, with a deviation starting at the non-planar level: where δ U is a non-trivial function of the 't Hooft coupling. The free energy and the vacuum expectation value of the BPS Wilson loop are examples of untwisted observables [8,22,39] but, as we will see, also the correlators of gauge invariant single-trace operators with even dimension in the E theory belong to the untwisted class.
In the N = 4 theory there are observables, which we call T 0 , that in the holographic correspondence are mapped to supergravity excitations of AdS 5 × S 5 which are odd under the Z 2 orbifold/orientifold parity. These excitations are removed by the projection and get replaced by states of the twisted sectors corresponding to the twisted observables T of the E theory. Therefore, one expects that T and T 0 differ already at the planar level: where δ T depends non-trivially on the 't Hooft coupling λ. As already discussed in [22], the correlators involving of gauge invariant single-trace operators with odd dimension in the E are of twisted type and behave as in (1.2). The main goal of this paper is to provide exact results for the functions δ U and δ T for a large class of observables, to exhibit their exact dependence on the 't Hooft coupling λ, and to find eventually their strong-coupling behavior when λ → ∞.

Results
By exploiting the simplicity of the matrix model and relying on the power of the full Lie algebra approach, we have analyzed in detail several observables of the E theory in the large-N limit in the untwisted and in the twisted sectors, and found the leading corrections with respect to the N = 4 theory using a simplified set of recursion relations [22]. As a consequence, we show that these quantities can be effectively evaluated in a Gaussian matrix model whose quadratic form is an infinite λ-dependent matrix 3) whose elements are convolutions of Bessel functions of the first kind. More precisely, we find that for the untwisted observables we have considered, the deviation δ U is given in terms of logarithmic λ-derivatives of log det 1 − X , (1.4) which is proportional to the free energy of the effective Gaussian model, while for the twisted variables the deviation δ T is expressed in closed form in terms of the propagator Since the matrix X is exactly known as a function of λ through the Bessel functions entering its definition, we can expand these expressions for small values of λ and obtain the weak-coupling expansions. Actually, proceeding in this way, we are able to push the perturbative calculations to very high orders with relatively low computational effort and generate very long expansions that can then be used for numerical investigations. Most importantly, we can exploit the exact knowledge of the λ-dependence of the X matrix and study the above deviations at strong coupling in an analytic way. For λ → ∞ we find that they behave as    √ λ for the untwisted observables , λ −1 for the twisted observables . (1.6) The precise numerical coefficients in front of λ −1 for the twisted observables we have analyzed can be fixed in full generality, while the coefficients in front of √ λ for the untwisted observables are more delicate and depend on the precise strong-coupling behavior of the free energy. However, once this is fixed, also these coefficients are fixed.
To substantiate our findings, we have compared the analytical strong-coupling predictions with numerical results obtained by using Padé approximants of the long perturbative series and by an independent Monte Carlo simulation. The agreement we have found is remarkable, especially in the twisted case.
Since the E theory has a holographic dual corresponding to a geometry of the type AdS 5 × S 5 /Z 2 , it would very interesting to retrieve the above results from a gravitational point of view. In particular, our results for the twisted observables provide a precise prediction for the twisted sectors of the AdS 5 ×S 5 /Z 2 geometry, which so far have remained largely unexplored. Furthermore, it would be interesting to go beyond the leading order in the strong-coupling expansion and study the properties of the resulting series for large λ, and also to investigate the role played by nonperturbative instanton corrections and/or resurgence effects, along the lines recently discussed in [46,47] for the N = 4 SYM theory.

Outline of the paper
In Section 2 we review the main properties of the E theory and of its associated matrix model arising from localization. In Section 3 we derive the exact expressions of the leading corrections in the large-N expansion for twisted and untwisted observables of the E theory in terms of the infinite matrix X. In Section 4 we work out the analytic behavior at strong coupling of these observables at leading order, which we check against numerical analyses and Monte Carlo simulations in Section 5. Finally, a few more technical details concerning the infinite matrix X and its properties are collected in Appendix A.

N = conformal SYM theories and matrix model
In this section we review some key properties of the N = 2 conformal theories with the aim of establishing the notation and describing the set-up in which we will perform our calculations. Most of this material has already been discussed in the literature (see in particular [22,25]) and thus we can be rather brief.

N = 2 conformal SYM theories
We consider N = 2 SYM theories in 4d with gauge group SU(N ) whose matter content consists of N F hypermultiplets in the fundamental representation, N S hypermultiplets in the symmetric representation and N A hypermultiplets in the anti-symmetric representation. These theories are conformal provided the 1-loop β-function coefficient vanishes. Among the five families of models that satisfy this condition (see [2,48,49] and also [13,20,22,25]), here we focus on the following two families, parametrized by N and characterized by D : Even though these theories have received less attention than SQCD, which has N F = 2N , N S = N A = 0, they are very interesting for various reasons. Most notably, they share many features with N = 4 SYM in the large-N limit, and admit a holographic dual description in terms of strings propagating in a geometry of the form AdS 5 × S 5 /Z 2 , where Z 2 is a suitable orientifold group (see for example [44]). For both D and E families, the ratio of the number of fundamental flavors with respect to (twice) the number of colors scales to zero in the planar limit, namely while in SQCD one has ν = 1. The vanishing of ν is a necessary condition for planar equivalence to N = 4 SYM, at least in a certain "even" sector [22,39]. Another interesting feature is that the difference between the c and a central charges is such that 1 24 for E .
On the other hand, the N = 4 SU(N ) SYM theory has ν = 0, since there is no fundamental matter, and α 0 = α 1 = α 2 = 0, since the two central charges are equal, c = a. We therefore see that in the large-N limit the D and E families share many features with the N = 4 SYM theory, and in particular the E theory is the one which goes closer. In this sense the E theory can be considered as the next-to-simplest SYM theory after the maximally supersymmetric one. From now on we will discuss only the E theory, even if most of our considerations can be applied with minor modifications to the D theory as well.
In the N = 2 SYM theories, an important set of local operators is provided by the multitraces where n = {n i } and ϕ is the (complex) scalar field in the adjoint vector multiplet. These operators are gauge-invariant, chiral, possess an R-charge equal to n = i n i and are automatically normalordered because of R-charge conservation. The anti-chiral operators O n (x) are constructed in the same way with the conjugate field ϕ(x). In the conformal case, the two-point functions between chiral and anti-chiral operators take the general form where G n,m is a non-trivial function of the gauge coupling g and of the number of colors N . To keep the presentation as simple as possible, in the following we restrict our attention to the single trace operators in which case the coefficients in the two-point functions (2.6) become diagonal and are just denoted by G n . These coefficients can be studied in the limit N → ∞ with the 't Hooft coupling λ = g 2 N held fixed, and compared with the corresponding coefficients in the N = 4 SYM theory, denoted by G n . For the E theory one finds where ∆ k and δ k are non-trivial functions of the 't Hooft coupling only. This means that when N → ∞ the two-point functions of even operators are identical to those in the N = 4 theory and differ only in sub-leading terms at order 1/N 2 , while the two-point functions of odd operators are different from the N = 4 ones already in the planar limit 2 . This structure is in full agreement with the expectations based on the AdS/CFT correspondence. Indeed, as discussed in the Introduction (see (1.1) and (1.2)), when the holographic dictionary is applied to the AdS 5 × S 5 /Z 2 orientifold that is dual to the E theory, the even chiral operators O 2k map to supergravity states of AdS 5 × S 5 which survive the projection. In other words, these excitations belong to the untwisted sectors of the orientifold and thus in the planar limit their correlators must be the same as in the maximally supersymmetric case. On the contrary, the AdS 5 × S 5 states that correspond to the odd chiral operators in the N = 4 theory are projected out in the orientifold construction, and the odd chiral operators O 2k+1 of the N = 2 theory are obtained from states of the twisted sectors of AdS 5 × S 5 /Z 2 . These supergravity excitations have no analogue in AdS 5 × S 5 , and thus it has to be expected that even in the planar limit their correlators be different from those in the N = 4 theory.
Our aim is to find an explicit expression of the functions ∆ k (λ) and δ k (λ) appearing in (2.8). One possibility to reach this goal is to use the standard perturbative field-theory methods and compute the Feynman diagrams that contribute to the two-point functions. This method is conceptually straightforward, but in practice it quickly becomes rather daunting as the number of loops increases. Even if some simplifications may occur in the large-N limit, still the number of diagrams that one has to compute is tremendously high (see [22,25] for some explicit examples at the lowest perturbative orders). Thus only the very first few perturbative terms can be explicitly computed in this way.
A much better approach is to exploit localization and map the superconformal gauge theory to an equivalent matrix model [1].

Matrix model
If one defines the N = 2 SYM theory on a sphere S 4 of unit radius, its partition function can be expressed in terms of a matrix model as follows [1] 3 : (2.9) Here we have adopted the so-called "full Lie algebra approach" [18,22,25] in which the integration is performed over all components of a Hermitean traceless matrix a belonging to the Lie algebra 4 of SU(N ). The factor Z 1−loop arises from a 1-loop computation, while Z inst is the non-perturbative instanton contribution. In the large-N 't Hooft limit, instantons are suppressed and thus we can put Z inst = 1. The 1-loop term depends on the representation R of the matter multiplets in the gauge theory and takes the following form Here with G being the Barnes G-function. In the N = 4 SYM theory, where R is the adjoint representation, S int clearly vanishes, while in the N = 2 SYM theories it accounts for the matter content of 3 As compared to [1], we have rescaled the matrix a with a factor of g 2 8π 2 in order to obtain a canonical normalization in the quadratic term, and have neglected the overall factor g 2 since it drops out in all calculations. 4 We write a = a b t b (b = 1, . . . , N 2 − 1) where t b are the generators of SU(N ) normalized in such a way that tr t b tc = 1 2 δ bc the so-called "difference theory" [3], in which one removes from the N = 4 result the contributions of the adjoint hypermultiplets and replaces them with the contributions of hypermultiplets in the representation R. In the perturbative regime when g → 0 we can use the expansion where γ E is the Euler-Mascheroni constant, and check that for the E theory S int is given by [22,25] S int (a) = 4 ∞ ,m=1 Therefore, the partition function (2.9) can be written as and S int (a) can be interpreted as an interaction action. Note that even for the E theory, which is the simplest N = 2 conformal gauge theory as we mentioned before, the corresponding matrix model contains an infinite number of interactions as is clear from (2.14). This is to be contrasted with the matrix model corresponding to the N = 4 SYM theory which is purely Gaussian.
Given any function f (a), its vacuum expectation value in the N = 2 matrix model is where in the second step we have denoted by 0 the expectation value in the free Gaussian model. In particular, we are interested in the matrix operators O n (a) that correspond to the chiral operators (2.7) of the gauge theory and in their correlators. Indeed, as shown in [9][10][11][12][14][15][16][17][18], the coefficients G n in the two-point functions of the gauge theory in flat space are given by (2.17) In this approach everything is then reduced to the calculation of vacuum expectation values in the free matrix model using (2.16). As discussed in [18,22,25] these vacuum expectation values can be efficiently computed using recursion relations, and explicit expressions can be easily generated even at very high perturbative orders.

Leading corrections in the 't Hooft limit
Our aim is to provide exact formal expressions of the leading corrections in the large-N limit for both twisted and untwisted operators in the E matrix model and for the corresponding observables in the gauge theory. In this section we show that these formal expressions can be encoded in an auxiliary theory of infinite free variables. In particular, we prove that the twisted corrections involve the propagator of this free theory, while the untwisted ones are expressed as derivatives with respect to λ of its free energy. From such exact expressions it is quite straightforward to expand in powers of λ and obtain the perturbative results at very high orders. From these long perturbative series it is then possible to compute Padé approximants which remain stable and accurate well beyond the radius of convergence of the perturbative series, thus allowing for an extrapolation towards strong coupling. In the next section we will describe how to extract the leading behavior of the corrections at large λ in an analytic way. To test the correctness of the results, in Section 5 we will make a comparison with the direct numerical evaluation of the vacuum expectation values in the matrix model by means of Monte Carlo simulations extrapolated to the large-N limit.

Gaussian representation of the E model at large N
Let us first recapitulate a few results of [22]. Consider the expectation value of a product of an even number of odd traces in the Gaussian theory: At the leading large-N order we have where H(k 1 , k 2 , k 3 , k 4 , . . .) represents the total Wick contraction computed with the "propagator" For instance, with 4 components one has Similarly, with 6 components one has the sum of the 15 possible ways to make a complete contraction, and so on. In order to exploit this Wick-like factorization property in the simplest way, it is convenient to rephrase it in terms of normal-ordered and normalized operators.
Normal-ordered operators: At the leading large-N order the normal ordered version O (0) n (a) of the single-trace operator tr a n contains only single traces of lower powers of a and is given by where the monic polynomial p n (a) is related the Chebyshev polynomials of the first kind T n (x) as follows [16]: For the first odd operators this amounts to while for the first even operators it gives These normal-ordered operators have diagonal two-point functions: The free-variable representation: Focusing on the odd sector, we introduce a normalized version of the normal-ordered operators, defining At large N , these objects have canonical two-point functions: The linear relation between the ω k (a)'s and the odd single traces ensures that, at the leading order for large N , the correlators of many ω's are computed using Wick's theorem just as it was the case for the correlators of many odd single traces (see (3.2)-(3.4)). Therefore, we can associate to each operator ω k (a) a real variable ω k and write where ω is the (infinite) column vector of components ω k and In this representation it is quite easy to include the interactions of the E matrix model, which are an infinite sum of double odd traces. By re-expressing these traces in terms of the operators ω k (a), the interaction action (2.14) takes the following form where X is an infinite symmetric matrix depending on λ originally introduced in [22]. The entries of X are initially expressed as series in λ since this is how they are generated but, remarkably, these series can be resummed in terms of an integral of Bessel functions of the first kind, according to [22] This expression can be easily re-expanded in series of λ at very high order, thus making it possible to obtain very long perturbative series in an efficient way. Moreover, as we will discuss in Section 4, the entries X k can also be expanded for large λ using Mellin transform techniques. This is the key to access the strong coupling behavior of all observables that can be given in terms of X.
A prime example of such observables is provided by the partition function (2.15) which can be rewritten as Since the Wick property exponentiates, so does the correspondence (3.12). Therefore, we can re-express Z in the free-variable representation as The corresponding free energy F is then given by Note that this F is actually the free energy in the "difference theory", namely

Twisted observables
The vacuum expectation value of any operator f (ω(a)) that can be written in terms of ω k (a) or, equivalently, purely in terms of the odd traces of a, can be realized using the free variables as follows From this expression, we see that the operators ω k (a), which were diagonal in the Gaussian theory, in the E model have the following propagator 5 This propagator D is an infinite matrix with a highly non-trivial dependence on λ. Multiple correlators ω k 1 (a) ω k 2 (a) · · · can just be treated as tree-level expressions and computed by applying Wick's theorem using D as propagator.
From these results it easy to deduce the expression of the two-point correlator of any two odd traces. In fact, it is sufficient to express these odd traces in terms of the ω k (a)'s by inverting the relations described in (3.5), (3.6) and (3.10). Explicitly, with the help of Eq. (5.6) of [22], one gets tr a 2k+1 = N 2 where the coefficients c k,i are Then, the correlator of any two odd traces in the E model is given by A few explicit examples are: Since each propagator D k, is known exactly in terms of λ through the matrix X, the above expressions represent a complete resummation of the perturbative results that have been previously considered in the literature. In turn, each D k, can be re-expanded in powers of λ and used to generate in an efficient way the perturbative series. This expansion can be done in two steps: firstly, by expanding in powers of X, then by expanding the matrix X using its definition (3.15). Thus, we first write D k, = δ k, + X k, + X 2 k, + X 3 k, + . . . , (3.25) and then we use the sum rule to obtain and so on, where we have defined for convenience . (3.28) These multiple integrals can be easily computed using the expansions of the Bessel functions.
In this way one finds that X n k, contains n powers of Riemann ζ-values with odd arguments and corresponds to the resummation of all perturbative contributions that in the gauge theory arise from diagrams involving n loops of matter hypermultiplets in the planar limit.
Normal-ordering in the E theory: While the above exact results for correlators in the matrix model are interesting by themselves, our main interest is in their bearing on the N = 2 SYM field theory of type E in flat space. The first step we have to carry out to reach this goal, is the Gram-Schmidt procedure to obtain the normal-ordered operators ω k (a) (3.29) that represent, at the leading large-N order, the chiral operators tr ϕ 2k+1 (x) of the N = 2 SYM theory, rescaled so as to have unit correlators in the Gaussian model according to (3.11). In other words, the ω k (a)'s are related to the operators O 2k+1 (a) introduced in Section 2 by By definition, these normal-ordered operators have diagonal two-point functions, which we write as where ∆ k (λ) vanishes for λ → 0. These operators capture precisely the twisted observables of the gauge theory defined in (2.8a), namely The key observation is that the Gram-Schmidt procedure expresses the two-point function of the normal ordered operators ω k (a) as a ratio of determinants of the correlation matrix (3.20), so that . (3.33) Here we have introduced the notation M (k) to indicate the upper left k × k block of a matrix M , with the convention that M (0) = 1. As shown in [22], the ratio of determinants in (3.33) can also be rewritten in a different and more suggestive way. To do so let us define M [k] as the sub-matrix which is obtained from M by removing its first (k − 1) rows and columns; in this way M [1] = M . Then, using the definition (3.20), one can prove that for any k We stress once again that this result is exact in λ. At weak coupling, it is quite straightforward to extract the perturbative expansion from the above expression. Indeed, one first expands (3.34) in powers of X [k] to get where X k,m X m,n X n,k , (3.36c) and so on, and then uses the sum rule (3.26) to reduce the calculation to multiple integrals of Bessel functions. This procedure can be easily automatized and the perturbative expansion of ∆ k (λ) can be obtained to very high orders with relatively short computational effort. For instance, at the first few orders for k = 1 and k = 2, one gets We have actually worked out these expansions up to order λ 160 , and constructed from them Padé approximants that turn out to be extremely stable far beyond the convergence radius of the perturbative expansion, as we will see in Section 5.

Untwisted observables
Let us now turn to the "untwisted" observables and consider an even-trace operator tr a 2k in the matrix model. Writing its expectation value as we realize that, at the first order in the interaction action, the deviation A 2k from the Gaussian result is given by As shown in Section 2, the interaction action S int (a) is a collection of double odd traces and can be recast in the form where the coefficients f p,q can be easily read from (2.14). Then, using the same notation as in (3.1), the correction (3.40) reads (3.42) As proved in [22], the even traces factorize in every Gaussian correlator at the planar level, meaning that t 2k,n,m,... = t 2k t n,m,... (3.43) at the planar level. Therefore, the first contribution to the connected correlator appearing in A 2k occurs at the sub-leading order and takes the form where the ellipses denote terms suppressed at large N . This contribution arises when just two propagators connect the even trace to the odd ones. Inserting this result into (3.42), we get 45) so that, at the lowest order in S int (a), we can rewrite (3.39) as follows The same pattern persists also when we consider higher powers of S int (a) and at order 1/N 2 we get contributions from connected diagrams where the even trace is linked to the odd ones. Thus we can express the result as tr a 2k tr a 2k where F is the free energy (3.18). If we consider multiple even traces, the contribution of order 1/N 2 to their vacuum expectation value arises when just one of them is linked to the odd traces coming from the interaction action. Thus, we simply have to sum over all such possibilities, getting tr a 2k 1 tr a 2k 2 · · · tr a 2k 1 tr a 2k 2 · · · 0 = 1 − l k l (k l + 1) Since the free energy F is known in terms of the matrix X as indicated in (3.18), also all these vacuum expectation values can be expressed in terms of X and thus the full λ-dependence of their first non-planar corrections is completely determined by (3.48). A similar reasoning can be used to find the correlators of the normal-ordered operators O 2k (a) that represent the single-trace observables of even weight in the gauge theory. To reach this goal, we have first to determine how O 2k (a) is written in terms of the even traces by means of the Gram-Schmidt diagonalization method. Once this step is completed, we can compute the correlators and see how they behave. This procedure is completely straightforward but rather tedious. Therefore, here we discuss only the simplest case corresponding to k = 1. In this case the normal-ordered operator associated to the chiral operator tr ϕ(x) 2 is just O 2 (a) = tr a 2 − tr a 2 , (3.49) and thus its two-point correlator is O 2 (a) O 2 (a) = tr a 2 tr a 2 − tr a 2 2 = tr a 2 tr a 2 0 − tr a 2 2 0 − A 2,2 . (3.50) At the first order in the interaction action, the deviation from the Gaussian result is A 2,2 = tr a 2 tr a 2 S int (a) 0 − tr a 2 tr a 2 0 S int (a) 0 − 2 tr a 2 0 tr a 2 S int (a) 0 + 2 tr a 2 2 0 S int (a) 0 . (3.51) Using the form of the interaction action given in (3.41), we can rewrite the previous expression as (3.52) The quantity in square brackets can be evaluated by exploiting the relations which were derived in [18] using the fusion/fission identities of the traces of SU(N ). In this way, after straightforward algebra, we find where the ellipses stand for terms suppressed at large N . This structure persists at higher orders in S int (a) and, as before, it exponentiates allowing us to express the final result in terms of the free energy F in the following way Mapping this to the chiral/anti-chiral correlator of the gauge theory, we find precisely the form anticipated in (2.8b), namely tr ϕ 2 (x) tr ϕ 2 (0) tr ϕ 2 (x) tr ϕ 2 (0) 0 We have also analyzed in this way the cases k = 2 and k = 3, and found leading us to conjecture that These expressions contain the exact dependence on λ since the free energy F can be written in closed form in terms of the matrix X as we have seen in (3.18). If one is interested in the perturbative expansion for λ → 0, this can be efficiently generated using These expansions can be easily pushed to very high orders without too much computational effort and then used for numerical investigations (see also [8] for similar calculations related to the vacuum expectation value of the BPS Wilson loop).

(4.3)
The asymptotic expansion of this expression for large λ receives contributions from the poles on the negative real axis and reads
It is quite remarkable that, as observed in [8], this is essentially the same matrix 6 appearing in the analysis of the zeros of Bessel functions performed in [50]. Note that the main diagonal and the super-and sub-diagonals of S converge to 0 when the size of the matrix increases. This is the key property that allows us to obtain finite results from this infinite matrix.

Twisted observables
In Section 3.2 we managed to express various twisted observables in terms of the propagator D of the free-variable representation. Using (4.5), we see that at strong coupling so that for large λ the twisted observables scale as 1/λ at the leading order.
In Appendix A, we show that it is possible to find an explicit form for the inverse of the infinite matrix S. In particular, from (A.38) we obtain Most importantly, also the quantities ∆ k (λ) of the E gauge theory admit an exact expression in terms of the matrix X, as we have shown in (3.33) and (3.34), which we rewrite here for convenience: For any given k we can use the expression in the middle, which is a ratio of determinants of finite matrices whose leading large-λ behavior can be deduced from (4.8). For example, for the lowest values of k we have Inserting this into the first equality of (4.10), it immediately follows that The second equality in (4.10) allows to to obtain the result in closed form for generic k. In fact, at large λ we can use (4.5) so that . (4.13) Note that in the last expression the ratio of the determinants of two infinite matrices appears. However, this ratio can be explicitly computed in closed form for any k, as we show in appendix A. In the end, using (A.23), we obtain (4.14) Therefore, from this analysis we find that the strong-coupling behavior of the two-point functions of chiral/anti-chiral operators in the E superconformal gauge theory is given by where the ellipses stand for sub-leading terms for large λ and large N .

Untwisted observables
In Section 3.3 we showed that, in the large-N limit, the correlators of even traces have a simple expression in terms of the free energy F given in (3.48). The strong coupling behavior of the free energy has been recently studied in [8], where it has been found that The scaling F ∼ √ λ is strongly supported by a Padé analysis of a long perturbative series for F obtained by exploiting the matrix X. Replacing the latter in (3.18) by its "leading order" (LO) approximation (4.5), one can analytically prove that µ LO = 1/(2π). A further discussion of the overall normalization µ will be presented in Section 5.2.2. Assuming the behavior in (4.16), we can use the exact formula (3.48) to express the correlators of multiple even traces at strong coupling as follows: tr a 2k 1 tr a 2k 2 · · · tr a 2k 1 tr a 2k 2 · · · 0 ∼ λ→∞ 1 − µ l k l (k l + 1) Likewise, for the correlators of the normal-ordered operators O 2k (a) that represent the chiral and anti-chiral operators of the gauge theory, we can use the exact formulas (3.58)-(3.61) and conclude that This expression has been checked for k = 1, 2, 3 and is conjectured for higher values of k.

Numerical analysis
In this section we discuss two independent numerical methods to evaluate the vacuum expectation value of a given observable in the E theory in the 't Hooft limit at strong coupling, i.e. for λ 1.
The first method relies on the fact that in the weak coupling regime these vacuum expectation values admit a perturbative expansion in λ that can be computed starting from their representation based on the matrix X defined in (3.15). These expansions may be pushed to very high order in λ and their numerical coefficients can be evaluated with high precision. Therefore, analytic continuation beyond the convergence radius can be done by standard Padé extrapolation [51,52].
A different approach is based on the direct numerical integration of the matrix model provided by supersymmetric localization (see the recent analysis in [37]). The vacuum expectation value of a given observable is computed by a finite N -dimensional integral with non-negative integrand. In the large-N limit, the relevant integration domain where the integrand is sizable shrinks around the saddle point configuration, i.e. around a zero-measure set. This forbids direct numerical quadrature methods. Nevertheless, this is a standard problem in statistical mechanics and lattice field theory and its solution requires dedicated Monte Carlo (MC) methods that are well established (see for example [53]).
As we will discuss in the following, the application of these two methods leads to strong coupling results in agreement with the analytical predictions previously discussed.

Padé approximants for ∆ 1 and ∆ 2
As we have seen in Section 3, it is possible to compute the quantity ∆ k defined in (3.31) to very high orders in perturbation theory, finding (5.1) For example, for ∆ 1 and ∆ 2 the first few coefficients of the above expansions are given in (3.37) and (3.38). As it was already shown in [22], the radius of convergence of the series for ∆ 1 is located at λ = π 2 . This is due to a branch point at λ = −π 2 . To reveal possible higher singularities it is convenient to make a conformal map before computing Padé approximants [54][55][56]. In particular, as in [8], we map the expansions (5.1) into the unit disk |z| ≤ 1 by the transformation More precisely, we make the replacement λ π 2 → 4z (1−z) 2 in the perturbative series expansion, then we construct the Padé approximants and finally we replace z → λ according to (5.2).
Since we expect ∆ k to be asymptotically constant with 1/λ corrections, a suitable choice for the Padé approximants is to take the diagonal ones. Hence, we shall extract information on the non-perturbative region λ > π 2 , from or by improving convergence using the above conformal map. Notice that M , the degree of the polynomials, must satisfy M < L/2. In the following we focus on the cases k = 1, 2 corresponding to ∆ 1 and ∆ 2 , respectively. In the z-plane, the poles of the diagonal Padé approximants are shown in Fig. 1 for ∆ 1 (left) and ∆ 2 (right). They tend to accumulate at z = −1 corresponding to the branch point λ = −π 2 . We can observe also two symmetric sub-sequences of zeroes that tend to a point on the unit circle corresponding to λ = −4 π 2 . The same analysis on longer perturbative expansions is expected to reveal collinear branch points at all λ = −n 2 π 2 with n = 1, 2, . . . .

Monte Carlo simulation
There are several options for the choice of a MC algorithm suitable to our problem. For our purposes, and with the level of precision we will need in the analysis, it is possible to use one of the simplest algorithms, i.e. the Metropolis-Hastings (see for instance [57]). This algorithm constructs a Markov chain of configurations {Y j } ∞ j=0 obeying a detailed balance and each state of this chain corresponds to a particular configuration Y of the eigenvalues {a i } of the localization matrix model. Moreover, each configuration Y j is weighted by a factor e −S(Y j ) , where S(Y j ) = S classical (Y j ) + S int (Y j ). The transition dynamics of the chain is built performing in an iterative way the two following steps: 1. Starting from the element Y j we generate the next element of the chain Y j+1 . This step is performed selecting in a random way one of the eigenvalues of the starting configuration Y j , let's denote it by a j . Then we randomly change the value of a j by a step of chosen size namely where u is a uniform random number in [− , ]. In this way we obtain a large set of configurations {Y j } with j = 1, · · · , n that are asymptotically distributed according to e −S in the limit n → ∞. The statistical estimator of the vacuum expectation value of a generic observable O(Y j ) is simply the arithmetic average over the elements of the chain, namely Evaluating (5.5) at finite n gives an estimate which has a non-negligible statistical error, but is exact for infinite statistics. For finite sampling sizes, the statistical error is determined from the variance of the estimator, corrected by the auto-correlation length of the measurements. Both these sources of errors have been measured and taken into account. As a common rule of thumb, we tune the size of the step in order to have an acceptance rate around 50-60%.

Twisted operators
We first consider the twisted observables ∆ 1 and ∆ 2 for which the analytic prediction reads (see (4.12)): Numerical results for these two quantities are reported in Fig. 2, where we compare the predictions (5.6) with the diagonal Padé approximants computed with M = 60 for ∆ 1 and M = 70 for ∆ 2 .
The agreement at large λ is excellent.
As we discussed, the conformal map (5.2) is effective in separating out different collinear branch points lying on the negative real axis. Actually, it is also expected to speed up convergence with a certain fixed number of terms in the series (5.1). This is a general feature of a broad class of divergent series, (see for example [55]), and remains valid here where the perturbative series has a finite radius of convergence with a branch point on the boundary of the convergence disk. 7 To appreciate the improvement associated with the conformal map (5.2), we also show in Fig. 3 the products λ π 2 γ 3 and λ π 2 γ 5 using this further refinement. In the first part of these curves we see how the Padé approximant tends to the right value according to (5.6), but does not reach the expected plateau at very large λ where it overshoots and eventually breaks down. A much better convergence is observed for the conformally improved Padé approximants that come very close to the analytical prediction before becoming unreliable. Figure 3: Comparison between the simple and the conformally improved Padé approximants for (λ/π 2 )γ 3 (left) and (λ/π 2 )γ 5 (right). This quantity is expected to tend to 24 for γ 3 and to 80 for γ 5 .
The same quantities have been analyzed by MC simulations at various pairs (N , λ). In principle, one should extrapolate the finite N results at fixed λ in the limit N → ∞. In practice, showing simultaneously our data at increasing values of N clearly shows convergence. Indeed, the data points very nicely tend to the Padé prediction as N increases. In the MC simulations we cannot reach very large values of λ and therefore in the comparison we use the simple Padé approximants which are equivalent to the conformally improved ones in this range of coupling. This analysis is shown in Figs. 4 and 5 for the quantities γ 3 = 1 + ∆ 1 (λ) and γ 5 = 1 + ∆ 2 (λ), respectively.

Untwisted operators
Let's now consider an application of these numerical methods to the untwisted operators. Although both the Padé approximants and the MC simulations can be applied to a vast set of untwisted  Figure 5: Padé curve (red curve), MC data and the large-λ theoretical prediction for γ 5 = 1 + ∆ 2 (λ) (black dashed curve) in the range 1000π 2 ≤ λ ≤ 4000π 2 . As N increases, the MC points systematically tend towards the Padé curve.
observables, for definiteness here we focus just on a particular example, namely the ratio tr a 4 tr a 4 0 = 1 + q 4 (λ) N 2 + . . . , (5.7) where, as can be seen from (3.47), the function q 4 (λ) is given in terms of the free energy F according to Therefore, in order to make an accurate prediction it is crucial to have a precise evaluation of F at strong coupling. This is far from trivial as we see in the following paragraph.
Strong coupling limit of F A first evaluation of the free energy has been recently performed in [8] by considering the LO expression of the matrix X at large λ given in (4.4). This gives the following result: As explained in [8], the comparison between (5.9) and a conformally improved Padé analysis is not fully satisfactory. In fact, the resummation clearly shows a ∼ √ λ behavior, but the overall coefficient is sizeably smaller then the prediction in (5.9). The apparent mismatch may be due to logarithmic corrections that affect convergence to the strong coupling regime but are hard to be visible in the explored range of λ accessible by numerical simulations. Including the next-toleading-order (NLO) term in the matrix elements of X is problematic. Indeed, the NLO correction requires to evaluate the ill-defined quantity where S is the three-diagonal matrix given in (4.6). One possibility is to evaluate (5.10) numerically by the empirical procedure described in Section 6.2. of [8]. In practice this means defining where, in the notation introduced after (3.33), S (k) denotes the upper left k ×k block of S, and then extrapolating the result for k → ∞. In this limit this provides an estimate for F NLO . To obtain this extrapolation, we first fix k and look for the value of λ that maximizes the ratio F NLO k / √ λ. We denote this maximum as Then, by successively increasing k, we generate a set of values of µ NLO k , which we fit as follows The value µ NLO = 0.113 , (5.14) provides an empirical estimate of the coefficient in front of √ λ at the NLO. This result is substantially smaller than the LO value 1/2π given in (5.9) by a factor of about √ 2. A further refinement consists in avoiding the use of the asymptotic expansion of the matrix elements X k and in using instead their exact expression as a definite integral given in (3.15). This is non-trivial because of the highly oscillatory behavior of the integrand associated with the Bessel functions. We can deal with it by splitting the integration region in sub-intervals between two zeroes of J 2k+1 J 2 +1 . This gives an alternating series for which it is possible to estimate the convergence error. Repeating the analysis in this way, we have obtained where "full" emphasizes that we used the exact (numerical) values of X k (λ). The extrapolation in (5.15) is done with smaller values of k than in (5.14), but the shown three digits appear to be stable.
In Fig. 6 we show the Padé resummation of F/ √ λ and its conformal improvement, taken from [8], together with the three estimates for the asymptotic value of this ratio, i.e. the LO prediction 1/2π from (5.9), the NLO value in (5.14), and the "full" extrapolation in (5.15). In the adopted logarithmic scale, we can reveal possible logarithmic contributions that may substantially slow down the convergence to the plateau. The conformally improved purple curve overshoots the NLO prediction, while it remains below µ full almost reaching it at λ ∼ 10 7 . In the range λ < 10 3 the curve is clearly not flat, but totally under control, given the agreement between the simple and the conformally improved Padé approximants.
Analysis of q 4 (λ) Let's now return to function q 4 (λ) appearing in (5.7). Using (5.9) for the strong-coupling behavior of the free energy, we obtain If instead we use the estimates (5.14) and (5.15), we have

(5.17)
We have compared these predictions with two independent numerical evaluations, namely a diagonal Padé approximant 8 and a MC simulation. Our results are collected in Fig. 7 where we observe equality between the Padé approximant and the MC data within error bars. Given the considered large values of λ, this is an important independent check of the proposed resummation based on the Padé representation. However, as it is clear from Fig. 6, we do not expect the expressions Figure 6: Comparison between various estimations of F/ √ λ. The highest horizontal line is the expected asymptotic constant using the LO analytical prediction. The other horizontal lines are the prediction from the truncated matrix analysis using the NLO expansion of the M matrix elements or their "full" integral representation. Also shown are the simple and conformally improved Padé resummations taken from [8]. Notice that a residual steady increase of the purple curve at quite large values of λ is clearly visible in the adopted logarithmic scale. . (5.16) and (5.17) to be a good approximation in the range explored in Fig. 7. Indeed, the asymptotic regime becomes visible at much larger values λ ∼ 10 6 , presumably due to slow logarithmic corrections. Nevertheless, it may be interesting to look at the behavior of the three analytical expressions in (5.16) and (5.17). We observe that the LO prediction (5.16) is systematically much lower than the MC data and the Padé approximant. The NLO curve starts to overshoot the Padé approximant at λ ∼ 350π 2 , while the "full" curve closely follows the Padé resummation keeping underneath it, but with similar slope.
We remark that the untwisted MC measurements that we could achieve still have sizeable relative errors around 10%. Besides, the data at higher values of λ are less reliable due to systematic errors in the large-N extrapolation needed to extract q 4 (λ) from the slope in (5.7). This fact is an indication that the untwisted observables, as the one we have considered, are much more difficult to be evaluated via a MC simulation than the twisted ones. In order to reduce the relative errors and increase the precision, a more intense computational effort would be necessary. However this is far beyond the scope of the present work. As a final comment, we stress the great flexibility of the MC approach. In this framework, in fact, more complicated observables may be studied at quite large values of gauge coupling with minor effort, even without resorting to the special role played by the X matrix in the considered model. (λ), q full 4 (λ) defined in (5.16) and (5.17). We observe the agreement between the Padé approximant P [70/70] (q 4 ) and the MC points (within error bars). Since, as shown in Fig. 6, the asymptotic regime becomes manifest only for larger values, i.e. λ ∼ 10 6 , we do not expect the expressions q LO 4 (λ), q NLO 4 (λ), q full 4 (λ) to be a good approximation in the range 0 ≤ λ/π 2 ≤ 550 analyzed here.
Pesando for helpful discussions. This research is partially supported by the INFN projects ST&FI "String Theory & Fundamental Interactions" and GSS "Gauge, Strings and Supergravity". The work of A.P. is supported by INFN with a"Borsa di studio post-doctoral per fisici teorici".

A Details on the derivation of the strong-coupling results
In Section 4 we have written in (4.5) the leading behavior for large λ of the infinite matrix X in terms of the infinite matrix S defined in (4.6). At strong coupling, the functions ∆ k (λ) that represent the deviation from the N = 4 result of the correlator of two twisted observables can then be expressed in terms of determinants of the sub-matrices of S, as shown in (4.13). Moreover, the expectation values of odd multi-traces are expressed in terms of the elements of the matrix D and thus, using (4.7), in terms of the elements of the inverse of S. In this appendix we give some details on how to treat these infinite matrices and compute the invariants connected with S. We first perform a similarity transformation and write and Y is not symmetric but has rational entries given by .
In the same way we can compute the determinants of the sub-matrices Y (k) corresponding to the upper (k × k) block of Y. Indeed, defining N (k) and Λ (k) in the same way, we have Y (k) = N (k) Λ (k) , (A.24) and det Y (k) = det N (k) det Λ (k) = 1 2 (k + 1)(k + 2) k =1 1 (2 + 1)(2 + 2) . (A.25) More generally, these techniques can be used to compute any element of the infinite matrix Y −1 . In fact, where the Y is the cofactor matrix whose entries are (A. 38) which has been used in Eq. (4.8) of the main text.