Loop Equation Analysis of the Circular $ \beta $ Ensembles

We construct a hierarchy of loop equations for invariant circular ensembles. These are valid for general classes of potentials and for arbitrary inverse temperatures $ {\rm Re}\,\beta>0 $ and number of eigenvalues $ N $. Using matching arguments for the resolvent functions of linear statistics $ f(\zeta)=(\zeta+z)/(\zeta-z) $ in a particular asymptotic regime, the global regime, we systematically develop the corresponding large $ N $ expansion and apply this solution scheme to the Dyson circular ensemble. Currently we can compute the second resolvent function to ten orders in this expansion and also its general Fourier coefficient or moment $ m_{k} $ to an equivalent length. The leading large $ N $, large $ k $, $ k/N $ fixed form of the moments can be related to the small wave-number expansion of the structure function in the bulk, scaled Dyson circular ensemble, known from earlier work. From the moment expansion we conjecture some exact partial fraction forms for the low $ k $ moments. For all of the forgoing results we have made a comparison with the exactly soluble cases of $ \beta = 1,2,4 $, general $ N $ and even, positive $ \beta $, $ N=2,3 $.


JHEP02(2015)173
1 Introduction Fundamental to random matrix theory is a set of equations known variously as Virasoro constraints, Ward identities, Schwinger-Dyson equation, Pastur equations or loop equations. We will use the latter terminology. These allow, in principle at least, the computation of the large N global scaled asymptotic expansion of correlation functions for eigenvalue probability density functions (PDFs) of the form |x k − x j | β , −∞ < x j < ∞ (j = 1, . . . , N ). (1.1) Here V (x) is referred to as the potential (for Gaussian ensembles V (x) is proportional to x 2 ), while β = 2κ > 0 is sometimes called the Dyson index, with β = 1, 2, 4 corresponding to matrices with orthogonal, unitary and symplectic symmetry respectively (see e.g. [22, chapter 1]). We recall that global scaling refers to a rescaling of the eigenvalues so that their support is a single finite interval (single cut), or a collection of finite intervals (multiple cuts). As a concrete example, consider the Gaussian orthogonal ensemble of real symmetric matrices, defined as the set of matrices of the form G = (X + X T )/2, where X is an N × N matrix with entries independent standard Gaussians. The eigenvalue PDF is given by (1.1) with V (x) = x 2 /2 and β = 1 (see e.g. [22, proposition 1.3.4]). By rescaling λ j → √ 2N λ j , the leading order support of the spectral density ρ (1) (λ; N ) is the interval (−1, 1).

JHEP02(2015)173
Let ρ (1) (λ; N ) denote the one-point function (eigenvalue density) with global scaling, normalised to integrate to unity. The loop equations allow the computation of the large N asymptotic expansion of the resolvent where up to and including terms O(N −6 ) [7,35,44]. The asymptotic expansion of the (smoothed) eigenvalue density follows from the inverse Cauchy transform and gives, with χ λ∈J = 1 for λ ∈ J and χ λ∈J = 0 otherwise, ρ (1) (λ; N ) = 1 π 1 − λ 2 χ λ∈(−1,1) Here the leading term is the celebrated Wigner semi-circle law. Our interest in this paper is in the loop equation formalism for generalised circular ensembles. The latter is the class of eigenvalue PDFs that extend (1.1) from the real line to unit circle, and are thus of the form |e iθ k − e iθ j | β , 0 ≤ θ j < 2π. (1.4) We are motivated by some earlier work of one of the present authors and collaborators [23]. That work relates to the bulk scaled limit of the two-point correlation function for the circular ensemble (1.4) with V (θ) independent of θ. This was first isolated by Dyson [18] in the study of unitary analogues of the Gaussian ensembles. Denoting the PDF by p N (θ 1 , . . . , θ N ), the two-point correlation function ρ (2) is specified by
To deduce (1.7), it was assumed that the small k expansion of the structure function is of the form πβ |k| S(k; β) = 1 + ∞ j=1 p j (κ)y j , (1.8) where p j (x) is a polynomial of degree j. Moreover, with f (k; β) := πβ |k| S(k; β), 0 < k < min (2π, πβ), (1.9) and f defined by analytic continuation for k < 0, it is a rigorous result that [23] f (k; β) = f − 2k β ; 4 β . (1.10) This applied to (1.8) requires that the polynomials in (1.8) have the reciprocal property p j (1/x) = (−1) j x −j p j (x). (1.11) Exact results for β → 0, β = 2, β = 4 and first order expansions about β = 2 and β = 4 were then used to determine the independent coefficients in the polynomial up to the highest order possible. We will show in the present paper that the loop equations provide a systematic approach to the generation of the expansion (1.7).
Our key results consist of two parts -a full and complete constructive proof of the hierarchy of loop equations for circular β ensembles in Propositions 3.1 and 3.3, and the application of this system of loop equations to the Dyson circular β ensemble upon specialisation of the forgoing theory. We have not seen this hierarchy written down in the literature and while it has resemblances with the system of loop equations on R it differs in many significant details. This resemblance is taken up in the discussion contained in section 6.
In addition we give the large N expansion of the m k for fixed but arbitrary k < O(N ) in a particular regime, which we call the global regime, in two ways -a direct one relating to exact forms above (see Corr. 4.1) and another through the generating function, the two-point connected resolvent function, or essentially the Riesz-Herglotz transform of the above two-point density where the leading terms are (see proposition 4.5) It is interesting to note the appearance of the Koebe function in our setting as the leading order and universal coefficient in W 2 , (see (4.12)). This function occupies an important role in the theory of univalent functions, [17,24,32], being the unique extremal example of such functions. However it is not clear how considerations arising from geometric function theory have interpretations in the context of the Dyson circular ensembles. We observe that the analytic properties of the circular ensembles differ markedly from Hermitian ensembles, in that the resolvent functions possess convergent expansions and not formal ones. This is related to the fact that under stereographic projection e iθ = 1+ix 1−ix the Dyson circular ensemble is equivalent to the Cauchy β ensemble with weight i.e. the potential has logarithmic growth and is not in the same universality class as say the Gaussian β ensembles.
As seen in the case of Hermitian matrices revised in the 2nd and 3rd paragraphs, and in the summary of some of the results to be derived for the circular ensembles, the loop equation analysis of correlation functions applies to the global scaling regime. In this JHEP02(2015)173 regime, the length scales are effectively macroscopic. Using different methods of analysis, typically based on Jack polynomial theory (see [22], chapter 12), correlations in local regimes on the length scales of the inter-eigenvalue spacings can be probed. References on that topic include [22], chapter 13 and [14,15,33].
The plan of our work is as follows: in section 2 we define the fundamental resolvent functions required in the theory and give some of their analytic and symmetry properties. The hierarchy of loop equations is derived in section 3 for a general class of potentials. A solution scheme to the loop equations is proposed for one of the large N regimes, based upon matching arguments in the decay of the resolvent functions, in section 4 and a solution scheme specialised to the Dyson ensemble is outlined. This is where our main results of the computer algebra calculations are given. As a reference point to the previous sections we augment the well-known results for the two point correlations for β = 1, 2, 4, general N in section 5.1-5.4 and for any even, positive β and low values of N = 2, 3 in section 5.5, and discuss the comparison of these special cases with those of general κ. In the final section of the present paper, section 6, we review earlier work on loop equations for circular ensembles so as to both contrast our contribution, and to put it in context.

Definitions for the general N , β circular ensembles
The unit circle is denoted T = {z ∈ C : |z| = 1}, the open unit disc is D = {z ∈ C : |z| < 1} and its exterior isD = {z ∈ C : |z| > 1}. The total number of particles in the system is N including the number of test particles. On the unit circle the co-ordinates are ζ = e iθ , and thus |ζ| = 1, with arguments θ ∈ [0, 2π]. The inverse temperature is β = 2κ and usually defined on C\{0}. Complex co-ordinates z, z 1 , . . . are generally defined on the Riemann sphere C . The measure dµ is taken to be absolutely continuous on T with density w of the form An example of the class of potentials that can be admitted are those drawn from the class of Laurent polynomials C[ζ, ζ −1 ] with the structure A vast literature studying the simplest case of the above example, M + = M − = 1, in the context of unitary matrix models was initiated in the works [6,26], which were known to arise as a one-plaquette lattice model of 2-D Yang-Mills theory. However, and we wish to emphasis this point, that we admit potentials with a finite number of isolated singularities at z s ∈ D or z s ∈D, and even on T however subject to additional restrictions. Due to the homotopical inequivalence of closed loops on the punctured Riemann sphere to those on the unpunctured sphere it will not be permissible in general to contract the integration contour T to an interval of the real line. Secondly, even when the forgoing contraction is permitted, unless there is additional symmetry (e.g. evenness with respect to θ = arg(ζ)) the projection of the lower and upper arcs onto the JHEP02(2015)173 interval I will lead to two, albeitly related, distinct weights w(x), x ∈ I. Further insight into this issue will be provided in the discussion contained in section 6.
As a minimum requirement on the potential we will henceforth assume the existence of all trigonometric moments of the form Furthermore we will generally require the winding number of w(ζ) about ζ = 0 to vanish however even this can be relaxed within our formalism, after the inclusion of additional boundary terms. Our ensemble is defined simply through the eigenvalue probability density function where the normalisation is specified by Averages of linear statistics of the eigenvalues are defined by with the implied normalisation 1 = 1. Defining ζ 1 = e iθ 1 , ζ 2 = e iθ 2 , the density and the two-point correlation function are given as the latter having been introduced in (1.5). Central to our theory are the resolvent functions which will serve as generating functions for the moments of the eigenvalues by virtue of the interior and exterior geometrical expansions of the Riesz-Herglotz kernel (2.8) In fact averages with this kernel for the linear statistic, while uncommon in applications of the loop equation method, are not novel in studies of unitary matrix models when one recognises that through ζ = e iθ , z = e iφ (see the remarks associated with eq. (3.2) of [36]). The Riesz-Herglotz kernel, or cotangent kernel, is particularly adapted to the circular case for another reason -it appears in the saddle point equations for the eigenvalue probability density functions of the form (see eq. (3.1) of [36]). The first of a sequence of resolvent functions, the Carathéodory function, is defined by Our first definition of a cumulant A 1 · · · A m+1 c for m ≥ 0 is given implicitly in terms of the average A 1 · · · as This definition differs from other authors, such as Mehta [34] by a factor of a sign, but our definition conforms to the more usual statistical conventions, see section 3.12 of [31] or section 15.10, or pg. 186 of [13]. In contrast to (2.10) Mehta's definition eq. (5.1.4) has sign factors. For example in section 5.1.1 of Mehta [34] the connected two-point correlation function is defined as the negative of (1.12). The unconnected resolvent function or moment of the linear statistic (2.8) is defined by whereas the connected resolvent function or cumulant is defined as In particular our study will focus on the second cumulant, which through a simple calculation is related to the first two densities by the integral formula

JHEP02(2015)173
We also require the potential resolvent functions, which are defined in their unconnected form by Q n+1 (z; z 1 , . . . , z n ) and their connected version by P n+1 (z; z 1 , . . . , z n ) In addition to the definition (2.10) the moments and cumulants are related through their formal exponential generating functions by an equivalent definition 11) and the related recursive relation However we will require a more refined recursive relation which properly recognises the arguments of the resolvents. In addition we will generally not assume symmetry in the arguments and therefore preserve their order, so that when combining sets of these we will perform a string concatenation operation, denoted , rather than the set union. Also I\I j will denote the excision of the variables in I j from those of I whilst retaining the original order. We state these generalised results without proof (these follow from the r 1 = · · · = r m+1 = 1 case of eq. (10) of [43]).
Theorem 2.1 ( [43]). Let I = (z 1 , . . . , z l ) and we designate z l+1 to be a distinguished variable. The moments U l and the cumulants W l satisfy the recursive relation, which is a generalisation of (2.12) The analogous result for the potential resolvents Q l and P l is the following recursive relation, where z is the distinguished variable

JHEP02(2015)173
The moments U n and therefore the cumulants W n are sectionally analytic with respect to z 1 , . . . , z n if the variables are strictly z j ∈ D or z j ∈D as one can see from simple bounds on the remainder terms for m ∈ N Thus there are at most 2 n distinct functions for each W n labelled by the string There are a number of trivial identities and properties satisfied by the cumulants (and moments) which we list for subsequent use: (i) re-labelling symmetry 1 ≤ i ≤ n and for all σ ∈ S n (ii) permutation symmetry within the subsets of variables in ∞ and 0 domains respectively

JHEP02(2015)173
3 Loop equations for general N , β circular ensembles with potential In this section we establish the set of loop equations from first principles for a general potential satisfying the assumptions (2.4) and (2.3), and for the parameters N ∈ N and Re(κ) > 0. We will assume these conditions henceforth. Our approach is an adaptation of Aomoto's method [1], which is also detailed in depth in Chapter 4.6 of [22].
Proposition 3.1. Under the above assumptions, z ∈ C * and z / ∈ T, the first Loop Equation is Proof. The Vandermonde determinant is defined in the standard way A key identity under the restriction ζ j = e iθ j , is the analytic re-expression of the squared Let us consider the following definition of J p and the rewriting of this using integration by parts Now we consider the various terms arising from the left-hand side of (3.3). Firstly we compute the derivative of the Vandermonde determinant Using this we next sum the left-hand side of (3.3) over all independent p and find N p=1 (3.5) Continuing we seek to express the terms on the right-hand side of (3.5) in terms of the connected resolvent functions. To this end we note the following averages have such evaluations -starting with , and also find 2z . This latter result gives the first term on the right-hand side of (3.5). Furthermore, for z, z / ∈ T, we compute

JHEP02(2015)173
Now we turn our attention to the third term on the right-hand side of (3.5). From the symmetry of the integral under p ↔ r we deduce The second term on the right-hand side of (3.5) is Lastly we can evaluate the average appearing above as . Such a limit exists given the analyticity of W 1 (z) for z ∈ D. Combining all of these results we arrive at (3.1).
Remark 3.1. Equation (3.1) of proposition 3.1 is directly comparable to the first loop equation of the hermitian matrix models, which can be found in numerous works. Amongst these works we mention eq. (2.13) and (2.17) of [3] in the a → ∞ limit and eq. (2.26) of [7].
Our next objective is to construct the hierarchy of loop equations, of which Proposition 3.1 is just the base or seed equation. To do this we will employ the insertion operator method [2,7] suitably adapted to the unit circle support. We rewrite potential given in (2.2) using the coefficients v k = kt k thus Employing this new parametrisation we define the insertion operator ζ ∈ C ∂ ∂V (ζ) which has the following properties:
Proceeding on with the task of constructing the higher loop equations we establish a number of preliminary Lemmas.
Lemma 3.1. The first resolvent function is given by or recursively with the convention W 0 := log Z N .
Proof. This, the first case (n = 1) of a sequence, is established by the computation Lemma 3.2. Let z ∈ C , and z 1 ∈ C , . . . , z m ∈ C be pair-wise distinct. The unconnected moment U m satisfies the recurrence relation for m ∈ N ∂ ∂V (z) Furthermore, with z ∈ C and distinct from the forgoing variables, the unconnected potential moment Q m+1 satisfies the recurrence relation Proof. Assume that z = z, z 1 , . . . , z n are all pair-wise distinct. Let us define the Riesz-Herglotz kernel sum A(z) := N l=1 ζ l +z ζ l −z , and the divided-difference potential analogue we compute that the action of the insertion operator on its configuration average is the sum of three parts, using (3.10), (3.7) and (3.8), Furthermore, employing (3.6) and (3.9), we compute that (3.14) Now we proceed to compute the action of the insertion operator on the product A 0 (z)A(z 1 ) · · · A(z n ) by applying the forgoing results. First we apply (3.13) to this particular product and note that ∂ ∂V (z) A(z j ) = 0. Next we substitute (3.14) into the appropriate term of the resulting expression and then deduce Both (3.11) and (3.12) now follow as applications of the above relation.
A key result is that the action of the insertion operator on a particular connected resolvent function generates the next connected resolvent function.
Proposition 3.2. Let us take the variables z 1 ∈ C , . . . , z n ∈ C pair-wise distinct. The resolvent functions W n , n ∈ N are computed from the generating function using the relation Proof. To establish this result we will prove it in its recursive form and then appeal to the initial relation (3.10). In order to prove the recursive form we consider the action of the insertion operator using (3.11) in two different ways, firstly in the form Now we compute the left-hand side of the above starting with the recursive moment- In the second step we have used (3.7); in the third (3.11); in the fourth we have noted that the two terms in the summand are just a division of a common term according to whether z l+2 is either in the argument of the W or the U factor; and the final step is a recognition of the sums involved. Upon comparing the two expressions we conclude In addition we require the action of the insertion operator on the potential resolvent functions. Lemma 3.3. Applying the insertion operator to P n gives, for n = 1 Proof. For (3.17) we apply Lemma 3.2 to the case P 1 (z) = A 0 (z) . Using (3.15) and deduce (3.17). In order to prove (3.18) we adopt a similar strategy to that employed in the proof of Proposition 3.2. Consider the action of the insertion operator on Q l+1 in two different ways, firstly in the form (3.12) Now we compute the left-hand side of the above starting with the recursive momentcumulant relation (2.14) (again I = (z 1 , . . . , z l )) in a sequence of steps Upon comparing the two expressions we arrive at (3.18).
Using Proposition 3.2 and Lemma 3.3 we can apply the action of the insertion operator repeatedly to the first Loop Equation.
Proposition 3.3. The second Loop Equation z = z 1 , z, z 1 / ∈ T is given by Let I denote the m-tuple of variables I = (z 1 , z 2 , . . . , z m ) and the string concatenation operation. In the general case the (m + 1)-th Loop Equation for m ≥ 2 is interchanging the z → 0 limit in the resulting expression and simplifying we deduce (3.19).
To prove the generic case (3.20), which applies for m + 1 ≥ 3, we are going to employ an induction argument and utilise all of our previous lemmas. We act on the left-hand side of (m + 1)-th loop equation (3.20) with the insertion operator ∂/∂V (z m+1 ) and note the following mappings of the terms (nowÎ = (z 1 , . . . , z m+1 )) From the fourth and fifth mappings in this list we note that where we recognise the two terms in the intermediate summation as arising from the latter as to whether z m+1 ∈ I j or not. Combining all these and sorting terms into appropriate categories we see that the the resulting expression is precisely the (m + 2)-th loop equation.  This regime requires a careful analysis of the jumps in W 2 (ζ) across the unit circle ζ ∈ T and of the densities on the unit circle which contain terms that are purely oscillatory with phases proportional to N , such as ζ N (in addition to the purely algebraic dependency on N ). This essentially implies a local analysis in the neighbourhood of distinguished or singular points on the unit circle and a new independent variable replacing ζ, depending on the details of the potential. The other regime is when either |ζ| < 1 or |ζ| > 1, i.e. bounded away from the unit circle, and thus ζ N is exponentially suppressed or dominant depending on the situation -we denote this the Global Regime. In this case the moment index k = O(1) is fixed or k = o(N ), and no information about the larger values of k ∼ O(N ) is apriori accessible. This is the only case we will study here. Nonetheless, by taking N, k → ∞ such that k/N is fixed in the resulting expressions, we can reclaim the expansion (1.7). This is consistent with f (k; β) as defined in (1.9) being analytic in k with radius of convergence min(2π, 2β).
For the Circular β Ensemble in the global regime it is possible to use elementary arguments to fix the algebraic growth of the cumulants, which we do in the following proposition. Proof. We will show this for the W l only as the arguments are identical in the case of the P l . For any z ∈ C such that ||z| − 1| > δ and ζ ∈ T we note the following bounds using the triangle inequality |1 − |z|| These bounds apply for all z ∈ C excluded from the annulus {z ∈ C : 1 − δ < |z| < 1 + δ} and thus we do not need to keep track of the configurations of the co-ordinates (z 1 , . . . , z l ). Applying these basic inequalities to the integral definition of U l , we have Therefore the U l have algebraic growth and because of the purely polynomial relationship with the W l (the inverse of (2.10)) the same conclusion can be drawn for them. However in order to refine the large N behaviour of the W l we will make an analysis of (3.20) using balancing arguments. Let us denote the leading order algebraic term by W l = O(N E l ) with the exponent E l . There are five types of terms in (3.20) with distinct exponents: Of the total number of matchings to apply the balancing conditions, the fifth Bell number B 5 = 52, a number are obviously logically inconsistent, such as B and D, of which there are sixteen of these. In addition a further eight are also inconsistent. The single case of no conditions can also be excluded. A further seven cases lead to E l = 0 which is just the original loop equation. A similar set are the eight neutral or fixed cases where E l is l independent however these are not relevant here. The remaining twelve have potential applications. Of these four are ascending E l+1 > E l , four are descending E l+1 < E l and another four are progressive E l+1 ≶ E l , depending on the sign of E 1 , E 2 , or E 2 − 1. In all these twelve cases the l dependence is linear. The descending cases are only of interest here and are: The last two cases are the same for E 1 = 1 and is the solution we are seeking as the others do not ensure the initial instance W 1 = O(N ). Taking E l = 2 − l we now seek the sub-leading term Matching the sub-leading terms from C, D, F j the only solution is δ l = −1, which also means that the remainder terms left over from the leading one come in at this level.  [19], which was further developed in [27], and culminating in the rigorous proof made for the general β case by Borot and Guionnet [4].
We now specialise all of the preceding theory to the Dyson circular ensemble case with V (z) = 0. In this work our focus will be on the two-point correlation function for the Dyson circular β ensemble analytically continued in the complex plane in the parameters β = 2κ and N . From its definition (1.5) one can readily deduce that for N ≥ 2 a (N −2)dimensional integral representation for this correlation function with the well-known form

JHEP02(2015)173
(see eq. (13.32) of [22]), where use has been made of the closed form evaluation of the normalisation as conjectured in Dyson's original paper [18], (see e.g. proposition 4.7.2 of [22]). Because V = 0 and thus P n = 0, n ≥ 1 there is rotational symmetry of the ensemble and the one-particle density is uniform Therefore we have All dependency of the higher n ≥ 2 resolvent functions on angles is via their differences and for n = 2 we denote θ = θ 2 − θ 1 . Let us define the Fourier coefficients of ρ (2)C (θ) through the trigonometric expansion They possess an evenness property m −k = m k . We can see, either from their definition or from the Loop Equations, that W 2 (z, z) = 0 for z ∈ D and z ∈D. The first Loop Equation (3.1) is satisfied by .
which has the solutions W We deduce that W : in general for the case of even orders N −2k , k ≥ 0 we find  : whereas for the odd orders N −2k−1 we have In order to solve this we have to consider the domains of z, z 1 in a particular orderfirstly we choose z ∈ D, z 1 ∈ D (denoted by 0, 0) which allows us to fix ∂ z W (0) 2 (0, z 1 ) = 0, and thus W (0) 2 (z, z 1 ) = 0. Next we consider z ∈D, z 1 ∈ D i.e. (∞, 0) and from the previous derivative evaluation we conclude (4.12) Proceeding we look at the domain 0, ∞, where z ∈ D, z 1 ∈D, and initially compute the derivative at the origin to be ∂ z W (0) 2 (0, z 1 ) = −4/κz 1 . This allows us to solve for W (0) 2 (z, z 1 ) and we obtain the same result as above. The reason why this is the same is because of the symmetry W 2 (z ∈ 0, z 1 ∈ ∞) = W 2 (z −1 ∈ ∞, z −1 1 ∈ 0). Lastly we examine the ∞, ∞ domain, and using the previous derivative evaluation we deduce that W : at the next order, N 0 , one can derive the equation Again this has a unique solution for W : next we come to the generic case at order In this case one solves for W (−k−1) 2 (z, z 1 ).

JHEP02(2015)173
For the third and higher Loop Equations a generic pattern has set in, so we only treat this general case. In addition to the simple and general statements about the initial coefficients we can give three known exact cases.
It is a non-trivial fact that the Koebe solution W (0) 2 satisfies the following functionaldifferential equation for all configurations of z, z 1 , z 2 subject to the initial values Proof. The (m + 1)-th loop equation resolved to the level N 1−m states From the previous theorem we know W and simplifying, we deduce (4.14).
The results of our calculations undertaken using the solution scheme in figure 1 are given, in the case of W 2 to order N −9 , in the following proposition.
Proposition 4.5. Let s 1 , s 2 be the first two elementary symmetric functions of z 1 , z 2 . Furthermore z 1 , z 2 are strictly bounded away from the unit circle. The two-point resolvent function W 2 in the ∞, 0 domain has the expansion as N → ∞ (4.16) Proof. Using a partial fraction decomposition of (4.15) with respect to ζ = z 2 /z 1 and then making the substitutions (ζ−1) −m → (−1) m (m) k k! for m = 2, 3, . . . we directly deduce (4.16).

JHEP02(2015)173
For m 1 + N − 1/κ we find the [1; 2] approximant yields a rational function of N which agrees with all terms in the expansion (4.16). Another signature of this fit is that higher approximants yield an indeterminate situation, i.e. vanishing of all subsequent Hankel determinants. This is (4.18). For m 2 + N − 2/κ we find the same situation in the case of the [3; 4] approximant, and [4; 5] and higher approximants are indeterminate. The result is (4.19). For k ≥ 3 we expect a [5; 6] approximant would be sufficient however we do not have enough terms in the expansion for this. In case the reader may doubt the veracity of the formulae (4.17)-(4.19) one can in fact directly prove these claims. If one takes the second Loop Equation (3.19) with z ∈D and z 1 ∈D then terms with W 2 (z, z 1 ), W 3 (z, z, z 1 ) vanish and W 1 (z) = W 1 (z 1 ) = −N , so that it reduces to 2N However from (4.7) ∂ z 0 W 2 (z 0 , z 1 )| z 0 =0 = −4(m 1 + N )/z 1 and we deduce (4.18).
Remark 4.5. The exact rational forms of the low moments given in proposition 4.6 are the analogue of the polynomial form of the Gaussian β ensemble moments, and these moments have been enumerated for low index on pg. 9 of [16], in eq. (24) of [35] and where the most extensive set is given in eqs. (3.27)-(3.33) and eqs. (A4)-(A7) of [44].
In the context of the circular Dyson ensembles we observe the following duality formulae are valid. Remark 4.6. Analogous dualities for the moments in the Gaussian β ensemble were established using Jack polynomial theory in the study of Dimitriu and Edelman [16], and the corresponding results for the generating functions were given in [44]. found for general β in the preceding section. These special cases also serve as illustrations of some key properties of the two-point resolvents for general β. At the same time we also highlight some of the differences between the exact results and those found within the global expansion regime, which will arise from contributions to m k (N, κ) when k = O(N ) and a failure of analyticity. We should point out that our formulation has differing normalisation conventions to those of Mehta [34], Chapter 10. Let us define [34], eq. (10.1.6) and (10.1.3) In addition we require the definition of the angular derivative [34], eq. (10.3.6) which can be expressed in terms of the first and second Chebyshev polynomials Furthermore we make the definition of the indefinite integral [34], eq.

β = 2 CUE
A standard result gives, see [34] pg. 196 eq. (10.1.13), From its Fourier decomposition the moments are and we note that this is not analytic at k = N . One can readily see that this has an exact large N continuum limit The two-point resolvent function is readily computed and found to be given by This differs from the leading, universal term of (4.15) by the term ζ N in the numerator, which is not accessible in the global regime. It is interesting to observe that W 2 (ζ) satisfies the Bieberbach property |m k + N | ≤ |k| for all k, N .

JHEP02(2015)173
or without loss of generality the partial fraction sum formula valid for k > 0 .
There are a couple of observations to note here -as well as terminating at |k| = 2N −2, m k + N has a maximum at k = N of 1 2 N + 1 4 N ψ(N + 1 2 ) − ψ( 1 2 ) . Thus W 2 , in this case, violates the Bieberbach inequality and fails to be univalent. The large N continuum limit is which exhibits a weak non-analyticity at x = 1.

β = 1 COE
Note that m k is non-zero for all k.
Proof. In this case we need to employ the formulae (5.2) and (5.4) from which we compute or in form of the partial-fraction sum, when k > 0 The moment m k does not have a maximum for finite k but approaches zero as k → ∞. The large N continuum limit of m k is now and is not analytic at x = 1 (very weakly though, as the difference between either side of x = 1 first appears at the third order).
In respect of the moments we can readily verify from (5.8) and (5.5) that they satisfy 11) and that in the global regime we have W 2 (ζ; −2N, 1 2 ) = 4W 2 (ζ; N, 2) as is evident by comparing (5.7) and (5.10). Thus there is consistency with proposition 4.7. However the exact forms (5.6) and (5.9) do not satisfy this latter relation because the symplectic moments terminate whilst the orthogonal ones do not even though the first set of 2N − 2 moments are related by (5.11).
Remark 5.2. In computing the large N expansions of the resolvent functions W 2 for these special β's we employed explicit and elementary function representations of the corresponding densities, without the need of other methods and in particular the use of skeworthogonal polynomials. The asymptotics of skew-orthogonal polynomials has been studied in [20] where one can find the leading order asymptotics for the skew-orthogonal polynomials for a polynomial potential and then applied to the kernels and the two-point correlations, for β = 1, 2, 4. However we have the exact forms from which the large N expansions are readily and systematically constructed to any order of approximation.

JHEP02(2015)173
the first parameter with β = 2M , M ∈ N. Alternatively this can be expressed as a 2k-sum m k = −9δ k,0 + 6(−1) k cos πκ 2 4κ Γ(κ + 1) 3 Γ(3κ + 1)Γ(2κ + 1 + k)Γ(2κ Proof. In order to compute the moments m k we expand the hypergeometric sum, integrate term-by-term and find an integral of the form (5.13), but with the replacements M → M + l/2 and k → k − M + l/2 for some l ∈ Z ≥0 . Resumming this again we have (5.19). This hypergeometric function with unit argument is an integer extension of a terminating Watson's Sum (see eq. (16.4.6) of [41]) for which alternative sums have recently become available -i.e. are 2k-fold sums rather than 2M -fold sums. From eq. (24) of [12] we see that what we seek is W 2k,0 (a, b, c) with a = −2M, b = 1 + β, c = −β (of course β = 2M but we only apply the termination through one parameter initially). Thus we can utilise Theorem 5, pg. 474, of that work with the above specialisations and employing a terminating form of Watson's Sum we arrive at (5.20).
Finally we display some of the low moments for the purposes of checking and comparison

A brief literature survey on loop equations for circular ensembles
In conclusion we have given a self contained and complete proof of a hierarchy of loop equations for circular β ensembles. We did this for the purpose of setting up a formalism to give a systematic derivation of the sequence of degree k polynomials in the coupling κ = β/2 occurring as the coefficients in the small k expansion of the bulk scaled structure function πβS(k; β)/|k| (1.7). To derive the loop equations, our starting point is an adaptation of what in the theory of Selberg integrals (see e.g. [22], chapter 4) is known as Aomoto's method, and in particular we work directly and specifically with the circular ensemble PDF. Our work differs from previous literature relating to loop equations for circular ensembles in its motivation, methodology and technical achievements, as we will indicate by giving a brief survey of some relevant literature. As remarked in section 2, the loop equation formalism for circular ensembles can be traced back to the study of (1.4) in the particular case β = 2 and V (θ) = t cos θ. The corresponding partition function is a special case of the so called Brezin-Gross-Witten where [dU ] is the normalised Haar measure on U (N ). The work [8] deduced that (6.1), upon the replacement g 2 → g 2 N satisfies the Schwinger-Dyson equation where the eigenvalues of J † J are {x j } N j=1 . With the exponent in (6.1) replaced by V = ∞ k=−∞ t k Tr U k , the studies [5,25] deduced the Schwinger-Dyson equations L ± n Z = 0, where L ± n are the Virasoro operators A highlight of this line of investigation, which includes [37][38][39][40]42], was the work of Hisakado [28][29][30]. In the latter, for the potential V (θ) = t cos θ, both the Toda lattice equation and Virasoro constraints were used to characterise the corresponding partition function in terms of a solution of the Painlevé III equation.

JHEP02(2015)173
In the Introduction, we recalled some results relating to the global scaling limit of the Gaussian β-ensembles, which for β = 2 is a particular Hermitian matrix model. In a paper published in 2005, Mizoguchi [36] made use of the Cayley transformation U = (I + iH)/(I − iH) between unitary and Hermitian matrices to initiate a study of unitary matrix integrals with the aim of obtaining genus expansions of the free energy. By way of motivation, he writes: "The recent use of matrix models for the study of gauge theory and string theory requires not only the knowledge of their critical behaviours but also their individual higher genus corrections away from criticality; the technology to compute them has been less developed in unitary one-matrix models than in Hermitian ones." In 2006 Chekhov and Eynard [10] undertook a loop equation analysis of a class of βgeneralised matrix models defined by (1.1), with V (x) analytic in x and the the eigenvalues restricted to a given contour in the complex plane. It is also required that the absolute value signs in the product of differences be removed. Formally at least, this includes a class of circular ensembles. However, the correlation functions are not based on the Riesz-Herglotz kernel (2.8), but rather the resolvent kernel 1/(ζ − z) familiar in the study of Hermitian matrix models. Various extensions of this study are given in [2,9,11]; none treat specifically the circular ensembles nor arrange for the correlation functions to be based on the kernel (2.8). Analytic features particular to circular ensembles (or more generally closed contours), such as the need to consider the domain inside, and the domain outside, the unit circle on equal footing do not show themselves.
Thus, by treating the circular ensemble directly, we have been able to make stronger analytic statements than hold for a β ensemble on a general curve. By way of application, we have been able to provide a computational scheme for the problem at hand, namely the systematic derivation of the polynomials in (1.7), and this in turn has lead to the discovery of some new rational function structures for the moments m k in expansion (4.6) as given in Proposition 4.6.