tcb@breakable

#  $q$-deformed Gaussian unitary ensemble:
spectral moments and genus-type expansions

Sung-Soo Byun  Department of Mathematical Sciences and Research Institute of
Mathematics, Seoul National University, Seoul 151-747, Republic of Korea
<EMAIL_ADDRESS>,  Peter J. Forrester  School of Mathematical and
Statistics, The University of Melbourne, Victoria 3010, Australia
<EMAIL_ADDRESS>and  Jaeseong Oh  June E Huh Center for Mathematical
Challenges, Korea Institute for Advanced Study, 85 Hoegiro, Dongdaemun-gu,
Seoul 02455, Republic of Korea<EMAIL_ADDRESS>
###### Abstract.

The eigenvalue probability density function of the Gaussian unitary ensemble
permits a $q$-extension related to the discrete $q$-Hermite weight and
corresponding $q$-orthogonal polynomials. A combinatorial counting method is
used to specify a positive sum formula for the spectral moments of this model.
The leading two terms of the scaled $1/N^{2}$ genus-type expansion of the
moments are evaluated explicitly in terms of the incomplete beta function.
Knowledge of these functional forms allows for the smoothed leading eigenvalue
density and its first correction to be determined analytically.

###### Key words and phrases:

Discrete orthogonal polynomial ensemble, $q$-Hermite weight, spectral moments,
Flajolet–Viennot theory, genus expansion

###### 2020 Mathematics Subject Classification:

Primary 60B20, 05A15; Secondary 33C20, 05A19, 05A30

Sung-Soo Byun was partially supported by the POSCO TJ Park Foundation (POSCO
Science Fellowship) and by the New Faculty Startup Fund at Seoul National
University. Funding support to Peter Forrester for this research was through
the Australian Research Council Discovery Project grant DP210102887. Jaeseong
Oh was partially supported by a KIAS Individual Grant (HP083401) via the
Center for Mathematical Challenges at Korea Institute for Advanced Study.

##  1\. Introduction

###  1.1. Unitary invariant ensembles and spectral moments

For a given weight $\omega,$ we consider the eigenvalues
$\boldsymbol{x}=\\{x_{j}\\}_{j=1}^{N}$ of an $N\times N$ unitary invariant
random matrix $X$ whose joint eigenvalue probability distribution follows

(1.1)
$d\mathbf{P}(\boldsymbol{x})\propto\prod_{j<k}(x_{j}-x_{k})^{2}\prod_{j=1}^{N}\omega(x_{j})\,dx_{j}.$

Such a random matrix is constructed according to its eigenvalue/eigenvector
decomposition $X=UDU^{\dagger}$, where the diagonal matrix of eigenvalues is
sampled according to (1.1), and $U$ is a Haar distributed random unitary
matrix corresponding to the eigenvectors. The $1$-point function of
$\boldsymbol{x}$ is defined by

(1.2)
$\rho_{N}(x)\equiv\rho_{N,\omega}(x):=N\int_{\mathbb{R}^{N-1}}\mathbf{P}(x,x_{2},\dots,x_{N})\,dx_{2}\,dx_{3}\dots
dx_{N}.$

This can be expressed as

(1.3)
$\rho_{N}(x)=\sum_{j=0}^{N-1}\frac{P_{j}(x)^{2}}{h_{j}}\omega(x)=\frac{P_{N}^{\prime}(x)P_{N-1}(x)-P_{N-1}^{\prime}(x)P_{N}(x)}{h_{N-1}}\omega(x),$

where $\\{P_{j}\\}_{j\geq 0}$ is the sequence of _monic_ orthogonal
polynomials satisfying the orthogonality condition

(1.4) $\int_{\mathbb{R}}P_{j}(x)P_{k}(x)\,\omega(x)\,dx=h_{j}\,\delta_{jk}.$

Here, the second expression of (1.3) follows from the Christoffel-Darboux
formula; see e.g. [25, Ch. 5].

The $p$-th spectral moment is defined by

(1.5)
$m_{N,p}:=\mathbb{E}\Big{[}\,\textup{Tr}X^{p}\Big{]}=\mathbb{E}\Big{[}\sum_{j=1}^{N}x_{j}^{p}\Big{]}=\int_{\mathbb{R}}x^{p}\,\rho_{N}(x)\,dx,$

where the expectation is taken with respect to the measure (1.1). We also
denote

(1.6)
$\mathfrak{m}_{p,j}:=\int_{\mathbb{R}}\frac{x^{p}P_{j}(x)^{2}}{h_{j}}\omega(x)\,dx.$

Then by (1.3), it follows that

(1.7)
$m_{N,p}=\sum_{j=0}^{N-1}\mathfrak{m}_{p,j},\qquad\mathfrak{m}_{p,j}=m_{j+1,p}-m_{j,p}.$

###  1.2. Spectral moments of the GUE

As a prominent example, we discuss the spectral statistics of the Gaussian
unitary ensemble (GUE). Let us write

(1.8) $\omega^{\rm(H)}(x)=\frac{e^{-x^{2}/2}}{\sqrt{2\pi}}$

for the standard Gaussian weight. With this weight, the measure (1.1)
coincides up to a scaling factor with the eigenvalue distribution of the GUE;
see e.g. [25, Definition 1.3.1] for the element-wise construction of such
Hermitian random matrices. The corresponding orthogonal polynomial is given by
the Hermite polynomial

(1.9) $He_{n}(x)=(-1)^{n}e^{x^{2}/2}\frac{d^{n}}{dx^{n}}e^{-x^{2}/2},$

which satisfies

(1.10)
$\int_{\mathbb{R}}He_{n}(x)He_{m}(x)\,\omega^{\rm(H)}(x)\,dx=n!\,\delta_{nm}.$

It satisfies the three-term recurrence relation

(1.11) $He_{n+1}(x)=x\,He_{n}(x)-n\,He_{n-1}(x).$

By (1.3), the $1$-point function $\rho_{N}^{\rm(H)}$ of the GUE is then given
by

(1.12)
$\displaystyle\begin{split}\rho_{N}^{\rm(H)}(x)&=\frac{e^{-x^{2}/2}}{\sqrt{2\pi}}\sum_{j=0}^{N-1}\frac{He_{j}(x)^{2}}{j!}=\frac{e^{-x^{2}/2}}{\sqrt{2\pi}(N-1)!}\Big{(}NHe_{N-1}(x)^{2}-(N-1)He_{N}(x)He_{N-2}(x)\Big{)}.\end{split}$

We also write

(1.13)
$\mathfrak{m}_{2p,j}^{\rm(H)}=\int_{\mathbb{R}}x^{2p}\frac{He_{j}(x)^{2}}{j!}\omega^{\rm(H)}(x)\,dx,\qquad
m_{N,2p}^{\rm(H)}=\sum_{j=0}^{N-1}\mathfrak{m}_{2p,j}^{\rm(H)}$

for the spectral moments of the GUE.

######  Proposition 1.1 (Spectral moments of the GUE, cf. [35, 36, 42, 30]).

The following evaluation formulas hold true:

  * (i)

(Positive sum formula) We have

(1.14) $\displaystyle\mathfrak{m}_{2p,j}^{\rm(H)}$
$\displaystyle=(2p-1)!!\sum_{l=0}^{p}\binom{j}{l}\binom{p}{l}2^{l},$ (1.15)
$\displaystyle m_{N,2p}^{\rm(H)}$
$\displaystyle=(2p-1)!!\sum_{l=0}^{p}\binom{N}{l+1}\binom{p}{l}2^{l}.$

  * (ii)

(Alternating sum formula) We have

(1.16) $\displaystyle\mathfrak{m}_{2p,j}^{\rm(H)}$
$\displaystyle=(2p-1)!!\sum_{r=0}^{p}(-1)^{r}\binom{j+2p-2r}{2p}\binom{p}{r},$
(1.17) $\displaystyle m_{N,2p}^{\rm(H)}$
$\displaystyle=(2p-1)!!\,\sum_{r=0}^{p-1}(-1)^{r}\bigg{[}\binom{N+2p-2r-1}{2p}+\binom{N+2p-2r-2}{2p}\bigg{]}\binom{p-1}{r}.$

######  Remark 1.2.

The formula (1.15) was shown in [35], where the three-term recurrence relation

(1.18)
$(p+1)m_{N,2p}^{\rm(H)}=(4p-2)N\,m_{N,2p-2}^{\rm(H)}+(p-1)(2p-1)(2p-3)\,m_{N,2p-4}^{\rm(H)}$

was also established. One notes too a combinatorial proof of (1.15) due to
Kerov, see [42, Theorems 1 and 2]. In particular, in [42, Lemma 1], the
quantity $\mathfrak{m}_{2p,j}^{\rm(H)}$ was identified with the number of rook
configurations on Ferrer boards; see also [34]. We also mention that the
spectral moment can be expressed as

(1.19)
$\displaystyle\begin{split}m_{N,2p}^{\rm(H)}=(2p-1)!!\,N\,{}_{2}F_{1}(-p,1-N\,;\,2\,;\,2);\end{split}$

see e.g. [22, Eqs.(4.2) and (4.3)] and [36, 68], where

(1.20)
${}_{2}F_{1}(a,b;c;z):=\frac{\Gamma(c)}{\Gamma(a)\Gamma(b)}\sum_{s=0}^{\infty}\frac{\Gamma(a+s)\Gamma(b+s)}{\Gamma(c+s)\
s!}z^{s}$

is the Gauss hypergeometric function.

######  Remark 1.3.

Notice that two different expressions (1.14) and (1.16) indicate the identity

(1.21)
$\sum_{l=0}^{p}\binom{j}{l}\binom{p}{l}2^{l}=\sum_{r=0}^{p}(-1)^{r}\binom{j+2p-2r}{2p}\binom{p}{r}.$

This identity seems to most naturally follow from a double counting; see
Subsection 3.2.

######  Remark 1.4.

The first few values of the spectral moments are given by

(1.22) $\displaystyle\begin{split}&m_{N,0}^{\rm(H)}=N,\qquad
m_{N,2}^{\rm(H)}=N^{2},\qquad m_{N,4}^{\rm(H)}=2N^{3}+N,\\\
&m_{N,6}^{\rm(H)}=5N^{4}+10N,\qquad
m_{N,8}^{\rm(H)}=14N^{5}+70N^{3}+21N.\end{split}$

On the other hand, it follows from the three-term recurrence (1.18) that

(1.23)
$\lim_{N\to\infty}\frac{1}{N^{p+1}}m_{N,2p}^{\rm(H)}=C_{p}:=\frac{1}{p+1}\binom{2p}{p},$

where $C_{p}$ is the $p$-th Catalan number. These limiting moments can be used
to reconstruct the limiting scaled spectral density (see the proof of
Proposition 4.4 below) and give rise to Wigner’s semi-circle law; see e.g.
[3].

######  Remark 1.5 (Genus expansion).

The Stirling numbers of the first kind [58, Section 26.8] are given by

(1.24) $s(n,k)=(-1)^{n-k}\sum_{1\leq b_{1}<\dots<b_{n-k}\leq
n-1}b_{1}b_{2}\dots b_{n-k},\qquad(n>k\geq 1).$

Their generating function is given by

(1.25) $\sum_{k=0}^{n}s(n,k)x^{k}=\frac{\Gamma(x+1)}{\Gamma(x-n+1)};$

see e.g. [58, Eq.(26.8.7)]. By using this and rearranging the terms in the
expression (1.15) of $m_{N,2p}^{\rm(H)}$, we have

$\displaystyle m_{N,2p}^{\rm(H)}$
$\displaystyle=(2p-1)!!\sum_{l=0}^{p}\bigg{(}\sum_{k=0}^{l+1}\frac{s(l+1,k)}{(l+1)!}N^{k}\bigg{)}\binom{p}{l}2^{l}$
$\displaystyle=(2p-1)!!\sum_{r=0}^{p+1}\bigg{(}\sum_{m=0}^{r}\frac{s(p+1-m,p+1-r)}{(p+1-m)!}\binom{p}{m}2^{p-m}\bigg{)}N^{p+1-r}.$

Furthermore, it follows from the identity

(1.26)
$\sum_{m=0}^{r}\frac{s(p+1-m,p+1-r)}{(p+1-m)!}\binom{p}{m}2^{p-m}=0,\qquad\textup{if
}r\textup{ is odd}$

that we have the $1/N^{2}$ expansion

(1.27)
$m_{N,2p}^{\rm(H)}=\sum_{g=0}^{\lfloor(p+1)/2\rfloor}\mathcal{E}_{g}(p)N^{p+1-2g},$

where

(1.28)
$\mathcal{E}_{g}(p)=(2p-1)!!\sum_{m=0}^{2g}\frac{s(p+1-m,p+1-2g)}{(p+1-m)!}\binom{p}{m}2^{p-m}.$

The coefficient $\mathcal{E}_{g}(p)$ corresponds to the count of pairings
among the edges of a $2p$-gon, which is dual to a map on a compact Riemann
surface of genus $g$; see [35]. Therefore, the expansion (1.26) is commonly
referred to as the genus expansion. In particular, the first few values of
$\mathcal{E}_{g}$ are given as

(1.29) $\displaystyle\mathcal{E}_{0}(p)$
$\displaystyle=C_{p},\qquad\mathcal{E}_{1}(p)=C_{p}\,\frac{(p+1)!}{(p-2)!}\frac{1}{12},\qquad\mathcal{E}_{2}(p)=C_{p}\,\frac{(p+1)!}{(p-4)!}\frac{5p-2}{1440},$
(1.30) $\displaystyle\mathcal{E}_{3}(p)$
$\displaystyle=C_{p}\,\frac{(p+1)!}{(p-6)!}\frac{35p^{2}-77p+12}{362880},\qquad\mathcal{E}_{4}(p)=C_{p}\,\frac{(p+1)!}{(p-8)!}\frac{175p^{3}-945p^{2}+1094p-72}{87091200}.$

The small order cases also follow from [68, Theorem 7], which were based on a
large $N$ analysis of (1.18). The form of $\mathcal{E}_{g}(p)$ given by (1.28)
appears to be new. We also mention that $\mathcal{E}_{g}(p)$ satisfies the two
parameter recurrence relation [35]

(1.31)
$(p+2)\mathcal{E}_{g}(p+1)=2(2p+1)\mathcal{E}_{g}(p)+p(2p+1)(2p-1)\mathcal{E}_{g-1}(p-1),$

which is also known as the topological recursion.

###  1.3. Discrete $q$-deformed GUE

We now discuss the discrete $q$-formed version of the GUE. For this purpose,
we shall frequently use the notations

(1.32)
$\displaystyle(a;q)_{n}:=(1-a)(1-aq)\dots(1-aq^{n-1}),\qquad(a_{1},\dots,a_{k};q)_{\infty}:=\prod_{l=0}^{\infty}(1-a_{1}q^{l})\dots(1-a_{k}q^{l}).$

The $q$-hypergeometric series ${}_{r}\phi_{s}$ is defined by

(1.33) ${}_{r}\phi_{s}\left(\begin{matrix}a_{1},\dots,a_{r}\vskip 3.0pt plus
1.0pt minus 1.0pt\\\
b_{1},\dots,b_{s}\end{matrix}\,;\,q,z\right)=\sum_{k=0}^{\infty}\frac{(a_{1};q)_{k}\dots(a_{r};q)_{k}}{(b_{1};q)_{k}\dots(b_{s};q)_{k}}\frac{z^{k}}{(q;q)_{k}}\Big{(}(-1)^{k}q^{k(k-1)/2}\Big{)}^{1+s-r};$

see e.g. [58, Chapter 17]. Note that if one of the parameters
$a_{1},\dots,a_{r}$ is a non-positive power of $q$, then the sum terminates.
We also mention that as $q\to 1$, it recovers the more familiar ${}_{r}F_{s}$
hypergeometric function:

(1.34) ${}_{r}\phi_{s}\left(\begin{matrix}q^{c_{1}},\dots,q^{c_{r}}\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
q^{d_{1}},\dots,q^{d_{s}}\end{matrix}\,;\,q,(q-1)^{1+s-r}\,z\right)\to{}_{r}F_{s}\left(\genfrac{.}{.}{0.0pt}{}{c_{1}{,}\mspace{8.0mu}\dots,\mspace{8.0mu}c_{r}}{d_{1}{,}\mspace{8.0mu}\dots,\mspace{8.0mu}d_{s}}\;\middle|\;z\right):=\sum_{k=0}^{\infty}\frac{(c_{1})_{k}\dots(c_{r})_{k}}{(d_{1})_{k}\dots(d_{s})_{k}}\frac{z^{k}}{k!}.$

Let us recall that the $q$-derivative $D_{q}$ is defined by

(1.35) $D_{q}f(x)=\frac{f(x)-f(qx)}{(1-q)x}.$

Recall also that its anti-derivative, the Jackson $q$-integral, is defined by

(1.36)
$\displaystyle\int_{0}^{\infty}f(x)\,d_{q}x=(1-q)\sum_{k=-\infty}^{\infty}q^{k}\,f(q^{k}),\qquad\int_{0}^{\alpha}f(x)\,d_{q}x=(1-q)\sum_{k=0}^{\infty}\alpha
q^{k}\,f(\alpha q^{k}),$

and

(1.37)
$\int_{a}^{b}f(x)\,d_{q}x=\int_{0}^{b}f(x)\,d_{q}x-\int_{0}^{a}f(x)\,d_{q}x.$

Note too that

(1.38)
$\int_{-1}^{1}f(x)\,d_{q}x=\frac{1}{a}\int_{-a}^{a}f\Big{(}\frac{x}{a}\Big{)}\,d_{q}x.$

The discrete $q$-Hermite weight and its scaled version are defined by

(1.39)
$\omega^{\rm(dH)}(x):=\frac{(qx,-qx;q)_{\infty}}{(q,-1,-q;q)_{\infty}},\qquad\widehat{\omega}^{\rm(dH)}(x):=\sqrt{1-q}\,\omega^{{\rm(dH)}}(\sqrt{1-q}\,x).$

The support of $\omega^{\rm(dH)}(x)$ is taken as the $q$-lattice $\pm q^{k}$,
$k=0,1,\dots$ (note that for these values of $x$, $\omega^{\rm(dH)}(x)>0$),
which in the continuum limit $q\to 1^{-}$ is the interval $[-1,1]$. In this
same limit, $\widehat{\omega}^{\rm(dH)}(x)$ has support of all the real line,
and moreover it recovers the Gaussian weight

(1.40) $\lim_{q\to
1^{-}}\widehat{\omega}^{\rm(dH)}(x)=\frac{1}{\sqrt{2\pi}}e^{-x^{2}/2};$

see e.g. [1]. For $\omega^{\rm(dH)}(x)$, one notes that the even moments are
given by

(1.41) $\int_{-1}^{1}x^{2p}\,\omega^{\rm(dH)}(x)\,d_{q}x=(1-q)(q;q^{2})_{p};$

see e.g. [55], whereas the odd moments vanish. The moments of
$\widehat{\omega}^{\rm(dH)}(x)$ follow from (1.41) by applying (1.38).

We shall use several different conventions for the $q$-Hermite polynomials.

  * •

The $q$-Hermite polynomial $H_{n}(x;q)$ is defined by

(1.42)
$\displaystyle\begin{split}H_{n}(x;q)&:=x^{n}\,{}_{2}\phi_{0}\left(\begin{matrix}q^{-n},q^{-n+1}\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
-\end{matrix}\,;\,q^{2},\frac{q^{2n-1}}{x^{2}}\right)\\\
&=\sum_{k=0}^{n/2}(-1)^{k}\,\frac{(q^{-n};q^{2})_{k}(q^{-n+1};q^{2})_{k}q^{k(2n-k)}}{(q^{2};q^{2})_{k}}x^{n-2k};\end{split}$

see [43, Section 14.28]. It satisfies the orthogonality

(1.43)
$\int_{-1}^{1}H_{m}(x;q)H_{n}(x;q)\omega^{\rm(dH)}(x)\,d_{q}x=(1-q)(q;q)_{n}q^{n(n-1)/2}\delta_{n,m}$

and the three-term recurrence relation

(1.44) $x\,H_{n}(x;q)=H_{n+1}(x;q)+q^{n-1}(1-q^{n})H_{n-1}(x;q).$

  * •

We also set

(1.45)
$\begin{split}\widehat{H}_{n}(x;q)&:=\frac{1}{(1-q)^{n/2}}H_{n}\Big{(}\sqrt{1-q}\,x;q\Big{)}\\\
&=\sum_{k=0}^{n/2}(-1)^{k}\,\frac{(q^{-n};q^{2})_{k}(q^{-n+1};q^{2})_{k}q^{k(2n-k)}}{(q^{2};q^{2})_{k}}\frac{1}{(1-q)^{k}}x^{n-2k}.\end{split}$

Then it follows that

(1.46)
$x\,\widehat{H}_{n}(x;q)=\widehat{H}_{n+1}(x;q)+q^{n-1}[n]_{q}\,\widehat{H}_{n-1}(x;q),$

where we used the convention that the $q$-analogue of an integer $n$ given by

$[n]_{q}:=\dfrac{1-q^{n}}{1-q}.$

Note that in the continuum limit, it recovers the classical Hermite polynomial

(1.47) $\lim_{q\to 1^{-}}\widehat{H}_{n}(x;q)=He_{n}(x).$

######  Remark 1.6.

In [38], the $q$-Hermite polynomial $\widetilde{H}_{n}(x;q)$ that satisfies
the recurrence relation

(1.48)
$x\,\widetilde{H}_{n}(x;q)=\widetilde{H}_{n+1}(x;q)+[n]_{q}\,\widetilde{H}_{n-1}(x;q)$

has been studied. We emphasize that while $H_{n}$ and $\widehat{H}_{n}$ are
related by simple scaling, $\widetilde{H}_{n}$ is not.

According to (1.3) the density $\rho_{N}^{\rm(dH)}$ of the discrete $q$-GUE is
given by

(1.49)
$\rho_{N}^{\rm(dH)}(x)=\frac{1}{1-q}\sum_{j=0}^{N-1}\frac{H_{j}(x,q)^{2}}{(q;q)_{j}q^{j(j-1)/2}}\omega^{\rm(dH)}(x).$

Similarly, we define

(1.50)
$\widehat{\rho}_{N}^{\rm(dH)}(x):=\sqrt{1-q}\,\rho_{N}^{\rm(dH)}(\sqrt{1-q}\,x).$

Let us define the corresponding spectral moments

(1.51) $\displaystyle m_{N,2p}^{\rm(dH)}$
$\displaystyle:=\int_{-1}^{1}|x|^{2p}\rho_{N}^{\rm(dH)}(x)\,d_{q}x,$ (1.52)
$\displaystyle\widehat{m}_{N,2p}^{\rm(dH)}$
$\displaystyle:=\int_{-1/\sqrt{1-q}}^{1/\sqrt{1-q}}|x|^{2p}\,\widehat{\rho}_{N}^{\rm(dH)}(x)\,d_{q}x=\frac{1}{(1-q)^{p}}\,m_{N,2p}^{\rm(dH)}.$

We also denote

(1.53)
$\mathfrak{m}_{2p,j}^{\rm(H)}=\int_{-1}^{1}\frac{|x|^{2p}H_{j}(x,q)^{2}\omega^{\rm(dH)}(x)}{(1-q)\,(q;q)_{j}q^{j(j-1)/2}}\,d_{q}x,\qquad\widehat{\mathfrak{m}}_{2p,j}^{\rm(H)}=\frac{1}{(1-q)^{p}}\mathfrak{m}_{2p,j}^{\rm(H)}$

so that

(1.54)
$m_{N,2p}^{\rm(dH)}=\sum_{j=0}^{N-1}\mathfrak{m}_{2p,j}^{\rm(H)},\qquad\widehat{m}_{N,2p}^{\rm(dH)}=\sum_{j=0}^{N-1}\widehat{\mathfrak{m}}_{2p,j}^{\rm(H)}.$

###  1.4. Objectives and previous literature

A celebrated result in the application of random matrices is that the
$1/N^{2}$ expansion of the GUE moments (1.27), albeit with the coefficients
defined implicitly via the two-variable recurrence relation (1.31), counts the
number of pairings of the sides of a $2p$-gon which are dual to a map on a
compact orientable Riemann surface of genus $g$ (equivalently, a surface of
genus $g$ results when the sides of the polygon are identified in pairs). This
line of work was initiated by Brézin et al. [11], with the recursive
specification of the coefficients obtained in [35]. As noted below (1.23), the
leading term in the expansion (1.27) was isolated much earlier in the
development of random matrix theory, specifically as a way of computing the
limiting spectral density [67].

In the ensuing years eigenvalue probability distributions of the general form
(1.1) have occurred in various applications, notably in matrix models of low-
dimensional quantum field theories [49], fermionic systems [23], quantum
transport [7], number theory [41], wireless communication [63] and tiling
models [8]. Consideration of observables in these settings has motivated a
number of studies into the corresponding spectral moments; references include
[59, 20, 37, 50, 66, 53, 32, 61, 12, 40, 4, 56, 6, 28, 51, 52].

An explanation of the three-term recurrence (1.18) in terms of differential
properties of the Laplace transform of $\rho_{N}^{\rm(H)}$, first uncovered in
[44] (see also [36]), has given motivation to the study of spectral moments in
other random matrix ensembles exhibiting this and related mathematical
structures [46, 68, 22, 5, 60, 31, 13, 14, 21, 33]. Considerations along these
lines relate to integrable properties of the corresponding ensembles.
Specifically, in [54] the ensemble (1.1) with discrete $q$-Hermite weight
(1.39) was isolated as possessing distinguished integrable properties, and
which furthermore has the status of being a generalisation of the GUE natural
from the viewpoint of $q$-orthogonal polynomial theory. The subsequent works
[30] and [16] furthered the study of the integrable properties of this model,
building on [17, 45]. We also refer to [9, 10, 29, 39, 47, 57] for related
studies on discrete orthogonal polynomial ensembles.

The aim of the present work is to study integrable properties of spectral
moments and the eigenvalue density of the discrete $q$-Hermite unitary
invariant random matrix ensemble beyond the results of the earlier works [54,
30, 16]. The details of our findings are specified in the next section. They
include a positive sum formula for the spectral moments (2.1), the leading two
terms of the scaled $1/N^{2}$ expansion of this formula (Theorem 2.3), and the
form of the smoothed leading eigenvalue density and its leading correction
(Theorem 2.6).

##  2\. Main results

###  2.1. Spectral moments

As a $q$-analogue of Proposition 1.1, we have the following.

######  Theorem 2.1 (Spectral moments of the $q$-deformed GUE).

The following hold true:

  * (i)

(Positive sum formula) We have

(2.1)
$\displaystyle\begin{split}\widehat{\mathfrak{m}}_{2p,j}^{\rm(H)}&=\sum_{l=0}^{p}q^{(j-l)(2p-l)+l(l-1)/2}\genfrac{[}{]}{0.0pt}{}{j}{l}_{q}\dfrac{[2p]_{q}!}{[2p-2l]_{q}!!\,[l]_{q}!},\end{split}$
(2.2) $\displaystyle\widehat{m}_{N,2p}^{\rm(dH)}$
$\displaystyle=\sum_{j=0}^{N-1}\sum_{l=0}^{p}q^{(j-l)(2p-l)+l(l-1)/2}\genfrac{[}{]}{0.0pt}{}{j}{l}_{q}\dfrac{[2p]_{q}!}{[2p-2l]_{q}!!\,[l]_{q}!}.$

Taking $q\to 1$, we reclaim (1.14) and (1.15).

  * (ii)

(Alternating sum formula) We have

(2.3)
$\displaystyle\begin{split}\widehat{\mathfrak{m}}_{2p,j}^{\rm(dH)}&=[2p-1]_{q}!!\sum_{r=0}^{p}(-1)^{r}q^{r(r-1)}\genfrac{[}{]}{0.0pt}{}{j+2p-2r}{2p}_{q}\genfrac{[}{]}{0.0pt}{}{p}{r}_{q^{2}},\end{split}$
(2.4) $\displaystyle\widehat{m}_{N,2p}^{\rm(dH)}$
$\displaystyle=[2p-1]_{q}!!\sum_{r=0}^{p-1}(-1)^{r}q^{r(r+1)}\bigg{(}\genfrac{[}{]}{0.0pt}{}{N+2p-2r-1}{2p}_{q}+\genfrac{[}{]}{0.0pt}{}{N+2p-2r-2}{2p}_{q}\bigg{)}\genfrac{[}{]}{0.0pt}{}{p-1}{r}_{q^{2}}.$

Taking $q\to 1$, we reclaim (1.16) and (1.17).

The positive formula (2.1) is new to our knowledge. The alternating formula
(2.4) was shown in [64] and [30, Eq.(2.39)]. Then (2.3) follows from

(2.5)
$\widehat{\mathfrak{m}}_{2p,j}^{\rm(dH)}=\widehat{m}_{j+1,2p}^{\rm(dH)}-\widehat{m}_{j,2p}^{\rm(dH)}.$

A linear recurrence for the spectral moment $\widehat{m}_{N,2p}^{\rm(dH)}$
extending (1.18) has not been discovered for general $q\in(0,1]$. However it
was found by Cohen [16, Theorem 4.4] that a certain weighted sum over $N$ of
$\widehat{\mathfrak{m}}_{2p,j}^{\rm(dH)}$ can be evaluated in terms of a
particular $q$-Hahn polynomial, from which a three-term recurrence of such
summed moments follows.

We mention that contrary to the $q=1$ case where one can apply Pascal’s
triangle identity consecutively, it seems not straightforward to simplify the
summation

$\displaystyle\sum_{j=0}^{N-1}q^{j(2p-l)}\genfrac{[}{]}{0.0pt}{}{j}{l}_{q}=\sum_{j=0}^{N-1}q^{j(2p-l)}\frac{(1-q^{j})(1-q^{j-1})\dots(1-q^{j-l+1})}{(1-q)(1-q^{2})\dots(1-q^{l})}$

in the expression (2.2).

######  Remark 2.2 (Other forms of alternating formulas).

In [30, Proposition 5.3], the following alternative formula has also been
derived

(2.6)
$\displaystyle\begin{split}\widehat{\mathfrak{m}}_{2p,j}^{\rm(dH)}=\sum_{k=0}^{p}\sum_{l=0}^{k}(-1)^{k}q^{k^{2}+2kj+l(4l-4k-2j-1)}(1-q)^{2l-p}\genfrac{[}{]}{0.0pt}{}{p}{k}_{q^{2}}\genfrac{[}{]}{0.0pt}{}{k}{k-l}_{q^{2}}^{2}\frac{[j]_{q}!}{[j-2l]_{q}!}.\end{split}$

It was shown in [16, Theorem 4.3] that the spectral moment
$\widehat{m}_{N,2p}^{\rm(dH)}$ can be expressed in terms of the hypergeometric
function as

(2.7)
$\displaystyle\begin{split}\widehat{m}_{N,2p}^{\rm(dH)}&=[2p-1]_{q}!!\,q^{-\frac{1}{2}(p^{2}+p)}\frac{(-q;q)_{p}}{(1-q^{p})(q;q)_{p}}{}_{3}\phi_{2}\left(\begin{matrix}-q^{p+1},q^{p},q^{-p}\vskip
3.0pt plus 1.0pt minus 1.0pt\\\ -q,q^{p+1}\end{matrix}\,;\,q,q\right)\\\
&\quad-[2p-1]_{q}!!\,q^{Np-\frac{1}{2}(p^{2}+p)}\frac{(-q;q)_{p}}{(1-q^{p})(q;q)_{p}}{}_{3}\phi_{2}\left(\begin{matrix}-q^{p+1},q^{p},q^{-p}\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
-q,q^{p+1}\end{matrix}\,;\,q,q^{N+1}\right).\end{split}$

Then by (1.33), we have

(2.8) $\displaystyle\widehat{m}_{N,2p}^{\rm(dH)}$
$\displaystyle=[2p-1]_{q}!!\,q^{-\frac{1}{2}(p^{2}+p)}\frac{(-q;q)_{p}}{(1-q^{p})(q;q)_{p}}\sum_{k=0}^{p}\frac{(-q^{p+1};q)_{k}\,(q^{p},q)_{k}\,(q^{-p};q)_{k}}{(-q;q)_{k}\,(q^{p+1};q)_{k}}\frac{q^{k}}{(q;q)_{k}}\Big{(}1-q^{(k+p)N}\Big{)},$
(2.9) $\displaystyle\widehat{\mathfrak{m}}_{2p,j}^{\rm(dH)}$
$\displaystyle=[2p-1]_{q}!!\,q^{-\frac{1}{2}(p^{2}+p)}\frac{(-q;q)_{p}}{(1-q^{p})(q;q)_{p}}\sum_{k=0}^{p}\frac{(-q^{p+1};q)_{k}\,(q^{p},q)_{k}\,(q^{-p};q)_{k}}{(-q;q)_{k}\,(q^{p+1};q)_{k}}\frac{q^{k}}{(q;q)_{k}}q^{(k+p)j}(1-q^{k+p}).$

###  2.2. Scaled moments and density on the exponential lattices

Following the scaling for the Stieltjes-Wigert random matrix ensemble [26], we
set

(2.10) $q=e^{-\lambda/N},\qquad\lambda\geq 0.$

Recall that the regularised incomplete beta function is defined by

(2.11)
$I_{x}(a,b)=\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\int_{0}^{x}t^{a-1}(1-t)^{b-1}\,dt.$

######  Theorem 2.3 (Genus expansion of the spectral moments).

Let $q$ be scaled as in (2.10). Then for all positive integers $p$ we have

(2.12)
$q^{p}m_{N,2p}^{\rm(dH)}=\Big{(}q(1-q)\Big{)}^{p}\widehat{m}_{N,2p}^{\rm(dH)}=\mathcal{M}_{2p,0}^{\rm(dH)}N+\frac{\mathcal{M}_{2p,1}^{\rm(dH)}}{N}+O(\frac{1}{N^{3}}),$

as $N\to\infty$, where

(2.13) $\displaystyle\mathcal{M}_{2p,0}^{\rm(dH)}$
$\displaystyle=\frac{1}{\lambda\,p}\,I_{1-e^{-\lambda}}(p+1,p),$ (2.14)
$\displaystyle\mathcal{M}_{2p,1}^{\rm(dH)}$
$\displaystyle=-\frac{\lambda\,p}{6}\bigg{(}I_{1-e^{-\lambda}}(p+1,p)+\frac{(2p-1)!}{p!(p-1)!}e^{-p\lambda}(1-e^{-\lambda})^{p-1}\Big{(}2+p-(2p+1)e^{-\lambda}\Big{)}\bigg{)}.$

(a) $(p,\lambda)\to\mathcal{M}_{2p,0}^{\rm(dH)}$

(b) $(p,\lambda)\to\mathcal{M}_{2p,1}^{\rm(dH)}$

(c) $p\to\mathcal{M}_{2p,1}^{\rm(dH)}$ with $\lambda=2$

(d) $\lambda\to\mathcal{M}_{2p,1}^{\rm(dH)}$ with $p=3$

Figure 1. The plots (A) and (B) show the graphs of
$\mathcal{M}_{2p,0}^{\rm(dH)}$ and $\mathcal{M}_{2p,1}^{\rm(dH)}$,
respectively. The plot (C) shows graph of $p\to\mathcal{M}_{2p,1}^{\rm(dH)}$
with $\lambda=2$ and its comparison with
$(q^{p}m_{N,2p}^{\rm(dH)}-\mathcal{M}_{2p,0}^{\rm(dH)}N)\,N$ for $N=10,20$,
and $N=40$ (green, red and blue dots, respectively). The plot (D) shows a
similar graph graph of $\lambda\to\mathcal{M}_{2p,1}^{\rm(dH)}$ with $p=3$.

######  Remark 2.4 (Continuum limit $q\to 1$).

As a consequence of Theorem 2.3, it follows that

(2.15) $\displaystyle\widehat{m}_{N,2p}^{\rm(dH)}$
$\displaystyle=\Big{(}\frac{N}{\lambda}\Big{)}^{p}\bigg{(}\mathcal{M}_{2p,0}^{\rm(dH)}N+\frac{3p}{2}\lambda\,\mathcal{M}_{2p,0}^{\rm(dH)}+\Big{(}\frac{(27p-1)p}{24}\lambda^{2}\,\mathcal{M}_{2p,0}^{\rm(dH)}+\mathcal{M}_{2p,1}^{\rm(dH)}\Big{)}\frac{1}{N}+O(\frac{1}{N^{2}})\bigg{)}.$

Then one can see that

(2.16) $\lim_{\lambda\to
0}\frac{1}{\lambda^{p}}\mathcal{M}_{2p,0}^{\rm(dH)}=\mathcal{E}_{0}(p)$

and

(2.17) $\lim_{\lambda\to
0}\frac{1}{\lambda^{p}}\Big{(}\frac{(27p-1)p}{24}\lambda^{2}\,\mathcal{M}_{2p,0}^{\rm(dH)}+\mathcal{M}_{2p,1}^{\rm(dH)}\Big{)}=\frac{(2p-1)!}{6(p-1)!(p-2)!}=\mathcal{E}_{1}(p).$

Thus the first two terms in the expansion (1.27) are reclaimed.

######  Remark 2.5.

In [30, Eq.(2.45)], the leading scaled moment is also derived using the
alternating formula, which takes the form

(2.18)
$\mathcal{M}_{2p,0}^{\rm(dH)}=\frac{1}{\lambda}\bigg{(}\frac{1}{p}+(-1)^{p-1}\sum_{k=p}^{2p}\frac{2}{k}\frac{(2k-2)!!\,(2p-1)!!}{(2k-2p)!!\,(k-1)!\,(2p-k)!}e^{-\lambda
k}\bigg{)}.$

Compared to this form, the formula (2.13) obtained from the positive formula
(2.2), which can be viewed as a power series in $1-e^{-\lambda}$, turns out to
be more convenient for further analysis.

We now discuss the spectral density. Wigner’s semi-circle law [3] gives that
the density $\rho_{N}^{\rm(H)}$ in (1.12) satisfies

(2.19) $\frac{1}{\sqrt{N}}\rho_{N}^{\rm(H)}(\sqrt{N}x)\to\rho_{\rm
sc}(x):=\frac{\sqrt{4-x^{2}}}{2\pi}\,\mathbbm{1}_{(-2,2)}(x),\qquad
N\to\infty.$

On the other hand, it follows from (1.38) and (1.50) that the maximal support
of $\widehat{\rho}_{N}^{\rm(dH)}(\sqrt{N}x)$ is

(2.20) $|x|<\frac{1}{\sqrt{N(1-q)}}=\frac{1}{\sqrt{\lambda}}+O(\frac{1}{N}).$

Thus the limiting scaled density is supported on $|x|<1/\sqrt{\lambda}.$
Another preliminary point is that the left-hand side of (2.12) is the moment
of

(2.21)
$\widetilde{\rho}_{N}^{\rm(dH)}(x):=\sqrt{q}\,\rho_{N}^{\rm(dH)}(\sqrt{q}x)=\sqrt{\frac{q}{1-q}}\,\widehat{\rho}_{N}^{\rm(dH)}\Big{(}\sqrt{\frac{q}{1-q}}x\Big{)}.$

######  Theorem 2.6 (Limiting spectral density and its finite-size
correction).

Let

(2.22) $b(\lambda)=2\sqrt{(1-e^{-\lambda})e^{-\lambda}}.$

Then as $N\to\infty$, in the sense of integration against a smooth function,
we have

(2.23)
$\frac{1}{N}\widetilde{\rho}_{N}^{\rm(dH)}(x)=\widetilde{\rho}_{(0)}^{\rm(dH)}(x)+\frac{\widetilde{\rho}_{(1)}^{\rm(dH)}(x)}{N^{2}}+O(\frac{1}{N^{4}}),\qquad|x|<1,$

where

(2.24)
$\displaystyle\begin{split}\widetilde{\rho}_{(0)}^{\rm(dH)}(x)&=\frac{1}{\pi\lambda
x}\arg\bigg{(}\frac{2e^{-\lambda}-x+i\sqrt{4(1-e^{-\lambda})e^{-\lambda}-x^{2}}}{2e^{-\lambda}+x+i\sqrt{4(1-e^{-\lambda})e^{-\lambda}-x^{2}}}\bigg{)}\mathbbm{1}_{(-b(\lambda),b(\lambda))}(x)\\\
&\quad+\begin{cases}0&\textup{if }\lambda\leq\log 2,\vskip 3.0pt plus 1.0pt
minus 1.0pt\\\
\displaystyle\frac{1}{\lambda|x|}\Big{(}\mathbbm{1}_{(-1,1)}(y)-\mathbbm{1}_{(-b(\lambda),b(\lambda))}(x)\Big{)}&\textup{if
}\lambda>\log 2\end{cases}\end{split}$

and

(2.25)
$\displaystyle\begin{split}\widetilde{\rho}_{(1)}^{\rm(dH)}(x)&=-\frac{1}{\pi}\,\lambda\frac{e^{-\lambda}(1-e^{-\lambda})}{2}(x^{2}-2e^{-\lambda})(b(\lambda)^{2}-x^{2})^{-5/2}\mathbbm{1}_{(-b(\lambda),b(\lambda))}(x)\\\
&\quad-\frac{1}{\pi}\frac{\lambda}{3}\int_{\tfrac{1}{2}-x^{*}}^{\min\\{1-e^{-\lambda},\tfrac{1}{2}+x^{*}\\}}t\Big{(}x^{2}+2t(1-t)\Big{)}\Big{(}4t(1-t)-x^{2}\Big{)}^{-5/2}\,dt.\end{split}$

Here, $x^{*}=\frac{1}{2}\sqrt{1-x^{2}}.$

As exhibited in Figure 2, the graph of $\widetilde{\rho}_{(0)}^{\rm(dH)}(x)$
as implied by (2.24) has a distinct shape depending on the parameter $\lambda$
being less or greater than $\log 2$. We mention that such a phase transition
in the spectral density for random matrix type models on a lattice has also
been observed in earlier works; see e.g. [48].

The density $\widetilde{\rho}_{(0)}^{\rm(dH)}(x)$ is characterised as the
probability density supported on $(-1,1)$ with even positive integer moments
$\mathcal{M}_{2p,0}^{\rm(dH)}$. Thus

(2.26)
$\int|x|^{2p}\widetilde{\rho}_{(0)}^{\rm(dH)}(x)\,dx=\mathcal{M}_{2p,0}^{\rm(dH)}.$

With $\mathcal{M}_{2p,0}^{\rm(dH)}$ given by (2.13), both sides are such that
Carlson’s theorem [62, Section 5.81] hold, telling us that they are uniquely
determined by their values at the non-negative integers. This allows the
validity of (2.26) to be extended to Re$(p)>-1/2$.

(a) $\lambda=0.3$

(b) $\lambda=\log 2$

(c) $\lambda=2$

(d) $\lambda=0.3$

(e) $\lambda=\log 2$

(f) $\lambda=2$

Figure 2. (A)–(C) shows the spectral density
$x\mapsto\widetilde{\rho}_{(0)}^{\rm(dH)}(x)$. Here the vertical line
indicates the support $|x|=1$. (D)–(F) shows the graph
$p\mapsto\mathcal{M}_{2p,0}^{\rm(dH)}$ (red full line) and its comparison with
$p\mapsto\int|x|^{2p}\widetilde{\rho}_{(0)}^{\rm(dH)}(x)\,dx$ (blue dots).

######  Remark 2.7.

For $\lambda<\log 2$, one can also write

(2.27) $\widetilde{\rho}^{\rm(dH)}_{(0)}(x)=\frac{1}{\pi\lambda
x}\bigg{[}\arctan\bigg{(}\frac{\sqrt{4e^{\lambda}-4-e^{2\lambda}x^{2}}}{2-e^{\lambda}x}\bigg{)}-\arctan\bigg{(}\frac{\sqrt{4e^{\lambda}-4-e^{2\lambda}x^{2}}}{2+e^{\lambda}x}\bigg{)}\bigg{]}\mathbbm{1}_{(-b(\lambda),b(\lambda))}(x).$

This form reveals a square-root decay at the edge, i.e.

(2.28) $\widetilde{{\rho}}^{\rm(dH)}_{(0)}(x)\asymp\sqrt{x-b(\lambda)},\qquad
x\to b(\lambda).$

Due to non-integrable singularities, the integral of
$\widetilde{{\rho}}_{(1)}^{\rm(dH)}$ should be interpreted as a meromorphic
continuation of a finite integral, referred to as the Hadamard regularised
form; see [68, Section III B] for a discussion in the context of the large $N$
expansion of the smoothed density for the Gaussian $\beta$ ensemble.

######  Remark 2.8 (Continuum limit $q\to 1$).

One can observe that

(2.29) $\sqrt{\lambda}\,\widetilde{\rho}_{(0)}(\sqrt{\lambda}x)\to\rho_{\rm
sc}(x),\qquad\lambda\to 0.$

Therefore $\widetilde{\rho}_{(0)}$ provides a $q$-extension of the Wigner
semi-circle law. Note also that for $|x|<2$,

(2.30)
$\sqrt{\lambda}\,\widetilde{\rho}_{(1)}(\sqrt{\lambda}x)\to\frac{1}{\pi}(4-x^{2})^{-5/2}\mathbbm{1}_{(-2,2)}(x),\qquad\lambda\to
0.$

Indeed, the first line in (2.25) remains and the second line involving
integral vanishes in the limit $\lambda\to 0$ for $|x|<2/\sqrt{\lambda}$. The
resulting $q\to 1$ limit (2.30) coincides with [68, Eq.(3.18)] with $h=0$ and
$x$ restricted to $(-2,2)$.

##  3\. Combinatorial derivation of the spectral moments

In this section, we review the _Flajolet-Viennot theory_ , a crucial
combinatorial framework in this paper. Exploiting this theory, we advance
towards computing spectral moments of the $q$-deformed GUE, establishing
Theorem 2.1.

Let $\\{P_{j}\\}_{j\geq 0}$ be the sequence of monic orthogonal polynomials
satisfying

(3.1) $\mathcal{L}(P_{j}(x)P_{k}(x))=h_{j}\,\delta_{jk},$

where $\mathcal{L}:\mathbb{R}[x]\to\mathbb{R}$ is the linear functional
associated with the weighted Lebesgue measure $\omega(x)\,dx$, i.e.

(3.2) $\mathcal{L}(P(x))=\int_{\mathbb{R}}P(x)\,\omega(x)\,dx.$

Note that $\mathfrak{m}_{p,j}$ in (1.6) can be written as

(3.3)
$\mathfrak{m}_{p,j}=\frac{\mathcal{L}(x^{p}\,P_{j}(x)^{2})}{\mathcal{L}(P_{j}(x)^{2})}.$

In general, the family of orthogonal polynomials satisfies the three-term
recurrence relation

(3.4) $P_{n+1}(x)=(x-b_{n})P_{n}(x)-\lambda_{n}P_{n-1}(x).$

Here, the orthogonal norm $h_{j}$ can be written in terms of the coefficients
$\lambda_{n}$ as

(3.5) $h_{j}=\lambda_{1}\dots\lambda_{j}.$

###  3.1. Flajolet–Viennot theory

Enumerating lattice paths is a classical subject in combinatorics, often
serving as a combinatorial interpretation for various invariants, with moments
of orthogonal polynomials being a prominent example. Flajolet [24] and Viennot
[65] provided a classical combinatorial interpretation of moments using
Motzkin paths. More recently, Corteel, Jonnadula, Keating and Kim [18] offered
another interpretation of moments using lecture hall graphs. In this section,
we review the Flajolet-Viennot theory. For further exploration of this topic,
we refer to [19] and references therein.

A (lattice) path $\omega$ is a sequence $\omega=(s_{0},\dots,s_{n})$, where
$s_{j}=(x_{j},y_{j})\in\mathbb{N}\times\mathbb{N}.$ Each step
$(s_{j},s_{j+1})$ is called _North-East_ if
$(x_{j+1},y_{j+1})=(x_{j}+1,y_{j}+1)$, _East_ if
$(x_{j+1},y_{j+1})=(x_{j}+1,y_{j})$, and _South-East_ if
$(x_{j+1},y_{j+1})=(x_{j}+1,y_{j}-1)$. We say a step $(s_{j},s_{j+1})$ is of
height $k$ if $y_{j}=k$.

######  Definition 3.1 (Motzkin path and weight).

A _Motzkin path_ is a path consisting of North-East, East, and South-East
steps which lie in the first quadrant. Let $\mathrm{Mot}_{n,k,l}$ be the set
of Motzkin paths from $(0,k)$ to $(n,l).$ For given sequences
$b=\\{b_{n}\\}_{n\geq 0}$ and $\lambda=\\{\lambda_{n}\\}_{n\geq 0}$, and a
Motzkin path $\omega,$ we define the associated weight
$\mathrm{wt}_{b,\lambda}(\omega)$ as a product of the weight of every step in
$\omega$, where the weight for each step is given as

  * •

A North-East step has weight $1$;

  * •

An East step of height $k$ has weight $b_{k}$;

  * •

A South-East step of height $k$ has weight $\lambda_{k}.$

See Figure 3 for an example of a Motzkin path of weight
$b_{0}b_{1}b_{3}\lambda_{1}\lambda_{2}\lambda_{3}^{2}$. For a random matrix
model associated with given orthogonal polynomials defined by the sequences
$b_{n}$ and $\lambda_{n}$, let $\mathfrak{m}_{p,j}$ and ${m}_{N,2p}$ be the
spectral moments. Then we have the following combinatorial interpretation of
spectral moments in terms of Motzkin paths.

$\lambda_{3}$$\lambda_{3}$$\lambda_{2}$$\lambda_{1}$$b_{1}$$b_{3}$$b_{0}$
Figure 3. A Mozkin path with weight
$b_{0}b_{1}b_{3}\lambda_{1}\lambda_{2}\lambda_{3}^{2}$

######  Proposition 3.2 (Combinatorial interpretation of the spectral moments,
cf. [65]).

Let $b,\lambda$,${\mathfrak{m}}_{p,j}$, and ${m}_{N,2p}$ as above. Then we
have

(3.6)
$\displaystyle\begin{split}{\mathfrak{m}}_{p,j}&=\sum_{\omega\in{\rm{Mot}}_{p,j,j}}\mathrm{wt}_{b,\lambda}(\omega),\end{split}$
(3.7) $\displaystyle{m}_{N,2p}$
$\displaystyle=\sum_{j=0}^{N-1}\sum_{\omega\in{\rm{Mot}}_{p,j,j}}\mathrm{wt}_{b,\lambda}(\omega).$

###  3.2. Spectral moments of the GUE revisited

Before going into the analysis of the spectral moments of the $q$-deformed
GUE, it is worth reviewing the original GUE scenario (at $q=1$). The idea of
using Flajolet–Viennot theory to derive the Harer-Zagier formula was
originally proposed by Kerov [42], employing concepts such as involutions and
rook placements. While we follow Kerov’s approach, we adopt the concept of
‘matchings’ to characterize the spectral moments.

A _matching_ on $[n]$ is a set partition of $[n]$ such that each part has size
1 or 2. These parts are termed _isolated vertices_ for those containing one
element and _arcs_ for those containing two elements. For an arc $(a,b)$ with
$a<b$, the vertex $a$ is termed an _opener_ , and the vertex $b$ is termed a
_closer_. We denote the set of matchings on $[a]$ with $b$ arcs as
$\operatorname{Mat}_{a,b}$. Typically, we represent a matching on $[n]$ as a
diagram with $n$ vertices arranged horizontally, connected by arcs.

###### Proof of Proposition 1.1.

We begin by revisiting a combinatorial interpretation of moments associated
with Hermite polynomials, in terms of matching. Let
$\operatorname{Mat}_{a,b}^{>c}\subseteq\operatorname{Mat}_{a,b}$ represent the
set of matchings on $[a]$ with $b$ arcs, ensuring that the first $c$ vertices
are either openers or isolated. Then we claim

(3.8) $\mathfrak{m}_{2p,j}^{\text{ (H)
}}=\left|\operatorname{Mat}_{2p+j,p}^{>j}\right|.$

This claim can be seen as follows. By Proposition 3.2 along with the three-
term recurrence (1.11), we obtain

(3.9) $\mathfrak{m}^{\text{
(H)}}_{2p,j}=\sum_{\omega\in\mathrm{Mot}_{2p,j,j}}\mathrm{wt}_{b,\lambda}(\omega),\qquad(b_{n}=0,\quad\lambda_{n}=n).$

Given $b_{n}=0$, our focus narrows down to the Motzkin paths without an East-
step, so-called the _Dyck paths_. We define a _Hermite history_ as a labeled
Dyck path from $(0,j)$ to $(2p,j)$ where the labeling of its South-East step
of height $k$ is labeled with an integer in $[k]$. An illustration of a
Hermite history is presented on the left side of the bottom of Figure 4. It is
evident that for a Dyck path $\omega$, there are
$\mathrm{wt}_{b,\lambda}(\omega)$-many _Hermite histories_ obtained by
labeling South-East steps of $\omega$.

Subsequently, we establish a bijective correspondence between the set of
Hermite histories from $(0,j)$ to $(2p,j)$ and
$\operatorname{Mat}_{2p+j,p}^{>j}$, in the following dynamic manner.

  * •

Start with $j$ vertices (depicted as blue dots in the right side of Figure 4).

  * •

Progressing through each step of the Dyck path from $(0,j)$, a vertex is added
(depicted as red dots in Figure 4). Consequently, a total of $2p$ vertices are
appended.

  * •

When the Hermite history advances with a South-East step labeled by $m$, we
introduce a new arc connecting the newly added vertex and the $m$-th isolated
point, counted from the left. If it advances with a North-East step, a
positive vertex is added without any accompanying arc.

See Figure 4 for a visual representation of the above bijection.

$2$$2$$2$$2$$2$$2$$1$$2$$2$$1$ Figure 4. A Hermite history from $(0,2)$ to
$(6,2)$ its corresponding matching

For the remainder of the proof, we enumerate
$|\operatorname{Mat}_{2p+j,p}^{>j}|$ in two distinct ways, leading to (1.14)
and (1.16), respectively.

We first show (1.14). Let $\operatorname{Mat}_{2p+j,p}^{>j}(i)$ be the set of
matchings in $\operatorname{Mat}_{2p+j,p}^{>j}$ such that there are $i$ arcs
connecting to the first $j$ vertices. Or equivalently, $j-i$ vertices among
the first $j$ vertices are isolated. We claim that

(3.10)
$\left|\operatorname{Mat}_{2p+j,p}^{>j}(i)\right|=(2p-1)!!\binom{j}{i}\binom{p}{i}2^{i},$

thereby completing the proof. To that end, we construct matchings in
$\operatorname{Mat}_{2p+j,p}^{>j}(i)$ step by step as follows.

  * •

Firstly, let us consider a matching in $\operatorname{Mat}_{2p-i,p-i}$. The
count of such matchings is given by

$\binom{2p-i}{i}(2p-2i-1)!!.$

Here, we select $2p-2i$ vertices out of $2p-i$ available vertices and form
perfect matchings with $p-i$ arcs on the chosen vertices.

  * •

For each matching $M\in\operatorname{Mat}_{2p-i,p-i}$, we introduce $2i$ new
vertices, with the first $i$ vertices placed in the leftmost part.
Subsequently, we create $i$ arcs among newly introduced $2i$ vertices,
connecting each of the first $i$ vertices with the remaining $i$ vertices. The
total count of such configurations is given by

$\binom{2p}{i}i!=\dfrac{(2p)!}{(2p-i)!}.$

  * •

Finally, we introduce $j-i$ isolated vertices in between the first $i$
vertices so that in the resulting matching, newly added vertices are in the
first $j$ vertices. The number of ways to achieve this is given by
$\binom{j}{i}$.

Figure 5 provides an illustrative example of constructing a matching in
$\operatorname{Mat}_{9,3}^{>2}(1)$ following the outlined steps. In the
initial figure on the left, a matching in $\operatorname{Mat}_{5,4}$ is
presented. Subsequently, in the middle figure, a vertex is added to the
leftmost part, followed by the addition of another vertex which are then
connected by a new arc, resulting in a matching in
$\operatorname{Mat}_{9,3}^{>1}(1)$. Finally, another vertex is added, yielding
a matching in $\operatorname{Mat}_{9,3}^{>2}(1)$. By counting the number of
ways in each step, we have

$\displaystyle\left|\operatorname{Mat}_{2p+j,p}^{>j}(i)\right|$
$\displaystyle=\binom{j}{i}\frac{(2p)!}{(2p-i)!}\binom{2p-i}{2i}(2p-2i-1)!!=(2p-1)!!\binom{j}{i}\binom{p}{i}2^{i}.$

This proves the claim, thereby concluding (1.14).

Figure 5. Constructing a matching in $\operatorname{Mat}_{9,3}^{>2}(1)$ step
by step, newly added vertices and arcs in each step are colored in red.

Next, we show (1.16). For an integer $0\leq r<p$ and $j$, we enumerate the
number $m_{j+2p,p,r}$ of matchings on $j+2p$ vertices with $p$ arcs, including
additional markings on $r$ arcs connecting negative vertices.

  * •

Firstly, we select $2r$ vertices from the set of $j$ negative vertices and
connect them with $r$ arcs.

  * •

Subsequently, we choose $2p-2r$ vertices from the remaining $2p+j-2r$ vertices
and connect them with $p-r$ arcs.

This procedural gives

(3.11)
$m_{j+2p,p,r}=\binom{j}{2r}(2r-1)!!\binom{j+2p-2r}{2p-2r}(2p-2r-1)!!=(2p-1)!!\binom{j+2p-2r}{2p}\binom{p}{r}.$

Let $m_{j+2p,p,r}^{i}$ denote the number of matchings on $j+2p$ with $p$ arcs,
where exactly $i$ arcs connect negative vertices, and $r$ of them are marked.
We then establish the following relationship:

(3.12)
$\displaystyle\begin{split}\sum_{r=0}^{p}(-1)^{r}m_{j+2p,p,r}&=(2p-1)!!\sum_{r=0}^{p}\binom{j+2p-2r}{2p}\binom{p}{r}\\\
&=\sum_{r=0}^{p}\sum_{i\geq
0}(-1)^{r}m_{j+2p,p,r}^{i}\binom{i}{r}.\end{split}$

The identity $\sum_{r=0}^{p}(-1)^{r}\binom{i}{r}=\delta_{i,0}$ implies (1.16).
∎

######  Remark 3.3.

The idea can be further extended to evaluate the case $j\not=l$ in general.
Then it follows that

(3.13)
$\displaystyle\frac{1}{(j+2m)!}\,\mathcal{L}\Big{(}x^{2p}He_{j}(x)He_{j+2m}(x)\Big{)}$
$\displaystyle=\sum_{l=0}^{p-m}\binom{j}{l}\frac{(2p)!}{(l+2m)!(p-l-m)!2^{p-l-m}},$
(3.14)
$\displaystyle\frac{1}{(j+2m+1)!}\,\mathcal{L}\Big{(}x^{2p+1}He_{j}(x)He_{j+2m+1}(x)\Big{)}$
$\displaystyle=\sum_{l=0}^{p-m}\binom{j}{l}\frac{(2p+1)!}{(l+2m+1)!(p-l-m)!2^{p-l-m}}.$

See [27, Proposition 13] and references therein for an extensive earlier
literature on such evaluations.

###  3.3. Proof of Theorem 2.1

We will prove Theorem 2.1 (i) using the path counting method. We first
describe the combinatorics of moments of $q$-Hermite polynomials. Let us
recall some statistics for matchings. A _crossing_ is either

  1. (1)

a pair of arcs $(a,b)$ and $(c,d)$ with $a<c<b<d$, or

  2. (2)

a pair consisting of an arc $(a,b)$ and an isolated vertex $c$ with $a<c<b$.

For a matching $M$, we denote the total number of crossings by
$\operatorname{cr}(M)$. Pictorially, the first two diagrams in Figure 6
illustrate two possible cases of crossings. For the second case, one might
conceptualize an isolated vertex as connected to an imaginary arc reaching
‘infinity’, thus forming a ‘crossing’. A _nesting_ is either

  1. (1)

a pair of arcs $(a,b)$ and $(c,d)$ with $a<c<d<b$ or $c<a<b<d$, or

  2. (2)

a pair consisting of an arc $(a,b)$ and an isolated vertex $c$ with $c<a<b$.

For a matching $M$, we denote the total number of nestings by
$\operatorname{ne}(M)$. Pictorially, the last two diagrams in Figure 6
illustrate two possible cases of nestings. For the second case, one might
again think of an isolated vertex as connected to an imaginary arc extending
to ‘infinity’, thus forming a ‘nesting’.

Figure 6. Crossings and nestings

As a $q$-analogue of (3.8), the moments of $q$-Hermite polynomials, defined by
the three-term recurrence in (1.46), admit a combinatorial interpretation in
terms of matchings, crossings, and nestings as follows.

######  Proposition 3.4.

Let $p,j$ be nonnegative integers. Then we have

(3.15)
$\widehat{\mathfrak{m}}_{2p,j}^{\rm(H)}=\sum_{\operatorname{Mat}_{2p+j,p}^{>j}}q^{\operatorname{stat}(M)},$

where $\operatorname{stat}(M)=\operatorname{cr}(M)+2\operatorname{ne}(M).$

###### Proof.

To analyze the statistic $\operatorname{stat}(M)$, let us begin with the
following observation: for a matching $M$ and an arc $\alpha$ of $M$, let
$\operatorname{stat}_{\alpha}(M)$ be defined as

$\operatorname{stat}_{\alpha}(M):=C_{1}+C_{2}+2(N_{1}+N_{2}),$

where $C_{1}$ is the number of arcs that form a crossing with $\alpha$ whose
opener is between the opener and closer of $\alpha$, $C_{2}$ is the number of
isolated vertices that form a crossing with $\alpha$ which are less than the
opener of $\alpha$, $N_{1}$ is the number of arcs that form a nesting with
$\alpha$ whose opener is less than the opener of $\alpha$, and $N_{2}$ is the
number of isolated vertices which are less than the opener of $\alpha$. Then,
it is straightforward to observe that

(3.16) $\operatorname{stat}(M)=\sum_{\alpha\in
M}\operatorname{stat}_{\alpha}(M).$

Recalling the proof of Proposition 1.1, we established a dynamic bijection
between Hermite histories and matchings. This bijection starts with $j$
vertices and introduces a new vertex and an arc to the $m$-th isolated vertex
for a South-East step labeled by $m$, or simply introduces a new vertex for a
North-East step. Following this bijection, let us consider the case
encountering a South-East step of height $k$ labeled by $m$. The corresponding
matching can be described as follows: for a given matching with $k$ isolated
vertices, we introduce a new vertex and connect it with the $m$-th isolated
vertex. Let $\alpha$ denote this additional arc. After the final stage of this
bijection to obtain a matching $M$, we have
$\operatorname{stat}_{\alpha}(M)={(k-1)+(m-1)}$. By considering all $k$
possible labels (1 through $m$) to obtain a new arc $\alpha$, we have the
additional statistic $q^{k-1}[k]_{q}$, which corresponds to the weight
$\lambda_{k}$. When we encounter a North-East step, the weight is 1, and on
the matching side, introducing a new vertex contributes 1 to
$q^{\operatorname{stat}(M)}$. For both cases, the weight of the Hermite
history for each step is recorded by $\operatorname{stat}_{\alpha}(M)$. This,
with (3.16) completes the proof of the claim. ∎

######  Remark 3.5.

For the other version of $q$-Hermite polynomial $\widetilde{H}_{n}(x;q)$,
given in (1.48), we have

(3.17)
$\widetilde{\mathfrak{m}}_{2p,j}^{\rm(H)}=\sum_{\operatorname{Mat}_{2p+j,p}^{>j}}q^{\operatorname{cr}(M)}.$

We are now prepared to establish one of the key theorems in this paper,
Theorem 2.1 (i), which provides a positive formula for the spectral moments of
the $q$-deformed GUE.

###### Proof of Theorem 2.1 (i).

According to Proposition 3.4, it is sufficient to establish the following
equality:

(3.18)
$\displaystyle\sum_{M\in\operatorname{Mat}_{2p+j,p}^{>j}}q^{\operatorname{cr}(M)+2\operatorname{ne}(M)}=\sum_{i\geq
0}q^{(j-i)(2p-i)+i(i-1)/2}\dfrac{[2p]_{q}!}{[2p-2i]_{q}!![i]_{q}!}\genfrac{[}{]}{0.0pt}{}{j}{i}_{q}.$

To achieve this, we extend the idea of proof employed in Proposition 1.1, the
$q=1$ case of Theorem 2.1 given in Section 3.2.

We first refine (3.18) as follows:

(3.19)
$\sum_{M\in\operatorname{Mat}_{2p+j,p}^{>j}(i)}q^{\operatorname{cr}(M)+2\operatorname{ne}(M)}=q^{(j-i)(2p-i)+i(i-1)/2}\dfrac{[2p]_{q}![j]_{q}!}{[2p-2i]_{q}!![i]_{q}![j-2i]_{q}!},$

where $\operatorname{Mat}_{2p+j,p}^{>j}(i)$ denotes the set of matchings in
$\operatorname{Mat}_{2p+j,p}^{>j}$ with $i$ arcs connected to the first $j$
vertices, as in the proof of Proposition 1.1. By summing both sides of (3.19)
over $1\leq i\leq j$, we establish (3.18). Therefore, it remains to validate
(3.19).

We start by claiming that

$\sum_{M\in\operatorname{Mat}_{a,b}}q^{\operatorname{cr}(M)+2\operatorname{ne}(M)}=\genfrac{[}{]}{0.0pt}{}{a}{b}_{q}[2b-1]_{q}!!.$

We establish this assertion through induction on $a$ and $b$. Let
$\operatorname{Mat}_{a,b}$ be partitioned as

$\operatorname{Mat}_{a,b}=\operatorname{Mat}_{a,b}^{(1)}\cup\bigcup_{x>1}\operatorname{Mat}_{a,b}^{(x)}$

Here, we let $\operatorname{Mat}_{a,b}^{(x)}$ be the set of matchings in
$\operatorname{Mat}_{a,b}$ where the first vertex $v=1$ forms an arc with the
vertex $x$. Notably, vertex $v$ is an isolated vertex in
$\operatorname{Mat}_{a,b}^{(1)}$.

For a matching $M\in\operatorname{Mat}_{a,b}^{(1)}$, since $v$ is isolated and
there are $b$ arcs after vertex $v$, vertex $v$ forms $b$ nestings and thus
contributes $2b$ to $\operatorname{stat}(M)$. Additionally, deleting the
vertex $v$ gives a bijection between $\operatorname{Mat}_{a,b}^{(1)}$ and
$\operatorname{Mat}_{a-1,b}$. By the induction hypothesis, we obtain

(3.20)
$\sum_{M\in\operatorname{Mat}_{a,b}^{(1)}}q^{\operatorname{stat}(M)}=q^{2b}\genfrac{[}{]}{0.0pt}{}{a-1}{b}[2b-1]_{q}!!.$

On the other hand, for a matching $M\in\operatorname{Mat}_{a,b}^{(x)}$, the
first arc $(1,x)$ contributes $x-1$ to $\operatorname{stat}(M)$. Furthermore,
deleting the arc $(1,x)$ gives a bijection between
$\operatorname{Mat}_{a,b}^{(x)}$ and $\operatorname{Mat}_{a-2,b-1}$. By the
induction hypothesis, we obtain

(3.21)
$\sum_{x=2}^{a}\sum_{M\in\operatorname{Mat}_{a,b}^{(x)}}q^{\operatorname{stat}(M)}=[a-1]_{q}\genfrac{[}{]}{0.0pt}{}{a-2}{2b-2}[2b-3]_{q}!!.$

Combining equations (3.20) and (3.21), we get

$\displaystyle\sum_{M\in\operatorname{Mat}_{a,b}}q^{\operatorname{stat}(M)}$
$\displaystyle=\sum_{x=1}^{a}\sum_{M\in\operatorname{Mat}_{a,b}^{(x)}}q^{\operatorname{stat}(M)}$
$\displaystyle=q^{2b}\genfrac{[}{]}{0.0pt}{}{a-1}{b}[2b-1]_{q}!!+[a-1]_{q}\genfrac{[}{]}{0.0pt}{}{a-2}{2b-2}[2b-3]_{q}!!$
$\displaystyle=q^{2b}\genfrac{[}{]}{0.0pt}{}{a-1}{b}[2b-1]_{q}!!+\genfrac{[}{]}{0.0pt}{}{a-1}{2b-1}[2b-1]_{q}!!=\genfrac{[}{]}{0.0pt}{}{a}{b}[2b-1]_{q}!!.$

In the third equation, we use the identity
$[n]_{q}\genfrac{[}{]}{0.0pt}{}{n-1}{k}_{q}=[n-k]_{q}\genfrac{[}{]}{0.0pt}{}{n}{k}_{q}$,
and in the last equation, we use the $q$-Pascal identity. In particular, by
taking $a=2p-i$ and $b=p-i$, we have

(3.22)
$\sum_{M\in\operatorname{Mat}_{2p-i,p-i}}q^{\operatorname{stat}(M)}=\genfrac{[}{]}{0.0pt}{}{2p-i}{i}[2p-2i-1]_{q}!!.$

We proceed to prove a special case of (3.19) where $j=i$,

(3.23)
$\sum_{M\in\operatorname{Mat}_{2p+i,p,i}^{<i}(i)}q^{\operatorname{stat}(M)}=q^{i(i-1)/2}\dfrac{[2p]_{q}!}{[2p-i]_{q}!}\genfrac{[}{]}{0.0pt}{}{2p-i}{2p-2i}[2p-2i-1]_{q}!!.$

To prove this, fix a matching $M\in\operatorname{Mat}_{2p-i,p-i}$. We
introduce $i$ new vertices to the left of $M$, placing the resulting matching
in $\operatorname{Mat}_{2p,p-i}^{<i}(0)$. Then, we introduce new arcs and
their closers larger than $i$ one by one by connecting vertex $1$ first,
followed by vertex $2$, and so forth. When introducing the first arc, subject
to the conditions mentioned above, there are $(2p-i+1)$ choices. These choices
contribute a factor $q^{i-1}[2p-i+1]_{q}$. Similarly, when introducing the
second arc, there are $(2p-i+2)$ choices, resulting in another factor
$q^{i-2}[2p-i+2]_{q}$. By inductively applying this argument, we establish the
claim.

Finally, consider a projection map

$\mathcal{P}:\operatorname{Mat}_{2p+j,p}^{<j}(i)\rightarrow\operatorname{Mat}_{2p+i,p}^{<i}(i)$

which is defined by deleting $j-i$ isolated vertices less than $j+1$. Suppose
we fix a matching $M\in\operatorname{Mat}_{2p+i,p}^{<i}(i)$. Then, there are
$\binom{j}{i}$ choices for introducing $j-i$ new isolated vertices to $M$ to
obtain a matching in $\operatorname{Mat}_{2p+j,p}^{<j}(i)$. For each such
choice, we observe that the additional contribution of an isolated vertex $v$
to $q^{\operatorname{stat}(M)}$ is $q^{(2p-i)+N}$, where $N$ is the count of
chosen vertices less than $j+1$ larger than $v$. Utilizing the standard
combinatorial interpretation of $q$-binomial coefficients, we conclude that

$\sum_{\tilde{M}\in\mathcal{P}^{-1}(M)}q^{\operatorname{stat}(\tilde{M})}=q^{\operatorname{stat}(M)+(j-1)(2p-i)}\genfrac{[}{]}{0.0pt}{}{j}{i}.$

Combining this with (3.23), we have

$\displaystyle\sum_{M\in\operatorname{Mat}_{2p+j,p}^{>j}(i)}q^{\operatorname{stat}(M)}$
$\displaystyle=q^{(j-1)(2p-i)+i(i-1)/2}\dfrac{[2p]_{q}!}{[2p-i]_{q}!}\genfrac{[}{]}{0.0pt}{}{2p-i}{2p-2i}[2p-2i-1]_{q}!!\genfrac{[}{]}{0.0pt}{}{j}{i}$
$\displaystyle=q^{(j-i)(2p-i)+i(i-1)/2}\dfrac{[2p]_{q}![j]_{q}!}{[2p-2i]_{q}!![i]_{q}![j-2i]_{q}!}.$

∎

######  Remark 3.6.

We expect that the alternating series formula (2.3) also has a combinatorial
proof. Let us rewrite (2.3) as

(3.24)
$\sum_{M\in\operatorname{Mat}_{2p+j,p}^{>j}}q^{\operatorname{stat}(M)}=\sum_{r=0}^{p}(-1)^{r}q^{r(r-1)}\genfrac{[}{]}{0.0pt}{}{j}{2r}_{q}[2r-1]_{q}!!\genfrac{[}{]}{0.0pt}{}{j+2p-2r}{2p-2r}_{q}[2p-2r-1]_{q}!!.$

We introduce further notation to refine the above identity. For a matching
$m\in\operatorname{Mat}_{2p+j,p}$, let $\operatorname{cl}^{\leq j}(m)$ denote
the count of closers less than or equal to $j$. Then, for $p,j\geq 0$ and
$r\leq p,\frac{j}{2}$, we expect that

(3.25)
$\sum_{M\in\operatorname{Mat}_{2p+j,p}}q^{\operatorname{stat}(M)}\genfrac{[}{]}{0.0pt}{}{\operatorname{cl}^{\leq
j}(m)}{r}_{q^{2}}=\genfrac{[}{]}{0.0pt}{}{j}{2r}_{q}[2r-1]_{q}!!\genfrac{[}{]}{0.0pt}{}{j+2p-2r}{2p-2r}_{q}[2p-2r-1]_{q}!!.$

This identity can be verified at small order using computer algebra. It is
worth noting that (3.24) follows from (3.25). For this purpose, let
$A_{p,j,k}$ be the generating function of
$q^{\operatorname{cr}(M)+2\operatorname{ne}(M)}$ over all matchings
$m\in\operatorname{Mat}_{2p+j,p}$ with exactly $k$ closers are less or equal
to $j$. Note that $A_{p,j,0}$ is equal to the left-hand side of (3.24). In
addition, the left-hand side of (3.25), denoted by $a_{p,j,r}$, is given by

$a_{p,j,r}=\sum_{k\geq r}A_{p,j,k}\genfrac{[}{]}{0.0pt}{}{k}{r}_{q^{2}}.$

Using the identity
$\sum_{r=0}^{k}(-1)^{r}q^{r(r-1)}\genfrac{[}{]}{0.0pt}{}{k}{r}_{q^{2}}=0$ for
$r\geq 1$, we have

$A_{p,j,0}=\sum_{r\geq 0}(-1)^{r}q^{r(r-1)}a_{p,j,r},$

which leads to (3.24).

##  4\. Scaled moments and density

###  4.1. Proof of Theorem 2.3

In this subsection, we prove Theorem 2.3.

Let us write

(4.1)
$F_{N}(l):=q^{p-l(2p-l)+l(l-1)/2}(1-q)^{p}\dfrac{[2p]_{q}!}{[2p-2l]_{q}!!\,[l]_{q}!}\sum_{j=0}^{N-1}q^{j(2p-l)}\genfrac{[}{]}{0.0pt}{}{j}{l}_{q}.$

Then we have

(4.2)
$\Big{(}q(1-q)\Big{)}^{p}\widehat{m}_{N,2p}^{\rm(dH)}=\sum_{l=0}^{p}F_{N}(l).$

We shall verify that for each $l$, $F_{N}(l)$ has a $1/N$-expansion with a
leading order of $O(N^{l+1-p})$. Hence, to prove the theorem, the specific
large-$N$ asymptotic behaviours of the three terms $F_{N}(p)$, $F_{N}(p-1)$,
and $F_{N}(p-2)$ are required. In particular, one needs to precisely analyze
the summation in (4.1), for which we will implement some ideas from recent
works [15, 2].

Recall that $s(n,k)$ denotes the Stirling numbers (1.24). We first provide the
asymptotic behaviours of the pre-factors in (4.1).

######  Lemma 4.1.

Let $q$ be scaled as in (2.10). Then as $N\to\infty,$ we have

$\displaystyle\quad(1-q)^{p}\dfrac{[2p]_{q}!}{[2p-2l]_{q}!!\,[l]_{q}!}=\Big{(}\frac{\lambda}{N}\Big{)}^{p}\frac{(2p)!}{(2p-2l)!!\,l!}\bigg{(}1+\frac{3l^{2}-4lp-2p^{2}-l}{4}\frac{\lambda}{N}+\mathfrak{a}(p,l)\frac{\lambda^{2}}{N^{2}}+O(\frac{1}{N^{3}})\bigg{)},$

where

$\displaystyle\mathfrak{a}(p,l)$
$\displaystyle=\frac{s(2p+1,2p-1)}{4}-s(p-l+1,p-l-1)-\frac{s(l+1,l-1)}{4}$
$\displaystyle\quad+\frac{1}{144}(12l-33l^{2}-30l^{3}+63l^{4}-8p+60lp+48l^{2}p-180l^{3}p+6lp^{2}+126l^{2}p^{2}-4p^{3}-36p^{4}).$

###### Proof.

Note that by definition,

$\displaystyle(1-q)^{p}\dfrac{[2p]_{q}!}{[2p-2l]_{q}!!\,[l]_{q}!}$
$\displaystyle=\frac{\prod_{j=1}^{2p}(1-q^{j})}{\prod_{j=1}^{p-l}(1-q^{2j})\prod_{j=1}^{l}(1-q^{j})}.$

Since

$1-q^{m}=1-e^{-m\lambda/N}=\frac{m\lambda}{N}\sum_{r=0}^{\infty}\frac{(-1)^{r}}{(r+1)!}\Big{(}\frac{m\lambda}{N}\Big{)}^{r},$

we have

$\displaystyle\quad\prod_{j=1}^{2p}(1-q^{j})=\Big{(}\frac{\lambda}{N}\Big{)}^{2p}(2p)!\prod_{m=1}^{2p}\bigg{(}1-\frac{m\lambda}{2N}+\frac{(m\lambda)^{2}}{6N^{2}}+O(\frac{1}{N^{3}})\bigg{)}$
$\displaystyle=\Big{(}\frac{\lambda}{N}\Big{)}^{2p}(2p)!\bigg{(}1-\frac{p(2p+1)}{2}\frac{\lambda}{N}+\Big{(}\frac{s(2p+1,2p-1)}{4}+\frac{p(2p+1)(4p+1)}{18}\Big{)}\frac{\lambda^{2}}{N^{2}}+O(\frac{1}{N^{3}})\bigg{)},$
$\displaystyle\prod_{j=1}^{p-l}(1-q^{2j})$
$\displaystyle=\Big{(}\frac{\lambda}{N}\Big{)}^{p-l}(2p-2l)!!\bigg{(}1-\frac{(p-l)(p-l+1)}{2}\frac{\lambda}{N}$
$\displaystyle\quad+\Big{(}s(p-l+1,p-l-1)+\frac{(p-l)(p-l+1)(2p-2l+1)}{9}\Big{)}\frac{\lambda^{2}}{N^{2}}+O(\frac{1}{N^{3}})\bigg{)},$

and

$\displaystyle\prod_{j=1}^{l}(1-q^{j})=\Big{(}\frac{\lambda}{N}\Big{)}^{l}l!\bigg{(}1-\frac{l(l+1)}{4}\frac{\lambda}{N}+\Big{(}\frac{s(l+1,l-1)}{4}+\frac{l(l+1)(2l+1)}{36}\Big{)}\frac{\lambda^{2}}{N^{2}}+O(\frac{1}{N^{3}})\bigg{)}.$

This completes the proof. ∎

We now begin with analyzing the summation in (4.1). For this, let

(4.3) $f_{l}(j)=q^{j(2p-l)}\genfrac{[}{]}{0.0pt}{}{j}{l}_{q}.$

######  Lemma 4.2.

Let $q$ be scaled as in (2.10). Then as $N\to\infty,$ we have

(4.4) $\displaystyle f_{l}(Nt)=\frac{N^{l}}{l!\,\lambda^{l}}e^{-(2p-l)\lambda
t}(1-e^{-\lambda
t})^{l}\Bigg{(}1+\mathcal{B}_{l,1}(t)\frac{\lambda}{N}+\mathcal{B}_{l,2}(t)\frac{\lambda^{2}}{N^{2}}+O(N^{-3})\bigg{)},$

where

(4.5) $\displaystyle\mathcal{B}_{l,1}(t)$
$\displaystyle=\frac{l(l+1)}{4}-\frac{e^{-\lambda t}}{1-e^{-\lambda
t}}\frac{l(l-1)}{2},$ (4.6)
$\displaystyle\begin{split}\mathcal{B}_{l,2}(t)&=\Big{(}\frac{e^{-\lambda
t}}{1-e^{-\lambda t}}\Big{)}^{2}s(l,l-2)-\frac{e^{-\lambda t}}{1-e^{-\lambda
t}}\frac{(-1+l)l(-2+7l+3l^{2})}{24}\\\
&\quad-\frac{s(l+1,l-1)}{4}+\frac{l(1+l)(-4+l+9l^{2})}{144}.\end{split}$

Furthermore, we have

(4.7) $\displaystyle
N\int_{0}^{1}f_{l}(Nt)\,dt=\Big{(}\frac{N}{\lambda}\Big{)}^{l+1}\bigg{(}\mathcal{A}_{l,0}+\mathcal{A}_{l,1}\Big{(}\frac{N}{\lambda}\Big{)}^{-1}+\mathcal{A}_{l,2}\Big{(}\frac{N}{\lambda}\Big{)}^{-2}+O(N^{-3})\bigg{)},$

where

(4.8) $\mathcal{A}_{l,0}=\frac{(2p-l-1)!}{(2p)!}I_{1-e^{-\lambda}}(l+1,2p-l),$
(4.9)
$\mathcal{A}_{l,1}=\frac{l(l+1)}{4}\frac{(2p-l-1)!}{(2p)!}I_{1-e^{-\lambda}}(l+1,2p-l)-\frac{(l-1)}{2}\frac{(2p-l)!}{(2p)!}I_{1-e^{-\lambda}}(l,2p-l+1),$

and

(4.10)
$\displaystyle\begin{split}\mathcal{A}_{l,2}&=\Big{(}-\frac{s(l+1,l-1)}{4}+\frac{l(1+l)(-4+l+9l^{2})}{144}\Big{)}\frac{(2p-l-1)!}{(2p)!}I_{1-e^{-\lambda}}(l+1,2p-l)\\\
&\quad-\frac{(l-1)(-2+7l+3l^{2})}{24}\frac{(2p-l)!}{(2p)!}I_{1-e^{-\lambda}}(l,2p-l+1)\\\
&\quad+s(l,l-2)\frac{1}{l(l-1)}\frac{(2p-l+1)!}{(2p)!}I_{1-e^{-\lambda}}(l-1,2p-l+2).\end{split}$

###### Proof.

With the scaling (2.10) and with $j=tN$,

$\displaystyle\genfrac{[}{]}{0.0pt}{}{j}{l}_{q}$
$\displaystyle=\frac{[j]_{q}\,!}{[l]_{q}!\,[j-l]_{q}!}=\frac{(1-q^{j})\dots(1-q^{j-l+1})}{(1-q^{l})\dots(1-q)}=\frac{N^{l}}{l!\,\lambda^{l}}(1-e^{-\lambda
t})^{l}\Big{(}1+O(N^{-1})\Big{)},$

as $N\to\infty.$ Here,

$\displaystyle\quad(1-q^{j})\dots(1-q^{j-l+1})=(1-e^{-\lambda
t})^{l}\prod_{m=0}^{l-1}\Big{(}1-\frac{e^{-\lambda t}}{1-e^{-\lambda
t}}\frac{m\lambda}{N}-\frac{e^{-\lambda t}}{1-e^{-\lambda
t}}\frac{m^{2}\lambda^{2}}{2N^{2}}+O(N^{-3})\Big{)}$
$\displaystyle=(1-e^{-\lambda t})^{l}\bigg{(}1-\frac{e^{-\lambda
t}}{1-e^{-\lambda t}}\frac{l(l-1)}{2}\frac{\lambda}{N}+\frac{e^{-\lambda
t}}{1-e^{-\lambda t}}\Big{(}\frac{e^{-\lambda t}}{1-e^{-\lambda
t}}s(l,l-2)-\frac{(l-1)l(2l-1)}{12}\Big{)}\frac{\lambda^{2}}{N^{2}}+O(N^{-3})\bigg{)}$

Since

$\displaystyle\quad\Big{(}(1-q^{l})(1-q^{l-1})\dots(1-q)\Big{)}^{-1}$
$\displaystyle=\frac{N^{l}}{l!\,\lambda^{l}}\bigg{(}1+\frac{l(l+1)}{4}\frac{\lambda}{N}+\Big{(}-\frac{s(l+1,l-1)}{4}+\frac{l^{2}(l+1)^{2}}{16}-\frac{l(l+1)(2l+1)}{36}\Big{)}\frac{\lambda^{2}}{N^{2}}+O(N^{-3})\bigg{)}$

we have

$\displaystyle\quad\genfrac{[}{]}{0.0pt}{}{j}{l}_{q}\Big{(}\frac{N^{l}}{l!\,\lambda^{l}}(1-e^{-\lambda
t})^{l}\Big{)}^{-1}$
$\displaystyle=1+\bigg{(}\frac{l(l+1)}{4}-\frac{e^{-\lambda t}}{1-e^{-\lambda
t}}\frac{l(l-1)}{2}\bigg{)}\frac{\lambda}{N}+\bigg{(}\Big{(}\frac{e^{-\lambda
t}}{1-e^{-\lambda t}}\Big{)}^{2}s(l,l-2)$
$\displaystyle\qquad-\frac{e^{-\lambda t}}{1-e^{-\lambda
t}}\frac{(-1+l)l(-2+7l+3l^{2})}{24}-\frac{s(l+1,l-1)}{4}+\frac{l(1+l)(-4+l+9l^{2})}{144}\bigg{)}\frac{\lambda^{2}}{N^{2}}+O(N^{-3}),$

which completes the proof of the first assertion. The second assertion follows
from straightforward computations using

$\displaystyle\quad\int_{0}^{1}e^{-\lambda t\,a}(1-e^{-\lambda
t})^{b}\,dt=\frac{1}{\lambda}\int_{0}^{1-e^{-\lambda}}s^{b}(1-s)^{a-1}\,ds$
$\displaystyle=\frac{1}{\lambda}B_{1-e^{-\lambda}}(b+1,a)=\frac{1}{\lambda}\frac{b!\,(a-1)!}{(a+b)!}I_{1-e^{-\lambda}}(b+1,a).$

∎

Next, we derive the asymptotic behaviour of $\sum_{j=0}^{N-1}f(j)$.

######  Lemma 4.3.

Let $q$ be scaled as in (2.10). Then as $N\to\infty,$ we have

(4.11) $\displaystyle\sum_{j=0}^{N-1}f_{l}(j)$
$\displaystyle=\mathcal{C}_{l,0}\Big{(}\frac{N}{\lambda}\Big{)}^{l+1}+\mathcal{C}_{l,1}\Big{(}\frac{N}{\lambda}\Big{)}^{l}+\mathcal{C}_{l,2}\Big{(}\frac{N}{\lambda}\Big{)}^{l-1}+O\Big{(}\frac{N}{\lambda}\Big{)}^{l-2},$

where

(4.12) $\displaystyle\mathcal{C}_{l,0}$ $\displaystyle=\mathcal{A}_{l,0},$
(4.13) $\displaystyle\mathcal{C}_{l,1}$
$\displaystyle=\mathcal{A}_{l,1}-\frac{1}{2}\frac{1}{l!}e^{-(2p-l)\lambda}(1-e^{-\lambda})^{l},$
(4.14) $\displaystyle\mathcal{C}_{l,2}$
$\displaystyle=\mathcal{A}_{l,2}+\frac{1}{12}\frac{1}{l!}e^{-(2p-l)\lambda}(1-e^{-\lambda})^{l}\,\Big{(}-\frac{3l^{2}+l+4p}{2}+(3l^{2}-2l)\frac{e^{-\lambda}}{1-e^{-\lambda}}\Big{)}.$

Here $\mathcal{A}_{l,0}$, $\mathcal{A}_{l,1}$ and $\mathcal{A}_{l,0}$ are
given by (4.8), (4.9) and (4.10), respectively.

###### Proof.

We apply the Euler–Maclaurin formula [58, Section 2.10]

(4.15)
$\begin{split}\sum_{j=m}^{n}f(j)=\int_{m}^{n}f(x)\,dx+\frac{f(m)+f(n)}{2}+\sum_{k=1}^{l-1}\frac{B_{2k}}{(2k)!}\Big{(}f^{(2k-1)}(n)-f^{(2k-1)}(m)\Big{)}+R_{l},\end{split}$

where $B_{n}$ is the Bernoulli number defined by the generating function

(4.16) $\frac{t}{e^{t}-1}=\sum_{n=0}^{\infty}B_{n}\frac{t^{n}}{n!}.$

Here the error term $R_{l}$ satisfies the estimate

$|R_{l}|\leq\frac{4\,\zeta(2l)}{(2\pi)^{2l}}\int_{m}^{n}|f^{(2l)}(x)|\,dx,$

where $\zeta$ is the Riemann zeta function. By using (4.15), we have

$\displaystyle\sum_{j=0}^{N-1}f_{l}(j)$
$\displaystyle=-f_{l}(N)+\sum_{j=0}^{N}f_{l}(j)$
$\displaystyle=\int_{0}^{N}f_{l}(x)\,dx+\frac{f_{l}(0)-f_{l}(N)}{2}+\sum_{k=1}^{m-1}\frac{B_{2k}}{(2k)!}\Big{(}f_{l}^{(2k-1)}(N)-f_{l}^{(2k-1)}(0)\Big{)}+R_{m}.$

Note that

$\displaystyle\int_{0}^{N}f(x)\,dx=N\int_{0}^{1}q^{2Nt(2p-l)}\genfrac{[}{]}{0.0pt}{}{Nt}{l}_{q}\,dt.$

Furthermore, as a consequence of Lemma 4.2, we have

$\displaystyle
f_{l}(N)=\frac{N^{l}}{l!\,\lambda^{l}}e^{-(2p-l)\lambda}(1-e^{-\lambda})^{l}\Bigg{(}1+\Big{(}\frac{l(l+1)}{4}-\frac{e^{-\lambda}}{1-e^{-\lambda}}\frac{l(l-1)}{2}\Big{)}\frac{\lambda}{N}++O(N^{-2})\bigg{)}$

and

$\displaystyle
f_{l}^{\prime}(N)=\frac{N^{l}}{l!\,\lambda^{l}}e^{-(2p-l)\lambda}(1-e^{-\lambda})^{l}\,\frac{\lambda}{N}\Big{(}-(2p-l)+l\frac{e^{-\lambda}}{1-e^{-\lambda}}\Big{)}.$

Combining all of the above with $B_{2}=1/6$, the lemma follows. ∎

We are now ready to prove Theorem 2.3.

###### Proof of Theorem 2.3.

By Lemma 4.1, we have

$\displaystyle\quad
q^{p-l(2p-l)+l(l-1)/2}(1-q)^{p}\dfrac{[2p]_{q}!}{[2p-2l]_{q}!!\,[l]_{q}!}\bigg{(}\Big{(}\frac{\lambda}{N}\Big{)}^{p}\frac{(2p)!}{(2p-2l)!!\,l!}\bigg{)}^{-1}\bigg{|}_{l=p}$
$\displaystyle=1-\frac{p(p+3)}{4}\frac{\lambda}{N}+\bigg{(}\frac{p(4+63p+20p^{2}-63p^{3})}{144}+\frac{s(2p+1,2p-1)}{4}-\frac{s(p+1,p-1)}{4}\bigg{)}\frac{\lambda^{2}}{N^{2}}+O(N^{-3}),$

and

$\displaystyle
q^{p-l(2p-l)+l(l-1)/2}(1-q)^{p}\dfrac{[2p]_{q}!}{[2p-2l]_{q}!!\,[l]_{q}!}\bigg{(}\Big{(}\frac{\lambda}{N}\Big{)}^{p}\frac{(2p)!}{(2p-2l)!!\,l!}\bigg{)}^{-1}\bigg{|}_{l=p-1}=1-\frac{p^{2}+p+4}{4}\frac{\lambda}{N}+O(N^{-2}).$

Therefore by (4.1), we have

$\displaystyle F_{N}(p)$
$\displaystyle=\frac{N}{\lambda}\frac{(2p)!}{p!}\mathcal{C}_{p,0}+\frac{(2p)!}{p!}\bigg{[}\mathcal{C}_{p,1}-\frac{p(p+3)}{4}\mathcal{C}_{p,0}\bigg{]}$
$\displaystyle\quad+\Big{(}\frac{N}{\lambda}\Big{)}^{-1}\frac{(2p)!}{p!}\bigg{[}\mathcal{C}_{p,2}-\frac{p(p+3)}{4}\mathcal{C}_{p,1}$
$\displaystyle\qquad+\bigg{(}\frac{p(4+63p+20p^{2}-63p^{3})}{144}+\frac{s(2p+1,2p-1)}{4}-\frac{s(p+1,p-1)}{4}\bigg{)}\mathcal{C}_{p,0}\bigg{]}+O(N^{-2}).$

Also we have

$\displaystyle F_{N}(p-1)$
$\displaystyle=\frac{1}{2}\frac{(2p)!}{(p-1)!}\mathcal{C}_{p-1,0}+\Big{(}\frac{N}{\lambda}\Big{)}^{-1}\frac{1}{2}\frac{(2p)!}{(p-1)!}\bigg{[}\mathcal{C}_{p-1,1}-\frac{p^{2}+p+4}{4}\mathcal{C}_{p-1,0}\bigg{]}+O(N^{-2}),$
$\displaystyle F_{N}(p-2)$
$\displaystyle=\Big{(}\frac{N}{\lambda}\Big{)}^{-1}\frac{1}{8}\frac{(2p)!}{(p-2)!}\mathcal{C}_{p-2,0}+O(N^{-2}).$

Combining these asymptotic behaviours, it follows that

(4.17) $\displaystyle
F_{N}(p)+F_{N}(p-1)+F_{N}(p-2)=\mathcal{D}_{0}\,\frac{N}{\lambda}+\mathcal{D}_{1}+\mathcal{D}_{2}\frac{\lambda}{N}+O(N^{-2}),$

where

$\displaystyle\mathcal{D}_{0}=\frac{(2p)!}{p!}\mathcal{C}_{p,0},\qquad\mathcal{D}_{1}=\frac{(2p)!}{p!}\bigg{[}\mathcal{C}_{p,1}-\frac{p(p+3)}{4}\mathcal{C}_{p,0}\bigg{]}+\frac{1}{2}\frac{(2p)!}{(p-1)!}\mathcal{C}_{p-1,0}$

and

$\displaystyle\begin{split}\mathcal{D}_{2}&=\frac{(2p)!}{p!}\bigg{[}\mathcal{C}_{p,2}-\frac{p(p+3)}{4}\mathcal{C}_{p,1}+\bigg{(}\frac{p(4+63p+20p^{2}-63p^{3})}{144}+\frac{s(2p+1,2p-1)}{4}-\frac{s(p+1,p-1)}{4}\bigg{)}\mathcal{C}_{p,0}\bigg{]}\\\
&\quad+\frac{1}{2}\frac{(2p)!}{(p-1)!}\bigg{[}\mathcal{C}_{p-1,1}-\frac{p^{2}+p+4}{4}\mathcal{C}_{p-1,0}\bigg{]}+\frac{1}{8}\frac{(2p)!}{(p-2)!}\mathcal{C}_{p-2,0}.\end{split}$

Write $e^{-\lambda}=x.$ By [58, Eq.(8.17.18)], we have

$\displaystyle I_{1-x}(p,p+1)$
$\displaystyle=I_{1-x}(p+1,p)+\frac{(2p)!}{(p!)^{2}}x^{p}(1-x)^{p},$
$\displaystyle I_{1-x}(p-1,p+2)$
$\displaystyle=I_{1-x}(p+1,p)+\frac{(2p)!}{(p-1)!(p+1)!}x^{p+1}(1-x)^{p-1}+\frac{(2p)!}{(p!)^{2}}x^{p}(1-x)^{p}.$

Using these, we have

$\displaystyle\mathcal{C}_{p,0}$
$\displaystyle=\frac{(p-1)!}{(2p)!}I_{1-x}(p+1,p),$
$\displaystyle\mathcal{C}_{p-1,0}$
$\displaystyle=\frac{p!}{(2p)!}I_{1-x}(p+1,p)+\frac{x^{p}(1-x)^{p}}{p!},$
$\displaystyle\mathcal{C}_{p-2,0}$
$\displaystyle=\frac{(p+1)!}{(2p)!}I_{1-x}(p+1,p)+\frac{x^{p+1}(1-x)^{p-1}}{(p-1)!}+\frac{p+1}{p!}x^{p}(1-x)^{p}$

and

$\displaystyle\mathcal{C}_{p,1}$
$\displaystyle=\frac{-p+3}{4}\frac{p!}{(2p)!}I_{1-x}(p+1,p)-\frac{1}{2}\frac{1}{(p-1)!}x^{p}(1-x)^{p},$
$\displaystyle\mathcal{C}_{p-1,1}$
$\displaystyle=\frac{-p^{2}+p+4}{4}\frac{p!}{(2p)!}I_{1-x}(p+1,p)+\frac{-p^{2}+p+4}{4}\frac{x^{p}(1-x)^{p}}{p!}-\frac{1}{2(p-2)!}x^{p+1}(1-x)^{p-1}.$

Also we have

$\displaystyle\mathcal{C}_{p,2}$
$\displaystyle=\Big{(}-\frac{s(p+1,p-1)}{4}+\frac{p+1}{p-1}s(p,p-2)-\frac{p(16-51p+14p^{2}+9p^{3})}{144}\Big{)}\frac{(p-1)!}{(2p)!}I_{1-x}(p+1,p)$
$\displaystyle\quad+s(p,p-2)\Big{(}\frac{1}{(p-1)\,p!}x^{p+1}(1-x)^{p-1}+\frac{(p+1)}{p(p-1)}\frac{1}{p!}x^{p}(1-x)^{p}\Big{)}$
$\displaystyle\quad+\frac{-2+4p-7p^{2}-3p^{3}}{24}\frac{1}{p!}x^{p}(1-x)^{p}+\frac{1}{12}\frac{1}{p!}(3p^{2}-2p)x^{p+1}(1-x)^{p-1}.$

Using these, after some computations, we obtain

(4.18)
$\mathcal{D}_{1}=0,\qquad\mathcal{D}_{2}=-\frac{p}{6}I_{1-x}(p+1,p)-\frac{(2p-1)!}{6(p-1)!^{2}}x^{p}(1-x)^{p-1}\Big{(}2+p-(2p+1)x\Big{)}.$

Here, we have used

(4.19) $s(n,n-2)=\frac{(n-2)(n-1)n(3n-1)}{24}.$

With the same strategy, one can also demonstrate that the next
$O(N^{-2})$-term vanishes, and we omit the details. This shows the theorem. ∎

###  4.2. Proof of Theorem 2.6

We first show the leading order of Theorem 2.6.

######  Proposition 4.4 (Limiting spectral density for the $q$-deformed GUE).

Let

(4.20)
$a(\lambda)=\frac{2\sqrt{(1-e^{-\lambda})e^{-\lambda}}}{\sqrt{\lambda}}$

Then for $|x|<1/\sqrt{\lambda}$, we have

(4.21)
$\frac{1}{\sqrt{N}}\widehat{\rho}_{N}^{\rm(dH)}(\sqrt{N}x)\to\widehat{\rho}^{\rm(dH)}(x),$

where $\widehat{\rho}^{\rm(dH)}$ is given as follows.

  * •

For $\lambda\leq\log 2$,

(4.22) $\widehat{\rho}^{\rm(dH)}(x)=\frac{1}{\pi\lambda
x}\arg\bigg{(}\frac{2e^{-\lambda}-\sqrt{\lambda}x+i\sqrt{4(1-e^{-\lambda})e^{-\lambda}-\lambda
x^{2}}}{2e^{-\lambda}+\sqrt{\lambda}x+i\sqrt{4(1-e^{-\lambda})e^{-\lambda}-\lambda
x^{2}}}\bigg{)}\mathbbm{1}_{(-a(\lambda),a(\lambda))}(x)$

  * •

For $\lambda\geq\log 2$,

(4.23)
$\displaystyle\begin{split}\widehat{\rho}^{\rm(dH)}(x)&=\frac{1}{\pi\lambda
x}\arg\bigg{(}\frac{2e^{-\lambda}-\sqrt{\lambda}x+i\sqrt{4(1-e^{-\lambda})e^{-\lambda}-\lambda
x^{2}}}{2e^{-\lambda}+\sqrt{\lambda}x+i\sqrt{4(1-e^{-\lambda})e^{-\lambda}-\lambda
x^{2}}}\bigg{)}\mathbbm{1}_{(-a(\lambda),a(\lambda))}(x)\\\
&\quad+\frac{1}{\lambda|x|}\Big{(}\mathbbm{1}_{(-1/\sqrt{\lambda},1/\sqrt{\lambda}}(y)-\mathbbm{1}_{(-a(\lambda),a(\lambda))}(y)\Big{)}.\end{split}$

###### Proof.

Let

(4.24)
$G(y):=\int_{-1/\sqrt{\lambda}}^{1/\sqrt{\lambda}}\frac{\widehat{\rho}^{\rm(dH)}(x)}{y-x}\,dx,\qquad
y\in\mathbb{C}\setminus\mathbb{R}$

be the resolvent. Note that it satisfies

(4.25)
$\displaystyle\begin{split}G(y)&=\sum_{p=0}^{\infty}\frac{1}{y^{2p+1}}\int_{-1/\sqrt{\lambda}}^{1/\sqrt{\lambda}}x^{2p}\,\widehat{\rho}^{\rm(dH)}(x)\,dx.\end{split}$

Then it follows that

$\displaystyle G(y)$
$\displaystyle=\frac{1}{y}+\frac{1}{y}\sum_{p=1}^{\infty}\frac{1}{y^{2p}}\frac{1}{\lambda^{p+1}}\binom{2p}{p}\,\int_{0}^{1-e^{-\lambda}}t^{p}(1-t)^{p-1}\,dt$
$\displaystyle=\frac{1}{y}+\frac{1}{\lambda
y}\int_{0}^{1-e^{-\lambda}}\frac{1}{1-t}\bigg{(}\sum_{p=1}^{\infty}\binom{2p}{p}\Big{(}\frac{t(1-t)}{\lambda
y^{2}}\Big{)}^{p}\bigg{)}\,dt.$

Now, the generalised binomial theorem tells us

(4.26)
$\sum_{p=0}^{\infty}\binom{2p}{p}z^{p}=\frac{1}{\sqrt{1-4z}},\qquad|z|<1/4.$

To use this, we note that the condition

(4.27) $\Big{|}\frac{t(1-t)}{\lambda y^{2}}\Big{|}<\frac{1}{4},\qquad
t\in(0,1-e^{-\lambda})$

is satisfied for

(4.28)
$|y|>\frac{1}{\sqrt{\lambda}}\begin{cases}2\sqrt{(1-e^{-\lambda})e^{-\lambda}}&\lambda<\log
2,\vskip 3.0pt plus 1.0pt minus 1.0pt\\\ 1&\lambda\geq\log 2.\end{cases}$

Consequently,

(4.29) $\displaystyle\begin{split}G(y)&=\frac{1}{y}+\frac{1}{\lambda
y}\int_{0}^{1-e^{-\lambda}}\frac{1}{1-t}\bigg{(}\frac{\sqrt{\lambda}y}{\sqrt{\lambda
y^{2}-4t(1-t)}}-1\bigg{)}\,dt\\\ &=\frac{1}{\lambda
y}\int_{0}^{1-e^{-\lambda}}\frac{1}{1-t}\Big{(}1-\frac{4t(1-t)}{\lambda
y^{2}}\Big{)}^{-\frac{1}{2}}\,dt.\end{split}$

To proceed further, suppose that $|y|<1/\sqrt{\lambda}$ is given. Then there
is range of $t$ values

(4.30)
$t\in\Big{(}\tfrac{1}{2}-y^{*},\min\\{1-e^{-\lambda},\tfrac{1}{2}+y^{*}\\}\Big{)},\qquad
y^{*}=\frac{1}{2}\sqrt{1-\lambda y^{2}}$

for which

(4.31) $1-\frac{4t(1-t)}{\lambda y^{2}}<0.$

By the Sokhotski–Plemelj formula applied to (4.24), we have

(4.32) $\widehat{\rho}^{\rm(dH)}(y)=\frac{1}{\pi}\lim_{\epsilon\to 0^{+}}{\rm
Im}\,G(y-i\epsilon).$

Making use of (4.29) and (4.31) then gives

(4.33) $\widehat{\rho}^{\rm(dH)}(y)=\frac{1}{\pi}\frac{1}{\lambda
y}\int_{\tfrac{1}{2}-y^{*}}^{\min\\{1-e^{-\lambda},\tfrac{1}{2}+y^{*}\\}}\frac{1}{1-t}\Big{(}1-\frac{4t(1-t)}{\lambda
y^{2}}\Big{)}^{-\frac{1}{2}}\,dt.$

It is possible to evaluate this integral. First, it follows from (4.33) that
if $\lambda<\log 2$,

(4.34) $\displaystyle\widehat{\rho}^{\rm(dH)}(y)$
$\displaystyle=\frac{1}{\pi\lambda
y}\int_{\tfrac{1}{2}-y^{*}}^{1-e^{-\lambda}}\frac{1}{1-t}\Big{(}\frac{4t(1-t)}{\lambda
y^{2}}-1\Big{)}^{-\frac{1}{2}}\,dt\,\mathbbm{1}_{(-a(\lambda),a(\lambda))}(y),$

while if $\lambda\geq\log 2$,

(4.35)
$\displaystyle\begin{split}\widehat{\rho}^{\rm(dH)}(y)&=\frac{1}{\pi\lambda
y}\int_{\tfrac{1}{2}-y^{*}}^{1-e^{-\lambda}}\frac{1}{1-t}\Big{(}\frac{4t(1-t)}{\lambda
y^{2}}-1\Big{)}^{-\frac{1}{2}}\,dt\,\mathbbm{1}_{(-a(\lambda),a(\lambda))}(y)\\\
&\quad+\frac{1}{\pi\lambda
y}\int_{\tfrac{1}{2}-y^{*}}^{\tfrac{1}{2}+y^{*}}\frac{1}{1-t}\Big{(}\frac{4t(1-t)}{\lambda
y^{2}}-1\Big{)}^{-\frac{1}{2}}\,dt\,\Big{(}\mathbbm{1}_{(-1/\sqrt{\lambda},1/\sqrt{\lambda}}(y)-\mathbbm{1}_{(-a(\lambda),a(\lambda))}(y)\Big{)}.\end{split}$

Now, using computer algebra we can check that

$\int\frac{1}{1-t}(At(1-t)-1)^{-\frac{1}{2}}\,dt=\frac{1}{i}\log\bigg{(}\frac{1+\sqrt{A}(1-t)+\sqrt{1-At(1-t)}}{1-\sqrt{A}(1-t)-\sqrt{1-At(1-t)}}\bigg{)}.$

Letting $A=4/(\lambda y^{2})$, we write

$\displaystyle
F(t):=\frac{1}{i}\log\bigg{(}\frac{1+\frac{2}{\sqrt{\lambda}y}(1-t)+i\sqrt{\frac{4}{\lambda
y^{2}}t(1-t)-1}}{1-\frac{2}{\sqrt{\lambda}y}(1-t)-i\sqrt{\frac{4}{\lambda
y^{2}}t(1-t)-1}}\bigg{)}.$

Substituting in (4.35) it follows that

(4.36)
$\displaystyle\begin{split}\widehat{\rho}^{\rm(dH)}(y)&=\frac{1}{\pi\lambda
y}\Big{(}F(1-e^{-\lambda})-F(\tfrac{1}{2}-y^{*})\Big{)}\mathbbm{1}_{(-a(\lambda),a(\lambda))}(y)\\\
&\quad+\begin{cases}0&\textup{if }\lambda\leq\log 2,\vskip 3.0pt plus 1.0pt
minus 1.0pt\\\ \displaystyle\frac{1}{\pi\lambda
y}\Big{(}F(\tfrac{1}{2}+y^{*})-F(\tfrac{1}{2}-y^{*})\Big{)}\Big{(}\mathbbm{1}_{(-1/\sqrt{\lambda},1/\sqrt{\lambda}}(y)-\mathbbm{1}_{(-a(\lambda),a(\lambda))}(y)\Big{)}&\textup{if
}\lambda\geq\log 2.\end{cases}\end{split}$

The expression (4.36) can be simplified. Note that by definition of $y^{*}$,

$F(\tfrac{1}{2}-y^{*})=\frac{1}{i}\log\bigg{(}\frac{\sqrt{\lambda}y+1+\sqrt{1-\lambda
y^{2}}}{\sqrt{\lambda}y-1-\sqrt{1-\lambda y^{2}}}\bigg{)},\qquad
F(\tfrac{1}{2}+y^{*})=\frac{1}{i}\log\bigg{(}\frac{\sqrt{\lambda}y+1-\sqrt{1-\lambda
y^{2}}}{\sqrt{\lambda}y-1+\sqrt{1-\lambda y^{2}}}\bigg{)}.$

Therefore it follows that

(4.37) $F(\tfrac{1}{2}+y^{*})-F(\tfrac{1}{2}-y^{*})=\pi.$

On the other hand, we have

$\displaystyle F(1-e^{-\lambda})$
$\displaystyle=\frac{1}{i}\log\bigg{(}\frac{\sqrt{\lambda}y+2e^{-\lambda}+i\sqrt{4(1-e^{-\lambda})e^{-\lambda}-\lambda
y^{2}}}{\sqrt{\lambda}y-2e^{-\lambda}-i\sqrt{4(1-e^{-\lambda})e^{-\lambda}-\lambda
y^{2}}}\bigg{)}$
$\displaystyle=\frac{1}{i}\log\bigg{(}\frac{1+y\sqrt{\lambda}}{1-y\sqrt{\lambda}}\bigg{)}^{1/2}+\arg\bigg{(}\frac{\sqrt{\lambda}y+2e^{-\lambda}+i\sqrt{4(1-e^{-\lambda})e^{-\lambda}-\lambda
y^{2}}}{\sqrt{\lambda}y-2e^{-\lambda}-i\sqrt{4(1-e^{-\lambda})e^{-\lambda}-\lambda
y^{2}}}\bigg{)}.$

Note here that

$\displaystyle\log\bigg{(}\frac{1+y\sqrt{\lambda}}{1-y\sqrt{\lambda}}\bigg{)}^{1/2}-\log\bigg{(}\frac{\sqrt{\lambda}y+1+\sqrt{1-\lambda
y^{2}}}{\sqrt{\lambda}y-1-\sqrt{1-\lambda y^{2}}}\bigg{)}=\pi i.$

Therefore it follows that

(4.38)
$F(1-e^{-\lambda})-F(\tfrac{1}{2}-y^{*})=\arg\bigg{(}\frac{2e^{-\lambda}-\sqrt{\lambda}y+i\sqrt{4(1-e^{-\lambda})e^{-\lambda}-\lambda
y^{2}}}{2e^{-\lambda}+\sqrt{\lambda}y+i\sqrt{4(1-e^{-\lambda})e^{-\lambda}-\lambda
y^{2}}}\bigg{)}.$

Combining all of the above in (4.36), the proof is complete. ∎

###### Proof of Theorem 2.6.

By Theorem 2.3, we have

(4.39)
$\frac{1}{N}\int_{1/\sqrt{q}}^{1/\sqrt{q}}\widetilde{\rho}_{N}^{\rm(dH)}(x)\,d_{q}x=\mathcal{M}_{2p,0}^{\rm(dH)}+\frac{\mathcal{M}_{2p,1}^{\rm(dH)}}{N^{2}}+O(N^{-4}).$

The leading order follows from Proposition 4.4 since

(4.40)
$\widetilde{\rho}_{(0)}(x)=\frac{1}{\sqrt{\lambda}}\,\widehat{\rho}^{\rm(dH)}\Big{(}\frac{x}{\sqrt{\lambda}}\Big{)}.$

Write $\eta=e^{-\lambda}$. By Theorem 2.3,

(4.41) $\displaystyle\mathcal{M}_{2p,1}^{\rm(dH)}$
$\displaystyle=-\log\eta\bigg{(}\frac{p}{6}I_{1-\eta}(p+1,p)+\frac{(2p-1)!}{6(p-1)!^{2}}\eta^{p}(1-\eta)^{p-1}\Big{(}2+p-(2p+1)\eta\Big{)}\bigg{)}.$

We need to compute

(4.42)
$\displaystyle\widetilde{G}_{(1)}(x):=\frac{1}{x}+\frac{1}{x}\sum_{p=1}^{\infty}\frac{\mathcal{M}_{2p,1}^{\rm(dH)}}{x^{2p}}.$

Note that

(4.43)
$\frac{1}{1-\eta}\,\sum_{p=1}^{\infty}\frac{(2p-1)!}{6(p-1)!^{2}}\Big{(}\frac{\eta(1-\eta)}{x^{2}}\Big{)}^{p}\Big{(}2+p-(2p+1)\eta\Big{)}=\frac{\eta(1-\eta)(x^{2}-2\eta)}{2x^{4}}\Big{(}1-\frac{4\eta(1-\eta)}{x^{2}}\Big{)}^{-5/2}.$

We also have

(4.44)
$\displaystyle\begin{split}\sum_{p=1}^{\infty}\frac{p}{6}\frac{1}{x^{2p}}I_{1-\eta}(p+1,p)&=\frac{1}{6}\int_{0}^{1-\eta}\frac{1}{1-t}\bigg{(}\sum_{p=1}^{\infty}\frac{(2p)!}{(p-1)!^{2}}\Big{(}\frac{t(1-t)}{x^{2}}\Big{)}^{p}\bigg{)}\,dt\\\
&=\frac{1}{3x^{4}}\int_{0}^{1-\eta}t\Big{(}x^{2}+2t(1-t)\Big{)}\Big{(}1-\frac{4t(1-t)}{x^{2}}\Big{)}^{-5/2}\,dt.\end{split}$

Then

(4.45)
$\displaystyle\begin{split}\widetilde{G}_{(1)}(x)&=\frac{1}{x}-\lambda\frac{e^{-\lambda}(1-e^{-\lambda})}{2}(x^{2}-2e^{-\lambda})\Big{(}x^{2}-4e^{-\lambda}(1-e^{-\lambda})\Big{)}^{-5/2}\\\
&\quad-\frac{\lambda}{3}\int_{0}^{1-e^{-\lambda}}t\Big{(}x^{2}+2t(1-t)\Big{)}\Big{(}x^{2}-4t(1-t)\Big{)}^{-5/2}\,dt.\end{split}$

Note that $x^{2}-4t(1-t)<0$ if

(4.46)
$t\in\Big{(}\tfrac{1}{2}-x^{*},\min\\{1-e^{-\lambda},\tfrac{1}{2}+x^{*}\\}\Big{)},\qquad
x^{*}=\frac{1}{2}\sqrt{1-x^{2}}.$

Then the Sokhotski–Plemelj (4.32) formula completes the proof. ∎

## References

  * [1] S. Alwhishi, R. S. Adigüzel and M. Turan, _On the limit of discrete $q$-Hermite I polynomials_, Commun. Fac. Sci. Univ. Ank. Ser. A1 Math. Stat. 68 (2019), 2272–2282.
  * [2] Y. Ameur, C. Charlier and J. Cronvall, _Free energy and fluctuations in the random normal matrix model with spectral gaps_ , arXiv:2312.13904.
  * [3] G. W. Anderson, A. Guionnet, and O. Zeitouni, _An introduction to random matrices_ , Cambridge University Press, Cambridge, 2009.
  * [4] T. Assiotis, B. Bedert, M.A. Gunes, and A. Soor, _On a distinguished family of random variables and Painlevé equations_ , Probab. Math. Phys. 2 (2021), 613–642.
  * [5] T. Assiotis, B. Bedert, M. Gunes and A. Soor, _Moments of generalized Cauchy random matrices and continuous-Hahn polynomials_ , Nonlinearity, 34 (2021), 4923–4943.
  * [6] T. Assiotis, J. P. Keating and J. Warren, _On the joint moments of the characteristic polynomials of random unitary matrices_ , Int. Math. Res. Not. 2022 (2022), 14564–14603.
  * [7] C. W. J. Beenakker _Random-matrix theory of quantum transport_ , Rev. Mod. Phys. 69 (1997), 731–808.
  * [8] A. Borodin and V. Gorin, _Lectures on integrable probability_. In: V. Sidoravicius and S. Smirnov, (eds) Probability and Statistical Physics in St. Petersburg, Proc. Sympos. Pure Math., vol. 91, pp. 155–214. Amer. Math. Soc., Providence, 2016.
  * [9] A. Borodin and E. Strahov, _Correlation kernels for discrete symplectic and orthogonal ensembles_ , Comm. Math. Phys. 286 (2009), 933–977.
  * [10] A. Borodin, V. Gorin and A. Guionnet, _Gaussian asymptotics of discrete $\beta$-ensembles_, Publ. math. IHES 125 (2017), 1–78.
  * [11] E. Brézin, C. Itzykson, G. Parisi and J. B. Zuber, _Planar diagrams_ , Comm. Math. Phys. 59 (1978), 35–51.
  * [12] A. Bufetov and A. Knizel, _Asymptotics of random domino tilings of rectangular Aztec diamonds_ , Ann. Inst. Henri Poincaré Probab. Stat. 54 (2018), 1250–1290.
  * [13] S.-S. Byun, _Harer-Zagier type recursion formula for the elliptic GinOE_ , arXiv:2309.11185.
  * [14] S.-S. Byun and P. J. Forrester, _Spectral moments of the real Ginibre ensemble_ , Ramanujan J. (to appear), arXiv:2312.08896.
  * [15] S.-S. Byun, N.-G. Kang and S.-M. Seo, _Partition functions of determinantal and Pfaffian Coulomb gases with radially symmetric potentials_ , Comm. Math. Phys. 401 (2023), 1627–1663.
  * [16] P. Cohen, _Moments of discrete classical $q$-orthogonal polynomial ensembles_, arXiv:2112.02064.
  * [17] P. Cohen, F. D. Cunden and N. O’Connell, _Moments of discrete orthogonal polynomial ensembles_ , Electron. J. Probab. 25 (2020), Paper No. 72, 19 pp.
  * [18] S. Corteel, B. Jonnadula, J. P. Keating and J. S. Kim, _Lecture hall graphs and the Askey scheme_ , arXiv:2311.12761.
  * [19] S. Corteel, J. S. Kim and D. Stanton, _Moments of orthogonal polynomials and combinatorics_ , Recent trends in combinatorics, 545–578, IMA Vol. Math. Appl., 159, Springer, 2016.
  * [20] J. B. Conrey, M. O. Rubinstein and N. C. Snaith, _Moments of the derivative of characteristic polynomials with an application to the Riemann zeta function_ , Comm. Math. Phys. 267 (2006), 611–629.
  * [21] F. D. Cunden, A. Dahlqvist and N. O’Connell, _Integer moments of complex Wishart matrices and Hurwitz numbers_ , Ann. Inst. Henri Poincaré Comb. Phys. Interact. 8 (2021), 243–268.
  * [22] F. D. Cunden, F. Mezzadri, N. O’Connell and N. Simm, _Moments of random matrices and hypergeometric orthogonal polynomials_ , Comm. Math. Phys. 369 (2019), 1091–1145.
  * [23] D. S. Dean, P. Le Doussal, S. N. Majumdar and G. Schehr, _Non-interacting fermions in a trap and random matrix theory_ , J. Phys. A 52 (2019), 144006.
  * [24] P. Flajolet, _Combinatorial aspects of continued fractions_ , Discrete Math. 32 (1980), 125–161.
  * [25] P. J. Forrester, _Log-gases and random matrices_ , Princeton University Press, Princeton, NJ, 2010.
  * [26] P. J. Forrester, _Global and local scaling limits for the $\beta=2$ Stieltjes-Wigert random matrix ensemble_, Random Matrices Theory Appl. 11 (2022), 2250020.
  * [27] P. J. Forrester, _Differential identities for the structure function of some random matrix ensembles_ , J. Stat. Phys. 183 (2021), 33.
  * [28] P. J. Forrester, _Joint moments of a characteristic polynomial and its derivative for the circular $\beta$ ensemble_, Probab. Math. Phys. 3 (2022), 145–170.
  * [29] P. J. Forrester and S.-H Li, _Rate of convergence at the hard edge for various Pólya ensembles of positive definite matrices_ , Integral Transforms Spec. Funct. 33 (2022), 466–484.
  * [30] P. J. Forrester, S.-H. Li, B.-J. Shen and G.-F. Yu, _$q$ -Pearson pair and moments in $q$-deformed ensembles_, Ramanujan J. 60 (2023), 195–235.
  * [31] P. J. Forrester and A. Rahman, _Relations between moments for the Jacobi and Cauchy random matrix ensembles_ , J. Math. Phys. 62 (2021), 073302.
  * [32] Y. V. Fyodorov and P. Le Doussal, _Moments of the position of the maximum for GUE characteristic polynomials and for log-correlated Gaussian processes_ , J. Stat. Phys. 164 (2016), 190–240.
  * [33] M. Gisonni, T. Grava and G. Ruzza, _Jacobi ensemble, Hurwitz numbers and Wilson polynomials_ , Lett. Math. Phys. 111 (2021), no. 3, Paper No. 67, 38 pp.
  * [34] P. J. Hanlon, R. P. Stanley and J. R. Stembridge, _Some combinatorial aspects of the spectra of normally distributed random matrices_ , Contemp. Math. 138 (1992), 151–174.
  * [35] J. Harer and D. Zagier, _The Euler characteristic of the moduli space of curves_ , Invent. Math. 85 (1986), 457–485.
  * [36] U. Haagerup and S. Thorbjørnsen, _Random matrices with complex Gaussian entries_ , Expo. Math. 21 (2003), 293–337.
  * [37] B. A. Khoruzhenko, D. V. Savin, and H.-J. Sommers, _Systematic approach to statistics of conductance and shot-noise in chaotic cavities_ , Phys. Rev. B 80 (2009), 125301.
  * [38] M. E. H. Ismail, D. Stanton and G. Viennot, _The combinatorics of $q$-Hermite polynomials and the Askey-Wilson integral_, European J. Combin. 8 (1987), 379–392.
  * [39] K. Johansson, _Discrete orthogonal polynomial ensembles and the Plancherel measure_ , Ann. Math. 153 (2001), 259–296.
  * [40] B. Jonnadula, J. P. Keating and F. Mezzadri, _Symmetry function theory and unitary invariant ensembles_ , J. Math. Phys. 62 (2021), 093512.
  * [41] J. P. Keating and N. C. Snaith, _Random matrices and L-functions_ , J. Phys. A 36 (2003), 2859–2881.
  * [42] S. V. Kerov, _Rooks on Ferrers boards and matrix integrals_ , J Math Sci 96 (1999), 3531–3536.
  * [43] R. Koekoek, P. A. Lesky and R. F. Swarttouw, _Hypergeometric Orthogonal Polynomials and Their $q$-Analogues_, Springer Monographs in Mathematics, Springer-Verlag Berlin Heidelberg, 2010.
  * [44] M. Ledoux. _Differential operators and spectral distributions of invariant ensembles from the classical orthogonal polynomials. The continuous case_ , Electron. J. Probab. 9 (2004) 177–208.
  * [45] M. Ledoux, _Differential operators and spectral distributions of invariant ensembles from the classical orthogonal polynomials: the discrete case_ , Electron. J. Probab. 10 (2005), 1116–1146.
  * [46] M. Ledoux, _A recursion formula for the moments of the Gaussian orthogonal ensemble_ , Ann. Inst. Henri Poincaré Probab. Stat. 45 (2009), 754–769.
  * [47] S.-H. Li, B.-J. Shen, G.-F. Yu and P. J. Forrester, _Discrete orthogonal ensemble on the exponential lattices_ , arXiv:2206.08633.
  * [48] K. Liechty, _Non-intersecting Brownian motions on the half-line and discrete Gaussian orthogonal polynomials_ , J. Stat. Phys. 147 (2012), 582–622.
  * [49] M. Marino, _Lectures on non-perturbative effects in large N gauge theories, matrix models and strings_ , Fortsch. Phys. 62 (2014) 455–540.
  * [50] F. Mezzadri and N. Simm, _Moments of the transmission eigenvalues, proper delay times, and random matrix theory I_ , J. Math. Phys. 52 (2011), 103511.
  * [51] F. Mezzadri and N. Simm, _Moments of the transmission eigenvalues, proper delay times and random matrix theory II_ , J. Math. Phys. 53 (2012), 053504.
  * [52] F. Mezzadri and N. Simm, _$\tau$ -function theory of quantum chaotic transport with $\beta=1,2,4$_, Comm. Math. Phys. 324 (2013), 465.
  * [53] F. Mezzadri, A.K. Reynolds and B. Winn, _Moments of the eigenvalue densities and of the secular coefficients of $\beta$-ensembles_, Nonlinearity 30 (2017), 1034–1057.
  * [54] A. Morozov, A. Popolitov and Shakirov, _Quantization of Harer-Zagier formulas_ , Phys. Lett. B 811 (2020), 135932.
  * [55] P. Njionou Sadjang, _Moments of classical orthogonal polynomials_ , Ph.D. Thesis, Universität Kassel, 2013.
  * [56] M. Novaes, _Time delay statistics for finite number of channels in all symmetry classes_ , Europhys. Lett. 139 (2022), 21001.
  * [57] G. Olshanski, _Macdonald-level extension of beta ensembles and large- $N$ limit transition_, Comm. Math. Phys. 385 (2021), 595–631.
  * [58] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark., eds. _NIST Handbook of Mathematical Functions_ , Cambridge: Cambridge University Press, 2010.
  * [59] E. Onofri, _SU(N) Lattice gauge theory with Villain’s action_ , Nuovo Cim. A 66 (1981), 293–318.
  * [60] A. Rahman and P. J. Forrester, _Linear differential equations for the resolvents of the classical matrix ensembles_ , Random Matrices Theory Appl. 10 (2021), 2250003.
  * [61] S. Sedighi, A. Taherpour, S. Gazor and T. Khattab, _Eigenvalue-based multiple antenna spectrum sensing: Higher order moments_ , IEEE Trans. Wireless Commun. 16 (2016), 1168–1184.
  * [62] E. C. Titchmarsh, _Theory of functions_ , Oxford University Press, London, 1939.
  * [63] A. M. Tulino and S. Verdu _Random matrix theory and wireless communications_ , Foundations and Trends in Communcations and Information Theory, vol. 1, Now Publisher, 2004, pp. 1–182.
  * [64] P. S. Venkataramana, _On $q$-analogs of some integrals over GUE_, arXiv:1407.79602.
  * [65] X. G. Viennot, _A combinatorial interpretation of the quotient-difference algorithm_ , Formal Power Series and Algebraic Combinatorics (Moscow, 2000), pages 379–390. Springer, Berlin, 2000.
  * [66] B. Winn, _Derivative moments for characteristic polynomials from CUE_ , Comm. Math. Phys. 315 (2012), 532–562.
  * [67] E. Wigner, _Characteristic vectors of bordered matrices with infinite dimensions_ , Ann. of Math. 62 (1955), 548–564.
  * [68] N. Witte and P. J. Forrester, _Moments of the Gaussian $\beta$ ensembles and the large $N$ expansion of the densities_, J. Math. Phys. 55 (2014), 083302.