Deformed Cauchy random matrix ensembles and large $N$ phase transitions

We study a new hermitian one-matrix model containing a logarithmic Penner's type term and another term, which can be obtained as a limit from logarithmic terms. For small coupling, the potential has an absolute minimum at the origin, but beyond a certain value of the coupling the potential develops a double well. For a higher critical value of the coupling, the system undergoes a large $N$ third-order phase transition.


Introduction
The statistical ensembles of random matrices have a vast number of applications in various domains [1,2,3,4]. Originally introduced by Wigner [5] -motivated by the spectral properties of nuclear resonance levels -they rose in the late 70's as a tool to elucidate aspects of the non-perturbative structure of QCD in the large N expansion [6]. Random matrices can be used to generate sums over random surfaces, a feature that led to non-perturbative formulations of two-dimensional gravity and non-critical strings [7]. More generally, matrix models are very efficient in computations of combinatorial graph generating functions [8,9]. This property has been exploited in different contexts, including classifying numbers of RNA complexes of an arbitrary topology [10,11]. The topological expansion of the free energy also provides an ideal framework to understand detailed aspects of resurgent analysis and large N instantons [12,13,14]. Large N expansions can be recursively generated by loop equations [15], which are encoded in a spectral curve.
Large N random matrices, in particular, multi-matrix ensembles, have also been used to model interesting statistical systems, including Ising [16,17] and Potts [18,19] models, O(N ) models [20,21], among others. Other important aspects of random matrices include the description in terms of a conformal field theory and relations with integrable hierarchies (see e.g. [22]). The Penner model [23] is a hermitian one-matrix model with a potential containing polynomial terms and logarithmic terms. It has been exhaustively investigated, including the classification of critical points to all genera [24]. In this paper we will explore random matrix models arising as a limit of models with potentials having logarithmic singularities in the complex plane. We will focus on a particular theory of a general class of models having large N quantum phase transition of the third order at a critical coupling. Large N phase transitions are familiar in random matrix models [25,26] and singularities typically reflect the finite convergence radius of the planar expansion. The model is introduced in section 2, where the large N limit is studied. The construction hints that the present matrix model might be exactly solvable at finite N , though solving the model is beyond the scope of this paper. A discussion of finite N partition functions is included in section 3.

The model
We consider an hermitian N × N matrix M with the dynamics governed by the potential For stability, we must require A > 0. However, convergence of the partition function will imply the stronger condition A > N − 1.

The 3-parameter deformed Cauchy model
Before studying the dynamics of the matrix theory (1) for its own sake, it is instructive to elucidate on the close relation to the familiar Penner random matrix models, where the potential is of the form where U 0 is a polynomial (for example, the Gaussian potential U 0 = c x 2 ). In a number of cases, the Penner model can be exactly solvable, with the partition function computed for any N , using Selberg's integral formula [27] and generalizations [1]. Alternatively, the partition function can be computed by recursion relations using the method of orthogonal polynomials [7,28] (see [29] for applications of this method to Penner models). We consider the potential This is in the class of Penner models (2) with the choice The particular case with A = N was recently studied in section 3 of [30], where it originated from a family of unitary matrix models through the map to the unit circle. This case is special, as we shall discuss, and in the large N theory corresponds to the marginal case for stability of the model. When = 0, the potential (3) describes a Cauchy ensemble, studied in [31,30].
The potential (3) has one absolute minimum at x = 0 for β ≤ A(1+ ) and two minima at x = ± β A − 1 − for β > A(1 + ). We will refer to the hermitian one-matrix theory with potential (3) as the 3-parameter deformed Cauchy model, or 3-parameter biphasic Cauchy model, since the theory has a large N phase transition on a critical line in parameter space. In this paper we will not study the phase transitions of this model, but instead focus on the matrix theory with potential (1). We will return to the 3-parameter biphasic Cauchy model in section 3.
To connect with theory (1), we take the limit β → ∞, → 0, with fixed β ≡ B. This exactly gives our model (1). The close relation with Penner models suggests that the model (1) could be exactly solvable at finite N (see section 3).

Partition function at large N
The partition function is given by After gauge fixing to a diagonal U (N ) matrix M = diag(a 1 , ..., a N ), the partition function becomes where We note that the integral is convergent provided A > N − 1. In the large N theory, this condition will naturally arise in order to avoid an instability in the distribution of eigenvalues.
The dynamics of the system can be intuitively understood as follows. There is the usual repulsive force between eigenvalues produced by the Vandermonde determinant. In addition, each eigenvalue is subject to the potential The first term, with A > 0, produces an attractive force that pushes the eigenvalue towards the origin. The second term produces an attractive force for B < 0 and a repulsive force for B > 0. The potential has an absolute minimum at the origin for B < A. However, when B > A, in the vicinity of the origin the repulsive force overcomes the attractive force and the vacuum at x = 0 becomes unstable: the potential develops a double well, with Figure 1: The potential has a minimum at x = 0 for B < A and two minima at x = We will study the large N limit with A, B → ∞ and fixed couplings λ, κ defined by We introduce the eigenvalue density It satisfies the normalization condition where L is the region of the complex x plane where eigenvalues condense.
In the large N limit, the partition function can be computed exactly by the saddle-point method. In this limit, the saddle-point equations reduce to the following singular integral equation Starting with small κ, one expects a one-cut eigenvalue distribution, as the potential has an absolute minimum at x = 0 and eigenvalues will be pushed towards the minimum. However, the potential develops a double well when κ > λ. As κ is further increased there should be critical point where the eigenvalues get separated into a symmetric two-cut distribution. This would imply that the system undergoes a phase transition at a critical κ c . In what follows we confirm this picture by explicitly solving the integral equation in the two phases.

The one-cut solution
For sufficiently small κ, the eigenvalues are expected to condense in one cut, −µ < x < µ, with a density satisfying the normalization condition (10) In order to solve the saddle-point equation, as usual one introduces the resolvent Then the density is determined by This leads to The integral can be computed by residues by considering a contour surrounding the cut. There is no pole at infinity, and the result is given by the residues at z = ±i. We obtain the following expression for the eigenvalue density The parameter µ representing the width of the eigenvalue distribution is determined by the normalization condition (12). Computing this integral, we obtain the condition This leads to a cubic equation for X ≡ µ 2 . A valid solution for a one-cut distribution requires that µ 2 is real and that ρ(x) is non-negative in the interval −µ < x < µ. For λ > 1, the normalization condition always has a real solution for µ, irrespective of the value of κ. If λ < 1, then the system is unstable at large N ; the repulsive (Vandermonde) force of eigenvalues dominate over the attractive force and eigenvalues spread out to infinity. The mathematical origin of the condition λ > 1 is the convergence condition of the (large N ) matrix integral (5). 2 Finally, the case λ = 1 is marginal and will be discussed separately in section 2.6. Next, let us consider the condition that ρ(x) is non-negative. From (16) we see that ρ becomes negative in an interval when where µ is determined in terms of λ and κ by the normalization condition (17). This condition indeed determines the critical value of the phase transition. Note that κ c /λ is greater than one, which implies that the transition occurs after the potential developed a double well. We can use (18) to find the critical µ: Substituting into (17), we find At the critical point, the eigenvalue density becomes The eigenvalue density is shown in fig. 2 for κ < κ c and the critical eigenvalue density is shown in fig. 3.
The critical coupling κ c lies in the interval λ < κ c < 2λ (in particular, this ensures that µ 2 c > 0). For example, for λ = 2, the potential has a double-well when κ > 2. However, eigenvalues are still distributed along one cut until the critical value κ c ≈ 3.73, where µ c ≈ 3.6. As shown below, for κ > κ c , eigenvalues get distributed symmetrically in a two-cut distribution.  Another interesting limit is the case κ = 0, where the density takes the form

The two-cut solution
Above the critical coupling, κ > κ c , one has to search for a two-cut solution. We assume that the eigenvalue density has support in two disconnected regions (−a, −b) and (b, a), with real a, b and 0 < b < a. The density is assumed to be of the form For a two-cut (Z 2 -symmetric) solution, the resolvent is given by The contour C is the union of two contours surrounding the two cuts. The integral can be computed by residues, and it is contributed by the poles at z = ±i. We find The parameters a and b representing the endpoints of the eigenvalue distribution can be computed by demanding two conditions: 1) normalization and 2) the asymptotic condition obeyed by the resolvent One can get the equivalent condition by substituting the above solution (23), (25) into the integral equation. The integral equation then implies that the residues at infinity coming from the two terms with coefficients λ and κ must cancel. This leads to the condition The normalization condition leads to the following relation between parameters Using condition (27), this simplifies to Equations (27) and (29) can be explicitly solved for a and b. One obtains Substituting these formulas in (25), the eigenvalue density dramatically simplifies, This represents the exact eigenvalue distribution in the supercritical phase in the stability regime λ > 1.
The critical point of the transition corresponds to b = 0. This gives in agreement with the critical point obtained in the subcritical regime. At the critical point, and ρ matches the critical density ρ c , given in (21), obtained from the subcritical phase. Finally, at large κ, the density takes the asymptotic shape shown in fig. 4(b), with Thus a 2 , b 2 grow linearly with κ with a/b tending to a fixed value for given λ.

Critical behavior
To understand the nature of the phase transition we now study the analytic properties of the free energy F = − ln Z at the critical point. The free energy can be computed in the subcritical and supercritical regimes by using the one-cut and two-cut eigenvalue densities obtained above. We shall fix λ and increase κ until it overcomes the critical value, where the system undergoes a phase transition. Instead of computing the free energy, it is more convenient to consider the first derivative Subcritical regime κ < κ c The integral (35) can be computed by residues using (16) and choosing a contour surrounding the cut from −µ to µ. We obtain where µ = µ(λ, κ) is implicitly defined by the condition (17) (we omit the explicit expression given in terms of a solution of a cubic equation).
Higher derivatives of the free energy in the subcritical phase can be computed by differentiating (36) with respect to κ. We need ∂µ 2 /∂κ, which is obtained by differentiating the normalization condition (17). We get .
Differentiating (36) and using (37) we find Supercritical regime κ > κ c Let us now compute the integral (35) by residues using (31). We choose a contour which is the union of two contours surrounding the two cuts from −a to −b and from b to a. We find the remarkably simple formula: By differentiating with respect to κ, we obtain second and third derivatives We can now examine the continuity properties of derivatives of the free energy at the critical point. For the first and second derivatives we obtain Therefore the first and second derivatives of the free energy are continuous at the transition point. For the third derivative, at the critical point we find and This is different from zero for any λ in the region of stability of the theory λ > 1 (it has a zero at λ = 1/2). Thus we conclude that the system undergoes a third-order large N phase transition. The susceptibility χ = − ∂ 2 F ∂κ 2 is continuous at the transition point, but its derivative has a jump. This is shown in fig. 5.

The marginal case λ = 1
The case λ = 1 is the marginal case for stability. It corresponds to having A = N + O(1) in the original coupling. The eigenvalue density can be obtained from the formulas of the previous subsections by taking the limit λ → 1.
In the subcritical case, λ → 1 gives µ → ∞, that is, eigenvalues are spread from −∞ to ∞. The critical coupling is In the supercritical case, a → ∞ and b 2 becomes The resulting eigenvalue densities in the subcritical and supercritical case are and One can check that they satisfy the normalization condition for any κ. The densities are shown in figs. 6(a),(b). At the critical point, the free energy exhibits the same non-analytic behavior as in the case λ > 1, with a discontinuous third derivative. The jump is obtained from (43) by setting λ = 1.

Finite N partition functions
Let us consider the 3-parameter deformed Cauchy matrix model with potential (3) at finite N , i.e. before the limit → 0. The partition function is given by The integral is a complicated generalization of the Selberg's integral, with singularities at different places. The particular case A = N was studied recently in [30].
In the case = 0, one has the Cauchy random matrix ensemble (see e.g. [31]). Note that for any given A all moments Tr M 2n with n ≥ 1 2 + A − N are divergent. This is of course the same pathology that is present in the Cauchy probability distribution.
This feature is also present in the above three-parameter deformed Cauchy matrix model. It is worth noting that the large N theory does not suffer from this pathology as long as A > N : then all moments x 2n are finite. The problem appears at large N when A ≤ N . In the marginal case A = N , discussed in section 2.6, already the moment x 2 is ill-defined. However, we stress that there are observables which are well defined, such as the free energy and its derivatives, computed in the previous section.
For any given N , the integral (48) can be carried out explicitly in terms of hypergeometric functions. In the U (1) case, we obtain where α ≡ A.
It is interesting to take the limit β → ∞, → 0, with fixed B = β , where the model reduces to the theory with potential (7) studied here. The partition function reduces to This integral can be computed by taking the → 0 limit on the result (49). For this, one first uses Kummer's transformations to bring the hypergeometric functions to hypergeometric functions with argument − . Then one considers the Taylor series expansion of the hypergeometric in powers of and uses the Stirling-de Moivre formula for the Gamma functions with large arguments in the coefficients of the series. Upon taking the → 0 limit with fixed B and resumming the resulting series, the hypergeometric functions become confluent hypergeometric functions and one finally obtains One can check numerically that this is indeed the exact formula for the integral (50). Thus we have computed the partition function for our theory in the U (1) case. For positive integer A, Z U (1) 0 is expressed in terms of Bessel functions. In the U (2) case, one obtains where J is the integral This is also expressed in terms of 2 F 1 hypergeometric functions.
A challenging problem is to derive a closed formula for arbitrary N . A standard method to compute the partition function is through orthogonal polynomials and recursion relations. In the particular case = 0, the partition function substantially simplifies and it can be computed using Romanovski polynomials [31,30].
The ensemble (48) with general A, β, and N seems to have been overlooked in the literature. It would be extremely interesting understand its different limits and phase structure.