Wilson loop correlators at strong coupling in N = 2 quiver gauge theories

We consider 4-dimensional N = 2 superconformal quiver theories with SU ( N ) M gauge group and bi-fundamental matter and we evaluate correlation functions of n coincident Wilson loops in the planar limit of the theory. Exploiting specific untwisted/twisted combinations of these operators and using supersymmetric localization, we are able to resum the whole perturbative expansion and find exact expressions for these correlators that are valid for all values of the ’t Hooft coupling. Moreover, we analytically derive the leading strong coupling behaviour of the correlators, showing that they obey a remarkable simple rule. Our analysis is complemented by numerical checks based on a Pad´e resummation of the perturbative series.

1 Introduction The study of the strong coupling regime in four dimensional gauge theories is a very challenging subject to tackle, yet, over the years, there have been considerable developments, especially in theories with a high amount of supersymmetry.Specifically many results have been achieved in the planar limit of the maximally supersymmetric theory in four dimensions, i.e.N = 4 Super Yang-Mills (SYM), by using supersymmetric localization, holography and integrability.However, more recently, mainly by exploiting supersymmetric localization techniques [1] also the strong coupling regime of some N = 2 gauge theories in the large-N limit has been investigated.In the context of theories with eight supercharges many BPS observables have been analysed using localization techniques both at weak and strong coupling, such as correlation functions between chiral/anti-chiral scalar operators [2][3][4][5][6][7][8][9][10][11][12], correlators between a chiral operator and a Wilson loop [13][14][15][16][17], the vacuum expectation value of BPS Wilson loops [18][19][20] and the free energy [21][22][23][24].
A particular class of N = 2 theories whose properties deserve to be analysed are the quiver gauge theories obtained by a Z M orbifold projection from N = 4 SYM.These are superconformal quiver theories with SU (N ) M gauge group and bi-fundamental matter hypermultiplets.Furthermore, they have a holographic dual realized by Type II B string theory on AdS 5 × S 5 /Z M space-time [25,26], thus they represent one of the simplest models in which to investigate the holographic correspondence when supersymmetry is not maximal.In the framework of these quiver theories many exact results for different correlation functions in the large-N limit have been found over the last few years [27][28][29][30][31][32].
Moreover these quiver gauge theories are the starting point for the construction of other theories by means of orientifold projections (see e.g.[33]).For example from the M = 2 quiver one obtains the so-called E-theory, which has SU (N ) gauge group and matter in the symmetric and antisymmetric representations, that establishes another very fruitful playground where to find exact results and explore the strong coupling regime of gauge theories (see e.g.[3,10,12,17]).
Most of the strong coupling results obtained for BPS-observables in the quiver gauge theories have been carried out at the orbifold fixed point, the most symmetric configuration where all the Yang-Mills couplings associated to the nodes of the quiver are taken to be equal, i.e. g I ≡ g.This way we introduce a unique 't Hooft coupling, namely λ ≡ g 2 N .
Indeed, at the orbifold fixed point, a very useful tool that allows to get exact results for every value of the 't Hooft coupling in the planar limit of the theory becomes available.This is the X-matrix, a semi-infinite matrix whose elements are a convolution of Bessel functions of the first kind.As a matter of fact, as shown in [11,28,29,31], it is possible to express the partition function, and also 2-and 3-point functions of chiral scalar operators, in terms of this semi-infinite matrix.Furthermore, in recent papers [11,31] the features of the X-matrix have been employed to work out a systematic strong coupling expansion in inverse powers of the't Hooft coupling.
In this article we focus on the study of correlators among coincident Wilson loops in the fundamental representation in Z M quiver gauge theories.Previous works on the same topic [15,[34][35][36] considered the more general configuration where all the couplings are different.Although very interesting, this general set up makes more complicated the analysis of the theory at strong coupling.As a matter of fact the authors of [34][35][36] only focus on the v.e.v. of a single Wilson loop.Instead, one of the main novelty of this work is that we only concentrate on the orbifold fixed point of the theory where, as mentioned above, the strong coupling regime of the theory can be successfully investigated using supersymmetric localization and the techniques of [11,31].Moreover, following the analysis carried out in [27] for chiral/anti-chiral operators, it turns out to be convenient to introduce a new basis of Wilson loop operators, namely, generalizing [15,34], we consider untwisted and twisted combinations of the Wilson loops associated to the different nodes of the quiver.
Specifically in this paper we derive an exact expression, valid for every value of the 't Hooft coupling, for the n-point Wilson loop correlator in the planar limit.Then, generalizing to the case of 4d N = 2 Z M quiver theory the techniques introduced in [37][38][39], we find the leading order of its strong coupling expansion, which agrees with the corresponding result in N = 4 SYM up to a numerical proportionality factor.To the best of our knowledge, this is the first example of an exact expression for a generic correlator among n coincident Wilson loops in the planar limit of a four dimensional N = 2 superconformal gauge theory for which also the strong coupling behaviour is analytically derived.This paper is organized as follows: in Section 2 we review the main features of the 4d N = 2 quiver gauge theory and we introduce the untwisted and twisted Wilson loop operators.In Section 3 we briefly recall the relevant aspects of supersymmetric localization and of the corresponding interacting matrix model, furthermore we review the notion of reducible correlator in the large-N limit, that will be extensively used in the following sections.Then, as a warm-up example, in Section 4 we consider the first non-trivial correlator between twisted Wilson loops, i.e. the 2-point function.We formally derive an exact expression for this correlator in the planar limit and we determine the leading term of its strong coupling expansion.For the simplest Z 2 quiver we cross check our analytic prediction through a Padé resummation of the perturbative series.In Section 5 we carry out a similar analysis for the 3-point function.Then in Section 6 we extend our results to the most generic reducible n-point Wilson loop correlator, determining its large-N expression and finding the leading term of its strong coupling expansion.Finally we wrap up in Section 7 with some closing remarks.More technical details concerning the properties of the Bessel functions and the derivation of the strong coupling expansion are collected in three appendices.

The 4d N = quiver gauge theory
In this section we introduce the theory and the observables that we are going to consider.Along the way we also fix the notation.
In this article we focus on the 4d N = 2 quiver gauge theory obtained by a Z M orbifold of N = 4 SYM with gauge group G = SU (N M ).Its matter content is summarized by the quiver diagram reported in Figure 1, where to each line between adjacent nodes is associated an N = 2 hypermultiplet transforming in the bi-fundamental representation (N, N) and, on the other hand, to each node of the quiver is associated an N = 2 vector multiplet transforming in the adjoint representation of the corresponding SU (N ) gauge group.In total there are M nodes labeled by the index I = 0, . . ., M −1. 1 This theory is conformal for arbitrary values of the coupling constants g I , since there are 2N fundamental flavours at each SU (N ) node.Henceforth we only focus on the orbifold fixed point configuration, where all the couplings are taken to be equal, namely Furthermore this theory has a non-trivial R-symmetry group given by SU (2) × U (1) R .
The gauge invariant operator that we are going to consider throughout this paper is the half-BPS Wilson loop along a circle C of radius R in the fundamental representation of the I-th node of the quiver [40,41], namely where P denotes the path-ordering, while A I µ and φ I the gauge field and the complex scalar of the I-th N = 2 vector multiplet.Finally, without any loss of generality, we place the circle C inside the R 2 ⊂ R 4 parameterized by (x 1 , x 2 ).Then the explicit expression of the function x µ (τ ) reads In this work we only consider circular Wilson loops with unitary radius, i.e. we set R = 1.
Starting from (2.2) we define where α = 1, . . ., M − 1 and ρ denotes the M -th root square of the identity, namely The operator (2.4a) is called untwisted Wilson loop, while the operators (2.4b) are called twisted Wilson loops.Of course the two definitions (2.4) can be joined writing where α = 0, 1, . . ., M − 1. 2 The linear combinations (2.6) turn out be very useful for the study of correlators among multiple coincident Wilson loops at the orbifold fixed point, that is 1 All over the paper the index I is taken modulo M , that is to say I ∼ I + M . 2 We recall that the notion of twisted/untwisted Wilson loop is not new.For instance the v.e.v. of a single Wilson loop at the orbifold fixed point was considered also in [15,34], where it was shown that Our goal is to examine the large-N limit of (2.8).Here we just notice that in order to get a non-trivial result one has to require Therefore, without loss of generality, we can always assume that Before addressing the most general case, we observe that the correlator of n coincident untwisted Wilson loops in the Z M quiver at the orbifold fixed point is planar equivalent to its counterpart in N = 4 SYM theory, namely in the 't Hooft limit it holds3 This expression can be regarded as an immediate extension of the results of [15,34,35].Since these correlators are planar equivalent to N = 4, in this work we focus on considering the cases of (2.8) with all α i ̸ = 0 or mixed correlators with both twisted and untwisted Wilson loops.
In the next section we explain how the computation of the large-N limit of (2.8) can be efficiently performed by exploiting supersymmetric localization.

The matrix model
Using supersymmetric localization we can recast the computation of the expectation value of an observable on R 4 in the evaluation of a finite dimensional matrix integral on the sphere S 4 .In this section we review the main properties of the matrix model for the 4d N = 2 quiver gauge theory introduced in Section 2 and we explain how, exploiting it, we can obtain exact expressions for the correlators (2.8) valid in the large-N limit and for any value of the 't Hooft coupling.

Explicit realization of the matrix model
Once the theory is placed on a 4-sphere S 4 of unitary radius, its partition function Z localizes and it can be written as a multiple integral over a set of M matrices a I taking values in the su(N ) I Lie-algebra where Z 1−loop contains the contributions due to the 1-loop determinants of the fluctuations around the localization locus, while Z inst encodes non-perturbative instanton corrections.However, since in the large-N limit instantons are exponentially suppressed, we can ignore this contribution and henceforth we set Z inst = 1.Here we use the so called "Full-Lie algebra approach" (introduced in [42]) and therefore the integrations in (3.1) are performed over all the elements of the matrices a I .We expand these matrices over the basis of the su(N ) I generators {T b } writing where b, c = 1, . . ., N 2 − 1.This way the integration measure in (3.1) can be written as where we choose the normalization factor such that Gaussian integration over each a I is equal to 1. Finally we observe that the 1-loop contribution can be rewritten in terms of an interaction action as where the explicit expression of S int is given by [28] with ζ 2m−1 denoting the odd Riemann ζ-values.Thus the vacuum expectation value of any observable f (a) can be written as where ⟨•⟩ 0 stands for the v.e.v. in the free matrix model.The expectation values in the Gaussian matrix model can be efficiently evaluated by exploiting the set of identities (see equation (2.43) of [3]) satisfied by where a denotes a generic matrix in the su(N ) Lie-algebra.Now we want to show how we can express correlators of untwisted and twisted Wilson loops in the matrix model representation.In order to do that, we introduce the operators with the understanding that Now, exploiting the definition of these operators, we can introduce the matrix model representation of the 1/2 BPS circular untwisted and twisted Wilson loop of unit radius in the fundamental representation, namely [1] It is worth mentioning that the matrix model representation of the Wilson loop in N = 4 SYM is Therefore, through relation (3.10), we connect the computation of Wilson loop correlators to that of correlation functions involving A α,k 's operators.From the definition (3.8) in general we have in the free matrix model where n i=1 α i = 0 mod M .Specifically, one case that will be relevant in the following sections is given by the subset of (3.12) with α i = α ∀ i = 1, . . ., M .We find where t c k 1 ,k 2 ,...,k M is the connected correlator 4 in the Gaussian matrix model, i.e. in N = 4 SYM.For instance, for connected 2-and 3-point correlation functions in the free matrix model we have where we have used the definition (3.7).However, more importantly, we are interested in the large-N behaviour of Wilson loop correlators and, for this purpose, as seen in [30], it is convenient to introduce the vevless basis of the A α,k 's operators, i.e.
As it has been shown in [28,30], it holds that 4 The generating functional of the connected correlators is the free energy F = −log ZN =4, where ZN =4 is the N = 4 partition function deformed by adding the sources {h k i ,I } for all the trace operators and is defined as From the free energy is then possible to obtain the connected n-point functions as Then the leading order of the large-N expansion of a generic higher point correlator can be factorized à la Wick among the product of 2-point functions (for correlators involving an even number of Âα,k 's) or the product of just one 3-point function and 2-point functions (for correlators involving an odd number of Âα,k 's).For example for a 4-point correlator we find where k 1 + k 2 + k 3 + k 4 is even.We notice that in general the coefficient of the leading term of the large-N expansion could vanish and then the correlator would scale as this turns out to be the case, the corresponding correlator is called irreducible since it cannot be factorized into the product of 2-point correlators, otherwise is called reducible.An example of irreducible correlator is provided by the 4-point function

The P-basis and its properties
Although correlation functions of Wilson loops can be naturally expressed using the Âα,k basis, in the interacting theory it is useful to use a new basis, the so called P α,k basis, that enormously simplifies the computations.This basis has been introduced in [30] and here we just review its main properties.In the large-N limit the A α,n basis and the P α,n basis are related as follows 5   Âα,k The interaction action in (3.5) has a remarkable simple form in the P-basis where P α is an infinite vector with entries P α,k , and X is a semi-infinite symmetric matrix, whose entries with opposite parity vanish 5 In principle, in the interacting theory one could define the Âα,n operators as i.e. subtracting the v.e.v of the operator in the interacting matrix model.However since our goal is to study the correlators in the planar limit of the theory, from the (3.5) and exploiting the properties of the v.e.v. of the traces proved in [12,28], one can show that ⟨Aα,n⟩ ≃ ⟨Aα,n⟩0 . (3.23) while its non-trivial elements read where k, ℓ ≥ 2.
Using the above expression for the interacting action (3.25) one can show that the interacting 2-point function reads [30] where We notice that for the untwisted sector (α = 0) D k,ℓ reduces to δ k,ℓ .The interacting 3-point function reads where Higher point functions in the interacting theory can be still computed using Wick's theorem.This is due to the fact that, at the leading term of the large-N expansion, the P α,k operators inherit the properties (3.20a)-(3.20b)valid for the Âα,k in the free theory.This way a correlation function among an even number of P α,k can be decomposed into the sum over all the possible Wick's contractions using the propagator (3.29).For example it holds that On the other hand a correlator among an odd number of P α,k can be decomposed into the sum over all the possible Wick's contractions using the 2-point functions (3.29) and just one 3-point function (3.31).For example We stress that the two relations (3.33)- (3.34) are valid in the large-N limit and are exact in λ.
4 The 2-point function Let us begin the computation of Wilson loop correlators in the Z M quiver gauge theory in the 't Hooft limit, starting from the simplest case, namely the 2-point function.The case of a 2-point function with one untwisted and one twisted Wilson loops is vanishing due to the Z M symmetry of the quiver; we focus on the case of a twisted 2-point function (α ̸ = 0).Using the expression (3.10) we get Since we are considering α ̸ = 0, we can replace the A α,n operators with their corresponding Âα,n ones and then we can exploit the relation (3.24) and (3.29) to write the leading term of the large-N expansion of the correlator as At this stage we consider the following redefinition of the indices therefore in the planar limit the expression (4.1) becomes where we have used the expansion of the modified Bessel functions of the first kind Equation (4.4) is the expression for the 2-point correlator valid for each value of the 't Hooft coupling in the planar limit of the quiver Z M .We notice that (4.4) is exponentially divergent for λ >> 1, for this reason we find convenient to consider its ratio with the corresponding observable in N = 4 SYM.This latter is given by where t c k,ℓ is defined in (3.16).Using the expression for the N = 4 Wilson loop in (3.11), we can rewrite this result as where W conn (λ) is the connected Wilson loop 2-point function that was previously considered also in [43].In the large-N expansion one obtains6 Therefore, we finally consider the ratio (4.9) We observe that when we turn off the interaction action S int the ratio (4.9) is just equal to 1.In the following, we aim to study the leading order of the large-λ expansion of ∆w (α) (M, λ) which, in the planar limit and for any value of the 't Hooft coupling, reads7

Strong coupling limit
In order to determine the strong coupling limit of (4.10) , we follow the same procedure applied in [17].As a first step, we divide the sums over k and ℓ considering separately the odd and even contributions.Then we expand the ratios between Bessel functions in (4.10) for large values of λ using polynomials Q (j) odd 2s (k) and Q (j) even 2s (k), that read Moreover we introduce the functions that, as explained in Appendix B, allows us to express the non-trivial matrix elements of the X-matrix (3.28) as where X denotes the operator acting on the Hilbert space. 8Hence, the strong coupling expansion of the expression (4.10) becomes where the coefficients are given by with and Based on the expansions for Φ (1,2) odd L (x) and Φ (1,2) even L (x) collected in Appendix C, we find that (4.20d) 8 We refer the reader to [31] for the definition of this quantity in terms of the Bessel operators At this point it is useful to introduce the coefficients with n, m ≥ 0 and where ϕ (ℓ) (x) ≡ J ℓ ( √ x) with ℓ = 1, 2. These coefficients are the generalization for every s α of the w n,m introduced in [39] and their main properties have been collected in Appendix B. Crucially for us, we can express the quantities S (P ) (s α ) appearing in (4.15) using the coefficients (4.21) and, for instance, for the first values of P we get S (1) 2,0 (s α )) + (w 2,0 (s α )) 3,0 (s α ) + w 3,0 (s α ) + w As shown in Appendix B, at strong coupling the leading term of the coefficients w n,m (s α ) scales as Moreover, we observe that the leading order contribution takes the generic form for P = n + m.Thus, since the leading term of w n,m (s α ) does not depend on ℓ, henceforth we can omit the dependence on ℓ and just write ω n,m (s α ).Based on these observations, we are now able to determine the leading order (LO) of the expansion of (4.15).We obtain where the overall factor of 2 is due to the fact that the coefficients ω n,m (s α ) are equal for both S (n,m) odd (s α ) and S (n,m) even (s α ).Therefore at the leading order we analytically find The generating function G (0) (s α , x, y) for the coefficients ω n,m (s α ) is defined in Appendix B, where it is also argued that where the integral I 0 (s α ) has been defined in (B.19).Therefore we find that This our final expression for the leading term of the large-λ expansion of the ratio (4.9).For example for the Z 2 quiver gauge theory s α = 1 and the expression (4.31) reads 1 + ∆w (1) (2, λ) 4.1.1Numerical checks for the Z 2 quiver gauge theory Here we provide some numerical checks for the Z 2 quiver gauge theory.We follow the same procedure described in Section 4 of [17].As a first step, using the perturbative expansion of the elements of the X-matrix (3.28), we have generated very long series for the coefficients D k,ℓ with k, ℓ ≥ 2. For example the first orders of the series expansion for where the dots stand for higher orders of the expansion.Then, using these expressions, we have obtained the perturbative expansion of the λ-dependent part of (4.10), namely where the C (p) Z 2 are numerical coefficients and the summation stops at some finite cut-off W .In turns this implies that we need to compute the expansion of the coefficients D (1) k,ℓ up to the order λ W .We choose to fix W = 110; this value represents a good compromise among the need to generate enough precise numerical results and the related computational cost.Using the ratio test we observe that the series (4.34) has a finite radius of convergence located at λ ≃ π 2 .Nevertheless we can extend it beyond this limit with a Padé resummation.For this reason we consider the diagonal Padé approximant As a last step, we find useful to consider a conformal Padé, that has the advantage to be very stable even for very high values of the coupling λ [20,44,45].Following [10] we perform the replacement then we construct the Padé approximant in the z variable inside the unit circle |z| ≤ 1. Finally we express the result as a function of λ using the inverse map The result of this analysis is reported in Figure 3.We observe that for very large values of λ the Padé curve tends to the constant value predicted by (4.32).We regard this numerical result as a strong confirmation of our theoretical strong coupling prediction.

The 3-point function
In Z M quiver gauge theories with M > 2 it appears a new irreducible correlator among twisted Wilson loops, namely the 3-point function In this section we first study the main properties of the 3-point correlator of twisted Wilson loops in the planar limit, focusing on its strong coupling regime, and then we will examine the mixed 3-point function with one untwisted and two twisted Wilson loops.Using (3.10) and (3.19), we obtain We utilise the change of basis (3.24) and in the planar limit we obtain where the Kronecker delta function is needed to enforce that, as discussed in [30], k 1 + k 2 + k 3 must be even.Then we use (3.31) and we write where d k was defined in (3.32)This way the planar limit of the correlator (5.2) factorizes in the product of three contributions, namely where we introduced Now we perform the sums over k p .We notice that it is convenient to treat separately the case of even k p and odd k p .For k p = 2k we have where in the second line we used the properties of the Bessel functions collected in Appendix A.
Similarly, for k p = 2k + 1, we obtain 2k+1,2ℓ+1 . (5.8) Finally using both (5.7) and (5.8) we get our final expression for the planar limit of the correlator (5.2) where the set of permutations reads (5.10) The structure of the correlator (5.9) can be efficiently summarized using the diagram reported in Figure 4. Also in this case it is useful to understand which is the corresponding observable in N = 4 SYM.For this reason we turn off the interaction action and we set M = 1, then the expression (5.9) becomes that corresponds to the planar limit of the N = 4 connected Wilson loop (5.12) Therefore we consider the ratio between (5.9) and (5.12), namely which is equal to 1 for S int → 0 and all the dependence on the coupling is encoded in the function ∆w (α,β) (M, λ).We find convenient to divide and multiply the r.h.s. of the above expression by Figure 4: Graphical representation of one of the "building blocks" of the 3-point function (5.9).To each line is associated a factor , while to the vertex (white circle) is associated the factor d Therefore, in the following, we evaluate the strong coupling behaviour of the expression ) . (5.14)

Strong coupling limit
In order to analyse the strong coupling limit of the ratio (5.14), we use the large-λ expansions (4.11) between ratios of Bessel functions and we introduce the following quantities where where Φ (1) odd P (x) and Φ (1) even P (x) have been defined by the expressions (4.19a)-(4.19b)and ϕ (1,2) (x) ≡ J 1,2 ( √ x).We can express both (5.16) and (5.17) in terms of the coefficients (4.21), this way, using the relations (4.20a), (4.20c) and (4.26), we determine the leading order (LO) of both S (5.19) In appendix B we argue that This way we determine the leading term of the strong coupling expansion of (5.13), which reads where α p are defined in (5.6).This is our final expression for the strong coupling limit of the ratio (5.13).For example for the Z 3 quiver gauge theory all the s αp = 3/4 and the expression (5.21) reads (5.22)

Mixed 3-point correlator
There is another non-trivial 3-point function which has to be considered, namely the correlator with one untwisted Wilson loop and two Wilson loops in conjugated twisted sectors.Hence we have Using (3.18), this expression is given by the sum of two contributions (5.24) Let us focus on the second contribution on the r.h.s. of (5.24).We set k 1 = 2ℓ (otherwise ⟨A 0,k 1 ⟩ 0 is vanishing) so that we get then, using (4.4), the expression (5.25) becomes in the planar limit Finally we observe that for the first contribution on the r.h.s of (5.24) we can perform the same steps as in the previous subsection, which lead to show that in the planar limit this term scales as O 1 N 4 (see (5.5)).Therefore the term in (5.26) is leading with respect to the other in the planar limit and we conclude that which means that the study of ⟨W 0 W α W † α ⟩ is reduced to the 2-point function that we analysed in Section 4.

Higher point correlators
In this section we perform the computation of higher point correlators.Specifically, we first deal with the case of correlation functions involving only twisted Wilson loops, secondly we conclude with the analysis of mixed correlators, i.e. correlators of both untwisted and twisted Wilson loops.As a warm-up example, we first consider the Z 2 quiver and then we generalize our analysis to the Z M quiver.

Higher correlators in the Z 2 quiver gauge theory
We start considering the case of observables involving only twisted Wilson loops in the Z 2 quiver gauge theory.We observe that a correlator among an odd number of twisted Wilson loops vanishes.On the other hand, in order to evaluate a correlator among an even number of twisted Wilson loops, we have to consider the following correlation function As shown in [30] and reviewed in Section 3, in the planar limit correlators of the form (6.1) are in general reducible and can be rewritten in terms of the product of 2-point correlators.In some cases the coefficient of the leading term of the correlator (6.1) vanishes, thus the corresponding correlator is no longer reducible and it is subleading with respect to the reducible ones.Since we just focus on the planar limit, in the following we only consider contributions arising from reducible Wilson loop correlators.Before addressing the general case, as a pedagogical example, we first examine with some details the case of the 4-point function.

The 4-point function ⟨W
We want to compute where we employed the change of basis (3.24) and we used the decomposition of the 4-point correlator in the planar limit as in (3.33).As in the previous sections, we consider the N = 4 observable that is obtained from (6.2) turning off the interaction action S int , namely Differently from the ⟨W 1 W 1 ⟩ case, we observe that the (6.3) does not correspond to an N = 4 connected correlator.Nevertheless, using the relation (3.21), we obtain where W (2) conn (λ) is given by (4.7).Finally we consider the ratio between ⟨W 1 W 1 W 1 W 1 ⟩ and W (4) (λ), which reads The expression (6.5) can be evaluated at strong coupling using (4.32).We find The study of the n-point correlator can be performed considering the ratio between a correlator involving 2n twisted Wilson loops and the N = 4 observable W (2n) (λ) ≃ (W conn (λ)) n obtained turning off the interaction action, i.e. (6.7) In the strong coupling limit, exploiting (4.32), we easily get the result

Higher correlators in the Z M quiver gauge theory
We extend the results obtained for the Z 2 quiver gauge theory to the most general Z M quiver.Unlike before, now we have to consider also correlators among an odd number of twisted Wilson loops.Since the derivation is slightly different, we split the discussion for correlators with even and odd number of twisted Wilson loops.

Even reducible correlators
Let us start with the case of a correlator among 2n Wilson loops, namely We use the (3.24) to expand the correlator (6.9) on the P α i ,k i basis.In the large-N limit it holds that We perform the sums over k i and the expression (6.9) in the planar limit becomes We observe that if we turn off the interaction action the expression (6.11) gets where we used the identity (A.5) and the overall numerical factor is given by N even (α 1 , α 2 , . . ., α 2n ) ≡ δ α 1 ,α 2 δ α 3 ,α 4 . . .δ α 2n−1 ,α 2n + other Wick contractions (6.13) and it holds 1 ≤ N even (α 1 , α 2 , . . ., α 2n ) ≤ (2n − 1)!! .This leads us to define the ratio

Odd reducible correlators
On the other hand, in the case of a correlator involving 2n + 1 twisted Wilson loops we have Since, as recalled in Section 3, in the large-N limit it holds that the correlator (6.17) becomes in the planar limit where the set of permutations Q 3 is given by expression (5.10).If we turn off the interaction action, the expression (6.19) becomes where and it counts the number of non-trivial Wick contractions.We consider the ratio where W conn (λ) is defined in (4.7) while W conn (λ) in (5.12).By construction the above ratio, when we turn off the interaction action, is equal to 1.We evaluate it for large values of λ and we find where we used the large-λ behaviour of the 2-and 3-point correlators given in (4.31) and (5.21), respectively.

Mixed correlators
In order to complete the analysis of the Wilson loop correlators in the planar limit of the Z M quiver gauge theory, let us finally consider the following correlation function with n untwisted Wilson loops and m twisted Wilson loops For each of the untwisted Wilson loops we perform the change of basis in (3.18), namely Then, exploiting the properties of the P α,k operators at large-N , generalizing the procedure shown in Section 5.2, it is easy to see that the leading contribution of (6.25) is due only to the term proportional to the product of the t k i .Therefore we find where the correlator involving only twisted Wilson loops is given by the expression (6.15) or the expression (6.24) for m even or odd respectively.

Conclusions
The main result of this paper is the systematic computation of Wilson loop n-point correlators in the planar limit of Z M quiver gauge theories, which led us to find an exact expression for these correlation functions valid for every value of the 't Hooft coupling.In particular, we considered the reducible Wilson loop correlators, since the irreducible ones vanish in the planar limit.It would be interesting to extend our analysis also beyond the leading planar approximation, thus including the discussion of irreducible correlators, but it seems to be very intricate to find a precise structure for these correlation functions in this regime (to this regard the formulas in Appendix E of [12] could be helpful).Non-planar corrections to other observables in the Z 2 quiver theory have been, for instance, recently investigated in [39].Furthermore, we managed to derive the leading order of the strong coupling expansion of these n-point correlators in expressions (6.15) (for even n) and (6.24) (for odd n), which show that the non-trivial contribution of each twisted Wilson loop is captured by a remarkable simple rule in the planar limit for λ → ∞, namely A similar pattern was observed in [30] for the 3-point functions among twisted chiral operators.
In particular, for the Z 2 case (i.e.only one twisted sector), we recover a large-N factorization analogous to the structure exhibited by correlators among Wilson loops of N = 4 SYM in the 't Hooft limit, i.e.
Z 2 theory: which shows that some properties of the planar limit of N = 4 are inherited also by a proper combination of twisted observables of the Z 2 quiver theory.In all our analysis a crucial role has been played by the relation (5.20) that provides an analytic expression for the the generating function G (0) (s α , −2π).From a purely mathematical point of view, this expression has been obtained by extending the analysis performed in [39] to a M -nodes quiver gauge theory.Moreover, since the generating function G (0) (s α , x) is not associated to any specific observable, the relation (5.20) could turn out to be useful also in the evaluation of the strong coupling regime of other type of correlators with operators belonging to conjugated twisted sectors, e.g.⟨W α T † k,α ⟩ (where T k,α denotes the chiral single trace twisted operators of [30]).We leave the study of these quantities for future work.
It would also be interesting to examine the subleading corrections in the coupling in the planar limit of the theory, as already done for 2-and 3-point functions of chiral scalar operators, the v.e.v of a Wilson loop and the free energy [11,31].In order to determine these contributions, the knowledge of the subleading orders of the generating function G(s α , ℓ, x), defined in (B.33), are required.
Finally, we stress the fact that these correlators, analysed in the large-N 't Hooft limit in the strong coupling regime, should be described and, possibly, computed from the dual supergravity theory through the AdS/CFT correspondence.Even though, to the best of our knowledge, a holographic prediction for these correlation functions is presently unavailable, the expressions that we found could be regarded as a QFT-side prediction and, therefore, constitute a natural starting point for a future holographic investigation.We plan to carry out this analysis in the future.
We begin from the first one.We rewrite the modified Bessel function I 2k+1 in terms of a Bessel function of the first kind J 2k+1 .Then, substituting the resulting expression in (A.1), we get At this point we can exploit the identity (A.20a) of [30], obtaining This concludes the proof of (A.1).The proof of (A.2) is analogous with the only difference that we have to use identity (A.20b) of [30].
We now formally derive another useful identity In order to prove it, we firstly consider and then we will take the particular limit As a first step we split the sum over even and odd contributions Now we exploit the following identities which can be proven by considering the analogous identities valid for the Bessel functions of the first kind Thus we are left with Now we use the following identity between the Bessel functions so that the previous expression becomes Thus we equal the arguments of the Bessel functions and we set them to √ λ obtaining B The coefficients w n,m (s α ) and their generating functions In this appendix we derive the strong coupling expansion of the coefficients w n,m (s α ) and some properties of the corresponding generating functions.Finally, using both analytical and numerical techniques, we provide an argument for the validity of the relation (5.20).

B.1 The X-matrix as an operator
We rewrite the matrix elements (3.28) as where we introduced the symbol function Then, with the aim to simplify the notation, we rescale the 't Hooft coupling introducing and we perform the change of variable 4) this way the matrix elements (B.1) become We then introduce the set of functions with k = 1, 2, . . .Starting from them we can construct two orthonormal basis, namely This way we can regard the non trivial elements9 of the X-matrix (B.5) as an integral representation of the operator X acting on the basis (B.7), namely In the following part of this Appendix it will be very useful to think about the X-matrix in a more abstract way, namely as an operator acting among infinite dimensional vector spaces.This in turn will allow to understand some of its properties at strong coupling.We warn the reader that we will employ the same symbol, i.e.X, to denote both the operator and its infinite dimensional representation.

B.2 The coefficients w (ℓ)
n,m (s α ) In this section we aim to study the main properties of the coefficients w n,m (s α ) introduced in (4.21) and whose explicit expression reads where ϕ (ℓ) (x) = J ℓ ( √ x).We observe that (B.9) can be obtained starting from the expression (C.1) of [39] simply by a rescaling of the symbol function (B.2), namely χ(x) → s α χ(x) .
(B.10) Therefore (B.9) constitutes the generalization, valid to each arbitrary value of s α , of the coefficients w nm introduced in [39], i.e. w n,m (1) ≡ w nm .Let us now move to analyse the properties of the coefficients (B.9) valid for any value of the coupling g.Following the same steps performed in Appendix C of [39], one can firstly show that the coefficients w which implies that only a subset of the initial coefficients (B.9) is independent, namely the ones given by the w (ℓ) 0,n (s α ) with n = 0, 2, 4, . . . .As a matter of fact all the others can be obtained by exploiting equation (B.11), for example 0,0 (s α ) , (B.12) Then one can also show that the coefficients (B.9) satisfy the equations where the function q (ℓ) 0 (z, g, s α ) is a solution of the differential equation subject to the boundary condition at weak coupling q (ℓ) 0 (z, g, s α ) = J ℓ (2gz)+O(g 2ℓ+1 ).We observe that applying several times (B.14b) we can express q (ℓ) n (z, g, s α ) in terms of q (ℓ) 0 (z, g, s α ) and the coefficients w Henceforth we focus on the strong coupling regime and we show how the iterative use of the equations (B.14a) and (B.16) permits to determine the strong coupling expansion of the coefficients w (ℓ) 0,n (s α ) with n = 0, 2, 4, . . . .Following [39] we assume that for large-g they scale as where the index i labels the order of the expansion.Although this analysis could be performed in full generality, for the scope of this article it is enough to determine the leading term of (B.17), i.e. the coefficient w 0,n (s α ).Therefore, in the following, we will mainly focus on this quantity and we will just briefly comment on the computation of the subleading terms of (B.17).
Let us start considering the n = 0 case.The coefficient ω (ℓ) 0,0 (s α ) can be computed applying the same procedure discussed in [37] with the only difference that the symbol χ(z) must be rescaled as in (B.10), this way we find where we introduced the function Then let us consider equation (B.15).A solution to this equation can be found applying the semiclassical methods of [37,38] and, at the leading order of the expansion, can be taken of the form C(z, s α )J ℓ (2gz), where C(z, s α ) is an arbitrary function of z and s α .Then, C(z, s α ) is fixed by evaluating for large-g the integral on the r.h.s. of (B.14a) and demanding that is equal to (B.18).This way, after some mathematical steps, we conclude that Let us now construct the solution q 0 (z, g, s α ) valid for any order of the large-g expansion.Firstly we observe that, for large gz, the Bessel function J ℓ (2gz) admits the following expansion which in turn suggests to consider the following ansatz for the general solution ).This way we construct a solution valid at any order of the strong coupling expansion.Then, the subleading coefficients of the expansion (B.17) can be determined substituting the expression (B.22) in (B.14a) and imposing the equality at each order of the large-g expansion.This leads to Now we consider the cases with n ≥ 1 and we just focus on the leading order of the large-g expansion.The expression (B.22) and the recursion relation (B.16) suggest the following ansatz for the functions q n (z, g, s α ) with n ≥ 1 Then, using (B.22) we determine the coefficients a 1 (z, g) and b 1 (z, g) appearing in the large-g expansion of q 1 (z, g, s α ), finding this way we completely fix the leading order expression of q 1 (z, g, s α ).Then we require that equation (B.14a) is satisfied for the case at hand and we find Let us also consider the n = 2 case, this time (B.16) reads where, using the leading term of equation (B.15), we wrote Similarly to the previous case, we determine the coefficients a 2 (z, g) and b 2 (z, g) and we demand that equation (B.14a) is satisfied.We finally find Applying n number of times this procedure we determine the general expressions The relations (B.32) permit to recursively determine the leading term of the w (ℓ) 0,n (s α ) coefficients, once are known the expressions of the integrals I −n (s α ).
B. 4 The integrals I n (s α ) and the generating function G (0) (s α , −2π) The coefficients ω (ℓ) 0,n (s α ), and consequently the ω n,m (s α ), can be written in terms of the function I n (s α ) defined in (B.19).According to equation (4.46) of [31] the expression of these integrals for n ≥ 1 is where M is the number of nodes of the quiver and ψ (m) denotes the m-th derivative of the digamma function.However, here we need the general expression of the above integrals also for n ≤ 0. For s α = 1 it was found that [39] I n (1) = 2 (−1) n−1 (1 − 2 2−2n )π 1−2n ζ(2n − 1) .(B.40) The case n = 1 require some attention and can be obtained taking the limit of (B.40), namely Then we evaluate it at x = −2π.We collect the results of this computation for the first values of s α in Table 1.As a first check of this numerical analysis we observe that our prediction for s α P [50/50] (G (0) (s α , −2π))   s α = 1 (last row of Table 1) is in remarkable good agreement with the results of [39], where it was conjectured that Furthermore, the numerical results found for s α = 3 4 , 1 2 , 1 4 are in agreement with the theoretical expressions (B.48).These checks permited us to test the accuracy of the Padé resummation (B.51) and strongly suggest that we can rely on this technique also for other values of s α .Therefore we performed a numerical evaluation of the generating function G (0) (s α , x) at x = −2π for all the values of s α up to M = 20, the results of this analysis are reported in Figure 5 and Figure 6.We observe that all the numerical points are compatible with the theoretical prediction within numerical errors.We regard these numerical results as a strong confirmation of our theoretical prediction (B.49).

Figure 2 :
Figure 2: Graphical representation of the two point function (4.4).To each line is associated a factor √ p N I p ( √ λ), while to the vertex (black dot) is associated the term D (α) k,ℓ .The sums over k and ℓ are understood.

B. 25 )
Then the iterative use of the relation (B.16) and of equation (B.14a) permits to determine the expressions of the remaining coefficients ω (ℓ) 0,n (s α ) with n ≥ 1.For example for n = 1 the relation (B.16) reads q (s α ) ones, while the other terms in the sums on the r.h.s. of the expressions (B.33)-(B.34)are the generating functions for the subleading contributions in the expansions of w (ℓ) 0,n (s α ) and w (ℓ) n,m (s α ).

1 GFigure 5 :P
Figure 5: We reported the comparison between the numerical evaluation of the function (B.49) for s α ∈ [0, 1] (black dashed line) and its numerical estimation via a Padé resummation of the perturbative series for fixed values of s α (red open circles) up to M = 20.We observe that all the points lie on the theoretical curve.

Figure 6 :
Figure 6: We reported the ratios (blue open circles) among the Padé points and the function (B.49) for all the different values of s α up to M = 20.Importantly we observe that all these ratios are compatible with the value one within numerical errors.

Table 1 :
Numerical evaluation of the generating function G (0) (s α , −2π) for the first values of s α via a Padé resummation of the corresponding series.