Exact equivalences and phase discrepancies between random matrix ensembles

We study two types of random matrix ensembles that emerge when considering the same probability measure on partitions. One is the Meixner ensemble with a hard wall and the other are two families of unitary matrix models, with weight functions that can be interpreted as characteristic polynomial insertions. We show that the models, while having the same exact evaluation for fixed values of the parameter, may present a different phase structure. We find phase transitions of the second and third order, depending on the model. Other relationships, via direct mapping, between the unitary matrix models and continuous random matrix ensembles on the real line, of Cauchy-Romanovski type, are presented and studied both exactly and asymptotically. The case of orthogonal and symplectic groups is studied as well and related to Wronskians of Chebyshev polynomials, that we evaluate at large $N$.


Introduction
Random matrix theory [1] has developed enormously, especially in the last two decades, attracting attention from researchers in a multitude of different fields [2,3,4,5,6,7,8]. Indeed, one of its most exciting aspects is an inherent interdisciplinarity in the study of random matrices. Among the myriad of connections with other mathematical and physical areas, we have the important relationship with the study of random partitions, which in turn is deeply linked with problems in statistical mechanics, as we shall see.
In this work, we will see that random matrix ensembles which are related in different ways, either via direct mapping or emerging as two random matrix descriptions of the same object in the study of random partitions, have the same analytical evaluations for fixed values of the parameters yet, in some cases, the consideration of the large N limit leads to a different phase structure.
Originally, we started by noting that the analysis of probabilities given by a finite ensemble with weight of the Meixner type, for example as in [9], could, in principle, equally be studied with a random unitary matrix model, with a distribution of eigenvalues supported on the unit circle. Indeed, we could check, as we shall see, that the corresponding unitary matrix model gives the same results. Furthermore, the equivalence holds for the calculation of correlations near the edge of the eigenvalue density, in a double-scaling limit. However, when studying the more general ensemble with Meixner weight and a hard wall [10] and the corresponding generalized unitary matrix model, while there is again an equivalence for fixed values of the parameter, a different asymptotic behaviour appears, with phase transitions of second order instead of third order. This discrepancy is a consequence of the complexification of the logarithm of the weight function of the more general unitary random matrix ensembles, studied in Section 2.4 below.
We further elaborate on this whole notion by studying other equivalences, with continuous random matrix ensembles on the real line, via direct mapping this time. We first introduce notions of random partitions and how either discrete random matrix ensembles, such as the Meixner ensemble with a hard wall, or continuous unitary random matrix ensembles appear.
In the present work, we are interested in a special case of Schur measure, which is a particular sub-case of what is known in the literature as z-measure [14,15]. We choose the parameters to be (1.2) t = (t, · · · , t N 1 times , 0, . . . ), t = (t, · · · , t N 2 times , 0, . . . ), 0 < t < 1, so that the Schur measure (1.1) becomes the particular instance of z-measure with z = N 1 , z = N 2 ∈ N, with the sequence of measures on partitions of n is assembled into the z-measure on all partitions (any n) through a negative binomial distribution with parameter t 2 . This choice of z, z is the degenerate case, i.e., non-negative measure, and positive definite when restricted to partitions of length at most N 1 . The probabilities where by Prob t we understand the probability taken with specialization of parameters (1.2), admit a determinantal representation in terms of the so-called hypergeometric kernel [14]. The properties of the hypergeometric function were used in [14] to show that, when z = N 1 , z = N 2 are integers, the hypergeometric kernel becomes proportional to the Meixner kernel. For a collection of results about the z-measure, see the survey [16], and in particular Proposition 6.1 therein for the connection with the Meixner ensemble.
In the present work we are interested in the quantities and (the sum is over all λ with λ 1 ≤ K) with choice of parameter understood to be as in (1.2), and |λ| = j λ j is the size of the partition λ, corresponding to the total number of boxes in its diagram, and dim λ is the dimension of the irreducible representation of the symmetric group S |λ| labelled by λ. The meaning of the subscripts H and E will be clear in a moment. Notice that, due to the property s λ (t) = 0 if length(λ) > length(t), the sums are effectively truncated to partitions of length at most N 1 . Furthermore, the ratio Z E (t) Z H (t) = Prob t (λ 1 ≤ K) is the probability that, picking a random partition λ with probability distribution as described above, the length of its rows is at most K. The z-measure induces a determinantal point process on Z + 1 2 , thus the correlation functions have determinantal form Prob t (X ⊂ S(λ)) = det (K) X , for X ⊂ Z+ 1 2 , where K is the operator whose kernel, the function K(x, y) on Z + 1 2 × Z + 1 2 , is the hypergeometric kernel. Therefore we have a Fredholm determinant representation of Z E (t): where 1 in the last line is the identity operator. Using the hook-length formula dim λ = j<k λ j − λ k − j + k k − j for the dimension in eq. (1.4) and changing variables h j = λ j −j +N , we arrive to the expression where we also used the symmetry in the h j variables to remove the restriction to the Weyl chamber λ 1 ≥ λ 2 ≥ · · · ≥ λ N 1 . G(·) is the Barnes G-function [17], which, for integer values of the argument is G(n) = n−2 j=0 j!. The expression (1.5) is a Meixner ensemble with summation restricted to 0 ≤ h j ≤ N 1 + K − 1 and coincides with the partition function of the dimer model studied in [9,10] 1 , and has also appeared in [19]. In the Coulomb gas picture, the restriction in the summation range corresponds to a hard wall for the charges (i.e. an infinite barrier) placed at N 1 + K − 1. We assumed N 1 ≤ N 2 , but a completely analogous expression can be easily obtained in the converse case.
Let us introduce the generating functions of, respectively, the complete homogeneous polynomials {h k } and the elementary symmetric polynomials {e k }, specialized at t: The partition functions (1.3)-(1.4) admit a determinantal representation in terms of Toeplitz determinants (see for example [20,21]), with symbol, respectively These Toeplitz determinants in turn admit a representation as unitary matrix integrals We will refer to (1.6) and (1.7) as the H-model and the E-model, respectively. Therefore we have two equivalent matrix model descriptions of the quantity Z E (t), or, equivalently, of the probability Prob t (λ 1 ≤ K): either as a discrete matrix model on the bounded subset {0, 1, . . . , N 1 + K − 1} ⊂ Z, or as a continuous matrix model on the unit circle. 1 The dimer model is associated with the random tiling of an Aztec diamond with a square [9] or a rectangle [10] cut off. See e.g. the recent work [18] for an overview on dimers, random tilings and random matrices. Table 1. Dictionary between the notation in [10], in [22] and the present work.

Random matrix ensembles on the unit circle
In the present Section, we consider the asymptotic behaviour of the unitary matrix models defined in eq. (1.6) and (1.7), when the rank K is large and N 1 , N 2 scale with K. Before that, we comment on the already known aspects of the exact solvability of some of the models above.
2.1. Prologue: Exact evaluation. Following [22], the authors of [10] gave an explicit evaluation of the discrete matrix model (1.5) at the limit value t = 1. The exact formula of [10, Proposition 3.1] was obtained thanks to the fact that at the limit value t = 1 the Meixner ensemble with a hard wall becomes a Hahn ensemble.
On the unitary matrix model side, at t = 1 (or more generally |t| = 1) the symbol σ E develops a Fisher-Hartwig singularity, and it can be evaluated exactly thanks to a formula by Böttcher and Silbermann [23] (see also [24] for another proof): This provides an exact check of the equivalence. We will provide a third independent derivation of this result later in Section 3.1.1. For ease of the reader, we report in table 1 the correspondence between the notation of [9,10], that in [22] and ours. Note that, if we choose the symmetric model N 1 = N 2 ≡ N and modify the symbol by inserting a monomial factor z s , with s ∈ Z, the exact evaluation of the corresponding matrix model at t = 1 is again given by the formula of [23,24], thanks to the simple identity when z = e iϕ , −π ≤ ϕ ≤ π. The effect of the monomial insertion is to shift the Fourier coefficients of the symbol σ E (z; 1) by s: the k th Fourier coefficient of the symbol z s σ E (z; 1) is the (k + s) th coefficient of σ E (z; 1). Therefore, we can evaluate exactly the partition function of the unitary ensemble with weight (2.2), and it is again given by formula (2.1), with It is worth mentioning that there are a wealth of analytical results on Toeplitz banded matrices, whose determinant is given by the E-models above [25]. More generally, there are many analytical results on determinants of Toeplitz matrices generated by a rational symbol, but we will explore complementary approaches here and indicate possible open problems at the end, in the Outlook Section.
We just quote here that for (1.7) with N 1 = N 2 = 1 we have the determinant of a symmetric tridiagonal Toeplitz matrix, known to be equal to a Chebyshev polynomial of the second type:

2.2.
Unitary matrix models. It is well known that, although no Cauchy-type identity exists for the E-model, its K → ∞ limit, with fixed N 1 , N 2 coincides with the H-model: In the scaled limit with N 1 , N 2 growing together with K, it was found in [9,10] that the presence of hard walls in the discrete matrix model (1.5) triggers a third order phase transition. For the case N 1 = N 2 , starting from a result by Baik [26] we will prove the phase transition from the point of view of the unitary matrix model. This is done in Section 2.3. For the general case N 1 = N 2 , however, the potential of the unitary matrix model is complex-valued, and the large K asymptotic becomes more involved. This topic is analyzed in Section 2.4. As the asymptotic behaviour at large K does not depend on N 1 , N 2 being integers, we consider the matrix models arising from Toeplitz determinants with more general symbols and, without loss of generality, we assume 0 ≤ β 1 ≤ β 2 . We also introduce the notation where β = (β 1 + β 2 )/2 is the average power and 0 ≤ v ≤ 1 measures the asymmetry in z ↔ z −1 . We then take the K → ∞ limit with 2 (2.4) γ := β K fixed.
As customary, when studying the large rank behaviour of matrix models, we introduce the density of eigenvalues which at large K becomes a continuous function of ϕ, with compact support and normalized so that π −π dϕρ(ϕ) = 1.
In each of the cases considered below, we will find the eigenvalue density ρ and use it to evaluate the free energy in this limit, defined as: Our calculations are based on standard saddle point techniques, and we omit them from the main text and refer to the Appendices A and B. We will show that all the models undergo a phase transition when a gap opens in the support of the eigenvalue density, as schematized in Figure 1. 2.3. Phase transition: symmetric case. We first focus on the symmetric case β 1 = β 2 ≡ β, while the analysis of the more general case β 1 = β 2 is undertaken later in Section 2.4.
Hence, the matrix models we analyze are: and we recall that Z sym.
u,H admits an exact solution through the Cauchy identity, whilst Z sym.
u,E does not. Nevertheless we have lim K→∞ Z sym. u,E = Z sym. u,H . In [26], Baik proved that the first system, described by the partition function (2.6), and which we will call for simplicity the H-model, undergoes a phase transition at large K. We prove that the second system (2.7), which we call E-model, undergoes the same phase transition. We prove it solving a singular integral equation in Appendix A, but in fact the result may also be directly obtained from [26], with minor changes.
2.3.1. The H-model. Consider the matrix integral in (2.6), and take the limit K → ∞ with β scaling as in (2.4). The leading contribution to Z sym.
u,H comes from the solution to the system of saddle point equations that, with the help of the eigenvalue density ρ H as defined in (2.5), can be rewritten as a single singular integral equation: where the symbol P means principal value of the integral, and ρ H is the eigenvalue density for the specific model considered presently. The details of the solution to eq. (2.8) are spelled in Appendices A.1 and A.2. Two phases exists, separated by the critical curve [26] γ = 1 + t 2t =: γ c,H (t).
The eigenvalue density, plotted in Figure 2 for various t and γ, reads: and allows to evaluate the free energy F sym. u,H , obtaining: , and therefore the system undergoes a third order phase transition at the critical curve γ c,H (t) = 1+t 2t .

2.3.2.
The E-model. We now turn to the second matrix model, defined in eq. (2.7). The leading contribution in the large K limit, with scaling (2.4), is obtained solving the saddle point equation We solve this singular integral equation in Appendices A.3 and A.4. From direct comparison of the integral representation of Z sym. u,H and Z sym. u,E in (2.6) and (2.7), one would expect that, the solution to the second model is related to the solution to the first model by Direct calculations prove that this is true and, in particular, the system undergoes a phase transition along the critical curve The eigenvalue density in the E-model, plotted in Figure 3 for different values of t and γ, is: which allows to compute the free energy F sym. u,E at large K (see Appendices B.1 and B.2), and the final result is: where C E (γ) is t-independent. Taking derivatives, one finds again that first and second derivatives are continuous functions, while thus the phase transition is of third order also in this case. . Eigenvalue density ρ E (φ). The blue curve is at γ = 1 2 γ c,E (t) and the red curve is at γ = 2γ c,E (t), for t = 0.1 (left), 0.5 (center), 0.9 (right).
2.3.3. Double-scaling limit. We have shown that the free energy, and consequently the phase structure and the critical curve, of the matrix model with symmetric weight of E type at large K and large N agrees with the equivalent descriptions as a Meixner ensemble with a hard wall studied in [9]. In fact, we find an even stronger agreement between the two pictures, as also the correlations among eigenvalues near the edge of the distributions in a double-scaling limit match. For the unitary ensemble, the double-scaling limit is described in [26], along the lines of [27]. On the other hand, in the discrete ensemble, if the critical region is approached from the small coupling phase, the hard wall is not active and the results of [22] directly apply. The doublescaling is the same in both cases, and leads to the Tracy-Widom law [28]. A third consistency check comes from the asymptotics of the hypergeometric kernel: standard computations using a steepest descend method (see for example [13,Sec. 4] for an introduction) show that the result in the double-scaling limit matches the Tracy-Widom behaviour.

2.4.
Phase transition: general case. While in the previous Section 2. 3 we focused on the analysis of two matrix models that are symmetric under z ↔ z −1 , we now allow the more general situation β 1 = β 2 , in which the symmetry is lost.
The matrix models we consider here are The weight functions are complex valued, thus we expect the eigenvalue densities at large K to be complex-valued functions. We rewrite (β 1 , β 2 ) in terms of the parameters (β, v) defined in (2.3), and consider the limit K → ∞ with scaling of β as introduced in (2.4).

2.4.1.
The H-model. We first focus on the behaviour of the partition function (2.13) in the limit described above, and show how the discussion of the symmetric case in Section 2.3 is modified. Details of the calculations can be retrieved in Appendices A.5 and A.6 for weak and strong coupling respectively.
The leading contributions at large K may be encoded in the eigenvalue density ρ H which solves the integral equation The asymmetry parameter v complexifies the left hand side of the latter equation, and the resulting eigenvalue density is complex (see Appendices A.5 and A.6): with critical value γ c,H = 1+t 2t , the same as in the symmetric case. The corresponding free energy F gen.
u,H = −K −2 log Z gen. u,H , computed in Appendix B.3, can be written as: . The introduction of the asymmetry reduces the order of the phase transition from third to second, with: The E-model. Consider now the second matrix model, the E-model of eq. (2.14), and take its scaled limit as in (2.3)-(2.4). The saddle point equation reads The eigenvalue density is complexified by the presence of the asymmetry parameter v (see Appendices A.7 and A.8): with the same critical value as for the symmetric case: γ c,E = 1−t 2t . The free energy is (see Appendix B.3): The firs derivative of the free energy is continuous, but the second derivative is not: thus the phase transition is second order. We plot the first derivative of free energy as a function of γ, at different values of t, in Figure 4.
2.5. The Gross-Witten limit. Consider the two matrix models of Section 2.3, defined in equations (2.6) and (2.7). The potential of those models can be written as: with + sign for the E and − for the H. We send t → 0 and β → ∞, keeping their product β GW := tβ fixed. This gives: which is the potential of the Gross-Witten matrix model [29,30]. Note that the limit is the same for both the E-and the H-model. In the more general, non-symmetric case, with potential we pass from (β 1 , β 2 ) to (β, v) as prescribed in (2.3), and define β GW := tβ. The same limit as above gives: which is the potential of the Gross-Witten model with a topological θ-term [31,32]. Since, by construction, we are in the regime |v| < 1, the Gross-Witten phase transition is still present, as discussed in [32].
Within the setting laid down in the Introduction, the limit above with β = N ∈ N of the normalization of the z-measure leads to the normalization of the Poissonized Plancherel measure on partitions [33,16]. This is consistent with what we mentioned above, as the Poissonized Plancherel measure is mapped to the Gross-Witten unitary matrix model [34]. On the other side, the Poissonized Plancherel measure is in the t → 1 limit of the Meixner ensemble [35] as well. The dashed back line is the symmetric case v = 0. The point at which the curve is continuous but with discontinuous derivative becomes more and more visible as v is increased.
2.6. Phase transition: from third to second order. For 0 < |v| < 1, the potential of the unitary matrix models (2.13)-(2.14) becomes complex-valued. We have studied the large K limit and showed that the model undergoes a second order phase transition, which becomes third order turning off the asymmetry parameter, v → 0. Nevertheless, the phase transition in the Gross-Witten matrix model remains third order if the potential is modified into cos ϕ − iv sin ϕ.
An alternative approach to face the matrix models with complexified potential is to analytically continue them in the sense of [36]. We relax the condition |z| 2 = 1 and deform the integration cycle S 1 C v in the complex plane so that the potential remains real along C v (see Figure 5). Of course, for the deformation to be smooth, C v should be a smooth Jordan curve, whose shape depends on v and which is the unit circle S 1 when v = 0.
As a warm up, we apply the analytic continuation to the Gross-Witten model: we seek a contour C v such that . We deform the integration contour from the unit circle to a Jordan curve along which the potential is real.
for all z ∈ C v . We find that C v is simply a circle of radius Hence, not only it is real-valued, but we also recover the original, symmetric Gross-Witten model sitting on a rescaled circle, see Figure 6. Furthermore, C v is sent to infinity or shrinks to a point as v → 1 or v → −1 respectively, and we cannot prolong beyond these values. This provides new perspective on the result of [31,32], where a drastic difference in the behaviour was observed crossing from 0 < |v| < 1 to |v| > 1.
For a generic polynomial potential the smooth integration cycle C v is determined as the locus in C that solves an algebraic equation.
For the models with weights of E and H type with 0 < t < 1, however, for every real v = 0 the unique contour such that the potential is real-valued is a half real line together with an open segment (the details depend on whether we consider the E-or the H-model), which is not a smooth deformation of S 1 .
2.6.1. Phase transition and universality. To sum up our conclusions in one sentence, the unitary matrix models with symmetric E and H weight undergo a third order GW phase transition, but their more general, non-symmetric extensions have complex potentials and therefore fall out of the GW universality class. This ought to be contrasted with what happens in the equivalent description as a discrete ensemble with a hard wall. In that case, increasing the parameter γ from zero, the system undergoes a third order phase transitions at the value γ = γ c,E when the hard wall becomes active. This holds both in the symmetric and the more general setting: they belong to the same universality class. See [37,38] for detailed discussion on the universality of the third order phase transition in presence of a hard wall. We also stress that the discrete topology further constrains the eigenvalue density, but the condition is satisfied for all values of the parameters [9,10]. Thus, the discrete nature of the ensemble plays no role in determining the phase structure of this model.

Random matrix ensembles on the real line
We now discuss two random matrix models on the real line that follow from a change of variables of the unitary matrix models above described. One family of random matrix ensembles is obtained through the 1d stereographic projection, leading to a Cauchy-Romanovski type of ensembles and the other one is based on the mapping x = cos θ which is more useful when the original matrix model has a symplectic or orthogonal symmetry (that is, corresponding with Toeplitz±Hankel determinants, instead of Toeplitz determinants), but that also lead to explicit expressions for the symmetric unitary matrix models, in terms of Wronskians of Chebyshev polynomials of the four types.
3.1. Cauchy ensembles and Romanovski orthogonal polynomials. We analyse the matrix models Z sym.
u,H and Z sym. u,E defined in eq. (2.6) and (2.7) and consider a change of topology, in which we remove a single point from the unit circle, to obtain S 1 \ {∞} ∼ = R. We expect the asymptotic behavior of these models at large K to correspond to the "gapped" phase of the unitary matrix model (see Figure 1), which appears at strong coupling. We pass from angular variables to the real line using the stereographic projection [2, Sec. 2.5]: as sketched in Figure 7. The Vandermonde determinant is mapped to: and, with the corresponding transformation of the symbols, the resulting matrix models are: Likewise, we can also use the stereographic projection in the more general, non-symmetric, matrix models (2.13) and (2.14). We obtain: where we used the redefinition of parameters β 1 = β(1 − v) and β 2 = β(1 + v).
3.1.1. Exact evaluation. In the limit case t = 1, the partition function of the symmetric E-model, after stereographic projection, takes the simple form: This random matrix ensemble has been studied as a Cauchy ensemble [39], Lorentz ensemble [40], and in this form it is a particular case of the classical ensemble with weight [2,41] σ This weight function satisfies the Pearson equation and hence the associated random matrix ensemble is classical [41,42,2], although in many references the listing of classical ensembles appears restricted to Hermite, Laguerre and Jacobi. See e.g. [42] for the expanded list of possible classical weights. The associated polynomials go under many names, including pseudo-Jacobi [43,44,45,46,47] due to their (non-trivial) relationship with Jacobi polynomials, and also appear in [48]. A proper name seems Romanovski polynomials R (α 1 ,α 2 ) n given [49,43], see [50] for a review. We evaluate now the matrix integral (3.5) using Romanovski polynomials. They satisfy (the dependence on (α 1 , α 2 ) is understood): where we have chosen a normalization such that the polynomials R n are monic, and their norm squared is [51] where the second square bracket is the change of normalization to obtain monic Romanovski polynomials from the normalization of [51]. Therefore, using with h n specialized to the case of interest α 1 = α 2 = K + β, we obtain where the Barnes G-function [17] is identified using The result is the same indeed as eq. (2.1). Therefore, the Z E (t = 1) is computed in three different ways: as a discrete ensemble on a finite set, as a unitary ensemble, and as the Cauchy ensemble on the real line. The tools used in each of the three approaches are, respectively: the Hahn polynomials [10], the Toeplitz determinant with symbol with a pure Fisher-Hartwig singularity [23], and the Romanovski polynomials, as represented in Figure 8.
Unit circle

Asymmetric Toeplitz matrix.
A much studied symbol in spectral analysis of Toeplitz matrices [52,53] is for integer s ∈ Z. As discussed in Section 2.1, this has the effect of shifting the Fourier coefficients of the symbol by s, thus it shifts the diagonals of the Toeplitz matrix upward (if s > 0) or downward (s < 0) by s and the Toeplitz banded matrix associated to our model becomes asymmetric, making it more similar to the generalized model in terms of this asymmetry of the associated Toeplitz matrix. The stereographic projection of this extra term is: For t = 1, the resulting weight function is This is the stereographic projection of the weight introduced in eq. (2.2), and it is still of the Romanovski form (see [41]) with the identification We can therefore give explicit evaluation of the determinant of the Toeplitz banded matrix with symbol (3.7) using again Romanovksi polynomials: Note that, if we consider the random matrix ensemble with weight (3.7) and keep the normalization such that it coincides with the stereographic projection of the unitary ensemble with weight (2.2), discussed in Section 2.1, the factor in bracket in (3.8) cancels exactly.
3.2. Large K limit. We now focus on Z sym. E,stereo. , in eq. (3.2), and study its large K limit, with γ = β/K fixed. The saddle point equation is: The parity symmetry of the matrix model guarantees that we can look for a symmetric solution with suppρ = [−A, A], A > 0. We expect that A will be back-projected to e iφ 0 of Section 2.3 undoing the stereographic projection.
We report the details of the computations in Appendix C.1. We arrive to the eigenvalue density: , with the boundary A fixed by normalization: This solution is well defined as long as and sending A → ∞, which would be back-projected to φ 0 → π, corresponds to the limit γ ↓ γ c,E = 1−t 2t . This matches the analysis on the circle. We can do the same with the more general matrix model Z gen.
E,stereo. of eq. (3.4). In this case, as happened on the circle, turning on the asymmetry parameter v complexifies the eigenvalue density, which becomes: , with same value of A as above. The derivation of the result is given in Appendix C.2. For completeness, we also report the details of the large K analysis of the matrix model with weight (3.6), corresponding to the determinant of an asymmetric Toeplitz matrix, in Appendix C.3.  π −π σ(e iθ )e −ijθ dθ, with 3 σ(e iθ ) = σ(e −iθ ). Then,

Orthogonal and symplectic symmetries and
We give a new proof of this result following [55]. First, we set the notation for the rest of the Section. Denote O ± (N ) the real orthogonal matrices with determinant ±1, and Note that the Haar measures on Sp(2K) and O − (2K + 2) coincide. Let us consider matrix integrals with weight σ(e iθ ) over the group G(K) (see [2,Sec. 2.6] for details on the Haar measure on orthogonal and symplectic groups), which we denote by Z G(K) σ . Proof of Lemma 3.1. The change of variable x j = cos θ j gives (3.10) with a, b ∈ {±1/2} depending on the symmetry of the original model, and the overall normalization depends on the specific G(K): Using Adreiéf identity, expression (3.10) gives the four Hankel determinants on the right hand sides in the Lemma, according to the choice of G(K). The Toeplitz±Hankel determinants on the left hand sides in the Lemma correspond with matrix integration over orthogonal and symplectic group [56,55]. Indeed, for σ(e iθ ) = σ(e −iθ ) Furthermore, as noted in [55], This concludes the proof of the four identities.
For the case of the symmetric E and H symbols considered in the paper, we have expressions of the general form: with ξ = 1 + t 2 /2t and a, b as in (3.11). The ± signs are always the upper sign for the E-model and the lower for the H-model. Moreover, we take β = N ∈ Z >0 a positive integer. Thus, Z G(K) E/H are averages of moments of characteristic polynomials, taken in a random matrix ensemble where the weight function is of the Jacobi type, ω (x) = (1 + x) a (1 − x) b , but always with the specific values of a and b above, corresponding to the Chebyshev subfamily. Hence, the corresponding orthogonal polynomials for are Chebyshev polynomials of first, second, third and fourth kind respectively, denoted by Their correspondence with group integration is (3.11). As mentioned, instead of interpreting the resulting model (3.12) as the partition function of a model with a semi-classical weight function is an interesting possibility and may lead to connections with Painlevé equations and other integrable models), we adopt the characteristic polynomial point of view [57,58], which, for the + sign case in (3.12) (below we discuss the other choice), gives that where denotes the Wronskian determinant and P n (·) denotes one of the four Chebyshev polynomials. We also have defined the coefficients .
(C H is for later use). The result follows immediately from Lemma 3.1 above (i.e. [54, Lemma 2.7]) together with the known characteristic polynomial results [57]. It has also been the specific subject of a later paper [59], where the same result is obtained, for a polynomial and, as above, even, symbol. Therefore it holds that: In addition, we can also study the case with symbol H in the same manner. For negative moments the characteristic polynomial description becomes involved, having to consider not the polynomials orthogonal w.r.t. the given weight, but the Cauchy transform of these polynomials [58]. In some cases this is a considerable complication, however, for Chebyshev polynomials, it is well known that the polynomials of first and second kind (and likewise the polynomials of third and fourth kind) are integral transforms of each other with respect to weighted Hilbert kernels [60]: with P the Cauchy principal value integral. The weighting is precisely the corresponding weight function, which is how it is required [58]. Note also that ξ > 1 for 0 < t < 1, thus ±ξ falls outside the integration region: in this case the Hilbert and Cauchy transform simply differ by a factor √ −1/2. Therefore, when 0 ≤ N < K we have: Thus, we obtain a duality between the E model with O − (2K + 1) and the H model with O + (2K + 1) and conversely. The parameters are mapped K → K − N , ξ → −ξ under the duality. We also obtain a tight relation between the E model with O + (2K) and the H model with Sp(2K) and conversely, although it is not really a duality in the sense that we cannot recover one result from the other via a bijective map of the parameters. So far, we have discussed the same models as above but with a symplectic or orthogonal symmetry (it would be interesting if such model had an interpretation in terms of a dimer problem, perhaps with boundaries, but we do not know of such an interpretation). The unitary case follows from the identity [61]: , which is used in matrix integral form in [55] and was used in [62] to precisely extend the Wronskian result [59] (again for a symbol of the E type) to the unitary case. Therefore, from this Wronskian point of view, the unitary matrix model is more complicated than the orthogonal and symplectic ones, because it is the product of two Wronskians. The formulas for the symmetric E symbol are 4 : and for the H symbol We obtain a duality for U (2K) with E and H symbol, with parameters mapped according to ξ → −ξ and K → K − N together with the factor (−1) N in the H-model. The relation between the model with E and H symbol and group U (2K + 1) is not a duality. It is worth mentioning that Wronskians of orthogonal polynomials in general satisfy an interesting identity which expresses them as Hankel determinants, but with the roles of the K and N swapped [63]. In the case of the Chebyshev Wronskian of the first type, the Hankel determinant is in terms of Legendre polynomials and for the second type one is in terms of a Hankel of Gegenbauer polynomials C > 0 and f (x) Lipschitz continuous, has been obtained by Basor and Chen [64] when a ≥ 0, b ≥ 0. They also showed that the resulting expression can be analytically continued to all a ≥ − 1 2 , b ≥ − 1 2 , and therefore it applies to the matrix ensemble (3.12) for each of the four symmetry groups. In that case, a, b ∈ + 1 2 , − 1 2 and f (x) = (ξ ± x) ±β , with upper sign for the E model and lower sign for the H model. Recall also that ξ = (1 + t 2 )/(2t) > 1, hence f (x) is admissible in both models. The formula of [64] for generic a, b and f is: for large K, where Z(f ) is independent of K and is given by: This expression directly gives the asymptotic behaviour of Z G(K) E/H for both symbols. First, we notice that the numerical prefactor in (3.12) simplifies with the one coming from formula (3.13), exactly for G(K) ∈ {Sp(2K), O ± (2K + 1)} and leaving behind a factor of 2 when G(K) = O + (2K). Then, we use and see that the contribution from the exponential in (3.13) cancels: We are left with the contribution from Z(f ). The ratio of products of Barnes G-functions and the Gamma function is trivial for G(K) ∈ {Sp(2K), O ± (2K + 1)}, and is 1 2 for G(K) = O + (2K), thus cancelling the left over factor from the numerical coefficient. The exponential of the single integral gives and in particular is trivial for odd rank symmetry, G(K) = O ± (2K + 1). The double integral which appears in the exponential in Z(f ) is universal for the four symmetries and takes the form The last integral admits an exact solution in terms of hypergeometric functions, albeit involved. Instead, we write ξ explicitly as (1 + t 2 )/2t and, using x = cos θ, we expand the logarithm: After integration the n = 0 term cancels, while for n ≥ 1 we use the Hilbert transform property of the Chebyshev polynomials to arrive at the final expressions with upper sign for the E-model and lower sign for the H-model. Note that we cannot recast the series into a logarithm because ξ > 1.

Outlook
It would be interesting if semiclassical polynomials can be used to obtain further analytical results on the model with orthogonal and symplectic symmetries, using the Wronskian representation presented here. Notice that computing these Wronskians is equivalent to computing determinants of banded matrices of the Toeplitz±Hankel type, which arise for example as Laplacian or difference matrices of problems with boundaries [65,66].
Other analytical approaches could have been followed as well. For example, it is immediate to interpret (3.2) as a characteristic polynomial average over a Cauchy ensemble 5 . It would be more interesting to use this approach for the case of the asymmetric symbol (3.6), which, when mapped on the real line, gives the full Cauchy-Romanovski weight. This could be studied in the context of asymmetrically banded Toeplitz matrices.
At large N , it could be interesting to consider the recent results [67] on correlators of characteristic polynomials for classical ensembles. If a confluent limit of these results on correlators can be taken, that would directly evaluate the Wronskian description of the orthogonal and symplectic ensembles discussed in Section 3.3. On the other hand, the asymptotic formulas in [67] already could be directly applied to a q-deformation of the E-model, such as the one studied in [21,Sec. 4.3].
It is worth mentioning that the Cauchy ensemble has interesting connections with the Laguerre ensemble and Painlevé equations [68,69] that may be further exploited in our context. We could also focus further on the exact equivalences part and see if other averages over the Meixner ensemble, such as the ones in [70], have a counterpart in terms of the E and H unitary random matrix ensembles.
With regards to the phase structure, it seems also worth exploring if the phase transitions obtained, which are phase transitions of the determinant in the Toeplitz matrix representation of the E and H models, are related to any spectral change of the Toeplitz matrices, which in the case of the E-model is a banded matrix. The spectra of such matrices is still of very much interest and a rich subject [52,53] and it could be a worthwhile direction to study the spectra in the natural scaling limit suggested by the random matrix representation, the one explored here.
We conclude by mentioning that passing from a third order to a second order phase transition opens up the possibility of intriguing interpretations. According to [71], the negative jump of the second derivative of the free energy when crossing the critical line is associated with the creation of a meta-stable state, and the second order phase transition corresponds to tunneling from the meta-stable vacuum to the stable one. This aspect may deserve further attention. Here we present the calculations to solve the saddle point equations obtained in the scaled large size limit of the unitary matrix models discussed in Section 2.
A.1. Solution for Z sym.
For the left hand side we use: where in the second equality we recognized the generating function of Chebyshev polynomials of second kind, U n (x). Exploiting the basic property U n (cos ϕ) = sin(n+1)ϕ sin ϕ , equation (A.1) is rewritten as: [a n sin(nϕ) − b n cos(nϕ)] , and we immediately get: (A.2) a n = γ π t n , b n = 0, for all n = 1, 2, . . . . The only yet undetermined coefficient a 0 is fixed by normalization to a 0 = 1/(2π). Putting all together: which can be further simplified recognizing the generating function of Chebyshev polynomials of first kind T n (x), using the property T n (cos(ϕ)) = cos(nϕ). We finally arrive to: The minima of this function are located at ϕ = ±π, and imposing the condition ρ(ϕ) ≥ 0 for all −π ≤ ϕ ≤ π, we see that the present solution holds in the regime Above the critical value γ c,H , solution (A.3) ceases to be valid and we must drop the assumption suppρ = [−π, π] and look for a different solution.
We introduce the following complex function, named resolvent: Here, L ⊂ S 1 is the arc of the unit circle from e −iφ 0 to e iφ 0 , oriented anti-clockwise, and the complex function ψ H (u) is the continuation of ρ H (ϕ) to the complex plane, with ψ H (e iϕ ) = ρ H (ϕ). The eigenvalue density is recovered from the resolvent through the relation with the subscript + (resp. −) meaning the limit taken from outside (resp. inside) the unit circle.
Following standard methods, we introduce two auxiliary complex functions: and a function Φ(z), to be determined, such that The saddle point equation (2.8) becomes a discontinuity equation for Φ(z). We will use this discontinuity equation, together with the fact that, by definition, Φ(z) decays as ∼ 1/z as |z| → ∞, to obtain an integral expression for Φ(z) (see [73,72]). Explicitly: where we denoted for shortness W (u) the left hand side of eq. (A.1), and the contour C encloses the branch cut L, the arc along the unit circle. See Figure 9. At this stage, we notice that the left hand side of eq. (A.1) has poles but no branch cuts, so we can manipulate the integration contour C in a convenient way. We arrive to: where the three contributions are: The sum in I 2 runs over the poles {z p } of the derivative of the potential (that is, the poles of the left hand side of eq. (A.1)), in this case z p = t ±1 , and C R is a large circle of radius R. We have I 3 (z) = 0 and I 1 (z)h(z) yields no discontinuity, so it is irrelevant for the evaluation of ψ(z). The unique relevant contribution comes from: where for the last equality we used h(t −1 ) = −t −1 h(t). In general, on the real axis h(x) > 0 if x > 1 and h(x) < 0 if x < 1. Therefore, since 0 < t < 1, we bare in mind that −h(t) > 0.
Plugging the expression for Φ(z) into Ψ(z) and taking its discontinuity along the arc L, we arrive to: The boundary φ 0 of the support is fixed by normalization: which provides (γ − 1) and hence, writing h(t) explicitly, we obtain: Also, plugging the expression for h(t) in (A.7), we arrive to the final expression A.3. Solution for Z sym. u,E at weak coupling. In this Appendix we solve the saddle point equation (2.11), which we rewrite here for completeness: The analysis will follow closely Appendix A.1 for small coupling and Appendix A.2 for strong coupling, and we omit most of the details. We again start with the ansatz suppρ = [−π, π], and manipulate the right hand side of (A.8) exactly as we did in Appendix A.1, and for the left hand side we use: Then, eq. (A.8) becomes [a n sin(nϕ) − b n cos(nϕ)] , with a n , b n the Fourier coefficients of ρ E (ϕ), and we immediately obtain: where, as usual, a 0 = 1/(2π) is fixed by normalization. We recognize the generating function of Chebyshev polynomials of the first kind, using the standard identity T n (cos(ϕ)) = cos(nϕ), and arrive to the final expression for the eigenvalue density ρ E (ϕ) = 1 2π 1 + 2γt cos ϕ + t 1 + t 2 + 2t cos ϕ .
The minima are placed at ϕ = ±π, and the non-negativity condition ρ E ≥ 0 holds as long as: We see that, as expected by naïve comparison of the integral representation of the matrix models (2.6) and (2.7), the eigenvalue densities ρ H and ρ E , as well as the respective critical points, are related through (γ, t) ↔ (−γ, −t).
A.4. Solution for Z sym. u,E at strong coupling. The procedure at strong coupling is as in Appendix A.2, and we avoid the technical details here. We assume a one-cut solution supported on [−φ 0 , φ 0 ], 0 < φ 0 < π, and introduce the resolvent with integration contour L meant to be the arc along the unit circle connecting e −iφ 0 to e iφ 0 , and the function ψ E (u) being the continuation of ρ E in C, with ψ(e iϕ ) = ρ E (ϕ). We introduce, as in Appendix A.2, the complex functions h(z) and Φ(z), and the saddle point equation (A.8) becomes a discontinuity equation for Φ(z). The only relevant contribution (for the computation of ρ E ) to Φ(z) is: where for the last equality we used h(−t −1 ) = t −1 h(−t), with both understood to be negative. We get: and sin(φ 0 /2) is fixed by the normalization: Thus, Combining this latter expression with eq. (A.9), we finally arrive to: A.5. Solution for Z gen. u,H at weak coupling. This Appendix is dedicated to the solution, in the weak coupling regime, to the saddle point equation (2.15), which we rewrite here for convenience: The procedure is similar to the one adopted in Appendix A.1, but now the parameter v makes the left hand side complex-valued. Therefore, we expand the eigenvalue density in the exponential form: while we also use 6 : The right hand side of (A.10) then becomes: For the left hand side we can expand the geometric series and obtain: Putting all together, and taking into account the normalization condition that fixes ρ 0 , we arrive to: We see that the parameter v, which controls the asymmetry, introduces an imaginary part in the eigenvalue density. Since −1 < v < 1, the relevant minima are those of the real part, thus the critical value is the same as in the symmetric case: and above this value the solution ceases to be valid.
A.6. Solution for Z gen. u,H at strong coupling. When γ > γ c,H we have to look for a different solution to the saddle point equation (A.10). At strong coupling, the procedure is exactly the same as in Appendix A.2, since we were already working in the complex plane, so the complexification of the left hand side of eq. (A.10) for v = 0 does not alter the strategy. We assume a gapped one-cut solution supported on [−φ 0 , φ 0 ] and introduce the resolvent Everything goes trough exactly as in Appendix A.2, except for the factors of (1 ± v) appearing in the left hand side of eq. (A.10). One easily finds: The normalization fixes sin (φ 0 /2) 2 , through: Notice that, due to the parity of the integral, the result is independent of v, in particular it is equal to the symmetric case (v = 0). The final expression for the eigenvalue density: with sin(φ 0 /2) the same as in [26] and Appendix A.2, that is A.7. Solution for Z gen. u,E at weak coupling. Here we solve the large K limit of the second matrix model, eq. (2.14), in the general case. The saddle point equation (2.17) is: The procedure is as in Appendix A.5. The right hand side reads: while, expanding the geometric series, the left hand side becomes: This gives: ρ E (ϕ) = 1 2π 1 + 2γt cos(ϕ) + t − iv sin(ϕ) 1 + t 2 + 2t cos(ϕ) .
Again, the validity of this solution extends up to: as in the symmetric case.
A.8. Solution for Z gen. u,E at strong coupling. When γ > γ c,E , we drop the assumption suppρ E = [−π, π] and look for a one-cut solution with a gap. The procedure is clear from the previous Appendices, and the result is: which, after imposing normalization and rewriting in terms of trigonometric functions, gives: supported on

Appendix B. Free energies
This Appendix contains the calculations to obtain the free energy of the unitary matrix models considered in the main text.
B.1. Free energy for Z sym.
u,H and Z sym. u,E at weak coupling. The simplest way to obtain the free energy is to evaluate its derivative with respect to the parameter t and then integrate. In the weak coupling phase, 0 ≤ γ ≤ γ c , we can use the boundary condition Z u (γ = 0) = 1, which follows immediately from the normalization of the Haar measure on U (N ). At strong coupling we use the continuity of log Z u (γ) at γ = γ c .
For what concerns F sym. u,H , at weak coupling we have: Integrating with boundary condition Z H u (γ = 0) = 1 gives For F sym. u,E at weak coupling we have: which, after integration, gives In particular, we see that the free energies of the two models are equal in the weak coupling phase. This is expected, since, in the weak coupling phase, the free energy equals the result provided by Szegő theorem, i.e. the limit without scaling (see e.g. [72] for a proof of this statement), and it is well known that the unscaled limits of the models (2.6) and (2.7) are equal.
B.2. Free energy for Z sym. u,H and Z sym. u,E at strong coupling. We now pass to the strong coupling phase, and use the form of the eigenvalue densities ρ H and ρ E at gamma > γ c .
Starting with F sym. u,H , we have: where we used the change of variables y = sin(ϕ/2) sin(φ 0 /2) and used the explicit form of sin(φ 0 /2) 2 . Immediate integration gives: with C H (γ) a t-independent integration constant fixed by continuity at γ = γ c,H . For F sym. u,E the calculations are almost the same, except for the change of sign in front of all factors of t and γ. We get: Notice that the final line here could not be inferred from the final line of F sym. u,H with reversed signs. After integration we obtain the result: where C E (γ) is t-independent and fixed by continuity at the critical value γ c,E .
B.3. Free energy for Z gen. u,H and Z gen. u,E . The free energies for the matrix models (2.13) and (2.14) are computed at weak coupling in exactly the same manner as in Appendix B.1. It is easy to check that, thanks to the parity of the integral, the free energies receive contribution from the terms proportional to v but they remain real. Direct computations give: u,E dt (γ < γ c,E ).
We proceed following standard methods for the analysis of Hermitian matrix models at large K, as reviewed for example in [74]. We introduce the resolvent where now, assuming a solution with symmetric support [−A, A], the path L is the segment [−A, A] on the real line. The eigenvalue density is recovered as the jump of the resolvent along L: with ± meaning the limit taken approaching the real line from the upper (resp. lower) half plane.
A solution for the resolvent is the following: write and the saddle point equation (C.1) becomes a discontinuity equation for Φ(z), with solution: where the contour C is a closed curve surrounding the cut L, as in Figure 10. W (u) is a meromorphic function, with simple poles at: We can thus deform the contour, to avoid the branch cut of the square root, and pick the poles of the integrand. This leaves a residual integral along an infinitely large circle. This latter contribution vanishes, since W (u) ∼ 1/u at large |u|. The pole at u = z generates the regular part of Ψ(z), thus yields no contribution to the eigenvalue density. Therefore, the relevant contributions to the ρ arise from the poles of W (u) in the complex plane. For the calculations, one has to be careful with the signs in front of the square roots, according to the definition of h(z). With our conventions, After simple computations, one arrives to: The eigenvalue density is obtained this function multiplied by the jump of h(z): .
The value of the boundary A can be fixed by normalization, or equivalently looking at the large z behaviour of the resolvent Ψ(z). From the definition, one has Ψ(z) ∼ 1/z, as z → ∞, while from the explicit evaluation, we have: This provides the solution (recall that t 0 = (1 + t)/(1 − t)) This is a positive quantity, and thus ρ stereo. is supported on the real line, as long as This is consistent with our general analysis: due to the change of topology, the stereographic projection should only provide the gapped phase of the unitary matrix model. We also notice that, undoing the stereographic projection, we obtain: thus the boundary we obtain on the real line is in fact the stereographic projection of the boundary of the model on the unit circle, as expected.
C.2. Solution for Z gen. E,stereo. . In this Appendix, we study the large K limit of the matrix model Z gen.
E,stereo. defined in eq. (3.4). The saddle point equation is: with W (x) in the present, more general case given by: . From this expression, it is clear that the asymmetry in the integrand introduces an imaginary part in the eigenvalue density, but does not introduce new poles of W (u) in the complex plane. We proceed as in the previous Appendix, and find: where now the residues at u = ±i, ±it 0 yield an extra term, proportional to the parameter v.
After some simple calculations, we see that this extra contribution is: The final expression for the eigenvalue density is: , Imposing normalization, from the parity of the integral we obtain the same value of A as above.
C.3. Solution for Z sym. E,stereo. with modified symbol. We now solve the saddle point equation (C.1) with, on the right hand side, the modified function The procedure follows closely that of Appendix C.1, but now we do not assume symmetry of the support for ρ stereo. , and let We therefore modify the definition of the auxiliary function h(z) = (−A − z)(B − z).
For later convenience, we also introduce the functions while from order z 1 , and using the previous expression to simplify the equation, we get the constraint γt 0 2h(t 0 )