O Marchal

Asymptotic expansions of some Toeplitz determinants via the topological recursion

In this article, we study the large n asymptotic expansions of n × n Toeplitz determinants whose symbols are indicator functions of unions of arc-intervals of the unit circle. In particular, we use a Hermitian matrix model reformulation of the problem to provide a rigorous derivation of the general form of the large n expansion when the symbol is an indicator function of either a single arc-interval or several arc-intervals with a discrete rotational symmetry. Moreover, we prove that the coefficients in the expansions can be reconstructed, up to some constants, from the Eynard-Orantin topological recursion applied to some explicit spectral curves. In addition, when the symbol is an indicator function of a single arc-interval, we provide the corresponding normalizing constants using a Selberg integral and thus we obtain a generalization of Widom's result.

A Toeplitz integral with symbol

f (e iθ k ) 1≤p<q≤n

where I is a union of d (d ≥ 1) intervals in [-π, π]:

I = d j=1 [α j , β j ] with -π ≤ α j ≤ β j ≤ π (1.2)
Note that in the case of a full support I = [-π, π] it is well known that (See [START_REF] Mehta | Random matrices[END_REF]):

Z n ([-π, π]) = 1 (1.3)
Thus, in the rest of the article, we assume that I = [ -π,π]. Since the integral is obviously invariant under a global angular translation (θ → θ -Cste 1), we may assume that no interval [α j , β j ] 1≤j≤d contains ±π. Integrals of type (1.1) can also be understood as the partition functions of a gas of particles restricted to the set T = {e it , t ∈ I} with 1 olivier.marchal@univ-st-etienne.fr interactions given by 1≤p<q≤n e iθp -e iθq 2 . Moreover, it is also well known that integrals of type (1.1) are Toeplitz integrals that can be reformulated as the determinant of a n × n Toeplitz matrix with a non-vanishing symbol (defined below) on some arc-intervals of the unit circle. Asymptotic expansions of Toeplitz determinants and integrals of the form (1.1) have been studied for a long time and many results already exist in the literature using different strategies. For example, studies using orthogonal polynomials on the unit circle, properties of powers of random unitary matrices, Fredholm determinants, Riemann-Hilbert problems, etc. have been used to tackle the problem. A non-exhaustive list of articles on the subject is [START_REF] Szegö | On certain Hermitian forms associated with the Fourier series of a positive function[END_REF][START_REF] Widom | Strong Szegö limit theorem on circular arcs[END_REF][START_REF] Widom | Strong Szegö limit theorem on circular arcs[END_REF][START_REF] Böttcher | Toeplitz determinants[END_REF][START_REF] Böttcher | Two Elementary Derivations of the Pure Fisher-Hartwig Determinant[END_REF][START_REF] Fisher | Toeplitz determinants, some applications, theorems and conjectures[END_REF][START_REF] Basor | The Fisher-Hartwig Conjecture and Toeplitz Eigenvalues[END_REF][START_REF] Borodin | A Fredholm determinant formula for Toeplitz determinants[END_REF][START_REF] Duits | Powers of large random unitary matrices and Toeplitz determinants[END_REF][START_REF] Deift | A Riemann-Hilbert approach to asymptotic problems arising in the theory of random matrix models, and also in the theory of integrable statistical mechanics[END_REF][START_REF] Krasovsky | Aspects of Toeplitz determinants[END_REF][START_REF] Deift | Asymptotics of Toeplitz, Hankel, and Toeplitz+Hankel determinants with Fisher-Hartwig singularities[END_REF][START_REF] Deift | On the asymptotics of a Toeplitz determinant with singularities[END_REF][START_REF] Deift | Toeplitz matrices and Toeplitz determinants under the impetus of the Ising model. Some history and some recent results[END_REF][START_REF] Duits | Relative Szegö asymptotics for Toeplitz determinants[END_REF]. The purpose of this article is to provide a rigorous large n expansion of Toeplitz integrals of the form (1.1) to all order in 1 n using the Eynard-Orantin topological recursion defined in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]. In particular we shall prove the following results:

• A complete large n expansion when the support is restricted to only one interval (i.e. d = 1) in Section 2.

• Results up to O(1) when the support is composed of d = 2r + 1 ≥ 3 intervals of the form [α j , β j ] = 2πj 2r+1 -π 2r+1 , 2πj 2r+1 + π 2r+1 with -r ≤ j ≤ r and 0 < < 1 in Section 3.3.

• Results up to O(1) when the support is composed of d = 2s ≥ 2 intervals of the form

[α j , β j ] = π(j-1 2 ) s -π 2s , π(j- 1 2 )
s + π 2s with -(s -1) ≤ j ≤ s and 0 < < 1 in Section 3.6.

The strategy used in this article is to reformulate the Toeplitz integrals in terms of some Hermitian matrix integrals (with some restrictions on the eigenvalues support). Then we compute the associated spectral curve and the corresponding limiting eigenvalues density. Using the theory developed in [START_REF] Borot | Asymptotic expansion of β matrix models in the one-cut regime[END_REF][START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF][START_REF] Borot | Large N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF] we are able to rigorously prove the general form of the large n expansions of the correlators and of the partition function, as well as relate them with quantities computed from the topological recursion. In the case of a single interval, we finally use a Selberg integral to fix the normalization issues of the partition function and thus provide the complete large n expansion of the Toeplitz integrals. We eventually compare our theoretical predictions with numeric simulations performed on the Toeplitz determinant reformulation (which is very convenient for numeric computations) up to o 1 n 4 .

Various reformulations of the problem

There are several useful rewritings of the integral (1.1). We list them in the following theorem: 2. The determinant of a n × n Toeplitz matrix:

Z n (I) = det (T p,q = t p-q ) 1≤p,q≤n (1.5) with discrete Fourier coefficients given by:

t 0 = 1 2π d j=1 (β j -α j ) = |I| 2π t k = 1 2π d j=1 e ik α j +β j 2 (β j -α j ) sin c k(β j -α j ) 2 , ∀ k = 0 (1.6)
where we denoted sin c (x) = sin x x the cardinal sine function. 3. A real n-dimensional integral with logarithmic potential and Vandermonde interactions:

Z n (I) = 2 n(n-1) (2π) n n! I n dθ 1 . . . dθ n 1≤p<q≤n sin 2 θ p -θ q 2 = 2 n 2 (2π) n n! J n dt 1 . . . dt n ∆(t 1 , . . . , t n ) 2 e -n n k=1 ln(1+t 2 k ) (1.7)
where ∆(t 1 , . . . , t n ) is the usual Vandermonde determinant ∆(t) = 1≤p<q≤n (t p -t q ) 2 .

4. A Hermitian matrix integral with prescribed eigenvalues support:

Z n (I) = c n Nn(J ) dM n (det(I n + M 2 n )) n (1.8)
where N n is the set of Hermitian matrices with eigenvalues belonging to J . The normalizing constant c n is related to the volume of the unitary group:

c n = 1 (2π) n n! 1 Vol U n = 1 (2π) n n! (n!) n-1 j=1 j! π n(n-1) 2 = 1 2 n π n(n+1) 2
n-1 j=1 j! 5. A complex n-dimensional integral over some segments of the unit circle with Vandermonde interactions:

Z n (I) = (-1)

n(n+1) 2 i n T n du 1 . . . du n ∆(u 1 , . . . , u n ) 2 e -n n k=1 ln u k (1.9)
proof:

The proof of the previous lemma is rather elementary. The reformulation in terms of Toeplitz determinant is standard [START_REF] Widom | Strong Szegö limit theorem on circular arcs[END_REF][START_REF] Böttcher | Toeplitz determinants[END_REF]. Indeed it is well known [START_REF] Böttcher | Toeplitz determinants[END_REF] that for a function (usually called "symbol" in the context of Toeplitz integrals) f measurable on the unit circle we have:

I n (f ) = 1 (2π) n n! I n dθ 1 . . . dθ n n k=1
f (e iθ k ) 1≤p<q≤n e iθp -e iθq 2

= det (T p,q = t p-q ) 1≤p,q≤n with t k = 1 2π

2π 0 f (e iθ )e ikθ dθ , ∀ k ∈ -(n -1), n -1
Thus, equality between (1.4) and (1.5) corresponds to the application of the last identity with f = 1 T . Equality between (1.4) and (1.7) follows from the change of variables θ j = tan

t j 2
which is allowed since the support of the angles is included into (-π, π). With this change of variables we get:

e iθp -e iθq 2 = 4 tan θp 2 -tan θq 2 2 1 + tan 2 θp 2 1 + tan 2 θq 2 Observing that dθ k provides a factor dθ k = 2 1+t 2 k
dt k and that we have:

1≤p<q≤n 1 1 + t 2 p 1 + t 2 q = n k=1 1 (1 + t 2 k ) n-1
immediately gives (1.7). Reformulating the real integral (1.7) in terms of a Hermitian matrix integral is standard (see [START_REF] Mehta | Random matrices[END_REF]) from diagonalization M = U ΛU † of normal matrices. We only note here that the support of eigenvalues is prescribed to J . Eventually the volume of the unitary group can be found in [START_REF] Zhang | Volumes of orthogonal groups and unitary groups[END_REF] and equality between (1.4) and (1.9) is straightforward from the change of variables u j = e iθ j .

As presented in the previous theorem, the complex (1.9) and real (1.7) integral reformulations of the problem share the important point that the interactions are given by a Vandermonde determinant ∆(x) 2 . We stress that this situation is rather exceptional since a (non-affine) change of variables in such integrals does not generally preserve the form of the interactions.

Remark 1.1 We inform the reader that part of the results proven in this article have already been presented in [START_REF]Matrix models, Toeplitz determinants and recurrence times for powers of random unitary matrices[END_REF] using the reformulation in terms of the complex integrals (1.9) in the context of return times for the eigenvalues of a random unitary matrix. Indeed, as one can obviously see from (1.7) and (1.9), the complex and real integral reformulations share the crucial fact that the interactions between the eigenvalues are of Vandermonde type: ∆(x) 2 . This implies that correlation functions of boths reformulations satisfy some loop equations and that the topological recursion may be applied to both models. This leads to the fact that the spectral curves found in this article are related by some symplectic transformations to the ones presented in [START_REF]Matrix models, Toeplitz determinants and recurrence times for powers of random unitary matrices[END_REF] and thus that they provide the same sets of free energies (that reconstruct up to some constants ln Z n (I)). However, we stress that several practical and theoretical issues were disregarded in [START_REF]Matrix models, Toeplitz determinants and recurrence times for powers of random unitary matrices[END_REF] that can be solved using the real integral reformulation:

1. In [START_REF]Matrix models, Toeplitz determinants and recurrence times for powers of random unitary matrices[END_REF], results from [START_REF] Borot | Large N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF] were used to justify the form of the large n expansions of the correlators and partition functions. However, results from [START_REF] Borot | Large N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF] have only been proved for integrals with support included in R, but not in the case of a generic closed curve in C as would require the unit circle. Though it is believed that the tools developed in [START_REF] Borot | Large N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF] should remain valid for certain non-real domains of integration, the mathematical proof is still missing.

2. In [START_REF]Matrix models, Toeplitz determinants and recurrence times for powers of random unitary matrices[END_REF], normalization issues related to the partition function were completely disregarded. Consequently, the large n asymptotic expansions of the Toeplitz determinants presented in [START_REF]Matrix models, Toeplitz determinants and recurrence times for powers of random unitary matrices[END_REF] lack some constant terms that we provide in this paper. Though the normalization issues were not particularly important for the physics problem studied in [START_REF]Matrix models, Toeplitz determinants and recurrence times for powers of random unitary matrices[END_REF], they become essential when one wants to compute exact probabilities and compare them with numeric simulations. It turns out that the reformulation (1.7) offers, at least in the one interval case, a simple way to deal with the normalization issues by connecting them to a Selberg integral. In the complex setting the connection to a known integral is less obvious and thus the normalization issues were disregarded in [START_REF]Matrix models, Toeplitz determinants and recurrence times for powers of random unitary matrices[END_REF].

3. Stronger results than those developed for general interactions [START_REF] Borot | Large N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF] regarding the large n expansion and the reconstruction by the topological recursion are available in the case of real integrals with Vandermonde interactions like (1.7) in [START_REF] Borot | Asymptotic expansion of β matrix models in the one-cut regime[END_REF][START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF][START_REF] Borot | Large N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF]. In particular these results allow a proper rigorous mathematical derivation of the full asymptotic expansion of the Toeplitz integrals (1.1) and provide the explicit expressions of the first orders including the proper normalizing factors.

4. Numerical simulations in [START_REF]Matrix models, Toeplitz determinants and recurrence times for powers of random unitary matrices[END_REF] were only performed to leading order O(n 2 ) while we are able to match the theoretical results with numeric simulations up to O (n -6 ) in this article.

2 Study of the one interval case

Known results

When d = 1, we can use the rotation invariance and take β 1 = -α 1 = π with 0 < < 1. For simplicity we shall denote in this section a = tan π 2 and a function of will equivalently be seen as a function of a and vice-versa depending on the relevance of the parameter in the discussion. In this section, we want to compute the large expansion of the Toeplitz integral:

Z n (a) notation ≡ Z n ( ) = 1 (2π) n n! [-π ,π ] n dθ 1 . . . dθ n 1≤p<q≤n e iθp -e iθq 2 (2.1)
The corresponding Toeplitz determinant reformulation is particularly easy:

Z n ( ) = det (T p,q = t p-q ) 1≤p,q≤n (2.2)
with the discrete Fourier coefficients given by:

t k = sin c (kπ ) , ∀ k ∈ -(n -1), n -1 (2.3)
The reformulation in terms of a n-fold integral corresponding to the diagonalized form of a Hermitian matrix integral is given by:

Z n (a) = 2 n 2 (2π) n n! [-a,a] n dt 1 . . . dt n ∆(t 1 , . . . , t n ) 2 e -n n k=1 ln(1+t 2 k ) (2.4)
Such integrals have been studied by Widom in [START_REF] Widom | Strong Szegö limit theorem on circular arcs[END_REF] and do not fall in the standard theory of Toeplitz integrals developed by Szegö. Indeed the standard case corresponds to a symbol f which is strictly positive and continuous on the unit circle. In that case, the standard theory of Toeplitz determinants can be applied and one would obtain the strong Szegö theorem [START_REF] Szegö | On certain Hermitian forms associated with the Fourier series of a positive function[END_REF]:

1 n ln det (T n (f )) n→∞ → 1 2π 2π 0 ln(f (e iθ ))dθ (2.5)
When the symbol is discontinuous but remains strictly positive, adaptations of the theory, known as Fisher-Hartwig singularities, have been developed and the convergence of

1 n ln det (T n (f )
) is still obtained though formulas get more involved. However, in our present case, the symbol is vanishing on several intervals of the unit circle and the convergence of 1 n ln det (T n (f )) does no longer hold. Indeed, in this case, Widom proved the following proposition [START_REF] Widom | Strong Szegö limit theorem on circular arcs[END_REF]:

Theorem 2.1 (Widom's result) Let 0 < θ 0 < π and define T (θ 0 ) = {e it , t ∈ [-θ 0 , θ 0 ]}.
Then we have:

ln det T n (1 T (θ 0 ) ) = n 2 ln sin θ 0 2 - 1 4 ln n - 1 4 ln cos θ 0 2 +3 ζ (-1) + 1 12 ln 2 + o(1)
where ζ denotes the Riemann zeta function.

We propose in this section to improve Widom's result by providing a mathematical proof of the form of the large n expansion of ln Z n ( ) as well as a general way to compute all sub-leading corrections.

Spectral curve and limiting eigenvalues density

Integral (2.4) may be seen as a gas of eigenvalues with Vandermonde interactions and evolving in a potential V (x) = ln(1 + x 2 ). Consequently, it falls into the category of integrals studied in [START_REF] Borot | Asymptotic expansion of β matrix models in the one-cut regime[END_REF][START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF][START_REF] Borot | Large N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF]. In particular, the potential is analytic on R and has only one minimum at x = 0. Moreover, since the support of the integration is restricted to a compact set [-a, a], the convergence issues are trivial. Under such conditions, it is proved in [START_REF] Wigner | Statistical properties of real symmetric matrices with many dimensions[END_REF][START_REF] Boutet De Monvel | On the statistical mechanics apporach in the random matrix theory: integrated density of states[END_REF] that the empirical eigenvalues density δ n converges almost surely towards an absolutely continuous limiting eigenvalues density:

δ n = 1 n n k=1 δ(x -t k ) → n→∞ dµ ∞ (x) = ν(x)dx (2.6)
whose support is a finite union of intervals in [-a, a]. Determining the limiting eigenvalues density can be done in many ways, but since we plan to use the topological recursion, it seems appropriate to derive the limiting eigenvalues density using the loop equations method. We define the following correlation functions: Definition 2.1 (Correlation functions) We define the correlation functions by:

W 1,a (x) = n k=1 1 x -t k W n.c. p,a (x 1 , . . . , x p ) = n i 1 ,...,ip=1 1 x 1 -t i 1 . . . 1 x p -t ip W p,a (x 1 , . . . , x p ) = n i 1 ,...,ip=1 1 x 1 -t i 1 . . . 1 x p -t ip c
where the average of a function of the eigenvalues is defined by:

g(t 1 , . . . , t n ) = 2 n 2 (2π) n n!Z n (a) [-a,a] n dt 1 . . . dt n g(t 1 , . . . , t n )∆(t) 2 e -n n k=1 ln(1+t 2 k )
The index c stands for "connected" or "cumulant" in the sense that:

W 2,a (x 1 , x 2 ) = W n.c. 2,a (x 1 , x 2 ) -W 1,a (x 1 )W 1,a (x 2 ) W 3,a (x 1 , x 2 , x 3 ) = W n.c. 3,a (x 1 , x 2 , x 3 ) -W 1,a (x 1 )W n.c. 2,a (x 2 , x 3 ) -W 1,a (x 2 )W n.c. 2,a (x 1 , x 3 ) -W 1,a (x 3 )W n.c. 2,a (x 1 , x 2 ) + 2W 1,a (x 1 )W 1,a (x 2 )W 1,a (x 3 ) etc.
or in a more general way by the inverse relation:

W n.c. p,a (x 1 , . . . , x p ) = µ {x 1 ,...,xp} l(µ) i=1 W |µ i |,a (µ i ) Integral (2.4
) is a Hermitian matrix integral with hard edges at t = ±a. Loop equations for Hermitian matrix integrals with hard edges have been written in many places [START_REF] Chekhov | Matrix models with hard walls: geometry and solutions[END_REF][START_REF] Eynard | Loop equations for the semiclassical 2-matrix model with hard edges[END_REF][START_REF] Chekhov | Solving matrix models in the 1/N -expansion[END_REF]. They can also be easily obtained with the integral:

1 (2π) n (n!)Z n (a) [-a,a] n dt 1 . . . dt n n j=1 d dt j 1 x -t j ∆(t 1 , . . . , t n ) 2 e -n n k=1 ln(1+t 2 k ) (2.7)
Indeed, the last integral is equivalent to:

W 2 1,a (x) + W 2,a (x, x) - 2nx 1 + x 2 W 1,a (x) + n n k=1 V (x) -V (t k ) x -t k = c 1 (n) x -a - c 2 (n) x + a (2.8)
Equation (2.8) is exact and the coefficients c 1 (n) and c 2 (n) are given by:

c 1 (n) = e -n ln(1+a 2 ) (2π) n (n!)Z n (a) n j=1 [-a,a] n-1 dt 1 . . . dt j-1 dt j+1 . . . dt n ∆(t 1 , . . . , t j-1 , a, t j+1 , . . . , t n ) 2 e -n n k =j ln(1+t 2 k ) c 2 (n) = e -n ln(1+a 2 ) (2π) n (n!)Z n (a) n j=1 [-a,a] n-1 dt 1 . . . dt j-1 dt j+1 . . . dt n ∆(t 1 , . . . , t j-1 , -a, t j+1 , . . . , t n ) 2 e -n n k =j ln(1+t 2 k )
(2.9)

Note that since the integral (2.4) is invariant under the change t → -t we automatically get c 1 (n) = c 2 (n). In order to obtain the spectral curve of the problem (which is the Stieltjes transform of the limiting eigenvalues density), we take the leading order in n of equation (2.8). Results from [START_REF] Borot | Asymptotic expansion of β matrix models in the one-cut regime[END_REF][START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF][START_REF] Borot | Large N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF] show that W 1,a (x)

∼ n→∞ nW (0) 1,a (x) and W 2,a (x 1 , x 2 ) = n→∞ O(1). Defining y(x) = W (0) 1,a (x) - 1 2 V (x) = W (0) 1,a (x) - x 1 + x 2
(2.10)

we end up with:

y(x) 2 = x 2 (1 + x 2 ) 2 - 2 1 + x 2 lim n→∞ 1 n n k=1 1 1 + t 2 k -x lim n→∞ 1 n n k=1 t k 1 + t 2 k + c(a) x -a - c(a)
x + a (2.11) where the constant c(a) is given by c(a) = lim

n→∞ 1 n c 1 (n). Since the integral (2.4) is invariant under t → -t we get that n k=1 t k 1+t 2 k = 0 so that: y(x) 2 = x 2 (1 + x 2 ) 2 - 2d(a) 1 + x 2 + c(a) x -a - c(a) x + a (2.12)
where d(a) and c(a) are so far undetermined constants (i.e. independent of x). Moreover, note that by definition we must have at large x:

W (0) 1,a (x) = 1 x -lim n→∞ 1 n n k=1 t k 1 x 2 + O 1 x 3 = 1 x + O 1 x 3 ⇒ y(x) = O 1 x 3 (2.13)
since the integral is invariant under t → -t. Using the fact that y 2 (x) = O 1 x 6 in (2.12) provides two independent equations satisfied by (c(a), d(a)) that can be explicitly solved. We find:

c(a) = 1 2a(1 + a 2 ) and d(a) = 2 + a 2 2(1 + a 2 ) (2.14)
Finally, we get:

y 2 (x) = 1 + a 2 (1 + x 2 ) 2 (x 2 -a 2 ) = 1 cos 2 ( π 2 )(1 + x 2 ) 2 (x 2 -tan 2 ( π 2 )) (2.15)
In other words since cos( π 2 ) > 0 for ∈ (0, 1):

y(x) = 1 cos( π 2 )(1 + x 2 ) x 2 -tan 2 ( π 2 ) (2.16)
This is equivalent to say that the limiting eigenvalues density is given by:

dµ ∞ (x) = dx π cos( π 2 )(1 + x 2 ) tan 2 ( π 2 ) -x 2 1 [-tan π 2 ,tan π 2 ] (x)
(2.17)

The last density can be verified numerically using Monte-Carlo simulations:

Fig. 1: Empirical eigenvalues density obtained from 100 independent Monte-Carlo simulations of the integral (2.4) in the case = 1 7 and n = 20. The black curve is the theoretical curve corresponding to (2.17).

The limiting eigenvalues density is supported on the whole interval [-tan π 2 , tan π 2 ] and exhibits inverse square-root behavior near the edges of the support:

dµ ∞ (x) x→tan π 2 = O dx x -tan π 2 and dµ ∞ (x) x→-tan π 2 = O dx x + tan π 2 (2.18) Moreover it is strictly positive inside [-tan π 2 , tan π 2 ].
Remark 2.1 The spectral curve (2.16) is related to the one found in [START_REF]Matrix models, Toeplitz determinants and recurrence times for powers of random unitary matrices[END_REF] (eq. C.14) using the complex integral reformulation rather than the real integral reformulation. Indeed, in [START_REF]Matrix models, Toeplitz determinants and recurrence times for powers of random unitary matrices[END_REF], the spectral curve in the one-interval case is ỹ2 =

(x+1) 2 4x 2 (
x-e iπ )(x+e -iπ ) . In fact, both curves are equivalent up to the symplectic transformation:

x = e 2iArctan(x) def = f (x) and ỹ = 1 f (x) y = 1 + x 2 2i e -2iArctan(x) y (2.19)
This transformation follows from the combination of x = tan θ 2 and x = e iθ corresponding to a simple reparametrization of the x-function (which is a trivial case of the symplectic invariance). In particular, as explained in Appendix B, both curves provide the same set of "symplectic invariants" F (g) g≥0 that reconstruct the large n expansion of ln Z n (a). However, the Eynard-Orantin differentials attached to both curves differ which is coherent with the fact that the correlation functions are different in both settings.

General form of the large n expansions

As soon as the limiting eigenvalue density (2.17) is determined, we may apply the main results of [START_REF] Borot | Asymptotic expansion of β matrix models in the one-cut regime[END_REF][START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF][START_REF] Borot | Large N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF] to obtain the general form of the large n expansions of the correlators and of the partition function. However, we first need to prove that Hypothesis 1.1 of [START_REF] Borot | Asymptotic expansion of β matrix models in the one-cut regime[END_REF] is satisfied so that we may apply the main results of [START_REF] Borot | Asymptotic expansion of β matrix models in the one-cut regime[END_REF]. We have: Proposition 2.1 The following conditions (Hypothesis 1.1 of [START_REF] Borot | Asymptotic expansion of β matrix models in the one-cut regime[END_REF]) are met for integral (2.4) (note that in our case, the potential V does not depend on n so that some conditions of [START_REF] Borot | Asymptotic expansion of β matrix models in the one-cut regime[END_REF] are trivially verified):

• (Regularity): The potential V is continuous on the integration domain [b -, b + ].
• (Confinement of the potential): If ±∞ belong to the integration domain, then the potential is assumed to be decaying sufficiently fast:

lim inf x→±∞ V (x) 2 ln |x| > 1
• (One cut regime): The support of the limiting eigenvalues density is a single interval [α -, α + ] not reduced to a point.

• (Control of large deviations): The function x → 1 2 V (x) + R |x -ξ|dµ ∞ (ξ) defined on [b -, b + ] \ (α -, α + ) achieves its minimum only in α -or α + .
• (Off-Criticality): The limiting eigenvalues density is off-critical in the sense that it is strictly positive inside the interior of its support and behaves like

O 1 √ x-b ± if b ± is a hard edge or like O ( √ x -α ± ) if α ± is a soft edge.
• (Analyticity): V can be extended to an analytic function inside a neighborhood of

[α -, α + ].
proof:

In our case, most of the points required are easily verified:

• (Regularity): x → ln(1 + x 2 ) is obviously continuous on [b -, b + ] = [-a, a].
• (Confinement of the potential): No confinement is required since the support is a compact set of R.

• (One cut regime): This condition directly follows from equation (2.17).

• (Control of large deviations): Since α -= b -= -a and α + = b + = a (the limiting density is supported on the whole integration domain) then

[b -, b + ] \ (α -, α + ) = {α -, α + }
so the condition is trivially realized.

• (Off-Criticality): We only have two hard edges and equation (2.18) provides the correct behavior. Moreover, we directly observe from its expression that dµ ∞ (x) is strictly positive inside its support.

• (Analyticity):

x → ln(1 + x 2 ) is trivially analytic in a neighborhood of [-a, a].
Therefore, we can apply the main result of [START_REF] Borot | Asymptotic expansion of β matrix models in the one-cut regime[END_REF] for β = 2 (as well as Theorems 1.3 and 1.4 of [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF] or results of [START_REF] Borot | Large N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF] for the partition function with hard edges) and we obtain that: Theorem 2.2 (Large n expansions) The correlators and the partition functions Z n (a) admit a large n expansion (usually called "topological expansion") of the form:

W p,a (x 1 , . . . , x p ) = ∞ g=0 W {2-p-2g} p,a (x 1 , . . . , x p )n 2-p-2g Z n (a) = n n+ 1 4 n! exp ∞ k=-2 F {k} (a)n -k ln Z n (a) = - 1 4 ln n + ∞ k=-2 F {k} (a)n -k (2.20)
The previous large n asymptotic expansions have the precise meaning that ∀ K ≥ 0:

W p,a (x 1 , . . . , x p ) = K g=0 W {2-p-2g} p,a (x 1 , . . . , x p )n 2-p-2g + o n 2-p-2K ln Z n (a) = - 1 4 ln n + K k=-2 F {k} (a)n -k + o n -K (2.21)
where the o n 2-p-2K and o n -K are uniform for x 1 , . . . , x n in any compact set of [-a, a] but are not uniform in n nor K. In the rest of the article, series like

∞ k=-2
f k n -k will always be considered as asymptotic expansions in the sense presented above.

Notice that the asymptotic expansion of a given correlation function only involves powers of n with the same parity and that the series expansion of W p,a starts at 1 n p-2 . On the contrary, the large n expansion of the partition function Z n (a) may involve all powers of n and has an extra factor n n+ 1 4 and n!. Indeed, the results of [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF][START_REF] Borot | Large N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF] providing the r.h.s. of (2.20) only apply directly to (2π) n n! 2 n 2 Z n (a). While, the factors (2π) n and 2 n 2 may be absorbed in the definition of the constants F {k} (a), the term n! may not. A direct corollary of the previous theorem is that the coefficients W {2-p-2g} p,a of the correlators are obtained from the topological recursion: Corollary 2.1 (Reconstruction of the correlators via the topological recursion) For all p ≥ 1 and g ≥ 0 we have:

W {2-p-2g} p,a (x 1 , . . . , x p )dx 1 . . . dx n = ω (g) p (x 1 , . . . , x p ) (2.22)
where ω (g) p (x 1 , . . . , x p ) is the (p, g) Eynard-Orantin differential (see Appendix B or [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]) computed from the application of the topological recursion to the (genus zero) spectral curve (2.16).

Details about the topological recursion are presented for completeness in Appendix B and more can be found in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] and [START_REF] Eynard | Counting surfaces[END_REF]. In particular since the spectral curve (2.16) is of genus zero, it can be parametrized globally on C = C ∪ {∞} via:

x(z) = 1 2 tan π 2 z + 1 z y(z) = 2 sin( π 2 ) 1 + 1 4 tan 2 ( π 2 ) z + 1 z z -1 z (2.23) Moreover, the normalized bi-differential ω (0) 2 (z 1 , z 2 ) required to initialize the topological re- cursion is ω (0) 2 (z 1 , z 2 ) = dz 1 dz 2 (z 1 -z 2 ) 2
for genus zero curves. The proof of corollary 2.1 is standard. Indeed, by construction the correlations functions W p,a (x 1 , . . . , x p ) satisfy the loop equations arising in Hermitian matrix models. Moreover, Theorem 2.2 ensures that they have a topological expansion and by definition these correlators may only have singularities at the branchpoints or at the edges. These properties are also satisfied by the Eynard-Orantin differentials (See [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]) and thus both sets must match since, under these conditions (topological expansion and location of the singularities), the loop equations admit a unique solution.

The situation is more complicated for the partition function Z n (a). Indeed, as discussed in [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF] and [START_REF] Borot | Large N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF], only ∂ ln Zn(a) ∂a can be matched with ∂ ln τ (a) ∂a where ln τ is the free energies generating series associated to the spectral curve (2.16) defined by:

ln τ (a) = - ∞ g=0 F (g) Top. Rec. (a) n 2-2g
(2.24)

The coefficients F are computed from the topological recursion and are commonly called "symplectic invariants" or "free energies" (See Appendix B or [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] for the formulas). Note that there is a change of convention regarding the indexes between the indexes of the topological recursion (noted (g) and corresponding to n 2-2g ) and the indexes of the asymptotic expansion of Z n (a) (noted {k} and corresponding to n -k ). We keep the notation of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] for the topological recursion side to avoid confusion (hence the notation is

F (0) Top. Rec. , F (1)
Top. Rec. for powers n 2 , n 0 and so on). Thus, the coefficients F {k} (a) k≥-2 match with F (g) Top. Rec. (a) g≥0 only up to some constants terms (in the sense independent of a). For applications in topological string theory and integrable systems, these constants are generally disregarded because normalization issues are mostly irrelevant. But in the context of Toeplitz determinants and probability one needs to find a way to compute them. Otherwise, one may only consider ratio of Toeplitz determinants as studied in [START_REF] Duits | Relative Szegö asymptotics for Toeplitz determinants[END_REF]. In our case, we may compute these constants from the knowledge of the asymptotic expansion of lim a 0 →0 ln Z n (a 0 ) and using results of [START_REF] Iwaki | Voros Coefficients for the Hypergeometric Differential Equations and Eynard-Orantin's Topological Recursion -Part II : For the Confluent Family of Hypergeometric Equations[END_REF] as we will see below.

Normalization and limit → 0

From [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF] we have for all (a, a 0

) ∈ R * + × R * + : ln Z n (a) -ln Z n (a 0 ) = - +∞ k=-2 F {k} (a)n -k + +∞ k=-2 F {k} (a 0 )n -k = - +∞ g=0 F (g) Top. Rec. (a) n 2-2g + +∞ g=0 F (g) Top. Rec. (a 0 ) n 2-2g (2.25)
Note that when a and a 0 are positive numbers, the corresponding potentials satisfy conditions of Proposition 2.1 (in particular the limiting eigenvalues distributions remain one-cut) thus allowing the use of results of [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF] that directly gives the previous result. In order to have convergent quantities when a 0 → 0, we then isolate a trivial divergent term for g = 0:

ln Z n (a) = - +∞ g=0 F (g) Top. Rec. (a) n 2-2g + ln (Z n (a 0 )) -n 2 ln a 0 + +∞ g=0 F (g) Top. Rec. (a 0 ) + δ g=0 ln a 0 n 2-2g
(2.26)

In Appendix A.3 we prove that (ln (Z n (a 0 )) -n 2 ln a 0 ) has a limit when a 0 → 0 which is given by a Selberg integral whose asymptotic expansion can be computed exactly. More precisely, we find:

ln Z n (a 0 ) -n 2 ln(a 0 ) → a 0 →0 4 ln(G(n + 1)) -ln(G(2n + 1)) + 2n 2 ln 2 -n ln(2π) = - 1 4 ln n + 3ζ (-1) + 1 12 ln 2 + ∞ g=1 (4 -2 -2g )B 2g+2 2g(2g + 2)n 2g (2.27)
where G is the G-Barnes function, ζ is the Riemann zeta function and (B k ) k≥0 are the Bernoulli numbers. Moreover, using invariance of the free energies under symplectic transformations (we only need a trivial symplectic transformation but not the difficult and controversial (x, y) ↔ (y, -x) exchange) we prove in Appendix A.4 that the limits when a 0 → 0 of the free energies F x 2 -1 whose values have been recently computed in [START_REF] Iwaki | Voros Coefficients for the Hypergeometric Differential Equations and Eynard-Orantin's Topological Recursion -Part II : For the Confluent Family of Hypergeometric Equations[END_REF] for all g ≥ 0. This term by term convergence is extended to a large n asymptotic expansion in Appendix A.5 using results of [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF]. We find:

+∞ g=0 F (g) Top. Rec. (a 0 ) + δ g=0 ln a 0 n 2-2g a 0 →0 → n 2 ln 2 - +∞ g=2 B 2g 4 -2 -(2g-2) 2g(2g -2)n 2g-2 (2.28)
Note in particular that we have the simplification (to be understood again in the sense of large n asymptotic expansions):

ln (Z n (a 0 )) -n 2 ln a 0 + +∞ g=0 F (g0)
Top. Rec. (a 0 ) + δ g=0 ln a 0 n

Top. Rec. ( ) given in Appendix A.2 (equations (A.9) and (A.11)). We remind here that a = tan π 2 . We obtain the first orders:

F (0) Top. Rec. ( ) = ln 2 -ln sin π 2 F (1) Top. Rec. ( ) = 1 4 ln cos π 2 F (2) Top. Rec. ( ) = 1 64 - 1 32 tan 2 π 2 F (3) Top. Rec. ( ) = - 1 256 - 1 128 tan 2 π 2 - 5

Final result

Using the rotation invariance of the problem, we may easily generalize the previous result for any interval [α, β] and we obtain the following theorem: where T (α, β) = {e it , t ∈ [α, β]} admits a large n asymptotic expansion of the form (with the same meaning as the one given in Theorem 2.2):

ln det T n (1 T (α,β) ) = n 2 ln sin |β -α| 4 - 1 4 ln n - 1 4 ln cos |β -α| 4 +3 ζ (-1) + 1 12 ln 2 - ∞ g=2 F (g) (a)n 2-2g
where a = tan |β-α| 4 and the coefficients F (g) (a) g≥2 are the Eynard-Orantin free energies (also called symplectic invariants) associated to the spectral curve

y 2 (x) = 1 cos 2 |β-α| 4 (1 + x 2 ) 2 x 2 -tan 2 |β-α| 4
The previous large n asymptotic expansion has the precise meaning given in Theorem 2.2. In particular, the first orders of the expansion are given by:

ln det T n (1 T (α,β) ) = n 2 ln sin |β -α| 4 - 1 4 ln n - 1 4 ln cos |β -α| 4 +3 ζ (-1) + 1 12 ln 2 + 1 64n 2 2 tan 2 |β -α| 4 -1 + 1 256n 4 1 + 2 tan 2 |β -α| 4 + 10 tan 4 |β -α| 4 + O 1 n 6

Numerical study

We can efficiently compute the Toeplitz determinants (2.2) numerically up to n = 35. This allows to compare the theoretical formula (2.32) with the numeric simulations up to order o 1 n 4 . We obtain the following picture: We obviously see on the last figure that the numeric simulations are compatible with the theoretical results up to order o 1 n 4 . This provides additional credit for the general formulas proved in Theorem 2.3 and the reconstruction of the expansion from the topological recursion. Note that the reformulation of Z n (a) in terms of the determinant of a symmetric Toeplitz matrix is particularly useful since it allows fast computations of Z n (a) even for relatively large values of n. We performed the computations using Maple software and could compute the values of Z n (a) from n = 2 to n = 35 in no more than a few minutes on a standard laptop.

3 Toeplitz determinants with a discrete rotational symmetry

Let r ≥ 0 be a given integer and let ∈ (0, 1) be a given number. In this section we consider the Toeplitz determinants:

Z n (I r ) = 1 (2π) n n! (Ir) n dθ 1 . . . dθ n 1≤p<q≤n
e iθp -e iθq 2 with

I r = r k=-r 2πk 2r + 1 - π 2r + 1 , 2πk 2r + 1 + π 2r + 1 (3.1)
For simplicity, we denote α

(r) k = 2πk 2r+1 -π 2r+1 , β (r)
k = 2πk 2r+1 + π 2r+1 and γ (r) k = 2πk 2r+1 for -r ≤ k ≤ r.
We also define:

a (r) k = tan α (r) k 2 and b (r) k = tan β (r) k 2 for -r ≤ k ≤ r
Note that Z n (I r ) can also be interpreted as the probability to obtain all angles (θ j ) 1≤j≤n in I r . Similarly to the last section, we also introduce the sets:

J r = r k=-r a (r) k , b (r) k
= r k=-r tan πk 2r + 1 - π 2(2r + 1) , tan πk 2r + 1 + π 2(2r + 1)
T r = e iθ , θ ∈ I r Therefore, from Lemma 1.1, Z n (I r ) can be reformulated as follow:

1. Z n (I r ) = det T (r) n with T (r) n i,j
= t i-j the n × n Toeplitz matrix given by:

t 0 = |I r | 2π = t k = sin c π k 2r + 1 δ k ≡ 0 [2r+1] for k = 0 (3.2)
Note in particular that the Toeplitz matrices are mostly empty since only bands with indexes multiple of 2r + 1 are non-zero.

2. A real integral with Vandermonde interactions:

Z n (I r ) = 2 n 2 (2π) n n! (Jr) n dt 1 . . . dt n ∆(t 1 , . . . , t n ) 2 e -n n k=1 ln(1+t 2 k )
(3.3)

A complex integral:

Z n (I r ) = (-1)

n(n+1) 2 i n (Tr) n du 1 . . . du n ∆(u 1 , . . . , u n ) 2 e -n n k=1 ln u k (3.4)
The situation can be illustrated as follow:

Fig. 3: Illustration (in green) of the set T r=2 for = 1 5 .

Computation of the spectral curve

We want to compute the spectral curve associated to the integral (3.3). This integral corresponds to a Hermitian matrix integral with hard edges at a .

Following the same method as in section 2.2, we define g(t) r as the average of the function g(t) relatively to the measure induced by (3.3). With the same arguments as in section 2.2 , we get a spectral curve of the form:

y(x) 2 = x 2 (1 + x 2 ) 2 - 2 1 + x 2   lim n→∞ 1 n n j=1 1 1 + t 2 j r -x lim n→∞ 1 n n j=1 t j 1 + t 2 j r   + r j=-r A k x -a (r) k + B k x -b (r) k (3.5)
where the constants A k (resp. B k ) are given by

A k = -lim n→∞ 1 n A k,n and B k = -lim n→∞ 1 n B k,n with: A k,n = - e -n ln(1+(a (r) k ) 2 ) (2π) n (n!)Z n (I r ) n j=1 (Jr) n-1 dt 1 . . . dt j-1 dt j+1 . . . dt n a (r) k -t j ∆(t 1 , . . . , t j-1 , a (r) k , t j+1 , . . . , t n ) 2 e -n n q =j ln(1+t 2 q ) B k,n = e -n ln(1+(b (r) k ) 2 ) (2π) n (n!)Z n (I r ) n j=1 (Jr) n-1 dt 1 . . . dt j-1 dt j+1 . . . dt n b (r) k -t j ∆(t 1 , . . . , t j-1 , b (r) k , t j+1 , . . . , t n ) 2 e -n n q =j ln(1+t 2 q ) (3.6)
Note that the integral (3.3) is invariant under t → -t, thus we get that

n j=1 t j 1+t 2 j r = 0.
Therefore we end up with:

y(x) 2 = x 2 (1 + x 2 ) 2 - 2d 1 + x 2 + r j=-r A k x -a (r) k + B k x -b (r) k (3.7)
where the coefficients d and (A k , B k ) -(r-1)≤k≤r-1 are so far undetermined and independent of x. Observe now that the choice of I r implies that we have:

a (r) -k = -b (r) k for all 0 ≤ k ≤ r
Using the invariance of the integral relatively to t → -t, we obtain that:

A -k = -B k and B -k = -A k for all 0 ≤ k ≤ r (3.8)
Consequently we can reduce the spectral curve to:

y(x) 2 = x 2 (1 + x 2 ) 2 - 2d 1 + x 2 + r j=-r 2b (r) k B k x 2 -(b (r) k ) 2 (3.9)
Similarly to the one interval case, the study of large x implies that y 2 (x) = O 1

x 6 giving some algebraic relations between the 2r + 2 undetermined coefficients (d, B -r , . . . , B r ). However unlike the one interval case, when r ≥ 1, these relations are not sufficient to determine completely the coefficients. Indeed, when r ≥ 1, the number of unknown coefficients is strictly larger than the number of relations obtained from the fact that y 2 (x) = O 1 x 6 . Therefore, some additional relations are required. In order to obtain them we perform the following change of variables (denoting 1 = (1, . . . , 1) t ∈ R n , t = (t 1 , . . . , t n ) t and t = ( t1 , . . . , tn ) t ) in the integral (3.3). Let -r ≤ j 0 ≤ r :

t = tan Arctan( t) + πj 0 2r + 1 1 = t + tan πj 0 2r+1 1 -t tan πj 0 2r+1 ⇔ t = tan Arctan(t) - πj 0 2r + 1 1 = t -tan πj 0 2r+1 1 + t tan πj 0 2r+1 (3.10)
Note that the domain of integration J r is invariant under the former change of variables since any interval a

(r) k , b (r) k is mapped to a (r) k-j 0 , b (r)
k-j 0 where the indexes k -j 0 are to be understood modulo 2r + 1. Then, straightforward computations show that:

dt p = 1 + tan 2 πj 0 2r+1 1 -tp tan πj 0 2r+1 2 d tp 1 + t 2 p = 1 + tan 2 πj 0 2r+1 (1 + t2 p ) 1 -tp tan πj 0 2r+1 2 (t p -t q ) 2 = 1 + tan 2 πj 0 2r+1 2 ( tp -tq ) 2 1 -tp tan πj 0 2r+1 2 1 -tq tan πj 0 2r+1 2 (3.11)
Therefore, the general form of the integral remains invariant:

Z n (I r ) = 2 n 2 (2π) n n! (Jr) n dt 1 . . . dt n ∆(t 1 , . . . , t n ) 2 e -n n k=1 ln(1+t 2 k ) = 2 n 2 (2π) n n! (Jr) n d t1 . . . d tn ∆( t1 , . . . , tn ) 2 e -n n k=1 ln(1+ t2 k )
and thus the function W 1 (x) is given by:

W 1 (x) = 2 n 2 (2π) n (n!)Z n (I r ) (Jr) n dt n k=1 1 x -t k ∆(t) 2 e -n n k=1 ln(1+t 2 k ) = 2 n 2 (2π) n (n!)Z n (I r ) (Jr) n d t n k=1 1 -tk tan πj 0 2r+1 x -tan πj 0 2r+1 -tk 1 + x tan πj 0 2r+1 ∆( t) 2 e -n n k=1 ln(1+ t2 k ) (3.12)
We now observe that the equation:

1 -at x -a -t(1 + ax) = µ + ν x -t + ρ t x 2 -t 2 (3.13)
where (µ, ν, ρ) are independent of t admits the following solution:

(µ, ν, ρ, a) = a 1 + ax , x(1 + a 2 ) (x -a)(1 + ax) , - (1 + a 2 )(1 + x 2 ) (x -a)(1 + ax) 2 , - 2x x 2 -1 (3.14)
In particular we always have µ 1-ν = x 1+x 2 = 1 2 V (x). Inserting these results with a = tan πj 0 2r+1 , t = tk and:

   x j 0 = tan πi 0 2r+1 -π 2(2r+1)
if j 0 = 2i 0 is even

x j 0 = tan πi 0 2r+1 + π 2(2r+1) if j 0 = 2i 0 + 1 is odd (3.15)
that are solutions of a(x 2 -1) + 2x = 0 located outside of J r , we end up in (3.12) with:

W 1 (x j 0 ) = µ + νW 1 (x j 0 ) + ρ n k=1 t x 2 j 0 -t 2 r (3.16)
Using the fact that

n k=1 t x 2 -t 2 r
= 0 from the symmetry t → -t, we get to: )) are 2r distinct zeros of the function x → y(x) located outside J r . This provides 2r distinct double zeros for y 2 (x). From (3.9) and the fact that y 2 (x) = x→∞ O 1 x 6 we get that the spectral curve must be of the form:

W 1 (x j 0 ) = µ + νW 1 (x j 0 ) ⇒ W 1 (x j 0 ) = µ 1 -ν = x j 0 1 + x 2
y 2 (x) = P 4r+4 (x) (1 + x 2 ) 2 r k=-r (x 2 -(b (r) k ) 2 )
with P 4r+4 a polynomial of degree 4r (3.18)

Since we have found 4r zeros (counted with their multiplicities), we have:

y 2 (x) = λ r r-1 k=0 x 2 -tan 2 πk 2r+1 + π 2(2r+1) 2 (1 + x 2 ) 2 r k=-r x 2 -tan 2 πk 2r+1 + π 2(2r+1) (3.19)
with the convention that for r = 0 we take empty products (like

-1 k=0 x 2 -tan 2 πk 2r+1 + π 2(2r+1)
2

) equal to 1. The constant λ r can be determined using the behavior of y 2 (x) around ±i. Indeed, by definition the function z → W 1 (z) is analytic in a neighborhood of ±i. Consequently, the poles of y 2 (x) at x = ±i only comes from the shift by -1 2 V (x) = x 1+x 2 . Therefore we should get y 2 (x) ∼ x→i 1 4(x-i) 2 . We get:

1 4 = λ r r-1 k=0 1 + tan 2 πk 2r+1 + π 2(2r+1) 2 4 r k=-r 1 + tan 2 πk 2r+1 + π 2(2r+1) ⇔ λ r = r-1 k=0 cos 4 πk 2r+1 + π 2(2r+1) r k=-r cos 2 πk 2r+1 + π 2(2r+1) (3.20)
Finally we get:

y 2 (x) = r-1 k=0 cos 4 πk 2r+1 + π 2(2r+1) r k=-r cos 2 πk 2r+1 + π 2(2r+1) r-1 k=0 x 2 -tan 2 πk 2r+1 + π 2(2r+1) 2 (1 + x 2 ) 2 r k=-r x 2 -tan 2 πk 2r+1 + π 2(2r+1) = r-1 k=0 x 2 cos 2 πk 2r+1 + π 2(2r+1) -sin 2 πk 2r+1 + π 2(2r+1) 2 (1 + x 2 ) 2 r k=-r x 2 cos 2 πk 2r+1 + π 2(2r+1) -sin 2 πk 2r+1 + π 2(2r+1) = r-1 k=0 (x 2 + 1) cos 2 πk 2r+1 + π 2(2r+1) -1 2 (1 + x 2 ) 2 r k=-r (x 2 + 1) cos 2 πk 2r+1 + π 2(2r+1) -1 (3.21)
The corresponding limiting eigenvalues density is therefore given by:

dµ ∞ (x) = r-1 k=0 cos 2 πk 2r+1 + π 2(2r+1) r k=-r cos πk 2r+1 + π 2(2r+1) r-1 k=0 tan 2 πk 2r+1 + π 2(2r+1) -x 2 π(1 + x 2 ) r k=-r tan 2 πk 2r+1 + π 2(2r+1) -x 2 1 Jr (x)dx = r-1 k=0 (1 + x 2 ) cos 2 πk 2r+1 + π 2(2r+1) -1 π(1 + x 2 ) r k=-r (1 + x 2 ) cos 2 πk 2r+1 + π 2(2r+1) -1 1 Jr (x)dx (3.22)
In particular, one can verify that it is properly normalized: Jr dµ ∞ (x) = 1. Note that in the case r = 0, we recover the spectral curve for a single interval (2.16). We can verify numerically the last limiting eigenvalues density with Monte-Carlo simulations:

Filling fractions

We want to determine the proportion of eigenvalues in each of the 2r + 1 intervals of the limiting eigenvalues density. We start with the previous result:

ydx = λ r r-1 j=-r x -tan jπ 2r+1 + π 2(2r+1) (1 + x 2 ) r k=-r x -tan πk 2r+1 -π 2(2r+1) x -tan πk 2r+1 + π 2(2r+1) dx (3.23)
and observe that the filling fraction corresponding to the k th 0 interval, with -r ≤ k 0 ≤ r, is given by:

k 0 +r+1 def = 1 iπ tan( πk 0 2r+1 + π 2(2r+1) ) tan( πk 0 2r+1 -π 2(2r+1) ) dµ ∞ (x) (3.24)
Let -r ≤ j 0 ≤ r be a given integer with j 0 = 0. We perform the change of vari-

ables x = x-tan πj 0 2r+1
1+x tan

πj 0 2r+1
described above in the previous integral. Observe in particular that we have dx 1+x 2 = dx 1+x 2 . Moreover under this change of variable, the interval

x ∈ tan πk 2r+1 -π 2(2r+1) , tan πk 2r+1 + π 2(2r+1) is directly mapped to the interval x ∈ tan π(k+j 0 ) 2r+1 -π 2(2r+1) , tan π(k+j 0 ) 2r+1 + π 2(2r+1)
. We also have the identities:

x -tan πk 2r + 1 ± π 2(2r + 1) = x -tan π(k+j0) 2r+1 ± π 2(2r+1) 1 -tan πk 2r+1 ± π 2(2r+1) tan πj0 2r+1 1 + x tan πj0 2r+1 x -tan i 0 π 2r + 1 + π 2(2r + 1) = x -tan (i0+j0)π 2r+1 + π 2(2r+1) 1 -tan πi0 2r+1 + π 2(2r+1) tan πj0 2r+1 1 + x tan πj0 2r+1 if i 0 + j 0 ≡ r[2r + 1] (3.25)
When i 0 + j 0 ≡ r[2r + 1] then we get instead:

x -tan i 0 π 2r + 1 + π 2(2r + 1) = -tan πj 0 2r+1 -tan i 0 π 2r+1 + π 2(2r+1) 1 + x tan πj 0 2r+1 = -tan πj 0 2r+1 - 1 tan( πj 0 2r+1 ) 1 + x tan πj 0 2r+1 = -1 cos πj 0 2r+1 sin πj 0 2r+1 1 + x tan πj 0 2r+1 (3.26)
Observe now that the powers of 1 + x tan πj 0 2r+1 produce a factor with power 1 at the numerator and that we can express it like:

1 + x tan πj 0 2r + 1 = tan πj 0 2r + 1 x -tan (r + j 0 )π 2r + 1 + π 2(2r + 1) (3.27)
Therefore collecting the terms depending on x we get:

r k=-r x -tan π(k + j 0 ) 2r + 1 - π 2(2r + 1)
x -tan π(k

+ j 0 ) 2r + 1 + π 2(2r + 1) = r k=-r x -tan πk 2r + 1 - π 2(2r + 1) x -tan πk 2r + 1 + π 2(2r + 1) x -tan (r + j 0 )π 2r + 1 + π 2(2r + 1) r-1 i 0 =-r,i 0 +j 0 ≡r[2r+1] x -tan (i 0 + j 0 )π 2r + 1 + π 2(2r + 1) = r i 0 =-r,i+j 0 ≡r[2r+1] x -tan (i 0 + j 0 )π 2r + 1 + π 2(2r + 1) = r-1 i 0 =-r x -tan i 0 π 2r + 1 + π 2(2r + 1) (3.28)
Thus we obtain:

k 0 +r+1 = C j 0 k 0 +j 0 +r+1 (3.29)
where the constant C j 0 is given by:

C j0 = - r-1 i0=-r,i0+j0 ≡r[2r+1] 1 -tan πi0 2r+1 + π 2(2r+1) tan πj0 2r+1 cos 2 πj0 2r+1 r k=-r 1 -tan πk 2r+1 -π 2(2r+1) tan πj0 2r+1 1 -tan πk 2r+1 + π 2(2r+1) tan πj0 2r+1 = - r-1 i0=-r,i0+j0 ≡r[2r+1] cos π(i 0 +j 0 ) 2r+1 + π 2(2r+1) cos( πi 0 2r+1 ) cos( πj 0 2r+1 ) cos 2 πj0 2r+1 r k=-r cos π(k+j 0 2r+1 + π 2r+1 cos π(k+j 0 ) 2r+1 -π 2r+1 cos( πk 2r+1 -π 2r+1 ) cos( πk 2r+1 + π 2r+1 ) cos 2 ( πj 0 2r+1 ) = cos πj 0 2r + 1 -2-(2r-1)+(2r+1) = 1 (3.30)
Hence we have just proved that all filling fractions are equal:

Proposition 3.1 The filling fractions are the same in each interval:

∀ 1 ≤ k ≤ 2r + 1 : k = 1 2r + 1 Remark 3.1
This result is compatible with the one found in [START_REF]Matrix models, Toeplitz determinants and recurrence times for powers of random unitary matrices[END_REF] where the filling fractions were also proved to be equal.

General form of the large n asymptotic

After determining the limiting eigenvalues density, the next step is to determine the general form of the asymptotic of the correlation functions and of the partition function. However, since the spectral curve is of strictly positive genus (g = 2r), the general form of the asymptotic is more complicated than in the one interval case. Using the main result of [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF] or [START_REF] Borot | Large N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF] we can prove the following: Theorem 3.1 We have the following large n expansion:

Z n (I r ) = n n+ 1 4 (2r+1) n! exp ∞ k=-2 n -k F {k}                  m≥0 l 1 ,...,lm≥1 k 1 ,...,km≥-2 m i=1 (l i +k i )>0 n - m i=1 (l i +k i ) m! m i=1 F {k i },(l i ) l i ! • ∇ ⊗ m i=1 l i ν                  Θ -n 0 F {-2}, (2) (3.31)
where Θ is the Siegel theta function:

Θ γ (ν, T) = m∈Z g exp - 1 2 (m + γ) t • T • (m + γ) + ν t • (m + γ)
and F {2k},(l) are defined as the l th derivative of the coefficient F {2k} relatively to the filling

fractions = ( 1 , . . . , 2r+1 ) t ∈ {u ∈ (Q + ) 2r+1 / n i=1 u i = 1}.
The "optimal" filling fractions corresponds to the vector of filling fractions (i.e. the proportion of eigenvalues in each interval of the support of the limiting eigenvalues density) of the limiting eigenvalues density:

= 1 2r + 1
, . . . , . . . , 1 2r + 1 t Again large n asymptotic expansions presented in this theorem are to be understood as asymptotic expansions up to any arbitrary large negative power of n as in Theorem 2.2.

The proof of the last theorem consists in verifying the conditions required to apply the main theorem of [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF] and [START_REF] Borot | Large N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF]. The conditions are very similar to the one-interval case but for completeness we summarize them here:

• (Regularity): The potential V is continuous on the integration domain r k=-r [b (k) -, b (k)
+ ]. In our case x → ln(1 + x 2 ) is obviously continuous on R.

• (Confinement of the potential): Not required since the integration domain is a compact set of R.

• (Genus 2r regime): The support of the limiting eigenvalues density is given by the union of 2r + 1 single intervals [α

+ ] not reduced to a point. This is trivial from the explicit expression (3.22).

• (Control of large deviations): The function x → 1 2 V (x) + R |x -ξ|dµ ∞ (ξ) defined on r k=-r [b (k) -, b (k)
+ ] \ r k=-r (α (k) -, α (k)
+ ) achieves its minimum only at the endpoints α

+ . In our case this condition is trivial since the limiting eigenvalues density spans the whole integration domain.

• (Off-Criticality): The limiting eigenvalues density is off-critical in the sense that it is strictly positive inside the interior of its support and behaves like

O 1 √ x-b ± if b ± is a hard edge or like O ( √ x -α ± ) if α ± is a soft edge.
In our case, we have 2(2r + 1) hard edges and the explicit expression (3.22) provides the correct behavior. Moreover it is obvious from (3.22) that dµ ∞ (x) is strictly positive inside its support.

• (Analyticity): V can be extended into an analytic function inside a neighborhood of the integration domain. In our case, x → ln(1 + x 2 ) is obviously analytic in a neighborhood of any compact set of R.

As one can see, the general form of the large n expansion in the multi-cut regime is much more complicated than in the one-cut regime since eigenvalues may move from an interval to another and thus a summation on the filling fractions is necessary and adds new terms (given by the Siegel Theta function). In particular, orders O(1) and beyond exhibit a far more complex structure than in the one-cut case. Indeed, equation 3.31 and the fact that the present situation has a discrete symmetry (in particular = 1 2r+1 , . . . , 1 2r+1 t ) imply that these orders depend on the remainder of the Euclidean division of n by 2r + 1.

The term exp

∞ k=-1 n -2k F {2k}
is usually called the "perturbative" or "formal" part of the expansion. It is connected to the symplectic invariants computed from the topological recursion by the following proposition (See [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF] for details):

Proposition 3.2 The coefficients F {k} k≥-2
are related to the symplectic invariants

F (g)
g≥0 computed from the topological recursion applied to the spectral curve (3.21) by:

∀ k ≥ -1 : F {2k} = -F (2k+2) + f 2k with f 2k independent of ∀ k ≥ -1 : F {2k+1} = f 2k+1 with f 2k+1 independent of
Note that the non-perturbative part of (3.31) starts at O(1). The main difficulty of the expansion (3.31) lies in the fact that the non-perturbative part requires the knowledge of the spectral curve in the case where the filling fractions are arbitrarily fixed in order to compute the derivatives relatively to the filling fractions. In our case, if we take arbitrary filling fractions, the symmetries (3.17) are lost and thus the determination of the spectral curve requires to solve fixed filling fractions conditions:

k = A k dµ ∞ (x) for all -r ≤ k ≤ r with A k a closed contour circling the interval tan πk 2r+1 -π 2(2r+1) , tan πk 2r+1 + π 2 (2r+1)
. Unfortunately, solving analytically these conditions remains an open challenge and therefore the non-perturbative part of (3.31) remains mostly out of reach for theoretical computations. Another possibility to bypass this difficulty may be to rewrite derivatives relatively to filling fractions of F {2k},(l) as some integrals over B-cycles of the correlation functions ω (k+1) l using results of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]. In this case, there appears to have no need for using arbitrary filling fractions but the strategy requires the computations of all correlation functions. However this approach seems more convenient than dealing with arbitrary filling fractions. In particular, the case m = 0 in (3.31) simply gives:

Term m = 0 in (3.31) : Θ -n 0 F {-2},(2) = Θ -n 0 -2iπτ | = * (3.32)
where τ is the Riemann matrix of periods of the spectral curve. In the case = * = 1 2r+1 , . . . , 1 2r+1 t the rotational symmetry of the problem may allow the computation of the Riemann matrix of periods even if the expression of the spectral curve is rather complicated (See (3.21)). Since we will have no use of this order, we let it as an open problem.

Normalization and limit → 0

Finally, we also note that even the perturbative part is challenging in the multi-cut regime. Indeed, proposition 3.2 determines the perturbative part up to some constants. The standard way to compute them is to choose the parameters (in our case ) in such a way that the initial integral is explicitly connected to a known case. In the one cut regime, we connected the integral to a Selberg integral after a proper rescaling. In the multi-cut case, we may also compute the limit of Z n (I r ) when → 0 since the integrals decouple in this regime. We obtain the following theorem:

Theorem 3.2 (Limit of Z n (I r ) when → 0) For a given r ≥ 0 and n ≥ 1 we have:

Z n (I r ) →0 ∼ (2π ) (2r+1) n 2r+1 2 +(2 n 2r+1 +1)mn (2π) -n n 2r+1 -1 k=1 k! 4(2r+1) n 2r+1 ! 4mn 2 n 2r+1 -1 k=1 k! 2r+1 2 n 2r+1 ! mn 2 n 2r+1 + 1 ! mn →0 ∼ (2π ) n 2 -m 2 n 2r+1 +mn (2π) -n   n-mn 2r+1 -1 k=1 k!   4(2r+1) n-mn 2r+1 ! 4mn   2n-2mn 2r+1 -1 k=1 k!   2r+1 2n-2mn 2r+1 ! mn 2n-2mn 2r+1 + 1 ! mn (3.33)
This formula may easily be verified numerically by computing the leading order in the → 0 series expansion of the Toeplitz determinants given by (3.2) for low values of n (up to 25) and low values of r (up to 10). Moreover, the formula reduces for r = 0 to (A.13). When r > 1, we observe as expected that the result depends both on the value of n but also on the value of n mod (2r + 1). This is of course in agreement with the asymptotic expansion given by (3.31).

proof:

Let us first decompose the partition function as a sum over all possible partitions of the eigenvalues in the 2r + 1 intervals:

Z n (I r ) = 1 (2π) n n! n∈Nn,r n n -r , . . . , n r [α (r) -r ,β (r)
-r ] dθ 1 . . . dθ n-r . . .

(n k )! [α (r) -r ,β (r)
-r ] dθ 1 . . . dθ n-r . . . stands for the partition function with filling fractions fixed to n = (n -r , . . . , n r ) t . Note that the normalization prefactor 1 (2π) n n! in the definition of the partition functions exactly cancels the combinatorial coefficients appearing in the decomposition.

In the limit → 0, interactions between eigenvalues in two different intervals reduce to interactions between the center of the two intervals. Thus, the integrals decouple and Z n -r ,...,nr (I r ) reduces to a product of one-cut partition functions whose values are given by (A.13) (with a 0 = tan π 2(2r+1) ): We now observe that the exponent in is minimal when the sum r k=-r

Z n -r ,...,nr (I r ) →0 ∼ r k=-r Z n k ( ) = (2π )
n 2 k is minimal.
Since the sum is constrained to r k=-r n k = n fixed, it is well-known that the previous sum is minimal when all (n k ) k∈ -r,r are equal. This is also coherent with the fact that the perturbative analysis leads to identical filling fractions (See Proposition 3.1). Indeed, let us

define f (n) = r k=-r n 2 k , then
f (n -r , . . . , n i 0 + 1, . . . , n j 0 -1, . . . , n r ) -f (n -r , . . . , n i 0 , . . . , n j 0 , . . . , n r ) = 2(n i 0 -n j 0 + 1)

   < 0 if n i 0 ≤ n j 0 -2 = 0 if n i 0 = n j 0 -1 > 0 if n i 0 ≥ n j 0
Therefore, to minimize f on N n,r , one needs to balance the vector n with differences between the components at most 1. This proves that if m n = n mod (2r + 1) = 0 the sum is minimal when we spread the additional m n eigenvalues in m n distinct intervals, i.e. we put n 2r+1 eigenvalues in 2r + 1 -m n intervals and n 2r+1 + 1 eigenvalues in the remaining m n intervals and that all other configurations provide a strictly greater value of the exponent. Since the sum contains a finite number of terms (n and r are fixed so N n,r is a finite set) and because each of them is strictly positive (thus the sum of terms with the same exponent in cannot give of a strictly lower order in ), we obtain from (3.34) the result presented in Theorem 3.2.

Final result and conjecture for F {0}

Note that the knowledge of the equivalent of Z n (I r ) when → 0 given by Theorem 3.2 is not sufficient to determine the constants f {k} k≥-2 of Proposition 3.2. Indeed, in order to obtain them, one needs to exchange the limit → 0 and the limit n → ∞ and possible entropy contributions may arise. It is also unclear if arguments similar to the one-cut case may be adapted to this context. The only known cases are F {-2} and F {-1} (given in section 1.4 of [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF] with the important reminder that we included a normalization factor 2 n 2 (2π) n n! in the definition of the partition function that is absent in [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF]) . We get:

F {-2} = ln 2 -F (0) ⇔ f -2 = ln 2 F {-1} = 0 ⇔ f -1 = 0
Remark 3.2 Note that the constants f -2 = ln 2 and f -1 = 0 differ from those given in [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF] because our choice of normalization of the partition function is different from the one used in [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF]. Indeed, our partition function is normalized with an additional factor 2 n 2 (2π) n n! compared to the one used in [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF]. In particular, the numerator 2 n 2 provides an additional factor n 2 ln 2 in the asymptotic expansion of ln Z n (I r ) thus giving f -2 = ln 2. Similarly, the denominator

(2π) n n! n→∞ ∼ (2π) n n n e -n √
2πn provides an additional factor n ln n already included in the general form of the asymptotic expansion of ln Z n (I r ) given in Proposition 3.1 as well as an additional factor -n(ln(2π) -1) in the asymptotic expansion of ln Z n (I r ) that precisely cancels the one obtained in [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF], hence leading to f -1 = 0. Thus, we end up with: Proposition 3.3 (First orders of the asymptotic expansion of Z n (I r )) For a given r ≥ 0, we have:

ln(Z n (I r )) = n 2 2r + 1 ln sin π 2 - 2r + 1 4 ln n + O(1)
As explained above, the next orders of the large n expansion, and in particular the O(1) term, exhibit a far more complex structure than in the one-cut case and depend on the remainder of the Euclidean division of n by 2r + 1. In fact, as discussed in [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF], for a given remainder m such that 0 ≤ m ≤ 2r, the sequence Z (2r+1)p+m p∈N admit a large p expansion taking a similar form as the one-cut case with n replaced by p (the precise form being given in Theorem 2.2) and with coefficients depending on the remainder m. However, because we could compute numerically the Toeplitz determinants up to sufficiently large n, we can propose some conjectures for the coefficient F {0} depending on the value of the remainder m. Conjecture 3.1 (Conjecture for F {0} ) We conjecture the following:

• In the case m = 0, F {0} has a dependence in given by -2r+1 4 ln cos π 2 + C 0 where C 0 is a constant. This generalizes the one-cut case that would correspond to r = 0.

• The cases m = j and m = 2r -1 -j for some 1 ≤ j ≤ r provides the same F {0} coefficient. However, the dependence in is more involved than for m = 0 and we conjecture that it is given by:

F {0} (j) = A j ln cos π 2 + B j ln tan π 2 + C j
with some non-zero constants (A j , B j , C j ) depending on the remainder.

• The constants (C j ) 1≤j≤r and C 0 are the same and correspond to the normalization issue of the partition function.

The case m = 0 exhibits an additional symmetry since the number of eigenvalues is precisely a multiple of the number of intervals so that they can spread evenly in all intervals. In particular when evaluating derivatives of the free energies at = 1 2r+1 1 cancellations are more likely to happen, thus explaining why the coefficient in front of ln tan π 2 vanishes only in this case. We illustrate our conjecture with numeric simulations: More precisely, the black curves correspond to the best matches with curves of the form:

f ( ) = α ln cos π 2 + β ln tan π 2 + γ (3.35)
with rational coefficients (α, β, γ) of the form i 256 with i ∈ Z. Note that we chose to express the coefficients as rational numbers with a specific denominator, but we have no evidence that the coefficients are indeed rational or that the denominator is a power of 2. However it seems that this particular choice provides very accurate results that are not numerically improved by adding multiplicative powers of 3 or 5 in the denominators or increasing the power of 2. Numerically, we obtain the best matches for the values:

• Case r = 1: m = 0 : (α, β, γ) = - 3 4 , 0, - 63 128 m ∈ {1, 2} : (α, β, γ) = - 11 128 , 85 128 , - 63 128
• Case r = 2:

m = 0 : (α, β, γ) = - 5 4 , 0, - 3 16 m ∈ {1, 4} : (α, β, γ) = - 15 32 , 203 256 , - 3 16 m
∈ {2, 3} : (α, β, γ) = - 1 16 , 6 5 , - 3 16
• Case r = 3:

m = 0 : (α, β, γ) = - 7

Even number of intervals

The method developed in the last section can be adapted in the case of an even number of intervals. However we need to be careful since in order to apply θ → tan θ 2 we need to avoid the angles θ = ±π. Therefore we use the invariance θ → θ + Cste 1 of the integral (3.36) to shift the intervals so that they do not contain ±π. We define for s ≥ 1: = t p-q the n × n Toeplitz matrix given by: The situation can be illustrated as follow: Fig. 7: Illustration (in green) of the set T s=3 for = 1 5 . The method developed in the case of an odd number of intervals can be easily adapted to this case. In particular we can prove the following: The O(1) order depends on the remainder m n of n modulo 2s. Numerically, its dependence in follows the conjecture 3.1 (with 2r + 1 replaced by 2s and the special case still given by m n = 0).

α (s) k = π(k -1 2 ) s - π 2s , β (s)
k = π(k -1 2 ) s + π 2s , γ (s)

5.

The limit → 0 of Z n (I s ) with n = m n mod 2s is given by:

Z n (I s ) →0 ∼ (

Conclusion and outlooks

In this article, we proposed a rigorous mathematical derivation of the asymptotic expansions of Toeplitz determinants with symbols given by f = 1 T (α,β) where T (α, β) = {e it , t ∈ [α, β]}. This generalizes Widom's result and completes the approach developed in [START_REF]Matrix models, Toeplitz determinants and recurrence times for powers of random unitary matrices[END_REF]. We also provided numerical simulations up to o 1 n 4 to illustrate these results. For symbols f = 1 T d with T d = d k=1 [α k , β k ] and d ≥ 2, the situation is more complex, but we were able to provide a rigorous derivation of the large n asymptotic of the corresponding Toeplitz determinants when the arc-intervals exhibit a discrete rotational symmetry on the unit circle. We also provided the first terms of the large n expansion up to O(1) and we proposed a conjecture for the O(1) term supported by numeric simulations. Moreover, the results presented in this article raise the following challenges:

• Prove conjecture 3.1 regarding the O(1) term in the symmetric multi-cut case. In particular, it would be interesting to find a way to obtain the normalization constants. This requires to connect the Toeplitz integral to a known integral for at least one value of the parameter (very likely → 0).

• In the multi-cut cases without discrete symmetry, it would be interesting to prove that, for generic choices of the edges, the corresponding spectral curves are regular and that the hypothesis required to prove the large n asymptotic of Theorem 3.1 are verified. Contrary to the case with a discrete rotational symmetry, it seems unlikely that we obtain an explicit formula for the spectral curves. However it may be possible to obtain sufficient information on the location of the zeros of y(x) in order to prove that the spectral curves are regular.

• The method developed in this article could also be tried for more general Toeplitz determinants. In particular, the transition from a symbol supported on a compact set of {e it , t ∈ (-π, π)} to a strictly positive symbol on the unit circle deserves some analysis. Indeed, in the case of a strictly positive symbol on the unit circle, Szegö's theorem implies that 1 n ln det T n admits a non-trivial limit (given by 1 2π 2π 0 ln(f (e iθ ))dθ) so that the large n expansion proposed in Theorem 3.1 requires some adaptations. However, it may happen that the sub-leading corrections may still be given by Theorem 3.1 and in particular that sub-leading corrections may be reconstructed from the topological recursion. This conjecture is supported by the fact that Toeplitz determinants can also be reformulated as Fredholm determinants (See [START_REF] Borodin | A Fredholm determinant formula for Toeplitz determinants[END_REF]) that are known to be deeply related with the topological recursion.

1

  Introduction: General setting and several reformulations of the problem 1.1 General setting In this article we are interested in the computation of Toeplitz integrals of the form: Z n (I) = 1 (2π) n n! I n dθ 1 . . . dθ n 1≤p<q≤n e iθp -e iθq 2 (1.1)

Lemma 1 . 1 ( 2 ]

 112 Various reformulations of the problem) Defining I = d j=1 [α j , β j ] ⊂ (-π, π), T = d j=1 {e it , t ∈ [α j , β j ]} and J = , the following quantities are equal to each other: = 1 (2π) n n! I n dθ 1 . . . dθ n 1≤p<q≤n e iθp -e iθq 2 (1.4)

  . (a) g≥0

  (g) Top. Rec. (a 0 ) + δ g=0 n 2 ln a 0 g≥0 are given by the free energies of Legendre's spectral curve y = 1 √

Theorem 2 . 3 (

 23 Asymptotic expansion of Toeplitz determinants in the one interval case) For (α, β) such that 0 < |β -α| < 2π, the Toeplitz determinant with symbol f = 1 T (α,β)

Fig. 2 : 2 , θ → -1 4 ln cos θ 2 +tan 2 θ 2 + 10 tan 4 θ 2 .

 22222 Fig. 2: Computations (from (2.2)) of the Toeplitz determinants θ → ln Z n ( θ π ) with 0 < θ < π for 2 ≤ n ≤ 35 with subtraction of the first coefficients of the large n expansion (2.32) (Colored dots: starting from orange to yellow, green and purple as n increases). The black curves are the theoretical predictions given by (from top to bottom and from left to right): θ → ln sin θ 2 , θ → -1 4 ln cos θ 2 + 3 ζ (-1) + 1 12 ln 2, θ → 1 64 2 tan 2 θ 2 -1 and θ → 1 256 1 + 2 tan 2 θ 2 + 10 tan 4 θ 2 .

Fig. 4 :Fig. 5 :

 45 Fig. 4: Histogram of 50 independent simulations of the eigenvalues density induced by (3.3) in the case r = 1, = 1 10 and n = 60. The black curve is the theoretical limiting eigenvalues density computed in equation (3.22)

  r )! . . . (n r )! Z n-r,...,nr (I r ) = n∈Nn,r Z n-r,...,nr (I r ) where N n,r = {(i -r , . . . , i r ) ∈ N 2r+1 / r k=-r i k = n}, n n -r ,...,nr stands for the multinomial coefficient and Z n-r,...,nr (I r )

  dθ n-nr+1 . . . dθ n 1≤p<q≤n e iθp -e iθq 2

Fig. 6 : 2 + 2 and ln tan π 2 as

 6222 Fig. 6: Representation of ln(Z n (I r )) -1 2r+1 ln sin π 2 + 2r+1 4 ln n for values of n ranging from 2 to 70 and classified by the remainder m (from left to right: m = 0, m ∈ {j, 2r -1 -j}, with j from 1 to r) of the Euclidean division of n by 2r + 1. The black curves correspond to the best numerical matches of an affine combination of ln cos π 2

2 T

 2 s = {e it , t ∈ I s } and the integral:Z n (I s ) = 1 (2π) n n! (Is) n dθ 1 . . . dθ n 1≤p<q≤n e iθp -e iθq 2 (3.36)which is also equal to det T

  0 [2s] for k = 0 Note again that the Toeplitz matrix is mostly empty since only bands with indexes multiple of 2s are non-zero. Moreover it is also equal to a Hermitian integral:Z n (I s ) = 2 n 2 (2π) n n! (Js) n dt 1 . . .dt n ∆(t 1 , . . . , t n ) 2 e n du 1 . . . du n ∆(u 1 , . . . , u n ) 2 e

Theorem 3 . 3 (x 2 -tan 2 πk 2s 2 ( 1 + x 2 ) 2 s 2 .x 2 -1 4 .

 332212224 Results for an even number of symmetric arc-intervals) In the case of an even number of symmetric arc-intervals (3.36) we have:1. The spectral curve attached to integral (3.36) is given by: The corresponding limiting eigenvalues distribution is given by:dµ ∞ (x) = tan 2 π(k-Js (x)dx3. The filling fractions are the same in all the 2s intervals: ∀ k ∈ -(s -1), s : k+s def The function ln Z n (I s ) admits a large n expansion given by Theorem 3.1 with partial reconstruction by the topological recursion given by proposition 3.2. In particular we have for s ≥ 1: ln(Z n (I s )

2π ) 2s n 2s 2 +

 2

)(z-1) 2 (z+1) 2 , then the Bergman tau-function in the present case is ln τ B =

8 ln cos π 2 .

Acknowledgments

I would like to thank A. Guionnet that indirectly gave me the idea to review in a rigorous way the tools and theorems used in this article. I also would like to thank G. Borot for very fruitful discussions about the normalizing constants and also an unknown referee for his suggestions. Eventually, I would like to thank Université de Lyon, Université Jean Monnet and Institut Camille Jordan for financial and material support.

A Topological recursion and normalization for the onecut case

In this appendix, we provide details about the computation of the free energies F

A.1 Computation of the topological recursion

In this section we present the computation of the topological recursion to the spectral curve:

The spectral curve is of genus 0 and the branchpoints are located at z = ±1 and we can define a global involution z = 1 z for which x(z) = x(z) and y(z) = -y(z). Note that y has simple poles at the branchpoints so that the one-form ydx given by:

is regular at the branchpoints.

We now need to compute the first free energies F (g) Top. Rec. ( ) g≥0 attached to this genus 0 spectral curve using the formulas presented in Appendix B. Specific formulas presented in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] are required for F (0) Top. Rec. ( ) and F [START_REF] Szegö | On certain Hermitian forms associated with the Fourier series of a positive function[END_REF] Top. Rec. ( ) and the corresponding computations are presented in Appendix A.2 (equations (A.9) and (A.11)). We find (remind that a = tan π 2 ):

p,a (x 1 , . . . , x p ) = 0 for all p ≥ 3 because y has a simple pole at the branchpoints. Moreover, from the definition and the form of the spectral curve, it is also easy to see that the coefficients W

A.2 Computation of F

Top.Rec and F

(1)

Top.Rec

We want to compute the first two free energies of the curve:

We first observe that the one-form ydx is given by:

with:

and t 1 = 1 2

and t 2 = -1 2

and

and

Thus the one form ydx has 4 simple poles that are not branchpoints but only poles of y(z).

The local coordinate around each point is given by:

Hence the local potential:

is trivially vanishing in all four cases since the poles are simple. Eventually we end up with the computation of:

we have:

Hence in the end, observing that

t j = 0 we get that the dependence in o vanishes (as claimed in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]) and we find:

The computation of F

Top.Rec. in the case of two hard edges is relatively straightforward. First, we observe that we have:

Then the computation of F

Top.Rec. can be performed with the formalism of [START_REF] Chekhov | Matrix models with hard walls: geometry and solutions[END_REF] that includes the possibility of hard edges. Straightforward, but lengthy, computations using equation [149] of [START_REF] Chekhov | Matrix models with hard walls: geometry and solutions[END_REF] gives:

An alternative approach with detailed computations using a refinement of the formalism of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] adapted to hard edges is developed in [START_REF]Matrix models, Toeplitz determinants and recurrence times for powers of random unitary matrices[END_REF]. However, we note that Bergman tau-function used in [START_REF]Matrix models, Toeplitz determinants and recurrence times for powers of random unitary matrices[END_REF] is slightly incorrect. Indeed since W

(1)

Normalization with a Selberg integral

Remark A.2 In this section, we will use the term "asymptotic expansion" and infinite series in the same sense as Theorem 2.2. In particular, infinite series like f

We have the following Selberg integral:

We get:

with:

S a 0 is a Hermitian matrix integral on I = [-1, 1] with potential V a 0 (x) = ln(1 + a 2 0 x 2 ). It is continuous in a 0 and in particular for a 0 = 0 we find:

Therefore, we have:

where the function G(z) is the G-Barnes function whose asymptotic expansion is:

Using this formula for N = n and N = 2n as well as Stirling's asymptotic formula:

we can compute the asymptotic expansion of the r.h.s. of (A. [START_REF] Zhang | Volumes of orthogonal groups and unitary groups[END_REF]). We obtain:

so that:

A. [START_REF]Matrix models, Toeplitz determinants and recurrence times for powers of random unitary matrices[END_REF] Term by term limit of the free energies when → 0

In this section, we compute the limit when goes to 0 of the free energies F

Top.Rec. ( )

of the spectral curve:

Using the symplectic transformation:

we obtain the equivalent spectral curve

In particular, since the symplectic transformation (A.22) is only a rescaling of x and y by a constant, the free energies computed from either (A.21) or (A.23) are the same for g ≥ 1.

The exceptional case g = 0 provides the same result up to a trivial factor ln tan π 2 . Finally, the new parametrization (A.23) admits a regular limit when → 0 given by:

which is exactly Legendre's spectral curve y 0 = 1 √

x 2 0 -1 whose free energies have recently been computed in [START_REF] Iwaki | Voros Coefficients for the Hypergeometric Differential Equations and Eynard-Orantin's Topological Recursion -Part II : For the Confluent Family of Hypergeometric Equations[END_REF] (Cf. section 2.3.4 and main theorem of [START_REF] Iwaki | Voros Coefficients for the Hypergeometric Differential Equations and Eynard-Orantin's Topological Recursion -Part II : For the Confluent Family of Hypergeometric Equations[END_REF]). Thus, we get:

Top.Rec. ( = 0) = 0

where (B k ) k≥0 are the Bernoulli numbers defined by:

A.5 Convergence of the asymptotic expansion

Top.Rec. ( )n 2-2g when → 0

In the previous section, we have proved that:

However, the previous results are only term by term convergence results and may only be used to study formal series in n -2g but not directly in asymptotic expansions that we deal with in the present article. In order to prove that the limit remains valid in the sense of large n asymptotic expansion as claimed in (2.28), we need to use results of [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF] rather than formal limit results found in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]. Let us consider the integrals defined by (we remind that a = tan π 2 ):

which is trivially connected to Z n (a) by equation (A.13). Note in particular that the integration domain is [-1, 1] n and thus is independent of a. On the contrary, the potential V a (u) = ln(1 + a 2 u 2 ) depends continuously in a. In particular, V a converges uniformly on [-1, 1] to V 0 = 0 when a goes to 0. Moreover, for all a ≥ 0, the potentials V a satisfy conditions of Proposition 2.1 and in particular the limiting eigenvalues density is always supported on a single interval. Therefore, we may apply results of [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF] stating that the large n asymptotic expansion of ln S a converges to the large n asymptotic expansion of ln S 0 :

Eventually, the connection between Z n (a) to S a given by (A.13) allows to extend the result to the large n asymptotic expansions of ln Z n (a) as claimed in (2.28).

Remark A. [START_REF] Widom | Asymptotic behavior of block Toeplitz matrices and determinants[END_REF] The case a = 0 only appears special at first sight because the definition of Z n (a) (See (2.4)) is taken as integrals over [-a, a] for which the value a = 0 looks ill-defined. But after a trivial transformation to integrals over [-1, 1] given by (A.13) (i.e. transforming Z n (a) to S a ), it becomes obvious that a = 0 is no different than any a > 0 regarding the application of results of [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF]. Therefore, the asymptotic expansion of ln S a converges to the asymptotic expansion of ln S 0 when a tends to 0 thus giving (A.29) and finally (2.28) using (A.13).

B The Eynard-Orantin topological recursion

In this section we briefly review the formalism of the topological recursion as presented in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]. More general versions of the topological recursion can be found in the literature but we restrict ourselves to the original simpler version of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] that is sufficient for the purposes of this article. Let us start by the definition of a spectral curve:

Definition B.1 (Spectral curve, branchpoints, normalized bi-differential) A spectral curve is the data of two meromorphic functions (x(z), y(z)) on a Riemann surface Σ of genus g. This is equivalent to the data of a polynomial P such that P (x, y) = 0 and therefore to an algebraic equation between x and y. When the genus g of Σ is strictly positive, we complete the data of the spectral curve with the choice of a basis of homology cycles (A j , B j ) 1≤j≤g such that:

Then it follows from standard results of algebraic geometry that there exists a unique symmetric bi-differential B(z 1 , z 2 ) (sometimes called "Bergman kernel") such that:

• B is holomorphic on Σ × Σ except at coinciding points where it behaves like:

• B is normalized on the basis of cycles (A j , B j ) 1≤j≤g in the following way:

The branchpoints (a j ) 1≤j≤R (with R ≥ 1) of the spectral curve are the points where dx vanishes. The spectral curve is said "regular" if the branchpoints are simple zeros of dx. When the spectral curve is regular, we can define locally around each branchpoint a holomorphic involution z → z such that x(z) = x(z).

In this paper we will only restrict ourselves to the case of regular spectral curves since the situation is much more complicated when the curve is not regular. We remark that when the spectral curve is regular and of genus 0, then there exists a global parametrization (x(z), y(z)) with z ∈ C ∪ {∞} of the spectral curve. Moreover, the holomorphic involution z → z is defined globally on the spectral curve and the normalized bi-differential B is explicitly given by B(z 1 , z 2 ) = dz 1 dz 2 (z 1 -z 2 ) 2 . We now have all the ingredients to define the correlators and free energies associated to a spectral curve. Definition B.2 (Definition 4.2 of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]) For g ≥ 0 and n ≥ 1, the Eynard-Orantin differentials (known also as "correlation functions" or "correlators") ω (g) n (z 1 , . . . , z n ) of type (g, n) associated to the spectral curve (x(z), y(z)) are defined by the following recursive relations:

Here

is called the recursion kernel, and the in the last line of (B.1) means that the cases (g 1 , I) = (0, ∅) and (g 2 , J) = (0, ∅) must be excluded from the sum.

The Eynard-Orantin differentials ω (g) n 's are meromorphic multi-differentials on Σ n and are known to be holomorphic except at the branchpoints if (g, n) = (0, 1), (0, 2). In [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF], the authors also introduced free energies (also called "symplectic invariants") F (g) g≥0 defined by: Definition B.3 (Definition 4.3 of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]) The g th symplectic invariant F (g) associated to the spectral curve (x(z), y(z)) is defined by: F (0) and F (1) are defined with specific formulas that can be found in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF].

Note that this definition extends to the case n = 0 (with the identification ω (g) 0 = F (g) ) the following property:

Res z→a i Φ(z)ω (g) n+1 (z, z 1 , . . . , z n ) for g ≥ 0 and n ≥ 0 Note also that the Eynard-Orantin differentials or the symplectic invariants do not depend on the choice of parametrization (x(z), y(z)). Eventually as suggested by their name, the symplectic invariants F (g) g≥0 are invariant under transformations of the spectral curve (x, y) → (x, ỹ) such that dx ∧ dy = dx ∧ dỹ, i.e. transformations that preserve the symplectic form dx ∧ dy. Note that the only non-trivial case is the exchange (x, y) → (y, -x) for which the invariance of the symplectic invariants is highly non-trivial and controversial. In this article we only used trivial symplectic transformations. We remind that this symplectic invariance property does not hold in general for the Eynard-Orantin differentials ω (g) n but only for the free energies.


BAKULA RINPOCHE"S FAMILY

In the middle of the 18 th century there was only a single heiress in the royal family of Zangla. Her name was Ladzes Wangmo (zla-mdzes dbang-mo) (Francke [1926(Francke [ ] 1982: 165): 165). She had to bring a makpa (mag-pa) i.e. to marry uxorilocally. According to N.T. Shakspo (2010: 147) her husband was Gyapo Tupchung (rgyal-po thub-chung ?) and according to Riaboff (1989: 205) he was Dubtse Wang Dpalde (grub-tse dbang dpal lde) from the royal family of Heniskut (or Henasku). She delivered at least three children: a daughter, Tsering Wangmo (tshe-ring dbang-mo), about whom I will speak later; a son, Lobzang Namgyal (blo-bzang rnam-rgyal), heir to the crown of Zangla; and the 18 th Bakula Rinpoche, Lobzang Yeshe Stanpa Gyaltsen (blo-bzang ye-shes bstan-pa rgyal-tshan) (Riaboff 1989: 205).

Before continuing on to the 19 th Bakula Rinpoche, I must introduce the Matho royal family.

The last king to actually ruled Ladakh was Tsespal Tondup Namgyal (tshe-dpal don-grub rnam-rgyal), the builder of the palace of Stok. Born in 1790, he ruled from 1820 till 1834. After the Dogra conquest he was dismissed and replaced by his Prime Minister Ngorub Standzin (ngos-grub bstan-"dzin), the Banka Kalhon or kalhon (bka'-blon, a kind of duke or marquis) of the Chimre-Sakti area. Then the Dogra general Zorawar Singh, suspecting that the Banka Kalhon was betraying him, dismissed him and put back Tsespal Tondup on the throne from 1839 to 1840 (Cunningham [1874] 1977: 258 and343-344). Both the king and the Banka Kalhon died from smallpox in 1840 (Datta 1972: 127). The crown prince was Tsewang Rabstan (tshe-dbang rab-brtan). He had three wives (cf Jamyang Gyaltshan).

-The first was Skalzang Dolma (bskal-bzang sgrol-ma), daughter of Tsewang Tondup (tshe-dbang don-grub), the powerful kalhon of Bazgo, and a minister of the king Tsespal Tondup Namgyal. According to Cunningham (p. 334-336), he was almost as powerful as the Banka Kalhon. Skalzang Dolma had a son, Jigmet Singge Namgyal (full name: "jigs-med chos-kyi sengge kun-ga" rnam-rgyal). -The second was Sonam Palkit (bsod-nams dpal-skyid). She was the daughter of the Banka Kalhon. She also delivered a son, Stansrung Urgyal (bstan-srung g.yul-rgyal). -The third was Zara Khatun, a Muslim and the daughter of the "prince" (cho, spelling jo) of Pashkyum, near Kargil. She had no issue.

The rule of primogeniture was not implemented. The son of the first wife received the palace of Stok, and the son of the second wife received land in Matho. Later this new younger branch built a palace in Matho.

Why was the rule of primogeniture not implemented? Why was the royal property shared between the two sons? One of my informants, Dawa Tondup (zla-ba don-grub), thinks that the Dogras wanted to treat the Banka Kalhon gently because he had become as powerful as the king. Of course the Dogras were victorious but they needed help and collaboration of notables to pacify the country.

The name of this first king of Matho was Stansrung Urgyal (bstan-srung g.yul-rgyal). He married Tsering Wangmo, the daughter of the queen of Zangla mentioned earlier. Hence the royal families of Matho and Zangla became allied for the first time.

The couple had two sons. One became a monk in Hemis monastery and the other one, Tashi Lawang (bkra-shis lha-dbang), succeeded his father.

Let"s come back to the Zangla royal family. The king, brother of the 18 th Bakula rinpoche, had three children: Tsetan Namgyal (tshe-brtan rnam-rgyal) who succeeded him as king of Zangla (Riaboff 1989: 205), Ishe Wangmo (ye-shes dbang-mo) about whom I will speak very soon, and Tsultim Nyima (tshul-khrims nyi-ma), the reincarnation of the founder of the monastery of Rizong (Kaplanian, 1981: 193-194). He is the second rinpoche in the royal family of Zangla. Tashi Lawang had five sons and one daughter. Of chief interest here are the elder son, Nangwa Taez (snang-ba mtha"-yas), and the fifth son, Chimed Dorje ("chi-med rdo-rje). Both married Ishe Wangmo. Here is the second alliance between the royal families of Matho and Zangla. From this polyandric trio the 19 th Bakula rinpoche was born. He was the fourth and last child of the trio.

Between the elder son and the 18 th Bakula there were two sisters. The second sister Norzin Wangmo (nor-"dzin dbang-mo) married the lhon-po (blon-po) of Shey, one of the most powerful nobles of the country. She gave birth to Thubstan Tsewang (thub-bstan tshe-dbang), also a rinpoche (Sera Rinpoche), who married a daughter of the royal family of Stok and is now (2017) member of the parliament in New Delhi. The elder son, Tashi Phuntsok Namgyal (bkra-shis phun-tshogs rnam-rgyal), succeeded his father as a king of Matho. He was married twice. His first wife died while delivering a child who was recognised as the reincarnation of Sras Rinpoche, the son of Tsultim Nyima, founder of Rizong monastery (Kaplanian, 1981: 194). Tashi Phuntsok Namgyal remarried and his second wife gave birth to a daughter, Namgyal Wangmo (rnam-rgyal dbang-mo). Till then things are simple: Namgyal Wangmo had to bring a makpa according to the rule given in the beginning of this article. But things became complicated because Chimed Dorje also remarried and also had a daughter, Tsering Dolkar (tshe-ring sgrol-dkar). Who was the legal heir, Namgyal Wangmo or Tsering Dolkar? A long trial took place between Namgyal Wangmo and Tsering Dolkar. Unfortunately I could not get any details of this trial and the informants give contradictory versions. Some say Namgyal Wangmo won the case. Some say she lost but Bakula Rinpoche reasoned with the Stok branch of the royal family (who was backing Tsering Dolkar) so that Tsering Dolkar renounced her claim and Matho was given to Namgyal Wangmo. They even say that that Bakula Rinpoche locked Tsering Dolkar as a nun in the Spituk Monastery for some years.

Reading this article one could think that rinpoches were always born in royal and noble families but it is not always the case.

  The 20 th Bakula Rinpoche was born in a very simple family in Nubra valley. 1 Ladakh, Physical, Statistical and Historical with notices of the surrounding areas. New-Delhi. Datta, C.L., 1975 Ladakh and Western Himalayan Politics 1819-1848. New-Delhi. Francke, A.H., [1926] 1982. Antiquities of Indian Tibet, vol. II. New-Delhi. Jamyang Gyal-tshan ("jam-byangs rgyal-tshan), s d, la-dwags kyi rgyal-rabs, Leh. Kaplanian, P., 1981, Les Ladakhi du Cachemire. Paris. Riaboff, I, 1989, Le Roi et le Moine, unpublished ph. D thesis. Paris. Shakspo, N.T., 2010, A cultural History of Ladakh. Sabu-Leh.

	Patrick Kaplanian
	BIBLIOGRAPHY
	Cunningham, A [1854] 1977

The author thanks N. T. Shakspo for inviting him to attend the international conference on Padma Bhushan Kushok Bakula Rinpoche and Kyle Gardner for correcting his bad English.