Higher-Order Linear Differential Equations for Unitary Matrix Integrals: Applications and Generalisations

Peter J. Forrester; Fei Wei (with an appendix by Folkmar Bornemann)

arXiv:2508.20797·math-ph·February 20, 2026

Higher-Order Linear Differential Equations for Unitary Matrix Integrals: Applications and Generalisations

Peter J. Forrester, Fei Wei (with an appendix by Folkmar Bornemann)

PDF

TL;DR

This paper develops higher-order linear differential equations to characterize unitary matrix integrals, enabling efficient computation and applications to permutation enumeration and Riemann zeta function moments.

Contribution

It introduces a novel approach using higher-order linear differential equations for unitary matrix integrals, extending to beta ensembles and linking to combinatorial and number theory problems.

Findings

01

Matrix differential equations efficiently compute matrix integral expansions.

02

Connections established between matrix integrals and permutation enumeration.

03

Application to moments of derivatives of the Riemann zeta function.

Abstract

In this paper, we consider characterisations of the class of unitary matrix integrals $⟨ (det U)^{q} e^{s^{1/2} Tr (U + U^{†})} ⟩_{U (l)}$ in terms of a first-order matrix linear differential equation for a vector function of size $l + 1$ , and in terms of a scalar linear differential equation of degree $l + 1$ . It will be shown that the latter follows from the former. The matrix linear differential equation provides an efficient way to compute the power series expansion of the matrix integrals, which with $q = 0$ and $q = l$ are of relevance to the enumeration of longest increasing subsequences for random permutations, and to the question of the moments of the first and second derivative of the Riemann zeta function on the critical line, respectively. This procedure is compared against that following from known characterisations involving the…

Tables1

Table 1. Table 1: Cost for computing T l ( N ) T_{l}(N) , 1 ≤ l ≤ N ≤ N ∗ 1\leq l\leq N\leq N_{*} ; single row means fixed l l ; θ \theta , c c in the last two columns are subject to 0 < θ 0<\theta , c < 1 c<1 .

method based on	total cost	single row $l$	single row $l = O (N_{*}^{θ})$	single row $l = c N_{*}$
Chazy-I (B.5)	$O (N_{*}^{3})$	$O ({(N_{} - l)}^{2}) + O (N_{} \log N_{*})$	$O (N_{*}^{2})$	${(1 - c)}^{2} \cdot O (N_{*}^{2})$
recursion (3.12)	$O (N_{*}^{3})$	$O (l \cdot N_{*})$	$O (N_{*}^{1 + θ})$	$c \cdot O (N_{*}^{2})$

Equations226

⟨ e^{z Tr (U + U^{†})} ⟩_{U (l)},

⟨ e^{z Tr (U + U^{†})} ⟩_{U (l)},

1, 2, 6, 24, 119, 694, 4582, 33324, 261808, 2190688, 19318688, 178108704, \dots .

1, 2, 6, 24, 119, 694, 4582, 33324, 261808, 2190688, 19318688, 178108704, \dots .

⟨ e^{z Tr (U + U^{†})} ⟩_{U (l)} = 1 + N = 1 \sum \infty \frac{T _{l} ( N )}{( N ! ) ^{2}} z^{2 N}

⟨ e^{z Tr (U + U^{†})} ⟩_{U (l)} = 1 + N = 1 \sum \infty \frac{T _{l} ( N )}{( N ! ) ^{2}} z^{2 N}

det [I_{j - k} (2 z)]_{j, k = 1, \dots, l} = 1 + N = 1 \sum \infty \frac{T _{l} ( N )}{( N ! ) ^{2}} z^{2 N},

det [I_{j - k} (2 z)]_{j, k = 1, \dots, l} = 1 + N = 1 \sum \infty \frac{T _{l} ( N )}{( N ! ) ^{2}} z^{2 N},

⟨ j = 1 \prod l a (θ_{j}) ⟩_{U (l)}

⟨ j = 1 \prod l a (θ_{j}) ⟩_{U (l)}

= det [\frac{1}{2 π} \int_{0}^{2 π} a (θ) e^{i (j - k) θ} d θ]_{j, k = 1, \dots, l} .

z^{4} Y^{(5)} + 20 z^{3} Y^{(4)} - 2 (10 z - 59) z^{2} Y^{(3)} - 2 (91 z - 110) z Y^{''}

z^{4} Y^{(5)} + 20 z^{3} Y^{(4)} - 2 (10 z - 59) z^{2} Y^{(3)} - 2 (91 z - 110) z Y^{''}

+ 4 (16 z^{2} - 87 z + 20) Y^{'} + 16 (8 z - 5) Y = 0

⟨ (det U)^{l} e^{z Tr (U + U^{†})} ⟩_{U (l)}

⟨ (det U)^{l} e^{z Tr (U + U^{†})} ⟩_{U (l)}

= (- 1)^{l (l - 1) /2} det [I_{j + k + 1} (2 z)]_{j, k = 0, \dots, l - 1}

F_{2} (l) = N \to \infty lim \frac{1}{N ^{l^{2} + 4 l}} ⟨ ∣ Z_{N}^{''} (1) ∣^{2 l} ⟩_{U (N)},

F_{2} (l) = N \to \infty lim \frac{1}{N ^{l^{2} + 4 l}} ⟨ ∣ Z_{N}^{''} (1) ∣^{2 l} ⟩_{U (N)},

Z_{N} (s) := e^{- π i N /2} e^{i \sum_{n = 1}^{N} θ_{n} /2} s^{- N /2} j = 1 \prod N (1 - s e^{- i θ_{j}}) .

Z_{N} (s) := e^{- π i N /2} e^{i \sum_{n = 1}^{N} θ_{n} /2} s^{- N /2} j = 1 \prod N (1 - s e^{- i θ_{j}}) .

F_{2} (l) = (- 1)^{\frac{l ( l - 1 )}{2}} m = 0 \sum 2 l (m 2 l) (\frac{d}{d x})^{4 l - 2 m} (e^{- \frac{x}{2}} x^{- m - \frac{l ^{2}}{2}} f_{m} (x))_{x = 0},

F_{2} (l) = (- 1)^{\frac{l ( l - 1 )}{2}} m = 0 \sum 2 l (m 2 l) (\frac{d}{d x})^{4 l - 2 m} (e^{- \frac{x}{2}} x^{- m - \frac{l ^{2}}{2}} f_{m} (x))_{x = 0},

det [I_{j + k + 1 + 2 μ_{j + 1}} (2 x)]_{j, k = 0, \dots, l - 1},

det [I_{j + k + 1 + 2 μ_{j + 1}} (2 x)]_{j, k = 0, \dots, l - 1},

f_{m} (x) = \frac{1}{x ^{m}} j = 0 \sum m x^{j} P_{j} (x) \frac{d ^{j} τ _{l} ( x )}{d x ^{j}} .

f_{m} (x) = \frac{1}{x ^{m}} j = 0 \sum m x^{j} P_{j} (x) \frac{d ^{j} τ _{l} ( x )}{d x ^{j}} .

E_{β = 2}^{hard} ((0, s); l) := N \to \infty lim E_{N} ((0, s / (4 N)); l) = e^{- s /4} ⟨ e^{(s^{1/2} /2) Tr (U + U^{†})} ⟩_{U (l)} .

E_{β = 2}^{hard} ((0, s); l) := N \to \infty lim E_{N} ((0, s / (4 N)); l) = e^{- s /4} ⟨ e^{(s^{1/2} /2) Tr (U + U^{†})} ⟩_{U (l)} .

j = 1 \sum n x_{j}^{2} \frac{\partial ^{2}}{\partial x _{j}^{2}} + \frac{2}{α} 1 \leq j < k \leq n \sum \frac{1}{x _{j} - x _{k}} (x_{j}^{2} \frac{\partial}{\partial x _{j}} - x_{k}^{2} \frac{\partial}{\partial x _{k}}) .

j = 1 \sum n x_{j}^{2} \frac{\partial ^{2}}{\partial x _{j}^{2}} + \frac{2}{α} 1 \leq j < k \leq n \sum \frac{1}{x _{j} - x _{k}} (x_{j}^{2} \frac{\partial}{\partial x _{j}} - x_{k}^{2} \frac{\partial}{\partial x _{k}}) .

[u]_{κ}^{(α)} := j = 1 \prod n \frac{Γ ( u - ( j - 1 ) / α + κ _{j} )}{Γ ( u - ( j - 1 ) / α )};

[u]_{κ}^{(α)} := j = 1 \prod n \frac{Γ ( u - ( j - 1 ) / α + κ _{j} )}{Γ ( u - ( j - 1 ) / α )};

F_{p}^{(} F_{q}^{(α)} (a_{1}, \dots, a_{p}; b_{1}, \dots, b_{q}; x) := κ \sum \frac{α ^{∣ κ ∣}}{h _{κ}^{'}} \frac{[ a _{1} ] _{κ}^{(α)} \dots [ a _{p} ] _{κ}^{(α)}}{[ b _{1} ] _{κ}^{(α)} \dots [ b _{q} ] _{κ}^{(α)}} P_{κ}^{(α)} (x) .

F_{p}^{(} F_{q}^{(α)} (a_{1}, \dots, a_{p}; b_{1}, \dots, b_{q}; x) := κ \sum \frac{α ^{∣ κ ∣}}{h _{κ}^{'}} \frac{[ a _{1} ] _{κ}^{(α)} \dots [ a _{p} ] _{κ}^{(α)}}{[ b _{1} ] _{κ}^{(α)} \dots [ b _{q} ] _{κ}^{(α)}} P_{κ}^{(α)} (x) .

F_{0}^{(} F_{1}^{(α)} (c + (n - 1) / α; (t)^{n}) =

F_{0}^{(} F_{1}^{(α)} (c + (n - 1) / α; (t)^{n}) =

\times 1 \leq j < k \leq n \prod e^{i θ_{k}} - e^{i θ_{j}}^{2/ α} d θ_{1} \dots d θ_{n},

B_{n} (c, α) = j = 1 \prod n \frac{Γ ( 1 + 1/ α ) Γ ( c + ( j - 1 ) / α )}{Γ ( 1 + j / α )} .

B_{n} (c, α) = j = 1 \prod n \frac{Γ ( 1 + 1/ α ) Γ ( c + ( j - 1 ) / α )}{Γ ( 1 + j / α )} .

⟨ e^{s^{1/2} Tr (U + U^{†})} ⟩_{U (l)} = F_{0}^{(} F_{1}^{(α)} (l; (s)^{l})_{α = 1};

⟨ e^{s^{1/2} Tr (U + U^{†})} ⟩_{U (l)} = F_{0}^{(} F_{1}^{(α)} (l; (s)^{l})_{α = 1};

⟨ (det U)^{q} e^{s^{1/2} Tr (U + U^{†})} ⟩_{U (l)} = \frac{s ^{q l /2}}{l ! B _{l} ( q + 1 , α )} F_{0}^{(} F_{1}^{(α)} (q + l; (s)^{l})_{α = 1};

⟨ (det U)^{q} e^{s^{1/2} Tr (U + U^{†})} ⟩_{U (l)} = \frac{s ^{q l /2}}{l ! B _{l} ( q + 1 , α )} F_{0}^{(} F_{1}^{(α)} (q + l; (s)^{l})_{α = 1};

\frac{1}{( 2 π ) ^{n} C _{β, n}} 1 \leq j < k \leq n \prod e^{i θ_{k}} - e^{i θ_{j}}^{β}, C_{β, n} = \frac{Γ ( β n /2 + 1 )}{Γ ( β /2 + 1 ) ^{n}};

\frac{1}{( 2 π ) ^{n} C _{β, n}} 1 \leq j < k \leq n \prod e^{i θ_{k}} - e^{i θ_{j}}^{β}, C_{β, n} = \frac{Γ ( β n /2 + 1 )}{Γ ( β /2 + 1 ) ^{n}};

⟨ (det U)^{q} e^{s^{1/2} Tr (U + U^{†})} ⟩_{CE_{β, l}} = \frac{s ^{q l /2}}{C _{β, l} B _{l} ( q + 1 , 2/ β )} F_{0}^{(} F_{1}^{(2/ β)} (q + 1 + (l - 1) β /2; (s)^{l}) .

⟨ (det U)^{q} e^{s^{1/2} Tr (U + U^{†})} ⟩_{CE_{β, l}} = \frac{s ^{q l /2}}{C _{β, l} B _{l} ( q + 1 , 2/ β )} F_{0}^{(} F_{1}^{(2/ β)} (q + 1 + (l - 1) β /2; (s)^{l}) .

l = 1 \prod N x_{l}^{a} e^{- β x_{l} /2} \mathbbm 1_{x_{l} > o} 1 \leq j < k \leq N \prod ∣ x_{k} - x_{j} ∣^{β};

l = 1 \prod N x_{l}^{a} e^{- β x_{l} /2} \mathbbm 1_{x_{l} > o} 1 \leq j < k \leq N \prod ∣ x_{k} - x_{j} ∣^{β};

E_{β}^{hard} ((0, s); a) := N \to \infty lim E_{N, β} ((0, s / (4 N)); a) = e^{- β s} F_{0}^{(} F_{1}^{(β /2)} (2 a / β; (s /4)^{a}) .

E_{β}^{hard} ((0, s); a) := N \to \infty lim E_{N, β} ((0, s / (4 N)); a) = e^{- β s} F_{0}^{(} F_{1}^{(β /2)} (2 a / β; (s /4)^{a}) .

M \to \infty lim F_{1}^{(} F_{1}^{(α)} (a M; b; x / M) = F_{0}^{(} F_{1}^{(α)} (b; x) .

M \to \infty lim F_{1}^{(} F_{1}^{(α)} (a M; b; x / M) = F_{0}^{(} F_{1}^{(α)} (b; x) .

F_{1}^{(} F_{1}^{(2/ β)} ((β /2) (n - 1) + a + 1; β (n - 1) + a + b + 2; (- s)^{n}) .

F_{1}^{(} F_{1}^{(2/ β)} ((β /2) (n - 1) + a + 1; β (n - 1) + a + b + 2; (- s)^{n}) .

= \frac{1}{S _{n} ( a , b , β )} \int_{0}^{1} d x_{1} \dots \int_{0}^{1} d x_{n} l = 1 \prod n x_{l}^{a} (1 - x_{l})^{b} e^{- s x_{l}} 1 \leq j < k \leq n \prod ∣ x_{k} - x_{j} ∣^{β},

H_{p} (x) = \frac{1}{C _{p}^{n}} \int_{0}^{1} d t_{1} \dots \int_{0}^{1} d t_{n} l = 1 \prod n t_{l}^{a} (1 - t_{l})^{b - 1} e^{- x \sum_{l = 1}^{n} t_{l}}

H_{p} (x) = \frac{1}{C _{p}^{n}} \int_{0}^{1} d t_{1} \dots \int_{0}^{1} d t_{n} l = 1 \prod n t_{l}^{a} (1 - t_{l})^{b - 1} e^{- x \sum_{l = 1}^{n} t_{l}}

H_{p} (x) = \times 1 \leq j < k \leq n \prod ∣ t_{k} - t_{j} ∣^{β} e_{p} (1 - t_{1}, \dots, 1 - t_{n}),

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

\FirstPageHeading

\ShortArticleName

Higher-Order Linear Differential Equations for Unitary Matrix Integrals

\ArticleName

Higher-Order Linear Differential Equations

for Unitary Matrix Integrals:

Applications and Generalisations

(with an Appendix by Folkmar Bornemann)

\Author

Peter J. FORRESTER a and Fei WEI b

\AuthorNameForHeading

P.J. Forrester and F. Wei

\Address

a) School of Mathematics and Statistics, The University of Melbourne, Victoria 3010, Australia \EmailD\[email protected]

\Address

b) Department of Mathematics, University of Sussex, Brighton, BN1 9RH, UK \EmailD\[email protected]

\ArticleDates

Received September 04, 2025, in final form January 26, 2026; Published online February 18, 2026

\Abstract

In this paper, we consider characterisations of the class of unitary matrix integrals $\big\langle(\det U)^{q}{\rm e}^{s^{1/2}\operatorname{Tr}(U+U^{\dagger})}\big\rangle_{U(l)}$ in terms of a first-order matrix linear differential equation for a vector function of size $l+1$ , and in terms of a scalar linear differential equation of degree ${l+1}$ . It will be shown that the latter follows from the former. The matrix linear differential equation provides an efficient way to compute the power series expansion of the matrix integrals, which with $q=0$ and $q=l$ are of relevance to the enumeration of longest increasing subsequences for random permutations, and to the question of the moments of the first and second derivative of the Riemann zeta function on the critical line, respectively. This procedure is compared against that following from known characterisations involving the $\sigma$ -Painlevé III*′* second-order nonlinear differential equation. We show too that the natural $\beta$ generalisation of the unitary group integral permits characterisation by the same classes of linear differential equations.

\Keywords

unitary matrix integral; matrix linear differential equation; random permutation; Riemann zeta function; Painlevé equation

\Classification

15B52; 37K10; 05A15

1 Introduction

1.1 A unitary matrix integral and the longest increasing

subsequence problem

The study of this article is the systematic derivation of linear differential equations satisfied by

[TABLE]

the same average with a further factor of $(\det U)^{q}$ , and their $\beta$ generalisation. In (1.1), $U(l)$ refers to the group of $l\times l$ unitary matrices sampled with Haar measure. This matrix integral first appeared in the theoretical physics literature on two-dimensional Yang–Mills gauge theory as a partition function [43]. However, it is not in this context that interest in differential equations for matrix integrals first arose. Rather, it came about due to the realisation that the coefficients $\{T_{l}(N)\}$ say of $z^{2N}/(N!)^{2}$ in the power series of (1.1) (note that this is even in $z$ ) are for a given $l$ an integer family, e.g., for $l=4$ and starting at $N=1$ [44, sequence A047889]

[TABLE]

Thus, according to [46], $T_{l}(N)$ counts the number or permutations of $\{1,2,\dots,N\}$ such that the length of the longest increasing subsequence is less than or equal to $l$ (this number is $N!$ for $N\leq l$ as illustrated in (1.2)). An alternative specification of $T_{l}(N)$ , used in, e.g., [11], is as the number of pairs of standard tableaux of content $N$ with length of the first row less than or equal to $l$ . For an explanation of these notions, and the equivalence of the two descriptions, see, e.g., [34, 47].

In a different guise, a formula equivalent to

[TABLE]

appeared in the earlier work [35]. There it was shown

[TABLE]

where $I_{\nu}(x)={1\over 2\pi}\int_{0}^{2\pi}{\rm e}^{x\cos\theta}{\rm e}^{{\rm i}\nu\theta}{\rm d}\theta$ , $\nu\in\mathbb{Z}$ , is the modified Bessel function of the first kind (i.e., the I-Bessel function). The reason why the Toeplitz determinant in (1.4) is equivalent to the matrix integral in (1.3) follows from the general identities

[TABLE]

Here the first equality follows from Weyl’s result on the explicit functional form of the eigenvalue probability density function (PDF) for matrices from $U(l)$ with Haar measure [21, 54], while the second may be regarded as a circular analogue of Andréief’s identity for the integral over the product of two determinants [26], with use too of the Vandermonde determinant formula to rewrite the product over the differences squared in this form.

Arguments were given in [35] that the generating function in (1.3) is $D$ -finite, meaning that it must satisfy a linear differential equation with polynomial coefficients; for more on this see Appendix B.1. Later, in [11], these differential equations where determined explicitly for $l\leq 7$ . For example, in the case $l=4$ and thus relating to (1.2), one has that (1.3) with $z\mapsto\sqrt{z}$ is a solution of the fifth-order equation [11, last collection of displayed equations]

[TABLE]

(there is also a known linear differential equation of degree four for the power series formed according to the rule that the coefficient of $z^{k}$ is given by the $k$ -th member of (1.2), and the coefficient of $z^{0}$ is unity [14, Section 3.3]). In the present paper, using the theory of Selberg correlation integrals (see [2, 19, 31] and [24, Section 4.6]), we will show how to systematically obtain a first-order $(l+1)\times(l+1)$ matrix linear differential equation, which has as the first component of its vector solution equal to (1.3). As a preliminary step, the matrix average is expressed in terms of a certain Jack polynomial based hypergeometric function of $l$ variables, specialised to the case that all the variables are equal; see Section 2. For a closely related class of hypergeometric functions, earlier work [28] has identified a matrix differential equation characterisation. It is shown in Section 3 how the Jack polynomial based hypergeometric functions of present interest can be deduced as confluent limits of those appearing in [28], allowing for these too to be similarly characterised. This implies a scalar linear differential equation of degree $l+1$ satisfied by (1.3) (Proposition 4.1 below), which with the use of computer algebra can be made explicit. As an illustration, we do this explicitly for $l=8$ , which is the first case beyond the listing in [11]; see Proposition 4.2 below.

Efficient methods for the computation of $\{T_{l}(N)\}$ have come to the fore in the recent works [13, 29] (see too [12]), where a primary aim was to study a scaling regime in the neighbourhood of $l=2\sqrt{N}$ for $N$ large. A precise limit law for the corresponding scaled distribution has been known since [7], and these recent works set about quantifying correction terms. With the numerical tabulations achieved in these works, clear functional forms of corrections were able to exhibited graphically. Such computations were carried out using the known relation for (1.1) and a particular $\sigma$ -Painlevé III*′* equation [32, 53], and also its Chazy-I equivalent form. In Appendix B.2, a comparative discussion of the computational complexity associated with the latter characterisation, and that based on the present matrix differential equation characterisation via the vector recurrence (3.12), is given.

1.2 A unitary matrix integral and the moments of the second derivative

of the Riemann zeta function on the critical line

We have a second main motivation for considering the characterisation of the matrix integral (1.1), and that of modifications of this integral, in terms of a linear differential equation equations satisfied by (1.1). This is the primary role played by the modification of (1.1) and the Toeplitz determinant in (1.4)

[TABLE]

in the recent studies [38, 39]. Here the second equality views the determinant in the Hankel class (i.e., indices a function of $j+k$ ), which is the structure emphasised in these studies. The theme of [38, 39] is the question of statistical properties, particularly moments, of higher-order derivatives of the Riemann zeta function on the critical line. Fundamentally, in keeping with the celebrated work [37], it is hypothesised that these zeta function moments coincide with the scaled limiting moments of the corresponding derivative of the characteristic polynomial of a Haar distributed random unitary matrix; see [3, 6, 8, 9, 18, 20, 27, 36, 55] for earlier work in this direction, and [1, 4, 5, 16, 50] for subsequent developments. In particular, in relation to the second derivative, this leads to the study of [38]

[TABLE]

where

[TABLE]

In [38], it was shown

[TABLE]

where each $f_{m}(x)$ is a combinatorial linear combination of shifted versions of the Hankel determinant appearing on the left-hand side of (1.7), which takes the form

[TABLE]

and where $\mu_{1}+\cdots+\mu_{l}=m$ . Moreover, in [39], it was shown to be given in the form

[TABLE]

Here, each $P_{j}(x)$ is a polynomial defined recursively, and $\tau_{l}(x)$ denotes the Hankel determinant appearing on the left-hand side of (1.7), with $z=2\sqrt{x}$ .

In [39], certain reasoning is given that leads to the conclusion that like (1.1), $\tau_{l}(x)$ satisfies a linear differential equation of degree $l+1$ (see Appendix A for an overview). As mentioned in the opening paragraph, our approach to systematically deriving a linear differential equation that (1.1) satisfies allows for an extra factor $(\det U)^{q}$ to be included for general $q$ , and thus includes (1.7). Thus in this regards our work complements [39]. We consider this further in Section 4.3, and discuss too the advantage of the linear differential equation characterisation over an alternative nonlinear differential equation characterisation put forward in [33].

2 Jack polynomial based hypergeometric function form

It has been revised in the introduction that the matrix integral (1.1) has the interpretations as a partition function for a two-dimensional gauge theory, and as the generating function for counting longest increasing subsequences. There is also an occurrence of (1.1) as a limiting gap probability at the hard edge (neighbourhood of the origin) of a positive definite matrix ensemble with unitary symmetry. To specify this, define complex Wishart matrices $W=X^{\dagger}X$ with $X$ an $n\times N$ $(n\geq N)$ standard complex Gaussian matrix (also referred to as a rectangular GinUE matrix [17]). Let $E_{N}((0,s);l)$ with $l=n-N$ denote the probability of no eigenvalues in $(0,s)$ in this class of complex Wishart ensembles. In [23], it was established that

[TABLE]

(Here the use of the symbol $\beta$ as a subscript on the left-hand side is in the context of the underlying Dyson index; see (2.6) below and subsequent text.) The most significant contribution in [23] for present purposes was the identification of the right-hand side of (2.1) (and generalisations) as a Jack polynomial based generalised hypergeometric function.

To specify this class of multivariable special function, we first recall that, with $\mathbf{x}=(x_{1},\dots,\allowbreak x_{n})$ the Jack polynomials $P_{\kappa}^{(\alpha)}(\mathbf{x})$ are symmetric polynomials of $n$ variables [42, Chapter VI.6], [24, Chapter 12], [41, Chapter 7]. They are dependent on a parameter $\alpha>0$ and labelled by a partition $\kappa:=(\kappa_{1},\dots,\kappa_{n})$ ; in the latter each part $\kappa_{i}\in\mathbb{Z}_{\geq 0}$ is ordered as $\kappa_{1}\geq\kappa_{2}\geq\cdots\geq\kappa_{n}$ . These polynomials are eigenfunctions of the differential operator

[TABLE]

Introducing the dominance partial ordering on partitions of the same length $\mu<\kappa$ , the Jack polynomial $P_{\kappa}^{(\alpha)}(\mathbf{x})$ is uniquely determined as the eigenfunction of (2.2) having the triangular structure with respect to the monomial basis $\{m_{\mu}(\mathbf{x})\}$ , $P_{\kappa}^{(\alpha)}(\mathbf{x})=m_{\kappa}(\mathbf{x})+\sum_{\mu<\kappa}c_{\kappa,\mu}^{(\alpha)}m_{\mu}(\mathbf{x})$ .

In the definition of Jack polynomial based hypergeometric functions, the Jack polynomials play the role of the monomials, although these must be scaled $P_{\kappa}^{(\alpha)}(\mathbf{x})\mapsto\alpha^{|\kappa|}P_{\kappa}^{(\alpha)}(\mathbf{x})/h_{\kappa}^{\prime}$ . Here the precise definition of $h_{\kappa}^{\prime}$ (which depends too on $\alpha$ ) can be found in, e.g., [24, equation (12.46)]. Playing the role of the Pochhammer symbol is the quantity

[TABLE]

see, e.g., [24, equation (12.46)]. With this notation we then have as the series definition of the generalised hypergeometric function in question (see, e.g., [24, Section 13.1])

[TABLE]

In the case $n=1$ , this is independent of $\alpha$ and reduces to its classical scalar counterpart ${}_{p}F_{q}$ .

We are now in a position to specify the result of [23] which allows the right-hand side of (2.1) to be identified as a Jack polynomial based generalised hypergeometric function. Thus, for general ${\alpha>0}$ ,

[TABLE]

(see too [24, equation (13.27)] and [22, equation (4.32)]), where $c\in\mathbb{Z}_{\geq 1}$ (this restriction can be removed by suitably deforming the implied circular contours in (2.3) to Hankel loops [25, proof of Proposition 2]) and

[TABLE]

The notation $(t)^{n}$ on the left-hand side denotes that $\mathbf{x}=(t,t,\dots,t)$ , and thus all $n$ components are equal to $t$ . Making use of (1.5), it follows from (2.3) that

[TABLE]

for a derivation of this result, and its generalisation for $\alpha=1/2$ and $\alpha=2$ , based on the theory of zonal polynomials, see [30]. Extending this, for any $q\in\mathbb{Z}_{\geq 0}$ we read off that

[TABLE]

recall (1.7) in relation to the case $q=l$ .

Specifying the circular $\beta$ ensemble CEβ,n is the eigenvalue PDF

[TABLE]

the case $\beta=2$ corresponds to $U(n)$ . This gives the $\beta$ generalisation of (2.5)

[TABLE]

For the $\beta$ -generalisation of the eigenvalue PDF for complex Wishart matrices, one has the functional form proportional to

[TABLE]

see, for example, [24, equation (3.16) with $a$ appropriately identified]. For this ensemble, let $E_{N,\beta}((0,s);a)$ denote the probability of there being no eigenvalues in $(0,s)$ . Generalising the equality implied by (2.1) and (2.4), it was shown in [23] that for $a\in\mathbb{Z}_{\geq 0}$

[TABLE]

This then relates to (2.7) with $\beta\mapsto 4/\beta$ , $l=a$ , $q=2/\beta-1$ .

3 Matrix differential equation characterisation

As in the classical case, there is the limit relation

[TABLE]

The significance of this is that for $\mathbf{x}=(t)^{n}$ , since the 2022 work [28] we have available a matrix differential equation characterisation of ${\vphantom{F}}_{1}^{\mathstrut}F_{1}^{(\alpha)}(a;b;\mathbf{x})$ . This comes about in the setting of the integral representation [24, Corollary 2.1]

[TABLE]

valid for Re $(a)>-1$ , Re $(b)>-1$ , Re $(\beta)\geq 0$ . Here $S_{n}(a,b,\beta)$ is the Selberg integral (see, e.g., [24, Chapter 4]), which normalises the right-hand side so that it is equal to unity for $s=0$ , and we further remark that the integral is well defined too in a range of negative $\beta$ values, although this detail is not used in our work. The matrix differential equation, which we will first present in the form of a differential-difference recurrence, relates to the family of multi-dimensional integrals

[TABLE]

for $p=0,1,\dots,n$ . The factor $e_{p}(1-t_{1},\dots,1-t_{n})$ in the integrand denotes the $p$ -th elementary symmetric polynomials in the variables $\{(1-t_{i})\}_{i=1}^{n}$ , while $C_{p}^{n}$ is the binomial coefficient $n$ choose $p$ . Setting $p=n$ gives the right-hand side of (3.2) and thus ${}_{1}F^{(2/\beta)}_{1}$ , up to the Selberg integral normalisation, and setting $p=0$ gives the same with $b$ replaced by $b-1$ .

Proposition 3.1 ([28, Corollary 2.5]).

The multi-dimensional integrals $\{H_{p}(x)\}_{p=0}^{n}$ satisfy the differential-difference system

[TABLE]

where

[TABLE]

Corollary 3.2.

Introduce the column vector $\mathbf{H}(x)=[H_{p}(x)]_{p=0}^{n}$ , and the bi-diagonal matrices

[TABLE]

Here $\operatorname{diag}^{+}$ ( $\operatorname{diag}^{-}$ ) denotes the leading diagonal above $($ below $)$ the diagonal itself. We have that the differential-difference equation (3.3) is equivalent to the matrix differential equation

[TABLE]

As a consequence [28, Remark 4.19], introducing the vector of power series solutions

[TABLE]

we have

[TABLE]

and, with $\mathbf{c}_{k}:=[c_{p,k}]_{p=0}^{n}$ , the vector recurrence relation

[TABLE]

where $\mathbf{I}_{n+1}$ denotes the $(n+1)\times(n+1)$ identity.

We now want to make use of (3.1) in Proposition 3.1, in view of both the first and last components of $\mathbf{H}(x)$ relating to ${}_{1}F^{(2/\beta)}_{1}$ .

Proposition 3.3.

Consider the differential-difference system111Here and in the remainder of the paper the functions $\{I_{p}(x)\}$ are not to be confused with the I-Bessel function as appearing in (1.4) and (1.7).

[TABLE]

where $B_{p}$ is as in (3.3). Introduce the column vector ${\mathbf{I}}(x)=[I_{p}(x)]_{p=0}^{n}$ , and the matrices

[TABLE]

We have that the differential-difference equation (3.6) is equivalent to the matrix differential equation

[TABLE]

Furthermore, introducing the vector of power series solutions

[TABLE]

normalised by choosing $d_{0,0}=1$ , we have

[TABLE]

and, provided all the $B_{s}$ are nonzero,

[TABLE]

The coefficient vector $\mathbf{d}_{k}:=[d_{p,k}]_{p=0}^{n}$ satisfies the recurrence relation

[TABLE]

Proof.

To deduce (3.6) and (3.7), we substitute in (3.3) and (3.4), respectively, $x\mapsto-x/a$ , $H_{p}(-x/a)=b^{p}I_{p}(x)$ and take the limit $a\to\infty$ with $a+b$ fixed. The evaluations (3.10) follow from the evaluations of $H_{n}(x)$ , $H_{0}(x)$ in terms of $\vphantom{F}{}_{1}^{\mathstrut}F_{1}^{(2/\beta)}$ noted in the sentence above Proposition 3.1. The evaluation (3.11) can be calculated as the eigenvector of ${\mathbf{Z}}$ corresponding to the eigenvalue 0, normalised to have its first component equal to unity. The vector recurrence relation (3.12) follows by substituting (3.8) in (3.7) and equating like powers of $x$ . ∎

Analogous to the situation with (3.5), since ${\mathbf{d}}_{0}$ is known from (3.11) the vector recurrence (3.14) uniquely determines ${\mathbf{d}}_{1},{\mathbf{d}}_{2},\dots$ provided $(l{\mathbf{I}}_{n+1}-\mathbf{Z})$ is invertible for all $l\in\mathbb{Z}^{+}$ . Since $\mathbf{Z}$ is lower triangular, the criterium for the latter is simply $-B_{p}\notin\mathbb{Z}^{+}$ , for each $p=1,\dots,n$ . Relevant to (2.7) are the choices in Proposition 3.3 of $n=l$ , and $a+b$ such that

[TABLE]

as follows from (3.10). This for $q\in\mathbb{Z}_{\geq 0}$ and $p=0,\dots,l$ is non-negative, allowing then for the applicability of (3.12) written in the form

[TABLE]

and with the initial condition $\mathbf{d}_{0}=[d_{p,0}]_{p=0,\dots,n}$ , where $d_{0,0}=1$ and $d_{p,0}$ ( $p\geq 1$ ) is as given by (3.11). Further specialising to $q=0$ , $\beta=2$ , we can use (3.14) to compute $\{T_{l}(N)\}_{N=1,,\dots,N^{*}}$ (here $N^{*}$ is some given largest value of $N$ ) according to the matrix product formula

[TABLE]

where the notation $(\cdot)_{1}$ denotes the first entry in the corresponding column vector and

[TABLE]

As a concrete example, implementation of (3.15) with $l=4$ , and thus involving products of ${5\times 5}$ matrices, allows for generation of the sequence (1.2).

Remark 3.4.

(1)

Suppose we replace $q$ by $q-1$ in (3.13). Then setting $q=0$ we observe that the column vector $[\delta_{p,l}]_{p=0}^{l}$ is an eigenvector of ${\mathbf{Z}}$ with $n=l$ , allowing for the choice of initial vector corresponding to $d_{p,0}=\delta_{p,l}$ ( $p=0,\dots,l$ ). Using this initial vector, an alternative formula to (3.15) is

[TABLE]

where the notation $(\cdot)_{l+1}$ denotes the $(l+1)$ -th entry (final entry) in the corresponding column vector. In the following section, when referring to $I_{n}(x)$ in the case $q=0$ as it relates to (2.4), it will be assumed that the initial vector has been similarly modified (and thus in (3.10) $d_{n,0}$ is replaced by unity — without this modification (3.10) is invalid since $d_{n,0}$ as specified by (3.11) involves a division by zero).

(2)

For efficient computation using the matrix differential equation, it is advantageous to forgo (3.15) in favour of a recursive computation based on (3.12). See Appendix B.2 for a complexity analysis when using this method to compute $\{T_{l}(N)\}$ in (1.3), which moreover is compared against the formalism used in [13, 29] for the same purpose.

4 Scalar differential equation

4.1 The case of general parameters

Next we would like to relate (3.7) to a scalar linear differential equation of degree $n+1$ for $I_{n}(x)$ (and thus by the replacement $b\mapsto b-1$ for $I_{0}(x)$ ; recall (3.9)).

Proposition 4.1.

Let $I_{n}(x)$ be specified via the differential-difference system (3.6). For some polynomials $\{q_{m}(x)\}_{m=0,1,\dots,n-1}$ , we have

[TABLE]

The degree in $x$ of the coefficient of each ${{\rm d}^{j}\over{\rm d}x^{j}}$ is $[(n-1+j)/2]$ , thus restricting the degrees of the $q_{m}(x)$ . Furthermore, when $\{B_{j}\}$ in (3.6) are all integers, the coefficients of these polynomials are also integers.

Proof.

Let222This notation is not to be confused with $D_{p}$ in (3.3). $D_{x}:=x{{\rm d}\over{\rm d}x}$ . From (3.6) in the case $p=n$ , we read off that

[TABLE]

Using this in the case $p=n-1$ of (3.6) then gives

[TABLE]

Using both (4.2) and (4.3) in the case $p=n-2$ of (3.6) then gives

[TABLE]

We see that the differential operators on the right-hand sides of (4.2), (4.3), (4.4) have the structure

[TABLE]

for certain coefficients.

Repeating this procedure down to the case $p=1$ of (3.6) allows each $I_{j}(x)$ ( $j=0,\dots,n)$ to be written as a differential operator of degree $n-j$ acting on $I_{n}(x)$ . In particular, $n!I_{1}(x)$ has the structure of a differential operator of the form (4.1) with $n\mapsto n-1$ and an extra factor of $x$ multiplying each of the terms involving a derivative. Similarly, $n!I_{0}(x)$ has the structure of a differential operator of the form in (4.1) as written, except for an extra factor of $x$ multiplying each of the terms involving a derivative.

The case $p=0$ of (3.6) gives

[TABLE]

From the established formulas for $I_{0}(x)$ , $I_{1}(x)$ , this implies (4.1). That the coefficients of the polynomials $\{q_{j}(x)\}$ are integers for $\{B_{j}\}$ integers is immediate from the structure of the right-hand sides of (4.2)–(4.4) as extends to $I_{0}(x)$ , $I_{1}(x)$ . ∎

4.2 Characterisation of the matrix integral in (1.3)

In the case $\beta=2$ , $n=l$ and $B_{p}=p(l-p)$ , $I_{n}(x)$ is according to (2.5) and (3.10) (the latter modified in keeping with Remark 3.4.1) equal to the matrix average in (1.3) with $z=x^{1/2}$ and thus is a generating function for the sequence $\{T_{l}(N)\}_{N=1,2,\dots}$ . Using different methods, the precise form of the corresponding differential operator in (4.1) has been computed for $l=2,3,\dots,7$ in [11]. From the substitution method given in the proof of Proposition 4.1, we can reclaim the listing of [11], and extend it to higher values of $l$ (see the work [45] for other applications in random matrix theory of this method in reducing first-order linear matrix differential equations to higher-order scalar differential equations). Here we make note of the case $l=8$ .

Proposition 4.2.

*The matrix integral in (1.3) with $z=x^{1/2}$ and $l=8$ satisfies the differential equation in (4.1) with $n=8$ *(and thus is of degree $9)$ , with

[TABLE]

In (1.3), let us replace $z$ by $x^{1/2}$ so that it reads

[TABLE]

As has just been illustrated in Proposition 4.2, we have from Proposition 4.1 that (4.6) satisfies a certain scalar linear differential equation of degree $l+1$ . As remarked in this context in [11], it is well known and classical that one can translate the latter into a recurrence for $\{T_{l}(N)\}_{N=1,2,\dots}$ with coefficients that are polynomial in $N$ . Previous literature contains the explicit form of the difference equations for small $l$ ; [13, footnote 22] gives references containing the explicit forms up to $l=5$ . As an explicit example, from [52], for $l=4$ we have

[TABLE]

valid for $N=1,2,\dots$ . Iterating this, starting with $N=1,T_{4}(0)=1$ indeed generates the sequence (1.2). Following the proof of Proposition 4.4 below, note is made of the explicit forms for $l=2,3$ and 5.

From [10], one has the conjecture that $\{T_{l}(N)\}_{N=1,2,\dots}$ satisfy a linear difference equation of the form

[TABLE]

with each $P_{k}(N)$ a polynomial in $N$ . However, already with $l=4$ we see from (4.7) that (4.8) is inaccurate, as (4.7) involves only $\{T_{4}(N),T_{4}(N-1),T_{2}(N-2)\}$ whereas (4.8) gives too a term involving $T_{2}(N-3)$ . Making use of Proposition 4.1, we can (by inspection) give the corrected form of (4.8).

Proposition 4.3.

Let $\{T_{l}(N)\}_{N=1,2,\dots}$ be as in (4.6). We have that this sequence satisfies a linear difference equation of the form (4.9) with the upper terminal replaced by $[(l-1)/2]+1$ .

Alternative to first determining the scalar differential equation and then converting it to a difference equation, it is possible to deduce the latter directly from the differential-difference system (3.6).333We thank Folkmar Bornemann for this observation.

Proposition 4.4.

Consider the differential-difference system (3.6). By substituting $I_{p}(x)=\sum_{N=0}^{\infty}d_{p,N}x^{N}$ as is consistent with (3.8), a difference equation for $\{d_{p,N}\}_{N=0,1,\dots}$ of the form

[TABLE]

with each $Q_{k}(p)$ a polynomial in $p$ , can be arrived at.

Proof.

Substituting the power series for $I_{n-1}(x),I_{n}(x)$ in (4.2) and equating coefficients of $x^{N}$ gives $d_{n-1,N}={1\over n}(B_{n}+N)d_{n,N}$ . Carrying out the analogous procedure in (3.6) with $p=n-j$ , $j=1,\dots,n-1$ gives

[TABLE]

After successive substitutions, at the completion of step $j$ , one has an expression for $d_{n-j-1,N}$ in terms of $\{d_{n,N^{\prime}}\}_{N^{\prime}=N-1-[j/2],\dots,N}$ . In relation to $d_{0,N}$ and $d_{1,N}$ there is the additional equation $(N+1)d_{0,N+1}-nd_{1,N}=0$ as follows from (4.5). Substituting the two previously obtained expressions for $d_{0,N}$ , $d_{1,N}$ in terms of $\{d_{n,N^{\prime}}\}$ , we deduce (4.9). ∎

As noted in the first line of this subsection, in the present context, $n=l,B_{p}=p(l-p)$ . Further, with $\{T_{l}(N)\}$ as in (4.9) we have $d_{l,N}=T_{l}(N)/(N!)^{2}$ . In terms of $\{T_{l}(N)\}$ , in low-order cases the procedure of the proof of Proposition 4.4 gives

[TABLE]

Note that the solution of the first recurrence with $T_{2}(0)=1$ is the Catalan number $T_{2}(N)={1\over N+1}\binom{2N}{N}$ ; see [52] for references. The procedure reclaims (4.7) for $l=4$ , and for $l=5$ gives the difference equation (4.9) with upper terminal on the sum $k=3$ and444This corrects the corresponding recurrence in [52, two displayed equations below (18)], which is reported with $P_{3}(N)$ of opposite sign.

[TABLE]

From the explicit forms (4.10), (4.7) and (4.11), one sees an obvious pattern in that the degree of the polynomial coefficients within each case is the same, and equal to $l-1$ . However, the cancellations which take place in going from the difference equation in Proposition 4.4 to the one for $\{T_{l}(N)\}$ makes this difficult to anticipate.

Remark 4.5.

(1)

We know that the matrix differential equation of Proposition 3.3 has a unique solution with $d_{p,0}=\delta_{p,l}(p=0,\dots,l)$ in (3.8); recall Remark 3.4 point 1. As a consequence, it must be that the scalar differential equation (4.1) has a unique power series solution subject only the normalisation $I_{n}(0)=1$ (i.e., that the constant term in the power series is fixed). For this latter property to hold, it must be that the corresponding difference equation (4.9) (corrected according to Proposition 4.3) must have the property that the coefficients of ${T_{l}(-1),T_{l}(-2),\dots,T_{l}(-[(l-1)/2])}$ , whenever they appear, must vanish. Note that this property is illustrated in (4.7) and (4.11).

(2)

In relation to computing the power series solution for general $l$ , there is no advantage in first reducing the matrix differential equation down to a scalar differential equation, rather use of the analogue of (3.15) (which is based on the matrix differential equation) is more efficient. As already mentioned in Remark 3.4, point 2, a quantification of the complexity when using this approach is given in Appendix B.

4.3 Characterisation of the matrix integral in (1.7)

The result of Proposition 3.3 with $\beta=2$ , $n=l$ , $B_{p}=p(2l-p)$ gives, according to (2.5) and (3.10), the matrix integral in (1.7) as proportional to $I_{n}(x)$ . Following the method of Proposition 4.2, we can similarly specify the polynomials in the scalar linear differential equation reduction of the matrix differential equation from Proposition 3.3 for small values of $l$ (without this latter restriction, the coefficients in the corresponding polynomials soon become unwieldy due to the number of digits involved and their increasing degree; recall the listing of Proposition 4.2).

Proposition 4.6.

Consider the matrix integral in (1.7) with $z=x^{1/2}$ , multiplied by $x^{-l^{2}/2}$ and up to a sign factor $(-1)^{\frac{l(l-1)}{2}}$ . For $l=2,\dots,5$ this quantity satisfies the linear differential equations

[TABLE]

respectively.

As revised in Appendix A, the work [39] has given an alternative way to derive the degree ${l+1}$ scalar linear differential equation satisfied by the matrix integral (1.7). In the case $l=2$ , this is made explicit in [39, equation (97)], and consistency with the third-order equation listed above can be checked. Implementing the recursive strategy leading to (A.10) of Appendix A, in the cases $l=3,4,5$ using computer algebra establishes consistency of the results obtained by the alternative method for all the differential equations listed in Proposition 4.6.

Motivating the investigation of [39] in this direction was the realisation that the formalism given in [33] relating to a characterisation in terms of a Painlevé III*′* $\sigma$ -form transcendent cannot be used to generate the full power series about the origin due to the stated boundary condition therein not uniquely determining the coefficients. Thus, considered in [33] was the quantity

[TABLE]

cf. (B.2). Using results from [32] it was noted that $v(s)$ satisfies the particular $\sigma$ -Painlevé III*′* equation

[TABLE]

(cf. (B.3)) subject to the small- $s$ boundary condition

[TABLE]

It was shown in [33] that this allowed for the generation of the first $l$ even powers recursively, as was the need of the application therein. However, missed in [33] was the fact made explicit in [39] that the boundary condition (4.12) does not uniquely determine a solution of (4.12), with there in fact being a family of solutions which require the specification of the coefficient of $s^{2l+1}$ to be made unique beyond degree $s^{2l}$ , and which was incorrectly set to zero (albeit without consequence to the application given there related to the moments of the first derivative of the Riemann zeta function).

The finding in [39] of a scalar differential equation characterisation of the matrix integral (1.7) overcame this latter point, and was used to obtain in particular the coefficients up to degree $s^{4l}$ as required for the application to the moments of the second derivative of the Riemann zeta function on the critical line as outlined in the text of the Introduction including (1.8). The present work complements such a computational objective, with the vector recurrence (3.12) as deduced from the matrix differential equation and analysed from a complexity viewpoint in the second subsection of Appendix B, providing an efficient scheme for this purpose.

Appendix A Alternative derivations from [39]

The purpose of this appendix is to outline the method given in [39] for the derivation of the degree $l+1$ linear differential equation satisfied by (1.7) (see Corollary A.10 below), which at a practical level allows for an alternative derivation of the results of Proposition 4.6. Moreover, the working of [39] will be extended to allow for consideration of the family of matrix averages in (2.5) with $q=l+\alpha$ (for $q\in\mathbb{Z}_{\geq 0}$ ), or more generally the family of Hankel determinants in (1.7) with the index of the $I$ -Bessel function therein generalised from $j+k+1$ to $j+k+\alpha$ (for ${\alpha\in\mathbb{R}_{\geq 0}}$ ). The starting point of the method is the connection between a shifted version of the Hankel determinant appearing in (1.7)

[TABLE]

and the derivative of the Hankel determinant itself. Specifically, we begin with the definition of the determinant of Hankel matrices whose column indices are shifted by Young diagrams.

Definition A.1.

Let $\alpha\in\mathbb{R}_{\geq 0}$ . Let $l\geq 1,h\geq 0$ be integers. Let $Y=(m_{1},\ldots,m_{h})$ be a Young diagram. Define $\tau_{l,Y,\alpha}(x):=\det\!\big[I_{j+k+\alpha+m_{l-k}}\bigl(2\sqrt{x}\bigr)\big]_{j,k=0}^{l-1}$ . Here we adopt the convention that $m_{h+1}=\cdots=m_{l}=0$ whenever $h<l$ . When $Y=\varnothing$ , we simply write $\tau_{l,\varnothing,\alpha}$ as $\tau_{l,\alpha}$ for short.

In particular, when $Y=\varnothing$ and $\alpha=1$ , the above definition reduces to $\tau_{l}(x)$ as given in (A.1). Using the trace of a square matrix, we define the following operations.

Definition A.2.

Let $\alpha\in\mathbb{R}_{\geq 0}$ . Let $l\geq 1$ , $h,n\geq 0$ be integers. Let $Y=(m_{1},\ldots,m_{h})$ be a Young diagram. Define

[TABLE]

and

[TABLE]

where for a matrix $M$ , $\mathrm{adj}(M)$ denotes its adjugate.

In deriving the higher-order linear differential equation for $\tau_{l,\alpha}$ , we mainly use a special form of Young diagram, namely

[TABLE]

of weight $m$ , where $m\geq j\geq 1$ . By applying [39, Theorem 3.1] to general Hankel determinants shifted by $Y_{m,j}$ and specializing to $\tau_{l,Y_{m,j},\alpha}$ , we obtain the following inverse linear system.

Proposition A.3.

Let $l,m\geq 1,m\geq j\geq 1$ . Let $Y_{m,j}$ be given as in (A.2). Then

[TABLE]

where $\mathbf{B}_{m}=(b_{j,k})_{j,k=1,\ldots,m}$ is the $m\times m$ matrix satisfying

[TABLE]

Using the recurrence relations of the modified Bessel function of the first kind,

[TABLE]

we can express $\tau_{l,Y_{1,1},\alpha}$ , $\tau_{l,Y_{2,1},\alpha}$ , and $\tau_{l,Y_{2,2},\alpha}$ in terms of the derivatives of $\tau_{l,\alpha}$ . The following lemmas for general $\alpha$ are obtained by a similar argument to that used in [39, Propositions 4.1 and 4.2] for $\alpha=1$ .

Lemma A.4.

Let $l\geq 1,m\geq 0$ and $Y$ be a Young diagram of length $m$ . Then

[TABLE]

In particular, if $Y=\varnothing$ , then

[TABLE]

Lemma A.5.

Let $l\geq 1$ . Then

[TABLE]

For a general weight $m$ of Young diagrams, we can recursively obtain an expression for $\tau_{l,Y_{m,j},\alpha}$ in terms of the derivatives of $\tau_{l,\alpha}$ . Thus, by identity (A.4), we can rewrite operations involving $T_{h}$ in Proposition A.3 in terms of $S_{h-1}$ and $T_{h-2}$ . Combining Lemmas A.4 and A.5, we then obtain the following formula for $\tau_{l,Y_{m,j},\alpha}$ , which depends recursively on $m$ .

Proposition A.6.

Let $l\geq 1$ and $m\geq 3$ . Let $m\geq j\geq 1$ . Let $Y_{m,j}$ be given as in (A.2). Then

[TABLE]

where $\mathbf{B}_{m}$ is given as in (A.3), $\mathbf{C}^{(1)}_{m}=\big(c_{j,k}^{(1)}\big)_{\begin{subarray}{c}j=1,\ldots,m\\ k=1,\ldots,m-1\end{subarray}}$ is an $m\times(m-1)$ matrix satisfying

[TABLE]

and $\mathbf{C}^{(2)}_{m}=\bigl(c_{j,k}^{(2)}\bigr)_{\begin{subarray}{c}j=1,\ldots,m\\ k=1,\ldots,m-2\end{subarray}}$ is an $m\times(m-2)$ matrix satisfying

[TABLE]

Remark A.7.

The above proposition generalizes [39, Proposition 4.3] from the case $\alpha=1$ to general $\alpha$ . It can be rewritten as a first-order matrix linear differential equation for the vector on the left-hand side analogous to (but distinct from) that implied by Proposition 3.3, with $n=m$ , $\beta=2$ and $B_{p}$ given by (3.13) with $q=m+\alpha$ . On this latter point, note that the entries in the corresponding vectors agree in the final entry but involve different functions otherwise.

Theorem A.8.

Let $l\geq 1$ . Suppose that $\tau_{l,Y_{l+1,l+1},\alpha}$ is obtained recursively by formula (A.8) with initial conditions (A.5), (A.6), and (A.7). Then

[TABLE]

where the coefficients $b_{j,n}^{(l+1)}(\alpha)$ are given recursively. Moreover, we have the following equation:

[TABLE]

Proof.

By induction, it is straightforward to verify that (A.9) holds. It now remains to prove (A.10). On one hand, by Proposition A.6,

[TABLE]

On the other hand, note that

[TABLE]

By the observation $\tau_{l,Y_{l,l},\alpha}(x)=\det\bigl(I_{j+k+\alpha+1}\bigl(2\sqrt{x}\bigr)\bigr)_{j,k=0,\ldots,l-1}$ and the recurrence relation

[TABLE]

we have

[TABLE]

Then by (A.11), we have the identity $\tau_{l,Y_{l+1,l+1},\alpha}\equiv 0$ . Hence, using above and (A.9), we obtain (A.10). ∎

Remark A.9.

Not only $\tau_{l,Y_{l+1,l+1},\alpha}$ , but also $\tau_{l,Y_{m,j},\alpha}$ for any $m,j$ with $l+1\leq j\leq m$ , which can be obtained recursively from formula (A.8) as explained in [39, Remark 3.5], should be [math]. Hence, this will also lead to a linear differential equation of a form similar to (A.10), but of degree $m$ .

As a consequence of Theorem A.8 with $\alpha=1$ , we obtain the following result, which provides a method to derive the degree $l+1$ linear differential equation satisfied by (1.7).

Corollary A.10.

Let $l\geq 1$ . Let $b_{j,n}^{(l+1)}(1)$ be given as in Theorem A.8 with $\alpha=1$ . Then

[TABLE]

Finally, applying Corollary A.10 and verifying with the code,555https://drive.google.com/file/d/1gd4kQO9ABKhALrzeYos-DbeLcelQmik6/view?usp=share_link we re-obtain Proposition 4.6.

Appendix B Computational aspects. Appendix by Folkmar Bornemann∗∗ $*$ Department of Mathematics, Technical University of Munich, 80290 München, Germany††E-mail: [email protected]

B.1 $\boldsymbol{D}$ -finiteness

We discuss the unitary matrix integrals (1.4) and (1.7) given in the form

[TABLE]

which are Toeplitz and Hankel determinants of $f_{n}(z):=I_{n}(2z)$ , where $I_{n}$ denotes the modified Bessel function of the first kind of integer order $n$ . By the symmetries of the unitary matrix integrals, $y_{l}$ and $w_{l}$ are even entire functions. By choosing the normalization $c_{l}=G(2l+1)/G(l+1)^{2}$ (where $G$ is the Barnes $G$ -function), their power series start as666The normalization of $w_{l}$ follows from (2.5) by noting that ${}_{0}^{\mathstrut}F_{1}^{(1)}\bigl(2l;\bigl(z^{2}\bigr)^{l}\bigr)=1+\frac{z^{2}}{2}+\cdots.$

[TABLE]

Bessel’s differential equation $z^{2}f_{n}^{\prime\prime}+zf_{n}^{\prime}-\bigl(4z^{2}+n^{2}\bigr)f_{n}=0$ shows that the $f_{n}$ are $D$ -finite (holonomic), that is, they satisfy a linear differential equation with polynomial coefficients. Since $D$ -finite power series constitute a subalgebra of the formal power series ${\mathbb{C}}\llbracket z\rrbracket$ (see [51, Theorem 2.3]), it is immediately clear that $y_{l}$ , $w_{l}$ are also $D$ -finite. A symmetry argument shows that the holonomic equations for the even functions $y_{l}$ and $w_{l}$ transform to holonomic equations of the same order for $Y_{l}(x):=y_{l}\bigl(\sqrt{x}\bigr)$ and $W_{l}(x):=w_{l}\bigl(\sqrt{x}\bigr)$ – so these entire functions are also $D$ -finite. $D$ -finite power series have $P$ -recursive (holonomic) coefficients (and vice versa, see [51, Theorem 1.5]), that is, the coefficients satisfy a linear recurrence equation with polynomial coefficients. Besides being of structural interest, such a recurrence allows for the effective calculation of the power series coefficients. Though the general theory of holonomic systems provides simple bounds on the order of the differential equations and recurrences – so that they can be constructed by Groebner bases elimination methods, see [56, Section 4] – those bounds are often too large to be of practical use. Since the Bessel differential equation is of order $2$ , the degree of the holonomic differential equation for $y_{l}$ , or $w_{l}$ , is superficially bounded by $2^{l}l!$ because the determinants are sums of $l!$ terms, each of which is the product of $l$ such Bessel functions.

However, by identifying the family $f_{n}(z)$ as a discrete-continuous holonomic system in $(n,z)$ (cf. [56, p. 322]), the order bound can be significantly reduced to the optimal $l+1$ (which is the order given in Proposition 4.1) and the coefficients of the higher-order linear differential equation can be computed by linear algebra. To this end, we supplement the Bessel differential equation with the holonomic relations between $f_{n+1}$ and $f_{n}$ (cf. (A.4)), namely

[TABLE]

Inserting this into the determinants, we get

[TABLE]

clearly we have $\dim_{{\mathbb{C}}(z)}V_{l}\leq l+1$ . A classical result of Siegel [49, Satz 9]777Since the Galois group of the differential equation for the Bessel function $J_{\nu}$ ( $\nu\not\in\frac{1}{2}+{\mathbb{Z}}$ ) is $\text{SL}_{2}({\mathbb{C}})$ (cf. [40, Appendix]), any fundamental system and their derivatives have transcendence degree $\dim_{{\mathbb{C}}}\text{SL}_{2}({\mathbb{C}})=3$ over ${\mathbb{C}}(z)$ . shows that $f_{0}$ , $f_{0}^{\prime}$ are algebraically independent over ${\mathbb{C}}(z)$ , so that $\dim_{{\mathbb{C}}(z)}V_{l}=l+1$ and $z$ , $f_{0}$ , $f_{0}^{\prime}$ can be treated as three independent indeterminates: calculations take place in the polynomial ring ${\mathbb{C}}[z,f_{0},f_{0}^{\prime}]$ . Since Bessel’s equation shows that $V_{l}$ is closed under differentiation, there is a unique set of relatively prime polynomials $a_{j}\in{\mathbb{C}}[z]$ , normalized to monic $a_{l+1}$ , such that

[TABLE]

the same holds for $w_{l}$ . The $a_{j}$ are computed by solving the homogeneous linear system of dimension $(l+1)\times(l+2)$ over ${\mathbb{C}}(z)$ of full rank that is obtained by setting the coefficients of the monomials $f_{0}^{j}(f_{0}^{\prime})^{l-j}$ on the left of (B.1) to zero. Using a computer algebra system, it is thus straightforward to reclaim the differential equations tabulated in Propositions 4.2 and 4.6.

Remark B.1.

In [11], using the order bound $l+1$ , the differential equations, such as (1.6) for ${l=4}$ , were obtained in a far less efficient way by first tabulating, using other means, the first $N$ entries of the power series of $Y_{l}(x)=y_{l}\bigl(\sqrt{x}\bigr)$ for sufficiently large $N$ . Using degree bounds for the coefficient polynomials $a_{j}$ , the presumed differential equation generates a set of linear equations for the coefficients of $a_{j}$ (normalized to monic $a_{l+1}$ ) being compatible with the tabulated values. If $N$ is large enough, and the bounds are chosen appropriately, the linear system has a unique solution. This procedure has been implemented in the Maple package gfun [48]; as an example, the differential equation (1.6) is obtained by the command listtodiffeq after adding three more terms to the list (1.2) as follows:

with(gfun): l := 4: L := [1, 1, 2, 6, 24, 119, 694, 4582, 33324, 261808, 2190688, 19318688, 178108704, ΨΨ 1705985883, 16891621166, 172188608886]: c := [seq(L[k + 1]/(k!)^2, k = 0 .. 15)]: gfun:-Parameters(‘maxordereqn’ = l + 1, ‘maxdegcoeff’ = l): listtodiffeq(c, Y(x));

The recursion (4.7) is thereafter obtained by executing the gfun-command listtorec as follows:

listtorec(L, T(N));

B.2 Complexity of tabulating $\boldsymbol{T_{l}(N)}$

As already noted, $T_{l}(N)=N!$ for $N\leq l$ , so compiling a table of those numbers of combinatorial interest can be restricted to the triangular array $1\leq l\leq N\leq N_{*}$ . Here we compare a tabulation as based on (3.12), to that used in the recent works [13, 29].888For a history of tabulation by combinatorial methods, see [13, p. 916]. A tabulation for $N_{*}=1000$ , computed with the method in [13, Section 3.2], is available for download at https://box-m3.ma.tum.de/f/7c4f8cb22f5d425f8cff/. There, instead of considering (1.1) directly, one first aims to provide a power series expansion of what is essentially the logarithmic derivative (with $z={1\over 2}\sqrt{s}$ ),

[TABLE]

The significance of (B.2) is the characterisation of $u(s)$ as the solution of the second-order second-degree $\sigma$ -Painlevé III*′* equation [32, 53]

[TABLE]

with boundary condition for small $s$

[TABLE]

Writing $u(s)$ as a power series consistent with (B.4),

[TABLE]

and substituting in (B.3) provides a recurrence for successive coefficients. However, since the differential equation (B.3) is cubic as a polynomial expression of $u$ , $u^{\prime}$ , $u^{\prime\prime}$ , the recursion is cubic in the coefficients $c_{n}$ , involving double sums for the evaluation of the Cauchy products. It was therefore suggested in [13, Section 3.2], to use the equivalent third-order first-degree Chazy-I equation, which is obtained from differentiating (B.3) and dividing the result by $2u^{\prime\prime}$

[TABLE]

The differential equation is now quadratic as a polynomial expression in $u$ , $u^{\prime}$ , $u^{\prime\prime}$ , $u^{\prime\prime\prime}$ which yields a much simpler recursive formula for the $c_{n}$ , $n=l+1,\ldots$ – quadratic in the $c_{n}$ with just a single sum for the evaluation of the Cauchy products

[TABLE]

In this way, after integrating and exponentiating the truncated power series of $u(s)$ as required by (B.2), we get the initial coefficients $T_{l}(N)$ in the power series expansion of the average in (1.3).

Counting the operations in ${\mathbb{Q}}$ , and noting that the fast exponentiation of truncated power series with $N_{*}$ terms takes just $O(N_{*}\log N_{*})$ operations with a small implied constant (a method originating with R. Brent, see [15]), computing a single row with a fixed $l$ accounts for

[TABLE]

operations. From this, we get the first row in Table B.2.999Note that the total cost when using the $\sigma$ -Painlevé III*′* equation as in [29] is $O\bigl(N_{*}^{4}\bigr)$ .

Bibliography56

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] Alvarez E., Conrey J.B., Rubinstein M.O., Snaith N.C., Moments of the derivative of the characteristic polynomial of unitary matrices, Random Matrices Theory Appl. 14 (2025), 2550002, 45 pages, ar Xiv:2407.13124 . · doi ↗
2[2] Aomoto K., Jacobi polynomials associated with Selberg integrals, SIAM J. Math. Anal. 18 (1987), 545–549. · doi ↗
3[3] Assiotis T., Bedert B., Gunes M.A., Soor A., On a distinguished family of random variables and Painlevé equations, Probab. Math. Phys. 2 (2021), 613–642, ar Xiv:2009.04760 . · doi ↗
4[4] Assiotis T., Keating J.P., Gunes M.A., Wei F., Exchangeable arrays and integrable systems for characteristic polynomials of random matrices, ar Xiv:2407.19233 .
5[5] Assiotis T., Keating J.P., Gunes M.A., Wei F., Joint moments for characteristic polynomials from the orthogonal and symplectic groups, ar Xiv:2508.09910 .
6[6] Assiotis T., Keating J.P., Warren J., On the joint moments of the characteristic polynomials of random unitary matrices, Int. Math. Res. Not. 2022 (2022), 14564–14603, ar Xiv:2005.13961 . · doi ↗
7[7] Baik J., Deift P., Johansson K., On the distribution of the length of the longest increasing subsequence of random permutations, J. Amer. Math. Soc. 12 (1999), 1119–1178, ar Xiv:math.CO/9810105 . · doi ↗
8[8] Bailey E.C., Bettin S., Blower G., Conrey J.B., Prokhorov A., Rubinstein M.O., Snaith N.C., Mixed moments of characteristic polynomials of random unitary matrices, J. Math. Phys. 60 (2019), 083509, 26 pages, ar Xiv:1901.07479 . · doi ↗

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

1 Introduction

1.1 A unitary matrix integral and the longest increasing

1.2 A unitary matrix integral and the moments of the second derivative

2 Jack polynomial based hypergeometric function form

3 Matrix differential equation characterisation

Proposition 3.1** ([28, Corollary 2.5]).**

Corollary 3.2**.**

Proposition 3.3**.**

Proof.

Remark 3.4**.**

4 Scalar differential equation

4.1 The case of general parameters

Proposition 4.1**.**

Proof.

4.2 Characterisation of the matrix integral in (1.3)

Proposition 4.2**.**

Proposition 4.3**.**

Proposition 4.4**.**

Proof.

Remark 4.5**.**

4.3 Characterisation of the matrix integral in (1.7)

Proposition 4.6**.**

Appendix A Alternative derivations from [39]

Definition A.1**.**

Definition A.2**.**

Proposition A.3**.**

Lemma A.4**.**

Lemma A.5**.**

Proposition A.6**.**

Remark A.7**.**

Theorem A.8**.**

Proof.

Remark A.9**.**

Corollary A.10**.**

Appendix B Computational aspects. Appendix by Folkmar Bornemann∗∗∗*∗Department of Mathematics, Technical University of Munich, 80290 München, Germany††E-mail: [email protected]

B.1 D\boldsymbol{D}D-finiteness

Remark B.1**.**

B.2 Complexity of tabulating Tl(N)\boldsymbol{T_{l}(N)}Tl​(N)

Proposition 3.1 ([28, Corollary 2.5]).

Corollary 3.2.

Proposition 3.3.

Remark 3.4.

Proposition 4.1.

Proposition 4.2.

Proposition 4.3.

Proposition 4.4.

Remark 4.5.

Proposition 4.6.

Definition A.1.

Definition A.2.

Proposition A.3.

Lemma A.4.

Lemma A.5.

Proposition A.6.

Remark A.7.

Theorem A.8.

Remark A.9.

Corollary A.10.

Appendix B Computational aspects. Appendix by Folkmar Bornemann∗∗ $*$ Department of Mathematics, Technical University of Munich, 80290 München, Germany††E-mail: [email protected]

B.1 $\boldsymbol{D}$ -finiteness

Remark B.1.

B.2 Complexity of tabulating $\boldsymbol{T_{l}(N)}$