√

) where Q stands for the number of independent randomly distributed eigenvalues of the system and f stands is the parameter controlling the size of the window defining that the system is close enough to its initial state. However the authors also indicate that their (independent) model is too simple and cannot be physically relevant since the eigenvalues cannot be completely random and uncorrelated in a general setting. In fact they suggest that a much proper setting would correspond to consider the

Matrix models, Toeplitz determinants and recurrence times for powers of random unitary matrices O. Marchal † † Université de Lyon, CNRS UMR 5208, Université Jean Monnet, Institut Camille Jordan, France 1 Abstract: The purpose of this article is to study the eigenvalues u t 1 = e itθ 1 , . . . , u t N = e itθ N of U t where U is a large N × N random unitary matrix and t > 0. In particular we are interested in the typical times t for which all the eigenvalues are simultaneously close to 1 in different ways thus corresponding to recurrence times in the issue of quantum measurements. Our strategy consists in rewriting the problem as a random matrix integral and use loop equations techniques to compute the first orders of the large N asymptotic. We also connect the problem to the computation of a large Toeplitz determinant whose symbol is the characteristic function of several arc segments of the unit circle. In particular in the case of a single arc segment we recover Widom's formula. Eventually we explain why the first return time is expected to converge towards an exponential distribution when N is large. Numerical simulations are provided along the paper to illustrate the results.

Summary of the article and main results

The purpose of the article is to study the probability P N (t) that the t th power of a large random unitary matrix U N of size N × N returns close to the identity matrix when the size of the matrix tends to infinity. In the following section, we present three usual definitions characterizing the notion of being "close" to the identity matrix. The first and second definitions consist in saying that respectively the real part (real-part return time) or modulus (weak return time) of the average of the eigenvalues is close to 1. The third alternative definition states that all eigenvalues belong to a small arc-interval around 1 (strong return time). We also present in this section the connection with quantum physics where these questions naturally occur. In section 3 we recall known results about the trace of powers of random unitary matrices and use them to estimate the probability of a real-part or weak return time. The rest of the paper is dedicated exclusively to the more interesting case of strong return times. Section 4 provides the connection with the asymptotic of large Toeplitz determinants that is useful for numerical simulations. The central aspects of the paper are then presented in section 5 where we use matrix models techniques to compute the probability of a strong return time. In particular our main results are:

• In theorem 5.1 we prove that the probability to have a strong return time satisfies the technical conditions required in the main result of [START_REF] Borot | Large-N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF] (briefly presented in proposition 5.1) so that the general form of the asymptotic expansion at large N is known explicitly.

• Since small times are shown to give a one-cut case (genus 0 spectral curve), we carry out the explicit computations of the first orders of 1 N 2 ln P N (t) up to O 1 N 6 for small times. These results are presented in theorem 5.3 and provide a new derivation (as well as the next orders) of Widom's results [START_REF] Widom | Strong Szegö limit theorem on circular arcs[END_REF] regarding the asymptotic of certain Toeplitz determinants.

• For integer times (i.e. integer powers of random unitary matrices) an additional symmetry allows the explicit computations of 1 N 2 ln P N (t) up to order O 1 N . Results are summarized in theorem 5.5

In section 6 we propose an Average Block Interaction Approximation to address the issue of non-integer times that turns to be very efficient numerically at leading order. Confrontation of our various results with numerical simulations are presented in section 7. Eventually in section 8 we propose two conjectures (conjectures 8.1 and 8.2) regarding the distribution of the first strong return time. More specifically we conjecture that the first strong return time obeys an exponential distribution whose parameter λ is also conjectured. We provide heuristic arguments as well as numerical simulations to support our conjectures.

2 Powers of random unitary matrices and definition of various return times

Recurrence time and random unitary matrices

In quantum physics one of the long standing problems is the issue of finding a proper interpretation of measurement that is compatible both with the probabilistic interpretation of quantum mechanics and the possibility to observe a unique outcome for each experiment. This issue has been discussed at length since the early days of quantum physics and we refer to [START_REF] Allahverdyan | Understanding quantum measurement from the solution of dynamical models[END_REF] for a recent review about the problems at stake. Although this paper is not directly aimed towards this problem, it deals with a question raised at the end of paragraph 6.1 (page 59) of [START_REF] Allahverdyan | Understanding quantum measurement from the solution of dynamical models[END_REF]. Indeed, in this paragraph, the authors explained that in the attempt to explain quantum measurements as a properly treated interaction between a quantum system S and a measurement apparatus A in a thermal bath, the discreteness of the spectrum of the operator m describing the pointer variable of the (finite but large) apparatus A gives rise to recurrence times when waiting long enough. Indeed if after a short time, the interaction between S and A destroys all correlations between the two systems (a phenomenon called truncation happening during a typical time τ trunc ), it happens that because the spectrum is discrete (and finite) there exists a Poincaré recurrence time τ rec for which the pointer variable of A comes back close to its original state. This means that even if the correlations between S and A are rapidly dephased (after a time of order τ trunc ), the memory of the correlations is not lost for good and emerges back periodically (with period τ rec ). In [START_REF] Allahverdyan | Understanding quantum measurement from the solution of dynamical models[END_REF] a simple model with independent and random eigenvalues for m is developed and leads to the conclusion that the recurrence time is inaccessibly large compared to the truncation time because the ratio is of order τrec τtrunc = π operator m as a random matrix whose eigenvalues are known to be correlated. Eventually, they claimed that using matrix models the recurrence time is expected to be smaller than the one found in their simple model but still sufficiently large so that is remains inaccessibly large in practice. In this spirit, this article can be seen as a confirmation of their argument in the case of a random unitary operator and direct comparison with their results is provided in section 8.2.

More specifically, we want to study the distribution of the eigenvalues of a unitary matrix U N of size N sampled uniformly according to the Haar measure of the unitary group U(N ).

Let us denote C the unit circle and u i = e iθ i ∈ C the i th eigenvalue of the matrix U N with θ i chosen in [-π, π]. It is well-known that the Haar measure on the unitary group is equivalent to the following distribution for the eigenvalues:

ZN = [-π,π] N dθ 1 . . . dθ N N i<j |e iθ i -e iθ j | 2 (2.1)
Equivalently, using the exponential variables u i = e iθ i the last partition function is equivalent to:

Z N = C N du 1 . . . du N ∆(u 1 , . . . , u N ) 2 e -N N k=1 ln u k = (-1) N (N -1) 2 i N ZN (2.2)
where ∆(u) = i<j (u i -u j ) 2 is the usual Vandermonde determinant. Indeed we have:

ZN = 1 i N C N du 1 u 1 . . . du N u N N i<j |u i -u j | 2 = 1 i N C N du 1 u 1 . . . du N u N N i<j (u i -u j ) 1 u i - 1 u j = 1 i N (-1) N (N -1) 2 C N du 1 u 1 . . . du N u N N i<j (u i -u j ) 2 u i u j = (-1) N (N +1) 2 i N C N du 1 . . . du N ∆(u 1 , . . . , u N ) 2 e -N N k=1 ln u k = (-1) N (N +1) 2 i N Z N (2.
3)

The partition function (2.2) can be seen as the diagonalized form of a matrix integral over a set of H N (Γ) of normal matrices where the eigenvalues are restricted to the contour Γ = C:

H N (Γ) = {M = ΛXΛ † ∈ M N (C) | ΛΛ † = I N and X = diag(x i )x i ∈ Γ} (2.4)
The measure on H N (Γ) is given by ∆(X) 2 dΛ dX with dΛ the Haar measure on U(N ) and

dX = N i=1
dx i is the curvilinear measure defined by: dx = f (t)dt for any parametrization x = f (t) , t ∈ [0, 1] of the curvilinear coordinate

In particular the measure dX is invariant under the choice of the parametrization f (t). The main interest of using (2.2) is that the loop equations and topological recursion method developed for hermitian matrix ensembles (i.e. Γ = R) is mostly applicable for normal matrices. The only difficulty lies in the determination of the spectral curve for which additional "hard edges" terms will appear. Note also that the potential in (2.2) is logarithmic but not polynomial as usually assumed in hermitian matrix models.

Time evolution and definition of various return times

Let (u i ) 1≤i≤N = e iθ i 1≤i≤N be the eigenvalues of the unitary matrix U N drawn from the measure given in (2.2) and define θ i (t) = tθ i corresponding to the eigenvalues of U t N . For convenience, we will always consider angles θ i (t) belonging to the interval [-π, π]. At time t = 0 it is obvious that all the eigenvalues are located at 1 (i.e. θ i = 0) and when t increases, the eigenvalues θ i (t) = tθ i start diluting around the unit circle. Note that since we are only interested in the eigenvalues we can extend the definition θ i (t) = tθ i for any real positive times t (and not only integers) even if the definition of U t N could be problematic. When waiting long enough, it may happen that all the eigenvalues come back close to 1 at the same time corresponding to a unitary matrix U t N close to the identity matrix I N . More precisely in this paper we are interested in the following quantities:

• Take > 0, we say that the eigenvalues have strongly returned around 1 at time t if

∀ 1 ≤ i ≤ N : u i (t) = e itθ i ∈ e -iπ , e iπ ⊂ C
where the last interval notation is to be understood as the arc-interval on the unit circle {e it , -π ≤ t ≤ π }.

• Denote U N (t) = U t N , we say that the eigenvalues have weakly returned around 1 at time t if:

S N (t) = 1 N Tr U N (t) = 1 N N i=1 e itθ i ≥ 1 -δ
where δ > 0 is a given fixed small constant. In the same way, we say that the eigenvalues have real-part returned at time t if:

R N (t) = 1 N N i=1 Re ( Tr U N (t)) = 1 N N i=1 cos (tθ i ) ≥ 1 -δ
This definition is the one used in [START_REF] Allahverdyan | Understanding quantum measurement from the solution of dynamical models[END_REF] to define the recurrence time.

• First return time: For a given sampling of U N , we define the first strong (resp. weak, real-part) return time to be the smallest time for which we have a new strong (resp. weak, real-part) return time:

τ N,strong = Inf t>0 (t strong return time and ∃ 0 < s < t | s is not a strong return time)
The last condition is necessary to avoid trivially that τ strong = 0. It can also be replaced by requiring t > 1 since we expect return times to be rare when N is large. Note in particular that the recurrence time τ rec introduced in [START_REF] Allahverdyan | Understanding quantum measurement from the solution of dynamical models[END_REF] corresponds to τ N,real-part .

We observe that we have the following inequalities:

τ weak (δ) ≤ τ real-part (δ) and τ real-part (δ = 1 -cos π ) ≤ τ strong ( ) (2.5)
In other words, it is more common to observe a weak return time than a real-part return than a strong return time after a trivial identification of the parameters used in the definitions. Eventually, we also introduce the probabilities of each return event:

P N,strong (t) = P(t is a strong return time)

P N,weak (t) = P(t is a weak return time) P N,real (t) = P(t is a real-part return time) (2.6) Note here that the probabilities are computed relatively to the Haar measure defined earlier for the matrix U N = U N (t = 1). Since we are only interested in the eigenvalues of the matrix, the problem can also be illustrated as follow: take N particles on the unit circle and make them rotate, starting at t = 0 at the point 1, with a constant velocity

v i = θ i ∈ [-π, π]
taken randomly according to the measure induced by (2.1). Then our problem consists in studying the times when all the particles come back close to the point 1 in the different ways presented above. For example the strong return time problem can be illustrated with the following picture: In this moving particles perspective, a more natural setting would correspond to draw the velocities uniformly and independently on [-π, π] rather than using the measure induced by (2.1). However since the two measures are different the results cannot directly be extended from one setting to the other. However we will see that for some quantities both problems can be thought as equivalent.

Weak and real-part return times

Eigenvalues of unitary matrices have been studied for quite a long time both at the macroscopic (distribution of the eigenvalues) and microscopic (local correlation between the eigenvalues) levels. Here we recall a few of the many results available in the literature. First since the works of Diaconis and Shashahani [START_REF] Diaconis | On the eigenvalues of random matrices[END_REF] we know that the following theorems hold: Proposition 3.1 Let U N be a sequence of N × N unitary matrices. Then the eigenvalues λ N = (λ 1,1 , . . . , λ 1,N ) tend almost surely to independent uniform variables on the circle.

Moreover we get from [START_REF] Duits | Powers of large random unitary matrices and Toeplitz determinants[END_REF]: Proposition 3.2 Let Z be a standard complex normal distributed random variable, and let k(N ) be an increasing sequence of integers. Then we have:

1 min(k(N ), N ) Tr U k(N ) N N →∞ → in distribution Z (3.1)
In particular there are two different regimes depending on whether the exponent increases more rapidly than the size of the matrix or not:

for l ≤ N : 1 √ l Tr U l N N →∞ → in distribution Z for l > N : 1 √ N Tr U l N N →∞ → in distribution Z (3.2)
Eventually, results regarding the eigenvalues for a finite N are known ( [START_REF] Rains | Images of eigenvalue distributions under power maps[END_REF][START_REF] Meckes | Spectral measures of powers of random matrices[END_REF]):

P N m,N - N θ 2π > x ≤ 2 exp -Min x 2 , x 2 4m 1 + ln m N (3.3)
Results 3.2 are interesting when studying the weak return times that involve the trace of the powers of the matrix U N . We will specify the corresponding results in the next section. However, these theorems do not provide direct information about the strong return times (which are rarer then weak return times) and do not specify the convergence rates or exact bounds at finite but large N . In particular, for a given time t, the size of the matrix will eventually becomes larger than t so that results from (3.3) will become useless.

Threshold results for weak and real-part return times at large N

We can use the two previous theorems in order to obtain some estimates of real-part and weak return times. Indeed, it is trivial that if the θ i 's are independent uniform variables on [-π, π] then the variables tθ i 's are also independent uniform variables on [-tπ, tπ]. Moreover an elementary computations shows that:

πt -πt cos(θ) dθ 2πt = sin πt πt = πt -πt e iθ dθ 2πt (3.4)
For convenience we define the cardinal sine function to be:

sin c (x) definition = sin x x (3.5)
Since it is known from 3.1 that our eigenvalues almost surely tend to independent uniform variables we get that ∀t > 0:

R N (t) where t c is the only positive time for which sin c (t c ) = 1 -δ when δ is small. This threshold result and in particular for t > t c means that for a large but finite N the weak or real-part return probability at a given time t should be extremely small and rapidly decreasing when N increases. In fact it seems natural to use theorem (3.2) to obtain an estimate of the return probability for a large but given N .

Estimates for the weak and real-part probability return for large N

Let us use 3.2 in order to get an estimate of R N (t) and S N (t) when t is a strictly positive integer. Theorem 3.2 implies that the leading order of P N,weak (t) as N → ∞ should be given by: for a given t > 0 : P N,weak (t) ∼

N →∞ P |Z| ≥ (1 -δ)N √ t for a sequence of times t N > N : P N,weak (t N ) ∼ N →∞ P |Z| ≥ (1 -δ) √ N
where Z is a complex standard normal distribution. The r.h.s. can be computed explicitly since the complex density of Z is given by f (z) = 1 π e -|z| 2 . Note that after integrating over the angle, |Z| has a Rayleigh distribution with parameter σ 2 = 1 2 whose density is given by g(r)dr = 2re -r 2 dr. A straightforward computation shows that:

P |Z| ≥ (1 -δ)N √ t = +∞ (1-δ)N √ t re -r 2 dr = e -(1-δ) 2 N 2 t (3.8)
In other words, we find that for a given large N :

for a given t > 0 :

P N,weak (t) ∼ N →∞ e -(1-δ) 2 N 2 t for a sequence of times t N ≥ N : P N,weak (t N ) ∼ N →∞ e -(1-δ) 2 N (3.9)
We can use similar techniques for the real-part return times. Note that the real part of Z is by definition a real normal distribution N (0, σ 2 = 1 2 ). We find:

P N,real (t) = P 1 N N i=1 Re Tr U (t) ≥ 1 -δ = P Re Tr U t N ≥ (1 -δ)N ∼ N →∞ P Re Z ≥ (1 -δ)N √ t = +∞ (1-δ)N √ t 1 √ π e -x 2 dx = 1 2 erfc (1 -δ)N √ t (3.10)
where erfc(x) is the complementary error function whose asymptotic expansion is given by: erfc

(x) = e -x 2 x √ π ∞ n=0 (-1) n (2n -1)!! (2x 2 ) n ∼ x→∞ e -x 2 x √ π (3.11)
Eventually we find:

for a given t > 0 :

P N,real (t) ∼ N →∞ √ t 2N (1 -δ) √ π e -(1-δ) 2 N 2 t
for a sequence of times t N ≥ N :

P N,real (t N ) ∼ N →∞ 1 2 √ N (1 -δ) √ π e -(1-δ) 2 N (3.12)
The general form of (3.9) and (3.12) for a given 0 < t < N is consistent with the expected form of the series expansion of the random matrix integrals like (2.2). Indeed, for random matrix integrals, it was proved in [START_REF] Hiai | A large deviation theorem for the empirical eigenvalue distribution of random unitary matrices[END_REF] that the leading order of the partition functions should be of order:

ln

Z N = N 2 F [-2] + O(ln N ) (3.13)
This result is also compatible with the full expansion provided by Borot, Guionnet and Kozlowski in the main theorem of [START_REF] Borot | Large-N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF] that we will present below.

The case of strong return times

It is tempting to apply the previous method to the strong return times and to consider the eigenvalues θ i 's as i.i.d. uniform variables on [-π, π] when N is large. Under this assumption a simple computation shows that we would get:

PN,strong (t) = I(t) 2π N (3.14)
where I(t) is a union of intervals characterizing the different possibilities that tθ ∈ [-π , π ] mod 2π. In particular the N dependence in this case is of the form PN,strong (t) = e N ln| I(t) 2π | which is no longer consistent with the expected e -N 2 F 0 (t) . As we will see in the next sections, this simple model does not give the right answer for P N,strong (t) that we find to be of the same form as P N,weak (t) and P N,real (t). The main difference with the other two cases is that for the real and weak return times, the presence of the trace (i.e. the sum over all eigenvalues) averages quantities sufficiently so that the simple independent and uniform variable model is a very good approximation. On the contrary strong return times do not involve averaged quantities and exhibit a peculiar behavior. It is the purpose of this article to study the large N asymptotic of these strong return times using first Toeplitz determinants and more importantly matrix models techniques.

4

Strong return times: Exact computations at finite N with Toeplitz determinants

Toeplitz determinants

In this section, we use the general theory of Toeplitz determinants (see [START_REF] Gray | Toeplitz and Circulant Matrices: A review[END_REF] for a simple review) to express the exact value of P N,strong (t). We start with the definition:

P N,strong (t) = 1 ZN Ĩ(t) N i<j |e iθ i -e iθ j | 2 dθ 1 . . . dθ N (4.1)
where (Cf (A.1) for the following computation) the normalization is given by: ZN

= [-π,π] N i<j |e iθ i -e iθ j | 2 dθ 1 . . . dθ N = (2π) N N ! (4.2)
The domain of integration Ĩ(t) is a union of intervals corresponding to the fact that e itθ i must belong to the segment of unit circle [e -iπ , e iπ ]. A straightforward computation shows that for any integer k ≥ 1:

∀t ∈ [2k + , 2(k + 1) -] : Ĩ(t) = k j=-k 2πj -π t , 2πj + π t ∀t ∈ [2k -, 2k + ] : Ĩ(t) = k-1 j=-k+1 2πj -π t , 2πj + π t ∪ 2πk -π t , π ∪ -π, - 2πk -π t (4.3)
In the first case, we get 2k + 1 complete intervals of length 2π t whereas in the first cases we get 2k -1 complete intervals and two pieces that reach the endpoints π and -π. We illustrate the situation with the following picture: Fig. 2: Evolution of the angles tθ i Integrals of type (4.1) are known as Toeplitz integrals. Toeplitz integrals can be rewritten as the determinant of a Toeplitz matrix (i.e. a matrix whose entries M i,j only depend on j -i) that can be obtained from the discrete Fourier transform of their symbol. We remind quickly the reader of the standard formalism. Define:

I f = [-π,π] N N i=1 f (e iθ i ) i<j |e iθ i -e iθ j | 2 dθ 1 . . . dθ N (4.4)
where f is an integrable function on the unit circle called the symbol of the integral I f . Then, the general theory of Toeplitz matrices implies that the integral I f can be expressed as the determinant of a Toeplitz matrix of size N × N whose entries are given by the Fourier coefficients of f . Namely:

I f = det (T i,j (f ) = t i-j ) 1≤i,j≤N with t k = 1 2π 2π 0
f (e iθ )e -ikθ dθ (4.5)

The matrix T N (f ) is a Toeplitz matrix and one can reconstruct f on the unit circle from its Fourier coefficients by the usual formula:

f (e iθ ) = N k=-N t k e ikθ (4.6)
making the correspondence one to one between a Toeplitz matrix and its symbol f on the unit circle. The general theory can first be applied for the normalizing constant ZN for which the symbol is f (x) = 1 and gives (2π) N N !. We rederive this result using permutation techniques in appendix. Moreover, the numerator in (4.1) is also a Toeplitz integral associated with the symbol f t (e iθ ) = 1 Ĩ(t) (θ) that is to say a piecewise constant function on the unit circle. Note that since Ĩ(t) can be decomposed into unions of disjoint intervals, one could rewrite f t as a product: f t (e iθ ) = j 1 Ĩj,t (θ). Unfortunately, the products of symbols is not known to behave nicely in the theory of Toeplitz determinants so we did not find a direct way to use this observation. However, since we know exactly the symbol f of our integral, we can rewrite our strong return probability as a Toeplitz determinant by computing the Fourier coefficients directly. The computation is relatively straightforward and we propose in appendix another equivalent way to recover the results by using permutations of S N . We find:

P N (t ∈ [2R + , 2(R + 1) -]) = det sin (j-i)(2R+1)π t sin (j-i)π t π(j -i) sin (j-i)π t 1≤i,j≤N P N (t ∈ [2R -, 2R + ]) = det δ j-i=0 - sin 2(j-i)πR t sin (1-)(j-i)π t π(j -i) sin (j-i)π t 1≤i,j≤N (4.7)
The problem is now to get the large N expansion of such determinants. It is also worth noticing that numerically such determinants are much easier to handle than the computation the initial integrals. Note also that in our case the Toeplitz matrices involved are symmetric.

Simplification at integer times

At integer times, the last formulas substantially simplify and we get:

P N (t = 2k + 1) = N det sin c (j -i)π 2k + 1 δ (j-i) ≡ 0 [2k+1] 1≤i,j≤N (4.8)
P N (t = 2k) = det δ j-i -(1 -) sin c (j -i)π(1 -) 2k δ (j-i) ≡ 0 [2k] 1≤i,j≤N (4.9)
Note that in the even integer case, the first delta function is necessary to force the diagonal entries of the matrix to be (and not -1 as the sin c would imply). As far as we know, there exists no general formula in the literature for such determinants although most of the entries are vanishing. Getting the asymptotic of such determinants is a difficult problem but as we will see in 5.4, matrix models techniques and the use of the specific rotation symmetry at these times give the possibility to actually compute the first terms of the series expansion when N → ∞.

Simple case when t > N

Let us consider an integer time t > N (or a sequence of times t N > N ). In this case the previous determinants simplify because the condition δ (j-i) ≡ 0 [2k] is only satisfied for i = j giving a diagonal determinant. It trivially gives:

∀ t > N ∈ N : P N (t) = N = e -N |ln | (4.10)
This formula is exact and recover the result of the simple independent and uniform model (3.14). We also note that this result has the same N dependence as the results (3.9) and (3.12). However in this case, the probability does no longer depend on t and remains constant. Result (4.10) is exact for integer times but can be adapted for non-integer times to get the leading order as N → ∞ but no longer an exact formula. Let us consider a time t > N (not necessarily in N). The exact return probability is given by the Toeplitz determinant (4.7).

In the first line of formula (4.7), we have: R = t- 2 while in the second we have R = t+ 2 . Observe also that in the determinants, the diagonal contributions are special. For large t, a simple expansion of the determinant shows that the diagonal terms represent the leading order since they are the only ones not to involve factors of 1 t . Hence we find that: The oscillations around the value N are easily recovered in numerical simulations for t ≈ N but are damped for very large time t >> N like 2 t .

∀ t > N : P N (t) ∼ N →∞ N 2 t t 2 N = e N ln(

General theorem for the asymptotic of a Toeplitz determinant

We now need to understand the situation when t ≤ N . Toeplitz determinants have been studied intensively [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 | Asymptotic behavior of block Toeplitz matrices and determinants[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] but as far as we know our general case is not part of the general results of the theory. Indeed, the first standard result about Toeplitz determinants is the following: Proposition 4.1 (Szegö) If f is strictly positive and continuous on the unit circle then we have:

1 N ln det T N N →∞ → 1 2π 2π 0 ln(f (λ))dλ (4.13)
Unfortunately this general result cannot apply to our case because our symbol f is not continuous nor strictly positive on the unit circle (it is the indicator function of a union of intervals). This is in fact not surprising because we expect the leading order in N to be of the form N 2 while the last theorem would only give a N dependence. Since the works of Szegö, many refinements of theorem 4.1 have been carried out. For example, Fisher-Hartwig singularities [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] can be introduced giving formulas for the asymptotic of a Toeplitz determinant whose symbol is a strictly positive function f having jump discontinuity or isolated zeros.

Unfortunately, these types of singularities are still not sufficient for our purpose since we have intervals over which our symbol function f is vanishing. Indeed, we stress here that the main difficulty does not lie in the discontinuity of the characteristic functions but rather in the fact that the functions identically vanishes on certain segments. To our knowledge, the closest known case regarding our problem has been carried out by Widom (See main theorem of [START_REF] Widom | Strong Szegö limit theorem on circular arcs[END_REF]) for a symbol f supported on the arc interval [α, 2π -α]. In this case the quantity 1

N 2 ln det T N (f ) is convergent in contrast with 4.1 for which 1 N ln det T N (f )
is convergent. This change of behavior in the asymptotic is typical of a symbol vanishing on a segment (and as we will see is logical from a matrix model perspective). In the next section, we will recover this result (and extend it to any interval [θ 0 , θ 1 ]) and provide a way to obtain all the subleading terms recursively by the use of the symplectic invariants developed by Eynard and Orantin in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]. Unfortunately the general case of a union of intervals is much more difficult even with the use of the symplectic invariants because the genus of the spectral curve increases rapidly. To our knowledge the problem of finding exact values or large N asymptotic of Toeplitz determinants whose symbol vanishes on multiple arc-intervals has never been answered in the literature. In the next section, we connect the strong return time probability problem with the series expansion of a suitable matrix model providing an algebraic characterization of the solution. In practice, this solution is mainly formal since when the spectral curve is of genus g we need to determine "optimal filling fractions" to get the complete spectral curve, a problem which is known to be transcendent and highly difficult to perform theoretically and even numerically. However, in the case of integer times, an additional rotation symmetry of the intervals allows a complete determination of the spectral curve so that we can get the asymptotic expansion of the corresponding Toeplitz determinant. Moreover, for general (non necessarily integer valued) t < N , we give an accurate way to approximate the leading order of the asymptotic of the Toeplitz determinant.

5 Probability of a strong return time using matrix models techniques

In the last section, we have formulated P N,strong (t) as a Toeplitz determinant with a piecewise symbol f vanishing on a union of intervals. Unfortunately, getting a large N asymptotic of such determinants is complicated and we propose to study the problem with random matrix models techniques. Let us observe that:

P N,strong (t) = 1 Z N I(t) N du 1 . . . du N ∆(u 1 , . . . , u N ) 2 e -N N k=1 ln u k (5.1)
(where Z N is given by (2.2)) is exactly an hermitian matrix integral in its diagonalized form and therefore can be tackled through matrix models techniques. Here the set I(t) is given by:

∀t ∈ [2k + , 2(k + 1) -] : I(t) = k j=-k [e i 2πj-π t , e i 2πj+π t ] ∀t ∈ [2k -, 2k + ] : I(t) = k-1 j=-k+1 [e i 2πj-π t , e i 2πj+π t ] ∪ [e i 2πk-π t , e iπ ] ∪ [e -iπ , e -i 2πk-π t ] (5.2)
In particular from matrix model considerations we expect P N,strong (t) to be of order exp(N 2 F [-2] (t)) when N is large. The last matrix integral exhibits two main differences from a traditional hermitian matrix model. First, the integration set is compact and thus in matrix models terminology our integral has several "hard edges" corresponding to the endpoints of the intervals of I(t). Secondly, the potential is only logarithmic but not polynomial (this would be impossible for hermitian matrix integrals because the integral would diverge).

We will see later that these issues must be taken into account in the analysis but do not prevent the usual theory to work properly at least at leading order. We also recall that the normalizing constant Z N is known explicitly:

Z N = (2π) N N !(-1) N (N +1) 2
i N and is expected to be subleading compared to the e -N 2 F [-2] term.

General form of the expansion at large N

Computations of integrals of the form:

ZN = ĨN i<j |e iθ i -e iθ j | 2 dθ 1 . . . dθ N (5.3)
given in (4.1) where Ĩ is a union of g + 1 intervals, falls into the category of mean field models with Coulomb gas interactions (i.e. a repulsive interaction equivalent to (θ i -θ j ) 2 when θ j → θ i ) whose large N asymptotic has been studied in [START_REF] Borot | Large-N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF]. In our context, we have β = 2 and only hard edges at the extremities of our intervals. The main theorem of [START_REF] Borot | Large-N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF] (theorem 1.1 page 6) states that:

Proposition 5.1 Main theorem of Borot-Guionnet-Kozlowski in [START_REF] Borot | Large-N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF] with β = 2 and hard edges. For a measure supported on g + 1 segments and under specific conditions (described and verified in appendix D) we have the following large N expansion:

ZN = N N + 1 4 (g+1) exp ∞ k=-1 N -2k F [2k]                  m≥0 l 1 ,...,lm≥1 k 1 ,...,km≥-1 m i=1 l i +2k i >0 N - m i=1 l i +2k i m! m i=1 F [2k i ],(l i ) l i ! • ∇ ⊗ m i=1 l i ν                  Θ -N 0 F [-2], (2)
(5.4)

where Θ is the Siegel theta function:

Θ γ (ν, T) = m∈Z g exp - 1 2 (m + γ) • T • (m + γ) + ν • (m + γ) (5.5)
and F

[2k],(l) are defined as the l th derivative of the coefficient F [2k] relatively to the filling fractions (See paragraph 8.1 of [START_REF] Borot | Large-N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF]).

Details of the quantities in the last proposition can be found in [START_REF] Borot | Large-N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF] but since we will not need the full expansion we only stress on the following points: [2k] is usually known in the literature as the perturbative (or sometimes called formal) series expansion of ln ZN . For β = 2, only even powers of N appears in the expansion and the numbers F [2k] are called symplectic invariants. They only depend on the spectral curve and can be computed using the so-called topological recursion described in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]. We stress here that there is a difference in the notation with [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] where powers are indexed F (k) with k ≥ 0. To avoid confusion we will keep the bracket notation F [k] for quantities referring to [START_REF] Borot | Large-N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF] formalism. For k ≥ 0, notation F (k) in Eynard-Orantin formalism [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] corresponds to -F [2k-2] in Borot-Guionnet-Kozlowski formalism.

• The exponential part ∞ k=-1 N -2k F
• For g = 0 (i.e. only one cut that is to say genus 0 spectral curves), the Theta function is identically zero and we recover a 1 N expansion of ln ZN . • Since the Siegel Theta function is evaluated at -N , we get a pseudo-periodic behavior at each order in 1 N including N 0 . • The g + 1 dimensional vector is the vector of optimal filling fractions i to spread over the various intervals. In general it is a solution of difficult equations and thus is usually impossible to determine analytically.

In our case, we are only interested in the first terms of the expansion and we get: Theorem 5.1 The probability to have a strong return time at times t ∈ N * (additional symmetry) or < t < 2 -(genus 0 curve) satisfies conditions of [START_REF] Borot | Large-N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF] for β = 2 and hard edges so in the g + 1 segments regime we have:

ln ZN = N 2 F [-2] + N ln N + 1 4 (g + 1) ln N + F [0] + ln Θ -N 0 F [-2], (2)
+ O 1 N (5.6)
where 2) is the matrix of the second derivatives of F [-2] relatively to the filling fractions, i.e. the g × g matrix whose entries are given by:

T = F [-2],(
T i,j = ∂ t 1 =0 ∂ t 2 =0 F [-2] +t 1 η i +t 2 η j with η k =   0, . . . , 0 k-1 , 1, 0 . . . , 0 g-k , -1   t ∈ R g+1
(5.7)

For general times t ∈ R * + , only the off-criticality condition of the equilibrium measure remains unproven though we conjecture that it should hold except at times 2k ± with k ∈ N * where we have an isolated point e iπ (for 2k -) or the splitting of an interval (2k + ).

The definitions and proofs of the technical conditions proposed by Borot, Guionnet and Kozlowski are presented in appendix D. Independently of their result, it is worth mentioning that the leading order N 2 F [-2] follows directly from potential theory (See [START_REF] Johansson | On fluctuations of eigenvalues of random hermitian matrices[END_REF] and [START_REF] Hiai | A large deviation theorem for the empirical eigenvalue distribution of random unitary matrices[END_REF]). In particular, this leading term does not require the off-criticality condition and is defined even for singular equilibrium measure. This is helpful to deal with singular times t = 2k ± for which conditions of Borot, Guionnet and Kozlowski may not apply. This is also in agreement with our numerical simulations (see section 7 and our Average Block Interaction Approximation (see section 6) where the leading order is a continuous function of t including at those times. Our strategy is now to compute the spectral curve of the problem in order to obtain the quantities F

[k] at least for k = -2 and k = -1. In particular for integer times, an additional symmetry gives the possibility to compute exactly even when we have multiple cuts. Since for non-integer times computing is analytically impossible, we propose instead an approximation to get back to genus 0 quantities.

Loop equations and spectral curve

We plan to apply the standard techniques of the topological recursion to obtain the asymptotic expansion of (5.1). We refer the reader to [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] for a review on this method. In random matrix models it is standard to introduce the following correlation functions:

Z N,t = I(t) N du 1 . . . du N ∆(u 1 , . . . , u N ) 2 e -N N k=1 ln u k W 1,t (x) = N i=1 1 x -u i t = 1 Z N,t I(t) N du 1 . . . du N N i=1 1 x -u i ∆(u 1 , . . . , u N ) 2 e -N N k=1 ln u k W p,t (x 1 , . . . , x p ) = = N i 1 ,...,ip=1 1 x 1 -u i 1 . . . 1 x p -u ip t c (5.8)
where the average t is taken relatively to the measure Z N,t . The subscript c means the connected or cumulant part. In general, the correlation functions are useful when taking the x → ∞ expansion that provide a generating series of the moments of the eigenvalues. These definitions recover the standard case of the one-matrix models integrals with hard edges at the endpoints of I(t). The first loop equation is therefore:

W 2 1,t (x) + W 2,t (x, x) - N x W 1,t (x) - N x N i=1 1 u i t = a k hard edges of I(t) α k x -a k (5.9)
where α k are undetermined constants depending on N and t (they can be written formally as a matrix integral of size N -1, Cf. (B.3)). The proof of (5.9) is standard in random matrix theory and is presented in appendix B for completeness. In random matrix theory, it is known that the so called 1 N expansion (also called perturbative expansion) of the correlation functions can be determined only with the spectral curve. The spectral curve is obtained from the projection of (5.9) at leading order. Since we only need the leading order of the expansion, results of [START_REF] Johansson | On fluctuations of eigenvalues of random hermitian matrices[END_REF] and [START_REF] Hiai | A large deviation theorem for the empirical eigenvalue distribution of random unitary matrices[END_REF] can be used even if the equilibrium measure is critical. We get:

ln Z N,t = F [-2] (t)N 2 + O(N ) W 1,t (x) = W [1] 1,t (x)N + O(1) W p,t (x 1 , . . . , x p ) = W [2-p] p,t (x 1 , . . . , x p )N 2-p + O(N 1-p ) (5.10)
As usual in matrix models integrals, it is convenient to translate the function W 1,t (x) to define:

y t (x) def = W [1] 1,t (x) - 1 2x (5.11)
Projecting the first loop equation (5.9) on N 2 gives the so-called spectral curve of the model:

y 2 t (x) = 1 4x 2 + c x + a k hard edges of I(t) α (0) k x -a k (5.12)
where c and α

k are defined as:

c = lim N →∞ 1 N N i=1 1 u i t , α (0)
k = lim N →∞ 1 N 2 α k (5.13)
Here, c, a k and α k as well as α

(0)
k are implicitly assumed to depend on t but for compactness we will omit to rewrite down the explicit t dependence.

From random matrix theory we know that the knowledge of this spectral curve should completely determine all functions

F [g] N,t , W [g] 1,t (x) and W [g]
p,t (x 1 , . . . , x p ) through the implementation of the topological recursion of this spectral curve as described in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]. Unfortunately, the knowledge of the spectral curve (5.12) is incomplete so far since we have not determined constants c and α

(0)
k in the formula. The general theory tells us that these constants can be fixed uniquely by specifying the following conditions:

• Asymptotic at infinity:

y t (x) = x→∞ 1 2x + c x 2 + O 1 x 3 ⇒ y 2 (x) = 1 4x 2 + c x 3 + O 1 x 4
(5.14)

• If we decompose the domain of integration in the form of union of arc intervals of the form

I(t) = q k=1 [a k , b k ] = q k=1
I k (t), the filling fractions 1 , . . . , q (that also depend on t) defined by:

I k (t) y t (x)dx = 2iπ k (5.15)
are fixed dynamically by minimizing the so-called chemical potential:

∀ 1 ≤ k ≤ q -1 : a k+1 b k y t (x)dx = 0 (5.16)
From these definitions it is clear that the filling fraction k represents the proportion of eigenvalues in the arc interval I k (t). The previous minimization physically represents the equilibrium achieved by transferring some eigenvalues from one interval to another (tunneling). Note here that all integrals are restricted on the unit circle and therefore can be also expressed on intervals for the corresponding angles.

In general, the first condition on the asymptotic is quite easy to deal with and provides simple relations between the unknown coefficients of the spectral curve. By contrast, the conditions on the filling fractions are known to be very difficult to handle on a computational perspective. The first condition on the asymptotic comes from the fact that by definition:

W 1,t (x) = N i=1 1 x + u i x 2 + O 1 x 3 t (5.17)
Now, since u i belongs to the unit circle we get that 1 u i = ūi and since the domain of integration I(t) is symmetric under complex conjugate, we get that

N i=1 u i t = N i=1 1 u i t ∈ R.
Projecting on the order of N gives the constant c introduced earlier and the asymptotic (5.14) for y(x). A straightforward computation combining (5.14) and (5.12) for y 2 (x), shows that the asymptotic is satisfied if and only if the following constraints are verified:

0 = c + a k hard edges of I(t) α (0) k (Order 1 x of y 2 t (x) ) 0 = a k hard edges of I(t) a k α (0) k (Order 1 x 2 of y 2 t (x) ) c = a k hard edges of I(t) a 2 k α (0) k (Order 1 x 3 of y 2 t (x) ) (5.18)
From (5.12), we know that the spectral curve should be of the form: ). Note that in the case of simple zeros, the equilibrium measure remains non-critical but the number of cuts is just strictly greater than d.

y t (x) = 1 2x P d (x) d k=1 (x -a k )(x -āk ) (5.
• P d (x) admits a double zero s 0 inside the cut [a k , b k ].
In that case the equilibrium measure is vanishing at this point and thus becomes critical. This phenomenon usually appears when two intervals are merging or one interval is splitting into two.

In conclusion, we see that the multiplicity and location of the zeros of P d (x) are crucial to determine the number of cuts and the criticality of the equilibrium measure. In theory imposing (5.16) or solving the equilibrium measure problem presented in (D.6) should determine completely y t (x) or its corresponding equilibrium density. However in practice solving analytically one of these problem is out of reach with current tools. We conjecture that for generic times t = 2N ± , y t (x) should only have double zeros outside the domain I(t) (one in each gap between intervals [a k , b k ]) and thus the corresponding equilibrium density is noncritical. As we will see in the following sections, there are two situations, namely when I(t) consists only in one interval (i.e. 0 < t < 2 -) or when t ∈ N * , where additional symmetry considerations are sufficient to determine the complete expression of y t (x). In particular we can verify our conjecture in those cases and the equilibrium density is non-critical. Indeed all angles θ i ∈ [-π, π] will be contracted to

tθ i ∈ [-tπ, tπ] ⊂ [-π , π ].
Let us now take t ∈ [ , 2 -] for which we have d = 1. In appendix C we prove that it gives a spectral curve of genus 0 for which computations can be carried out completely. In particular, conditions of proposition 5.1 are met. In appendix C we study the more general case (we will need it below) given by :

Z N, 0 ,θ 0 ,θ 1 = [e iθ 0 ,e iθ 1 ] N 0 du 1 . . . du N 0 ∆(u 1 , . . . , u N 0 ) 2 e -N N 0 i=1 ln u i (5.21)
for which we find the leading orders of ln Z N, 0 ,θ 0 ,θ 1 to be:

-F [-2] (N, 0 , θ 0 , θ 1 ) = -(2 0 -1) 2 ln sin |θ 1 -θ 0 | 4 -F [0] (N, 0 , θ 0 , θ 1 ) = - 1 24 ln 2 - 1 12 ln 0 + 1 24 ln tan |θ 1 -θ 0 | 4 - 1 8 ln sin |θ 1 -θ 0 | 2 -F [2] (N, 0 , θ 0 , θ 1 ) = 3 cos θ 1 -θ 0 2 -1 128 2 0 cos 2 θ 1 -θ 0 4 (5.22)
Note here that it is important to select the angles within the same interval of length 2π (for example [-π, π]). As expected from rotation invariance the result only involves the length of the interval θ 1 -θ 0 but not θ 0 or θ 1 directly and only real quantities are involved (since the partition function Z N, 0 ,θ 0 ,θ 1 is real, it is natural to check that its expansion are large N only involves real quantities). The case t ∈ [ , 2 -] is obtained from (5.22) by taking 0 = 1 and θ 1 = -θ 0 = π t . Moreover since there is only one cut, no quasi-periodic terms are present in (5.4) so we get:

Theorem 5.3 For ∀ t ∈ [ , 2
-], the probability to have a strong return time is given by:

1 N 2 ln P N,strong (t) + 1 N 2 ln((2π) N N !) = N →∞ ln sin π 2t + ln N N + 1 4 ln N N 2 + 1 24N 2 ln 2 sin 3 π t tan π 2t + 1 64N 4 1 -3 cos π t 1 + cos( π t ) + O 1 N 6 (5.23)
Details of the computations can be found in C. As claimed earlier, this result is also in agreement with Widom's result [START_REF] Widom | Strong Szegö limit theorem on circular arcs[END_REF] on the asymptotic expansion of a Toeplitz determinant with a symbol f supported on an arc interval [α, 2π -α]. In other words, we have just provided here a new derivation of Widom's result with topological recursion techniques as well as the next order in the expansion. In particular we observe that sin π 2t is smaller than 1 so that the probability is indeed exponentially small when N gets large. Moreover at t = we recover that P N,strong (t = ) ∼ 1 which is consistent with (5.20). We can verify numerically that the last formula is correct by computing the exact value of P N,strong (t) for low values of N using the Toeplitz determinant (4.7). We find the following picture: We clearly observe the convergence to formula (5.23). The simulations validate also the rate of convergence of order O 1 N 2 . Since the spectral curve is of genus 0 and explicit (See (C.14)) one could apply the topological recursion on this curve to obtain all F [2k] (N, 0 , θ 0 , θ 1 ) terms and therefore reconstruct up to any order the 1 N expansion of ln P N,strong (t).

Exact Results for integer times

At integer times the problem exhibits an additional discrete rotation symmetry giving extra information that is precious to solve the filling fractions problem. In particular it gives the opportunity to compute exactly the spectral curve and the optimal filling fractions.

Consequently we can then compute the symplectic invariants associated to this spectral curve giving the 1 N series expansion of the corresponding Toeplitz determinant. Let us illustrate what happens in the case of an odd integer t = 2k + 1. For an odd integer time, we have 2k + 1 intervals of size 2π 2k+1 that are distributed along the unit circle with a discrete rotation invariance of angle θ = 2π 2k+1 . For example for k = 4, the situation is illustrated by: Fig. 5: Situation at time t = 2k + 1 with k = 4. The 9 distinct arc intervals are represented in green.

As explained earlier the general form of the spectral curve coming from loop equations should be of the form: and P 2k+1 (x) is a monic polynomial of degree 2(2k + 1). The rotation symmetry of the problem should apply also to our spectral curve. Since the denominator is already invariant under such rotations, the numerator should also share this symmetry. This implies that:

y 2 2k+1 (x) = P 2k+1 (x) 4x 2 k j=-k (x -a j )(x -b j ) = P 2k+1 (x) 4x 2 (x 2k+1 -e -iπ ) (x 2k+1 -e iπ ) ( 5
P 2k+1 (x) = k j=-k (x -e 2iπj 2k+1 a)(x -e 2iπj 2k+1 b) (5.25)
Since the situation is also invariant by complex conjugation, it means that b = ā so that:

P 2k+1 (x) = k j=-k (x -e 2iπj 2k+1 a)(x -e 2iπj 2k+1 ā) (5.26)
We can then use a similar symmetry arguments as the one developed in (C.8). Let us take -k ≤ j ≤ k, we can compute (defining ũi = e -i 2πj 2k+1 u i ):

W 1 -e 2iπj 2k+1
= -e - = 0 (5.28)

Consequently we must have:

P 2k+1 (x) = k j=-k (x + e 2iπj 2k+1 ) 2 ) (5.29)
and the spectral curve reads:

y 2 2k+1 (x) = x 2k+1 + 1 2 4x 2 (x 2k+1 -e -iπ ) (x 2k+1 -e iπ ) (5.30)
In particular this shows that the equilibrium density at t k = 2k + 1 is supported on the the whole domain I(t k ) and is non-critical (zeros of y(x) fall into C \ I(t k )). This spectral curve looks very similar to the case of one interval with the exception that x is now replaced by x 2k+1 . A natural Zhukowsky parametrization would be:

x 2k+1 (z) = cos π + 1 2 sin π z -1 z y(z) = 1+x 2k+1 (z) x(z)(z+ 1 z ) sin π (5.31)
One must be careful here since x(z) is not well defined with the first equation. In fact around each interval, one should specify the choice of the branch used. In other words, the spectral curve (5.30) defines a Riemann surface that is obtained as the gluing of 2k + 1 identical genus 0 Riemann surfaces. Define X(z) = x 2k+1 (z) and Y (z) = y(z) (2k+1)x 2k (z) satisfying dx ∧ dy = dX ∧ dY . Then (5.30) is equivalent to 2k + 1 copies of the following spectral curve:

X(z) = cos π + 1 2 sin π z -1 z Y (z) = 1+X(z) (2k+1)X(z)(z+ 1 z ) sin π (5.32)
The spectral curve (5.32) defines a genus 0 Riemann surface of the same nature as the one studied in appendix C with identical filling fractions 0 = 1 2k+1 (as symmetry implies) and θ 1 = -θ 0 = π . Thus, we can directly use results given in appendix C (keeping in mind that we have 2k + 1 identical copies) to compute F [-2] and F [0] . Moreover in expansion (5.4) we have contributions from the Siegel Theta function terms with = 1 2k+1 , . . . , 1 2k+1 ∈ R 2k+1 . Eventually we get: Theorem 5.4 At integer times t = 2k + 1 with k ∈ N, the probability to have a strong return time is given by:

1 N 2 ln P N, strong (2k + 1) + 1 N 2 ln((2π) N N !) = 1 2k + 1 ln sin π 2 + ln N N + 2k + 1 4 ln N N 2 - 2k + 1 24N 2 2 ln(2k + 1) + ln 4 tan π 2 + 1 N 2 ln Θ -N 0 F [-2], (2)
+ o 1 N 2 with = 1 2k+1 , . . . , 1 2k+1 ∈ R 2k+1
We stress that this result is exact and that one could obtain the next orders F

[g] by applying the topological recursion to the curve (5.32) (see for example (C.48)). Moreover, these results are consistent with the simulations presented in section 7 as well as the computations for general times t / ∈ N presented in section 6. Nevertheless it is worth mentioning that our method does not provide a way to compute explicitly F [-2], (2) in theorem 5.4. Indeed, it would require to compute explicitly the energy functional (and its Hessian) or the spectral curve (and its first symplectic invariant) in a neighborhood of . Since our method relies mostly on the symmetry at it cannot be extended easily in a neighborhood where the symmetry of the polynomial P 2k+1 will be lost. On the bright side, we do not expect the support of the equilibrium density to change, only the position of the zeros of P 2k+1 should change a little to accommodate the change in filling fractions.

At even integer times t = 2k ∈ N, the situation is very similar and the same method can be applied. The spectral curve is given by:

y 2 2k (x) = x 2k + 1 2 4x 2 (x 2k -e -iπ ) (x 2k -e iπ ) (5.33)
Computations for this curve provides similar results as the odd case. Eventually we find that the last theorem can be extended to all integers:

Theorem 5.5 At any integer time t = k ∈ N * , the probability to have a strong return time is given by:

1 N 2 ln P N, strong (k) + 1 N 2 ln((2π) N N !) = 1 k ln sin π 2 + ln N N + k 4 ln N N 2 - k 24N 2 2 ln(k) + ln 4 tan π 2 + 1 N 2 ln Θ -N 0 F [-2], (2)
+ o 1 N 2 with = 1 k , . . . , 1 k ∈ R k .
We stress again that this result is exact and that one could apply the topological recursion to (5.30) and 5.1 to get the next orders of the series expansion of ln P N, strong (k).

Non-integer times: Average Block Interaction Approximation

When t increases, the number of cuts also increases and therefore we have to deal with a spectral curve with a strictly positive genus. We have seen in last section that for integer times t ∈ N then an additional symmetry (rotation) gives us the opportunity to compute the exact spectral curve and thus the series expansion of ln P N,strong (t). In theory, one could compute the corresponding F [-2] in a similar way than the one presented in the genus 0 case. As mentioned earlier, the main difficulty for non-zero genus curves is that some coefficients of the spectral curve are determined by conditions (5.16) on the filling fractions. In practice, these conditions are impossible to solve except in some rare and exceptional cases with additional symmetries. A possible way to avoid this difficulty may consist in computing F [-2] for a spectral curve with unknown coefficients 1 , . . . d for example using formulas developed in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]. However this seems difficult since at some point we need to integrate ydx, that is to say the Abel map of the Riemann surface between two points. Such computations are usually very hard for a given Riemann surface and become unreasonable when the Riemann surface depends on some formal parameters 1 , . . . , d . Even if this could be performed and one could get an explicit formula for F [-2] ( 1 , . . . , d ), we would still have to minimize this function relatively to ( 1 , . . . , d ) in order to obtain the correct dynamical filling fractions.

Since the function F [-2] ( 1 , . . . , d ) is not expected to be simple, this last step is likely to be impossible analytically. In order to avoid such complications we propose to obtain some suitable approximations of 1 N 2 ln P N (t) where only genus 0 computations appear. Of course the price to pay is that we will only get an approximation (suitable for computation) but not the exact value.

General Average Block Interaction Approximation

In order to have a strong-return time at time t, we need to draw the eigenvalues (θ 1 , . . . , θ N ) from a union of r(t) intervals (with r(t) increasing with t). Let us denote generically θ m i (t), θ M i (t) 1≤i≤r(t) the r(t) arc-intervals given by (4.3). We label them so that:

-π < θ m 1 (t) < θ M 1 (t) < θ m 2 (t) < • • • < θ m r (t) < θ M r (t) ≤ π.
The probability to have a strong return time at time t is thus given by:

P N,strong (t) = 1 Z N r(t) K=1 [θ m K (t),θ M K (t)] N dθ 1 . . . dθ N i<j e iθ i -e iθ j 2 (6.1)
The idea is then the following: we impose for each 1 ≤ i ≤ r(t) to have N i eigenvalues in the arc interval [θ m i (t), θ M i (t)] (To be precise since we need the number of eigenvalues to be an integer we should take N i eigenvalues in the interval [θ m i (t), θ M i (t)] but to avoid unnecessary notational complication we will omit the symbols) We denote by M ( 1 , . . . , r , N ) the number of ways to select the eigenvalues according to this restriction. Hence we define:

P N,strong (t, 1 , . . . , r ) = M ( 1 , . . . , r(t) , N ) PN,strong (t, 1 , . . . , r(t) ) M ( 1 , . . . , r(t) , N ) = N ! (N 1 )!(N 2 )! . . . (N r(t) )! (6.2)
Here PN,strong (t, 1 , . . . , r(t) ) is computed with the first θ 1 , . . . , θ

N 1 eigenvalues in [θ m 1 (t), θ M 1 (t)], the next θ N 1 +1 , . . . θ N ( 2 -1 ) eigenvalues in [θ m 2 (t), θ M 2 (t)
] and so on. Now computing exactly PN,strong (t, 1 , . . . , r(t) ) is still as hard as the initial problem therefore we need to find a good approximation to compute it in practice. We denote by θ i,K the i th eigenvalue in the K th interval. The idea is to decompose the interaction in the following form: first we have interactions between eigenvalues belonging to the same interval i<j e iθ i,K -e iθ j,K 2 . Then we also have interactions between eigenvalues belonging to different intervals:

k<k N k i=1 N k j=1
e iθ i,k -e iθ j,k 2 . Since the intervals are small compared to the distance between them (they are of length 2π t and distant from 2π t ), we can approximate in the last term every eigenvalue by a constant which we take as the center of its corresponding interval: Definition 6.1 We call the Average Block Interaction Approximation (ABIA) the following approximation:

k<k N k i=1 N k j=1 e iθ i,k -e iθ j,k 2 ≈ k<k |c k (t) -c k (t)| 2 k k N 2 (6.3) where c k (t) = e i θ m k (t)+θ M k (t) 2
is the center of the arc interval [e iθ m k (t), e iθ M k (t) ].

This assumption is equivalent to average the interactions between the eigenvalues of different intervals (hence the name average block interaction approximation). The initial problem is now significantly simplified. Indeed, we have N j eigenvalues in each interval submitted to the usual interaction between themselves (but not the other ones) and an additional potential term coming from the averaged interaction with the other intervals. More precisely the probability is given by:

P N,strong (t, 1 , . . . , r(t) ) ABIA ≈ M ( 1 , . . . , r(t) , N ) Z N   r(t) k<k =1 |c k (t) -c k (t)| 2 k k N 2     r(t) k=1 [θ m k (t),θ M k (t)] N k dθ 1 . . . dθ N k 1≤i<j≤N k e iθ i -e iθ j 2   = M ( 1 , . . . , r(t) , N ) Z N r k<k =1 |c k (t) -c k (t)| 2 k k N 2   r(t) k=1 [a k (t),b k (t)] N k du 1 . . . du N k ∆(u 1 , . . . , u N k ) 2 e -N k N k i=1 ln u i   (6.4)
where we have defined a ( t) = e iθ m k (t) ∈ C and b k (t) = e iθ M k (t) ∈ C. The last quantity is a one interval (genus 0 curve) computation with only N k eigenvalues and is computed in appendix C. At leading order in N → ∞ the quantities Z N and M ( 1 , . . . , r(t) , N ) do not contribute. Indeed a Stirling expansion in (6.2) with the fact that Z N = (2π) N N ! easily shows that the contributions are subleading compared to a e -N 2 F [-2] term. Hence we are left at leading order with:

1 N 2 ln P N,strong (t, 1 , . . . , r ) ≈ 2 k k r(t) k<k =1 ln |c k (t)-c k (t)|- r k=1 F [-2] (a k (t), b k (t), k , N )+O 1 N (6.5)
We remind the reader here that to obtain P N,strong (t) one needs to minimize the previous result relatively to the filling fractions. As we will see below, (6.5) is a quadratic form in the filling fractions ( 1 , . . . , r(t) ) and hence can easily be minimized by standard techniques.

Case when t ∈ [2k + , 2(k + 1) -]

In this case, we have a union of arc segments given by

I(t) = k j=-k [e i(2πj-π ) t

, e i(2πj+π ) t

]. There are 2k + 1 intervals:

a j = e i(2πj-π ) t and b j = e i(2πj+π ) t with j ∈ [-k, k].
In particular we observe here that for every interval we have θ 1 -θ 0 2 = π t which does not depend on the interval. Therefore from appendix C we get that at leading order each interval contributes with a factor:

F [-2] j = 2 j ln sin π 2t + O 1 N (6.6)
Then we take c j = e 2iπj t as the central point of each interval. Therefore we have a contribution of the form:

k j<j =-k e 2iπj t -e 2iπj t 2 j j N 2 = 2 N 2 k j<j =-k j j k j<j =-k 1 -cos 2π(j -j ) t N 2 j j = 4 N 2 k j<j =-k j j k j<j =-k sin 2 π(j -j ) t N 2 j j (6.7)
In the end we find (we denote = ( -k , . . . , k ) t ):

1 N 2 ln P N (t, ) ≈ ln sin π 2t k j=-k 2 j + 2 ln 2 k j<j =-k j j + k j<j =-k j j ln sin 2 π(j -j ) t + O 1 N (6.8)
We stress here that the filling fractions vector = ( -k , . . . , k ) should not be confused with the initial scalar parameter corresponding to the size of the window allowed around θ = 0 to declare a strong return time. As explained earlier we now need to optimize the previous function relatively to the filling fractions . Moreover from symmetry consideration (the problem is invariant by complex conjugation) we have -k = k so the dimension of the problem can be lowered from 2k + 1 down to k + 1. Let us introduce:

g( ) = ln sin π 2t k j=-k 2 j + k j<j =-k j j
2 ln 2 + ln sin 2 π(j -j ) t (6.9)

It is clearly a quadratic form in the filling fractions. However by definition of the filling fractions we have an affine constraint k j=-k j = 1 that from symmetry can be recast as

1 -0 -2 k j=1 j = 0.
Therefore we introduce a symmetric matrix A of size (k + 2) × (k + 2):

A 1,1 = 2 ln sin π 2t A 1,k+2 = A k+2,1 = -1 A j,k+2 = A k+2,j = -2 , ∀ 2 ≤ j ≤ k + 1 A 1,j = A j,1 = 4 ln 2 + 2 ln sin 2 π(j -1) t , ∀ 2 ≤ j ≤ k + 1 A j,j = 4 ln sin π 2t + 4 ln 2 + 2 ln sin 2 2π(j -1) t , ∀ 2 ≤ j ≤ k + 1 A k+2,k+2 = 0 A i,j = 8 ln 2 + 2 ln sin 2 π(j -i) t + 2 ln sin 2 π(j + i -2) t , ∀ 2 ≤ i = j ≤ k + 1 (6.10)
We also introduce the vector b of size k + 2 given by: b = (0, . . . , 0, -1) t (6.11)

and finally we introduce the vector of filling fractions and Lagrange multiplier x of size k + 2:

x = ( 0 , . . . , k , λ) t (6.12)

So that we get:

h( x) = g( ) + λ 1 - k j=-k j = 1 2 x t A x -x t b (6.13)
Our problem is now equivalent to find the extreme values of h( x). The general theory of quadratic forms gives:

x extr = A -1 b and h( x extr ) = - 1 2 b t A -1 b = - 1 2 A -1 k+2,k+2 (6.14)
Note that we can also compute x extr by ( x extr ) i = (-1) i+k+1 det( Ãi,k+2)

2 det A
where Ãi,k+2 is the submatrix of A with the i th line and (k + 2) th column removed. (The last expression is interesting because determinants are more stable and faster to compute than the inverse of a matrix) In the end we find:

∀t ∈ [2k + , 2(k + 1) -] : 1 N 2 ln P N,strong (t) ABIA ≈ - 1 2 b t A -1 b = - 1 2 A -1 k+2,k+2 = - 1 2 det Ãk+2,k+2 det A (6.15)
with A and b defined in (6.10) and (6.11) and Ãk+2,k+2 being the submatrix of A where we have excluded the last line and column.

We would like now to see if the former result is compatible with exact results obtained at integer times in theorem 5.4. At t = 2k + 1 ∈ N, symmetry considerations imply that all filling fractions are identical. In that specific case, our ABIA approximation gives:

1 N 2 ln P N,Strong (2k + 1) ABIA ≈ 1 2k + 1 ln sin π 2(2k + 1) + 2k 2k + 1 ln 2 + 1 (2k + 1) 2 2k j=1 j ln sin 2 πj 2k + 1 (6.16)
From formula (E.3) we can compute the last sum and we obtain:

1 N 2 ln P N,Strong (2k + 1) ABIA ≈ 1 2k + 1 ln (2k + 1) sin π 2(2k + 1) (6.17)
We note here that the ABIA fails to reproduce the exact formula of theorem 5.4 but is very close to it. Indeed, as soon as the approximation sin π 2(2k+1) ∼ π 2(2k+1) can be made we recover theorem 5.4. In particular this is the case for → 0 or k → ∞. Remark 6.1 A more drastic approximation would correspond to remove all interactions between different intervals instead of averaging them. In this case we would get:

∀ j ∈ [-k, k] : j = 1 2k + 1 and 1 N 2 ln P N,strong (t) ≈ 1 2k + 1 ln sin π 2t (6.18)
However this strong approximation provides a completely wrong answer since at large k (and hence large t) its behavior is very different from theorem 5.4. Consequently removing interactions between different intervals fails to reproduce the refined structure of the situation whereas averaging them is a much better approximation.

Case when t ∈ [2k -, 2k + ]

The case when t ∈ [2k -, 2k + ] is trickier than the previous one since we have to deal with an incomplete interval around e iπ = -1 which gives a contribution ln cos 2πk-π 2t . Indeed, we remind the reader that we must choose a suitable determination of the angles that is compatible with the cut in order to obtain F [-2] and thus it would be wrong to take . Then, the approximation of the interaction between intervals is given by:

θ 0 = -2πk-π
1 N 2 ln k j<j =-(k-1) |c j -c j | 2 = 2 ln 2 k-1 j<j =-(k-1) j j + k-1 j<j =-(k-1) j j ln sin 2 π(j -j ) t +2 ln 2 k-1 j=-(k-1) j k + k-1 j=-(k-1) j k ln cos 2 πj t (6.19)
In the end we find the following approximation:

1 N 2 ln P N (t, ) ABIA ≈ ln sin π 2t k-1 j=-(k-1) 2 j + 2 k ln cos πk t - π 2t + k-1 j<j =-(k-1) j j
2 ln 2 + ln sin 2 π(j -j ) t

+ k-1 j=-(k-1) j k 2 ln 2 + ln cos 2 πj t + O 1 N (6.20)
This is again a quadratic form in . Moreover, since the problem is invariant under complex conjugation, the filling fractions must satisfy -j = j for j ∈ [1, k -1]. This observation only leaves k + 1 unknown filling fractions. Moreover, the Lagrange multiplier looks like

λ 1 -0 -k -2 k-1 j=1 j .
Therefore we introduce the following (k + 2) × (k + 2) symmetric matrix A:

A k+2,k+2 = 0 A 1,k+2 = A k+2,1 = -1 A k+1,k+2 = A k+2,k+1 = -1 A i,k+2 = A k+2,i = -2 , ∀ 2 ≤ i ≤ k A 1,1 = 2 ln sin π 2t A 1,j = A j,1 = 4 ln 2 + 2 ln sin 2 π(j -1) t , ∀ 2 ≤ j ≤ k A 1,k+1 = A k+1,1 = 2 ln 2 A j,j = 4 ln sin π 2t + 4 ln 2 + 2 ln sin 2 2π(j -1) t , ∀ 2 ≤ j ≤ k A k+1,k+1 = 2 ln cos πk t - π 2t A i,j = 8 ln 2 + 2 ln sin 2 π(j -i) t + 2 ln sin 2 π(j + i -2) t , ∀ 2 ≤ i = j ≤ k A j,k+1 = 4 ln 2 + 2 ln cos 2 π(j -1) t , ∀ 2 ≤ j ≤ k (6.21)
The vector b of size k + 2 is now given by: b = (0, . . . , 0, -1) (6.22)

and finally we introduce the reduced (by symmetry) vector of filling fractions and Lagrange multiplier x of size k + 2:

x = ( 0 , . . . , k , λ) t (6.23) so that we get again the following function to extremize:

h( x) = 1 2 x t A x -x t b (6.24)
We find:

∀t ∈ [2k -, 2k + ] : 1 N 2 ln P N,strong (t) ABIA ≈ - 1 2 A -1 k+2,k+2 = - 1 2 det Ãk+2,k+2 det A (6.25)
where A and b are defined in (6.21) and (6.22) and Ãk+2,k+2 is the submatrix of A where we have excluded the last line and column. At time t = 2k, it is easy to observe by symmetry that all filling fractions are again identical (and take the value 1 2k ). Therefore we can compute explicitly an approximation of the probability for these times. We get from (6.20):

1 N 2 ln P N,strong (t = 2k) ABIA ≈ 1 2k ln sin π 4k + 2(2k -1)(k -1) (2k) 2 ln 2 + 1 (2k) 2 k-1 j<j =-(k-1) ln sin 2 π(j -j ) 2k + 2(2k -1) (2k) 2 ln 2 + 1 (2k) 2 (ln k -(2k -2) ln 2) = 1 2k ln sin π 4k + 2(2k -1)k -(2k -2) (2k) 2 ln 2 + 1 (2k) 2 ln k + 1 (2k) 2 k-1 j<j =-(k-1) ln sin 2 π(j -j ) 2k = 1 2k ln sin π 4k + 2(2k -1)k -(2k -2) (2k) 2 ln 2 + 1 (2k) 2 ln k + 1 (2k) 2 2k-1 j=1 (j -1) ln sin 2 πj 2k (6.26)
Eventually using (E.3) we find:

1 N 2 ln P N,strong (2k) ABIA ≈ 1 2k ln 2k sin π 4k (6.27)
As in the case of odd integer times, we do not recover completely the exact value of theorem 5.5 but in the limit where sin π 4k ∼ π 4k (i.e. → 0 or k → ∞) we recover it properly.

General remarks about our approximation

Our average block interaction approximation satisfies various interesting aspects:

• The filling fractions and the probability computed with our approximation are both continuous functions of t including at t = 2k -where the dimension of the matrix increases. This is coherent with the fact that we do not expect any specific singular behavior at t = 2k -.

• At t ∈ N we recover identically distributed filling fractions and our approximation almost gives the right formula. Indeed, our approximation gives 1 t ln t sin π 2t while the exact computation gives 1 t ln sin π 2 . As one can see, the difference is small and vanishes at first order in or t → ∞.

• At t = 2 -the probability is continuous and we know that the left limit is exact because it is a genus 0 computation. Therefore it means that our approximation is consistent at this point.

• Numerical simulations (Cf. 7) show very good agreement between the approximate curve and the exact values of the probability computed at finite large N by Toeplitz determinants.

• When t = N we recover that the leading order of the approximation is given by (4.10)

We also mention that our approximation could be improved in the following way: in order to carry out the approximation we had to take N i eigenvalues in each interval. However in the rest of the analysis, we considered i to be a real parameter in [0, 1] whereas it can of the form p N ∈ Q. This aspect is not problematic in the matrix model perspective developed in appendix C since it is valid for any i ∈ R but it may be important in the optimization process. However optimizing a quadratic form on a discrete space { p N , p ∈ [0, N ]} N seems much more complicated than on a continuous space. In particular this discrete/continuous issue will be important when t (which basically gives the number of cuts) is of order N . In that case, we can only have a few eigenvalues (and possibly none) in each interval and the analysis developed should be adapted.

Numerical simulations for strong return times

We want to compare numerical evaluations of the probability for finite values of N with the theoretical solutions or approximations developed in the previous sections. We take = 1 5 and we are able to compute exact values of the probability up to N = 9 for t < 8. We have several following possibilities:

• Exact N computations given by the determinantal formulas (4.7) that allow numerical computations of exact finite N values of the probability at any time t. Since the matrices involved are Toeplitz matrices, the determinant computation is much faster than usual. We computed the exact value up to N = 35.

• Specific exact cases for finite N at integer times given by the reduction of the determinantal formulas at these points in (4.9) and (4.8).

• The exact limiting curve for t < 2 -given by (5.23)

• For t > 2 -we can use our average block interaction approximation. The results are given in (6.15) and (6.25) and can be carried out numerically for small values of N .

Numerical simulations are presented here and show that our average block interaction approximation matches perfectly the large N exact computations. The simulations give the following results for = 1 5 (the choice of 1 5 is purely conventional but the picture is the same with other values of ): (t) computed exactly with the formula (A.3) for values of N from 2 to 35 (colored points). The black curve is our average block interaction approximation given by (6.15) and (6.25).

For a given N (i.e. one of the dot-curve presented above), we clearly observe that at time t = N the behavior of the curve changes and is given by (4.11). In particular it exhibits a quasi-periodic behavior of period 2 and the probability at integer times remains constant as claimed in (4.10).

First strong return time

In the previous section we have computed the probability that the eigenvalues come back into the interval [e -iπ , e iπ ] at time t. A natural question is now to find information about the distribution of the first strong return time τ N,strong when N becomes large. When N is large, having all the eigenvalues back into [e -iπ , e iπ ] is a very rare event since each eigenvalue only spends a proportion of time in the interval and that they are not synchronized a priori. Therefore, we expect the first-return time to be much larger than the size of the matrix N (a typical time is expected to be of order -N ). Moreover it is known from 3.3 that for times t > N the eigenvalues can be considered independent and uniformly distributed on the unit circle. Thus we expect the first-return time τ N,strong to be equivalent to the first-return time for the following problem: Definition 8.1 Denote a the distance to the nearest integer of a real number a. Take X 1 , . . . , X N i.i.d. uniform variables on [-1 2 , 1 2 ] (it would correspond to take X i = θ i 2π in the unitary matrix interpretation). Define

S t = 1 tX 1 ≤ 2 ,..., tX N ≤ 2 (8.1)
and the first return time τN, :

τN, = Inf t>0 {t such that S t = 1 and ∃ s < t such that S s = 0} (8.2)
(It is the exact analogue of our first strong return time in this i.i.d. setting)

Then we conjecture that:

Conjecture 8.1 In the large N limit we should have:

N τ N,strong 4 -(N -1) Law → N →∞ E(1) and N τN, 4 -(N -1) Law → N →∞ E(1) (8.3)
where E(1) stands for an exponential distribution with parameter λ = 1.

The heuristic idea behind the previous proposition is the following: First observe that for a given sampling of the initial angles θ i the only possible first return times are given by:

t i,k = 2πk |θ i | - 2 Sign(θ i ) , ∀ 1 ≤ i ≤ N and k ∈ N * (8.4)
Indeed a simple continuity argument of the trajectories implies that one of the eigenvalues must be reentering the interval at τ N,strong (otherwise there would be a immediate earlier return time). Then let us consider the variables:

S i,k = 1 t i,k θ 1 2π ≤ 2 ,..., t i,k θ N 2π ≤ 2 (8.5)
Since the θ i 's are independent and uniform variables, it is well known that for a given t > 1 the variables t θ i 2π are independent and uniformly distributed on 0, 1 2 . Hence the variables S i,k are Bernoulli variables with parameter p N, = N -1 (note that by construction of the times t i,k we know that the i th eigenvalue is in the success zone and thus only N -1 eigenvalues remain to be tested). Unfortunately, the variables S i,k are not independent because knowing the position of an eigenvalue at time t provides information on its location at other times. However we note that S i,k = 1 is a very rare event when N is large. Indeed, conjecture 8.1 states that the average strong return time is expected to be of order 4 N N -1 . Additionally, it is known that rotating around the circle is a mixing system (in the sense that a little change of the initial velocity makes the position very different after some time, thus making the position of a particle unpredictable after a sufficiently long time if the velocity is known approximately). Hence it seems reasonable that the dependence between the positions of an eigenvalue at time t 0 and at time t 1 should rapidly decrease when the times t 0 and t 1 become distant (typically ∆t = t 1 -t 0 should be greater than the so-called mixing time of the system). Consequently it seems reasonable that a memoryless exponential distribution should appear if the typical mixing time is lower than the typical strong-return time. Anyway, keeping the discussion rigorous, we have to compute the first success time for dependent Bernoulli variables corresponding to rare events. Recent results using Stein's method about the convergence of dependent Bernoulli variables for rare events towards Poisson variables (and exponential variables regarding the first success time) can be found in [START_REF] Arratia | Poisson Approximation and the Chen-Stein Method[END_REF][START_REF] Peköz | Stein's method for geometric approximation[END_REF][START_REF] Chen | Approximating dependent rare events[END_REF] (for example theorem 4.1 of [START_REF] Chen | Approximating dependent rare events[END_REF]). In particular these theorems explain why is seems natural to recover an exponential distribution in the limit. The last difficulty is to get an estimate of the average first success which we claim to be 4 N -(N -1) . Considering temporarily that the variables S i,k are independent trivially gives that the mean first success should be:

Mean First Success = 1 (∆ N, t) p N, (8.6)
where ∆ N, t is the average time between two consecutive possible first return times t i,k :

∆ N, t = E Min p∈N * (s p+1 -s p ) (8.7)
where the time sequence (s p ) p∈N * is a reordering in increasing order of the possible return times (t i,k ) i,k (with 1 ≤ i ≤ N and k ∈ N * ). When N increases the times (s p ) n∈N * become denser hence the typical difference between two of them decreases. A direct computation shows that we have ∆ N, t = 4 N hence giving conjecture 8.1. The proposition is supported by the following numerical simulations: Similar histograms can be obtained for the continuous i.i.d. case confirming that the two problems are intimately related: However we mention here that the results involving Stein's method found in [26, 27, 28] (and especially theorem 4.1 of [START_REF] Chen | Approximating dependent rare events[END_REF]) cannot be directly used in our situation. Indeed, these results are very efficient when the range of dependence between the Bernoulli variables is bounded (or is rapidly decreasing with the distance between the variables) which is not the case here. They are also very efficient in the case where we have an explicit strong bound on the level of dependence on the Bernoulli variables which would require a precise analysis here.

The discrete version of the first return time

We mention here a discrete version of the first return time problem by only considering integer powers of the matrix U N . Let us consider n 0 (N, ) the smaller integer n ≥ 1 for which U n N has all its eigenvalues in [e -iπ , e iπ ]. Since we also expect n 0 (N, ) to be much larger than N , we can recast this problem into the following one: Take X 1 , . . . , X N i.i.d. uniform variables on -1 2 , 1 2 and define ñ0 (N, ) the smallest integer n ≥ 1 such that:

∀ 1 ≤ i ≤ N : nX i ≤ 2 (8.8)
Then similar arguments as the one developed in the continuous time setting give:

Conjecture 8.2 In the large N limit we have: . This is because the difference between two consecutive testing times is now trivially 1 (we test all integer times). Moreover for a given integer time n we are no longer sure that one of the eigenvalues is located inside the success zone therefore giving p N, = N as the success probability for variables S n defined similarly to (8.5). The discrete version of the first-return time is supported by the following simulations: Note here that the simplified i.i.d. version of the discrete first-return time is a pure number theoretic problem. It can be recast into simultaneous Diophantine approximation of N random independent real numbers and may have some interests for people in this field. In this spirit it seems possible to tackle similar number theory inspired problems (like the lonely runner conjecture, see [START_REF] Goddyn | Tight Instances of the Lonely Runner[END_REF]) by extending them to unitary matrix models in first approximation (of course loosing the independence property since the eigenvalues are correlated) and apply known results there (like Toeplitz determinant or matrix models tools). Numerically we also observe that in the integer time setting always underestimate the value of λ while the continuous time setting always overestimate the value of λ. This probably means that the subleading term of the large N expansion is of different signs in the two settings.

n 0 (N, ) -N Law → N →∞ E (

Connection with the recurrence time quantum measurement

Let us try to interpret the results of the last section in relation with the recurrence time developed in [START_REF] Allahverdyan | Understanding quantum measurement from the solution of dynamical models[END_REF]. Conjecture 8.1 suggests that the average first strong return time is of order

τ strong ≈ 4 -(N -1)

N

when N is large. The first difficulty to compare this estimate with results from [START_REF] Allahverdyan | Understanding quantum measurement from the solution of dynamical models[END_REF] is that τ rec is related to the first real-part return time τ real-part and not the first strong return time τ strong . In this paper we focused on the strong return time because it is easier to handle from the matrix models perspective as well as arguments presented in the previous section. However it is easy to see that strong return times will be less frequent than real-part return times whose definition is weaker. Hence τ strong ≈ 4 -(N -1) N can be seen as an upper bound on the recurrence time. Let us now translate into our context the notation of [START_REF] Allahverdyan | Understanding quantum measurement from the solution of dynamical models[END_REF]. In our case, we have Q = N eigenvalues and the threshold required to have a strong (or real-part) return time is f = 1 -δ = cos π . Hence results from from [START_REF] Allahverdyan | Understanding quantum measurement from the solution of dynamical models[END_REF] states that we should have:

τ rec τ trunc = π √ 2e N cos π (8.10)
Our conjecture for τ strong supports such result. Indeed even we did not compute the truncation time in this paper, authors in [START_REF] Allahverdyan | Understanding quantum measurement from the solution of dynamical models[END_REF] claim that it should be of order τ trunc = 1 ∆ω where ∆ω stands for the average difference between two consecutive eigenvalues around the initial time. A trivial lower bound in our situation is to take τ trunc of order 1, though as mentioned in [START_REF] Allahverdyan | Understanding quantum measurement from the solution of dynamical models[END_REF] a more natural estimate is given by the uniform case: τ trunc = 1 ∆ω = N 2π . In this last case we find that:

τ rec τ trunc = 8π -(N -1) N 2 = 8π N 2 e (N -1)| ln | (8.11)
Thus we recover the exponential dependence in N of the ratio and only the dependence in the window parameter is modified from cos(π ) to | ln | (both being decreasing functions of in (0, 1]). In particular this makes the recurrence time rapidly inaccessible when N increases. We believe that the typical time τ real-part should not be very different from τ strong because of the following argument: Let us consider a first strong return time τ a in an arcinterval [e -ia , e ia ]. At τ a one of the eigenvalues (say the N th ) is re-entering the arc-interval and therefore its real-part is cos a. Since we expect τ a >> N , we know that the eigenvalues can be considered independent and uniformly distributed on the unit circle. Thus it seems reasonable at time τ a to consider that the N -1 remaining eigenvalues θ a i are i.i.d. on the arc interval [e -ia , e ia ]. We now consider the random variable:

S a = Re 1 N N i=1 e iθ i = cos a N + 1 N N -1 i=1 cos θ a i (8.12)
where θ a i are i.i.d. uniform random variables on [-a, a]. In the large N limit, we can approximate the last sum with a normal distribution using the central limit theorem. A simple computation shows that (we denote sin c x = sin x x ):

P (S > 1 -δ) = 1 -F 0 1 -δ -cos a N -N -1 N sin c a 1 2 + 1 2 sin c (2 ) -sin 2 c √ N -1 (8.13)
where F 0 is the cumulative distribution function of the standard normal distribution. Hence when N is large (we discard the cos a N term and take N -1 N ≈ 1) we observe the following threshold: when 1 -δ -sin c a < 0 we get that in the large N limit the eigenvalues have also real-part returned in the zone 1 -δ at time τ a while for 1 -δ -sin c a > 0 they have not real-part returned. Thus with this argument we can get an estimate of the first real-part return time with the knowledge of the first strong return time but evaluated with a suitable parameter a. We find that:

When N is large, a typical first real-part return time (corresponding to a recurrence time in [START_REF] Allahverdyan | Understanding quantum measurement from the solution of dynamical models[END_REF]) should be of order 4a -(N -1)

N

where a is the unique solution of 1 -δ -sin c a = 0 (i.e. for small values a ∼ √ 6δ)

In this approach we see that the inaccessibility claim in [START_REF] Allahverdyan | Understanding quantum measurement from the solution of dynamical models[END_REF] is still correct since the last estimate is still exponentially large in N .

Conclusion and outlooks

In this article we presented general results regarding the eigenvalues of powers of large random unitary matrices using different classical tools like Toeplitz determinants and matrix models. In particular we focused our analysis on determining various return times into the arc-interval [e -iπ , e iπ ] for a fixed . We supported the work with numerical simulations and recover some known results with different methods (like Widom's formula (5.23)). We also motivated our work with the recurrence time problem occurring in quantum measurement developed in [START_REF] Allahverdyan | Understanding quantum measurement from the solution of dynamical models[END_REF]. An interesting extension of our work could be to vary the size of the return zone according to the power n (or t in the continuous setting) of the matrix or the size of the matrix N . We believe that the methods presented here could be used to deal with these problems. We would be very happy to work on these matters with anyone interested in these problems.

A Recovering the Toeplitz determinants with permutations

In this section, we derive the standard results of Toeplitz integrals and determinants for our case. We use a permutation approach that is different from the standard Fourier approach known for the general case.

A.1 Computing the normalization factors

We have the following:

ZN = [-π,π] N i<j |e iθ i -e iθ j | 2 dθ 1 . . . dθ N = [-π,π] N dθ 1 . . . dθ N det e i(j-1)θ i i,j=1..N det e -i(j-1)θ i i,j=1..N = [-π,π] N dθ 1 . . . dθ N σ,τ ∈S N (-1) |σ|+|τ | N k=1 e i(σ(k)-1)θ k N r=1
e -i(τ (r)-1)θr

B First loop equation

The first loop equation presented in (5.9) is obtained from the following observation. Define:

A t (x) = 1 Z N I(t) N du 1 . . . du N N j=1 d du j 1 x -u j ∆(u 1 , . . . , u N ) 2 e -N N i=1 ln u i (B.1)
Then we can compute this quantity in two different ways. First, since it is a total derivative, it can be directly integrated and since we have hard edges, we get a contribution of the form:

A t (x) = a k hard edge of I(t) α k x -a k (B.2)
where we can write explicitly:

α k = ± 1 Z N I(t) N -1 i =k du i ∆ N -1 (u 1 , . . . , u k-1 , u k+1 , . . . , u N ) 2 i =k (u i -a k ) 2 e -N ln a k e -N N i =k ln u i (B.3)
Here the sign is +1 if the hard edge is an upper bound of an interval or is -1 if the hard edge is a lower bound of an interval. The second way to evaluate A t is to make the derivative acts on every term and identify the terms with their corresponding correlation functions. We find:

A t (x) = k=1 1 (x -u k ) 2 t + 2 i,j =i 1 x -u i 1 x -u j t -N N j=1 1 (x -u j )u j t = -W 1,t (x) + W 2,t (x, x) + W 1,t (x) 2 + W 1,t (x) - N x W 1,t (x) - N x N j=1 1 u j t = W 2,t (x, x) + W 1,t (x) 2 - N x W 1,t (x) - N x N j=1 1 u j t (B.4)
thus proving (5.9).

C The one-cut case with two arbitrary points on the unit circle and an arbitrary filling fraction C.1 Spectral curve

In this section, we are interested in computing:

I(a, b, 0 ) = [a,b] N 0 du 1 . . . du N 0 ∆(u 1 , . . . , u N 0 ) 2 e -N 0 N 0 i=1 ln u i (C.1)
We also define the correlation functions as well as their 1 N expansion similarly to (5.8) and (5.10). (but with all sums going from 1 to 0 N instead of N ). We restrict to the case where the endpoints are chosen on the unit circle and write a = e iθ 0 and b = e iθ 1 with θ 0 and θ 1 chosen in [-π, π]. Let us study:

A t,a,b, 0 (x) = 1 I(a, b, 0 ) [a,b] N 0 du 1 . . . du N 0 N 0 j=1 d du j   1 x -u j ∆(u 1 , . . . , u N 0 ) 2 e -N 0 N 0 i=1 ln u i   (C.
2) Since it is a total derivative, it can be evaluated at the endpoints of the integral and we get:

A t,a,b, 0 (x) = µ x -a + ν x -b (C.3)
Alternatively, we can apply the derivative to each term and we get:

A t,a,b, 0 (x) = N 0 i=1 1 (x -u i ) 2 + 2 N 0 i,j =i=1 1 x -u i 1 x -u j -N 0 N 0 i=1 1 (x -i )u i = -W 1,a,b, 0 (x) + W 2,a,b, 0 (x, x) + W 1,a,b, 0 (x) 2 + W 1,a,b, 0 (x) - N 0 x W 1,a,b, 0 (x) - N 0 x N j=1 1 u j = W 2,a,b, 0 (x, x) + W 1,a,b, 0 (x) 2 - N 0 x W 1,a,b, 0 (x) - N c x (C.4)
We thus introduce the following shift:

y(x) = W [1]
1,a,b, 0 (x) -0 2x (C.5) so that the spectral curve is given by:

y 2 (x) = A t,a,b, 0 (x) + c x + 2 0 4x 2 = 2 0 4x 2 + c x + µ x -a + ν x -b (C.6)
By definition of W 1,a,b, 0 (x) the asymptotic at infinity must be of the form:

W [1] 1,a,b, 0 ∼ 0 x + O 1 x 2 ⇒ y(x) ∼ 0 2x + O 1 x 2 (C.7)
So far we get:

y 2 (x) = 2 0 (x -α)(x -γ) 4x 2 (x -a)(x -b) = 2 0 4x 2 + c x + µ x -a + ν x -b (C.8)
We want to prove that α = γ, i.e. that the spectral curve does not admit simple zeros. Indeed, simple zeros would imply that the equilibrium density is supported on two intervals included in [θ 0 , θ 1 ] on not the whole [θ 0 , θ 1 ] segment. Let us compute (we define ũi = e -i θ 0 +θ 1 2 u i ):

W 1,a,b, 0 -e i θ 0 +θ 1 2 = -e -i θ 0 +θ 1 2 N 0 i=1 1 1 + e -i θ 0 +θ 1 2 u i = -e -i θ 0 +θ 1 2 N 0 i=1 1 1 + ũi ũi ∈ e -i θ 1 -θ 0 2 ,e i θ 1 -θ 0 2 (C.9)
If we translate the origin of angles at e i θ 0 +θ 1 2 we get by symmetry under complex conjugation (with 1 u = ū on the unit circle) that:

N 0 i=1 1 1 + ũi ũi ∈ e -i θ 1 -θ 0 2 ,e i θ 1 -θ 0 2 = N 0 i=1 1 1 + 1 ũi ũi ∈ e -i θ 1 -θ 0 2 ,e i θ 1 -θ 0 2 (C.10) Hence: 2W 1,a,b, 0 -e i θ 0 +θ 1 2 = -e -i θ 0 +θ 1 2 N 0 i=1 1 1 + ũi + ũi ũi + 1 ũi ∈ e -i θ 1 -θ 0 2 ,e i θ 1 -θ 0 2 = -N 0 e -i θ 0 +θ 1 2 (C.11)
In particular we must have:

W [1] 1,a,b, 0 -e i θ 0 +θ 1 2 = - 0 2 e -i θ 0 +θ 1 2 ⇒ y -e i θ 0 +θ 1 2 = 0 (C.12)
Inserting this result into (C.8) gives:

y 2 (x) = 2 0 (x + e i θ 0 +θ 1 2 ) 2 4x 2 (x -a)(x -b) = 2 0 4x 2 + c x + µ x -a + ν x -b (C.13)
or equivalently:

y(x) = 0 (x -α) 2x (x -a)(x -b)
where α = -e i (θ 0 +θ 1 )

(C.14)

This shows that the equilibrium measure is supported on the whole interval [θ 0 , θ 1 ] and is non-critical (the double zero of y(x) is outside [θ 0 , θ 1 ]). We can then use a Zhukowsky parametrization:

x(z) = a + b 2 + (b -a) 4i z - 1 z = e i (θ 0 +θ 1 ) 2 cos θ 1 -θ 0 2 + 1 2 sin θ 1 -θ 0 2 z - 1 z y(z) = 0 x(z)e -i (θ 0 +θ 1 ) 2 + 1 x(z) sin θ 1 -θ 0 2 (z + 1 z ) = 0 1 + cos θ 1 -θ 0 2 + 1 2 sin θ 1 -θ 0 2 z -1 z x(z) sin θ 1 -θ 0 2 (z + 1 z ) (C.15)
In particular the one form ydx is:

ydx(z) = 0 x(z) + e i θ 0 +θ 1 2 2zx(z) dz = 0 1 2z + 1 2(z -z + ) - 1 2(z -z -) (C.16)
We observe again that the branchpoints are not singularities of ydx. Eventually we denote:

z + = 1 -cos θ 1 -θ 0
the two zeros of x(z). Moreover we have:

dx(z) x(z) = - 1 z + 1 z -z + + 1 z -z- (C.18)
Eventually we observe that the factor e i θ 0 +θ 1 2 can be factorized using symplectic invariance. Indeed, if we take x = e -i θ 0 +θ 1 2 x and ỹ = e i θ 0 +θ 1 2 y then the spectral curve reads:

x(z) = cos θ 1 -θ 0 2 + 1 2 sin θ 1 -θ 0 2 z - 1 z ỹ(z) = 0 (x(z) + 1) x(z) sin θ 1 -θ 0 2 (z + 1 z ) (C.19)
from which we can deduce that all F [g] will only depend on the difference θ 1 -θ 0 but not directly from (θ 0 , θ 1 ). It would be interesting to see if this spectral curve comes from a specific integrable system.

C.2 Computation of F [-2]

We want to determine F [-2] (noted -F (0) in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]) associated to the spectral curve (C.15). In order to do so, we need to identify the type of singularities arising in the one-form ydx. Here we have 4 singularities:

• z = 0 is a pole of x(z) of order 1. In [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF], it corresponds to a type 1 singularity

• z = ∞ is a pole of x(z) of order 1. In [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF], it corresponds to a type 1 singularity

• z + and z -that satisfy x(z ± ) = 0 are simple poles of ydx without being singularities of x(z) nor y(z). In [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF], it corresponds to a type 2 singularity.

Computation of the temperatures: The first step to determine F [0] is to compute the so-called temperature of the singularities defined as t α = Res z→α ydx. We find:

Singularity Temperature t 0 0 2 t ∞ -0 2 t + 0 2 t - -0 2
For example, we get for z + the following:

t + = Res z→z + ydx = Res z→z + (x(z) + 1) 0 2zx(z) dz = 0 2 Res z→z + dz zx(z) = 0 2z + x (z + ) = 0 2 (C.20)
We observe that the sum of the temperatures equals 0 as expected from the theory.

Computation around z = 0: We now need to define some local coordinate to compute the contribution of each singularity. We start with z = 0. It is a pole of order 1 of x(z) and therefore we need to choose in Eynard-Orantin language z 0 (z) = x(z). The potential V 0 (z) should be given by:

V 0 (z) = Res q→0 y(q)dx(q) ln 1 - x(z) x(q) (C.21)
But since ydx(q) = 1 2q + O(q) then we find:

V 0 (z) = 0 (C.22)
Therefore we must compute:

µ 0 = K 0 ydx -dV 0 + t 0 dx(z) x(z) +V 0 (K)-t 0 ln x(K) = K 0 ydx + 0 2 dx(z)
x(z) -

ydx(z) + t 0 dx(z) x(z) = ydx(z) + 0 2 dx(z) x(z) = 0 z -z + (C.24)
Therefore the integral gives:

µ 0 = 0 ln (K -z + ) -0 ln(-z + ) - 0 2 ln x(K) (C.25)
Computation around z = ∞: It is a pole of order 1 of x(z) and therefore we need to choose in Eynard-Orantin language z 0 (z) = x(z). The potential V 0 (z) should be given by:

V ∞ (z) = Res q→∞ y(q)dx(q) ln 1 - x(z) x(q) (C.26)
But since ydx(q) = 0 2q + O(q) then we find:

V ∞ (z) = 0 (C.27)
A straightforward computation shows that:

ydx(z) + t ∞ dx(z) x(z) = ydx(z) - 0 2 dx(z) x(z) = 0 z - 0 z -z - (C.28)
Therefore the integral gives:

µ ∞ = 0 ln K -0 ln (K -z -) + 0 2 ln x(K) (C.29)
Computation around z = z + : This time we have to deal with a zero of x(z) which is not a branchpoint. The local parameter in Eynard-Orantin formalism is z +

(z) = 1 x(z)-x(z + ) = 1 x(z)
since x(z + ) = 0 by definition. Therefore we should compute:

V + (z) = Res q→z + y(q)dx(q) ln 1 - x(q) x(z) (C.30)
We have y(q)dx(q) of the form α q-z + while ln 1 -x(q)

x(z)

∼ q→z + -x (z + ) x(z) (q -z + ) + O ((q -z + ) 2
). Therefore there is no residue here since the function is regular at q = z + and we find:

V + (z) = 0 (C.31) Now note that d dz ( 1 x(z) ) 1 x(z) = -x (z)
x(z) so that it gives a negative sign so that the integrand of Eynard-Orantin becomes:

ydx(z) + t + dz + (z) z + (z) = ydx(z) - 0 2 x (z) x(z) = 0 z - 0 z -z - (C.32)
and therefore we find:

µ + = 0 ln K -0 ln(K -z -) -0 ln(z + ) + 0 ln(z + -z -) + x(z) and as before we find:

V -(z) = Res q→z - y(q)dx(q) ln 1 -x(q) x(z) = 0 (C.34)

Now since the temperature of the pole is the opposite of z + , we find:

ydx(z) + t - dz -(z) z -(z) = ydx(z) + 0 2 x (z) x(z) = 0 z -z + (C.35)
Therefore we find:

µ -= 0 ln(K -z + ) -0 ln(z --z + ) - 0 2 ln x(K) (C.36)
Computation of F [-2] : We can now insert all the results in order to compute F [-2] . It is given in general by the formula:

-F [-2] = 1 2 α Res α V α ydx + 1 2 α t α µ α - 1 4iπ i A i ydx B i ydx (C.37)
Here since the curve is of genus 0 there are no cycles so the last term is 0. Moreover, all potential V α are 0 so we are left with the central term. Eynard-Orantin theory proves that F [-2] should be independent of the choice of K. We can check this fact here by the explicit computation: Indeed we have: In this section we compute F [0] (noted -F (1) in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]) associated to the spectral curve (C.15). Since our curve is of genus 0, we can compute F [0] (see [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]):

F [-2] θ 0 ,θ 1 , 0 = -
t 0 µ 0 + t ∞ µ ∞ + t + µ + + t -µ -= 2
-F [0] = -1 2 ln τ Bx -1 24 ln a i branchpoint y (a i ) with y (a i ) = dy(a i ) dz i (a i ) and z i (p) = x(p)

x(p) -x(a i ) (C.40) Note here a slight difference in our formula compared to definition 4.3 of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]. Indeed, in definition 4.3 of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF], the local coordinate is z i (p) = x(p) -x(a i ) because the branchpoints appear traditionally in matrix models with polynomial potentials at the numerator of y(x). However in our case, the branchpoints are poles of y(z) (and not zeros) and the standard definition z i (p) =

x(p) -x(a i ) is not consistent with definition of the moduli of poles defined in section 3.4.2 of the same paper [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]. In the case where a branchpoint is a pole of y(z) the local coordinate must be taken as z i (p) =

x(p) √ x(p)-x(a i )

. In particular the numerator is necessary to maintain F [0] invariant under symplectic transformations. In our spectral curve, we have 2 branchpoints located at z = ±i. We have x(i) = e iθ 1 while x(-i) = e iθ 0 and we remind the reader that by definition of a branchpoint x (i) = x (-i) = 0. It is then a straightforward computation from (C.15) to compute x (i) = -1 2 e iθ 1 -e iθ 0 and x (-i) = 1 2 e iθ 1 -e iθ 0 . In the end we have: 1 + e i θ 1 -θ 0 2 8x(i)x(-i) x (i)x (-i)

dz z=i (z) = - √ 2x(i) x (i)(z -i) 2 + O 1 z -i
and the cardinal sine function is sin c (x) = sin x x . Note that in our case, the so-called 2linear potential T is indeed symmetric and does not have any large N expansion (thus only T [0] (x 1 , x 2 ) = T (x 1 , x 2 ) is present in [START_REF] Borot | Large-N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF]). Moreover, we only have to deal with compact subsets of R since A ⊂ [-π, π]. Hence confinement issues are trivial in our case. We remind now the reader of the technical conditions to satisfy in order to have proposition 5.1: Proposition D.1 If the following conditions are satisfied:

• A is the union of g + 1 disjoint segments A i not reduced to a point: A = g i=0 A i

• (Regularity) T is a continuous bounded function on A 2

• (Uniqueness of the minimum) The energy functional E defined by:

E[µ] = - A 2 T (x 1 , x 2 ) 2 + β 2 ln |x 1 -x 2 | dµ(x 1 )dµ(x 2 ) (D.4)
has a unique minimum on the space of probability measure on A noted M 1 (A) in the unconstrained model.

• T is holomorphic in a neighborhood of A 2

• (Non-criticality) The equilibrium measure µ eq (i.e. the unique measure minimizing E) is off-critical in the sense that its support S is a union of g 0 + 1 segments not reduced to a point and: We now discuss the validity of the previous conditions in our situation:

• First condition is obvious for times t / ∈ {2k -, k ∈ N * }. Indeed except at those times we always have disjoint intervals described in (5.2) that are not reduced to a single points. At times t k = 2k -with k ∈ N * , we have an isolated point at θ = π so the first condition fails.

• The regularity condition and the holomorphic condition for the function T (x, y) are obvious since T (x, y) = sin 2 c

x-y 2

Eventually we see that only the conditions regarding the uniqueness of the minimum and the non-criticality are non trivial in our situation. In our case, the energy functional looks like: Hence following [START_REF] Borot | Abstract loop equations, topological recursion, and applications[END_REF], a sufficient condition to obtain uniqueness of the minimum is to show that L(s) is negative so that the interaction is strictly convex on R and thus also on A. This can be done by direct computation of the Fourier transform L(s) =

E[µ] = - 1

-2π

+∞ m=1 1 m δ(s -m 2π ) + δ(s + m 2π ) or using equations A.3 and A.4 of [START_REF] Borot | Abstract loop equations, topological recursion, and applications[END_REF]. We stress here that this result is valid at any time t > 0 in our situation. Eventually the last remaining point deals with the off-criticality condition. Off-criticality is never easy to prove and in our case we expect it to fail for times of the form t k = 2N ± . Indeed at those times we have:

• Times of the form t = 2k -with k ∈ N * . In that case, there is an isolated point at ±π in the domain A = I(t). Hence the first condition of D.1 fails.

• Times of the form t = 2k + with k ∈ N * . In that case, an interval is splitting into two intervals (at θ = ±π). Such situations are known to be problematic and in particular we expect the equilibrium density to exhibit a singular behavior at those edge points. We expect the edge to become soft and no longer hard (the polynomial P d (x) in the spectral curve is expected to vanish there) making (D.5) fail.

Nevertheless there are special cases where the direct computation of the spectral curve allows us to prove the off-criticality condition. In particular when the domain A = [θ 0 , θ 1 ] is supported on a unique interval, we proved in (C.13) (using symmetry arguments) that the equilibrium measure is supported on the whole interval [θ 0 , θ 1 ] and that the equilibrium measure is non-critical. Using similar arguments, we proved that for integer times t ∈ N * , the problem can be recast into several copies of the previous problem hence giving an explicit and off-critical spectral curve corresponding to an equilibrium measure supported on the whole domain A t . However for general times t > 2 -and t / ∈ N, off-criticality remains a challenge since the initial domain is supported on several intervals and no symmetry arguments can be used. On the energy functional perspective, the situation seems difficult too since moving eigenvalues from a place to another increases the energy contribution of some of the eigenvalues while it decreases the energy contribution of others. In the end, it is even unclear why a gap in one of the intervals [a i , b i ] of the domain A t may not be more beneficial to the energy than an equilibrium density supported on the whole domain A t . Additionally Jensen inequalities seem difficult to use in (D.6) because the interaction -ln | sin ∆θ 2 | is repulsive at short distances (∆θ ≤ π) but attractive at long distances (∆θ ≥ π).

E A few computations with trigonometric sums

It is known that the following identities hold [START_REF] Weisstein | Sine." from MathWorld-A Wolfram Web Resource[END_REF]:

Fig. 1 :

 1 Fig. 1: Illustration of the strong return problem for N = 4. Dashed lines give a strong return time interval.

Proposition 3 . 3

 33 For any m > N , the eigenvalues u 1,m , . . . , u N,m of U m N are independent and uniformly distributed along the unit circle Proposition 3.4 Let m ≤ N and θ ∈ [0, π] and denote N m,N (θ) the number of eigenvalues of U m N inside the arc-interval 1, e iθ . Then we have ∀ x > 0:

Fig. 3 :

 3 Fig. 3: Plot of the curves t → 2 t

Fig. 4 :

 4 Fig. 4: Plot of t → 1 N 2 ln P N,strong (t) for ≤ t ≤ 2 -with = 1 5 for N ∈ {1, 2, . . . , 7} computed exactly with the formula (A.3). The black curve is the expected value ln sin π 2t

2t and θ 1 =

 1 2πk-π 2t . In the present situation we should take θ 0 = 2πk-π 2t and θ 1 = 2π -2πk-π 2t for which we find sin π 2 -2πk-π 2t = cos 2πk-π 2t

Fig. 6 :

 6 Fig. 6: Plot of t → 1 N 2 ln P N,strong, = 1 5

Fig. 7 :

 7 Fig. 7: Histograms of N τ N,strong 4 -(N -1) for N = 6 and ∈ {0.15, 0.2, 0.25, 0.3} from top-left to right-bottom obtained for n = 10 3 independent samples. The black curve is the normalized exponential distribution E(1). Empirical estimation of parameter λ ranges from 1.021 ( = 0.3) to 1.002 ( = 0.20).

Fig. 8 :

 8 Fig. 8: Histograms of N τN,strong 4 -(N -1) for N = 6 and ∈ {0.15, 0.2, 0.25, 0.3} from top-left to right-bottom obtained for n = 10 3 independent samples. The black curve is the normalized exponential distribution E(1). Empirical estimation of parameter λ ranges from 1.1 ( = 0.3) to 1.049 ( = 0.15).

9 )

 9 Note that the normalizing factor is now -N and no longer N 4 -(N -1)

Fig. 9 :Fig. 10 :

 910 Fig. 9: Histograms of n 0 (N, ) -N for N = 6 and ∈ {0.15, 0.2, 0.25, 0.3} from top-left to right-bottom obtained for n = 10 3 independent samples. The black curve is the normalized

  from (C.18) and (C.15) shows that:

  around z = z -: The computation are very similar to the z = z + case. The local parameter is z -(z) = 1

1 -

 1 cos θ 1 -θ 0 2 (C.39) C.3 Computation of F [0]

+ e i θ 0 +θ 1 2 2 2ix+ e i θ 0 +θ 1 20 1 + e i θ 1 -θ 0 2 2ix0 1 + e -i θ 1 -θ 0 2 2ix 2 √ 2 √ 2 0 1 +

 12112122221 2ix (i)x(i)(z -i) + O(1) = 0 1 + e -i θ 1 -θ 0 2ix (-i)x(-i)(z + i) + O(1) = -(i)(z -i) 2 + O(1) dy(z) = 0 1 + e i θ 1 -θ 0 2 2ix (-i)(z + i) 2 + O(1) (C.44)Eventually we can compute:j dy(a j ) dz j (a j ) = -0 1 + e -i θ 1 -θ 0 2 2x(i)i x (i) 0 1 + e i θ 1 -θ 0 2 2x(-i)i x (-i) = e -i θ 1 -θ 0 2

-1 2 (D. 5 )

 25 with M (x) > 0 on S. then proposition 5.1 is valid.

2 A 2 ln sin 2 c θ 1 6 ) 2 gives: R 2 dµ

 221622 -θ 2 2 + 2 ln |θ 1 -θ 2 | dµ(θ 1 )dµ(θ 2 ) = -1 2 A 2 ln 4 sin 2 θ 1 -θ 2 2 dµ(θ 1 )dµ(θ 2 ) = -ln 2 -Taking the Fourier transform L(s) of the function L(z) = ln sin z (x)dµ(y)L(x -y) = R L(s)|μ(s)| 2 (D.7)

  5.[START_REF] Szegö | On certain Hermitian forms associated with the Fourier series of a positive function[END_REF] The one cut case: Exact results for t ≤ 2 -

	Let us start with the following observation:
	Theorem 5.2
	∀ 0 ≤ t ≤ , t is a strong return time and Z N,t = 1	(5.20)

Acknowledgements

The author would like to thank Université Lyon 1 and particularly Université Jean Monnet and Institut Camille Jordan for the opportunity to make this research possible. O. Marchal would also like to thank his family and friends for moral support during the preparation of this article as well as B. Eynard for fruitful discussions and R. Balian for references and explanations about the underlying quantum measurement physics problem. The author also thanks G. Borot for useful discussions about the existence of a unique minimizer and an unknown referee for suggesting the use of proposition 5.1.

dθ k e i(σ(k)-τ (k))θ k = (2π) N N ! (A.1) Indeed, the integrals are only non zero when σ = τ for which the result is trivial. In particular we get by change of variables u i = e iθ i that:

A.2 Case for a single interval

We can apply the same method for the following probability: P N (b) = Prob(All eigenvalues are in the arc segment [-b, b]) where b is a fixed real angle in [0, π]. We find:

where the previous expression of Z N was used. Moreover the N ×N symmetric matrix C[b, N ] is given by:

In particular, the last formula is useful when dealing with a unique interval that is to say for t ∈ [ , 2 -]. Indeed, in that case we only need to take b = π t in order to recover P N,strong (t).

A.3 Several intervals case and exact computation at time t When the domain of integration is on the unit circle of the form θ ∈ R r=-R [a r , b r ] with a -r = -b r , b -r = -a r for r > 0 and a 0 = -b 0 we can apply a similar computation. We find:

where the matrix

In particular the previous result can be applied to compute the strong returning time probability. There are two cases depending on whether t ∈ [2k -, 2k + ] or not. Indeed, in this case, one of the interval is incomplete and one has to be careful with the determination of the angles used. We have:

In the first case we find a r = 2πr-π t and b r = 2πr+π t so that:

so that we get:

t (A.9) which can be rewritten with a determinant:

Note here that in the last formula the case when j = i is special since the sine product have to be replaced π (2R+1) t which corresponds to its Taylor expansion and is coherent with the integral formulation.

The second case is a little different since the two extremal intervals are limited by -π and π. Observe that in the sum this two terms simplify so that we find:

so we get:

12) Here one has to be careful with the case σ(k) = k for which the last formula does not make direct sense. In the integral formulation, the integral is equal to |I| = 2π 2R t + 2π(1 -2R t ) whereas the Taylor expansion of the sine functions does not contain the factor 2π. This is why we have inserted the factor 2πδ σ(k)-k in order to recover the right factor when taking the usual Taylor expansion. The determinant representation is:

In both cases, we have been able to rewrite an exact formula to compute the strong return probability as the determinant of a symmetric N × N matrix (4.7).

= -

Since the curve is of genus 0, the Bergmann tau-function is trivial and equals τ Bx ∝ (b -a)

In the end we find:

-

Computations of F [2] (noted -F (2) in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]) and higher orders fall into the general case of the topological recursion presented in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]. We find for the curve (C.15):

In the case θ 0 = -θ 1 and 0 = 1 it reduces to:

-θ 1 ,θ 1 ,1 = -3 cos (θ 1 ) -1 64 (1 + cos θ 1 ) (C.49) D Verifying conditions for proposition 5.1

In this section, we provide details on the formalism developed by Borot, Guionnet and Kozlowski and prove that it can be applied to our situation. We want to compute the large N asymptotic of the following integral:

where

A h is a union of segments described in (5.2). In our case, the intervals (as well as how many we have) depend on a time parameter t that we will consider fixed in this appendix so we omit it everywhere to simplify notations. In the formalism of [START_REF] Borot | Large-N asymptotic expansion for mean field models with Coulomb gas interaction[END_REF], our situation corresponds to pairwise interactions, i.e. r = 2 that can be rewritten like:

From them we can deduce:

and then:

F Information about numerical simulations

Numerical simulations have been carried out on a standard laptop using Maple 11 (or more recent versions) software and can easily be reproduced. Perhaps the most difficult step in all the simulations is to sample a random unitary matrix properly according to the Haar measure. We took procedures as described in [START_REF] Mezzadri | How to generate random matrices from the classical compact groups[END_REF] to sample them. Note also that exact computations for the various probabilities can directly be computed from Toeplitz determinants and turn out to be much more efficient than using matrix integrals. We can provide on demand our code to any interested reader.