Orthogonal polynomials for a class of measures with discrete rotational symmetries in the complex plane

We obtain the strong asymptotics of polynomials $p_n(\lambda)$, $\lambda\in\mathbb{C}$, orthogonal with respect to measures in the complex plane of the form $$ e^{-N(|\lambda|^{2s}-t\lambda^s-\overline{t\lambda}^s)}dA(\lambda), $$ where $s$ is a positive integer, $t$ is a complex parameter and $dA$ stands for the area measure in the plane. Such problem has its origin from normal matrix models. We study the asymptotic behaviour of $p_n(\lambda)$ in the limit $n,N\to\infty$ in such a way that $n/N\to T$ constant. Such asymptotic behaviour has two distinguished regimes according to the topology of the limiting support of the eigenvalue distribution of the normal matrix model. If $0<|t|^2T/s$ the eigenvalue distribution support consists of $s$ connected components. Correspondingly the support of the limiting zero distribution of the orthogonal polynomials consists of a closed contour contained in each connected component. Our asymptotic analysis is obtained by reducing the planar orthogonality conditions of the polynomials to an equivalent system of contour integral orthogonality conditions. The strong asymptotics for the orthogonal polynomials is obtained from the corresponding Riemann--Hilbert problem by the Deift--Zhou nonlinear steepest descent method.


Introduction
We study the asymptotics of orthogonal polynomials with respect to a family of measures supported on the whole complex plane. To set up the notation for the general case, let p n (λ) denote the monic orthogonal polynomials of degree n such that C p n (λ)p m (λ)e −N W (λ) dA(λ) = h n,N δ nm , n, m = 0, 1, 2, . . . where W : C → R is called the external potential, N is a positive parameter 1 , dA(λ) is the area measure in the complex plane and h n,N is the norming constant. The external potential is assumed to have sufficient growth at infinity so that the integrals in (1.1) are bounded. Planar orthogonal polynomials satisfying (1.1) appear naturally in the context of normal matrix models [12] where one studies probability distributions of the form M → 1 Z n,N e −N Tr(W (M )) dM, Z n,N = Nn e −N Tr(W (M )) dM, (1.2) where N N is the algebraic variety of n × n normal matrices Since normal matrices are diagonalizable by unitary transformations, the probability density (1.2) can be reduced to the form [43] 1 Z n,N i<j |λ i − λ j | 2 e −N n j=1 Tr W (λi) dA(λ 1 ) · · · dA(λ n ), where λ j are the complex eigenvalues of the normal matrix M and the normalizing factor Z n,N , called partition function, is given by The statistical quantities related to eigenvalues can be expressed in terms of the orthogonal polynomials p n (λ) defined in (1.1). In particular, the average density of eigenvalues is ρ n,N (λ) = 1 n e −N W (λ) The partition function Z n,N can be written as a product of the normalizing constants Z n,N = n j=0 h j,N .
While the asymptotic density of eigenvalues can be studied using an approach from potential theory [49], the zero distribution of orthogonal polynomials remains an open issue for general potential weights despite general results in [51]. The density of eigenvalues ρ n,N (λ) converges (in the sense of measures) in the limit to the unique probability measure µ * in the plane which minimizes the functional [25], [35] I(µ) = log |λ − η| −1 dµ(λ)dµ(η) + 1 T W (λ)dµ(λ). (1.4) The functional I(µ) in (1.4) is the Coulomb energy functional in two dimensions and the existence of a unique minimizer is a well-established fact under mild assumptions on the potential W (λ) [49]. If W is twice continuously differentiable and its Laplacian ∆W is non-negative, the equilibrium measure is given by dµ * (λ) = ∆W (λ)χ D (λ)dA(λ), where χ D is the characteristic function of the compact support set D = supp(µ * ). Sub-leading order corrections to the behaviour of the eigenvalues distribution ρ n,N (λ) as n, N → ∞ and fluctuations in the bulk and at the boundary of the support D have been considered in [9], [3], [42] [4], [5].
The measure µ * can be also uniquely characterized by the Euler-Lagrange conditions for all values λ ∈ C with equality in (1.5) on the support 2 of µ * . The Lagrange multiplier D is called the (generalized) Robin constant. It is a non-trivial problem to determine the shape of the support set D. In some cases this problem is called Laplacian growth. When the potential W (λ) is real analytic the boundary ∂D is a finite union of analytic arcs with at most a finite number of singularities [48], see also [35]. There is only a handful of potentials W (λ) for which the polynomials p n (λ) can be explicitly computed. The simplest example is W (λ) = |λ| 2 for which the orthogonal polynomials p n (λ) are monomials of degree n, the constants h n,N and the average density of eigenvalues ρ n,N can be computed explicitly in terms of the Gamma function. The matrix model associated to this potential is known as the Ginibre ensemble [31,30] and the density of eigenvalues ρ n,N converges to the normalized area measure on the disk of radius √ T centered at the origin. In general, for radially symmetric potentials W = W (|λ|), the orthogonal polynomials are always monomials and in the limit (1.3), the eigenvalue distribution is supported either on a disk or an annulus by the Single-ring Theorem of Feinberg and Zee [28], whose rigorous proof can be found in [32]. The correlation functions in this case have been studied in [12].
The harmonic deformation of the Gaussian case W (λ) = |λ| 2 +(harmonic) has been intensively studied. In particular the potential W (λ) = |λ| 2 − t(λ 2 +λ 2 ) is associated to the Hermite polynomials for |t| < 1/2. In the limit (1.3) the distribution of eigenvalues is the normalized area measure on an ellipse, while the distribution of the zeros of the orthogonal polynomials is given by the (rescaled) Wigner semicircle law with support between the two foci of the ellipse [23]. The normal matrix model with a general deformation W (λ) = |λ| 2 + Re(P (λ)) where P (λ) is a polynomial of a fixed degree has first been considered in the seminal paper [47], [54], (see also the review article [56]) where the connection with the Hele-Shaw problem and integrable structure in conformal dynamic has been pointed out. More general potentials have been considered later in [55]. For such potentials, however, the matrix integrals have convergence issues in the complex plane and therefore a natural cut-off has been introduced in the work of Elbau and Felder [25]. In [24] the polynomials associated to such deformed Gaussian potentials have been studied and it was argued that the Cauchy transform of the limiting zero distribution of the orthogonal polynomials coincides with the Cauchy transform of the limiting eigenvalue distribution of the matrix model outside the support of the eigenvalues. Moreover, it was also conjectured in [24] that the zero distribution of the polynomials p n (λ) is supported on tree-like segments (the mother body, see definition below) inside the compact set (the droplet) that attracts the eigenvalues of the normal matrix model.
For the external potential W (λ) = |λ| 2 + Re(tλ 3 ), Bleher and Kuijlaars [11] defined polynomials orthogonal with respect to a system of unbounded contours on the complex plane, without any cut-off and which satisfy the same recurrence relation that is asymptotically valid for the orthogonal polynomials of Elbau and Felder. They then study the asymptotic distribution of the zeros of such polynomials confirming the predictions of [24]. Similar results were obtained for the more general external potential [38] W (λ) = |λ| 2 − Re(t|λ| k ), k ≥ 2 and |t| sufficiently small so that the eigenvalue distribution of the matrix model has an analytic simply connected support. Cases in which the eigenvalue support has singularities were analyzed in [39] and [6]. In particular in the work [6] the external potential W (λ) = |λ| 2 −2c log |λ−a| with c and a positive constants, has been studied and the strong asymptotics of the corresponding orthogonal polynomials has been derived both in the case in which the support of the eigenvalues distribution is simply connected (pre-critical case) or multiply connected (post-critical case) and critical transition was observed (see also [53,7]). We remark that in the work [6], differently form the previous works, the zeros of the orthogonal polynomials do not always accumulate on a curve that corresponds to the mother-body of the domain where the eigenvalues of the normal matrix models are distributed.
In this work we study the strong asymptotic of the polynomials p n (λ) orthogonal with respect to a density e −N W (λ) where the external potential is of the form W (λ) = |λ| 2s − tλ s −tλ s with λ ∈ C, s a positive integer and t ∈ C \ {0}. By a simple rotation of the variable λ the analysis can be reduced to the case of real and positive t. Therefore, without loss of generality, we may and do assume that t ∈ R + , that is, (1.6) and the associated orthogonality measure is Note that the potential W (λ) has a discrete rotational Z s -symmetry. It was observed in [8] (see also [26]) that if a potential W (λ) can be written in the form the equilibrium measure for W can be obtained from the equilibrium measure of Q by an unfolding procedure. In our particular case corresponds to the Ginibre ensemble so that the equilibrium measure for the potential Q is the normalized area measure of the disk where T has been defined in (1.3). The equilibrium measure for W turns out the be equal to where χ D is the characteristic function of the support set We observe that for t < t c the equation (1.9) describes a simply connected domain in the complex plane with uniformizing map from the exterior of the unit disk in the ξ-plane to the exterior of D given by For t > t c the domain defined by the equation (1.9) consists of s connected components which have a discrete rotational symmetry. For s = 2 the domain D is called Cassini oval. The boundary of D, namely ∂D can be identified with the real ovals of the Riemann surface S defined by S : c } which has genus (s − 1) 2 . Such Riemann surface does not coincide with the Schottky double of D for s > 2.
The boundary of the domain D can be also expressed by the equation The function S(λ) is analytic in a neighbourhood of ∂D and it is called the Schwarz function associated to ∂D (see e.g. [18]).
Remark 1.1. The domain D is a quadrature domain [33] with respect to the measure dµ * . Indeed for any function h(λ) analytic in a neighbourhood of D one has, by applying Stokes' theorem and the Residue Theorem, where S(λ) is the Schwarz function (1.10), c k = 1 s and λ k = t 1 s e 2πik s .

Statement of results
The goal of this manuscript is to determine pointwise asymptotics of the polynomials p n (λ) defined in (1.1) orthogonal with respect to the weight (1.7) in the two cases • pre-critical: t < t c ; • post-critical: t > t c .
The Z s -symmetry of the orthogonality measure (1.7) is inherited by the corresponding orthogonal polynomials. Indeed the non-trivial orthogonality relations are where k and l are such that n = ks + l, 0 ≤ l ≤ s − 1, 5 i.e., the n-th monic orthogonal polynomial satisfies the relation It follows that there exists a monic polynomial q (l) k of degree k such that Therefore the sequence of orthogonal polynomials {p n (λ)} ∞ n=0 can be split into s subsequences labelled by the remainder l ≡ n mod s, and the asymptotics along the different subsequences can be studied via the sequences of reduced polynomials By a simple change of coordinates it is easy to see that the monic polynomials in the sequence {q are orthogonal with respect to the measure namely they satisfy the orthogonality relations As a result of this symmetry reduction, starting from the class of measures (1.7), it is sufficient to consider the orthogonal polynomials with respect to the family of measures (1.12). It is clear from the above relation that for l = s − 1 one has γ = 0 and the polynomials q It follows that the monic polynomials p ks+s−1 (λ) have the form Remark 1.2. Observe that the weight in the orthogonality relation (1.13) can be written in the form and it is similar to the weight e −N W (u) with W (u) = |u| 2 − c log |u − a| with c > 0 studied in [6]. However in our case it turns out that c = −2γ/N < 0, so the point interaction near a = 0 is repulsive and the asymptotic distribution of the zeros of the polynomials (1.13) turns out to be substantially different from the one in [6].

6
Let us consider the level curveĈ r The level curvesĈ r consist of s closed contours contained in the set D, where D has been defined in (1.9). For these curves we consider the usual counter-clockwise orientation. Define the measureν associated with this family of curves given by dν = 1 2πis dφ r (λ), (1.16) and supported onĈ r .
The a-priori complex measure dν in (1.16) is a probability measure on the contourĈ r defined in (1.15) for 0 < r ≤ t tc . Let us denote by ν(p n ) the zero counting measure associated with the polynomials p n , namely where δ λ is point distribution with total mass one at the point λ. • for n = ks + l, l = 0, . . . , s − 2 the polynomial p n (λ) has a zero in λ = 0 with multiplicity l and the remaining zeros in the limit n, N → ∞ such that N = n − l T accumulates on the level curvesĈ r as in (1.15) with r = 1 for t < t c and r = t 2 c t 2 for t > t c . Namely the curveĈ on which the zeros accumulate is given bŷ with |λ s − t| ≤ z 0 t. The measureν in (1.16) is the weak-star limit of the normalized zero counting measure ν n of the polynomials p n for n = sk + l, l = 0, . . . , s − 2.  with z 0 = t 2 c t 2 and |z| ≤ z 0 . The curve C is similar to the Szegő curve {z ∈ C : |ze 1−z | = 1, |z| ≤ 1} that was first observed in relation to the zeros of the Taylor polynomials of the exponential function [52], and coincides exactly with such curve in the critical case z 0 = 1. The Szegö curve also appeared in the asymptotic analysis of the generalized Laguerre polynomials, see e.g. [45], [13], [40]. The curve (1.18) is the limiting curve for the zeros of the polynomials where q   for t > t c , and the contourĈ defined in (1.17), the equation defines the boundary of a domain which coincides with D defined in (1.9). The measure µ * in (1.8) of the eigenvalue distribution of the normal matrix model and the measureν of the zero distribution of the 8 orthogonal polynomials are related by (1.20) Remark 1.7. The identities (1.19) and (1.20) in Lemma 1.6 are expected to hold in general for a large class of normal matrix models. It has been verified for several other potentials (see, for example [2,24,6,11,38]). We also observe that for the orthogonal polynomials appearing in random matrices, in some cases, the asymptotic distribution of the zeros is supported on the so-called mother body or potential theoretic skeleton of the support D of the eigenvalue distribution. We recall that a measure ν is a strong mother body for a domain D with respect to a measure µ if [34] 1) log |λ − η|dµ(η) ≤ log |λ − η|dν(η), for λ ∈ C with equality for λ outside D; 2) ν ≥ 0 and supp ν ⊂ supp µ; 3) the support of ν has zero area measure; 4) the support of ν does not disconnect any part of D fromD c .
If the measure ν has only the property 1), 2) and 3) it is called weak mother body. The problem of constructing mother bodies is not always solvable and the solution is not always unique [50].
Concerning the explicit examples appearing in the random matrix literature, for the exponential weight W (λ) = |λ| 2 + Re(P (λ)) where P (λ) is polynomial, the support of the zero distribution of the orthogonal polynomials is a strong mother body of the domain that corresponds to the eigenvalue distribution of the matrix model (see, e.g., [11,38,31,25,54]). In contrast, in the model studied in [6] and also in the present case, the support of the zero distribution of the orthogonal polynomials does not have property 4) and therefore it is a weak mother body of the set D.
The proof of Theorem 1.4 is obtained from the strong and uniform asymptotics of the polynomials p n (λ) in the whole complex plane which is obtained by characterising the orthogonal polynomials p n (λ) via a Riemann-Hilbert method. Uniform Asymptotics. In the next theorem we describe the strong and uniform asymptotic of the polynomials p n (λ) in the complex plane. We distinguish the pre-critical case t < t c and post-critical case t > t c . We first define the function whereφ r (λ) is given by (1.14).
We observe that in compact subsets of the exterior ofĈ there are no zeros of the polynomials p n (λ). The only possible zeros are located in λ = 0 and in the region where the second term in parenthesis in the expression (1.23) is of order one. Since Reφ(λ) is negative insideĈ and positive outsideĈ it follows that the possible zeros of p n (λ) lie insideĈ and are determined by the condition The expansion (1.23) shows that the zeros of the polynomials p n (λ) are within a distance O 1/k 2 from the level curve (1.25). This level curve approachesĈ defined in (1.17) at a rate O (log k/k). For t > t c the polynomials p n (λ) with n = ks + l, l = 0, . . . , s − 2, have the following asymptotic behaviour. Theorem 1.9 (Post-critical case). For t > t c the polynomials p n (λ) with n = ks + l, l = 0, . . . , s − 2, and γ = s−l−1 s ∈ (0, 1) have the following behaviour when n, N → ∞ in such a way that N T = n − l (1) for λ in compact subsets of the exterior ofĈ one has (2) for λ in the region nearĈ and away from the points λ s = t(1 − z 0 ) one has withφ(λ) defined in (1.21) and the constant c defined in (4.24); (3) for λ in compact subsets of the interior region ofĈ one has (4) in the neighbourhood of each of the points that solve the equation where U(a; ξ) is the parabolic cylinder function (Chap. 13, [1]) satisfying the equation We observe that in compact subsets of the exterior ofĈ the polynomials p n (λ) have zero at λ = 0 with multiplicity l. The other possible zeros are located in the region where the second term in parenthesis in the expression (1.26) is of order one. This happens in the region where Reφ(λ) < 0 and Reφ(λ) = O (log k/k), namely in a region inside the contourĈ and defined by the equation (1.27) The expansion (1.26) shows that the zeros of the polynomials p n (λ) are within a distance O 1/k 2 from the level curve (1.27). This level curve approachesĈ defined in (1.17) at a rate O (log k/k).
The proofs of Theorem 1.8 and Theorem 1.9 are obtained by reducing the planar orthogonality relations of the polynomials p n (λ) to orthogonality relations with respect to a complex density on a contour. More precisely, the sequence of polynomials p n (λ) can be reduced to s families of polynomials q (l) k (λ s ), l = 0, . . . , s − 1 as in (1.11) with n = ks + l. The orthogonality relations of the polynomials q (l) (u) are reduced to orthogonality relations on a contour. We then reformulate such orthogonality relation as a Riemann-Hilbert problem. We perform the asymptotic analysis of the polynomials q (l) k (λ) using the nonlinear steepest descent/stationary phase method introduced by Deift and Zhou [22] and successfully applied in Hermitian random matrices and orthogonal polynomials on the line in the seminal papers, [19], [20]. See [21] for an introduction to this method, with a special emphasis on the applications to orthogonal polynomials and random matrices. It would be also interesting to explore the asymptotic for the polynomials p n (λ) using the∂-problem introduced in [36].
The zeros of p n (λ) accumulate along an open contour as shown in Figure 1. The determination of this contour is a first step in the analysis. Riemann-Hilbert analysis in which a contour selection method was required, also appeared in the papers [6], [40].
The Riemann-Hilbert method gives strong and uniform asymptotics of the polynomials p n (λ) in the whole complex plane. The asymptotic behaviour of the polynomials p n (λ) in the leading and sub-leading order can be expressed in terms of elementary functions in the pre-critical case t < t c , while in the postcritical case t > t c we have used, in some regions of the complex plane, parabolic cylinder functions as in [22], [37], [46]. The proof of Theorem 1.4 is then deduced from the strong asymptotic of the orthogonal polynomials.
The paper is organized as follows: • In Sect. 2, by using the symmetry properties of the external potential, we reduce the orthogonality relations of the polynomials p n (λ) on a contour and recall how to formulate a Riemann-Hilbert problem for such orthogonal polynomials [27].
• In Sect. 3 we obtain the strong and uniform asymptotics for the polynomials p n (λ) in the pre-critical case t < t c .
• In Sect. 4 we find the strong and uniform asymptotics for the polynomials p n (λ) in the post-critical case t > t c and we prove Theorem 1.4.
In the critical case t = t c , the set D in (1.9) that supports the eigenvalue distribution of the normal matrix models has a singularity in λ = 0. The corresponding asymptotics of the orthogonal polynomials seems to be described in terms of Painlevé IV equation as in [17] and this situation is quite different from the generic singularity that have been described in terms of the Painlevé I equation [41], [39]. Such problem will be investigated in a subsequent publication.
Furthermore, from the norming constants of the orthogonal polynomials, one may consider the problem of determining the asymptotic expansion of the partition function Z n,N in the spirit of [16]. Also this problem will be investigated in a subsequent publication.

The associated Riemann-Hilbert problem
In this section we set up the Riemann-Hilbert problem to study the asymptotic behaviour of the orthogonal polynomials p n (λ). As seen above, the analysis of p n (λ) can be reduced to that of the polynomials q

Reduction to contour integral orthogonality
The crucial step in the present analysis is to replace the two-dimensional integral conditions (2.1) by an equivalent set of linear constraints in terms of contour integrals. In what follows it will be advantageous to perform the change of coordinate and characterize q in the new variable z. The polynomial π k (z) is also a monic polynomial of degree k and it can be characterized as follows.
Theorem 2.1. The polynomial π k (z) is characterized by the non-hermitian orthogonality relations where Σ is a simple positively oriented contour encircling z = 0 and z = 1 and the function z z−1 γ is analytic in C\[0, 1] and tends to one for |z| → ∞. Proof. In order to prove the theorem we first show that the orthogonality relation for the polynomials q l k (u) on the plane can be reduced to an orthogonality relation on a contour. To this end we seek a function χ j (u,ū) that solves the∂-problem Having such a function, for any polynomial q(u) one has where d denotes the operation of exterior differentiation. If such χ j (u, u) exists, one can use Stokes' theorem and reduce the planar orthogonality relation to an orthogonality relation on a suitable contour. The equation (2.3) has a contour integral solution |u| → ∞ It follows that for any polynomial q(u) the following integral identity holds: where R and R 0 are sufficiently large and G j (u) = Γ(j−γ+1) . So it follows that for any polynomial q(u) the following identity is satisfied: where γ ∈ (0, 1), j is an arbitrary non-negative integer, andΣ is a positively oriented simple closed loop enclosing u = 0 and u = t. Making the change of coordinate u = −t(z − 1) one arrives to the statement of the theorem. Q.E.D.

The Riemann-Hilbert problem
Our aim is to study the behaviour of the polynomials π k (z) in the limit k → ∞ and N → ∞ in such a way that for n = ks + l one has Let In terms of the weight function In the limit k → ∞, we distinguish two different cases: • pre-critical case z 0 > 1, corresponding to 0 < t < t c , • post-critical case z 0 < 1, corresponding to t c < t .
Our goal now is to characterize the polynomial π k (z) as a particular entry of the unique solution of a matrix-valued Riemann-Hilbert problem. Let us first define the complex moments where, for simplicity, the dependence on k is suppressed in the notation. Introduce the auxiliary polynomial Note that Π k−1 is not necessarily monic and its degree may be less than k − 1: its existence is guaranteed just by requiring that the determinant in the denominator does not vanish.
Proof. We have where the last identity has been obtained by the reflection of the column index j → k − 1 − j. Due to Theorem 2.1 we have an hence the second determinant is given by (2.6) Finally, the determinant on the right-hand side is strictly positive because where the equality follows from the fact that the columns of the two matrices are related by a unimodular triangular matrix, while the inequality follows from the positivity of the measure. Finally, since Γ(z) has no zeros (and no poles since j − γ + 1 > 0), the non-vanishing follows from (2.6).
Q.E.D. Define the matrix It is easy to verify that the matrix Y (z) is the unique solution of the following Riemann-Hilbert problem of Fokas-Its-Kitaev-type [27]: 1. Piecewise analyticity: 2. Jump on Σ: 3. Behaviour at z = ∞: The last relation is obtained by noticing that if Π k−1 is given by (2.5) then the entry

Initial undressing step
In order to simplify the subsequent analysis, we define the following modified matrix: with V (z) as (2.4).

Asymptotic analysis in the pre-critical case
In order to analyse the large k behaviour ofỸ (z) we use the Deift-Zhou nonlinear steepest descent method [22]. The first step to study the large k behaviour of the matrix functionỸ (z) is to make a transformatioñ Y (z) → U (z) so that the Riemann-Hilbert problem for U (z) is normalised to the identity as |z| → ∞. For this purpose we introduce a contour C homotopically equivalent to Σ in C \ [0, 1] and a function g(z) analytic off C. Both the contour C and the function g(z) will be determined later. We assume that the function g(z) is of the form where dν(s) is a positive measure with support on C such that C dν(s) = 1.
With this assumption clearly one has where the logarithm is branched on the positive real axis.

First transformationỸ → U
Since C is homotopically equivalent to Σ in C \ [0, 1] we can deform the contour Σ appearing in the Riemann-Hilbert problem forỸ to C. Define the modified matrix where is a real number, to be determined below. Then U (z) solves the following RHP problem 1. Piecewise analyticity: 2. Jump discontinuity on C: 3. Jump discontinuity on (0, 1): 4. Endpoint behaviour at z = 0 and z = 1:

Large z boundary behaviour:
The polynomial π k (z) is determined from U (z) by (3.5) 18

The choice of g-function
In order to determine the function g and the contour C we impose that the jump matrix (3.3) becomes purely oscillatory for large k. This is accomplished if the following conditions are satisfied: Next we show that we can find a function g and a contour C that satisfy the conditions (3.6). For the purpose we use the following elementary result in the theory of boundary value problems.
can be represented in the form The boundary values of the function Φ on the two sides of the contour L then satisfy Now we apply Lemma 3.1 to the function g (z) that satisfies the differentiated boundary condition The functions ψ ± (z) of Lemma 3.1 are ψ + (z) = 1 z 0 and ψ − (z) = − 1 z and therefore the function g (z) takes the form so that the measure dν in (3.1) is given by Integrating the relation (3.7) and using (3.2) one has where is an integration constant and log z is analytic in C\R + . Performing the integral in (3.1) for a specific value of z ∈ C, say z = 0 and deforming C to a circle of radius r one can determine the value of : The total integral of dν in (3.8) is normalized to one on any closed contour containing the point z = 0. However we have to define the contour C so that dν(z) is a real and positive measure along C. To this end we introduce the function (3.10) Observe that φ r (z) is analytic across C, namely φ r+ (z) = φ r− (z), z ∈ C and that Figure 4: The family of contours C r for three values of r: r = z 0 , r = 1 and r < 1. For r = z 0 , the region |z| > z 0 is plotted as well.
The next identity follows in a trivial way dν(z) = 1 2πi dφ r (z).
This equation defines a family of contours which are closed for |z| ≤ r (easy to verify). Since the function − log r + r z 0 for r > 0 has a single minimum at r = z 0 and diverge to +∞ for r → 0 and r → ∞, it is sufficient to consider only the values 0 < r ≤ z 0 . We define the contour C r associated to r as (see Figure 4) Note that the point z = 0 lies inside C r since Re φ(z) → −∞ for |z| → 0. Furthermore C r intersects the real line in the points z = r and z ∈ (−r, 0). Indeed Re φ(−r) > 0 which shows that the point z = −r lies outside C r .
Lemma 3.2. The a-priori complex measure ν in (3.8) is a probability measure on the contour C r defined in (3.12) for 0 < r ≤ z 0 .
Proof. By the Residue Theorem the measure dν has total mass 1 on any contour C r . Since Re φ r = 0 on the contour C r it follows that the measure dν = 1 2πi dφ r is real on the contour C r . In order to show that the measure is positive on the contour C r we introduce the variable ψ r = e φr , On the contour C r we have that |ψ r | = 1, so that in the ψ r -plane the contour C r is mapped to a circle of radius one and the map ψ r = e φr is a univalent conformal map from the interior of C r to the interior of a circle [52]. We have which shows that the measure dν in the variable ψ r is a uniform measure on the circle and therefore dν is a positive measure on each contour C r . Q.E.D. We are now ready to prove Lemma 1.3 and Lemma 1.6 announced in the Introduction.
Proof of Lemma 1.3. By using the Residue Theorem, it is straightforward to check that dν is a probability measure on any contourĈ r defined in (1.15). In the post-critical case such curveĈ r has s connected componentsĈ j r , j = 0, . . . , s − 1, each encircling exactly one root of the equation λ s = t. In the pre-critical case we denote with the same symbolĈ j r the connected components ofĈ r \ {0} together with {0}.
This relation, together with Lemma 3.2, where we proved that dν(z) is real and positive on C r , implies that also dν is real and positive on anyĈ j r and hence on the wholeĈ r . Q.E.D.
Proof of Lemma 1.6. By using the Residue Theorem we first calculate the r.h.s. of (1.19) obtaining which is equivalent to the equation for the boundary of D defined in (1.9). Next we show that the l.h.s. of (1.20) is equal to its r.h.s. By using Stokes' Theorem, the relation (1.10) and the Residue Theorem one obtains which coincides with (3.13).
Q.E.D. Figure 5: The contour C r<1 ∪C that is homotopic to Σ. The region where Re φ(z) < 0 is coloured in pink. Color figure online.

Choice of the contour
We now argue that the relevant contour on which the zeros of the orthogonal polynomials π k (z) accumulate in the precritical case is given by the level r = 1. The family of contours 0 < r < 1 can be immediately ignored because in this case Σ has to be deformed to two contours Σ C r<1 ∪C as shown in Figure 5. On the contourC we have Re(φ) > 0 and it is not possible to perform any contour deformation to get exponentially small terms in the jump matrix (3.3) as k → ∞.
For 1 < r ≤ z 0 the original contour Σ can be deformed homotopically to the contours C r . In this case the asymptotic expansion of the matrix entry Y 11 (z) = π k (z) as k → ∞, does not give any sub-leading contribution. For this reason we discard this case and we omit the corresponding asymptotic analysis. The only possible case remains r = 1. In this case the analysis as k → ∞ to the matrix entry Y 11 (z) = π k (z) which will performed in subsections 3.2 to 3.5, gives leading and subleading terms. The comparison of these two terms enable us to locate the zeros of the polynomial π k (z) on a contour that lies within a distance O (log k/k) from the contour C r=1 . So for the reasons explained above, we are going to perform the asymptotic analysis of the RH-problem (3.3)-(3.4) by deforming the contour Σ to the contour C r=1 . For simplicity we denote this contour by C: Also, for simplicity, the function φ r=1 will be referred to as φ below: where log z is analytic in C\(−∞, 0]. With this choice of contour, the jump matrices of the RH problem (3.3)-(3.4) for the matrix U (z) are summarized in Figure 6. Figure 7: The jump matrices for T (z) with con-

The second transformation
Consider two extra loops C i and C e as shown in Figure 7. These define new domains Ω 0 , Ω 1 , Ω 2 and Ω ∞ . Define the new matrix-valued function Then T (z) satisfies the following Riemann-Hilbert problem: 1. Piecewise analyticity: 3. Endpoint behaviour at z = 0 and z = 1:

Large z boundary behaviour:
Note that e kφ(z) is analytic in C \ {0} . The important feature of this Riemann-Hilbert problem is that the jumps are either constant or they tend to the identity matrix as k → ∞ at an exponential rate. as k → ∞ uniformly for z ∈ C e ∪ C i \U 1 , where U 1 is a small neighbourhood of z = 1.
The proof of this proposition follows immediately from the fact that Re φ(z) < 0 for z ∈ Int(C)\U 1 and Re φ(z) > 0 for z ∈ Ext(C)\U 1 . Here Int(C) is the interior region bounded by C and Ext(C) is the exterior region bounded by C.
From the above proposition it follows that e −γπiσ3 as z ∈ (0, 1) I as z ∈ C e ∪ C i .

The outer parametrix for large z
We need to find a matrix-valued function P ∞ (z) analytic in C \ (C ∪ [0, 1]) such that it has jump discontinuities given by the matrix v ∞ (z) in (3.16), namely with endpoint behaviour and large z boundary behaviour (3.17) The k-independent matrixP ∞ (z) has no jump on C and it satisfies the following Riemann-Hilbert problem: 1. Piecewise analyticity:P ∞ is holomorphic in C \ (0, 1) 2. Jump discontinuity on (0, 1): 3. Boundary behaviour at z = ∞:P The solution to this Riemann-Hilbert problem is given bỹ which leads to the particular solution z ∈ Ext(C) .

The local parametrix at z = 1
The aim of this section is to construct a local parametrix P 0 (z) in a small neighbourhood U 1 of z = 1 having the same jump property as T for z near 1 and matching the outer parametric P ∞ (z) in the limit k → ∞ and z ∈ ∂U 1 . Then the Riemann-Hilbert problem for P 0 (z) is given by and P 0 (z) = P ∞ (z)(I + o(1)) as k → ∞ and z ∈ ∂U 1 . (3.19) In order to build such local parametrix near the point z = 1 we first construct a new matrix function B(z) from P 0 (z): The matrix B(z) satisfies the following jump relations in a neighbourhood of 1: Next we construct a solution to the so called model problem, namely a Riemann-Hilbert problem that has the same jumps as B.

Model problem
Consider the model problem for the 2×2 matrix function Ψ(ξ) analytic in C\R − with boundary behaviour we obtain the following Riemann-Hilbert problem forΨ(ξ): where the function ξ γ is analytic in C\(−∞, 0] and (ξ γ ) + denotes the boundary from the right with respect to the contour (−∞, 0) oriented in the positive direction. This is an abelian Riemann-Hilbert problem, and therefore easily solvable by the Sokhotski-Plemelj formula: In particular . (3.24)

Construction of the parametrix
We are now ready to specify the parametric P 0 (z) which follows from combining (3.23) with (3.20), that is where Q(z) and χ(z) are defined in (3.21)) and (3.17) respectively and E(z) is an analytic matrix in a neighbourhood of U 1 and w(z) is a conformal mapping from a neighbourhood of 1 to a neighbourhood of 0. The conformal map w(z) is specified by We observe that The matrix E(z) is obtained from condition (3.19) which, when combined with (3.24) gives where we use the fact that Q(z) → I exponentially fast as k → ∞, and z ∈ ∂U 1 . From the above expression it turns out that the matrix E(z) takes the form We observe that the function E(z) is single valued in a neighbourhood of U 1 . Indeed the jumps of (z−1) and w(z) − γ 2 σ3 cancel each other. From the expression (3.25) and (3.26) the matching between P 0 (z) and P ∞ (z) takes the form for z ∈ ∂U 1 , where γ ∈ (0, 1). Figure 8: The jump contour structure C R We now define the error matrix R in two regions of the plane, using our approximations to the matrix T . Set

Riemann-Hilbert problem for the error matrix
The matrix R is piecewise analytic in C with a jump across the contour C R = C i ∪ C e ∪ ∂U 1 given in Fig. 8.
For z ∈ C R , we have The jump matrices across the contour C R \ ∂U 1 are all exponentially close to I for large k because v T converges exponentially fast to v ∞ defined in (3.16) and the product (3.18). The only jump that is not exponentially small is the one on ∂U 1 since from (3.27) for any integer M > 2 we have with the shift matrix Λ defined in (3.23). By employing the notation v (j) we can rewrite the matrix v R (z) in the form where, in particular, (3.24) implies v (1) By a standard perturbation theory argument one has the expansion which gives, using (3.29), (3.30) and (3.32)

R
(1) where we observe that the function v By a simple residue calculation we obtain (3.33) In general, given the structure of the jump matrix (3.30) the error matrix R (j) (z) in (3.32) is of the form namely, only the (1, 2)-entry of the matrices R (j) (z) is non-zero. From (3.28) and (3.32) one has where, in particular, the first term R (1) (z) is given in (3.33).
3.5 Proof of Theorem 1.8: asymptotics for p n (λ) for z 0 > 1 In order to obtain the asymptotic expansion of the polynomials p n (λ) for n → ∞, and N T = n − l we first derive the asymptotic expansion of the reduced polynomials π k (z) for k → ∞ with n = sk + l. For the purpose we use (3.5), (3.15), (3.33) and (3.35) to obtain From the above relation we obtain the expansions in the following regions.
The exterior region Ω ∞ . In this region, from (3.34) and (3.35), for any integer M ≥ 2 we have on any compact subset of Ω ∞ . Therefore there are no zeros accumulating in this region. The interior region Ω 0 \ U 1 .
The leading term of the above expansion is of order O k −1−γ , so there are no zeros in this region. The interesting region Ω 1 \ U 1 .
The region U 1 : for z ∈ U 1 ∩ Ext(C) whereΨ 12 is the (1, 2)-entry of the matrixΨ defined in (3.22) and v (1) R (z) has been defined in (3.31). Here Int(C) and Ext(C) is the interior and exterior of C respectively.
In order to obtain the asymptotic behaviour of the polynomials p n (λ) we use the substitution where the function φ(z) has been defined in (3.14). The curves in (3.39) approach C at the rate O (log k/k) and lies in Int(C). The normalised counting measure of the zeros of π k (z) converges to the probability measure ν defined in (3.8).
Proof. Observing the asymptotic expansion (3.36) of π k (z) in Ω ∞ \ U 1 it is clear that π k (z) does not have any zeros in that region, since z = 0 and z = 1 do not belong to Ω ∞ \ U 1 . The same reasoning applies to the region Ω 0 \ U 1 where there are no zeros of π k (z) for k sufficiently large.
From the relations (3.37) and (3.38) one has that in Ω 1 ∪ Ω 2 using the explicit expression of g(z) defined in (3.9) The zeros of π k (z) may only lie asymptotically where the expression is equal to zero. Since Ω 2 ⊂ {Re(φ) ≥ 0} and Ω 1 ⊂ {Re(φ) ≤ 0}, it follows that the zeros of π k (z) may lie only in the region Ω 1 and such that Re φ(z) = O (log k/k). Namely the zeros of the polynomials π k (z) lie on the curve given by (3.39) with an error of order O 1/k 2 . Such curves converge to the curve C defined (1.18) at a rate O (log k/k) (see Fig. 9). k=65 k=80 k=95 Figure 9: The zeros of π k (z) for s = 3, l = 0, t < t c , and k = 65, 80, 95. The red contour is C while the yellow contour is the curve (3.39). Color figure online. We now know that the zeros of π k (z) accumulate on the curve C. We still need to determine the asymptotic zero distribution. Due to the strong asymptotic of π k (z) we have from (3.36) uniformly on compact subsets of the exterior of C. Furthermore C is the boundary of its polynomial convex hull. Then it follows that the measure dν in (3.8) is the weak-star limit of the zeros distribution of the polynomials π k (z) (see [44], Theorem 2.3 and [49] Chapter 3). Q.E.D.

Post-critical case
In this section we assume that 0 < z 0 < 1 .
To perform the analysis similar to the pre-critical case, one has to choose the appropriate contour from the homotopy class of the contour Σ that encircles the branch cut [0, 1]. Since the point z 0 is lying in (0, 1), the family of curves C r , 0 < r ≤ z 0 defined in (3.12) are loops encircling z = 0 and crossing the real line in z = r. The regions where Re φ r (z) < 0 with φ r (z) defined in (3.10), are depicted in the figures below in red for two values of the parameter r. Figure 10: The contour C r for r = z 0 on the left and for 0 < r < z 0 on the right. The region where Re φ r (z) < 0, with φ r (z) defined in (3.10), is pink. Color figure online. Figure 11: The contour C z0 ∪ C o . The region Re φ z0 (z) < 0 is plotted in pink. Color figure online.
In order to perform the asymptotic analysis of the Riemann-Hilbert problem (2.8), (2.9) and (2.10), we need to deform Σ to a homotopic contour C r ∪ C o in such a way that the Re φ r (z) is negative on C o .
It is clear from Figure 10 that the only possibility is to deform the contour Σ to the contour C z0 ∪C o with C o as depicted in Fig. 11. For simplifying the notation we define the following.
According to (3.11) on the contour C we have

First transformationỸ → S
Since C ∪ C o is homotopically equivalent to Σ in C \ [0, 1] when z 0 < 1, we can deform the contour Σ appearing in the Riemann-Hilbert problem (2.8)-(2.10) for the matrixỸ (z) to C ∪ C o . Define the modified matrix whereỸ (z) has been defined in (2.11). Then the matrix S(z) is the unique solution of the following Riemann-Hilbert problem with standard large z behaviour at z = ∞: 1 e kφ(z) 0 1 Figure 12: The jump matrices for S(z). The function φ(z) is defined in (4.1).

Large z boundary behaviour:
The orthogonal polynomials π k (z) are recovered from the matrix S(z) using the relation

The second transformation S → T : opening of the lenses
Consider two extra loops C i and C e as shown in Figure 13. These define new domains Ω 0 , Ω 1 , Ω 2 and Ω ∞ . Define Then this matrix-valued function has the following jump discontinuities: where Σ T is the contour defined in the Figure 13. Figure 13: The jump matrices for T (z) and the contour Σ T = C ∪ C i ∪ C e ∪ C o ∪ (0, 1). Endpoint behaviour: Large z boundary behaviour: The proof of the proposition follows in a straightforward way by observing that by contruction (see Figure 11) Re φ(z) > 0 for z ∈ C e \U z0 and Re φ(z) < 0 for z ∈ C i ∪ C o \U z0 . Therefore we have as z ∈ (0, 1), In order to build such local parametrix near the point z = z 0 we first construct a new matrix function B(z) from P 0 (z). Let us first define ∆(z) = I Im(z) < 0 e −γπiσ3 Im(z) > 0, (4.6) and the matrix Q as follows (4.7) Then the matrix B(z) is defined from P 0 (z) by the relation The matrix B(z) satisfies the jump relations specified in Figure 14 in a neighbourhood of z 0 . In the next section we construct the solution of the so called model problem, namely a 2 × 2 matrix Ψ that has the same jumps as the matrix B. Figure 14: The jumps for the matrix B

Model problem
Consider the model problem for the 2×2 matrix function Ψ(ξ) analytic in C\{R ∪ iR} with jumps and boundary behaviour with Ψ 1 Ψ 2 and Ψ 3 constant matrices (independent from ξ) and where the matrix v Ψ (ξ) is defined as The solution of the Riemann-Hilbert problem (4.9) and (4.10) is obtained in the following way, [22,37,46]. Recall the parabolic cylinder equation

Construction of the local parametrix
From (4.8) and the previous section we are now ready to specify the form of the local parametrix at z = z 0 , by where E(z) is an analytic matrix in a neighbourhood of U z0 , the matrix Ψ has been defined in (4.16) and ∆(z) and Q(z) have been defined in (4.6) and (4.7) respectively. The function w(z) is a conformal mapping from a neighbourhood of z 0 to a neighbourhood of 0 and it is specified by We observe that The matrix E(z) is obtained from condition (4.5) which, when combined with (4.18), gives where we used the fact that Q(z) → χ(z) −1 exponentially fast as k → ∞, and z ∈ ∂U z0 . From the above expression and (4.5), it turns out that the matrix E(z) takes the form We observe that the function E(z) is single valued in a neighbourhood of U z0 . Indeed the boundary values of (z − z 0 ) γσ3 and w −γσ3 cancel each other. The boundary value of (z(z − 1)) γ + = (z(z − 1)) γ − e 2πiγ , so that ∆(z)(z(z − 1)) γ 2 σ3 remains single valued in a neighbourhood of z 0 . From the expression (4.21) and (4.22) the matching between P 0 (z) and P ∞ (z) takes the form P ∞ (z)(P 0 (z)) −1 = E(z) I − γ 2 where γ ∈ (0, 1). It is clear from the above expression that the (2, 1)-entry of the matrix above is not small as k → ∞. For this reason we need to introduce an improvement of the parametrix.

Improvement of the local parametrix
In order to have a uniformly small error for k → ∞ we have to modify the parametrices as in [10], [15], [14]:P ∞ (z) := I + C z − z 0 P ∞ (z) ,

40
where C is a nilpotent matrix to be determined and where the matrix Ψ 1 has been defined in (4.10) and With those improved definitions of the parametrices,P ∞ (z) andP 0 (z) have the same Riemann-Hilbert jump discontinuities as before but they might have poles at z = z 0 . Note thatÊ(z) also has a pole at z = z 0 . However we can choose C in such a way thatP 0 (z) is bounded in z 0 . This is accomplished by The improved parametrix gives the following matching betweenP ∞ (z) andP 0 (z) as k → ∞ and z 0 ∈ ∂U z0 : P ∞ (z)(P 0 (z)) −1 which shows thatP ∞ (z)(P 0 (z)) −1 = I + O (1/k α ), α > 0, as k → ∞ when z 0 ∈ ∂U z0 for γ ∈ [0, 1). Figure 15: The contour Σ R = C i ∪ C e ∪ C o ∪ ∂U z0 where C i , C e and C o are defined only in C\U z0 .

Riemann-Hilbert problem for the error matrix
We now define the error matrix R in two regions of the plane, using our approximations to the matrix T . Set (4.26) The matrix R is piecewise analytic in C with a jump across Σ R (see Fig. 15). Riemann-Hilbert problem for R

R is analytic in
3. As z → ∞, we have R(z) = I + O z −1 .
The jump matrices across the contour Σ R \∂U z0 are all exponentially close to I for large k because v T (z) converges exponentially fast to v ∞ (z) defined in (4.2) for z ∈ C\U z0 and the product with P ∞ (z) defined in (4.4) and C in (4.23). The only jump that is not exponentially small is the one on ∂U z0 . Indeed one has from (4.25) v R (z) =P ∞ (z)(P 0 (z)) −1 = I + γe(z) 2 2k 1 2 +γ w(z)β 21 where e(z) = E 11 (z)k γ/2 (4.28) and we have substituted the explicit expressions of the matrix Ψ 1 , Ψ 2 and Ψ 3 as given by (4.19).
We have the following two cases depending on the value of γ ∈ (0, 1).
Compatibility of (4.27), (4.31) and (4.34) and the jump condition R + = R −vR on ∂U z0 gives the following relations. R (i) with v (i) R (z) defined in (4.32). In addition R (i) (z) is analytic in C\∂U z0 and R (i) (∞) = 0. Since v (1) R (z) has a third order pole at z = z 0 , the unique function that satisfies those conditions is given by