A matrix model with a singular weight and Painleve III

The purpose of this article is to study the asymptotics as N → ∞ of the partition function where z ∈ R\{0} and 0 ≤ t < ∞; the particular case z = t = 0 corresponds to the partition function of the Gaussian Unitary Ensemble (GUE). This multiple integral belongs to a general class where w(x) = exp -V (x) is the weighting function and V (x) is known as the potential. These integrals have been the subject of extensive investigations in Random Matrix Theory (RMT), because they contain all the information on the correlations of the eigenvalues and are the starting point to study their linear statistics as well as global fluctuations of the spectra. The potential of the GUE is V (x) = -x 2 /2; therefore, E N (z, t) can be thought of as the partition function of a matrix model obtained by perturbing the GUE potential with a first and a second order pole, with t and z measuring the strength of the perturbations. Then, a natural question arises: what happens as the average (1.1) approaches the GUE partition function? The main result of this paper is that as N → ∞ while N 1 2 t and N 1 2 z converge to finite constants, a phase transition emerges where the asymptotics of E N (z, t) is characterized by a solution of the Painlevé III (PIII) equation. In most applications the potential V (x) is required to have some regularity properties; for example, imposing that V (x) should be real analytic together with appropriate boundary conditions suffices in many cases, as it guarantees that the limiting mean density of the eigenvalues, known as equilibrium measure, is supported on a finite union of intervals. However, recently matrix models whose weight w(x) has an essential singularity have appeared in several area of mathematics and physics, like number theory, quantum transport and finite-temperature field theory (see, e.g., [4,9,10,30,34]). For such singular potentials the asymptotic analysis of the partition function, or of any statistics of the spectra, becomes substantially involved and no studies of the double scaling limits are available. The eigenvalues of a matrix ensemble with partition function (1.2) form a determinantal point process. It is remarkable that with an appropriate choice of the scaling limit, as N → ∞ the kernel of this process becomes universal and depends only on the local properties of the equilibrium measure (see, e.g., [15] and references therein). Finding the universal kernel can be reduced to the asymptotics analysis in the same scaling limit of the polynomials orthogonal with respect to w(x). By standard theory of orthogonal polynomials [35,Chap. 2.2] and by the results in [5], the asymptotics of E N (z, t) can expressed in terms of the same system of orthogonal polynomials. In order for the integral (1.1) to converge either z = 0 or z = t = 0, when it is simply the partition function of the GUE whose equilibrium measure is the semicircle law, which is supported in an interval symmetric with respect to the origin. In the limit N → ∞ the parameter t does not contribute to the equilibrium measure because with the correct scaling it becomes asymptotically negligible. Thus, the limiting values z = 0 is bifurcation point: away from it, the second order pole in the potential splits the support of the equilibrium measure in two intervals symmetric with respect to the origin [33]. Previous studies of double scaling limits in matrix models have concerned various types of critical points of the equilibrium measure. When it vanishes quadratically inside the support, the universality of the spectra correlations is characterized in the double scaling limit by the second Painlevé equation [7,8,12,13]. In this case too, two intervals in the support of the equilibrium measure coalesce into one at the critical point; however, the main difference from the model of the present paper is that the potential is real analytic and has not any singularity at the critical point. This changes the nature of the problem entirely. When the equilibrium measure vanishes as the power 5/2 at an endpoint of its support, then in the double scaling limit universality is identified by a fourth order analogue of Painlevé one [14]. This is the first time that a critical phenomenon associated to an essential singularity in the weighting function has been studied. It seems that the characterization of the double scaling limit by PIII is specific to a pole emerging in the background of a smooth potential and appears to be a new universal feature of the spectra of unitary matrix ensembles. This property is likely to be shared by other matrix models whose weighting function has the same type of singularity. The average (1.1) was introduced by Berry and Shukla [4] in their study of the random function where Λ N (x) := N j=1 (x-x j ) and x 1 , . . . , x N are a set of random variables -in general not independent. The value distribution P (Q N ) is important mainly for two reasons: firstly, it is a sensitive indicator of the degree of the repulsion of two neighbouring x j 's, in the sense that the rate of decay of P (Q N ) is a measure of the rigidity of x 1 , . . . , x N ; secondly, if Λ N (x) is replaced by the Riemann zeta function (or more precisely the Hardy function, which is real on the critical line), then the Riemann hypothesis implies that Q N (x) > 0. Therefore, P (Q N ) provides valuable information on the statistics of the zeros of the Riemann zeta function. The connection between RMT and the theory of the Riemann zeta function suggests to replace Λ N (x) with the characteristic polynomial of a random matrix from the GUE. In this case an explicit formula for the probability density P (Q N ) is very difficult to find. However, all the information on P (Q N ) is contained in the integral (1.1). It is straightforward to see from the definition that E N (z, t) is a real analytic function of z and t and is even. This symmetry suggests the introduction of the power expansion which for z = 0 converges uniformly in t in any closed subset of R + including the origin. Therefore, E N (z, t) can be interpreted as the generating function of the coefficients E N 2m (z). Berry and Shukla [4] showed that the moments of P (Q N ) are given by In a previous article [33] we showed that in the range where B N is the ensemble average These limits, however, are not uniform in z: in the regime where zN 1 2 remains finite they are determined by PIII and change drastically. Since the integral (1.5) extends to zero, the fact that the asymptotics are not uniform cannot be ignored when computing the moments of P (Q N ). As already mentioned, recently ensembles whose weighting function have essential singularities have appeared in several applications. Usually, such models depend on external parameters beside the matrix dimension N ; therefore, double scaling limits play a fundamental role in their asymptotic behaviour. Since the universal properties of the spectrum and the associated Painlevé transcendents are identified by the singular points of the equilibrium measure, which in turn are determined by the analytical properties of V (x), we would expect the same critical behaviour whenever V (x) has the same poles as the potential appearing in E N (z, t). Therefore, the results in this paper have implications for other matrix models, beside its direct application to the theory of the Riemann zeta function suggested by Berry and Shukla [4], and provides a pathway to tackle similar asymptotic problems. As an example consider the weight introduced by Chen and Its [10], This is a singular perturbation of the Laguerre Unitary Ensemble. For the unperturbed weight, i.e. when s = 0, Forrester and Witte [24,25] discovered that the generating function of the probability that an interval contains k eigenvalues in the hard edge scaling limit can be evaluated in terms of a PIII transcendent in σ-form. Chen and Its [10] studied the polynomials orthogonal with respect to the weight (1.9) and showed that for finite N the partition function can be written as an integral involving PIII. In unpublished work [11] 1 they also investigated the asymptotics of these orthogonal polynomials. When t = 0 in Eq. (1.1) the system of monic polynomials orthogonal with respect to the weight of the partition function E N (z, 0) can be mapped into that orthogonal with respect to W ± 1 2 (x) by a change of variables -the respective partition functions, however, would still be different. The pole of order one in the exponent of the weighting function of E N (z, t) does not contribute to the equilibrium measure, as it is asymptotically negligible [33]. Therefore, in the general setting we would expect that the matrix model with weighting W α (x) should manifest the same critical behaviour discovered in this article. There exist other matrix models whose weighting function has essential singularities and whose understanding could shed light on important unsolved problems. By setting α = 3n and s = n, the partition function associated to W 3n (x) becomes the moment generating function of the probability density of the Wigner delay time [34]. This is the average time that an electron spends when scattered by an open cavity and plays a fundamental role in the theory of mesoscopic quantum dots. The distribution of the Wigner delay time is far from being understood [34,36]. An other example of a similar problem is provided by the distribution of the roots of the derivative of the characteristic polynomials of a random unitary matrix. It has an integral representation that reduces to a matrix average over the unitary group; the weighting function of this average has essential singularities analogous to those in the partition function E N (z, t) [17,31]. This distribution is very elusive; the calculation of more refined formulae than those known at the moment would improve present results on the percentage of the number of zeros of the Riemann zeta function on the critical line. The asymptotics of the system of polynomials orthogonal with respect to the weighting function exp (-V z,t (x)), where is characterized by two distinct regimes. It is well known that such asymptotics is determined, after appropriate rescaling, by the behaviour of the equilibrium measure, which is the solution of a particular variational problem. If z = t = 0, then the equilibrium measure is the semicircle law, which is supported in [-2, 2] and is defined by (1.11) In the limit N → ∞ the simple pole in the potential (1.10) does not contribute to the equilibrium measure, which depends only on z; in [33] we showed that in the interval , where c 1 and c 2 are two positive constants, it is supported on two disjoint intervals symmetric with respect to the origin. As z → 0 the quadratic singularity vanishes and the gap between the intervals closes. If this migration is fast enough, i.e. when zN 1/2 remains finite as N → ∞, formulae (1.6) and (1.7) cease to be valid and a phase transition emerges. The main result of this paper is that in this region of the phase space (z, N ) the asymptotics of E N (z, t) is expressed in terms of a solution of the PIII equation. As we shall see in Sec. 3, for technical reasons we will use the equilibrium measure (1.11) and then deform the relevant Riemann-Hilbert problem (RHP) by incorporating the essential singularities in the jump matrix. Let us introduce the scaling where u 1 ≥ 0 and u 2 > 0 are finite. The rate of convergence of (u 1,N , u 2,N ) is chosen so that it does not affect the leading order asymptotics of E N (z, t) (see Remark 2.1). Consider the rescaled weight where y is rescaled as in (1.11). Let π j (y) denote the monic polynomials orthogonal with respect to w N (y), i.e. w N (y)π j (y)π k (y)dy = h j δ jk , j, k = 1, 2, . . . , (1.15) where the subscript in π j (y) refers to the degree of the polynomial. This change of variables turns the integral (1.1) into where where w N (y)y j dy. j = 0, 1, . . . Our approach to the analysis of the integral (1.1) is based on the steepest descent method to compute the asymptotics of orthogonal polynomials introduced by Deift and Zhou [20] and further developed by Deift et al. [18,19] (see also [7] in connection to the double scaling limit). The starting point of this technique is the characterization of the orthogonal polynomials in terms of the solution of a RHP due to Fokas et al. [22,23]. Define the matrix valued function where κ N -1 = -2πi/h N -1 . It solves the following RHP: where Y + (y) and Y -(y) denotes the limiting value of Y (y) as it approaches the left-and right-hand side of the real axis. The jump matrix in the second condition is continuous, as the weight (1.14) can be uniformly bounded on R; therefore, there is no need to specify any special behaviour of Y (y) near the origin in the definition (1.21). Bertola et al. [5] showed that the logarithmic derivatives of G N (u 1,N , u 2,N ) admit integral representations involving Y (y), which in our context can be phrased as follows [33]: Lemma 1.1 (Bertola, Eynard and Harnad [5]). The following differential identities hold: where the contour of integration is a small loop around y = 0 oriented counter-clockwise and σ 3 = ( 1 0 0 -1 ). Note that, although the function Y (y) has a jump discontinuity on the real axis, it has a uniform asymptotic expansion near the origin, which is used to compute the residues in the right hand sides of (1.22). Thus, the differential identities (1.22) allow us to compute the asymptotics of G N (u 1,N , u 2,N ) in terms of that of Y (y). A standard technique in the Deift-Zhou steepest descent is the opening of the lens. The g-function is used to construct a sequence of transformations that turn (1.21) into a RHP which is amenable to an asymptotic analysis. As N → ∞ the solution of the modified RHP outside small neighbours of the critical points of the equilibrium measure can be approximated in the same way as in [15,Chap. 7.3] and [18]. In our problem such critical points are the edges of support of the equilibrium measure and the origin, where the essential singularity appears. These aspects of the Riemann-Hilbert analysis are discussed in Sec 3. The challenge for the weight (1.14) comes from the essential singularity. Inside a small disc around the origin the asymptotic approximation of the RHP breaks down; thus, it is replaced by an exactly solvable "model" RHP, whose solution Φ, known as the local parametrix, matches the asymptotic solution outside such a disc. A widely used technique in investigations of double scaling limits in RMT (see, e.g., [8,[12][13][14]) is to express Φ in terms of a special solution of a nonlinear differential equation, which in turn gives the compatibility conditions for the following system of linear ODEs: where ζ ∈ C and u is a parameter measuring the strength of the perturbation from the critical point. The matrices A(ζ, u) and B(ζ, u) are rational functions of ζ. In our case the parameter u is replaced by the vector (u 1 , u 2 ). The compatibility conditions for (1.23) are among isomonodromic deformations which include the Painlevé equations. In previous studies of double scaling limits in one matrix models, one parameter was sufficient to represent the whole family of perturbations from the critical point. In the RHP that we study, however, the dependence of the local parametrix on two perturbation parameters instead of one increases the technical difficulties considerably. Indeed, the compatibility conditions for the linear system analogous to (1.23) give us a nonlinear PDE and not an ODE. In Sec. 4 we use the Hamiltonian theory of isomonodromic deformations [29] to express such a PDE in terms of a time-dependent Hamiltonian system ODEs in the variables u 1 and u 2 . This simplifies the problem substantially from the computational point of view. The differential identities (1.22) give a link between the partition function E N (z, t) and the RHP (1.21); in Secs. 5 and 6 we compute their asymptotics in terms of the solution of the Hamiltonian system. This completes the proof of Theorem 2.1. The boundary conditions of the system of ODEs are provided by the asymptotics of the local parametrix as u 2 → 0, which we compute in Sec. 7. Sec. 8 concerns the reduction of the Hamiltonian system to the PIII equation. In order to achieve this goal we need to make two observations. Firstly, when t = 0, and therefore by Eqs. (1.12) and (1.13) u 1 = 0 too, the symmetries of the partition function E N (t, z) and the initial conditions computed in Sec. 7 simplify the system of ODEs, reducing it to a nonlinear second order ODE. The Taylor coefficients in Eq. (1.4), and in turn the partition function (1.1), can be expressed in terms of the solution of this nonlinear ODE. Secondly, Chen and Its [10] showed that the solution of RHP for the orthogonal polynomials associated to the weighting function (1.9) is described by the PIII equation. When α = ±1/2 in the weight (1.9) and t = 0 in the potential (1.10), the two respective systems of orthogonal polynomials are related by a change of variables. Thus, the corresponding RHPs can be mapped into each other. These two facts and a lengthy calculation allow us to identify the second order ODE that we discovered with a scaling limit of the PIII equation in [10]. The ODE associated to the problem in this paper belongs to the PIII family too, but it is different, as the parameters that define it are not those of the ODE in [10]. The differential identities (1.22) are the starting point to compute an asymptotic formula for the partition function E N (z, t). When the perturbation parameters z and t are O(N -1/2 ), Eqs. (1.22a) and (1.22b) can be expressed at leading order in terms of a particular solution of non-linear Hamiltonian equations, which are isomonodromic deformations of the compatibility conditions for the ODEs (1.23). The phase space of this Hamiltonian system is four-dimensional and the time variables are the parameters (u 1 , u 2 ) defined in Eq. (1.13). Denote by Z(u 1 , u 2 ) = (P 1 , Q 1 , P 2 , Q 2 ) such a solution. Then, Z(u 1 , u 2 ) is determined uniquely by the monodromy data of the local parametrix at the origin of the RHP (1.21), which in turn gives the behaviour of Z(u 1 , u 2 ) as u 2 → 0, thereby providing the initial conditions. The asymptotic formulae of Z(u 1 , u 2 ) in the limit as u 2 → 0 are rather involved and are given by Theorem 7.1. Theorem 2.1. Let u 1 and u 2 be the parameters defined in Eq. (1.13). At leading order as ) with Hamiltonian functions ) as well as Remark 2.1. The definitions of (u 1,N , u 2,N ) and of their limits (u 1 , u 2 ) (see Eqs. (1.12) and (1.13)) imply that the error terms in Eqs. (2.1a) and (2.1b) are not affected by the rate of convergence of (u 1,N , u 2,N ). All the information on the partition function E N (z, t) is contained in the coefficients of the generating function (1.4), which can be studied by looking at the projection at u 1 = 0. The expressions of the initial conditions of this subsystem simplifies considerably compared to the formulae in Theorem 7.1. Corollary 2.1. Let u 1 = 0. As u 2 → 0 we have 2 , (2.5a) ) 2 . (2.5d) The subsystem (2.4) is equivalent to a fourth order ODE. However, from the integral (1.17) we see that the average G N (u 1,N , u 2,N ) is an even function in u 1,N , and hence the derivatives ∂ 2k+1 u 1,N log G N are zero when u 1,N = 0. Then, Eq. (2.1a) gives The first two equations in (2.6) provide two possible sets of relations between the variables P 1 , Q 1 , P 2 and Q 2 , namely or The condition ∂ 4 u 1 H 1 u 1 =0 = 0 does not yield any new constraints; when k > 2 Eqs. (2.6) involve high powers of P 1 and Q 1 and cannot be solved analytically. Equation (2.7a) is compatible with the behaviour of the coordinates P 1 and Q 1 as u 2 → 0 in Corollary 2.1, while Eq. (2.7b) is not. Therefore, at u 1 = 0, the correct algebraic relations among the canonical coordinates are (2.7a). For convenience let us make the change of variable r := √ u 2 . Then, Eqs. (2.7a) reduce the Hamilton equations (2.2) in Theorems 2.1 to the system ) where The following theorem is the main result of this article. Theorem 2.2. Consider the PIII equation with initial conditions when N is even, and when N is odd. Then, the trajectory ) solves the system of ODEs (2.8), where the solution v Y (r) of (2.9) is specified uniquely by the monodromy data of the local parametrix at the origin of the RHP (1.21). Finally, there exist two closed intervals where and H(r) is computed along the trajectory (2.12). Remark 2.2. Recall that PIII is the following second order nonlinear ODE Hence, Eq. (2.9) is a special case of (2.14) with parameters a = (-1) ) is an odd function of u 1 , the series in the right-hand side of Eq. (2.13) can be expressed as as it is also apparent from Eq. (2.1a). The series expansion, however, provides a direct way of evaluating the integral in the right-hand side of (2.15), as the derivatives can be computed explicitly by combining the Hamilton equations (2.2) with the algebraic relations (2.7a). For example, we have In this section we apply the Deift-Zhou steepest descent analysis to the RHP (1.21). Define the g-function where dµ(q) is the equilibrium measure for the potential V (y) = y 2 /2, defined in Eq. (1.11). for some constant l. Define where Y (y) is given in Eq. (1.20). Then, T (y) solves the RHP The jump matrix J T (y) is given by J T (y) := e -N (g + (y)-g -(y)) e -N (g + (y)+g -(y))+2F (y) 0 e N (g + (y)-g -(y)) , y ∈ R. Note that since the essential singularities can be controlled on the real line, the jump matrix J T (y) is continuous on R. We now perform a standard technique in the steepest decent method (see, e.g., [7,18,19]). Take a point y 0 sufficiently close to the origin and define the lens contours as in Fig. 1. The precise location of y 0 will be specified in Sec. 4.1, where we study the local parametrix near zero. The boundaries of the lens regions L ±,j are ∂L ±,j = Ξ ±,j ∪ Ξ R,j and we denote by Ξ the union of the lens contours in Fig. 1. Write T (y), y outside the lenses, Then, S(y) satisfies the RHP The jump matrix J S (y) is defined piecewise: Equation (3.6) combined with the inequalities (3.2) imply that outside of some discs D ±2 and D 0 of sufficiently small radius δ centred at the points ±2 and zero, J S (y) can be approximated by for y ∈ [-2, 2] and the identity on the rest of the contour. This suggests the following approximation to S(y) outside of D 0 and D ±2 : The outer parametrix that solves (3.8) can be constructed as in [15,Chap. 7.6] and [18]: where the branch cut is chosen to be on [-2, 2] and γ ∼ 1 as y → ∞. The approximation of S(y) by S ∞ (y) fails near the points ±2. Therefore, we must find exact solutions to the RHP for S(y) and match them with S ∞ (y) as y moves away from the edges of the interval [-2, 2]. These local parametrices can be constructed using Airy functions as in [15,Chap. 7.6] and in [7,18,19] and we shall not repeat the derivation here. The essential singularity in the weight (1.14) means that the approximation of S(y) by S ∞ (y) breaks down near the origin too, as becomes apparent from Eq. (3.7). Therefore, we must solve the RHP for S(y) in a neighbourhood of the origin exactly and match it with S ∞ (y) away from the singularity. Take a small neighborhood of the origin D 0 containing the interval (-y 0 , y 0 ). We want to solve the following RHP inside D 0 : where the jump matrix is , y ∈ (-y 0 , y 0 ). In order to study the model problem (4.1), we introduce a conformal map ζ from D 0 into an open neighbourhood of the origin that maps the lens contours Ξ as in Fig. 2 and sends Note that g+ (y) ∈ iR for y ∈ R ∩ D 0 . For small but finite δ we denote by ζ the analytic continuation of (4.3) to the complex plane for all y ∈ D 0 . In the limit N → ∞, ζ maps the disc D 0 to the whole complex plane. Equation (4.3) implies that, in order for the images ζ(±y 0 ) to remain at a finite distance from the essential singularity as N → ∞, we need to choose y 0 = O(1/N ). For small enough δ, U 1 (y) and U 2 (y) are conformal inside D 0 and Let Γ be the union of the contours in the right-hand side of Fig. 2. As N → ∞ let us introduce the following RHP in the ζ-plane: The jump matrices are (4.7) The local parametrix that solves the RHP (4.1) can be obtained from the solution of (4.6) through the relations S (0) (y) = S ∞ (y)e -N g+ (0)σ 3 0 -1 1 0 P (ζ)e N g+ (0)σ 3 , Im(y) > 0, (4.8a) with N g+ (0) ∈ iR. Before proceeding we need to prove the following existence theorem. Theorem 4.1. The solution to the RHP (4.6) exists and is unique. Then, P (ζ) solves the RHP where the jump matrices J P (ζ) are given by (4.10) Let J P ,k (ζ) denote the jump matrix on Γ k . In order to prove the existence and uniqueness of P (ζ) we refer to the theory developed by Zhou [38], which for our purposes is equivalent to the following statements: a. the RHP (4.9) and the jump matrices (4.10) do not have singularities at the selfintersection points ±y 0 of the contours; b. the jump matrices J P ,k (ζ) satisfy the cyclic conditions J -1 P ,2 J P ,6 J P ,0 J P ,1 = I and J -1 P ,2 J P ,5 J P ,4 J P ,3 = I at the points ±y 0 ; c. there is no non-trivial function P0 (ζ) that satisfies properties 1. and 2. of the RHP (4.9) and behaves as The conditions a. and b. can be easily verified by direct inspection. In order to prove c., suppose that such a P0 (ζ) exists and denote the region between Γ 1 (Γ 3 ) and Γ 0 (Γ 4 ) by Ω 0,+ (Ω 4,+ ) and the region between Γ 0 (Γ 4 ) and Γ 6 (Γ 5 ) by Ω 0,-(Ω 4,-) (see Fig. 3). Then, define The matrix X(ζ) has the following properties [18]: Let C be a close contour in the upper-half plane consisting of the real line and a large semicircle S. Consider the integral The second condition in Eq. (4.11) gives The essential singularity in the jump matrix does not affect this integral, since the exponential can be uniformly bounded in a small interval on the real line containing the origin; therefore, Cauchy's theorem implies that integral (4.13) is zero. By adding this matrix to its Hermitian conjugate, we see that This implies that the first column of X -(ζ) is identically zero. From the jump conditions in (4.11), it follows that the second column of X + (ζ) is identically zero too. Now, let the first column of X be X 1 and the second column of X be X 2 and define the vector From Eq. (4.11) we see that f (ζ) is analytic in the whole complex plane and behaves as ζ -1 as ζ → ∞. Hence, by Liouville's theorem, we have f (ζ) = 0. This shows that there is no non-trivial function P0 that satisfies 1. and 2. in (4.9) and such that P0 = O(ζ -1 ) as ζ → ∞. Therefore, the RHP (4.6) is uniquely solvable. We will now transform the RHP (4.6) into one with constant jumps but with essential singularities at ζ = 0 and ζ = ∞. Then, theory of isomonodromy deformations developed by Jimbo, Miwa, and Ueno [29] can be applied to derive differential equations that give the solution to (4.6). Let us deform the jump contours of (4.6) as in Fig. 4 and let Γ ± be the semicircles in the uppe/lower half planes. Note that ∂Ω 2,± = Γ ± ∪ Γ 2 . Then, define the function By direct inspection we can see that the matrix Φ(ζ) satisfies the RHP where the jump matrices J Φ (ζ) are given by (4.16) The function satisfies a RHP with constant jumps on Γ + ∪ Γ -∪ Γ 2 : Furthermore, it has essential singularities at the points ζ = 0 and ζ = ∞: where Ψ0 (ζ) and Ψ∞ (ζ) have the asymptotic expansions Note that the determinants of Ψ0 and Ψ∞ are both equal to one. From (4.4) it follows that U j (y), j = 1, 2, behaves as where u j,N and u j were defined in (1.12) and (1.13), respectively. Let Φ c (ζ) and Ψ c (ζ) denote the functions (4.14) and (4.17) with U j (y) replaced by u j . Since the jump matrices J Φ (ζ) are all bounded, the matrix Φ(ζ)Φ -1 c (ζ) will have jump discontinuities on the contours Γ 2 and Γ ± of order I + O(N -2 ) and behaves as I + O(ζ -1 ) as ζ → ∞. By the general theory of the RHP (see e.g., [18,Sec. 7] and also Eq. (5.4)), we have Let us consider the RHP with jumps (4.18) and asymptotic behaviour (4.19). Since the matrix Ψ(ζ) has constant jumps on the complex plane, the matrices are rational functions with poles at ζ = 0 and ζ = ∞ only. From (4. 19) we see that their explicit expressions are where (•) pp denotes the singular part at ζ = 0. We can write the matrix A as follows: The coefficient of ζ -1 can be extracted by expanding (4.22 From (4.23) we see that the coefficients of ζ -6 , ζ -5 and ζ -4 in the determinant of A(ζ) coincide with those in det 1 2 This implies Finally, the compatibility of the linear differential systems (4.22) implies where ∂ j is the derivative with respect to u j and ∂ ζ the derivative with respect to ζ. The compatibility conditions (4.26) constitute a system of PDEs in the variables u 1 and u 2 . The theory of isomonodromic deformations of linear ODEs [29] states that such PDEs are equivalent to an integrable Hamiltonian system of ODEs, which can be studied much more easily. The main idea is to express the terms involving the commutators in Eqs. as well as In addition, require that the coefficients of the expansion (4.27) should be of the form Note that the A j 's are traceless; therefore, they belong to the Lie algebra sl 2 (C). We shall see that the deformations Ȧ = [A, B] (4.30) define a Hamiltonian flow. The ODEs (4.30) do not identify the pair of matrices A and B uniquely. If we conjugate both A and B by a constant matrix G ∈ sl 2 (C), which is either diagonal or a scalar multiple of A 1 , then the form of the coefficients (4.29) and the trajectories of the flow (4.30) remain invariant. Thus, rather than the set of matrices (4.27), we look at their conjugacy classes. This imposes two conditions on the matrix entries of the coefficients A j and reduces the independent parameters from eight to six. In the particular example (4.23a), these extra relations are equivalent to Eqs. (4.25). For simplicity, and without loss of generality, we shall not distinguish between the set of matrices defined in (4.27) and their equivalent classes. To begin with we need to equip this six-dimensional manifold with a Poisson structure. Let the loop algebra g be the set of smooth maps f : S 1 → sl 2 (C), where S 1 is the unit circle. Then, split g into the subalgebras g + and g -, where g + (g -) is the set of maps that admits holomorphic extension to the inside (outside) of the unit circle; the maps in g -should also vanish at infinity. Using the pairing we can identify g -with the dual of g + . This structure allows us to define Poisson brackets on g -(see, e.g., [21]) where f and g are functions on g -= g * + . Using the method of the moment map [2,3,26,28], this Poisson structure can be restricted to the space of matrices of the form (4.27). The Poisson brackets between the matrix entries a j , b j and c j are degenerate and are given by All the other brackets vanish. Equations (4.32) define a Poisson manifold. Remark 4.1. Although, strictly speaking, the matrices (4.27) do not belong to the algebra g -, because of the term iσ 3 , one can think of the space of matrices (4.27) as being parameterised by the entries a j , b j and c j and use Eqs. (4.32) to define the Poisson brackets. In fact, the original construction in [2,3] applies to much more general matrices. The interested reader is invited to consult these references. The symplectic leaves in this Poisson manifold are the co-adjoint orbits of the matrix A, specified by det(A) = const. A set of canonical coordinates on these symplectic leaves can be found as follows [2,37]. Let ψ(ζ) be an eigenvector of A(ζ) with eigenvalue λ(ζ). Then, the poles q j of ψ(ζ) and the values p j = λ(q j ) of the eigenvalue λ(ζ) at q j satisfy the relations {p j , q k } = δ jk , {p j , p k } = {q j , q k } = 0. For example, let ψ(ζ) be an eigenvector of A(ζ) normalized by The poles q j of ψ(ζ) are the zeros of the polynomial and the eigenvalues at these points are Therefore, the symplectic form on these symplectic leaves are given by Hence, are a set of canonical coordinates on the symplectic leaves and the parameters where v 1 and v 2 are functions of the matrix elements a j , b j and c j , constitute a set of coordinates on the Poisson manifold. The equivalence between the deformations (4.30) and Hamiltonian systems is provided by the following proposition, which is a particular case of the flows studied in [1][2][3]. whose Hamiltonians are Equations (4.33) and (4.25) allow us to express the matrix entries a j , b j and c j in terms of the canonical variables P 1 , Q 1 , P 2 , Q 2 and the parameters u 1 , u 2 , which together provide a set of coordinates for the Poisson manifold. In turn, the Hamiltonians H 1 and H 2 can be written in terms of these coordinates. The outcomes of this calculation are formulae (2.3a) and (2.3b). The terms ∂ ζ B j in Eq. (4.26a) imply that deformations (4.34), and hence the Hamiltonian flows generated by H 1 and H 2 , are not sufficient to describe the compatibility conditions (4.26a). However, using a technique developed by Mazzocco and Mo [32], we can show that Eq. (4.26a) is equivalent to a time-dependent Hamiltonian system. We are now in a position to prove the following. Theorem 4.2. Let H j and h j , j = 1, 2 be the Hamiltonians given in Theorem 2.1, Eqs. (2.3). The deformations of the matrix A(ζ) defined by are equivalent to the time-dependent Hamiltonian equations Proof. Let ν 1 , . . . , ν 4 , u 1 , u 2 be a new set of coordinates on the Poisson manifold and denote by ∂ ν 1 and ∂ ν 2 the partial derivatives that keep the variables ν 1 , . . . , ν 4 fixed. We can choose these coordinates so that they are related to the matrix entries a j , b j and c j by the equations Therefore, we can write Eq. (4.36) as time-dependent Hamiltonian equations: In order the find the evolution equations in the form for suitable Hamiltonians H j , we must replace ∂ ν j by the partial differentiation ∂ can j that keeps the canonical coordinates We have From (4.38) and (4.33), we can compute the derivatives ∂ ν j P k and ∂ ν j Q k : Combining these formulae with we arrive at and H j = H j + h j . Finally, Eq. ( 4.40) becomes Remark 4.2. The Hamilton equations (4.37) are ODEs in the variables u 1 and u 2 . A special solution of these ODEs, specified by the jump conditions (4.18), will give the local parametrix via (4.22), (4.17), (4.14) and (4.8). In particular, we will show that the Hamiltonians H 1 and H 2 provide the leading order terms of the differential identities (1.22). Here we prove that the outer and local parametrices that we constructed in Secs. 3.2 and 4, respectively, are indeed good approximations to the solution S(y) of the RHP (3.5). where S (p) (y) are the local parametrices inside the neighborhoods D ±2 and D 0 . The function R(y) has jump discontinuities on the contour Γ R shown in Fig. 5. In particular, R(y) From the definition of R(y), it follows that the order of magnitude of the jump matrix where Γ 0 is the subset of [-2, 2] that lies outside of D ±2 and D 0 . The general theory of the RHP (see, e.g. [18,Sec. 7], and also Eq. (4.21)) asserts that the third statement of RHP (5.2) and Eq. ( 5.1) lead to uniformly in C. Therefore, at leading order as N → ∞ the solution of the RHP (3.5) is given by Combining these expressions with equations (3.4) and (3.3) we obtain an asymptotic formula for the solution of the original RHP (1.21). In this section we complete the proof of Theorem 2.1 by computing the leading order terms of the differential identities (1.22). We shall see that their asymptotics are determined by the behaviour of the local parametrix near the origin. In turn, such formulae can be identified with the Hamiltonians (2.3a) and (2.3b). Lemma 6.1. Let G N (u 1,N , u 2,N ) be the partition function (1.17) and Ψ0 (ζ) the matrix introduced in (4.19a). At leading order as N → ∞ we have In addition, from Eq. (5.1) near the origin we have S(y) = R(y)S (0) (y). Then, Eqs. (4.8) give Now, from Eqs. (4.14) and (4.17) we have Therefore, P (ζ) is related to Ψ0 (ζ) by Define By combining Eqs. (6.3) and (6.5) we can write The final step consists in inserting this expression into into the differential identities (1.22). Note that the definitions of the jump matrices of the RHPs (3.5) and (4.1) combined with Eq. (5.1) imply that the problem (5.2) does not have any oscillatory factors. Therefore, we can safely differentiate the right-hand side of Eq. (5.4) and see that R (y) = O(1/N ). Finally, the conformal map (4.3) leads to (6.1). The leading order terms of the identities (1.22) are the derivatives of the τ -function for the equations (4.26) (see, e.g., [29]). In unpublished work, Bertola, Harnad, Hurtubise and Putsai [6,27] have shown that they are given by the Hamiltonians (4.35). For completeness we repeat the proof of this statement for the system that we study. Proposition 6.1. The Hamiltonians (4.35) can be expressed as series in √ u 2 , whose coefficients can be computed using the technique in [18,Sec. 7.2]. The first four terms of this asymptotic expansion provide the initial conditions. Before we discuss the proof of Theorem 7.1 we need to introduce preliminary definitions and results. Throughout this section we shall assume that the parameters u 1 and u 2 scale as u 1 = ũ√ u 2 , where ũ remains finite. Let and consider the coefficients of the expansion of φ(x) as x → ∞ in the upper/lower half planes, namely Here E + (E -) consists of the union of the two intervals (-∞, -) and ( , ∞), oriented from -∞ to ∞, together with the semi circle of radius around the origin in the upper (lower) half plane. The exact value of is not important, as the integrand is analytic away from zero; therefore, by Cauchy's theorem, integrals with different values of will yield the same result. The residue theorem leads to the identity This formula gives us a set of relations between the φ +,j 's and φ -,j 's. For example, for the first few coefficients we have By integrating by parts, we can express the coefficients φ ±,j in terms of φ ±,0 . The integrals (7.2) are rather cumbersome to compute directly. The following lemma provides an efficient tool to this purpose. Lemma 7.1. The coefficients φ ±,j satisfy the recurrence relation with initial conditions 2 dq, (7.6a) Proof. The proof of Eq. (7.5) is an immediate consequence of integration by parts: When j = 1 this formula gives To compute φ ±,0 , we expand the factor e ũq in the integrand in (7.2) in a power series: To evaluate the integral E ± q -1 e -q 2 2 dq, note that if we change the integration variable from s to -s, then the path of integration will change from E ± to -E ∓ . Therefore, we have Now, the residue theorem gives This equation and Eq. (7.9) imply that i 2π E ± q -1 e -q 2 2 dq = ± 1 2 . Inserting this value into (7.8) and combing it with (7.6) completes the proof. We are now in a position to state the main result of this section. 2 and suppose ũ exists and is finite as u 2 → 0. Let the φ -,j 's be the coefficients introduced in Eq. (7.2). Furthermore, define and denote by φ 1 be the coefficient of x in the expansion of φ(x) in a neighbourhood of x = 0. Then, the solution Φ(ζ) of the RHP (4.15) determines the following initial conditions for the Hamilton equations in Theorem 2.1: 2 . Remark 7.1. These formulae simplify considerably at ũ1 = 0 and become the initial conditions in Corollary 2.1. Divide the ζ-plane in the three regions Ω 2,± and Ω = C \ (Ω 2,+ ∪ Ω 2,-) as in Fig. 6. Recall that ∂Ω 2,± = Γ ± ∪ Γ 2 (see Fig. 4). Define Let ∆ 0 be a small disc around the origin of fixed radius and let ξ = u From Eqs. (7.4), (7.11) and the property of the Cauchy transform, one can check that Φ (p) (ζ) satisfies the jump condition Hence, the matrix where Γ R Φ is the contour in Fig. 7. The jump matrices are Indeed, a series expansion in √ u 2 of the matrix R Φ (ζ) can be computed using the method in [18,Sec. 7.2]. In order to find the necessary initial conditions to solve the ODEs in Theorem 2.1, we will need the expansion of R Φ (ζ) up to order O u Then, the matrix R(ζ) satisfies the RHP (7.17) It follows from Eq. ( 7.3), that the jump matrix J R(ζ ) is given by These relations and the requirement that Rk = O(ζ -1 ) as ζ → ∞ define RHPs for the Rk 's that we can solve. From Eqs. (7.14) and (7.4), we see that J R,1 is given by Let us write By using (7.4) and (7.19), one can verify that the solution to the RHP for R1 is Let us now denote the jump matrix of the RHP for R2 by Ĵ R,2 := R-,1 J R,1 +J R,2 . Then, we have . This expression gives Proceeding in the same way, we can compute the matrix elements of the jump matrix We obtain where α 0 , α ±,1 and ι ±,0 are given by Then, introduce the matrix (c.f. (1.20)) where κ Using standard arguments in the theory of integrable systems, similar to those that we outlined is Sec. 4.3, Chen and Its [10,Sec. 5] showed that the Lax pair associated to X(x) is that for the PIII equation. More precisely, in our notation their result is the following. Theorem 8.1 (Chen and Its [10]). Define the function Then, v X (ω) satisfies the PIII equation When u 1,N = 0 the system of monic polynomials orthogonal with respect to the weight (1.14) can be mapped to the polynomials P n (x) with α = ±1/2 by a change of variables. In what follows, we shall use this observation to relate X(x) to the solution Y (y) of the RHP (1.21). Upon taking the double scaling limit, the matrix Y (y) can be expressed in terms of a solution of the Hamiltonian system (2.8). ) when N is even, and if N is odd. Then, the trajectory solve the Hamilton equations (2.8). This is the main result of this section and completes the proof of Theorem 2.2. Proposition 8.1 is proved in Secs. 8.1 and 8.2. Note that, although both ODEs in Theorem 8.1 and Proposition 8.1 belong to the family of PIII equations, they are different in that they have different parameters. Here we derive the relations between the entries of the matrix X(x) in (8.2) and those of Y (y) in (1.20). We will consider the case where N is odd and even separately. Let π j (y) be the orthogonal polynomials for the weight w N (y) in (1.14) and P n (x) be the orthogonal polynomials for the weight (8.1) with α = -1 2 . Notice that w N (y) = w N (-y) so by the orthogonality conditions, R π j (y)y j w N (y) dy = (-1) j R π j (-y)y j w N (y) dy = 0; therefore, π j (-y) = (-1) j π j (y). (8.9) From (8.9), we see that π N (y) = p N 2 (y 2 ), where p N 2 is a polynomial of degree N/2. Let us make the change of variables As the difference between u 2,N and its limit u 2 does not affect the asymptotic results, we shall replace u 2,N by u 2 . By the orthogonality of π N (y), we can see that (omitting the constant) Now, let us look at X 12 (0) and Y 12 (0). As N is even, since it is simply the integration of an odd function over a symmetric interval. But if we expand Y 12 (y) near the origin, we see that The term with j = 2 in the above is given by We can make the same change of variables, q 2 = 2q 1 /N , to obtain Expanding Y (y) near the origin, we can express X 12 (0) in terms of Y 12 (0), i.e. Y 12 (0). (8.12) In order to complete the mapping Y (y) → X(x), we are only left to determine the relation between the orthogonality constants. We have Finally, the change of variable and so When N is odd, we set α = 1/2 and replace n = (N -1)/2 in the transformation (8.10). Thus, we have where p N -1 2 is a polynomial of degree (N -1)/2. The orthogonality conditions for π N (y) now give As previously, the monic polynomials orthogonal with respect to W 1 2 (x) are The relation between Y 11 (y) and X 11 (x) follows immediately: Differentiating both sides with respect to y gives Therefore, Y 11 (0). (8.17) When N is odd the integral is different from zero. Equation (8.14) and the change of variables q 2 = 2q 1 /N give .18) The relation between the orthogonality constants can be found in the same way as for N even; we obtain Now that we have obtained the relations between entries in X(x) and Y (y), we can express the elements of Y (y) in terms of v X (ω). In turn, as N → ∞ and N 1 2 z → √ u 2 such relations will allow us to derive formulae for the canonical coordinates P 1 and Q 1 involving the PIII transcendent v Y (r). This completes the proof of Theorem 2.2. In addition, the canonical coordinates P 1 and Q 1 can be written as ) When Im(ζ) < 0, formula (3.9) allows write the outer parametrix S ∞ (y) in terms of ζ: where we have used the expansion In order to prove Eqs. (8.24), note that from (4.23a) and (4.24) we have . Since u 1 = 0, P 2 = a 3 = 0; therefore, we can write We are now in a position to write P 1 and Q 1 in terms of the transcendent v Y (r). We shall outline the calculation for N even; the derivation when N is odd is almost identical. Using (8.33) The initial conditions (2.5a) and (2.5b) provide the asymptotic behaviour of v -1 Y (r) as r → 0 and prove (8.7). Higher order terms can be computed; however, they are rather complicated and we shall not write them down explicitly. Differentiating v -1 Y (r) and using the relations (2.7a) between P 1 , Q 1 and the Hamilton equations (2.8) gives dv Y dr = iv 2 Y (r) 1 (r 2 )r 4 Q 2 1 (r 2 ) 2r 3 + v 2 Y (r). (8.40) Equation (8.39) can be rearranged as Inserting these expressions into Eq. (8.40) we find P 1 and Q 1 in terms of v -1 Y , v Y and r. Namely, we have ) This completes the proof of Proposition 8.1. We are grateful to Professors Chen and Its for making their manuscript available to us. For the sake of simplicity, we use the same notation for the Lax pair (4.22) and for the generic matrices (4.27) and (4.28). It will be evident from the context which ones we are referring to.