N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N} $$\end{document} = 2 Conformal SYM theories at large N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N} $$\end{document}

We consider a class of N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N} $$\end{document} = 2 conformal SU(N) SYM theories in four dimensions with matter in the fundamental, two-index symmetric and anti-symmetric representations, and study the corresponding matrix model provided by localization on a sphere S4, which also encodes information on flat-space observables involving chiral operators and circular BPS Wilson loops. We review and improve known techniques for studying the matrix model in the large-N limit, deriving explicit expressions in perturbation theory for these observables. We exploit both recursive methods in the so-called full Lie algebra approach and the more standard Cartan sub-algebra approach based on the eigenvalue distribution. The sub-class of conformal theories for which the number of fundamental hypermultiplets does not scale with N differs in the planar limit from the N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N} $$\end{document} = 4 SYM theory only in observables involving chiral operators of odd dimension. In this case we are able to derive compact expressions which allow to push the small ’t Hooft coupling expansion to very high orders. We argue that the perturbative series have a finite radius of convergence and extrapolate them numerically to intermediate couplings. This is preliminary to an analytic investigation of the strong coupling behavior, which would be very interesting given that for such theories holographic duals have been proposed.


Introduction
One ambitious but necessary goal in theoretical physics is to understand the dynamics of interacting Quantum Field Theories (QFTs) at strong coupling. Many ideas have been proposed and investigated, often involving the use of duality relations. Among them, a prominent role is played by the Anti-de Sitter/Conformal Field Theory (AdS/CFT) correspondence [1][2][3]. N = 4 Super Yang-Mills (SYM) theory represents a benchmark for exact computations in QFTs and for an explicit realization of the AdS/CFT correspondence. In a way, this theory is the simplest interacting four-dimensional QFT, since it enjoys the highest possible amount of superconformal symmetry [4] for a theory with at most spin-one fields, and Sduality. This large symmetry constrains so much of its dynamics that many sectors can be described in an exact way by resorting to powerful techniques, among which we can mention localization and the relation to integrable models. Furthermore, the N = 4 SU(N ) SYM theory admits an holographic description as type II B strings on AdS 5 × S 5 which is the prototype of all AdS/CFT relations. In this context, the 't Hooft large-N limit singles out the planar diagrams on the field theory side and suppresses string loop effects.
Important achievements have been obtained in this highly symmetric context also in presence of extended objects, like the BPS Wilson loops, that preserve part of the N = 4 superconformal symmetry [5][6][7][8] and are examples of conformal defects [9][10][11][12][13][14][15][16]. Many of these results can be efficiently derived using supersymmetric localization [17,18], which allows to reduce the calculation of the partition function on a sphere S 4 and of other observables to a computation in a Gaussian matrix model.
Many efforts have been devoted over the years to extend as much as possible these results to less symmetric theories. In this perspective, the N = 2 SYM theories play an outstanding role. For such theories the localization procedure is available [17]. It expresses a class of observables on S 4 -including chiral operators and BPS Wilson loops -in terms of a matrix model. This matrix model is no longer Gaussian as in the N = 4 case, and contains both perturbative and non-perturbative contributions. When the N = 2 theory is conformal, 1 from these localization results it is possible to obtain information also about the analogous observables in flat space 2 [22][23][24][25].
It is obviously of great importance to study N = 2 SYM theories in the large N limit and at strong coupling and to understand if and how some analogue of the AdS/CFT duality applies [26][27][28][29]. In this paper we provide some contributions to this long term goal by exploiting the localization matrix model to extract in a rather efficient way the expression of protected observables in the large-N limit. We do this for a certain class of superconformal theories with matter in the fundamental, symmetric and anti-symmetric representations. 1 Large classes of N = 2 superconformal theories were early found in [19]. 2 In the case of chiral/anti-chiral two-point functions, it has been argued in [20] that the matrix model reproduces to a large extent the flat-space result also in theories with a non-vanishing β-function. On the contrary, when conformal invariance is explicitly broken by mass terms as in the N = 2 * theories, the observables computed from the matrix model differ from those computed in flat space [21].

JHEP09(2020)116
For a sub-class of these theories a dual holographic description, built out as an appropriate orientifold projection of the AdS 5 × S 5 geometry, has been proposed in [26]. These theories are extremely close to the N = 4 SYM theory: many observables coincide at large N with the N = 4 results. However, observables involving chiral operators built with traces of odd powers of the scalar fields do not; in the holographic correspondence of [26] these odd-dimensional observables are related to twisted sectors of the orientifold. If we regard the N = 4 SYM theory as the simplest non-trivial QFT, these theories represent the next-to-simplest ones. One could hope to be able to explicitly study them beyond the perturbative regime, at least in the 't Hooft large-N limit, and to match the field-theoretic description with its holographic dual. In this paper, from the matrix model we obtain, through an effective description in terms of free variables, closed forms of the perturbative series for the odd observables which appear to have a finite radius of convergence. Although we are not yet able at this stage to infer analytically the strong coupling behavior, these expressions allow for reliable numerical extrapolation to the intermediate coupling regime and indicate that the strong coupling regime might not be out of reach in a nearby future.
Let us now be more specific about the content of this work which is divided in two parts. In the first part we review the matrix model methods for N = 2 conformal models at large values of the rank N of the gauge group SU(N ). We distinguish two matrix model approaches. In the original Pestun derivation [17] the matrix model was written as an integral over the Cartan sub-algebra variables, i.e. over the matrix eigenvalues. In this framework the large-N limit is described by the asymptotic eigenvalue distribution which satisfies an integral equation obtained with a saddle-point approximation [22-25, 28, 30-38]. For the N = 2 conformal theories we are considering, this integral equation depends on the matter content only through a single parameter ν, which counts the fraction of hypermultiplets transforming in the fundamental representation: . (1.1) The second approach has been developed more recently, and it has been named the "full Lie algebra approach" in [39]. It consists of keeping the matrix integral over the full Lie algebra and developing a series of recursive rules [40][41][42] to evaluate correlation functions. These techniques are very efficient in the perturbative regime at finite N . They allow to explore different sectors of the gauge theory [43][44][45][46][47][48] and can be used also in non-conformal cases [20]. Similar methods have been used also in [39,49]. We shall see that in fact also the large-N regime is easily accessible within the full Lie algebra approach by exploiting the recursion relations in a suitable way. We present the main technical points of the two approaches and present a thorough evaluation of the following observables: • the vacuum expectation value of the 1/2 BPS circular Wilson loop; • the two-point functions of chiral/anti-chiral operators; • the one-point function of chiral operators in presence of a Wilson loop.

JHEP09(2020)116
In doing so, we address also some specific issues related to the computation of correlators with operators of odd dimensions, which play a crucial role throughout the paper.
In the second part of the paper we concentrate on a set of N = 2 theories whose fundamental matter content does not scale with N , so that they have ν = 0. As noted above, these models have a holographic dual [26] and are very close to the N = 4 SYM theory, as confirmed by the fact that some observables, such as the vacuum expectation value of the Wilson loop [42] and the Bremsstrahlung function [23,50], do not deviate from the N = 4 result in the large-N limit. In the present paper we compute the set of observables listed above for the ν = 0 theories, using matrix model techniques, and clarify which observables are different with respect to N = 4 in the planar limit. It turns out that the correlation functions involving only chiral operators made of traces of even order have the same behavior as in the N = 4 SYM theory; this applies to both chiral/anti-chiral correlators and one-point functions in presence of the Wilson loop. Correlation functions with traces of odd order, instead, do deviate from the N = 4 results through an infinite perturbative series. This analysis allows to identify a sort of "twisted" sector of operators that, in the holographic correspondence, feel the presence of the orientifold, consistently with the analysis of [26].
We also investigate the difference between even and odd correlators with a perturbative expansion of the N = 2 field theory directly in flat space. Using the N = 1 super-space formalism and the diagrammatic difference between N = 2 and N = 4 [24,[40][41][42][51][52][53][54] we are able to perform an explicit large-N analysis of the two-point functions of operators with low dimensions up to three-loops in the ν = 0 models. This allows us to understand at the diagrammatic level the origin of the different behavior of the even and odd correlators in the planar limit, and to infer the general structure of the leading term in the two-point correlator of operators with arbitrary odd dimension.
Finally, we develop some new matrix model techniques which are particularly efficient when ν = 0 and make it very easy to obtain for any odd correlator expansions at any desired order in perturbation theory. By applying some numerical resummation technique to these long expansions, we produce a first attempt in going beyond perturbation theory. Even if this still falls short of an analytic treatment of the strong coupling regime, we think that our results for this special class of N = 2 theories may represent a first step towards this important goal.
Several technical details, which are useful to reproduce our results, and some explicit high-order expansions are collected in the appendices.

Part I
In the first part of this paper we study N = 2 conformal SYM theories and several of their physical observables in the planar limit using matrix model techniques.

JHEP09(2020)116
2 N = 2 CFT theories We consider N = 2 SYM theories in 4d with gauge group SU(N ). 3 The field content of these theories consists of one N = 2 vector multiplet, which contains the gauge vector A µ (x) and a complex scalar field ϕ(x) plus their fermionic partners, all transforming in the adjoint representation, and several N = 2 matter hypermultiplets, each containing two complex scalars plus their fermionic partners transforming in a representation R. Given this field content, the gauge coupling constant g receives contributions only at one loop, and the coefficient β 0 of the β-function is where i R is the index of R. In the following we will focus on conformal theories, for which β 0 vanishes. This condition is clearly satisfied if R is the adjoint representation, in which case we have a SYM theory with N = 4 supersymmetry. An important set of local operators in these theories is provided by the multitraces where n = {n i }. These operators are chiral, i.e. they are annihilated by half of the supercharges. Their R-charge is n = i n i and they are automatically normal-ordered because of R-charge conservation. The analogous anti-chiral operators, constructed with the conjugate field ϕ(x), are denoted by O n (x). We will study the two-point functions between chiral and anti-chiral operators, which in conformal theories take the general form In the following we will often consider diagonal cases, for which we will employ the streamlined notation G n ≡ G n,n . Another operator that we will consider is the half-BPS Wilson loop in the fundamental representation on a circle C of radius R: where P denotes the path-ordering. In particular, we will study its vacuum expectation value (v.e.v.) and its one-point functions with the chiral operators, whose form is fixed by conformal invariance to be [9,41] Table 1. The five families of N = 2 superconformal theories with SU(N ) gauge group and matter in fundamental, symmetric and anti-symmetric representations. These theories have been identified long ago in [55], and more recently they have been considered in [23,44].
The sets of functions G n,m (g, N ), w(g, N ) and w n (g, N ) are the main subject of our analysis. In particular we will study these functions in the large-N 't Hooft limit in which N → ∞ with λ = g 2 N held fixed. Notice that all these observables involve operators constructed entirely with fields of the gauge multiplet. Therefore, in their perturbative evaluation, the matter hypermultiplets run only inside the loops.

The ABCDE theories
To be specific, we will consider theories whose matter fields transform in the following representation of SU(N ): corresponding to N F fundamental, N S symmetric and N A anti-symmetric hypermultiplets. For these theories the β-function coefficient (2.1) reads The condition β 0 = 0 leads to five families of N = 2 superconformal theories, whose field content is displayed in table 1. Theory A is the N = 2 conformal QCD. Theories D and E are superconformal models for which a holographic dual of the form AdS 5 × S 5 /Γ has been identified [26]. In the last column of the table we have written the values of the variable ν defined in (1.1), which in the specific case becomes As mentioned in the Introduction, this quantity determines the large-N behavior of the theory, so that theories B and C become equivalent in the large-N limit, and the same is true for theories D and E.

Matrix model from localization
Exploiting localization, it is possible to prove that certain protected observables of the N = 2 SYM theories can be exactly reduced to a matrix model computation [17,18]. Among JHEP09(2020)116 these observables there are the partition function on a four-sphere S 4 , the correlators between chiral and anti-chiral operators, as well as the v.e.v. of a circular BPS Wilson loop and the one-point functions of chiral operators in presence of the loop.
The sphere partition function: the partition function of a N = 2 SYM theory with gauge group SU(N ) on a sphere S 4 of unit radius can be expressed as follows [17]: Here m u are the eigenvalues of a Hermitean traceless (N × N ) matrix M and ∆(m) is their Vandermonde determinant Moreover Z(im, g) is the partition function for the theory on R 4 with gauge coupling g evaluated at a point in the Coulomb moduli-space parametrized by the eigenvalues m u . It consists of a classical, a one-loop and an instanton factor: The classical part is simply The instanton part can be neglected when working in perturbation theory; moreover it does not contribute in the large-N 't Hooft limit, and thus in the following we set Z inst = 1. The one-loop part depends on the matter representation R. Denoting by m the Ndimensional vector of components m u , by W (R) the set of the weights w of the representation R and by W (adj) that of the adjoint representation, we have where with G being the Barnes G-function. Writing we deduce from (2.14) that the interacting action is [42] where we have introduced the combination of traces

JHEP09(2020)116
In the N = 4 theory, where R is the adjoint representation, this combination clearly vanishes while in the N = 2 theories it accounts for the matter content of the so-called "difference theory" which is often used in field theory computations [51], where 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.
For the class of theories listed in table 1, by combining the tree-level and one-loop factors, we obtain In terms of the eigenvalues m u , the matrix model action S is explicitly given by The BPS Wilson loop: in [17] it was shown that also the v.e.v. of the supersymmetric circular Wilson loop (2.4) can evaluated exactly using localization. The result, which was anticipated by direct diagrammatic computations in [6,7], is when the Wilson loop is in the fundamental representation. This means that in the matrix model W C is represented by the following operator For conformal theories the v.e.v. of W C on S 4 coincides with the flat-space observable w(g, N ) defined in (2.5).

The full Lie algebra approach
The sphere partition function (2.10) is expressed in terms of the eigenvalues m u of M . The integration measure over these eigenvalues with the Vandermonde determinant is precisely the one that arises from diagonalizing the flat measure over the entire matrix M ; moreover, the exponential weight, e −S , can be written in terms of the matrix M using (2.13) and (2.17). Thus, the partition function can be recast in the form namely as an integral over all elements of M . In the following we will use the conventions of [40,41] and rescale the matrix M so as to get a tree-level term in the matrix model action with unit weight. We then introduce the matrix

JHEP09(2020)116
understanding, from now on, that it is traceless. The partition function (2.23) reads where in the second step we used the notation for any function of a. This shows that Z S 4 can be regarded as the expectation value of e −S int (a) in the free Gaussian model. We now consider a basis of su(N ) generators T b , with b = 1, . . . , N 2 − 1, normalized as and write a = a b T b . Then, the flat integration measure appearing above becomes where the normalization has been chosen in such a way that (0) = 1 and the "propagator" for the components of a is simply Let us now discuss the interacting part of the matrix model in this full Lie algebra approach. From (2.17) and (2.24) we see that If g is small, we can use the expansion where γ E is the Euler-Mascheroni constant, and obtain Notice that the quadratic term, proportional to (1 + γ E ), drops out since it is proportional to the coefficient β 0 which vanishes for the theories we are considering; indeed Tr R a 2 = Tr R a 2 − Tr adj a 2 = −β 0 tr a 2 = 0 . (2.33)

JHEP09(2020)116
The higher traces Tr R a 2k with k > 1 appearing in the interacting action can be re-expressed in terms of traces in the fundamental representation as follows [42]: (2.34) Therefore, in the full Lie algebra approach the action S int of the interacting matrix model is a linear combination of single traces and double traces of powers of a in the fundamental representation with coefficients depending on the gauge coupling g and on ζ-values. Given this structure, the expectation value of a generic function f (a) in the N = 2 matrix model is defined by . (2.35) Evaluating this v.e.v. in perturbation theory is therefore just a matter of computing expectation values in the free Gaussian model using the propagator (2.29).
The BPS Wilson loop: in the full Lie algebra approach, the field theoretic BPS Wilson loop expectation value is exactly captured by the v.e.v. of the following operator in the matrix model which simply follows from (2.22) upon using the redefinition (2.24).
Chiral operators and normal ordering: as shown in [22, 25, 34-36, 40, 56-58] the matrix model encodes information also about the flat space correlators of chiral/anti-chiral operators in conformal theories. 4 To obtain this information one maps the multitrace operators O n (x) introduced in (2.2) to suitable matrix operators O n (a) such that their two-point functions correspond to the quantities G n,m defined in (2.3), namely where the v.e.v. in the right hand side is computed according to (2.35). Naïvely one would associate O n (x) to matrix operators with the same trace structure, that is Ω n (a) = tr a n 1 tr a n 2 . . . . are not diagonal. Therefore, as argued in [25,40], one has apply the Gram-Schmidt orthogonalization procedure with respect to C n,m , and construct the normal-ordered version of Ω n (a), which we denote by O n (a). If C (n) is the matrix of two-point functions among the operators of dimension lower than n, then one finds (2.40) With this definition, the correlator of O n (a) with any operator of lower dimensions is zero and the two-point functions G n,m vanish for n = m, as required. Using (2.40), these twopoint functions can be expressed in terms of the correlators C n,m . It turns out that they have a particularly simple expression when there is a single independent operator for each dimension n, as in the SU(2) case [25]. We will see that in the large-N limit also the set of single-trace operators is closed under normal ordering and we will restrict our attention to this set. Since there is one single-trace operator for each dimension n, then the matrix G has a simple expression in terms of the matrix C.
Recursion relations: from what we have reviewed above, it is clear that the basic ingredients for the calculation of the various observables in the N = 2 matrix model are the expectation values of the multitrace operators in the Gaussian theory, which we denote as t n = Ω n (a) (0) .

(2.41)
One obvious fact is that t n = 0 for n odd . For n even, instead, they are non-vanishing and one can explicitly evaluate their expressions starting from the initial condition t 0 = N and a set of recursion relations of the form and so on. These relations follow [40] from the fusion/fission identities satisfied by the su(N ) generators T b , namely for two arbitrary (N × N ) matrices A and B. In fact, using these identities one can recursively relate any correlator t n to the combination of correlators obtained after a single Wick contraction with the propagator (2.29).

JHEP09(2020)116
3 Large-N limit from the recursion relations We now consider the 't Hooft limit in which N → ∞ with λ = g 2 N (3.1) kept fixed. As argued in the previous section, all relevant observables can be expressed in terms of the quantities t n defined in (2.41). Therefore, as a preliminary step, we study the large-N limit of the latter.

Basic ingredients
Single traces: eq. (2.42) implies that the odd single traces have an identically vanishing v.e.v.: The v.e.v. of the even single traces, t 2k , can be computed by applying Wick's theorem and taking into account that all contractions in which some propagators cross, are suppressed in the large-N limit. In other words, only rainbow diagrams count [6] and, up to subleading terms in the 1/N expansion, one finds where C k are the Catalan numbers which enumerate the distinct rainbow diagrams 5 .
Double traces: the mixed even/odd double traces are clearly vanishing due to (2.42): Thus we have to consider only the even and the odd double traces. At the leading order in the large-N limit, the even double traces factorize: To show this, one can consider the relation (2.43b) and observe that the term in the righthand side proportional to n 1 , obtained when a propagator connects the first component with the second component, is subleading with respect to the first term in which the propagator remains within the first component. This fact leads to the factorized result (3.6). 5 The generating function of the Catalan numbers is

JHEP09(2020)116
The next-to-leading terms determine the connected part of the double trace v.e.v.'s. They are defined as and at large N they behave as Comparing (3.6) and (3.8), we see that the connected v.e.v. t c 2k 1 ,2k 2 is suppressed by a factor of 1/N 2 with respect to t 2k 1 ,2k 2 .
Also in the odd case we can use the recursion relations (2.43), but we cannot discard the terms that superficially look subleading. Clearly, the odd double traces coincide with their connected part because of (2.42), and at large N they are given by 6 In what follows it will turn out to be convenient to write the odd double trace v.e.v. as is closely related to the so-called Hilbert matrix.
Multitraces: let us now consider the v.e.v. of the multitraces. When there is an even trace component, at large N this factorizes as follows t 2k 1 ,n 2 ,n 3 ,... = t 2k 1 t n 2 ,n 3 ,... , (3.12) and the subleading contributions are suppressed by a factor of 1/N 2 . When all components are odd (with a total even number of factors), the recursion relations at large N imply that the full result is obtained by pairing all components in all possible ways, and replacing each pair by its expectation value (3.9). In other words, we have where H(k 1 , k 2 , k 3 , k 4 , . . .) represents the total Wick contraction computed with the "propagator" H k i ,k j . For instance, with 4 components we have Similarly, with 6 components we have the sum of the 15 possible ways to make a complete contraction, and so on. Repeatedly using the above results, one can obtain the large-N expansion of all observables in a quite explicit and detailed form.

The N = 4 SYM theory
Let us begin by illustrating the procedure for the N = 4 SYM theory. In this case the matrix model is purely Gaussian and thus, following the convention introduced above, we use a label (0) to distinguish its observables. The results we present in this subsection are well-known, but it is convenient to briefly review them before moving to the N = 2 cases of our interest.
The BPS Wilson loop: the v.e.v. of the Wilson loop (2.36) in the N = 4 theory is (3.15) Using (3.3), in the large-N limit one gets: where I is the modified Bessel function of the first kind. This well-known result was first obtained in [6] by resumming ladder diagrams in perturbation theory.
Single-trace operators and their mixing: let us consider the single-trace operators Ω n (a) = tr a n (3.17) with n ≥ 2, and subtract their v.e.v., so as to work with a basis formed by the vevless operators This is convenient because, after this redefinition, which is trivial when n is odd, the even and odd cases are on the same footing. Indeed, in both cases their two-point functions coincide with the connected correlators we introduced above: whose expression in the large-N limit is given in (3.8) and (3.9). Our goal is to apply the Gram-Schmidt diagonalization procedure with respect to the pairing (3.19) and find the normal-ordered version of the single-trace operators Ω n (a), which we denote by O n (a). The set of operators of dimension lower than n comprises both Ω p (a) with p < n, and multitrace operators with a total dimension less than n. However, in the large-N limit these multitrace operators do not play any role, because the factorization and the Wick-like expansion properties of the expectation values (3.12) and (3.13) imply that if one imposes the vanishing of the correlators between O  is enough to run the normal-ordering procedure considering only the single-trace operators. In this way, (2.40) reduces to Given the structure of the two-point functions (3.19), the even operators O (3.21) Actually, the result can be given in closed form for any n as follows: with the understanding that Ω 1 (a) = 0, which is true in the SU(N ) theory. Using in agreement with the results of [35].
Two-point functions: the two-point functions of the normal-ordered operators are indeed diagonal: as it can be proven by using (3.22) and (3.19). We note that the Gram-Schmidt procedure yields an expression of G (0) n directly in terms of the elements of C (0) : . (3.26)

JHEP09(2020)116
One-point functions in presence of the Wilson loop: also the one-point functions of the operators O (0) n (a) with the Wilson loop W C can be easily computed in the large-N limit. Indeed, using (3.22) and expanding the Wilson loop operator (2.36), we find Inserting the large-N behavior of the connected correlators given in (3.8) and (3.9) and retaining the leading contributions for N → ∞, after simple algebra we can recast the above sum as an expansion in powers of λ which can be resummed into a modified Bessel function I n . More precisely, we have This result was originally obtained in [8] by resumming rainbow diagrams in the planar limit.

The ABCDE theories
The matrix model for the N = 2 conformal theories of table 1 contains an interacting action S int (a) given in (2.32), and the v.e.v. of the various observables are computed according to (2.35). Using (2.34), S int (a) can be written as a sum of terms that are either quadratic or linear in the single-trace operators (3.17). We find convenient to split this sum into three parts as follows: (3.29) Here S odd (a) contains odd double traces:

30)
S even (a) contains even double traces: (a) contains single-trace operators, all even: Finally, we have introduced the rescaled 't Hooft coupling to make the end results more compact.

JHEP09(2020)116
Single-trace operators and their mixing: to obtain the normal-ordered operators O n (a) we have to repeat the Gram-Schmidt procedure and diagonalize the matrix of the two-point functions of the single-trace operators Ω n (a) computed in the interacting matrix model. To take advantage of the calculations already performed, we proceed in two steps and start from the operators O n,m is the term of order k in S int . In particular, at the first order we have The corrections terms G (k) n,m make the two-point functions (3.34) non-diagonal. We have therefore to rerun the Gram-Schmidt procedure and redefine the operators. At large N several important simplifications occur. We illustrate them by considering the N dependence of the various terms contributing to G (1) n,m , but these arguments can be readily extended also to the higher correction terms G m (a) decomposes into a sum of terms which contain a couple of single-trace operators Ω r (a) Ω s (a). The tree-level contribution of any such term to the two-point function G n,m is proportional to t r,s . Its contribution to the first correction G (1) n,m depends on which part of the interacting action one considers. The single-trace part S s.t. given in (3.32) corresponds to an insertion of Ω 2p+2 (a) accompanied by a factor of 1/N p+1 , so that the contribution of Ω r (a) Ω s (a) is proportional to where the second step follows from (3.12) and (3.3). This correction is therefore subleading in N with respect to the tree-level result t r,s . This fact means that all contributions arising from the single-trace part of the interacting action can be ignored in the large-N limit. Since S s.t. is the only part of S int that does not depend on the parameter ν, it follows that in the large-N limit the two-point correlators will only depend on ν and not on the more specific matter content of the N = 2 theory. Thus, G n,m will be equal for the B and C models, and for the D and E models.
Let us now consider the corrections coming from the even double-trace part S even of the interaction action given in (3.31). In this case, the interaction produces an insertion of Ω 2k (a) Ω 2p−2k+2 (a) accompanied by a factor of 1/N p+1 , so that the typical term goes like 1 N p+1 t r,s,2k,2p−2k+2 − t r,s t 2k,2p−2k+2 ∝ t r,s . (3.37)

JHEP09(2020)116
Again we have used again the factorization property (3.12) and the expression of the onepoint function (3.3). We see that now the correction scales in the large-N limit just like the tree level result. So the even part S even , which is proportional to ν, will always contribute to the corrections of the two-point functions.
Finally, we consider the odd double-trace part S odd defined in (3.30). Here we have to distinguish two cases: when the operators Ω r (a) Ω s (a) are both even and when they are both odd. In the first case, which occurs in the two-point functions G n,m of two even operators, the typical contribution due to S odd has the following behavior at large N for r and s even. This result follows again upon using (3.12) and (3.9). Thus, being subleading with respect to the tree-level term, the odd part S odd has no effect on the correlators of two even operators in the large-N limit. Instead, it corrects the correlators of two odd operators. Indeed, in this case the typical contribution due to S odd is of the type for r and s odd. Here we have taken into account the structure of odd multitraces described in (3.13) and (3.14). This correction has the same N -dependence of the tree-level term and survives in the large-N limit. To summarize, we have shown that the correlators of two even operators in the large-N limit only receive corrections from S even , while those of two odd operators are corrected both by S even and by S odd . This means that the N = 2 correlators of two even operators differ from the N = 4 ones by terms which are proportional to ν. In particular, for the D and E theories, for which ν = 0, this difference vanishes. Thus, in the sector of the even operators, these two N = 2 models are indistinguishable from the N = 4 SYM theory in the planar limit. On the contrary, the correlators of odd operators have corrections proportional to (2 − ν) arising from S odd and corrections proportional to ν arising from S even , and these are non-trivial even for the D and E theories.
Some explicit examples: we have at our disposal all elements to carry out explicitly the Gram-Schmidt procedure starting from the operators O while, applying the analogue of (3.26), we obtain the diagonal two-point functions which we parametrize as follows:

JHEP09(2020)116
Here we have factorized the N = 4 result (3.25), so that γ n , which is a function of the 't Hooft coupling, reduces to 1 in the limit λ → 0. This procedure is entirely algorithmic and it is easy to reach quite high orders in the coupling and in transcendentality. Here we report the explicit results for the lowest dimensional operators, showing only the lowest terms in their expansion to avoid excessive clutter.
• At dimension 2 the operator O (0) 2 (a) can only mix with the identity and we find that the corresponding correction is where the ellipses stand for terms of higher order in λ. Computing the two-point function of O 2 (a), we find that the correction factor γ 2 is • The operator O The corresponding two-point coefficient γ 3 is (3.47) The two-point function is captured by 3 (a) and the mixing term is given by

JHEP09(2020)116
while the two-point function coefficient is 7 All these results explicitly show the pattern discussed above, namely that the even operators and their two-point functions are not corrected with respect to the N = 4 theory when ν = 0. On the contrary the odd operators and their two-point functions are different from the corresponding ones in the N = 4 theory even when ν = 0. In particular we point out that in the D and E theories, the first correction to γ 3 is proportional to ζ(5) λ 3 , while the first correction to γ 5 is proportional to ζ(9) λ 5 . In section 7 we will give a diagrammatic explanation of this fact and of its generalization to the two-point functions γ 2k+1 with k > 2.
The BPS Wilson loop: we now consider the Wilson loop operator (2.36). Its v.e.v. in the N = 2 theories is given by .
Expanding the exponentials and using the explicit form of the interacting action (3.29), one gets a non-vanishing contribution only when k is even. Thus, the difference with respect to the N = 4 result (3.15) can be written as (3.52) This quantity greatly simplifies in the large-N limit. In fact, using the same arguments explained above, one can show that the single-trace part of the interacting action (3.32) does not contribute at leading order and that the same thing happens for S odd . Thus, one is left only with the contributions arising from the even piece of the action, S even , which can be evaluated using the large-N factorization property (3.12) together with (3.8). Collecting all terms contributing to a given ζ-value, one finds a series in λ which can be resummed into a combination of Bessel functions I n . Explicitly, the very first few terms are 8 (3.53) 7 We report this result to a higher order in λ with respect to the previous ones to exhibit the fact that it does not vanish at ν = 0. 8 The term proportional to ζ(3) was already obtained in [41,59] for the A theory, corresponding to ν = 1.

JHEP09(2020)116
where for brevity we have defined Terms with higher powers of λ for which the ellipses in (3.53) stand, can be systematically computed without any difficulty.
One-point functions in presence of the Wilson loop: in presence of the Wilson loop W C , the normal-ordered operators O n (a) have a non-trivial one-point function which is given by .
To compute the difference of these quantities with respect to the N = 4 theory, we have to expand the exponentials and also to take into account that the operators bear a dependence on the interacting action because of the normal ordering that we discussed in the previous part of this section. We then find When we take the large-N limit, drastic simplification occur in these quantities and the results can be compactly written in terms of the rescaled Bessel functions (3.54). For example, for n = 2, 3 we find (3.58) Similar expressions for higher values of n, as well as the contributions at higher orders in λ, can be obtained without any problem since the whole procedure is purely algebraic.

Large-N limit from the eigenvalue distribution
We now discuss the Cartan algebra approach based on the integration over the matrix model eigenvalues in the large-N limit.

JHEP09(2020)116
A specific application of this approach to the N = 2 ABCDE theories has already been presented in [23] (see also [34]), where it is shown that the ratio ν defined in (1.1) is the unique relevant parameter at large N . The study in [23] focused mainly on the Wilson loop v.e.v. and its eventual extrapolation at strong coupling. Here, instead, we discuss how this method can be applied to the calculation of the two-point functions of chiral primaries and of their one-point functions in presence of a Wilson loop. In doing so, we also provide new results about the large-N mixing of single-trace operators in terms of N = 2 deformed orthogonal polynomials, setting in a broader perspective the findings of [35]. Our analysis also leads to an alternative and more compact derivation of the results obtained by the full Lie algebra methods described in the previous section which, at least in principle, may be a convenient starting point for a non-perturbative investigation.

The large-N universal integral equation
(4.1) Here the function K(x) is defined by where ψ(x) is the digamma function. Notice that in N = 2 superconformal theories, the linear part of K(x) drops out in the saddle-point equation and thus it may be conveniently subtracted from the start. 9 Moreover, exploiting the expansion (2.31), for small x we have 3) The large-N limit of (4.1) is captured by an integral equation for the continuum limit of the discrete density Assuming that for N → ∞ the eigenvalues m u condense on a single 10 segment [µ − , µ + ] ⊂ R, one finds that ρ and µ ± are determined by the following equation [23]:

JHEP09(2020)116
where we have understood the prescription for taking the Cauchy principal value of the integral and have introduced the parameter ν as in (2.9), together with the normalization condition From the solution to (4.5) we can readily compute the v.e.v. of single-trace operators. For example, the integral representation of the half-BPS Wilson loop v.e.v. is The right hand side of (4.5) is odd under x → −x and thus it is consistent to assume a symmetric density ρ(x) = ρ(−x) supported on a symmetric cut with µ + = −µ − ≡ µ. With these assumptions, we may replace K(x + y) with K(x − y) under integration and rewrite (4.5) in the much simpler form This implies that, whenever (4.8) may be used, any ν = 0 model gives the same results as the N = 4 SYM theory where the exact density is the Wigner semi-circle distribution . (4.9) Using this distribution in (4.7), one finds in agreement with (3.16). Thus, for the ν = 0 theories we have ∆w = 0 in agreement with (3.53). However, the symmetric one-cut assumption is too restrictive. Indeed, to compute observables that are more general than (4.7), like for instance two-point functions or one-point functions in presence of a Wilson loop, we need to extend the action by suitable sources coupled to the operators that are inserted in the correlators. This leads to asymmetries in the cut when one considers "odd" operators, i.e. operators involving traces of odd powers of the matrix model variable. Moreover, intermediate calculations require to consider the unrestricted integral equation (4.5) that does not reproduce the N = 4 results for ν = 0. In other words, one cannot expect that such "odd" observables in the ν = 0 theories are equal to those in the N = 4 SYM theory. Several examples of this claim will be illustrated and discussed later.

Single-trace mixing at large N and two-point functions
As we remarked in the previous section, an important bottleneck is the operator mixing that seems to require an ad hoc treatment depending on the operators under study. We now JHEP09(2020)116 discuss how to deal in general with the mixing in the single-trace sector starting from the integral equation (4.5). Then, as an application, we recompute the two-point functions and the one-point functions in presence of the half-BPS Wilson loop. As a technical point, note that we shall work in the U(N ) matrix model which is expected to give the same results as the SU(N ) model in the large-N limit, since the difference between enforcing or not the tracelessness condition turns out to be subleading at large N [60]; see also appendix A. Indeed, we will perfectly reproduce the previous results obtained in the SU(N ) matrix model following the full Lie algebra approach.
To evaluate the two-point functions of single-trace operators and eventually impose the orthogonality conditions, we need to add a set of source terms to the effective action of the matrix model. We thus modify the action S given in (2.20) according to where σ n is the source term for a monic polynomial P n in the matrix M of degree n: P n = M n + . . . . with respect to the pairing n, m ≡ tr P n tr P m − tr P n tr P m . (4.14) In presence of the sources σ n , the eigenvalue distribution satisfies a modified integral equation, that we will consider below.
Since the normal-ordered operators have to be orthogonal with respect to the identity, P n must also satisfy tr P n = 0. This can be left to the end. This is because the connected two-point function in the above pairing is unchanged if we replace tr P n by tr P n + c n . As a result we can impose (4.13) and determine P n . The normal-ordered operator is then tr P n − tr P n .

The N = 4 SYM theory
Let us begin with the simplest case of the N = 4 SYM theory. In previous sections, we have adopted the convention that quantities in the N = 4 case are distinguished by a (0) super/sub-script. Here, however, to avoid excessive clutter, we suppress this index where there is no risk of confusion.
We start by considering even sources σ n with n ∈ 2N. They modify the integral equation (4.8) with ν = 0 as follows: Both the density and the cut edge now depend on σ n .

JHEP09(2020)116
In terms of this modified density, the pairing is expressed as 11 (4. 16) Taking into account that the density vanishes at the source-dependent edges of the cut, this becomes where µ 0 = µ(0) and To determine ρ m , we take a derivative of (4.15) with respect to σ m , obtaining Since the differentiated density is expected to be unbounded at the cut edges, the general solution to (4.19) is [61] ρ m (x) = C m where C m is an arbitrary constant 12 . Inserting (4.20) into (4.17), we get (4.21) The orthogonality condition n, m ∝ δ n,m is realized if we take P n to be proportional to the Chebyshev polynomial of the second kind T n . Indeed, in such a case the contribution proportional to C m is identically zero (so that we can safely set C m = 0) and the double integral in (4.21) vanishes for n = m.
To see this, let us observe that integrating first over x gives a result proportional to the Chebyshev polynomial of the first kind U n−1 (y) (see (B.3) in appendix B). On the other hand, we have T n (x) = n U n−1 (x), and thus (4.21) reproduces the orthogonality relation of polynomials U n with respect to the weight µ 2 0 − y 2 . In conclusion, in the N = 4 SYM theory, taking into account the normalization choice in (4.12), we have 11 Note that the density integral provides a trace divided by a factor of N which is compensated by differentiation with respect to σ which inserts N times the trace of the associated field. 12 Indeed, (µ 2 0 − x 2 ) −1/2 is a zero-mode of the Cauchy integral kernel 1/(x − y) in (4.19).

JHEP09(2020)116
where in the second step we used (4.9). Plugging this into (4.20), we get while the pairing becomes n, m = G (0) n δ n,m (4.24) with where the last step follows from (4.22) and the orthogonality properties of the Chebyshev polynomials.
This result is in full agreement with [36] where it is shown that (4.22) and (4.25) hold also for odd n. It is also in complete agreement with the results of the full Lie Algebra approach given in section 3.2. Indeed, comparing (4.22) with (3.24), and taking into account the rescaling (2.24) between the matrices M and a, we see that for n ≥ 2  In a perfectly consistent way, the two-point functions (4.25) and (3.25) are related as follows:

The ABCDE theories
In a N = 2 theory of the ABCDE series with parameter ν and generic (even or odd) sources, the modified integral equation that determines the eigenvalue density reads +µ(σ)+c(σ) −µ(σ)+c(σ) where we have allowed for a non zero cut center c(σ). Of course c(σ) = 0 if all sources are even. The orthogonality condition has the same form as in (4.13) and the pairing is expressed in terms of the eigenvalue distribution by the analogue of (4.16). Thus we must find polynomials P n such that where again µ 0 ≡ µ(0). Taking a derivative of (4.28) with respect to σ m yields +µ 0 we obtain, due to the restrictions in (4.31), the diagonal values The coefficients a  13 It would be interesting to understand more deeply such a deformation, as it happened in a different context for the hypergeometric polynomials [62]. 14 Convolutions are conveniently evaluated by The next step is finding the expression of µ 0 .
Determination of µ 0 : to obtain explicit expansions in powers of λ we need to determine the expression of the cut endpoint µ 0 (λ). To this aim, we restart from (4.28) without sources, i.e. from (4.8) with µ → µ 0 . Since we are interested in the contributions coming from a finite set of terms in the expansion (4.3), we can conveniently represent the density, which is an even function, in the form where the sum is over a finite number of terms and a 0 is fixed by the normalization condition Then, we can use the formula which follows from the second relation in (B.4). In this way, replacing (4.36) in (4.8) gives an algebraic equation for µ 0 that can be easily expanded to any desired order (see appendix (C) for an efficient algorithm and details). With this procedure we obtain the following expansion in terms of the rescaled 't Hooft coupling λ:

JHEP09(2020)116
where γ n encodes the deviation of the N = 2 result from the N = 4 one given in (4.25).
For the first few values of n we explicitly find Additional data for γ n with n = 5, . . . , 11 are collected in appendix D.

One-point functions in presence of the Wilson loop
The previous analysis provides us with all necessary ingredients to compute also the v.e.v. of the half-BPS Wilson loop and the one-point functions defined in (3.55).
The v.e.v. of the Wilson loop may be obtained from the density in (4.36) and reads where the last step follows form the second equation in (B.5). When a chiral operator is inserted we can proceed as for the two-point functions in (4.16) and find

JHEP09(2020)116
where the differentiated density is given in (4.32). Exploiting the first relation in (B.5), it is easy to obtain This result is perfectly consistent with what we obtained within the full Lie algebra approach in section 3 (see in particular (3.55)). In fact, the two results differ just by an overall factor due to the rescaling factor (2.24), namely We now give some explicit results.

The N = 4 SYM theory
In N = 4 SYM we just take µ 0 = √ λ 2π , a k = 0 for k = 0 and a 0 given by the normalization condition (4.37). With this input, it is immediate to find from (4.44) the well-known N = 4 result (3.16). On the other hand, from (4.46) using a (n) n = 1 for n > 1, we get which is consistent with (3.28) and the results of [8].

The ABCDE theories
When µ 0 is non-trivial, the result (4.46) may be formally expanded in the Riemann ζ-values and each term captures the exact dependence in λ. Proceeding in this way, we find that the difference of the Wilson loop v.e.v. with respect to the N = 4 result exactly matches the expression in (3.53).
In the case of insertions, the above procedure yields where for n = 2, 3 the quantities ∆w 2 and ∆w 3 precisely match the differences (3.57) and (3.58) computed in the full Lie algebra approach. Additional data for n = 4, 5, 6, 7 may be found in appendix D.

Asymptotic correlators of operators with large dimension
Inspection of the explicit expressions for γ n and ∆w n shows that each transcendentality structure depends on n in a simple way when n is large enough. This is because beyond some point, i.e. for n greater than a certain n(ζ) depending on the specific Riemann ζ-values, the coefficients a The detailed analysis presented in this section shows the complete equivalence of the methods based on the use of the distribution of the matrix eigenvalues in the large-N limit and those of the full Lie algebra approach based on the recursion relations satisfied by the v.e.v.'s of multitrace operators, that we discussed in section 3.

Part II
In the second part of this paper we study in detail the D and E superconformal theories in the large-N limit. These models, which both have ν = 0, are known to have a holographic dual description in terms of an orbifold of AdS 5 × S 5 [26], and in some sense they can be regarded as the next-to-simplest theories after the N = 4 SYM theory. It would therefore be extremely interesting to be able to extrapolate the various observables which so far we have discussed in a weak-coupling perturbative approach and see what one can say in the planar limit at finite or strong coupling. As we shall see in the following, remarkable simplifications occur when ν = 0 and some resummation methods can be applied to our perturbative expansion. 5 The full Lie algebra approach for the ν = 0 theories As discussed in section 3, when ν = 0 the interaction action of the matrix model in the large-N limit simply reduces to the odd part S odd (a) given in (3.30). Moreover, the v.e.v. of the Wilson loop and the two-point functions of even single-trace operators are not corrected with respect to their N = 4 values. We concentrate therefore on the odd single trace JHEP09(2020)116 operators and analyze their two-point functions and their one-point functions in presence of the Wilson loop. The key property we will exploit is the Wick-like factorization (3.13) of the expectation values of such odd operators in the Gaussian model at large N . The linear relation (3.22) ensures that this property also applies to correlators of the operators O (0) 2i+1 (a), which are normal-ordered with respect to the Gaussian measure. Their correlators are therefore diagonal, as one can see from (3.25).
Let us rescale them by setting (3.25)). At large N , the operators ω i (a) have a canonical two-point function and the correlators of many of them are again computed using Wick's theorem. In other words, we can regard the matrix operators ω i (a) as a set of real variables ω i normally distributed. Indeed we can write where we have denoted by ω the (infinite) column vector of components ω i and have defined The rules (5.2) and (5.3) exponentiate, so that for any constant matrix A we have In order to efficiently compute the observables involving the odd single-trace operators, it is convenient to rewrite the interaction action S odd (a) in terms of the quantities ω i (a) we have just introduced. To this aim, we first invert the relation (3.22), which for odd operators leads to we then rescale the operators according to (5.1) and plug everything in the action (3.30). After carrying out the algebra, the resulting expression takes the form where X is an infinite numerical matrix whose entries are given by It is interesting to observe that the matrix elements X ij can be given an integral representation in terms of the Bessel functions of the first kind J (x); indeed one can show that This expression resums the perturbative expansion in λ given (5.8), which in turn can be recovered by Taylor expanding the Bessel functions and performing the resulting t-integral using The definition (5.2) and the property (5.5) allow us to rewrite this expression in terms of ordinary gaussian integrals as follows: Two-point functions: using this free-field formalism, we now reconsider the computation of the observables γ 2i+1 that parametrize, according to (3.41), the two-point functions of the normal ordered odd operators O 2i+1 (a). The starting point is finding the correlators of the operators O (0) 2i+1 in the N = 2 theories with ν = 0, to which we have to apply the Gram-Schmidt procedure. Let us thus define the matrix G with elements

JHEP09(2020)116
Applying the Gram-Schmidt procedure one can express the observables γ 2i+1 in terms of the correlators G ij by means of the analogue of (3.26). Taking into account (3.41), we find Here G (i+1) is the submatrix of G comprising its first i rows and columns, namely the matrix of correlators of the operators O (0) 2k+1 with k ≤ i. It then follows that This ratio of determinants can be rewritten in a different way taking into account the expansion ( − X) −1 = + X + X 2 + X 3 + . . . (5.19) and introducing the infinite matrices X (i) ≡ X with the first i − 1 rows and columns deleted. (5.20) Of course, X (1) = X. Then one has This expression, together with the form of the matrix X given in (5.10), is very powerful. In fact, it readily reproduces the formulae in (3.46), (3.50) and (D.1)-(D.7) for γ 3 up to γ 11 for ν = 0, and also it yields in a quite straightforward way the expansions to very high orders in λ for any desired value of the index i and any power of Riemann ζ-values. In the next section we will show how one can obtain equivalent expressions using the eigenvalue approach discussed in section 4 and provide also a few explicit examples of expansions to high orders. We conclude by observing that it is also possible to exploit the Gaussian variables ω i to express in a very efficient way the observables w 2i+1 defined (3.55), which are related to the one-point functions of odd operators in presence of the BPS Wilson loop. To avoid redundancy, however we will discuss them only in the eigenvalue approach in the next section.
6 The eigenvalue distribution approach for the ν = 0 theories In the previous section 5, we have discussed the ν = 0 theories D and E in the "full Lie algebra" approach. Here, consider them in the "Cartan algebra approach". We discuss in full details both the calculation of the two-point unctions and that of the one-point functions in presence of the Wilson loop.

JHEP09(2020)116
Two-point functions: in the D and E models, the cut edge µ 0 given in (4.39) simply reduces to µ 0 = 2 λ = √ λ 2π . The two-point functions (4.34) read therefore Comparing with (4.25) and (4.40), we see that the coefficients a (n) n yield in this case directly the quantities γ n which encode the deviation of the N = 2 result from to the N = 4 one: Let us recall that the coefficients a (n) m appear in the differentiated density (4.31), which in the present case is determined by the integral equation (4.30) with ν = 0, namely Since we are interested in odd chiral primaries with n = 2k + 1, according to (4.32), the differentiated density ρ n has an expansion in Chebyshev T -polynomials of odd degree. After inserting the Ansatz (4.32) in the left hand side of (6.3), the part of the integral operator that depends on the function K can be written as We can now use the following integral representation for K(x) valid in the strip | Im(x)| < 1, given by and obtain This expression can be expanded in the even Chebyshev U -polynomials as follows where (6.10) With some work, it is possible to evaluate this expression in closed form using the Bessel functions of the first kind and to show that The matrix Y is related to the matrix X introduced in (5.10) in a simple way: Inserting this result into (6.8) and (6.4) we can solve the integral equation (6.3) for the coefficients a (2k+1) 2j+1 , taking into account the ansatz (4.31) for the polynomials P 2k+1 . After some straightforward manipulations, we find that the vector 2k+3 , · · · is given by where C = (1, 0, 0, . . . ) and Y (k) is the matrix obtained from Y by removing its first k − 1 rows and columns. Using this into (6.2), we conclude that the two-point coefficients γ 2i+1 are given by This result is fully equivalent to (5.21). For the first few values of i we get explicitly + 368235 2 ζ(7) ζ(9) + 1000 ζ(5) 3 λ 9 + · · · , (6.15a) ζ(17) λ 9 + · · · , (6.15b) We have reported here only the first few terms of the expansions but, as we will discuss in section 8, using (6.14) we can push the expansions to very high orders with minor effort.

JHEP09(2020)116
One-point functions in presence of the Wilson loop: the one-point functions of the chiral odd operators in presence of the circular Wilson loop, given in (4.46), can be written for the ν = 0 theories as follows ("t" denotes matrix transposition) where B k is the k-th component of the vector with I being the modified Bessel functions of the first kind. Taking the difference with respect to the N = 4 results given in (4.48), we obtain the corrections ∆ w 2i+1 which are related to quantities ∆w 2i+1 computed in the full Lie algebra approach by a simple rescaling due to the relation (2.24), namely We report below the explicit expressions of ∆w 2i+1 we obtain from the above procedure for the first few values of i: difference between correlators with even and odd operators and in particular to trace the diagrammatic origin of the fact that the even correlators do not receive corrections in the large-N limit with respect to the N = 4 case, while the odd correlators do. We mainly concentrate on the two-point function (2.3) of single-trace operators O n (x) = tr ϕ n (x), namely Superconformal symmetry guarantees that the space-time dependence 1/(4π 2 x 2 ) 2n is preserved at the quantum level, so we have simply to compute the coefficient G n and compare it with the matrix model results. In particular our achievements are twofold: • We provide a systematic diagrammatic analysis up to three loops for the two-point function coefficients G 2 , G 3 , G 4 and G 5 .
• We infer the general contribution for the first deviation from the N = 4 theory of the coefficient G n for generic n.
We use the same tools which are frequently used in supersymmetric contexts [20,33,[40][41][42]53], namely the N = 1 superspace formalism and the diagrammatic difference between N = 2 and N = 4 theories. We refer to [42] for a detailed account of the Lagrangian and the Feynman rules that are needed for the computations.

Field theory calculations up to 3-loops
We start from the operator insertions represented in figure 1. The tree-level contribution to the two-point function is obtained by contracting the legs in this diagram in all possible ways using the tree-level propagators. Doing so, we obtain In the N = 4 SYM theory all higher-loop corrections cancel and the tree-level result (7.2) represents the full answer to the correlator in the large-N limit, matching the matrix model result (3.25). The N = 2 theories, instead, contain a full tower of quantum corrections.

JHEP09(2020)116
a2 a1 Here we follow the same notation introduced in (3.42) and factorize the N = 4 result (7.2), so that any N = 2 two-point function coefficient is uniquely described by γ n which is a function of the 't Hooft coupling. The loop corrections of N = 2 theories can be effectively worked out in the diagrammatic difference with N = 4. This procedure allows us to discard all diagrams which only contain fields from the N = 2 vector multiplet, since they are in common with the N = 4 SYM theory, and to isolate the genuine N = 2 contributions, corresponding to the diagrams where hypermultiplet are present. As we see from figure 1, the operator insertions only contain fields belonging to the vector multiplet, hence the hypermultiplets appear only inside loops. The diagrams surviving in the difference can then be organized in terms of some building blocks, as displayed in figure 2. This procedure reduces the number of Feynman diagrams to be computed and allows a rapid comparison with the matrix model. Indeed each building block carries a color factor C (n) , where n counts the number of adjoint fields attached to hypermultiplet loops. In the difference theory this color factor is given by the combination C a 1 ...an = (Tr R − Tr adj ) T a 1 . . . T an (7.4) which exactly reproduces what we found in the matrix model (see (2.18)). We now systematically classify the Feynman diagrams by looking at their color factors. Many of them vanish because of the matter content R of the ν = 0 theories or can be discarded because they are subleading in the large N limit.  One loop: at order g 2 there is a unique way to insert a hypermultiplet loop inside figure 1, namely by using the one-loop correction to the scalar propagator. This diagram corresponds to the first building block of figure 2, but its color factor C (2) vanishes because of conformal invariance due to the condition (2.33).

JHEP09(2020)116
Two loops: since C (2) always vanishes, the diagrams at order g 4 can only be built out of the C (3) and C (4) building blocks. Any diagram proportional to C (3) is again proportional to the coefficient β 0 of the β-function, and so it vanishes in any superconformal theory. There are only two diagrams which can be built out of the third diagram of figure 2 with a C (4) color factor. They have been fully analyzed in the conformal SQCD, which is theory A, in [40,41] and are represented in figure 3. These diagrams do give a non-trivial contribution to correlation functions for ν = 0 theories, and actually account for the terms linear in ζ(3) in γ n . However, if we analyze the color factor C (4) in the D and E theories by expanding the trace combination in (2.34) and using FormTrace [63], we obtain that C (4) is always subleading in the large-N limit. This fact explains why the ν = 0 theories have no two-loop order term proportional to ζ(3) in the large-N limit.
Three loops: at order g 6 we analyze all diagrams that could provide a ζ(5) term (we do not consider diagrams made of ζ(3)-subdiagrams). They can arise only from the C (4) , C (5) and C (6) building blocks. Such a list is rather long as one can see from figures 4, 5 and 6. However, in the ν = 0 theories all these diagrams have a color factor which is either zero or subleading in N .
The unique special diagram that deserves to be analyzed in detail is the one depicted in figure 7, which we call the "hexagon diagram" since it gives rise to an hexagon when it is unfolded. We analyze its contribution in the following subsection. JHEP09(2020)116

Hexagon diagram in even and odd correlators
First of all, it is obvious that the hexagon diagram cannot be inserted inside G 2 . This fact confirms that for the ν = 0 theories in agreement with the matrix model result.
We now evaluate the contribution of the hexagon diagram to G 3 and G 4 . Remembering that each vertex brings a power of √ 2g and computing the corresponding color factor in the large-N limit, we obtain the following prefactors This shows that the hexagon diagram is planar only inside G 3 whereas it is subleading, and hence non-planar, inside G 4 . Therefore we conclude that for ν = 0 the factor γ 3 receives a three-loop correction, whereas γ 4 remains 1, at least up to four loops, namely

JHEP09(2020)116
This diagrammatic analysis highlights the difference between even and odd correlators. For the two-point functions of even operators we expect the same cancellations in the planar limit also at higher loops, so that the total result for any even correlator in the ν = 0 theories is equal to the N = 4 one, confirming the matrix model predictions.
Let us now return to G 3 and γ 3 . The relevant space-time integral W 6 (x) contributing to the two-point function of O 3 has been computed in appendix C of [40] and is (7.8) The space-time dependence is the expected one because of conformal invariance as prescribed by (7.1), and so we can just focus on the remaining terms. Multiplying by K 3 given in (7.6) and taking into account a symmetry factor of 3! corresponding to all possible ways of inserting the hexagon inside the correlator, we get: This result confirms the matrix model prediction (6.15a) up to three loops. Notice that if we had inserted the hexagon diagram inside correlators with odd operators with higher order, we would have obtained a non-planar result. This implies that the coefficient G 2k+1 with k > 1 starts deviating from the N = 4 result at higher-loop orders.

Higher loops for higher order correlators
We are able to generalize the previous reasoning to higher correlators. In particular we show that for each correlator in ν = 0 theories the first deviation from N = 4 theory is given by the insertion of the following diagram: This is a hypermultiplet loop with 2k adjoint scalar legs, and represents a generalization of the hexagon diagram. The result (7.10) is obtained by exploiting a map with the k-loop ladder diagrams contributing to the four-point function of a φ 3 -scalar theory, which was computed in [64]. Following the reasoning of figure 2, the color factor arising from the insertion of the diagram with 2k legs is the totally symmetric tensor C (a 1 ...a 2k ) = Tr R − Tr adj T (a 1 . . . T a 2k ) . (7.11)

JHEP09(2020)116
We again evaluate this color factor using FormTrace [63], and find that at large N it is given by This shows that this 2k-leg diagram can contribute to correlation functions in ν = 0 theories only for odd k.
Proceeding as before, we see that for each G k with k = 2i + 1, the first non-trivial deviation from N = 4 appears at O(λ 2i+1 ), and its full contribution is given by the insertion of the diagram (7.10) with (4i + 2) legs. We have explicitly computed the i = 2 and i = 3 cases corresponding to k = 5 and k = 7 respectively, finding These perfectly match the first perturbative corrections in the matrix model results (6.15).
Then it is quite easy to infer the general behavior for the ν = 0 theories: For each correlator we recognize the contribution coming from the integral (7.10) with 4i + 2 legs, while the coefficient 1 2 2i−1 is due to the multiplicity and color factors. This diagrammatic analysis can be readily generalized to the one-point functions of chiral primaries in presence of the Wilson loop, and also in this case one nicely recovers the first perturbative terms in perfect agreement with the matrix model results.

Resummation in the ν = 0 theories
The computational tools that we have developed in the previous sections are particularly efficient for the D, E models and allow us to generate perturbative expansions to a very high order without too much effort. In this section we try to analyze these long perturbative expansions in order to get some preliminary non-perturbative information. This is clearly a very important issue given the difficulties that are notoriously encountered in the strong coupling analysis of these models, at both numerical and analytical level [23,30,34,36].
We begin with a short recap of the available results. The strong-coupling scaling of the two-point function G 2 for SQCD (i.e. theory A with ν = 1) was considered in [36]. The same scaling was reproduced in [34] which extended the analysis to G 4 . For the correction factors γ 2 and γ 4 the results are that for large λ → ∞,

JHEP09(2020)116
Similarly the Wilson loop scaling was first considered in [30] for SQCD and the analysis was extended to general ν in [23], which obtained for large λ the following behaviors 16 Here, we point out two shortcomings of these results: • Apart from the case of the Wilson loop without insertions, it has not been possible to figure out the coefficient that should go in front of these scaling factors, which have been obtained primarily using the Wiener-Hopf method which is quite hard to control, see e.g. appendix D of [34].
• All these results are for the even sector of observables. This means in particular that they don't shed any light on the difference between the strong coupling dynamics of the D and E models compared to the N = 4 SYM theory.
Although the first point is beyond the reach of the numerical analysis of the long perturbative expansions we employ, we are going to explore the structure that can be seen from such an approach. We also take a first step addressing the second point, presenting strong numerical evidence that the two-point correction γ 3 and one-point correction ∆w 3 have a power-law growth.

The two-point function correction factor γ 3
We begin our analysis with the discussion of the correction coefficient γ 3 . In order to compute a long expansion, we use the method described in sections 5 and 6, and write (see for instance (5.21) for i = 1) where the matrix X is defined in (5.10). The powers X k can be computed using the sum rule

(8.4)
This relation gives the following pattern of iterated integrals 16 The Wilson loop v.e.v. w in the ν = 0 models is equal to the one in N = 4 SYM. The power correction, which is well-known from the matrix model solution, has been a hard test of AdS/CFT correspondence and has been recovered only quite recently in [65].   or its analyticity structure. This prevents us to perform any rigorous resummation of the perturbative expansion. Nevertheless, we can try to analytically continue beyond the convergence radius by applying the Padé-Borel resummation technique [66]. We evaluate a Padé approximant P [M/K] (λ) of the Borel-improved series, namely and then transform back according to: In the left panel of figure 9 we show what is obtained by considering three diagonal Padé approximations. In all cases, the reconstructed function coincides with the convergent perturbative sum when λ < λ c . Beyond λ c , the resummation appears to be well defined in the sense that its value at fixed λ stabilizes at increasing Padé degree M . Our data suggest that we can safely trust this reconstruction up to the rather large values λ   behavior appears to be ∆γ 3 ∼ C λ 1/4 for moderate λ, as shown in the right panel of figure 9. The exponent 1/4 is just a qualitative estimate in the considered range of coupling, since one can expect also logarithmic corrections as discussed in [34]. It would be very interesting to match such numerical indications by an analytic strong-coupling calculation.

The one-point function correction ∆w 3
The same numerical investigation can be worked out for ∆w 3 . From the perturbative expansion up to order λ 60 we obtain a finite convergence radius consistent with that of ∆γ 3 and a pattern which is very similar to that of figure 8. The smallest singularity of the Padé approximants is now shown in table 3 which strongly suggests a singularity at the same position as in ∆γ 3 (λ), namely at λ c /π 2 1. Finally, in figure 10, we present the results from the Padé-Borel resummation. We emphasize again that this kind of resummation is just a numerical exploration, given the lack of theoretical control on the asymptotic properties of the perturbative expansion. Nevertheless, the analytic continuation of our numerical data turns out to be reliable at least up to λ/π 2 10 (left panel). In this region, the asymptotic behavior appears to be ∆w 3 ∼ C λ 8 . The specific value of the exponent should not be interpreted as an analytic prediction that we don't have. It is just an estimate valid in the range suggested by the numerical data. To guide the eye, in the plot we have also included a dashed linear fit of |∆w 3 | 1/8 to the rightmost part of the data. Again, we remind that such asymptotic representation may well be a crude approximation to the actual answer, due to possible logarithmic corrections that are natural in this context, but which are out of the reach of the current numerical analysis.

JHEP09(2020)116 9 Conclusions and perspectives
In this paper we have considered a set of Lagrangian N = 2 conformal theories with SU(N ) gauge group. Following the localization procedure, we have reviewed the existing matrix model techniques for large values of the rank N . Developing both the full Lie algebra approach and the Cartan sub-algebra approach, we have computed a large set of observables in the planar limit. We have generated efficient algorithms to produce perturbative series up to very high loop orders, without any limitations due to the conformal dimension, extending some results already present in the literature.
In the second part we have concentrated on N = 2 theories whose fundamental matter content does not scale with N . We have shown several reasons why they can be considered as "the next-to-simplest" gauge theories. At the diagrammatic level it is immediate to reach the three-loops order in perturbation theory, and we can easily go even beyond for certain classes of observables. Moreover we have showed that many observables (like the Wilson loop v.e.v. and chiral correlators with even dimensions) are equivalent to those of the N = 4 SYM theory in the large-N limit. Only odd correlators feel the difference with N = 4 and represent a set of interesting observables to explore the gravity dual of this special N = 2 class of theories.
From this perspective, a natural continuation of the present analysis of the field theory side could be the study of these special observables in the strong coupling regime λ 1. To this aim, one should capture the evolution of the matrix model eigenvalue density by solving at large λ the associated ν = 0 integral equation. As we remarked in the main text, this requires dealing with various deviations from the N = 4 SYM case, like extra (odd) sources and the associated cut asymmetry, that play a role when computing the odd observables. At the moment, it is unclear whether this can be done analitically, for instance by Wiener-Hopf methods. The preliminary results presented in section 8 may be useful in this respect.
To give an example, we presented strong support for a finite convergence radius of the perturbative expansion at |λ c |/π 2 = 1. This non-perturbative feature of the exact solution in the intermediate coupling range is not completely unexpected. Indeed, it resembles what happens in the N = 4 SYM theory where branch-point singularities are present on the negative real axis of the 't Hooft coupling, with the leading one located at λ = −π 2 (see the discussion in [29]). The origin of this singularity is in the structure of the single-magnon dispersion relation, which was derived in the past in perturbation theory [67,68] and using superconformal symmetry [69,70]. Nowadays, it is well-understood in the more modern quantum algebraic curve treatment of the N = 4 SYM theory (see for instance [71]). The fact that the same singularity could also be present in the N = 2 theories considered here was already noticed in the study of mass deformations of N = 4 SYM, where |λ c |/π 2 = 1 is indeed the radius of convergence of the free energy [72]. This sort of universality could be ascribed to the similar structure of the combinatorics of planar diagrams; from a deeper perspective, it might be related to integrability structures yet to be fully clarified [73].

JHEP09(2020)116
We finally remark that in the U(N ) matrix model the operator Ω 1 (a) = tr a is non-zero and that it mixes with all operators of odd dimensions. This fact implies that normalordered version of these operators always contains a component proportional to Ω 1 (a). For example, in the N = 4 U(N ) SYM theory one finds that the normal-ordered single trace operator of dimension 3 at large N , instead of being simply Ω 3 = tr a 3 , is It is interesting to remark that if we compute the two-point function O 3 (a)O 3 (a) in the large-N limit we obtain the same result as in the SU(N ) case for any ν, namely the function γ 3 is still given by the expression given in (3.46) or (4.42). This means that the mixing term (A.9) does not give any contribution in the planar limit. We have explicitly verified that the same thing occurs also for the mixing terms involving Ω 1 (a) in the single trace operators with odd dimensions up to n = 7. These findings confirm the expectation that for the observables that exist in both theories, the SU(N ) and U(N ) models are indistinguishable at large N .

B Chebyshev polynomials of the first and second kind
The Chebyshev polynomials of the first and second kind T n (x) and U n (x) can be defined as

JHEP09(2020)116
They obey the orthogonality relations where I n is the modified Bessel function of the first kind. These relations are easily proved by changing variable to x = cos θ and using the well-known integral representation of the Bessel function I n (a) = 1 π π 0 dθ e a cos θ cos(nθ).

C On the determination of the cut edge µ(λ)
The solution of (4.8) with the Ansatz (4.36) cannot be given in closed form. Nevertheless, there exists a simple iterative scheme that allows us to determine µ 0 in an efficient way. After writing (4.8) in the form The second term in the left hand side above can be expanded in the Chebyshev polynomials of the T type as follows Then the condition (C.2) becomes k ,k T 2k+1 x µ 0 δ k ,k + ν E k, k a 2k = 8π λ T 1 x µ 0 .
(C.4) Let us denote by E the matrix of elements E k ,k ; notice that in the conventions we are using, the index labels start from 0. The solution to (C.4) can then be written as In particular, The coefficients E k ,k can be determined using the orthogonality relation (B.3) in (C.3); they are given by Repeating the analysis discussed after (6.7), we can show that for k > 0 we have E k ,k = 4 (−1) k+k (2k + 1) On the other hand, a 0 is fixed by the normalization condition (4.37) to a 0 = 2/(πµ 2 0 ). This turns the above relation into an equation for µ 0 which can be solved perturbatively, obtaining the result reported in (4.39) of the main text. Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.