Exact 1/N expansion of Wilson loop correlators in $\mathcal{N}=4$ Super-Yang-Mills theory

Supersymmetric circular Wilson loops in $\mathcal{N}=4$ Super-Yang-Mills theory are discussed starting from their Gaussian matrix model representations. Previous results on the generating functions of Wilson loops are reviewed and extended to the more general case of two different loop contours, which is necessary to discuss coincident loops with opposite orientations. A combinatorial formula representing the connected correlators of multiply wound Wilson loops in terms of the matrix model solution is derived. Two new results are obtained on the expectation value of the circular Wilson loop, the expansion of which into a series in $1/N$ and to all orders in the 't~Hooft coupling $\lambda$ was derived by Drukker and Gross about twenty years ago. The connected correlators of two multiply wound Wilson loops with arbitrary winding numbers are calculated as a series in $1/N$. The coefficient functions are derived not only as power series in $\lambda$, but also to all orders in $\lambda$ by expressing them in terms of the coefficients of the Drukker and Gross series. This provides an efficient way to calculate the $1/N$ series, which can probably be generalized to higher-point correlators.

A useful approach to Wilson loops is to consider suitable generating functions. Quite generally, Wilson loop generating functions are most elegantly formulated in the language of symmetric functions [44], which allows to encode the information on Wilson loops in arbitrary representations of the gauge group and to translate between different sets of basis correlators using combinatorial identities. The generating functions for higher rank Wilson loops introduced in [31] are contained in this language as special cases. In particular, the connected correlators of multiply would Wilson loops turn out to be a natural basis to work with and are a key ingredient in the proof of an interesting involution property [38,40,45,46]. In the case of N = 4 SYM theory with gauge group U(N ), the expression of these correlators in terms of the matrix model solution has been worked out in [38,45,46], and it will be one of the aims of this paper to further elaborate on this relation.
At a time when localization had not been established yet as a theorem [25] and the relation between supersymmetric Wilson loops in N = 4 SYM theory and the Gaussian matrix model had only been conjectured [28], Drukker and Gross [29] calculated the circular Wilson loop expectation value exactly, as a series in 1/N and to all orders in λ. To date, their result remains a rare example of a full series in 1/N that can be obtained from the exact matrix model solution.
The aim of this paper is to develop methods that allow to extract as many results as possible about coincident circular Wilson loops from the exact matrix model expression of the general Wilson loop generating functions. In particular, we shall be interested in the correlators of multiply wound Wilson loops, Tr U k i (1.1) and their connected variants, where the non-zero integers k 1 , k 2 , . . . k h represent the winding numbers. Following Okuyama's slight abuse of nomenclature [38], we will also call such correlators "h-point functions".
The structure of the paper is as follows. We will start by reviewing, in section 2, the general theory of Wilson loop generating functions in the language of symmetric functions [44,46]. This review will end with a generalization to the case of two independent contours, which is necessary for treating coincident circular loops winding with opposite orientations. In section 3, we discuss these Wilson loop generating functions in the case of N = 4 SYM theory, in which they can be obtained exactly by solving a Gaussian matrix model. The result of this discussion will be a general formula for the connected correlators of multiply wound Wilson loops in terms of the traces of symmetrized matrix products, generalizing the results of [38,46]. Section 4 is dedicated to the manipulation of the matrix model result using harmonic oscillator algebra. The deep relation between the Gaussian matrix model and the harmonic oscillator is well known, but the very elegant treatment by Okuyama [47] has gone nearly unnoticed. Therefore, we shall review it here. In section 5, we will consider the one-point functions W (k 1 ) starting with a review of Drukker and Gross' exact series in 1/N . We will be able to add two new results here. First, we find a recursive set of differential equations, from which the series can be constructed and, second, we provide a new derivation of the series in 1/N , which results in an explicit combinatorial formula for the numerical coefficients in this series that have been defined only recursively. These two results find their analogues in the treatment of the connected two-point functions, which is carried out in section 6.
The direct approach will result in the full series in 1/N , but the individual terms are given only as series in λ. It will be checked that the leading term in 1/N coincides with known results. Using the differential equation approach, however, we will show that the series can be constructed from the knowledge of the Drukker and Gross' series for the one-point functions, resulting in expressions to all order in λ. A procedure for this construction, which can be coded on the computer, will be given. Finally, we will conclude in section 7 and add two appendices for technical details.

Generating functions of Wilson loops and correlators 2.1 Brief review of combinatorics and symmetric functions
In this subsection, we will recall some basic combinatorial notions and introduce the symmetric functions. Readers not familiar with them should consult a standard reference such as [48] or the lecture notes [49].
A partition λ ⊢ n is a weakly decreasing (or weakly increasing) set of positive integers λ i The numbers λ i are called the parts of λ, and the number of the parts of λ is denoted by l(λ). Obviously, it holds that l(λ) ≤ |λ|. In order to avoid ambiguities, we will assume l(λ) > 0 throughout the paper, i.e., we exclude the empty partition.
Flipping the diagram along its diagonal defines the transpose partition λ † .
Sometimes, the notation is used, meaning that the integer i is contained in λ a i times. Then, we have (2. 2) The notation (2.1) is particularly useful when relating partitions to permutations. A partition λ is associated with the cycle type of a permutation, if the permutation contains a i cycles of length i.
Thus, λ defines a conjugacy class C λ of the permutation group S n . Defining the centralizer size by we have that |C λ | = |λ|!/z λ is the size of the conjugacy class, i.e., the number of permutations of cycle type λ.
A composition K is a sequence of positive integers, K = (k 1 , k 2 , . . .), which are called the parts of K. The length l(K) is the number of the parts of K, and the weight of K is the sum of its parts, Writing the parts of K in weakly decreasing (or increasing) order uniquely associates K with a partition. Therefore, two sequences of positive integers that differ in the order of their parts define different compositions, while they are considered to define the same partition. As for the partitions, we will assume l(K) > 0.
Given a set X, a set partition ν = (ν 1 , ν 2 , . . .) of X is a sequence of disjoint sets ν i , called the parts of ν, the union of which is X. The length l(ν) is the number of its parts. Of particular interest in this paper are the set partitions of X = [n], where n is some positive integer, and [n] = (1, 2, . . . , n). In particular, given a composition K of weight |K| and length l(K), let ν be a set partition of [l(K)]. This implies that l(ν) ≤ l(K). Then, let K ν be defined as the following composition, Clearly, |K ν | = |K| and l(K ν ) = l(ν).
Let us now introduce the symmetric functions [48,49]. Let e n , h n , and p n be the elementary, complete homogeneous and power-sum polynomials of degree n, respectively. For a ∈ {e, h, p}, given a composition K, we define a K = i a k i . Because a composition K is uniquely associated with a partition λ, we can identify a λ = a K . 1 These functions form bases of symmetric functions (on some countably infinite alphabet). There are three additional classical bases, the monomials, m λ , the Schur basis, s λ , and the "forgotten" basis, f λ . Their role is captured best by considering the Hall inner product, ·, · , or the Cauchy kernel. The monomial basis is the adjoint of the complete homogeneous basis, m λ , h ν = δ λν , the forgotten basis is the adjoint of the elementary basis, f λ , e ν = δ λν , whereas the power-sum basis and the Schur basis satisfy p λ , p ν = z λ δ λν and s λ , s ν = δ λν , respectively. The Schur functions are related to the monomials by the Kostka matrix [48]. 2

Wilson loop generating functions
Generating functions for Wilson loops in arbitrary representations of the gauge group can be formulated elegantly [44,46] using the language of symmetric functions. We will restrict our treatment to unitary gauge groups.
Let U be the holonomy of the gauge connection for a single Wilson loop, an "open" Wilson loop, so to say. We can take U diagonal, U = diag(u 1 , u 2 , . . .) and denote by u = (u 1 , u 2 , . . .) the 1 Most often, these functions are defined with reference to a partition. For commuting variables, the two definitions are clearly equivalent. For non-commuting variables [50], using compositions is more appropriate. 2 The Kostka matrix was used in [32] to obtain the Wilson loops in irreducible representations (Schur basis) from the matrix model solution (monomial basis), but we will not use it here.
alphabet of its eigenvalues. 3 Then, it is obvious that the n-fold multiply-wound Wilson loop is simply given by the power-sum symmetric polynomial of degree n in the eigenvalues. Introduce an alphabet of real parameters y = (y 1 , y 2 , . . .) and define the two generating functions 4 Because these two generating functions contain the same information, we shall work with E(y) in what follows.
Expanding E(y) as a formal power series in the parameters y, one obtains where λ denotes the sum over all non-empty partitions. 5 One may call m λ (u) a monomial representation of the Wilson loop. A Wilson loop in some irreducible representation of U(N ) is given by a Schur function Tr λ (U ) = s λ (u) . (2.8) Using the Cauchy identity [48] in (2.7) one has, Products of multiply wound Wilson loops are given by power-sum functions. More precisely, given a composition K, we have Recall that p K (u) = p λ (u), with λ the partition that is associated with K, because the traces commute in the product. In terms of the power-sum basis, (2.7) reads where ε λ is a shorthand for Wilson loop expectation values . are now encoded in E(y) , and connected correlators . are defined in terms of ln E(y) . For example, expanding ln E(y) in the power-sum basis, yields the connected correlators of multiply wound Wilson loops. Note that the unity term is absent in (2.13).
The above argument can be extended to two or more different Wilson loops. 6 Let us consider two of them, with gauge holonomies U andŪ , respectively. We now need two alphabets of parameters, y and x, and two generating functions, E(y) andĒ(x). Correlators between the two Wilson loops are encoded in the generating function E(y)Ē(x) . Moreover, taking the logarithm defines the connected correlators. For example, 14) The first two terms on the right hand side arise from the unity term in (2.11), and we recall that the sums include only non-empty partitions. Specializing to x = 0 reduces (2.14) to (2.13). Moreover, one may consider the special case of equal Wilson loops, U =Ū , which impliesĒ(x) = E(x). In this case, we simply have E(y)E(x) = E(y ⊕ x), where ⊕ is the operation of alphabet addition [49].
Thus, one would not get any new information from the product.
3 Wilson loop generating functions in N = 4 SYM theory

Review of the single loop case
In this subsection, we briefly review the results of [46]. In N = 4 SYM theory, 1 2 -BPS circular Wilson loops can be mapped by localization to a Gaussian matrix model [28][29][30][31][32]. Our conventions for the matrix model are 6 By different we mean following different contours.
We shall consider only the case of U(N ) as gauge group, in which the matrix integral is over hermitian matrices X. 7 In terms of the Gaussian matrix model, the Wilson loop generating function E(y) defined in (2.6) is given by where g = 1 2 g YM . The 't Hooft coupling is then λ = 4N g 2 . The matrix model integral can be done with standard techniques [52], and the result is [41] where A k represents the N × N matrix [38,41] (A k ) m,n = (A k ) n,m = n! m! e Taking the logarithm of (3.3) yields which tells us, from (2.7), that the traces of the symmetrized matrix products 8 are related to the connected Wilson loop expectation value in the monomial basis, Using purely combinatorial relations, the monomials can be translated into any other basis. For the power-sum basis, which represents the correlators of multiply wound loops, the relation is [49] Tr For SU(N ), the matrices must be hermitian and traceless. 8 Symmetrization includes a normalization factor 1/l(λ)!.
In the above equations, all λ i are strictly positive, because λ always denotes a proper partition.
This implies that all of the loops in the correlator (2.10) have the same orientation. For unitary groups, inverting the orientation of the loop maps each representation to its complex conjugate.
In the above calculation, this amounts to swapping the sign of the gauge coupling, g. Therefore, denoting the generating function of the Wilson loops in the complex conjugate representations bȳ E(y), one has As is obvious from (3.4), we have A −k = A k , the difference arising from the term (kg) m−n .
However, because this term cancels in matrix products, it holds that Tr ). In turn, together with (3.5), this implies Ē (y) = E(y) . This result is of course expected, because the choice of the common orientation of all the loops is irrelevant.
In order to discuss correlators of loops with opposite orientations, we need to consider the case of two different loops discussed at the end of subsection 2.2. This is what we will do next.

Oppositely wound loops
In this subsection, our aim is to generalize the results reviewed in the previous subsection to products of multiply would Wilson loops with arbitrary orientation. The loops are still spacially overlapping to ensure that the configuration remains 1 2 -BPS, so that it can be mapped to the Gaussian matrix model. This is the general case considered by Okuyama [38]. Specifically, we are interested in correlators of the form where k 1 , k 2 , . . . , k h are non-zero integers. 10 The results of the previous subsection can be applied to (3.10), if either all k i are positive, or all k i are negative. These two cases are equivalent, because U −k = (U −1 ) k , but we have seen that E(y) = Ē (y) holds for representations that are complex conjugates of each other. In the general case, we can collect the positive and negative integers into two sets using the commutativity of the traces and rewrite (3.10) as These correlators belong to the generic two-loop case discussed at the end of subsection 2.2.
Consider the two-loop generating function E(y)Ē(x) , where y and x are two independent alphabets of parameters. The matrix model expression of this generating function is There are different ways to proceed from here. The first way is to rewrite (3.12) as One could now take the logarithm in (3.13), which would yield the generating function of the connected correlators (2.14). While this would reproduce nicely the first two terms on the right hand side of (2.14), expanding the remaining term in y and x would seem dreadfully complicated, because of the matrix inverses. Let alone the conversion to the power-sum basis.
10 Zeros lead to trivial modifications, because Tr U 0 = N .
Another way of manipulating (3.12) is to reorder the double sum and collect the terms with equal A k . This yields where the functions α k (y, x) are defined by This time, after taking the logarithm, the expansion in powers of α k /α 0 is straightforward, but the conversion into the power sum basis still seems dreadful. Therefore, we shall proceed differently and follow Okuyama [38].
Okuyama considered the generating function of multiply wound Wilson loops, Evaluating the matrix model expectation value in (3.16) results in where the sum runs over all non-empty subsets of [h], and we defined Then, taking the logarithm in (3.17) yields In (3.19), the terms containing the maximum-rank elementary polynomial e h (y) are precisely those in which (µ 1 , µ 2 , . . . , µ c ) constitutes a set partition of [h] (for the definition of a set partition, see section 2.1), such that each parameter y i appears exactly once in the product. Therefore, where the sum is over all set partitions of [h]. In analogy with (2.4), K ν denotes the set 11 and we introduced Eqn. (3.20) is our final result of this subsection, which generalizes the expressions given by Okuyama [38] to arbitrary h. It is also equivalent to formula (4.18) of [46], which it generalizes to arbitrary integers k i . 12 To end this section, let us consider the special (trivial) case when at least one of the integers k i is zero. We simply have and Here, we cannot call Kν a composition, because the integers ki are not necessarily positive. 12 Formula (4.18) of [46] is expressed in terms of a partition λ and contains explicit symmetrizations both over the kis and over the matrix products. To establish the equivalence with (3.20), one can use the unique association of a partition λ with a given set partition, as explained in [46]. The factor |P λ | stems from the multiplicity of set partitions associated to the same partition λ. The other factor l(λ)! is the normalization factor in the symmetrized product of matrices.
Eqn. (3.24) is a consequence of 1O = 1 O for any operator O, which means that the connected part of the correlator is trivial. To see this explicitly in (3.20), take k 1 = 0 and separate the set partitions into two groups. In the first group, the integer 1 sits alone in a set, in the second group not. Consider first a set partition ν with l(ν) = c parts, which belongs to the first group, i.e., in which one part is ν i = {1}, and the remainder ν ′ = ν/ν i is a set partition of {2, 3, . . . , h} with l(ν ′ ) = c − 1 parts. There are c equivalent choices for i, so that these set partitions contribute to (3.20). Compare this to the contribution of the set partitions ν ′′ of length c − 1 that belong to the second group. These set partitions are obtained by adding the number 1 to one of the c − 1 parts of ν ′ defined above. Their contribution to (3.20) is Thus, the two contributions cancel, which proves (3.24) after iterating through all c.
4 Exact results from the matrix model

Matrix model results from harmonic oscillator
In this section, we will review and elaborate on exact results that can be obtained from the Gaussian matrix model (3.1). Our analysis will be based on the very elegant treatment of Okuyama [47] and exploits the intricate relation between the hermitian matrix model and the simple harmonic oscillator quantum mechanics. Mathematically, this relation appears, because the Vandermonde determinant, which is introduced in the matrix integral by the matrix diagonalization, is most conveniently expressed in terms of Hermite polynomials [52], which represent the eigenfunctions of the harmonic oscillator [38,41]. Therefore, the matrices A k given in (3.4) are nothing but the matrix elements [47] ( where a and a † are the oscillator lowering and raising operators satisfying and the states |i are the normalized eigenstates of the number operator, Before going on, let us slightly change notation by introducing which can be treated as a continuous (real or complex) variable. Because k appears in A k only within z, we will also let A z ≡ A k , so that (4.1) reads Whereas the harmonic oscillator eigenstates are given by i, j = 0, 1, . . . , ∞, the matrix model involves only the elements i, j = 0, 1, . . . , N − 1. 13 The clever insight of Okuyama [47] is that one can work in the infinite-dimensional Hilbert space, if one truncates any sum over the eigenstates by inserting the projector The elegance of this approach can already be seen in the calculation of the matrix element (4.5), which we wish to report here from [47]. One starts with rewriting (4.5) as The relation between the matrix Az and the algebra of the truncated harmonic oscillator was made explicit in [41].
The truncated harmonic oscillator is defined by N × N matrix lowering and raising operators b and b † satisfying where PN−1 is the projector onto the highest eigenstate. This is required by the fact that the trace of any commutator must vanish in a finite-dimensional system, in contrast to the infinite-dimensional system of the standard harmonic oscillator. Nevertheless, the number operator n b = b † b retains the standard commutators After introducing w = z s , the generating function of the Laguerre polynomials [53] can be recognized, so that (4.7) becomes reproducing (4.5).

One-point functions
Although the calculation of the one-point function is very easy by tracing over the matrix (4.8), it is instructive to do the calculation in the infinitedimensional system [47]. One starts with where Tr ∞ denotes the trace in the infinite-dimensional Hilbert space. Then, one exploits the the cyclic property of the trace, as well as the commutators Thus, one can write which reproduces the known result (4.14)

Two-point functions
Let us extend the procedure of the previous subsection to the two-point function Tr(A z 1 A z 2 ). A similar calculation yields [47] ( Furthermore, one can write (abbreviating ∂ n = ∂ zn for n = 1, 2) This is where the calculation stops in [47], and we will take it from there. First, using the sum [53, 8.974.1], (4.16) can be rewritten as which gives Then, using [53, 8.976.4], (4.18) is equal to After reordering the summations, one obtains (4.20) We will further manipulate this expression and integrate it in section 6.

One-point functions
In this section, we revisit the one-point functions W (k) = Tr U k , with k of order unity (as opposed to order N or √ N , for example). Without loss of generality, we can set k = 1, keeping in mind that k appears in the exact result (4.14) only in the combination z 2 = k 2 λ 4N . The general case can be recovered from the case k = 1 by scaling λ → k 2 λ. The 1/N expansion of W (1) was obtained by Drukker and Gross [29]. We will review their solution and provide a very simple check of it. The solution contains certain numerical coefficients, which are defined through a recursive procedure.
Then, we will present an explicit construction, which results in a direct combinatorial formula for these coefficients.

Drukker and Gross' series in 1/N
The expansion into a series in 1/N of the exact result (4.14), where the genus-m contributions W m are given by 14 Here, I α (z) are the modified Bessel functions of the first kind, and the coefficients B(m, k) are determined by the recurrence relation (for m, k > 0) We find it useful to express W m in terms of generalized hypergeometric series, see (B.7). In particular, for m = 0, For m > 0, we have where we have introduced the new coefficients C(m, k) = 4 m+1 (2m + k + 3)B(m + 1, k + 1) . (5.9) The recurrence relation for C(m, k) is easily found from (5.4) and reads  Table 2.
There is a slightly different form of (5.8), which we wish to derive for later purposes. Let us first use a recurrence relation for the modified Bessel function to write (5.3) as Next, consider the second term in the bracket in (5.11). For k = 0, this term does not contribute to the sum (B(m, 0) = 0 for m > 0), nor would it for k = m + 1, so we can safely shift the summation index by one for this term. After this, the two terms can be combined using (5.4) to give The expression (5.12) is valid for all m ≥ 0, but it has one summand more than (5.8) for m > 0.
To the best of my knowledge, the easiest way to check the series of Drukker and Gross (and find it, as one might say with hindsight) is as follows. Consider the exact solution (5.1). Laguerre polynomials satisfy the differential equation [53] x With the help of (5.13) it is straightforward to verify that (5.1) satisfies the differential equation (5.14) Knowing that W (1) has an expansion of the form (5.2), (5.14) implies a differential recurrence relation for the genus-m contributions W m , The leading term W 0 is the unique homogeneous solution of (5.15) (up to a normalization constant) that can be written as a power series in λ.

Explicit construction
In this subsection, we shall present a new explicit construction of W (1) as a series in 1/N . Our immediate aim is to find non-recursive expressions for the coefficients B(m, k) and C(m, k). The calculation will also serve as a blueprint for the analogous calculation in the case of the two-point function, which we consider in the next section.
Let us start with the exact expression (5.1), which can be rewritten in terms of a Whittaker function as [54, 18.11.2] . where we have reordered the double summation, and the sum over k is over all values for which the summand is non-vanishing. 15 We denote by s(n, k) the signed Stirling numbers of the first kind, the unsigned ones being simply (−1) n−k s(n, k).
We remark that an alternative representation of A(n, m) can be found by a similar calculation that starts with the hypergeometric series 2 F 1 (−n, 1 − N ; 2; 2). It yields A(n, m) = 1 (n + 2m)! 2m l=0 2 −l n + 2m l s(n + 2m − l + 1, n + 1) (n + 2m − l + 1)! . (5.31) In this approach, however, the vanishing of the terms with even powers of N is not obvious from the explicit expression and must be checked by other means. Moreover, to establish the equivalence of (5.30) and (5.31), some convolution formula of the Stirling numbers [55] might be employed. In any case, we have verified using computer algebra [56] that the two formulae give the same values.
Another remark is that one can use the recurrence relations of the Stirling numbers and the binomial coefficients to show that the coefficients A(n, m) satisfy the recurrence relation (n + 2m)(n + 2m + 1)A(n, m) = A(n − 1, m) + 1 4 A(n, m − 1) .
At this point, we can make contact with the Drukker-Gross series. Taking (5.12) and substituting the generalized hypergeometric series, one finds (5.33) Because C(m, k) = 0 for k > m, we can extend the sum over k to infinity and reorder the two sums by setting n = k + l. This yields

Connected two-point functions
In this section, we will evaluate the 1/N expansion of the connected two-point functions 17 In subsection 6.1, we will perform an explicit calculation along the same line as we did for the one-point function in subsection 5.2. This will result in a series in 1/N 2 with coefficients that are series in λ and functions of k 1 and k 2 . In subsection 6.2, the leading term will be compared to known expressions from the literature. In subsection 6.3 we develop a new procedure. It will be shown how the connected correlators (6.1) can be constructed from the knowledge of the one-point functions W (k) and develop a procedure by which the 1/N series can be construced. This will be the main new result of the paper.

Explicit construction
Our starting point is the exact expression (4.20), Unfortunately, there does not appear to be an easy way to integrate (6.2) in this form, but we can proceed to expand it as we did with the one-point function in subsection 5.2. First, we use [54, 18.11.2] to express the Laguerre polynomial in terms of a Whittaker function, where we have formally extended the summation over k to ∞, which is safe, because of the binomial coefficient. Then, expanding the Whittaker function into a series [54, 13.14.6], we get where we have introduced the coefficients S(n, k; N ) = 1 (k!) 2 n! N + k 2k + 1 2 F 1 (−n, k + 1 − N ; 2k + 2; 2) . (6.5) 17 We recall our definitions z = kg and Az ≡ A k .
It is helpful to express (6.4) in terms of the variables z ± = 1 2 (z 1 ± z 2 ) , (6.6) and expand it in powers of z 2 + − z 2 − , which gives Equation (6.7) can be readily integrated. The result is The integration constant F (z + ) is determined uniquely by considering the special case z 2 = 0, in which z + = z − and Therefore, F (z + ) = Tr A z 1 +z 2 . Comparing this with (6.1) reveals that the rest of (6.8) represents the connected two-point function, Our next aim is to rewrite (6.10) as a series in 1/N . First, let us return to using z 1 and z 2 , where by ∆ we denote the N -and λ-independent combination 18 We show in appendix A that S(n, k; N ) has the form N 2k+n+1−2m σ(n, k, m) . (6.13) Therefore, substituting z = kg and g 2 = λ 4N into (6.11) yields Then, pulling the sum over m in front, we get where the coefficients are given by Here, we have performed a sequence of sum rearrangements and introduced the coefficients

Leading term
In the case of the leading term, W (k 1 ,k 2 ) 0 , the sum (6.15) simplifies significantly. Using (A.6), the coefficient A(n, j, 0) becomes Then, relabelling p → p − k and exchanging the order of summation gives Substituting (6.17) into (6.15) yields Using hypergeometric function identities, this can be written in several equivalent forms. In particular, The above expressions do not simplify further in terms of generalized hypergeometric series, except for the special cases |k 1 | = |k 2 |. Setting, without loss of generality, |k 1 | = |k 2 | = 1, we have Let us compare these expressions with the results of Beccaria and Tseytlin [42]. They have calculated the genus expansion of the correlators and Tr U Tr In order to find the connected contributions to (6.23) and (6.24), we need to subtract Tr U 2 . From (5.2) and (5.3) we have To square this, we can use the product formula (B.8), in which, for our parameters, the 2 F 3 ()'s simplify to 1 F 2 ()'s. Furthermore, one can use the contiguous function relations of the generalized hypergeometric functions, which we review in appendix B. This gives Thus, after subtracting (6.26) from (6.23) and (6.24) and using again the contiguous function relations of appendix B, one finds (6.21) and (6.22), respectively.
To prove the equivalence with our result, let us first expand the modified Bessel functions in (6.27) into series. After some rearrangement of the two infinite sums one gets Thus, to show that (6.18) is equal to (6.27), we have to establish that n j=0 n j n + 1 j Consider first the left hand side of (6.29). For simplicity, we shall omit the summation limits using the convention to sum over all possible non-zero summands. Using the identity and letting j → n − j in the term with even powers of k 2 in the numerator, we find lhs. = To manipulate the right hand side of (6.29), we first use the hypergeometric function identity [53, 9.137.16], which leads to rhs. = (2∆k 1 k 2 ) n 2 F 1 −n − 1, By means of the quadratic transformation law [53, 9.134.2] and recalling the definition of ∆ (6.12), this is equal to After writing out the hypergeometric series, this becomes Finally, with the help of the identity (6.30) one can show that which is just (6.31). Thus, we have proven (6.29).

Recursive construction
The result (6.15) for the genus-m contribution to the connected two-point function, although exact as a power series in λ, is extremely unwieldy. Beyond the leading order term, a general pattern is not apparent, and operations such as finding the large-λ behaviour would require further work.
Therefore, we shall abandon this explicit approach. Which alternatives do we have for making progress? A look at the one-point function can help. As explained in subsection 5.1, the easiest way to find the genus expansion of the one-point function W (k) is to use the differential equation (5.14) to construct a recursive series of differential equations for W m , (5.15). Although these are second-order differential equations, the physically relevant solutions are unique once the leading order solution W 0 is taken as the start of the recursion. So, the question is whether a similar technique exists for the two-point functions. In this section we shall see that the answer to this question is indeed affirmative.
To start, let us return to the exact expression (4.16), Tr(A z 1 A z 2 ) differs from the connected correlator W (k 1 ,k 2 ) conn by a minus sign and an additional term that depends only on k 1 + k 2 , cf. (6.1). Therefore, (6.34) implies that Let us also recall the exact one-point function (4.14) which satisfies the differential equation Equation (6.37) can be established either by direct calculation or by changing the independent variable in (5.15). One can show by direct comparison with (6.36) and using some Laguerre polynomial identities that (6.35) is nothing but To continue, let us rewrite (6.37) and (6.38) in terms of k, k 1 and k 2 as independent variables, recalling that z = kg = k λ 4N . Therefore, (6.37) becomes Similarly, (6.38) takes the form where, here and henceforth, ∂ n is a shorthand for ∂ n ≡ ∂ kn , andD Next, consider the operator Applying (∂ 1 − ∂ 2 ) from the left yields When acting with this on W (k 1 ) W (k 2 ) , one can use (6.39) to replace the second derivatives, which where we have introduced the new operator Clearly, (6.40) and (6.44) imply that This suggests the following recursive procedure. Let D (k 1 ,k 2 ) n andD (k 1 ,k 2 ) n be operators independent of N and containing at most first derivatives with respect to k 1 or k 2 (the only allowed second derivative is the mixed ∂ 1 ∂ 2 ). They are defined in a recursive fashion by and by fixing the integration constant in D from which we can read off the N −2m coefficient It is reassuring to verify that the leading order term W (k 1 ,k 2 ) 0 is just (6.27).
Let us flesh out this procedure. We start by writing where a n , b ± n , c n , as well asb ± n andc n are functions of k 1 , k 2 and λ. After inserting (6.50) and (6.51) into (6.47) and eliminating the second derivatives by means of (6.39), the terms of order N 0 give rise to the following system of equations, (∂ 1 − ∂ 2 )a n − 3 k 1 a n + 3 k 2 a n − 2b − n = 0 , (6.52a) Moreover, the terms of order N −2 determine the coefficients inD The functions corresponding toD In order to make progress for n > 0, let us introduce the variables and let a n = 1 where the new variablesâ n ,b ± n andĉ n are functions of y and z. The various factors of λ serve the purpose of removing it from the system. Moreover, we note from (6.53a) and (6.53b) that which can be used to form a homogeneous equation from (6.52b) and (6.52c) (for n > 0 only). Then, the system (6.52), with the right hand sides determined by (6.53), can be transformed into the following recursive system for the hatted variables, The case n = 0 is special. In that case, the right hand sides of (6.59) are to be replaced by 0, y 2 8z 2 , y 16z 2 and 1 4 , respectively. The start values are given bŷ The system (6.59) can be solved recursively for n > 0. In each step, one must impose thatâ n , b ± n andĉ n vanish for y = 0, which gives a unique solution and ensures that it corresponds to the connected 2-point function. This essentially implies that one can construct a particular solution of the inhomogeneous equation in the form of polynomials in y. More precisely, one can make a solution ansatz in which the functionsâ n ,b ± n andĉ n are given by y 3 times polynomials of degree n − 1, with coefficients that depend algebraically on z. Then, finding the coefficients amounts to solving a system of linear algebraic equations. This can be easily coded. The solutions until n = 4, obtained with the help of [56], are listed in table 3.
6.4 Special cases |k 1 | = |k 2 | For completeness, we shall provide the explicit expressions for a few subleading terms in the special cases k 1 = k 2 and k 1 = −k 2 . In these cases, it is possible to simplify the general expressions that result from (6.49) by applying the product formula (B.8) and the contiguous function relations listed in appendix B. As before, we can limit the discussion to |k 1 | = |k 2 | = 1, because the general case can be recovered by rescaling λ. Without further details, the first sub-leading terms in the

Conclusions
In this paper, we have discussed various aspects of 1 2 -BPS Wilson loops in N = 4 SYM theory, focusing one exact results that can be obtained starting from the Gaussian matrix model representation. build. Although we have not considered the three-and higher-point correlators in this paper, it is reasonable to believe that the harmonic oscillator formulation carries a lot of potential for further progress also in these cases.
Fourth, we have reviewed the Drukker and Gross expansion of the one-point function W (1) as a series in 1/N 2 and added two new approaches. We have shown that the entire series can be constructed also from a recursive system of differential equations, cf. (5.15). Furthermore, the direct approach, in which the power series in λ is reordered such as to give a series in 1/N 2 , results in several explicit expressions of the numerical (rational) coefficients of the Drukker and Gross series, which were originally defined only in terms of a recursion.
Last, we have considered the 1/N 2 expansion of the connected two-point functions, W (k 1 ,k 2 ) conn , using two different approaches, both of which start from the exact result of the harmonic oscillator approach mentioned above. The direct approach of reordering the series in λ into a series in 1/N 2 results in an exact, although unwieldy, solution to this problem. However, we have also shown how the connected two-point functions, W (k 1 ,k 2 ) conn , are related to the product of one-point functions, W (k 1 ) W (k 2 ) , and constructed a systematic procedure to calculate the series coefficients W m , cf. (6.49). This construction is perhaps the main result of the paper. It is possible that this result is related to other methods that exploit the integrability of the Gaussian matrix model, such as the Toda integrability structure, and it would be very interesting to investigate this. Finally, a generalization of our construction to three-and higher-point functions is left for the future.

Acknowledgements
This work was supported in part by the INFN, research initiative STEFI.
A Explicit forms of S(n, k; N ) and σ(n, k, m) In this appendix, we will derive several forms of the coefficients S(n, k; N ) defined in (6.5). Using a hypergeometric function identity, (6.5) can be also written as Depending on which expression we start with, we shall find different, but non-trivially equivalent, results. This is similar to the expressions (5.30) and (5.31) for the coefficients A(n, m). We shall start by considering (A.1), which will directly show that S(n, k; N ) has an expansion in 1/N 2 .

B Some properties of generalized hypergeometric functions
In this appendix, we will review some properties of generalized hypergeometric series, with special regard to relations between contiguous functions. The main sources of this material are [57] and the earlier [58], as well as [54]. The series (B.3) is convergent for all z if p ≤ q and diverges for all z in the case p > q + 1. In the case p = q + 1, it is convergent for |z| < 1, divergent for |z| > 1 and convergent for |z| = 1 if Re( i b i − i a i ) > 0. If one of the coefficients a i is a negative integer or zero, the series terminates, in which case the above generic statements of divergence or convergence are irrelevant.
If an element of a coincides with an element of b, then this pair of parameters can be omitted. For example, p F q (a, a 2 , . . . , a p ; a, b 2 , . . . , b q ; z) = p−1 F q−1 (a 2 , . . . , a p ; b 2 , . . . , b q ; z) . where by a + 1 we intend that every element of a is increased by unity.
Let θ = z d dz . (B.10) Then, using θz k = kz k one can easily derive that the generalized hypergeometric function satisfies the following differential equation of degree q + 1, (θ + a i ) p F q (a; b; z) = 0 . (B.11) Two generalized hypergeometric functions are said to be contiguous, if their parameters differ by integers. As in the case of the standard hypergeometric functions, there exist a number of linear relations between contiguous functions. The differential equation (B.11), together with (B.9), is an example, but there are others. Following [58], we shall introduce some shorter notation, F = p F q (a; b; z) , (B.12) F(a 1 ±) = p F q (a 1 ± 1, a 2 , . . . ; b; z) , (B.13) F(b 1 ±) = p F q (a; b 1 ± 1, b 2 , . . . ; z) , (B.14) and so on. Then, one has (θ + a i ) F = a i F(a i +) , (B.15a) from which follow the contiguous function relations (a i − a j ) F = a i F(a i +) − a j F(a j +) , (B.16a) Other relations can be found in [57], but will not be used here.