ar
X
iv
:1
00
5.
12
98
v1
[
ma
th.
CA
]
7 M
ay
20
10
THE LOWEST EIGENVALUE OF JACOBI RANDOM MATRIX
ENSEMBLES AND PAINLEVE´ VI
EDUARDO DUEN˜EZ, DUC KHIEM HUYNH, JON P. KEATING, STEVEN J. MILLER,
AND NINA C. SNAITH
Abstract. We present techniques for the numerical calculation of the distribution
of the lowest eigenvalue of random matrices in a Jacobi ensemble. Our method is
to solve numerically the associated Painleve´ VI nonlinear differential equation with
suitable initial conditions that we determine. Alternatively, we compute a power-series
expansion of the same Painleve´ VI function. Our results are applied in a forthcoming
paper [4] in which we model the distribution of the first zero above the central point
of elliptic curve L-function families of finite conductor and of conjecturally orthogonal
symmetry.
Contents
1. Introduction 1
2. Jacobi Ensembles 2
3. Distribution of the first eigenvalue and Painleve´ VI 3
4. Initial Conditions for the Painleve´ Equation 6
5. Numerical Painleve´ VI solver 17
6. Alternate numerical solution using power series 20
7. Numerical integration vs. series expansion 25
8. Summary 26
References 27
1. Introduction
We present techniques for calculating numerically the distribution of the lowest eigen-
value (or synonymously, we say the ‘first eigenvalue’) of random matrices in Jacobi en-
sembles J (a,b)N . We proceed as follows. We introduce the Jacobi ensemble JN = J
(a,b)
N of
N ×N random matrices. We relate the distribution of the lowest eigenvalue of matrices
in JN to the probability E(α,β)N (0; I) of having no eigenvalues in some interval I. The
method we use to obtain a numerical approximation for E(α,β)N is through its interpreta-
tion as the Okamoto τ -function of a Painleve´ VI system, for whose auxiliary Hamiltonian
Date: May 11, 2010.
2010 Mathematics Subject Classification. 34M55 (primary), 15B52, 33C45, 65F15 (secondary).
Key words and phrases. Jacobi ensembles, Painleve´ VI, Selberg-Aomoto’s integral.
The second-named author was partially supported by EPSRC, a CRM postdoctoral fellowship and
NSF grant DMS-0757627. The fourth-named author was partially supported by NSF grant DMS-
0855257. The final author was supported by funding from EPSRC. The first, second and fourth authors
thank the University of Bristol for its hospitality, where much of the write-up was done.
1

2 DUEN˜EZ, HUYNH, KEATING, MILLER, AND SNAITH
function h(t) (the values of which in the interval 0 ≤ t ≤ 1 are the ones of interest) For-
rester and Witte [6] have established an explicit differential equation of Painleve´ VI
type. Using their result together with the Selberg-Aomoto integral, we obtain explicit
initial conditions valid close to the edge t = 1 for the Painleve´ equation characterizing
h(t). With these in hand we provide code using MATLAB’s built-in ordinary differential
solver (which uses a Runge-Kutta method) to evaluate h(t) numerically.
As an independent alternate approach, we use the symbolic power series manipula-
tion abilities of SAGE and Maxima to obtain a power series expansion for h(t) at the
other edge t = 0 by solving an explicit recurrence relation for each successive, hitherto
unknown, coefficient in terms of the prior ones. The recurrence relation is obtained
through the Painleve´ equation.
Finally, we analyze the range of parameters a, b,N where both methods are stable in
the sense that both numerically-computed solutions agree in some subinterval [u, v] of
(0, 1), which implies that the numerical solver is robust for this range of parameters.
It is important to note that the methods we use apply to non-integer values of N .
Painleve´ differential equations have played a role in many problems in random matrix
theory, ranging from the distribution of the eigenvalues in the bulk to the largest and
smallest eigenvalues, and have been extensively studied; for our purposes, the most
relevant are the investigations of solutions to Painleve´ VI. We briefly mention some of
the literature. We refer the reader to the special edition of the Journal of Physics A
(Volume 39, Number 39, 2006), which celebrates 100 years of Painleve´ VI, especially
the historical introduction and survey [15] and the article by Forrester and Witte [7] on
connections with random matrix theory; see also the recent works by Dai and Zhang [2]
and Chen and Zhang [1] for determinantal formulas obtained from ladder operators.
The main contribution of this paper is the derivation of an algorithm to compute
numerically the distribution of the lowest eigenvalue in the Jacobi ensembles, and a
discussion of its implementation and accuracy. The motivation for this project comes
from attempts to understand the observations in [14] on the distribution of the first zero
above the central point in families of elliptic curve L-functions when the conductors
are small. The Katz-Sarnak conjectures [9, 10] predict that as the conductors of the
elliptic curves tend to infinity their zero statistics should agree with the N → ∞ scaling
limits of the corresponding statistics of the eigenvalues of matrices from a classical
compact group. For suitable test functions this was proved in [13,18]; however, for finite
conductors the numerical data in [14] is in sharp disagreement with the limiting behavior
of these random matrix ensembles. In particular, the first zero above the central point
is repelled, with the repulsion decreasing as the conductors increase. In the companion
paper [4] we complete the study of the low lying zeros of elliptic curve L-functions, and
obtain a model which describes the behavior of these zeros for finite conductors. One
of the key ingredients in our model is the lowest eigenvalue of these Jacobi ensembles
of N ×N matrices, often requiring non-integer values of N , which is the main result of
this paper.
2. Jacobi Ensembles
Let JN = J (a,b)N denote the Jacobi ensemble on N levels 0 ≤ xj ≤ 1, j = 1, 2, . . . , N ,
with real parameters a, b > −1. Explicitly, the N -level (joint) probability density of

THE LOWEST EIGENVALUE IN JACOBI ENSEMBLES AND PAINLEVE´ VI 3
levels of JN on [0, 1]N with respect to its Lebesgue measure dx1 dx2 · · · dxN is given by
(2.1) C˜(a,b)N
N∏
j=1
W (xj)
∏
1≤j<k≤N
(xk − xj)2
where the weight function W = W (a,b) on [0, 1] is given by
(2.2) W (x) = xb(1− x)a
and C˜(a,b)N is the ensemble’s normalization constant. Jacobi ensembles as described
above correspond to suitable ensembles of self-dual random matrices via the angular
variables φj defined by
(2.3) xj =
1 + cosφj
2 , 0 ≤ φj ≤ pi.
Note that the edges x = 0, x = 1 correspond respectively to φ = pi, φ = 0 under this
change of variables. We refer the reader to [3] and the forthcoming monograph [8] for
details regarding the matrix realizations of Jacobi ensembles, for which we will otherwise
have no direct use. In what follows we will go back and forth between the abscissæ xj
and the angular variables φj, but will in any case refer to the associated ensemble by the
Jacobi name and denote it by JN . In terms of the angular variables, and with respect
to Lebesgue measure on (0, pi)N , the N -level (joint) probability density for JN is given
by
(2.4) C(α,β)N
N∏
j=1
w(cos φj)
∏
1≤j<k≤N
(cosφk − cosφj)2
with the weight function w = w(α,β) on (0, pi),
(2.5) w(cos φ) = (1− cosφ)α(1 + cosφ)β .
The parameters α, β > −12 are related to a, b above by α = a + 1/2, β = b + 1/2,
and C(α,β)N is the appropriate normalization constant, namely C
(α,β)
N = 2−(N+a+b+2)N
C˜(α−1/2,β−1/2)N for the constant C˜
(a,b)
N of (2.1). For suitable choices for α and β we obtain
the joint probability density of the N independent eigenphases for the classical groups
of matrices SO(2N), SO(2N + 1) and USp(2N), when the latter are endowed with an
invariant (Haar) probability measure and regarded as random matrix ensembles. The
case α = β = 0 corresponds to SO(2N), α = 1 and β = 0 corresponds to SO(2N + 1),
and α = β = 1 to USp(2N). This is explained in detail in [3]. Below in (4.23) we give an
explicit expression for the normalization constant C˜(a,b)N . (Jacobi-distributed pseudoran-
dom sequences of levels can be generated from a uniform pseudorandom sequence using
only the Jacobi joint probability density via for instance the Accept-Reject Algorithm
whose applicability is quite broad; see for instance [16].)
3. Distribution of the first eigenvalue and Painleve´ VI
As remarked above, the Jacobi ensemble JN describes the eigenvalue statistics in
suitable ensembles of self-dual random matrices having N pairs of eigenvalues e±iφj ,
j = 1, 2, . . . , N ; we call φ the eigenphase of the eigenvalue eiφ. Let E(α,β)N (n; I) denote
the probability that a random matrix A ∈ JN has exactly n eigenphases in the inter-
val I = [0, φ]. As shorthand notation we write E(α,β)N (φ) for E
(α,β)
N (0; [0, φ]), namely

4 DUEN˜EZ, HUYNH, KEATING, MILLER, AND SNAITH
the probability of having no eigenphases in the interval [0, φ]. The probability density
function ν(α,β)N (φ) of the distribution of the first eigenphase is related to E
(α,β)
N (φ) by
(3.1) ν(α,β)N (φ) = −
d
dφE
(α,β)
N (φ).
We can deduce the relation (3.1) as follows: assume that the interval [0, φ] contains no
eigenvalues. Then a small increment ε > 0 of the interval to [0, φ + ε] has two possible
outcomes. Either the interval [0, φ+ε] contains no eigenvalues, or it contains some. The
probability of the first event is E(α,β)N (φ + ε). It follows that E
(α,β)
N (φ) − E
(α,β)
N (φ + ε)
is the probability that the interval [φ, φ + ε] contains at least one eigenvalue; as ε → 0
there can be only one, namely the first eigenphase in [φ, φ + ε]. Thus
(3.2) lim
ε→0
E(α,β)N (φ)− E
(α,β)
N (φ + ε)
ε = −
d
dφE
(α,β)
N (φ)
indeed yields the probability density function ν(α,β)N (φ) of the first eigenphase. An
alternative way to prove (3.1) is to observe that 1 − E(α,β)N (φ) is the probability that
[0, φ] contains at least one eigenphase, hence that the first eigenphase φmin is at most φ;
otherwise said 1−E(α,β)N (φ) is the cumulative distribution function of the first eigenphase
φmin, so its derivative is equal to the probability density function ν(α,β)N (φ). In Section 5
we shall need to scale the angular variable φ by a factor of N/pi in order to consider
eigenphases of mean unit spacing on [0, N ].
We have
E(α,β)N (φ) = C
(α,β)
N
∫ pi
φ
· · ·
∫ pi
φ
N∏
j=1
(1− cosφj)α(1 + cosφj)β
×
∏
1≤j<k≤N
(cosφj − cosφk)2 dφ1 · · · dφN
(3.3)
for fixed α, β > −1/2 and the normalization constant C(α,β)N of (2.4). There is no known
method to evaluate the multiple integral in equation (3.3) exactly. E(α,β)N (φ) is related
to a Painleve´ VI transcendental function h(t), namely a certain solution to a second-
order nonlinear ordinary differential equation. In Proposition 4.4 we provide the first
few terms of a power-series expansion of E(α,β)N (φ) for φ close to 0; these provide the
initial conditions for the differential equation we aim to solve. Our main reference for
the theory is the work of Forrester and Witte [6]. Their result is stated for the abscissal
counterpart to the function E(α,β)N (φ) of (3.3), namely the function E˜
(a,b)
N (t) defined by
(3.4) E˜(a,b)N (t) = C˜
(a,b)
N
∫ t
0
· · ·
∫ t
0
N∏
j=1
xbj(1− xj)a
∏
1≤j<k≤N
(xj − xk)2dx1 · · · dxN
with the normalization constant C˜(a,b)N of (2.1). The functions (3.4) and (3.3) are related
by the change of variables
(3.5) a = α− 1/2 and b = β − 1/2
along with
(3.6) t = 1 + cosφ2 and xj =
1 + cosφj
2

THE LOWEST EIGENVALUE IN JACOBI ENSEMBLES AND PAINLEVE´ VI 5
where 0 ≤ t, xj ≤ 1. Explicitly,
(3.7) E(α,β)N (φ) = E˜
(a,b)
N
(1 + cosφ
2
)
.
The method we use to obtain a numerical approximation for E˜(a,b)N is via its interpre-
tation as an Okamoto τ -function and ensuing relation to a Painleve´ system with associ-
ated auxiliary Hamiltonian h(t); this Hamiltonian arises as the solution of a Painleve´ VI
equation with the exact parameters determined by Forrester and Witte in Proposition
13 of [6] as follows.
Proposition 3.1. Let a, b > −1 and N be a positive integer. The auxiliary Hamiltonian
(3.8) h(t) = t · e′2[b]− 12e2[b] + t(t− 1) ddt log E˜
(a,b)
N (t)
where
b = (b1, b2, b3, b4) =
(
N + a+ b2 ,
a− b
2 ,−
a + b
2 ,−N −
a+ b
2
)
,(3.9)
e2[b] = b1b2 + b1b3 + b1b4 + b2b3 + b2b4 + b3b4,(3.10)
e′2[b] = b1b3 + b1b4 + b3b4,(3.11)
satisfies the following Painleve´ VI equation in Jimbo-Miwa-Okamoto σ-form
(3.12) h′(t)
(
t(1− t)h′′(t)
)2 +
(
h′(t)[2h − (2t− 1)h′(t)] + b1b2b3b4
)2 =
4∏
k=1
(h′(t) + b2k).
Furthermore, we have the boundary condition (as t → 0)
(3.13) h(t) =
(
−12e2[b]−N(b +N)
)
+
(
e′2[b] +
N(N + b)(2N + a+ b)
2N + b
)
t+O(t2).
Note that besides simplifying the notation in Proposition 13 of [6] we also swap the a’s
and b’s therein. (Note that a determines the order of vanishing of the Jacobi eigenphase
density at the edge φ = 0.)
Following Edelman and Persson [5], we seek to compute simultaneously E˜(a,b)N (t) and
the Hamiltonian h(t) via a (non-autonomous) differential equation for the triple of func-
tions
(3.14) H(t) =


E˜(a,b)N (t)
h(t)
h′(t)


of the form
(3.15) dHdt = Ft(H)
with initial conditions
(3.16) H0 = H(t0) =


E˜(a,b)N (t0)
h(t0)
h′(t0)


where t0 = 1 − ε for small ε > 0. (The singularity of the Painleve´ equation and its
solutions at t0 = 1 preclude taking simply ε = 0.)

6 DUEN˜EZ, HUYNH, KEATING, MILLER, AND SNAITH
From (3.8) we obtain
(3.17) ddt E˜
(a,b)
N (t) =
h(t)− te′2[b] + 12e2[b]
t(t− 1) E˜
(a,b)
N (t),
and likewise from (3.12)
(3.18) h′′(t) = 1t(1− t)
√∏4
j=1(h′(t) + b2j)− (h′(t)[2h − (2t− 1)h′(t)] + b1b2b3b4)2
h′(t) .
Therefore, (3.15) holds with
(3.19)
Ft


h1(t)
h2(t)
h3(t)


=


h2(t)− te′2[b] + 12e2[b]
t(t− 1) h1(t)
h3(t)
1
t(1− t)
√∏4
j=1(h3(t) + b2j )− (h3(t)[2h2(t)− (2t− 1)h3(t)] + b1b2b3b4)2
h2(t)


.
The MATLAB code to compute Ft(H) is described in Section 5 and is also available for
download at http://www.maths.bris.ac.uk/∼mancs/publications.html. We em-
ploy the built-in ordinary differential solver ode45 from MATLAB, which implements a
Runge-Kutta method giving an approximate solution of (3.15) (and thus of the sought
density ν(α,β)N (φ) = − ddφE
(α,β)
N (φ) of the distribution of the first eigenphase). It remains
still to determine the initial condition H0 = H(t0) for (3.14) as follows. We shall find
the small-ε asymptotic behavior of E˜(a,b)N (1 − ) (equivalently, what we will actually do
is find the small-φ asymptotics of E(α,β)N (φ)). By differentiation we then find
(3.20) ddtE˜
(a,b)
N (t) and
d2
dt2 E˜
(a,b)
N (t)
and thus (through its definition (3.8)) we obtain asymptotically good approximations
to h(t0) and h′(t0) for any t0 close to 1. This gives the triple H0 of initial conditions
for (3.14). In the following section we compute the asymptotic behavior of E(α,β)N (φ) for
φ close to 0.
4. Initial Conditions for the Painleve´ Equation
In this section we compute the probability E(α,β)N (φ) that a random matrix from a
Jacobi ensemble JN as defined in Section 2 has no eigenphase in the interval [0, φ] for
small φ > 0. As described at the end of Section 3, this enables us to derive the initial
conditions for the system of differential equations (3.15) which gives the distribution
of the first eigenphase (3.1). Forrester and Witte [6] consider this same limit in their
equation (1.38), but here we derive a further term in the approximation. Our result is
stated in Proposition 4.4.

THE LOWEST EIGENVALUE IN JACOBI ENSEMBLES AND PAINLEVE´ VI 7
We require some notation. For n = 1, . . . , N and α, β > −1/2 we define the integral
I(n) := C(α,β)N
∫ pi
0
· · ·
∫ pi
0︸ ︷︷ ︸
N−n times
∫ φ
0
· · ·
∫ φ
0︸ ︷︷ ︸
n times
N∏
j=1
(1− cosφj)α(1 + cosφj)β
×
∏
1≤j<k≤N
(cos φj − cosφk)2dφ1 · · · dφN .
(4.1)
Then we have the following lemmata:
Lemma 4.1. For E(α,β)N (φ) as given in (3.3) and for I(n) as defined in (4.1) we have
(4.2) E(α,β)N (φ) = 1−N · I(1) +
(N
2
)
I(2)−
(N
3
)
I(3) + · · · + (−1)N I(N).
Lemma 4.2. For I(1) as defined in (4.1) we have
(4.3) I(1) = H1
φ2α+1
2α + 1 −
[
(N − 1)H2 +
( α
12 +
β
4
)
H1
] φ2α+3
2α + 3 +O(φ
2α+4)
where
(4.4) H1 :=
Γ(α +N + 1/2)Γ(α + β +N)
22αΓ(α + 1/2)Γ(α + 3/2)Γ(N + 1)Γ(β +N − 1/2)
and
(4.5) H2 :=
Γ(α +N + 1/2)Γ(α + β +N + 1)
22α+1Γ(N + 1)Γ(β +N − 1/2)Γ(α + 1/2)Γ(α + 5/2) .
Lemma 4.3. For I(n) as defined in (4.1) we have for n ≥ 2 and α > −1/2
(4.6) I(n)  φ2αn+2n2−n.
We postpone the proof of the above lemmata for a moment in order to state the
desired Taylor series expansion of E(α,β)N (φ) for small φ > 0.
Proposition 4.4. For E(α,β)N (φ) as given in (3.3) we have
E(α,β)N (φ) = 1−N
(
H1
φ2α+1
2α + 1 −
[
(N − 1)H2 +
( α
12 +
β
4
)
H1
] φ2α+3
2α + 3
)
+ O(φ2α+4)
(4.7)
where H1 is defined in (4.4) and H2 in (4.5).
Proof. Let n ≥ 2. We can apply Lemma 4.3 to every term I(n) in (4.2) of Lemma 4.1.
Thus, I(n)  φ2αn+2n2−n. As n ≥ 2, we have 2αn + 2n2 − n ≥ 2α + 4 for every
α > −1/2, hence I(n)  φ2α+4. Thus every I(n) in Lemma 4.1 can be absorbed into
the error term O(φ2α+4) of I(1) in (4.3) of Lemma 4.2. 
Now we prove Lemmata 4.1, 4.2 and 4.3.
Proof of Lemma 4.1. Recall that E(α,β)N (φ) is the probability of having no eigenphases
in the interval [0, φ]. We denote the event that there is at least one eigenphase in [0, φ]
by B. Then its complement {B is the event of having no eigenphases in [0, φ]. Hence
(4.8) E(α,β)N (φ) = P ({B) = 1− P (B).

8 DUEN˜EZ, HUYNH, KEATING, MILLER, AND SNAITH
We now focus on P (B). Let
(4.9) Bk := {(φ1, . . . , φk, . . . , φN ) ∈ [0, pi]N where φk ∈ [0, φ]},
which is the event that φk lies in [0, φ] and the remaining eigenphases lie anywhere in
[0, pi]. Note that Bj and Bk, j 6= k, are not necessarily disjoint.
We can write the event of having at least one eigenphase in [0, φ] as
(4.10) B =
N⋃
k=1
Bk.
For the probability of the event Bk we have
P (Bk) = C(α,β)N
∫
Bk
N∏
j=1
(1− cosφj)α(1 + cosφj)β
∏
1≤j<k≤N
(cosφj − cosφk)2 dφj .(4.11)
For any k and any j 6= k we have
(4.12) P (Bk) = I(1) and P (Bj ∩Bk) = I(2),
and in general, for 1 ≤ i1 < · · · < ik ≤ N ,
(4.13) P (Bi1 ∩ · · · ∩Bik) = I(k).
By the inclusion-exclusion principle and the symmetry above,
P (B) = P
( N⋃
k=1
Bk
)
=
N∑
k=1
(−1)k+1
∑
1≤i1<···<ik≤N
P (Bi1 ∩ · · · ∩Bik)
= N · P (B1)−
(N
2
)
P (B1 ∩B2) +
(N
3
)
P (B1 ∩B2 ∩B3)
+ · · · + (−1)N+1P (B1 ∩ · · · ∩BN ).
(4.14)
Thus
P (B) = N · I(1)−
(N
2
)
I(2) +
(N
3
)
I(3) + · · · + (−1)N+1I(N).(4.15)
Finally, the probability of having no eigenphases in [0, φ]N is given by
P ({B) = C(α,β)N
∫ pi
φ
· · ·
∫ pi
φ
N∏
j=1
(1− cosφj)α(1 + cosφj)β
×
∏
1≤j<k≤N
(cosφj − cosφk)2dφ1 · · · dφN
= 1− P (B)
= 1−N · I(1) +
(N
2
)
I(2) −
(N
3
)
I(3) + · · ·+ (−1)N I(N).(4.16)


THE LOWEST EIGENVALUE IN JACOBI ENSEMBLES AND PAINLEVE´ VI 9
Proof of Lemma 4.2. First we recall Selberg’s integral (see Chapter 17 of [12]):
SN (ρ, η; γ) :=
∫ 1
0
dx1 · · ·
∫ 1
0
dxN
N∏
l=1
xρ−1l (1− xl)η−1
∏
1≤j<k≤N
|xk − xj |2γ
=
N−1∏
j=0
Γ(1 + γ + jγ)Γ(ρ + jγ)Γ(η + jγ)
Γ(1 + γ)Γ(ρ+ η + [N + j − 1]γ)(4.17)
and Aomoto’s extension of Selberg’s integral for 1 ≤ R ≤ N
∫ 1
0
dx1 · · ·
∫ 1
0
dxN
R∏
j=1
xj
N∏
l=1
xρ−1l (1− xl)η−1
∏
1≤j<k≤N
|(xj − xk)|2γ
=
R∏
j=1
ρ+ (N − j)γ
ρ+ η + (2N − j − 1)γ
N−1∏
j=0
Γ(1 + γ + jγ)Γ(ρ + jγ)Γ(η + jγ)
Γ(1 + γ)Γ(ρ + η + [N + j − 1]γ) ,(4.18)
both valid for integer N and complex ρ, η, γ with
Re(ρ) > 0, Re(η) > 0, Re(γ) > −min
( 1
N ,
Re(ρ)
(N − 1) ,
Re(η)
(N − 1)
)
.(4.19)
The version of Selberg’s integral we are interested in is related to (4.17) by
ρ = s+ 1/2, η = r + 1/2, γ = 1,
and the change of variables
yj =
1 + cosφj
2
as follows (note Γ(2) = 1):
∫ pi
0
dφ1 · · ·
∫ pi
0
dφN
N∏
l=1
(1− cosφl)r(1 + cosφl)s
∏
1≤j<k≤N
(cosφj − cosφk)2
= 2N (N+r+s−1)
∫ 1
0
dy1 · · ·
∫ 1
0
dyN
N∏
l=1
(1− yl)r−1/2ys−1/2l
∏
1≤j<k≤N
(yj − yk)2
= 2N (N+r+s−1)
N−1∏
j=0
Γ(2 + j)Γ(s + 1/2 + j)Γ(r + 1/2 + j)
Γ(s+ r +N + j) .
(4.20)
The version of Aomoto’s extension of our interest has parameters
ρ = r + 1/2, η = s+ 1/2, γ = 1,
and we change variables
zj =
1− cosφj
2

10 DUEN˜EZ, HUYNH, KEATING, MILLER, AND SNAITH
in (4.18) to obtain
∫ pi
0
dφ1 · · ·
∫ pi
0
dφN
R∏
k=1
(1− cosφk)
N∏
l=1
(1− cosφl)r(1 + cosφl)s
∏
1≤j<k≤N
(cosφj − cosφk)2
= 2R+N (N+r+s−1)
∫ 1
0
dz1 · · ·
∫ 1
0
dzN
R∏
k=1
zk
N∏
l=1
zr−1/2l (1− zl)s−1/2
∏
1≤j<k≤N
(zj − zk)2
= 2R+N (N+r+s−1)
R∏
j=1
r + 1/2 +N − j
r + s+ 2N − j
N−1∏
j=0
Γ(2 + j)Γ(r + 1/2 + j)Γ(s + 1/2 + j)
Γ(r + s+N + j) .
(4.21)
We can now determine the normalization constant C(α,β)N in (4.11). We have
(4.22)
C(α,β)N
−1
=
∫ pi
0
dφ1 · · ·
∫ pi
0
dφN
N∏
j=1
(1− cosφj)α(1 + cosφj)β
∏
1≤j<k≤N
(cosφj − cosφk)2.
By setting N = N, r = α, and s = β in (4.20) we obtain
(4.23) C(α,β)N
−1
= 2N(N+α+β−1)
N−1∏
j=0
Γ(2 + j)Γ(β + 1/2 + j)Γ(α + 1/2 + j)
Γ(α+ β +N + j) .
Now, for small φ > 0 and using Selberg’s integral (4.20) we wish to evaluate
(4.24) I(1) = C(α,β)N I˜(1)
where
I˜(1) :=
∫ pi
0
· · ·
∫ pi
0
∫ φ
0
K(φ2, . . . , φN )H(φ1, . . . , φN )dφ1dφ2 · · · dφN
with
(4.25) K(φ2, . . . , φN ) :=
N∏
j=2
(1− cosφj)α(1 + cosφj)β
∏
2≤j<k≤N
(cosφj − cosφk)2,
and
H(φ1, . . . , φN ) := (1− cosφ1)α(1 + cosφ1)β
N∏
k=2
(cosφ1 − cosφk)2.(4.26)
Now we evaluate H(φ1, . . . , φN ). The Taylor expansion around φ1 = 0 of the first
factor of H(φ1, . . . , φN ) in (4.26) is
(1− cosφ1)α(1 + cosφ1)β =
[φ2α1
2α − α
φ2α+21
2α · 12 +O(φ
2α+4
1 )
] [
2β − β2β−2φ21 +O(φ41)
]
= φ
2α
1
2α−β −
( α
2α−β · 12 +
β
2α−β+2
)
φ2α+21 +O(φ2α+41 ).
(4.27)

THE LOWEST EIGENVALUE IN JACOBI ENSEMBLES AND PAINLEVE´ VI 11
The Taylor expansion of the terms in the second factor defining H(φ1, . . . , φN ) in (4.26)
is
N∏
k=2
(cosφ1 − cosφk)2 =
N∏
k=2
[(
1− φ
2
1
2! +O(φ
4
1)
)
− cosφk
]2
=
N∏
k=2
[
(1− cosφk)2 − φ21(1− cosφk) +O(φ41)
]
=
N∏
k=2
(1− cosφk)2 − φ21


N∑
j=2
N∏
k=2
(1− cosφk)2
1− cosφj

+O(φ41).
(4.28)
Using (4.27) and (4.28) in (4.26) gives
H(φ1, . . . , φN ) =
[ φ2α1
2α−β −
( α
2α−β · 12 +
β
2α−β+2
)
φ2α+21 +O(φ2α+41 )
]
×


N∏
k=2
(1− cosφk)2 − φ21


N∑
j=2
N∏
k=2
(1− cosφk)2
1− cosφj

+O(φ41)

 .
(4.29)
Multiplying out and collecting terms by powers of φ1 gives
H(φ1, . . . , φN ) =
φ2α1
2α−β
N∏
k=2
(1− cosφk)2 −
φ2α+21
2α−β


N∑
j=2
N∏
k=2
(1− cosφk)2
1− cosφj


−
( α
2α−β · 12 +
β
2α−β+2
)
φ2α+21
N∏
k=2
(1− cosφk)2 +O(φ2α+41 ).
(4.30)
Integration of the last expression (4.30) gives
∫ φ
0
H(φ1, . . . , φN )dφ1 =
N∏
k=2
(1− cosφk)2
φ2α+1
(2α + 1)2α−β
−


N∑
j=2
N∏
k=2
(1− cosφk)2
1− cosφj

 φ
2α+3
(2α + 3)2α−β
−
N∏
k=2
(1− cosφk)2
( α
2α−β · 12 +
β
2α−β+2
) φ2α+3
2α + 3 +O(φ
2α+5).
(4.31)
Hence, to evaluate I˜(1) =
∫ pi
0 · · ·
∫ pi
0
∫ φ
0 K(φ2, . . . , φN )H(φ1, . . . , φN )dφ1 · · · dφN , we have
to compute
∫ pi
0
· · ·
∫ pi
0
K(φ2, . . . , φN )
N∏
k=2
(1− cosφk)2dφ2 · · · dφN and(4.32)
∫ pi
0
· · ·
∫ pi
0
K(φ2, . . . , φN )


N∑
j=2
N∏
k=2
(1− cosφk)2
1− cosφj

 dφ2 · · · dφN .(4.33)

12 DUEN˜EZ, HUYNH, KEATING, MILLER, AND SNAITH
Observe that the integrand of (4.33) is symmetric in its variables φ2, . . . , φN . There-
fore we have
(4.34)
∫ pi
0
· · ·
∫ pi
0
K(φ2, . . . , φN )


N∑
j=2
N∏
k=2
(1− cosφk)2
1− cosφj

 dφ2 · · · dφN
= (N − 1)
∫ pi
0
· · ·
∫ pi
0
K(φ2, . . . , φN )(1 − cosφ2)
N∏
k=3
(1− cosφk)2dφ2 · · · dφN .
Evaluating the integral (4.32) yields
(4.35)
∫ pi
0
· · ·
∫ pi
0
K(φ2, . . . , φN )
N∏
k=2
(1− cosφk)2dφ2 · · · dφN
=
∫ pi
0
· · ·
∫ pi
0
N∏
j=2
(1− cosφj)α+2(1 + cosφj)β
∏
2≤j<k≤N
(cosφj − cosφk)2dφ2 · · · dφN ,
and using Selberg’s integral (4.20) with N = N − 1, r = α + 2, and s = β gives
(4.36)
∫ pi
0
· · ·
∫ pi
0
K(φ2, . . . , φN )
N∏
k=2
(1− cosφk)2dφ2 · · · dφN
= 2(N−1)(N+α+β)
N−2∏
j=0
Γ(2 + j)Γ(β + 1/2 + j)Γ(α + 5/2 + j)
Γ(α + β +N + 1 + j) .
Normalizing the last expression (4.36) with C(α,β)N from (4.22) we obtain
(4.37) C(α,β)N
∫ pi
0
· · ·
∫ pi
0
K(φ2, . . . , φN )
N∏
k=2
(1− cosφk)2dφ2 · · · dφN
= 12α+β
Γ(α+N + 1/2)Γ(α + β +N)
Γ(α + 1/2)Γ(α + 3/2)Γ(N + 1)Γ(β +N − 1/2) .
We now evaluate the integral in (4.34). For this we note that
(4.38) (1− cosφ2)
N∏
k=3
(1− cosφk)2 =
N∏
k=2
(1− cosφk)
N∏
k=3
(1− cosφk).

THE LOWEST EIGENVALUE IN JACOBI ENSEMBLES AND PAINLEVE´ VI 13
So
∫ pi
0
· · ·
∫ pi
0
K(φ2, . . . , φN )(1 − cosφ2)
N∏
k=3
(1− cosφk)2 dφ2 · · · dφN
=
∫ pi
0
· · ·
∫ pi
0
K(φ2, . . . , φN )
N∏
k=2
(1− cosφk)
N∏
k=3
(1− cosφk)dφ2 · · · dφN
=
∫ pi
0
· · ·
∫ pi
0
N∏
j=2
(1− cosφj)α(1 + cosφj)β
∏
2≤j<k≤N
(cosφj − cosφk)2
×
N∏
k=2
(1− cosφk)
M∏
k=3
(1− cosφk)dφ2 · · · dφN
=
∫ pi
0
· · ·
∫ pi
0
N∏
k=3
(1− cosφk)
N∏
j=2
(1− cosφj)α+1(1 + cosφj)β
×
∏
2≤j<k≤N
(cosφj − cosφk)2dφ2 · · · dφN .
(4.39)
With R = N − 2,N = N − 1, r = α + 1, and s = β in Aomoto’s integral (4.21) this
evaluates to
∫ pi
0
· · ·
∫ pi
0
K(φ2, . . . , φN )(1− cosφ2)
N∏
k=3
(1− cosφk)2dφ2 · · · dφN
= 2(N−2)+(N−1)(N+α+β−1)
N−2∏
j=1
α +N + 1/2 − j
α + β + 2N − 1− j
×
N−2∏
j=0
Γ(2 + j)Γ(α + 3/2 + j)Γ(β + 1/2 + j)
Γ(α + β +N + j)
= 2(N−2)+(N−1)(N+α+β−1) Γ(α +N + 1/2)Γ(α + 5/2)
Γ(α + β +N + 1)
Γ(α + β + 2N − 1)
Γ(α +N − 1/2)
Γ(α + 1/2)
×
N−2∏
j=0
Γ(2 + j)Γ(α + 3/2 + j)Γ(β + 1/2 + j)
Γ(α + β +N + j) .
(4.40)
Normalizing the last expression (4.40) with C(α,β)N from (4.22) we obtain
C(α,β)N
∫ pi
0
· · ·
∫ pi
0
K(φ2, . . . , φN )(1 − cosφ2)
N∏
k=3
(1− cosφk)2dφ1 · · · dφN
= 2−α−β−1 Γ(α +N + 1/2)Γ(α + β +N + 1)Γ(N + 1)Γ(β +N − 1/2)Γ(α + 1/2)Γ(α + 5/2) .
(4.41)

14 DUEN˜EZ, HUYNH, KEATING, MILLER, AND SNAITH
Putting (4.37) and (4.41) into (4.24) we obtain
I(1) = H˜1
φ2α+1
(2α + 1)2α−β − (N − 1)H˜2
φ2α+3
(2α + 3)2α−β
−
( α
2α−β · 12 +
β
2α−β+2
)
H˜1
φ2α+3
2α + 3 +O(φ
2α+5)
(4.42)
where
(4.43) H˜1 :=
1
2α+β
Γ(α +N + 1/2)Γ(α + β +N)
Γ(α + 1/2)Γ(α + 3/2)Γ(N + 1)Γ(β +N − 1/2)
and
(4.44) H˜2 :=
1
2α+β+1
Γ(α +N + 1/2)Γ(α + β +N + 1)
Γ(N + 1)Γ(β +N − 1/2)Γ(α + 1/2)Γ(α + 5/2) .
Finally, we rewrite the result slightly and get
I(1) = H1
φ2α+1
2α + 1 − (N − 1)H2
φ2α+3
2α + 3 −
( α
12 +
β
4
)
H1
φ2α+3
2α + 3 +O(φ
2α+5)
= H1
φ2α+1
2α + 1 −
[
(N − 1)H2 +
( α
12 +
β
4
)
H1
] φ2α+3
2α + 3 +O(φ
2α+5)
(4.45)
where
(4.46) H1 :=
Γ(α +N + 1/2)Γ(α + β +N)
22αΓ(α + 1/2)Γ(α + 3/2)Γ(N + 1)Γ(β +N − 1/2)
and
(4.47) H2 :=
Γ(α +N + 1/2)Γ(α + β +N + 1)
22α+1Γ(N + 1)Γ(β +N − 1/2)Γ(α + 1/2)Γ(α + 5/2) .
This completes the proof of Lemma 4.2.

Proof of Lemma 4.3. The definition of I(n) is
I(n) = C(α,β)N
∫ pi
0
· · ·
∫ pi
0︸ ︷︷ ︸
N−n times
∫ φ
0
· · ·
∫ φ
0︸ ︷︷ ︸
n times
N∏
j=1
(1− cosφj)α(1− cosφj)β
×
∏
1≤j<k≤N
(cosφj − cosφk)2dφ1 · · · dφN .
(4.48)
Here we are only interested in the size of I(n) in terms of n, so we can disregard the
normalization constant C(α,β)N . For n ≥ 2 we consider
I˜(n) := C(α,β)N
−1
I(n).

THE LOWEST EIGENVALUE IN JACOBI ENSEMBLES AND PAINLEVE´ VI 15
Then we have
I˜(n) =
∫ pi
0
dφn+1 · · ·
∫ pi
0
dφN
N∏
j=n+1
(1− cosφj)α(1 + cosφj)β
×
∏
n+1≤j<k≤N
(cosφj − cosφk)2
×
∫ φ
0
{
(1− cosφn)α(1 + cosφn)β
N∏
j=n+1
(cos φn − cosφj)2
×
∫ φ
0
(1− cosφn−1)α(1 + cosφn−1)β
N∏
j=n+1
(cos φn−1 − cosφj)2
× (cos φn−1 − cosφn)2dφn−1
...
×
∫ φ
0
(1− cosφ1)α(1 + cosφ1)β
N∏
j=n+1
(cosφ1 − cosφj)2
× (cos φ1 − cosφ2)2 × (cosφ1 − cosφ3)2 × · · · × (cosφ1 − cosφn)2dφ1
}
dφn.
(4.49)
Now for j, k ≤ n, j 6= k and φj , φk ∈ [0, φ] for small φ > 0 we have
(4.50) (cosφj − cosφk)2 ≤ (1− cosφ)2.
There are (n − 1) + · · · + 1 = n(n − 1)/2 terms of the form (cos φj − cosφk)2 (with
j = 1, . . . , n − 1, and k = j + 1, . . . , n) occurring in (4.49), each of which we bound
using (4.50).
From equation (4.31) in the proof of Lemma 4.2 we derive for k = 1, . . . , n that
∫ φ
0
(1− cosφk)α(1 + cosφk)β
N∏
j=n+1
(cos φk − cosφj)2dφk
=
N∏
j=n+1
(1− cosφj)2
φ2α+1
2α−β(2α + 1) +O(φ
2α+3).
(4.51)

16 DUEN˜EZ, HUYNH, KEATING, MILLER, AND SNAITH
Using (4.50) and (4.51) we obtain
I˜(n) ≤
∫ pi
0
dφn+1 · · ·
∫ pi
0
dφN
N∏
j=n+1
(1− cosφj)α(1 + cosφj)β
×
∏
n+1≤j<k≤N
(cos φj − cosφk)2(1− cosφ)n(n−1)
×


N∏
j=n+1
(1− cosφj)2 ×
φ2α+1
2α−β(2α + 1) +O(φ
2α+3)


n
=
∫ pi
0
dφn+1 · · ·
∫ pi
0
dφN
N∏
j=n+1
(1− cosφj)α(1 + cosφj)β
×
∏
n+1≤j<k≤N
(cos φj − cosφk)2(1− cosφ)n(n−1)
×


N∏
j=n+1
(1− cosφj)2n ×
φ(2α+1)n
[2α−β(2α + 1)]n +O(φ
n(2α+1)+2)


=
∫ pi
0
dφn+1 · · ·
∫ pi
0
dφN
N∏
j=n+1
(1− cosφj)α+2n(1 + cosφj)β
×
∏
n+1≤j<k≤N
(cos φj − cosφk)2(1− cosφ)n(n−1)
×
[
φ(2α+1)n
[2α−β(2α + 1)]n +O(φ
n(2α+1)+2)
]
.
(4.52)
Hence for some constants c1 and c2
I˜(n) ≤ c1(1− cosφ)n(n−1)
[
φ(2α+1)n +O(φn(2α+1)+2)
]
= c2
[
φ2 +O(φ4)
]n(n−1) [φ(2α+1)n +O(φn(2α+1)+2)
]
= c2
[
φ2n(n−1) +O(φ2n(n−1)+2)
] [
φ(2α+1)n +O(φn(2α+1)+2)
]
= c2φ(2α+1)n+2n(n−1) +O(φ(2α+1)n+2n(n−1)+2)
= c2φ2αn+2n
2−n +O(φ(2α+1)n+2n(n−1)+2).
(4.53)
Hence I˜(n)  φ2αn+2n2−n, which implies I(n)  φ2αn+2n2−n as required to finish the
proof. 
Remark 4.5. It is natural to ask whether we can determine further terms of the Taylor
expansion of E(α,β)N (φ) than provided in Proposition 4.4. By Lemma 4.1 this requires
us to determine I(n) with n ≥ 2. However, for I(n) with n ≥ 2 we encounter multiple
integrals which cannot be dealt with using Selberg’s or Aomoto’s integral. To our
knowledge the required generalization of Selberg’s integral does not exist.

THE LOWEST EIGENVALUE IN JACOBI ENSEMBLES AND PAINLEVE´ VI 17
5. Numerical Painleve´ VI solver
We combine the results from the previous sections and provide a MATLAB algorithm
computing a numerical approximation to the probability density function ν(α,β)N (φ) of
the first eigenphase of random matrices in the Jacobi ensemble JN = J (α,β)N . We follow
some of the ideas in Edelman and Persson [5] for the implementation.
The code snippet below shows the MATLAB definition of the function p6diff which
computes Ft(H) as given in (3.19); p6diff takes as inputs the parameters b1, b2, b3,
b4 (corresponding to b1, b2, b3, b4 from (3.9)), e2 (corresponding to e2[b] from (3.10)),
and e2p (corresponding to e′2[b] from (3.11)).
function dH = p6diff(t,H,b1,b2,b3,b4,e2p,e2)
dH = zeros(3,1);
dH(1) = H(1)/t/(t-1)*(H(2)-e2p*t+e2/2);
dH(2) = H(3);
dH(3) = sqrt(((H(3)+b1^2)*(H(3)+b2^2)*(H(3)+b3^2)*(H(3)+b4^2)...
-(H(3)*(2*H(2)-(2*t-1)*H(3))+b1*b2*b3*b4)^2)/H(2))/t/(1-t);
To determine the initial conditions (3.16) for some t0 = 1− ε with small ε > 0 we use
Proposition 4.4, according to which
(5.1)
E(α,β)N (φ) = 1−N
(
H1
φ2α+1
2α + 1 −
[
(N − 1)H2 +
( α
12 +
β
4
)
H1
] φ2α+3
2α + 3
)
+O(φ2α+4),
with
(5.2) H1 =
Γ(α +N + 1/2)Γ(α + β +N)
22αΓ(α + 1/2)Γ(α + 3/2)Γ(N + 1)Γ(β +N − 1/2)
and
(5.3) H2 =
Γ(α +N + 1/2)Γ(α + β +N + 1)
22α+1Γ(N + 1)Γ(β +N − 1/2)Γ(α + 1/2)Γ(α + 5/2) .
Via the substitution φ = cos−1(2t−1), equation (5.1) provides a good approximation
to E˜(a,b)N (t) for t small, which we code in MATLAB as the function EN (note that r,s
correspond to α, β):
function E = EN(r,s,N,phi)
E=1-N*Hone(r,s,N)*phi.^(2*r+1)/(2*r+1)...
+N*(N-1)*Htwo(r,s,N)*phi.^(2*r+3)/(2*r+3)...
+N*(r/12.0+s/4.0)*Hone(r,s,N)*phi.^(2*r+3)/(2*r+3);
Calling EN(r,s,N,phi) with phi = arccos(2t0 − 1) for small t0 gives a good numerical
approximation to the first component of the initial condition (3.16). The MATLAB code
for the auxiliary functions Hone(r,s,N) and Htwo(r,s,N) defined by (5.2) and (5.3) is
function H1 = Hone(r,s,N)
H1=gamma(r + N + 1/2)/gamma(N + 1)*gamma(r + s + N)/...
gamma(r + 3/2)/gamma(r + 1/2)/gamma(s + N - 1/2)/2^(2*r);
function H2 = Htwo(r,s,N)
H2=gamma(r + N + 1/2)/gamma(N + 1)*gamma(r +s +N + 1)/...
gamma(r + 5/2)/gamma(r + 1/2)/gamma(s + N - 1/2)/2^(2*r+1);
In view of the definition (3.8) of the auxiliary Hamiltonian h, we need to differentiate
E˜(a,b)N (t) twice with respect to t in order to obtain the initial conditions h(t0), h′(t0) (the

18 DUEN˜EZ, HUYNH, KEATING, MILLER, AND SNAITH
second and third components of (3.16)). We call ENderiv1 (resp., ENderiv2) the imple-
mentation of (an asymptotic approximation to) the function E(α,β)N
′
(φ) = ddφE
(α,β)
N (φ)
(resp., of E(α,β)N
′′
(φ) = d2dφ2E
(α,β)
N (φ)) obtained by differentiation of (5.1). The MATLAB
code is
function Ed = ENderiv1(r,s,N,phi)
Ed = -N*Hone(r,s,N)*phi.^(2*r) ...
+ N*(N-1)*Htwo(r,s,N)*phi.^(2*r+2) ...
+ N*(r/12.0+s/4.0)*Hone(r,s,N)*phi.^(2*r+2);
function Edd = ENderiv2(r,s,N,phi)
Edd = -2*r*N*Hone(r,s,N)*phi.^(2*r-1) ...
+ N*(2*r+2)*(N-1)*Htwo(r,s,N)*phi.^(2*r+1) ...
+ N*(r/12.0+s/4.0)*(2*r+2)*Hone(r,s,N)*phi.^(2*r+1);
Since φ = cos−1(2t − 1) implies dφ/dt = −1/
√
t(1− t) and E˜(a,b)N (t) = E
(α,β)
N (φ), the
chain rule and equation (3.8) give
h(t) = t · e′2[b]− 12e2[b] + t(t− 1)
d
dtE˜
(a,b)
N (t)
E˜(a,b)N (t)
= t · e′2[b]− 12e2[b] +
√
t(1− t)E
(α,β)
N
′
(cos−1(2t− 1))
E(α,β)N (cos−1(2t− 1))
(5.4)
and
h′(t) = e′2[b] +
1− 2t
2
√
t(1− t)
E(α,β)N
′
(cos−1(2t− 1))
E(α,β)N (cos−1(2t− 1))
+
√
t(1− t)
E(α,β)N (cos−1(2t− 1))
(
−E(α,β)N
′′
(cos−1(2t− 1))√
t(1− t)
)
+
√
t(1− t)E(α,β)N
′
(cos−1(2t− 1))
(
E(α,β)N
′
(cos−1(2t− 1))
√
t(1− t)E(α,β)N (cos−1(2t− 1))2
)
= e′2[b] +
1− 2t
2
√
t(1− t)
E(α,β)N
′
(cos−1(2t− 1))
E(α,β)N (cos−1(2t− 1))
− E
(α,β)
N
′′
(cos−1(2t− 1))
E(α,β)N (cos−1(2t− 1))
+ E
(α,β)
N
′
(cos−1(2t− 1))2
E(α,β)N (cos−1(2t− 1))2
.
(5.5)
With (5.1), (5.4) and (5.5) we have the initial conditions and code them in MATLAB
where t0 is close to 1 (e.g. t0 = 1 - e-5) as
r = a+0.5;
s = b+0.5;
phi0 = acos(2*t0-1);
root0 = sqrt(t0*(1-t0));
E0 = EN(r,s,N,phi0);
E0d = ENderiv1(r,s,N,phi0);
E0dd = ENderiv2(r,s,N,phi0);
H0 = [ E0; ...
t0*e2p - 0.5*e2 + root0*E0d/E0; ...

THE LOWEST EIGENVALUE IN JACOBI ENSEMBLES AND PAINLEVE´ VI 19
e2p + (0.5-t0)/root0 * E0d/E0 - E0dd/E0 + (E0d/E0)^2 ];
We now have all ingredients for the main program painleve6(a,b,N) which computes
an approximated solution for the system (3.15). Typing
(5.6) function [t,H,theta,Fp] = painleve6(a,b,N)
defines our main program which returns the variables in the square brackets on the
left-hand side. Here t, H correspond to t and H, and theta corresponds to the rescaled
angular variable θ = Nφ/pi. The output variable Fp in (5.6) is a vector of values of the
rescaled distribution
(5.7) pi2N sin
(piθ
N
)E˜(a,b)N (t)
t(t− 1)
[
h(t)− te′2[b] + 12e2[b]
]
of the first eigenphase (the rescaling achieves mean unit spacing of the N eigenphases
θ1, . . . , θN on [0, N ]), as obtained from (3.17) via the rescaled variable θ (= theta),
namely
Fp = (pi/2/N)*sin(pi*theta/N).*F./t./(t-1).*(h-e2p*t+e2/2);
where F is the first component (row) of H (so the entries of F provide values of E(α,β)N (φ)),
and h, the second component of H, has as entries the values of the Painleve´ function h,
thus
F = H(:,1);
h = H(:,2);
Below follows some self-explanatory code for (3.9), (3.10), and (3.11)
b1 = (a+b)/2+N;
b2 = (a-b)/2;
b3 = -(a+b)/2;
b4 = -(a+b)/2-N;
e2p = b1*b3 + b1*b4 + b3*b4;
e2 = e2p + b2*(b1+b3+b4);
The syntax for MATLAB’s solver ode45 is
(5.8) ode45(function, domain, initial conditions, options, parameters).
The ‘options’ argument gives the user limited control over the MATLAB algorithm by
specifying the admissible relative and absolute tolerances. At the i-th step MATLAB
ensures that the approximated error e(i) satisfies
(5.9) |e(i)| ≤ max(RelTol*abs(y(i)),AbsTol(i))
where y(i) is the approximation for the exact value. Obviously the smaller we choose
the relative and absolute tolerance the better the approximation becomes at the cost of
running time and memory. We set the relative and absolute tolerance as follows
opts=odeset(’reltol’,1e-5,’abstol’,1e-6).(5.10)
Finally we call the solver ode45 with the syntax given in (5.8):
[t,H]=ode45(@p6diff,[t0,0.01],H0,opts,b1,b2,b3,b4,e2p,e2);
Here we apply ode45 to the function p6diff with domain from t0 = 1 - 1e-5 down
to 0.01 (which means that theta ranges approximately from 0 to 1) with initial con-
ditions H0 defined in (3.16), options opts set in (5.10) and parameters b1, b2, b3, b4,
e2p, e2.
At long last, the command
[t H theta Fp] = painleve6(a,b,N)

20 DUEN˜EZ, HUYNH, KEATING, MILLER, AND SNAITH
computes the numerical solution to the system of differential equations (3.15) for the
Jacobi ensemble J (a,b)N . The subsequent command
plot(theta, Fp)
plots the distribution Fp (as defined in (5.7)) of the first rescaled eigenphase θmin for JN .
Notice that a=-0.5 and b=-0.5 corresponds to SO(2N). Likewise setting a=0.5 and
b=-0.5 gives SO(2N +1) and finally USp(2N) would correspond to choosing a=0.5 and
b=0.5.
The full MATLAB code can be obtained from the authors or downloading from our
webpages, such as
http://www.maths.bris.ac.uk/∼mancs/publications.html.
6. Alternate numerical solution using power series
In this section we describe an algorithm to compute the power series expansion of
the Painleve´ function h(t) at t = 0, leading to the numerical computation of E(α,β)N (φ)
for φ close to pi. It is rather unfortunate (for our intended application) that t = 1 is
a branch-point singularity of h(t), and as a consequence a power-series expansion of h
about t = 1 is a Puisseaux series (i. e., a series in fractional powers of t), at least if
the parameters a, b (and, eventually, N) are rational numbers; for arbitrary values of
the parameters the situation would be even more complicated. Therefore, we content
ourselves for the time being with the power series expansion about t = 0.
The idea of the algorithm is very simple: The coefficients h0, h1 of the expansion
h(t) = h0 + h1t + h2t2 + . . . are given in equation (3.13). These are used to bootstrap
a recursive search for the higher coefficients h2, h3, . . . , regarding each unknown hk as
implicitly defined by the earlier coefficients h0, h1, h2, . . . , hk−1 through the Painleve´
equation (3.12). Given the complicated nonlinear nature of the Painleve´ equation it
is not immediately obvious that this approach will work in practice, but fortunately
it does, and each successive coefficient hk is expressible as a rational function of the
previous ones, hence ultimately as a rational function of the parameters a, b,N . Once
many terms are computed, the power series for h(t) can be used to evaluate E˜(a,b)N (t) by
solving for the latter in equation (3.8); then we use the change of variables t → φ and
differentiation to compute the density ν(α,β)N (φ) of the distribution of the first eigenvalue,
at least for φ relatively close to pi.
We seek to find the coefficients hk as exact rational numbers; for this reason, we will
work with exact (truncated) power series with rational coefficients. This approach has
the enormous advantage that the complicated operations needed to evaluate both sides
of the Painleve´ equation (3.12) introduce no numerical errors at all in the evaluation of
successive higher coefficients. The price paid is a more expensive calculation compared
to one done using exclusively floating-point arithmetic. Remarks on the choice of the
number of needed terms to reach adequate numerical precision follow below in section 7.
The algorithm is written in SAGE [17], although Maxima [11] plays an important role
behind the scenes.
As pointed out above, we take the parameters a, b and N to be (fixed) rational
numbers, and regard the function
(6.1) h(t) = h0 + h1t+ h2t2 + · · ·

THE LOWEST EIGENVALUE IN JACOBI ENSEMBLES AND PAINLEVE´ VI 21
as having coefficients which are rational numbers. In particular, from equation (3.13)
we have
h0 = −
1
2e2[b]−N(b +N)(6.2)
h1 = e′2[b] +
N(N + b)(2N + a+ b)
2N + b(6.3)
with e2[b], e′2[b] given in terms of a, b,N by (3.9), (3.10) and (3.11).
The following lines of SAGE code take care of some preliminaries:
a=-1/2
b=-1/2
n=4
b1=(a+b)/2+n
b2=(a-b)/2
b3=-(a+b)/2
b4=-(a+b)/2-n
e2p = b1*b3 + b1*b4 + b3*b4
e2 = e2p + b2*(b1+b3+b4)
Note that the algorithm’s implementation depends crucially on inputting rational values
for the parameters a, b, n on the lines above. With a limitation explained below, the
algorithm can handle integral as well as rational values of the parameter n = N .
When, at any given point, the coefficients h0, . . . , hk−1 are known (hence are rational
numbers from SAGE’s perspective), but hk is still unknown, we wish to regard the latter
as an indeterminate, say hk = X. Let
(6.4) h(t) = h0 + h1t+ · · ·+ hk−1tk−1 +Xtk.
Then h lies in the ring Q[X, t] of polynomials in X, t with rational coefficients. Let
us denote by LHS(h) and RHS(h) the polynomials obtained by substitution of (6.4)
in the Painleve´ equation (3.12), and let PEZ(h) = RHS(h) − LHS(h) (“Painleve´-
Equal-to-Z ero”). The natural hope is that if h0, h1, . . . , hk−1 are chosen correctly, then
PEZ(h) = pk(X)tk + O(tk+1) for a non-constant polynomial pk ∈ Q[x] and some
(unspecified) polynomial O(tk+1) divisible by tk+1. Then hk should be chosen to be a
root of pk(X).
Performing exact polynomial arithmetic in Q[X, t] to compute PEZ(h) is quite ex-
pensive, especially since we only need to determine the coefficient p(X) of the lowest
power of t. For computational purposes, however, it is enough to regard h as a finite trun-
cation of an infinite power series in t, systematically neglecting any higher-order terms
not needed for the immediate purpose at hand—namely the determination of pk(X).
Fortunately, SAGE can do algebra in power series rings (with help from Maxima).
Henceforth we work in the ring S = Q[X][[t]] of power series in t whose coefficients
are in Q[X] (polynomials in X with rational coefficients). In order to determine pk(X)
it is enough to know h up to O(tk+2) terms. (It is not quite enough to work modulo tk+1
because pk(X) a priori depends also hk+1, not just on the currently unknown X = hk;
luckily, the solution found a posteriori shows this dependence to be fictitious.)
We thus set
(6.5) h(t) = h0 + h1t+ · · ·+ hk−1tk−1 +Xtk +O(tk+2)
and use SAGE to compute PEZ(h) = pk(X)tk + O(tk+1) for a linear polynomial pk
(the only exception is p2, which is quadratic with a trivial root X = 0). Therefore, for

22 DUEN˜EZ, HUYNH, KEATING, MILLER, AND SNAITH
each k = 2, 3, . . . it suffices to take hk as the unique (nontrivial) root of pk(X), which is
a rational number.
DEGREE=50 # degree (plus 1) of the truncated power series
R.<x> = QQ[’x’] # R = R[x]
S.<t> = PowerSeriesRing(R, default_prec=DEGREE) # S = R[x][[t]]
We store the coefficients hk of the power series h in a SAGE list g, so we have initially
leadexp = n*(n+b)
g = [-e2/2-leadexp, e2p+leadexp*(1+a/(b+2*n))]
h = g[0] + t*g[1]
in accordance with (6.2) and (6.3), where we have introduced the auxiliary constant
leadexp (the reason for the name leadexp will be explained later).
The following loop is the heart of the algorithm. It computes the first DEGREE terms
of the power series expansion of h(t). Note that PEZ is defined according to the Painleve´
equation (3.12).
X = var(’X’) # symbolic Maxima variable
for i in xrange (2,DEGREE):
h = h + x*t^i + O(t^(i+2))
hp = h.derivative()
hpp = hp.derivative()
PEZ = hp*(hpp*t*(1-t))^2 + (hp*(2*h-hp*(2*t-1))+b1*b2*b3*b4)^2 -\
(hp+b1^2)*(hp+b2^2)*(hp+b3^2)*(hp+b4^2)
LHS = PEZ[i](X) # Substitute x=X in p_i(x)
solucion = solve (LHS==0, X) # Find root of p_i(X)
g.append(QQ(solucion[0].rhs())) # The root becomes the next coefficient
# Now reconstruct h up to degree i using the newly-found g[i]
h = g[0]+g[1]*t+g[2]*t^2
for j in xrange (3,i+1):
h += g[j]*t^j
# End of main loop
h = h + O(t^DEGREE)
A couple of remarks are in order. Maxima does not seem to understand that we wish to
interpret the variable x of the power series ring F to be a “symbolic” variable with respect
to which the equation PEZ[i] == 0 (i. e., pi(X) = 0) is to be solved. We therefore need
explicitly to replace the ring’s variable x with the symbolic variable X before finding the
root of PEZ[i], which is then appended at the end of the coefficient list g. We also found
it simpler to reconstruct h from scratch using the coefficient list g rather than figuring
out a way both to (i) increase the precision of h, and (ii) substitute the newly-found
rational coefficient g[i] for x. (The wrapper QQ(...) around the argument of g is used
to convert (“coerce” in SAGE lingo) the root found by Maxima to a bona fide SAGE
rational number.)

THE LOWEST EIGENVALUE IN JACOBI ENSEMBLES AND PAINLEVE´ VI 23
Solving for E˜(a,b)N (t) in the definition (3.8) of the auxiliary Hamiltonian h(t) yields
E˜(a,b)N (t) = exp
(∫ h(t)− t · e′2[b] + 12e2[b]
t(t− 1) dt+ c
)
for a suitable constant c
= C exp
(∫ h0 + 12e2[b] + t · (h(1)(t)− e′2[b])
t(t− 1) dt
)
where C = ec and h(1)(t) = (h(t)− h0)/t = h1 + h2t+ h3t2 + . . .
= C exp
(∫ −N(N + b)
t(t− 1) dt +
∫ h(1)(t)− e′2[b]
t− 1 dt
)
since h0 = −
1
2e2[b]−N(N + b) from (3.13)
= C exp
(
−N(N + b)
∫ ( 1
t− 1 −
1
t
)
dt +
∫ h(1)(t)− e′2[b]
t− 1 dt
)
= C exp
(
N(N + b) log t+
∫ h(1)(t)− e′2[b]−N(N + b)
t− 1 dt
)
= CtN(N+b) exp
(∫ h(1)(t)− e′2[b]−N(N + b)
t− 1 dt
)
.
(6.6)
Note that the integrand in the last expression above is regular at t = 0. We define
(6.7) F(t) = exp
(∫ t
0
h(1)(τ)− e′2[b]−N(N + b)
τ − 1 dτ
)
.
Then F(t) has a series expansion (in integral powers) about 0 with F(0) = 1, and
equation (6.6) reads
(6.8) E˜(a,b)N (t) = CtN(N+b)F(t).
It follows that the leading-order term of the power series expansion of E˜(a,b)N (t) about
t = 0 is CtN(N+b), explaining the nomenclature leadexp (“leading exponent”) for the
auxiliary SAGE variable leadexp = n*(n+b).
The following SAGE code evaluates DEGREE terms of the power series expansion
of F(t) and stores it in F
h1 = (h-h[0])/t
F = h1
F -= e2p+leadexp
F /= t-1+O(t^DEGREE)
F = F.integral()
F = F.exp()
F = F.truncate(DEGREE+1)
The truncate method returns the “polynomial part” of F, in other words, it coerces a
power series into a polynomial by setting all higher-order terms equal to zero.
The value of the constant C in (6.8) can be read off from equation (3.5) in [6] as
the quotient of the normalization constants C˜(a,b)N = 1/SN (a+1, b+1; 1) for the Jacobi

24 DUEN˜EZ, HUYNH, KEATING, MILLER, AND SNAITH
ensembles J (a,b)N and C˜
(0,b)
N = 1/SN (1, b + 1; 1) for J
(0,b)
N (recall that we swap the role of
a and b relative to Forrester and Witte; moreover, C˜(a,b)N = C˜
(b,a)
N ). Explicitly,
(6.9) C = C˜
(a,b)
N
C˜(0,b)N
= SN (1, b + 1; 1)SN (a+ 1, b + 1; 1)
.
The SAGE function Jac(a,b,n) evaluates the Selberg integral SN (a + 1, b + 1; 1) for
integer values of N = n using formula (4.17); the value of C is then computed and
stored in leadcoef.
def Jac(a,b,n):
jac=1.
for i in xrange(n):
jac *= real(gamma(a+i+1.)) * real(gamma(b+i+1.))\
* real(gamma(i+2.)) / real(gamma(a+b+n+i+1.))
return jac
leadcoef = Jac(0,b,n)/Jac(a,b,n)
This part of the code naturally depends on N = n being a positive integer. However,
the calculation of leadcoef is the only dependence of the SAGE code on the integrality
of n. In particular, it would be possible to evaluate Jac(a,b,n) in closed form using
the Barnes G-function (which, unfortunately, is not yet implemented in SAGE) using
the following pseudocode
# *NOT* VALID SAGE CODE UNTIL BARNES G IS IMPLEMENTED!!!
def Jac(a,b,n):
return G(n+a+b+2.)/G(2*n+a+b+1.) \
*G(n+a+1.)/G(a+2.) \
*G(n+b+1.)/G(b+2.) \
*G(n+2.)
Given a proper implementation for the function G, the above would correctly compute
leadcoef and the program would be capable of handling general rational values of n.
We now rewrite (3.17) in the form
(6.10) E˜(a,b)′N (t) =
(
h(1)(t)− e′2[b]
t− 1 −
N(N + b)
t(t− 1)
)
E˜(a,b)N (t),
and compute E˜(a,b)N (t) and its derivative E˜
(a,b)′
N (t) using (6.8) and the above by
htilde = h1 - e2p
htilde = htilde.truncate(DEGREE-1)
def E(u):
return leadcoef * u^leadexp * F(u)
def Ep(u):
return (htilde(u)/(u-1) + leadexp/u/(1-u)) * E(u)
Finally, we compute the cumulative distribution function
1− E(α,β)N (φ) = 1− E˜
(a,b)
N
(1 + cosφ
2
)
and its derivative (cf. (3.1) and (6.10))
ν(α,β)N (φ) = −
d
dφE
(α,β)
N (φ) =
sinφ
2 E˜
(a,b)′
N
(1 + cosφ
2
)
.

THE LOWEST EIGENVALUE IN JACOBI ENSEMBLES AND PAINLEVE´ VI 25
using the SAGE code below, where pcummul computes 1−E˜(a,b)N and nu computes ν
(α,β)
N .
(Note that 12 sinφ =
√
t(1− t) if t = 1+cosφ2 .)
def pcummul(theta):
return 1-E((1+cos(theta))/2)
def nu(theta):
if theta < 1.e-6 or theta > pi - 1.e-6:
return 0
t = (1+cos(theta))/2
return sqrt(t*(1-t)) * Ep(t)
The entire SAGE code can be obtained from the authors and is also available for down-
load at http://www.maths.bris.ac.uk/∼mancs/publications.html.
7. Numerical integration vs. series expansion
As might be expected, our numerical implementation in Matlab of the Painleve´ solver
does not work equally well in all parameter regimes. In this section we describe the tests
we have carried out on the code and the conclusions about the parameter regimes where
a robust solution can be obtained. The Matlab (Runge-Kutta) solver starts from initial
conditions near φ = 0 (that is, t = 1) and numerically extends the computed solution
towards φ = pi (or θ = Npi φ = N in scaled units). The power series, on the other hand,
is an expansion around φ = pi (or, equivalently, t = 0); it is therefore accurate at the
opposite end of the interval on which we are solving. Hence, if the tail of our numerical
solution matches the initial behaviour of the series solution, we are confident that the
numerical solver has worked correctly.
There are three parameters to vary in the input to the numerical solver: a = α− 1/2,
b = β − 1/2 and N . There are also three variables we can adjust in the Matlab code to
try to coax a solution: t0, the starting point near t = 1 (φ = 0) for the numerical solver;
reltol and abstol, which control the accuracy of the numerical solution.
After testing the code for various values of N (integer and non-integer) from about 1
up to 100, it appears that a good solution can be found on a standard desktop machine in
a few seconds with t0 = 1−10−7, reltol= 10−5 and abstol= 10−6 for any −0.5 ≤ a ≤ 0
and −0.5 ≤ b ≤ 0.5 (the range for b that is relevant to the classical groups). In Figure
1 we see examples of code that runs efficiently and matches the series expansion in the
tail of the distribution.
For values of a > 0 the Matlab solver breaks down. Trials with a = 0.001 still work,
but already a = 0.01 fails to produce a good solution. Moving t0 away from 1, for
example to 1 − 10−2, helps achieve a better curve, but as can be seen from Figure 2,
the initial conditions are not close enough to the true curve at this point to produce a
valid solution. Decreasing reltol and abstol, even by a factor of 1000, does not make
a visible difference to the curve. We note that unfortunately this means that while
the Matlab solver works very well for the group SO(2N), we cannot use it to produce
solutions for SO(2N + 1) and USp(2N).
For a > 0 and small values of N a solution valid over the whole interval can still be
glued together by matching the series solution (with a sufficient number of coefficients)
for the tail and bulk of the curve with the known asymptotic behavior near φ = 0. The
right-hand plot in Figure 2 shows that this is certainly possible for N = 5.
The series solution produces a very accurate curve with only 50 terms when N = 2,
but the number of terms needed to obtain a solution which is meaningful over a large

28 DUEN˜EZ, HUYNH, KEATING, MILLER, AND SNAITH
[11] Maxima, a Computer Algebra System. Open source software. http://maxima-project.org.
[12] M.L. Mehta. Random Matrices. Academic Press, 2nd edition, 1991.
[13] S.J. Miller. One- and two-level densities for rational families of elliptic curves: evidence for the
underlying group symmetries. Compos. Math., 140:952–992, 2004.
[14] S.J. Miller. Investigations of zeros near the central point of elliptic curve L-functions. Experiment.
Math., 15 (3):257–279, 2006. (E-print: math.NT/0508150).
[15] Clarkson P.A., Mazzocco M. Joshi N., Nijhoff. F. W., and Noumi M. One hundred years of pVI,
the fuchs-Painleve´ equation. Journal of Physics A, 39, 2006.
[16] Christian Robert and George Casella. Monte Carlo Statistical Methods (Springer Texts in Statistics).
Springer-Verlag, second edition edition, 2004.
[17] SAGE. Open source mathematics software. http://www.sagemath.org.
[18] M. P. Young. Low-lying zeros of families of elliptic curves. J. Amer. Math. Soc., 19:205–250, 2006.
E-mail address: eduenez@math.utsa.edu
Department of Mathematics, University of Texas at San Antonio, San Antonio, TX
78249
E-mail address: dkhuynhms@googlemail.com
Department of Pure Mathematics, University of Waterloo, Waterloo, ON, N2L 3G1,
Canada
E-mail address: J.P.Keating@bristol.ac.uk
School of Mathematics, University of Bristol, Bristol BS8 1TW, UK
E-mail address: Steven.J.Miller@williams.edu
Department of Mathematics and Statistics, Williams College, Williamstown, MA 01267
E-mail address: N.C.Snaith@bristol.ac.uk
School of Mathematics, University of Bristol, Bristol BS8 1TW, UK