Orthogonal Polynomials for a Class of Measures with Discrete Rotational Symmetries in the Complex Plane

We obtain the strong asymptotics of polynomials pn(λ)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_n(\lambda )$$\end{document}, λ∈C\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda \in {\mathbb {C}}$$\end{document}, orthogonal with respect to measures in the complex plane of the form e-N(|λ|2s-tλs-t¯λ¯s)dA(λ),\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \hbox {e}^{-N(|\lambda |^{2s}-t\lambda ^s-\overline{t}\overline{\lambda }^s)}\hbox {d}A(\lambda ), \end{aligned}$$\end{document}where s is a positive integer, t is a complex parameter, and dA\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {d}A$$\end{document} stands for the area measure in the plane. This problem has its origin in normal matrix models. We study the asymptotic behavior of pn(λ)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_n(\lambda )$$\end{document} in the limit n,N→∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n,N\rightarrow \infty $$\end{document} in such a way that n/N→T\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n/N\rightarrow T$$\end{document} constant. Such asymptotic behavior has two distinguished regimes according to the topology of the limiting support of the eigenvalues distribution of the normal matrix model. If 0<|t|2T/s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|t|^2>T/s$$\end{document}, 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 equivalent 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.

Keywords Orthogonal polynomials on the plane · Normal matrix model · Riemann-Hilbert problem · Logarithmic potential theory Mathematics Subject Classification 31A15 · 30E25 · 60K35

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, . . . , (1.1) 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 where λ j are the complex eigenvalues of the normal matrix M and the normalizing factor Z n,N , called the 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 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 nonnegative, the equilibrium measure is given by where χ D is the characteristic function of the compact support set D = supp(μ * ). Sub-leading order corrections to the behavior 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 [3][4][5]9,42]. The measure μ * can also be 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 nontrivial 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 are 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, and 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 with this potential is known as the Ginibre ensemble [30,31], 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 with 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 papers [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 with 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 [7,53]).
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 to be equal to where χ D is the characteristic function of the support set We observe that for t < t c , the Eq. (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 Eq. (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 also be expressed by the equation The function S(λ) is analytic in a neighborhood of ∂ D, and it is called the Schwarz function associated with ∂ 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: The Z s -symmetry of the orthogonality measure (1.7) is inherited by the corresponding orthogonal polynomials. Indeed, the nontrivial orthogonality relations are where k and l are such that i.e., the n-th monic orthogonal polynomial satisfies the relation p n (e 2πi s λ) = e 2πin s p n (λ).
It follows that there exists a monic polynomial q 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 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 (s−1) 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].
and for r > 0, the function (1.14) 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 counterclockwise orientation. Define the measureν associated with this family of curves given by and supported onĈ r .

Lemma 1.3
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 t c .
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 λ.

Theorem 1.4
The zeros of the polynomials p n (λ) defined in (1.1) behave as follows: for t > t c . We define asĈ the curve on which the zeros accumulate that is given bŷ (1.17) 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.
Remark 1. 5 We observe that the curve (1.17) in the rescaled variable z = 1 − λ s /t takes the form 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., [13,40,45]. The curve (1.18) is the limiting curve for the zeros of the polynomials where q (l) k (u) have been defined in (1.11). From Theorem 1.4 and (1.5), the following identity follows immediately.   [2,6,11,24,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 properties 1), 2), and 3), it is called weak mother body. The problem of constructing strong 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,25,31,38,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 characterizing 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 precritical case t < t c and post-critical case t > t c . We first define the function whereφ r (λ) is given by (1.14). (1) for λ in compact subsets of the exterior ofĈ, one has for any integer M ≥ 2, (2) for λ nearĈ and away from λ = 0, (1.23) whereφ(λ) has been defined in (1.21) and (z) = (z − 1)!; (3) for λ in compact subsets of the interior ofĈ and away from λ = 0, where the (1, 2)-entry of the matrix˜ is defined in (3.22).
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 parentheses 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) (Fig. 2). For t > t c , the polynomials p n (λ) with n = ks + l, l = 0, . . . , s − 2, have the following asymptotic behavior. 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 behavior 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 neighborhood of each of the points that solve the equation where U(a; ξ) is the parabolic cylinder function (Chap. 13, [1]) satisfying the 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 parentheses 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 also be interesting to explore the asymptotic for the polynomials p n (λ) using the∂-problem introduced in [36].
The zeros of p n (λ) accumulate along a contour in the complex plane as shown in Fig. 1 or Fig. 2. 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 behavior 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 post-critical 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 the Painlevé IV equation as in [17], and this situation is quite different from the generic singularity that has been described in terms of the Painlevé I equation [39,41]. That 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 behavior 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 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| → ∞ (Fig. 3).
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 Eq. (2.3) has a contour integral solution It follows that for any polynomial q(u), the following integral identity holds: where R and R 0 are sufficiently large and j+1 . So it follows that for any polynomial q(u), the following identity is satisfied: where γ ∈ (0, 1), j is an arbitrary nonnegative 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 at the statement of the theorem.

The Riemann-Hilbert Problem
Our aim is to study the behavior of the polynomials π k (z) in the limit k → ∞ and N → ∞ in such a way that for n = ks + l, one has Let and introduce the function In terms of the weight function the orthogonality relations (2.2) can be written in the form In the limit k → ∞, we distinguish two different cases: 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 (2.5) 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.
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 and hence the second determinant is given by 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 nonvanishing follows from (2.6).
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: Y (z) is analytic in C\ and both limits Y ± (z) exist along , 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 Y 22

Asymptotic Analysis in the Pre-critical Case
In order to analyze the large k behavior ofỸ (z), we use the Deift-Zhou nonlinear steepest descent method [22]. The first step to study the large k behavior of the matrix functionỸ (z) is to make a transformationỸ (z) → U (z) so that the Riemann-Hilbert problem for U (z) is normalized 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 behavior at z = 0 and z = 1: 5. Large z boundary behavior: (3.5)

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 this purpose, we use the following elementary result in the theory of boundary value problems.
where ψ + (ζ ) is the boundary value of a function ψ + (z) analytic for z ∈ D + and ψ − (ζ ) is the boundary value of a function ψ − (z) analytic for z ∈ D − and such that ψ − (∞) = 0. Then the Cauchy integral 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 log z z ∈ Ext(C), (3.9) 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) Fig. 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 Observe that φ r (z) is analytic across C, namely φ r + (z) = φ r − (z), z ∈ C, and that The next identity follows in a trivial way: Imposing the relation (3.6) on (3.11), one has This equation defines a family of contours that 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 diverges 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 with r as (see Fig. 4) Since ∂ x Re φ r | y=0 = |x| 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 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 .
We are now ready to prove Lemma 1.3 and Lemma 1.6 given 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Ĉ This relation, together with Lemma 3.2, where we proved that dν(z) is real and positive on C r , implies that dν is also real and positive on anyĈ j r and hence on the wholeĈ r .
Proof of Lemma 1.6 By using the residue theorem, we first calculate the r.h.s. of (1.19), obtaining (3.13)

Therefore (1.19) is equivalent to
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).

Choice of the Contour
We now argue that the relevant contour on which the zeros of the orthogonal polynomials π k (z) accumulate in the pre-critical 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 Fig. 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 be performed in Sects. 3.2 to 3.5, gives leading and be sub-leading 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 Fig. 6.

The Second Transformation U → T
Consider two extra loops C i and C e as shown in Fig. 7. These define new domains 0 , 1 , 2 , and ∞ . Define the new matrix-valued function (3.15) Then T (z) satisfies the following Riemann-Hilbert problem: 1. Piecewise analyticity: Fig. 7 The jump matrices for

Large z boundary behavior:
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. 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.

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 behavior 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 = ∞: The solution to this Riemann-Hilbert problem is given bỹ which leads to the particular solution (3.18)

The Local Parametrix at z = 1
The aim of this section is to construct a local parametrix P 0 (z) in a small neighborhood 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): (3.20) where The matrix B(z) satisfies the following jump relations in a neighborhood 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 behavior 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,

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 neighborhood of U 1 and w(z) is a conformal mapping from a neighborhood of 1 to a neighborhood 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 neighborhood of U 1 . Indeed, the jumps of (z − 1) γ 2 σ 3 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).

Riemann-Hilbert Problem for the Error Matrix R
We now define the error matrix R in two regions of the plane, using our approximations to the matrix T . Set 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.

RH Problem for R
1. R is analytic in C\C R , 2. 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) (3.18). The only jump that is not exponentially small is the one on ∂U 1 , since from (3.27), for any integer 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 Therefore, where we observe that the function v (1) R (z) has a simple pole in z = 1 with expansion v (1) 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 35) where, in particular, the first term R (1) (z) is given in (3.33).

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 this purpose, we use (3.5), (3.15), (3.33), and (3.35) to obtain 11 , z ∈ 2 ∩ U 1 .
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 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) are the interior and exterior of C, respectively.
In order to obtain the asymptotic behavior of the polynomials p n (λ), we use the substitution which gives the relations (1.22)-(1.24) in Theorem 1.8.

Proposition 3.4
The support of the counting measure of the zeros of the polynomials π k (z) outside an arbitrary small disk U 1 surrounding the point z = 1 tends uniformly to the curve C defined in (1.18). The zeros are within a distance o(1/k) from the curve defined by 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 lie in Int(C). The normalized 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), (3.40) 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 in (1.18) at a rate O (log k/k) (see Fig. 9). 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 weakstar limit of the zeros distribution of the polynomials π k (z) (see [44], Theorem 2.3 and [49] Chapter 3).

Post-critical Case
In this section, we assume that 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 Fig. 10 in red for two values of the parameter r . 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 Fig. 10 that the only possibility is to deform the contour to the contour C z 0 ∪ C o with C o as depicted in Fig. 11. To simplify the notation, we define the following. with φ z 0 (z) defined in (3.10) and C z 0 defined in (3.12).
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 Y (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 behavior at z = ∞: 2. Jump discontinuities (with φ(z) as in (4.1)) ( Fig. 12): 3. Endpoint behavior at z = 0 and z = 1: 4. Large z boundary behavior: 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 Fig. 13. These define new domains 0 , 1 , 2 , and ∞ . Define Then this matrix-valued function has the following jump discontinuities: where T is contour defined in Fig. 13. Fig. 13 The jump matrices for T (z) and the contour Endpoint behavior: Large z boundary behaviour:

Proposition 4.2 There exists a constant c
The proof of the proposition follows in a straightforward way by observing that by contruction (see Fig. 11), Re φ(z) > 0 for z ∈ C e \U z 0 and Re φ(z) < 0 for e −πiγ σ 3 as z ∈ (0, 1), 3. Large z boundary behavior: A particular solution is given bỹ which leads to the particular solution (4.4)

The Local Parametrix at z = z 0
The aim of this section is to construct a local parametrix P 0 (z) in a small neighborhood U z 0 of z 0 having the same jump properties as T (z) for z near z 0 and matching the outer parametric P ∞ (z) in the limit k → ∞ for z ∈ ∂U z 0 .

RH Problem for P
as k → ∞ and z ∈ ∂U z 0 . (4.5) In order to build such a 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 Fig. 14 in a neighborhood 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.

Model Problem
Consider the model problem for the 2×2 matrix function (ξ ) analytic in C\{R∪iR} with jumps and boundary behavior as ξ → ∞ , (4.10) 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 (see, e.g., Chapter 19, [1]). This equation has a nontrivial solution U(a, ξ) associated with any a ∈ C specified by the asymptotic behavior which is an entire analytic function of ξ . Three other solutions can be obtained by symmetry: U(a, −ξ), U(−a, iξ), U(−a, −iξ).
The relations among the above four solutions are (4.14) Moreover, solutions for a and a + 1 are connected by By using (4.13)-(4.15), the solution of the Riemann-Hilbert problem (4.9) and (4.10) takes the form Furthermore, from (4.12) one obtains the extra terms of the asymptotic expansion of (ξ ) as ξ → ∞:

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 neighborhood of U z 0 , 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 neighborhood of z 0 to a neighborhood of 0, and it is specified by We observe that (4.20) 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 z 0 . 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 neighborhood of U z 0 . 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 neighborhood of z 0 .
From the expression (4.21) and (4.22), the matching between P 0 (z) and P ∞ (z) takes the form 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,14,15]:P 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 From the above relation, the matrix C takes the form The improved parametrix gives the following matching betweenP ∞ (z) andP 0 (z) as k → ∞ and z 0 ∈ ∂U z 0 : which shows thatP ∞ (z)(P 0 (z)) −1 = I + O (1/k α ), α > 0, as k → ∞ when z 0 ∈ ∂U z 0 for γ ∈ [0, 1).

Riemann-Hilbert Problem for the Error Matrix R
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
1. R is analytic in C\ R , 2. For z ∈ R , we have with 3. As z → ∞, we have R(z) = I + O z −1 .
The jump matrices across the contour R \∂U z 0 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 z 0 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 z 0 . 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): (a) 0 < γ < 1 2 : v R (z) =P ∞ (z)(P 0 (z)) −1 = I + v 1 R (z) k (4.31) where v (1) By the standard theory of small norm Riemann-Hilbert problems, one has a similar expansion for R(z) in the large k limit.

R
(1) R (z) defined in (4.30). In addition R (1) (z) is analytic in C\∂U z 0 and R (1) (∞) = 0. The unique function that satisfies those conditions is given by where the integral is taken along ∂U z 0 . Using the definition of β 21 , w (z), and e(z) given in (4.17), (4.20), and (4.28), respectively, one has with the constant c defined in (4.24).
Since v (1) R (z) has a third order pole at z = z 0 , the unique function that satisfies those conditions is given by From the comments after the statement of Theorem 1.8 and Theorem 1.9, it is clear from (1.25) and (1.27) that the zeros of the polynomials p n (λ) accumulate in the limit n = ks + l → ∞, with s and l given and l = 0, . . . , s − 2, along the curveĈ defined in (1.17). In order to show that the measureν is the weak-star limit of the zero density ν n , we will show that uniformly for λ in compact subsets of unbounded component of C\Ĉ. Furthermore, C is the boundary of its polynomial convex hull. Then it follows that the measureν is the weak-star limit of the zero counting measure ν n of the polynomials p n (λ) for n = kd + l, l = 0, . . . , s − 2 (see [44] and [49] Chapter 3). The proof of Theorem 1.4 is then completed.