Finite-Range Coulomb Gas Models I: Some Analytical Results

TL;DR

This paper introduces finite-range Coulomb gas models that interpolate between Poisson and classical random matrix statistics, revealing new universality classes and providing insights into banded matrices and related dynamical systems.

Contribution

It generalizes Dyson's models by incorporating finite-range interactions, establishing a framework for new universality classes in random matrix theory.

Findings

Transition from Poisson to classical RMT statistics with increasing interaction range

Introduction of new universality classes in random matrix ensembles

Application to banded random matrices and dynamical systems

Abstract

Dyson has shown an equivalence between infinite-range Coulomb gas models and classical random matrix ensembles for the study of eigenvalue statistics. In this paper, we introduce finite-range Coulomb gas (FRCG) models as a generalization of the Dyson models with a finite range of eigenvalue interactions. As the range of interaction increases, there is a transition from Poisson statistics to classical random matrix statistics. These models yield new universality classes of random matrix ensembles. They also provide a theoretical framework to study banded random matrices, and dynamical systems whose matrix representation can be written in the form of banded matrices.

Equations

$$P(A)=C\exp (-Tr{A}^2/4v^2),$$

$$P(A)=C\exp (-\beta TrV(A)),$$

$$P(U)=C\exp (-\beta TrV(U)),$$

$$\frac{dM(\tau )}{d\tau }=-\beta V^{\prime }(M(\tau ))+\xi (\tau ).$$

$$\overline{{\xi }_{jk}^{(\gamma )}(\tau )}=0,$$

$$\overline{{\xi }_{ij}^{(\gamma )}(\tau ){\xi }_{kl}^{(\eta )}({\tau }^{\prime })}={\delta }_{\gamma \eta }({\delta }_{ik}{\delta }_{jl}+(2{\delta }_{\gamma 0}-1){\delta }_{il}{\delta }_{jk})\delta (\tau -{\tau }^{\prime }).$$

$$\overline{\delta M}=-\beta V^{\prime }(M)\delta \tau .$$

$$\overline{\delta M_{ij}^{(\gamma )}\delta M_{kl}^{(\eta )}}={\delta }_{\gamma \eta }\left[{\delta }_{ik}{\delta }_{jl}+(2{\delta }_{\gamma 0}-1){\delta }_{il}{\delta }_{jk}\right]\delta \tau .$$

$$\overline{\delta x_j}=\beta \left[-V^{\prime }(x_j)+\sum _{k\ne j}(x_k-x_j{)}^{-1}\right]\delta \tau =-\beta \frac{\partial W}{\partial x_j}\delta \tau ,$$

$$W=-\sum _{j<k}\log |x_k-x_j|+\sum _{j}V(x_j),$$

$$\overline{\delta x_j\delta x_k}=2{\delta }_{jk}\delta \tau .$$

$$\frac{\partial p}{\partial \tau }=\sum _{j=1}^N\left[\frac{{\partial }^{2}p}{\partial x_j^2}-\beta \frac{\partial }{\partial x_j}\left(p\frac{\partial W}{\partial x_j}\right)\right].$$

$$\displaystyle\begin{split}p(x_{1},\cdots,x_{N})={}&C\exp(-\beta W)\\ =&C\prod_{j<k}|x_{j}-x_{k}|^{\beta}\exp\Big{(}-\beta\sum_{j}V(x_{j})\Big{)}.\end{split}$$

$$\overline{\delta\theta_{j}}=\beta\left[-V^{\prime}(\theta_{j})+\frac{1}{2}\sum_{k\neq j}\cot\Big{(}\frac{\theta_{k}-\theta_{j}}{2}\Big{)}\right]\delta\tau=-\beta\frac{\partial W}{\partial\theta_{j}}\delta\tau,$$

$$W(\theta_{j})=-\sum_{j<k}\log\Big{|}\sin\Big{(}\frac{\theta_{k}-\theta_{j}}{2}\Big{)}\Big{|}+\sum_{j}V(\theta_{j}),$$

$$\overline{\delta {\theta }_j\delta {\theta }_k}=2{\delta }_{jk}\delta \tau .$$

$$\displaystyle\begin{split}p(\theta_{1},\cdots,\theta_{N})=&C\exp\big{(}-\beta W\big{)}\\ =&C\prod_{j<k}|\sin(\theta_{j}-\theta_{k})|^{\beta}\exp\Big{(}-\beta\sum_{j}V(\theta_{j})\Big{)}.\end{split}$$

$$W=-{\sum _{j<k}}^{\prime }\log |x_k-x_j|+\sum _{j}V(x_j).$$

$$p(x_1,\cdots ,x_N)=$$

$$=$$

$$p({\theta }_1,\cdots ,{\theta }_N)=C\exp (-\beta W),$$

$$W=-{\sum_{j<k}}^{\prime}\log\Big{|}\sin\Big{(}\frac{\theta_{k}-\theta_{j}}{2}\Big{)}\Big{|}+\sum_{j}V(\theta_{j}).$$

$$p(\theta_{1},\cdots,\theta_{N})=C{\prod_{j<k}}^{\prime}|\sin(\theta_{j}-\theta_{k})|^{\beta}\exp\Big{(}-\beta\sum_{j}V(\theta_{j})\Big{)}.$$

$$\rho (x_1)={\int }_0^{\infty }\cdots {\int }_0^{\infty }p(x_1,\cdots ,x_N)d{x}_2\cdots d{x}_N.$$

$$M_p$$

$$⟨F⟩={\int }_0^{\infty }\cdots {\int }_0^{\infty }F(x_1,\cdots ,x_N)p(x_1,\cdots ,x_N)d{x}_1\cdots d{x}_N,$$

$$M_p=$$

$$=$$

$$M_p=$$

$$=$$

$$C{\int }_0^{\infty }{x}^p\exp \left(-\frac{\beta V(x)}{\beta d+1}\right)dx=\frac{C\beta }{(p+1)(\beta d+1)}{\int }_0^{\infty }{x}_j^{p+1}V{}^{{}^{\prime }}(x_j)\exp \left(-\frac{\beta V(x)}{\beta d+1}\right),$$

Full text

1 Introduction

There has been extensive use of random matrices in many branches of physics as well as in other disciplines. For example, they have found applications in quantum chaotic systems, most significantly complex nuclei, atoms, molecules and mesoscopic systems [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. In recent years, novel applications have emerged in biology [14, 15], economics [16, 17], and communication engineering [18, 19, 20]. In these applications, classical random matrix ensembles (uniform circular and Gaussian ensembles) have provided a framework for understanding complex spectra.

Dyson demonstrated that the joint probability distribution (jpd) of eigenvalues of these ensembles are equilibrium states of a Brownian motion model of Coulombic particles interacting via a logarithmic potential [21]. The positions of the Brownian particles are identified as the eigenvalues of the random matrix ensemble. In his original papers, Dyson considered (a) Coulombic particles moving on a real line with a harmonic binding (referred to as Gaussian ensembles (GE)), and (b) Coulombic particles moving on the unit circle (referred to as circular ensembles (CE)). In related work, Calogero and Sutherland [22] have demonstrated that the corresponding quantum Hamiltonians are exactly solvable.

In an important generalization, Dyson also considered non-harmonic confining potentials on the real line [23]. In this paper we will refer to these as linear ensembles, viz., ensembles of hermitian matrices. Similar generalizations can be done for circular ensembles, viz., ensembles of unitary matrices. The equilibrium properties of these non-Gaussian ensembles have been studied in detail [24, 25]. Similarly non-uniform circular ensembles have also been studied [26]. Further, non-equilibrium ensembles have also been studied extensively [26, 27, 28].

In the Dyson model and the above extensions, all particles have pairwise Coulombic interactions (i.e., each particle interacts with all other particles). In this paper, we consider the natural generalization to the case where particles have finite-range interactions. We will refer to these as finite-range Coulomb gas (FRCG) models. Such ensembles have important applications in the study of banded random matrices (BRM). Two specific examples are systems of quantum kicked rotors (QKR) and embedded Gaussian ensembles (EGE). Surprisingly, these generalizations have received very limited attention in the literature. In this (paper I) and its companion paper II, we present detailed analytical and numerical results for FRCG models. A brief account of papers I and II has been published in [29]. Some of the analytical results reported in this paper were obtained earlier by one of the authors [30] and have been used in [31, 32, 33].

This paper is organized as follows. In Sec. 2, we review Dyson’s Brownian motion models. In Sec. 3, we generalize the Dyson models to obtain FRCG models. In Secs. 4 and 5, we derive the level density of short-range models for linear and circular ensembles respectively. In Sec. 6, we derive the level density of long-range models for both linear and circular ensembles. In Sec. 7, we show the universality of spectral fluctuations with respect to different binding potentials in the circular and linear cases. In Sec. 8, we give some exact results for spectral properties of short-range models ($d=0,1$), where $d$ is a parameter quantifying the range of the interaction. In Sec. 9, we present a mean-field approximation for $d>1$. In Sec. 10, we present a detailed study of spectral properties for $d=2$ via an integral-equation approach. In Sec. 11, we describe the corresponding integral equation for $d>2$. We conclude with a summary and discussion in Sec. 12.

2 Dyson’s Brownian Motion Models

In this section we will introduce various models used in this paper. We consider $N$-dimensional random hermitian matrices $A$, which can be symmetric hermitian, complex hermitian, quaternion real hermitian or quaternion real self-dual. These are labeled by the Dyson parameter $\beta=1,2,4$ respectively [1, 6]. The matrices are $N$-dimensional in real, complex, and quaternion space. In Secs. 2.1, 2.2 and 2.3, we will review Dyson’s formulation of Brownian matrix ensembles.

2.1 Linear and Circular Ensembles

The joint probability density (jpd) of matrices $A$ in the case of Gaussian random matrix ensembles is given by

[TABLE]

where $C$ is the normalization constant (we will generally use $C$, $\tilde{C}$ etc. to denote the normalization constants for various distributions). Each matrix element $A_{jk}^{(\gamma)}$ has $\beta$ distinct ‘sites’, labeled by $\gamma=0,....,\beta-1$. The matrix elements of $A$ follow independent Gaussian distributions and have mean zero and variance $v^{2}$ for each of the $\beta$ distinct sites. The ensembles corresponding to Eq. (1) are Gaussian orthogonal ensemble (GOE), Gaussian unitary ensemble (GUE) and Gaussian symplectic ensemble (GSE) for $\beta=1,2,4$ respectively. A natural generalization of Eq. (1) is a non-Gaussian ensemble defined by

[TABLE]

where $V(A)$ is a positive-definite function of the matrix $A$.

In a similar fashion, one can define circular ensembles of unitary matrices $U$ (symmetric for $\beta=1$, general for $\beta=2$, and self-dual for $\beta=4$). The jpd of $U$ is

[TABLE]

where $V(U)$ is a positive-definite function of $U$ and $U^{{\dagger}}$, e.g., $V(U)=U+U^{{\dagger}}$. The case $V=0$ corresponds to the usual circular ensemble introduced by Dyson [21]. We will refer to this as the “uniform circular ensemble”, as opposed to $V\neq 0$ for the “non-uniform circular ensemble”. These ensembles are referred to as circular orthogonal ensemble (COE), circular unitary ensemble (CUE), and circular symplectic ensemble (CSE) for $\beta=1,2,4$, respectively.

2.2 Brownian Motion Model for Matrices

Dyson [21] introduced a Brownian matrix process $M(\tau)$ (over-damped case, i.e., Smoluchowski process [34]), in fictitious time $\tau$, yielding Eqs. (1) and (2) as equilibrium densities. For linear ensembles, $M(\tau)$ represents a Hamiltonian operator. Remarkably this leads to a Brownian process for the eigenvalues of $\{A\}\equiv\{M(\infty)\}$, which interact via a Coulomb gas (logarithmic) potential. In a similar fashion, a Brownian model can be written for the circular case, which yields the equilibrium density in Eq. (3) [21].

We review Dyson’s formulation, starting with the Langevin equation for the matrix variable $M(\tau)$ in the linear ensembles:

[TABLE]

Here, $\xi(\tau)$ is a Gaussian white-noise matrix which is similar in structure to $M$. Note that $M(\tau)$ is a matrix variable executing Brownian motion in the space of matrices $M$. The first two moments of the matrix elements $\xi_{jk}^{(\gamma)}$ are [21]

[TABLE]

[TABLE]

Here, the bar denotes the ensemble average.

We discretize the Langevin equation over an infinitesimal time interval $\delta\tau$. The matrix increment, $\delta M\equiv M(\tau+\delta\tau)-M(\tau)$, has a drift term (first term of RHS in Eq. (4)) and a diffusive term (second term of RHS in Eq. (4)). The first term contributes to the average and the second term contributes to the covariance of the matrix elements. The average of $\delta M$ is given by

[TABLE]

The second moments of elements of $\delta M$ are given by

[TABLE]

In the Gaussian case, $\delta M$ is an infinitesimal GOE, GUE, GSE respectively for $\beta=1,2,4$. The off-diagonal variances are $v^{2}\delta\tau$ in each case. The corresponding Fokker-Planck (FP) equation for the probability distribution of the matrix elements is given in [23]. The equilibrium jpd of $M(\infty)$ is the same as in Eq. (2).

For circular ensembles of matrices $U$, the infinitesimal increment is $\delta U=iS\delta MS^{D}$. Here, $S$ is unitary and $S^{D}$ is the transpose of $S$ for $\beta=1$, hermitian adjoint for $\beta=2$, and quaternion dual for $\beta=4$.

2.3 Brownian Motion Formulation for Eigenvalues

We consider the eigenvalues $\{x_{j};j=1,\cdots,N\}$ of matrices $M(\tau)$, which obey the Brownian process on the real line. The increment in eigenvalue $\delta x_{j}=x_{j}(\tau+\delta\tau)-x_{j}(\tau)$ for the matrix increment $\delta M$ can be computed by using second-order perturbation theory. Using Eqs. (7)-(8), the first two moments of $\delta x_{j}$ at fixed time $\tau$, correct upto first-order in $\delta\tau$ are

[TABLE]

where

[TABLE]

and

[TABLE]

(See Appendix C in [28] for the derivation of moments in Eqs. (9) and (11).) These moments are also known as “conditional moments”.

The Fokker-Planck equation for the jpd $p(x_{1},\cdots,x_{N})$ is

[TABLE]

The equilibrium jpd of eigenvalues is

[TABLE]

Here, $V(x)$= $x^{2}/4v^{2}$ gives the Gaussian ensemble results.

For circular ensembles, one deals with the eigenangles $\theta_{j}$ instead of real variables $x_{j}$. Using second-order perturbation theory for unitary matrices, Eqs. (9)-(11) are replaced by

[TABLE]

where

[TABLE]

and

[TABLE]

In Eqs. (14) and (15), the potential $V(\theta_{j})=V\left(e^{i\theta_{j}}\right)$.

The Fokker-Planck equation for the jpd $p(\theta_{1},\cdots,\theta_{N})$ is given by Eq. (12), with $x_{j}$ replaced by $\theta_{j}$ and $W$ as in Eq. (15). The equilibrium jpd is

[TABLE]

Note that the logarithmic potential terms in Eqs. (10) and (15) correspond to the $2$-dimensional Coulomb potential. This is the reason for the usage of the term “Coulomb gas”. In the next section, we will generalize the above formulation for the FRCG.

3 Finite-Range Coulomb Gas Models

We now introduce FRCG models as a generalization of the above FP equation for the eigenvalue dynamics. In the FP equation (12), we restrict the eigenvalue interaction to the range $d$, i.e., $W$ is given by

[TABLE]

Here ${\sum}^{\prime}$ denotes the sum over all $|j-k|\leq d$ with $j\neq k$. In Eq. (18), the logarithmic terms represent the finite-range repulsive two-dimensional Coulomb gas potential, and $V$ is a one-body binding potential. In equilibrium, the jpd of eigenvalues is

[TABLE]

where the prime denotes that the product is restricted to $|j-k|\leq d$. We will refer to these ensembles as linear ensembles. Note that we have considered the case of arbitrary potential $V(x)$ in the above discussion. The Gaussian case corresponds to $V(x)=x^{2}/{4v^{2}}$.

At this stage, it is appropriate to make some remarks about these FRCG models.

1. In the $d=0$ case, the interaction term in Eq. (18) is absent and the particles $\{x_{i}\}$ move independently. This corresponds to the Poisson limit. In the $d=N-1$ case, all particles interact with each other. This corresponds to the Wigner-Dyson classical ensembles. Thus, as $d$ changes from 0 to $N-1$, there is a crossover from the Poisson limit to the Wigner-Dyson limit. We would like to understand the nature of this crossover.

2. The term “finite range” refers to the range in eigenvalue indices, not actual distances. In principle, nearest-neighbor eigenvalues could lie far apart on the real number line.

3. It is relevant to ask whether the above crossover is realized in physical systems. In paper II, we will demonstrate that the FRCG models (and their extension to non-integer $d$) provide a framework to understand transitions in QKR and BRM.

4. The FRCG models were first proposed by one of the authors [30]. Their formal properties were studied by Pandey [30], Bogomolny et al. [31, 33], and Jain-Khare [32]. We will discuss connections to earlier work at appropriate places in papers I and II. The present work constitutes the first detailed application of FRCG models to study transitions in physical systems.

5. It is tempting to ask whether the FRCG models can be motivated from a Brownian matrix evolution with banded noise matrices. Then, the $d=N-1$ limit would correspond to the usual Dyson prescription in Secs. 2.2 and 2.3. We caution the reader that such a connection is not straightforward. At present, we will treat the FRCG models as generalizations of Dyson’s Brownian motion models for eigenvalue spectra. Clearly, an important direction for future work is to identify the matrix ensembles which yield the jpd of eigenvalues in Eq. (3).

The corresponding banded version of the circular ensemble follows again from Dyson’s prescription. The jpd of the resultant equilibrium ensemble is

[TABLE]

where $\{\theta_{j}\}$ are eigenangles in ascending order. The term $W$ is given by an appropriate generalization of Eq. (15):

[TABLE]

In this case, $V(\theta)$ is a potential periodic on the unit circle. We can also write $p$ as in Eq. (3):

[TABLE]

Note that $V(\theta)=0$ and $d=N-1$ corresponds to the usual Dyson’s circular ensemble, which we also refer to as the uniform circular ensemble. The introduction of $V$ renders it a non-uniform circular ensemble [26]. We will concentrate primarily on Gaussian ensembles in the linear case, and uniform ensembles in the circular case. However, we will consider some non-Gaussian and non-uniform circular ensembles also. We will derive analytic results in this paper, and supplement them with Monte Carlo (MC) calculations in paper II. In paper II, we will also discuss applications of FRCG ensembles.

4 Level Density: Linear Ensembles with $d=O(1)$

The level density in the linear case is defined by

[TABLE]

The corresponding $p^{\rm th}$ moment is given by

[TABLE]

Here, the angular brackets denote

[TABLE]

for a function $F(x_{1},\cdots,x_{N})$. After a partial integration, we get

[TABLE]

In the second step of the above equation, we have differentiated $W$ and written the double sum in the symmetrized form.

For $d=O(1)$, $x_{j}$ and $x_{k}$ in the double sum of Eq. (4) can be taken to be equal to each other. For a given $j$, $2d$ values of $k$ contribute to the sum. Thus

[TABLE]

In the last step of the above equation, $j$ is any of the $N$ indices. Now, using

[TABLE]

we identify the level density as

[TABLE]

Here, $C$ is the normalization constant. Note that, in the Gaussian ensemble, the density has a Gaussian form, but with a larger variance. In paper II, we will verify Eq. (29) by MC simulations. We will consider the following two potentials:

(a) Harmonic potential

[TABLE]

(b) Quartic potential

[TABLE]

In the harmonic and quartic cases, $\kappa$ sets the scale of the $V$-axis. In the quartic potential, $\alpha$ determines whether the potential is single well $(\alpha<0)$ or double well $(\alpha>0)$.

5 Level Density: Circular Ensembles with $d=O(1)$

For the circular ensembles, we follow the method described in the previous section. Consider the moment $M_{p}$ of the density $\rho(\theta)$ defined by

[TABLE]

Here, a partial integration has been used in the last step. Using $W$ from Eq. (21), and using angular brackets to denote the ensemble averages as in Eq. (25), we have

[TABLE]

In the first step, we have used

[TABLE]

In the second step of Eq. (5), symmetrization has been done in the double sum. Now, as in the linear case for $d=O(1)$, we take $\theta_{j}$ and $\theta_{k}$ to be equal in the double sum, leading to

[TABLE]

This yields

[TABLE]

In subsequent discussion, we will consider the following two potentials:

(a) Uniform potential: $V(\theta)=0$,

(b) Cosine potential: $V(\theta)=\kappa\cos(\theta)$.

6 Level Density: Linear and Circular Ensembles with $d=O(N)$

We first consider linear ensembles in the long-range case $(d=O(N))$. The typical pair of $x_{j}$ and $x_{k}$ in the last step of Eq. (4) can be taken to be independent. Thus, $M_{p}$ can be written as

[TABLE]

The number of distinct terms in ${\sum}^{\prime}$ can be shown to be $N_{d}=Nd-d(d+1)/2$. Thus, for $d=N-1$, we have $N_{d}=N(N-1)/2$ (i.e., all distinct pairs). For $d=(N-1)/2$, we have $N_{d}=(N-1)(3N-1)/8$. We assume that each term in the double sum in $M_{p}$ contributes equally. Thus, since there are $2N_{d}$ pairs of $(j,k)$, we obtain

[TABLE]

Here,

[TABLE]

with $d=\gamma N$. In long-range models,

[TABLE]

gives a density with finite support (i.e., independent of $N$). Dividing both sides of Eq. (6) by $N$ and ignoring $O(N^{-1})$ terms in the limit $N\rightarrow\infty$, we get

[TABLE]

Note that

[TABLE]

Comparing the final steps of Eqs. (41) and (42), we obtain

[TABLE]

where the principal value of the integral is implied. Integral equations of this type have been considered earlier; see for example [24, 35].

As a first example, we consider the above equation for the Gaussian case $u(x)=x^{2}/4$. Then, Eq. (43) yields the semi-circular density

[TABLE]

Notice that the semicircle radius $\sigma^{2}=1$ for $d=N-1$, which is the usual Wigner result. For $d\leq N-1$, the semicircle applies once again but with reduced radius. Consider also the Jacobi ensemble, with $V(x)=a\log(1-x)+b\log(1+x)$ where $a$ and $b$ are of order $1$. Since $V$ is independent of $N$, $u(x)=0$ for $|x|<1$. In this case, the density is independent of $\sigma^{2}$, given by

[TABLE]

For the quartic potential, $u(x)=(\kappa/N)(x^{4}/4-\alpha x^{2}/2)$. For $d=N-1$, it is known that the density undergoes a transition from one-band to two-band at the critical point $\alpha_{c}=\sqrt{2N/\kappa}$ [24]. For $d=O(N)$, the pre-factor becomes $\tilde{\kappa}=\kappa/N\sigma^{2}$, and $\alpha_{c}$ changes to $\tilde{\alpha_{c}}=\sigma\alpha_{c}$. The results for the density are

[TABLE]

and

[TABLE]

The support of these densities is single-band in Eq. (46) and two-band in Eq. (47) respectively.

For $d=O(N)$, with $V(\theta)=Nu(\theta)$ in the circular case, we take $\theta_{j}$ and $\theta_{k}$ to be independent as in the linear case, giving thereby

[TABLE]

Here, again, the principal value of the integral is implied. $\sigma^{2}$ is defined as in Eq. (39).

When $u(\theta)=0$, the density is uniform, given by

[TABLE]

Note that this simple result is a consequence of the circular symmetry, and is valid for all $d$ and $N$. This will be useful in computing the fluctuation results for large $d$.

For the non-uniform potential $u(\theta)=\kappa\cos(\theta)$, $\kappa$ changes to $\tilde{\kappa}=\kappa/\sigma^{2}$ for $d=O(N)$. The density is given by [26, 36]

[TABLE]

and

[TABLE]

The density in Eq. (50) corresponds to the banded case with a peak at $\theta=\pi$. On the other hand, the density in Eq. (51) describes the non-banded case, valid for the entire range of $\theta$ (again peaked at $\theta=\pi$).

Note that the critical point $\tilde{\alpha_{c}}$ decreases with $d$ in the quartic case, whereas the effective transition parameter $\kappa$ in the circular case increases with $\sigma^{2}$. Thus, in both the cases, the onset of the transition occurs more rapidly with decreasing $d$.

7 Equivalence of Fluctuations for Linear and Circular Ensembles

As we have mentioned earlier, the well-known classical ensembles (Gaussian linear as well as uniform circular) arise for $d=(N-1)$ [6]. For these ensembles, it has been proven analytically that the fluctuation properties are identical for each $\beta$ after proper unfolding of the spectra, i.e., there is a distinct universality class for each $\beta$. Further, these properties are also the same for an arbitrary smooth potential $V$ [29, 35, 37].

Therefore, it is natural to investigate the nature of universality for short-range models. In this section, we demonstrate that

(a) For $d=O(1)$: Each ($d,\beta$) gives rise to a distinct universality class.

(b) For $d=O(N)$: For each $\beta$, the corresponding universality class of classical ensembles applies.

Let us first focus on $d=O(1)$ cases for the circular ensemble. We define

[TABLE]

where $D(\theta)$ is the average spacing at $\theta$. The $\theta_{j}$’s are defined modulo $2\pi$. Note that the spacing $s_{j}$ has an average equal to $1$ for all $j$. This is commonly referred to as unfolding of the spectra [3]. $D(\theta)$ is given by

[TABLE]

where $\rho(\theta)$ is the normalized level density given in Eq. (36). For example, in the case of zero potential, $D(\theta)=2\pi/N$. Note that, for large $N$ and $j>k$,

[TABLE]

Thus, using Eq. (22), we obtain the jpd of the nearest neighbor spacings $s_{j}$ as

[TABLE]

valid for $N\gg d\geq 1$. This can also be written as

[TABLE]

The $\delta$-function term in Eqs. (55), (7) arises because the eigenvalues lie on the unit circle. Note that, in deriving Eqs. (55) and (7), we have used Eq. (36) to obtain

[TABLE]

Thus, the entire $V$-dependence is eliminated in the $\{\theta_{j}\}$ $\rightarrow$ $\{s_{j}\}$ transformation, including the contribution from the Jacobian of the transformation. Therefore, the fluctuation properties are independent of $V(\theta)$ for $d=O(1)$ in the circular ensemble.

In a similar way, for the linear ensembles we use Eqs. (52) and (53) with $\theta$ replaced by $x$. In this case, the density $\rho(x)$ is given by Eq. (29). Now, using Eq. (3), we obtain Eqs. (55) and (7) without the $\delta$-function term. For large $N$, the mean of $\sum s_{j}$ is $N$ with relative fluctuations of $O(N^{-1/2})$. Therefore, the $\delta$-function term can be reinstated. This proves the equivalence of circular and linear ensembles for the fluctuation properties for large $N$ with $d=O(1)$. This result is independent of the potential $V$.

Let us next turn to the case of $d=O(N)$. Eqs. (55) and (7) are valid for large $N$. To take the $N\rightarrow\infty$ limit, we need to consider the jpd of $n$ consecutive spacings where $n\ll N$. Thus,

[TABLE]

is the jpd for each $\beta$ and $d$, and independent of $V$. We will deal with explicit forms of $P^{(n)}_{d}$ in sections 8, 9 and 10.

The potential $V$ is important in determining the density; see Eqs. (29) and (36). On the other hand, the fluctuations on a scale $n\ll N$ are governed by the logarithmic potential. When $d=O(N)$ it is natural to expect that, after proper unfolding, the fluctuations on the scale of order $n$ will be independent of $d$ and therefore universal for each $\beta$.

The above discussion demonstrates the equivalence of linear and circular ensembles. Therefore, in our subsequent discussion, we will focus on the case of circular ensembles with potential $V(\theta)=0$.

We now introduce the statistical measures used to characterize eigenvalue fluctuations. These are as follows [3, 6].

1. $(n-1)^{\rm th}$ nearest neighbor spacing distribution $(n=1,2,3,\cdots)$

[TABLE]

A related quantity is

[TABLE]

which is the jpd of two separated spacings $s_{1}$ and $s_{n}$.

2. Spacing variance

[TABLE]

3. Two level correlation function $R_{2}$ and cluster function $Y_{2}$.

[TABLE]

4. Number variance

[TABLE]

Note that the results for $d=(N-1)$ with $N\rightarrow\infty$ correspond to the standard results given by Dyson, Mehta and others [1, 6]. In Wigner-Dyson statistics, all the results for eigenvalue and eigenvector fluctuations correspond to the $N\rightarrow\infty$ limit. In paper II, we will demonstrate from MC calculations that these results are applicable for $d=O(N)$. Here, we briefly review the main analytical results for fluctuation properties. These will be used for comparison with (a) results from numerical integration presented in this paper, and (b) the MC results presented in paper II.

The exact spacing distributions have complicated analytical forms deriving from integral equations. However, $p_{0}(s)$ is well-approximated by the following results [6]:

[TABLE]

The number variance for $r\gtrsim 1$ is given by

[TABLE]

where

[TABLE]

Here, $\gamma=0.5772156$ is the Euler constant. The spacing variance (for $\beta=1,2,4$) up to a good approximation is given by [38, 3]

[TABLE]

which is exact for $k\rightarrow\infty$.

The two-level cluster functions are known exactly [6]:

[TABLE]

Here,

[TABLE]

We will see in paper II that the Dyson-Mehta results already arise for FRCG models with large values of $d$, even for $d\ll N$.

8 Fluctuation Properties for $d=0,1$

We start with the simplest case, viz., $d=0$. From Eq. (58)

[TABLE]

In the first step of Eq. (8), the integration is facilitated by the change of variables

[TABLE]

and using $\delta(ax)=|a|^{-1}\delta(x)$. In the last step, the normalization constant is unity. This is the standard Poisson result as there is no interaction between the “particles”. The spacings $s_{j}$ are independent and have exponential distributions. Note that the results in this case are independent of $\beta$.

In this case, the fluctuation measures given in Eqs. (7)-(63) yield the following well-known results:

[TABLE]

and

[TABLE]

Next, we consider the $d=1$ case. Now, each particle (i.e., eigenvalue) interacts with its nearest neighbor. From Eq. (7) it is clear that there will be an extra factor $s_{j}^{\beta}$ in the jpd $P_{1}(s_{1}\cdots s_{N})$. Again, from Eqs. (55) and (58), we find

[TABLE]

In this case, the spacings are again independent, and

[TABLE]

Further, the $(n-1)^{\rm th}$ spacing distribution is given by

[TABLE]

The proof of this also follows directly from Eq. (7) by change of variables $s_{j}=sx_{j}$. The spacing variance is obtained as

[TABLE]

To obtain the number variance $\Sigma^{2}$, we first consider the Laplace transform of the two-level cluster function. Using Eqs. (62) and (82), and introducing $\alpha$ as the Laplace variable, we find

[TABLE]

Comparing the constant term, and the $\alpha$-order term, we have

[TABLE]

[TABLE]

Since the two $Y_{2}$ integrals in Eq. (63) converge rapidly, we can use Eqs. (85) and (86) in Eq. (63) to get

[TABLE]

The rapid convergence of integrals comes from the fact that the inverse Laplace transform gives exponential terms in $Y_{2}$. For example, for $\beta=1$ we find

[TABLE]

For other $\beta$-values, $Y_{2}(s)$ contains several exponential terms and is given in Eq. (95) for $d=1,\beta=2$.

9 Mean-Field Approximation for $d>1$

In this section, we present an approximation which reduces the arbitrary $d$ case to an effective $d=1$ case. This reduction is analogous to the mean-field (MF) approximation in statistical mechanics. The MF approximation will give good estimates for the above fluctuation measures. In Eq. (55), we set $s_{j+1},s_{j+2},\cdots,s_{j+d-1}\simeq s_{j}$ in each of the factors under the product sign, neglecting fluctuations in neighboring spacings. This yields

[TABLE]

Then, integrating over variables $s_{n+1},\cdots,s_{N}$ as in the $d=0,1$ cases, we obtain

[TABLE]

i.e., the spacings are independent for large $N$. This approximation is valid for $n\gtrsim d$. Note that Eq. (9) is exact for $d=0,1$.

Substituting Eq. (9) in Eq. (7), we obtain the spacing distribution for $n\geq 1$:

[TABLE]

The spacing variance is given by

[TABLE]

The leading term in Eq. (92) (i.e., the linear term in $n$) is determined by the MF result in Eq. (91). The constant term $\gamma(\beta,d)$ does not arise in the MF approximation, and its introduction is motivated by our exact $d=2$ result below. For $d>2$, we estimate this constant from our MC calculations as discussed in paper II.

To calculate the number variance, we use Eqs. (9) and (91) to obtain the Laplace transform of the two-level correlation function $R_{2}(s)$ in Eq. (62). We find for $\alpha\geq 0$,

[TABLE]

Here, the final expression results from the sum of an infinite geometric series. The cluster function $Y_{2}$ falls off exponentially for all $d=O(1)$. For example, for $d=1,\beta=2$ (i.e., $\xi=3$), we have from Eq. (9) after taking the inverse Laplace transform:

[TABLE]

This gives

[TABLE]

Using Eq. (9), we can evaluate the number variance $\Sigma^{2}(r)$ in Eq. (63) as

[TABLE]

For $d>2$, one can improve the MF approximation by considering a reduction to an effective $d=2$ case. This process gives an estimate for $\gamma(\beta,d)$ in Eq. (92), but the corresponding calculation requires solutions of the integral equation discussed below for $d=2$.

10 Fluctuation Properties for $d=2$

Let us now consider FRCG ensembles with next-nearest-neighbor interactions, i.e., $d=2$. This case will involve significant complications as compared to the $d=1$ case. This is due to the correlations introduced between the consecutive level spacings. From Eq. (55), the jpd of the nearest-neighbor spacings for $d=2$ is

[TABLE]

The corresponding jpd of $n\geq 2$ consecutive spacings is

[TABLE]

or

[TABLE]

Here, the exponential factor arises from the $\delta$-function. The factor $(2\beta+1)$ becomes $(\beta d+1)$ for general $d$, as there are $d$ multiplying factors of $s$-terms. The term $G(s_{1},s_{n})$ arises because there are two extra factors in Eq. (97) which depend on $s_{1}$ and $s_{n}$, viz., $(s_{1}+s_{N})^{\beta}$ and $(s_{n}+s_{n+1})^{\beta}$.

As $s_{1}$ and $s_{n}$ are well-separated by the integration variables $(s_{n+1},\cdots,s_{N})$, and the interaction is short-ranged, we assume that $G(s_{1},s_{n})$ is factorisable:

[TABLE]

From Eq. (99), a further integration on $s_{n}$ yields $P_{2}^{(n-1)}(s_{1},\cdots,s_{n-1})$. In this process, the function $F(s)$ does not change, but the constant changes. This implies that $F(s)$ is an eigenfunction of the integral equation:

[TABLE]

Here, $\lambda_{\mu}$ and $f_{\mu}$ are eigenvalues and eigenfunctions of this equation. The kernel $(t+s)^{\beta}$ is a polynomial of order $\beta$ in $t$, implying that there are $(\beta+1)$ eigenvalues. We order them as $\lambda_{0}<\lambda_{1}\cdots<\lambda_{\beta}$. The eigenfunction $f_{\mu}$ is a polynomial of order $\beta$, and this enables us to solve Eq. (101) numerically by matrix diagonalization. The eigenfunction $f_{\beta}$ corresponding to the largest eigenvalue will have only positive coefficients. Thus, for $f_{\beta}$ normalized as in Eq. (103) below, we obtain

[TABLE]

valid for $n=2,3,\cdots$.

It is useful to describe some properties of the integral equation (101). The kernel $(t+s)^{\beta}$ is symmetric in $t$ and $s$, so the integral equation is hermitian. Therefore, the eigenvalues $\lambda_{\mu}$ are real and the eigenfunctions $f_{\mu}$ are orthogonal. The eigenfunctions are appropriately normalized, giving the orthonormality relation:

[TABLE]

The properties of hermitian operators can be exploited to obtain sum rules involving the eigenvalues and eigenfunctions. The kernel $I(s,t)$ of the identity operator can be expressed as

[TABLE]

The spectral decomposition property can be written as

[TABLE]

This property can be generalized, using Eqs. (101) and (104), as

[TABLE]

Notice that, for $n=1$ in Eq. (10), we recover Eq. (105). Finally, we have the trace relation

[TABLE]

In Eq. (10), we integrate over all variables except $s_{1}$ and $s_{n}$ to obtain

[TABLE]

For $n=1$, this gives the nearest-neighbor spacing distribution

[TABLE]

We point out that the distributions in Eqs. (108) and (109) are normalized to unity. Further, because of unfolding, the average spacings are all unity.

For higher-order spacing distributions $(p_{n-1}$ or $p_{k}$ for $k=n-1,n\geq 1)$, we substitute Eq. (10) in Eq. (7) to obtain

[TABLE]

Thus, for example,

[TABLE]

We have performed a numerical integration of Eq. (10), and plot $p_{k}(s)$ for different values of $k$ for $d=2$ and $\beta=1,2,4$ in Fig. 1. We can also calculate $R_{2}(s)$ and $Y_{2}(s)$ from Eq. (62) with $p_{n-1}(s)$ from Eq. (10). We will show numerical results for $R_{2}(s)$ and $Y_{2}(s)$ for $d=2$ later.

It is convenient to define an average over level spacings:

[TABLE]

The variance of the nearest-neighbor spacing distribution is given by

[TABLE]

The covariance between $s_{1}$ and $s_{n}$ can be derived from Eq. (108) as

[TABLE]

valid for $n\geq 2$. Note that the covariance goes to zero exponentially as $n$ increases.

Using Eqs. (113) and (114) in Eq. (61), we obtain the variance of the $n^{\rm th}$ spacing distribution as

[TABLE]

The result in Eq. (10) consists of three terms: linear in $n$, constant, and exponentially decaying in $n$. In the large-$n$ limit, the spacing variance becomes

[TABLE]

where $\gamma(\beta,2)$ is a constant. This will be explicitly verified below for $\beta=1,2,4$. We will use the notation $\gamma(\beta,d)$ to denote this constant for general $d$.

For $\beta=1$, one can derive the above quantities explicitly. The two eigenvalues and their eigenfunctions are as follows:

[TABLE]

The other statistical quantities are

[TABLE]

and

[TABLE]

The equalities in Eqs. (120) and (121) apply for $n\geq 2$.

For $\beta=2,4$ we have numerically computed the eigenvalues and eigenfunctions from the integral equation (101). These are tabulated in the Appendix of this paper. Figs. 1-6 show different statistical measures obtained from numerical integrations for $d=2$ and $\beta=1,2,4$. In Fig. 1, we show the higher-order spacing distributions $p_{k}(s)$ for $\beta=1,2,4$. In Fig. 2, we make a comparison between the spacing distributions $p_{0},p_{1},p_{2}$ for all three $\beta$-values. In Fig. 3, we show a comparison of $p_{0}(s)$ for $d=0,1,2$ for each $\beta$. In Figs. 3(a), (b), (c), we have also included the corresponding classical-ensemble results for GOE, GUE and GSE, respectively. (As mentioned earlier, these arise for $d=N-1$.) These plots demonstrate the crossover from Poisson to classical results in the FRCG models as $d$ is increased from 0 to $N-1$. In Fig. 4, we show a comparison between our exact results and the MF approximation. We have compared results for $p_{k_{m}}(s)$ vs. $s$, where $k_{m}$ is the smallest value for which the exact results are numerically indistinguishable from the MF results on this scale. The quality of agreement improves further for higher values of $k$. Notice that $k_{m}$ is smaller for higher values of $\beta$. Fig. 5 and Fig. 6 correspond to the two-point correlation function $R_{2}(s)$ and the two-point cluster function $Y_{2}(s)$, respectively.

We have also used the eigenvalues and eigenfunctions to compute $\sigma^{2}(n-1)$ and verified Eq. (116) with $\gamma(2,2)=0.0451566$ and $\gamma(4,2)=0.0305985$. For $\beta=1,\gamma(1,2)=1/18$ as in Eq. (121).

Following the above discussion regarding $\sigma^{2}(n-1)$, and from the result for $\Sigma^{2}(n)$ in the $d=1$ case in Eq. (87), we rewrite the result of number variance for $d=2$ as

[TABLE]

This will be verified by MC calculations in paper II.

11 Fluctuation Properties for $d>2$

The $d>2$ case is analytically less tractable, as expected. However, the basic structure of many of the equations is formally similar to that for the $d=2$ case. Let us start with the jpd in Eq. (55) for $d>2$. We follow the calculations in Sec. 10. The first major change is encountered in the integral equation (101), which now becomes

[TABLE]

In this case, there are $\zeta+1$ eigenvalues and eigenfunctions, where $\zeta$ can be calculated from the multinomial form of $f_{\mu}$:

[TABLE]

Thus $\mu=0,1,\cdots,\zeta$. We consider the eigenvalues in ascending order with $\lambda_{\zeta}$ as the highest.

The appropriate generalization of Eq. (10) is

[TABLE]

This expression can be used to calculate various fluctuation properties. However, the subsequent calculations become more involved as $d$ increases. Thus, we use MC techniques (described in paper II) to obtain the complete picture for arbitrary $d$.

As we have mentioned in the introductory section, significant work on FRCG models was also done by Bogomolny et al. (BGS) [31, 33]. It is useful to compare and contrast our approach with that of BGS. The cases $d=0,1$ are simple, so we will focus on cases with $d\geq 2$.

1. We use the integral equation to directly calculate the jpds and the corresponding fluctuation properties. In this context, we adopt a scaling approach to simplify the complicated integral equations. This should be contrasted with Ref. [33], where BGS use a saddle-point approximation to make the equations tractable.

2. Due to the simplicity of our approach, we are able to explicitly calculate important physical quantities (inclusive of prefactors), e.g., level densities, spacing variance, number variance, spacing distributions, etc.

3. We present a MF approximation for $d\geq 2$, which enables the straightforward calculation of all statistical quantities for arbitrary $d$. Our MF approximation is validated by comparison with analytical results for $d=2$, and MC results for $d>2$ (presented in paper II).

4. In paper II, we present detailed MC results for the case with arbitrary $d$, enabling us to characterize the crossover from the Poisson limit to the Wigner-Dyson limit.

5. We present detailed numerical results in paper II to demonstrate that FRCG models provide a good framework to understand the spectral statistics of QKR and BRM.

6. In the context of point 5), we formulate FRCG models with fractional $d$ to cover the entire parameter range for QKR and BRM.

12 Summary and Discussion

In this paper, we have generalized Dyson’s Brownian motion model for eigenvalues of random matrix ensembles to introduce finite-range Coulomb gas (FRCG) models. These FRCG models are parametrized by the range $d$, which characterizes the extent of the particle-particle interaction in the Coulomb gas model. The FRCG models are solvable, and we have presented detailed analytical results. In this context, we calculate various fluctuation properties, e.g., spacing distributions, spacing variance, two-level correlation functions, etc. Further, we have presented an approximate mean-field (MF) solution which works very well for $d\geq 2$. This MF solution is validated by comparison with exact results for $d=2$.

From our analytical and numerical results, we observe that FRCG models provide an elegant route for transition from Poisson to classical ensembles as the interaction range increases from $d=0$ to $d=N-1$. Further, the onset of this transition is rapid. There are several interesting features of this crossover. For example, the number variance $\Sigma^{2}(r)$ exhibits a linear dependence on $r$ for FRCG models with $d=O(1)$, whereas this quantity grows logarithmically ($\sim\ln r$) for the classical ensembles. Moreover, $Y_{2}(s)$ decays exponentially (or faster) for small $d$, but decays algebraically for the classical ensembles. It is clearly relevant to quantify this crossover as a function of $d$, and we undertake this task in paper II of this two-part exposition.

In paper II, we present Monte Carlo (MC) results for FRCG models. The MC approach will provide “exact” results for FRCG models at intermediate values of $d$, where the framework described in paper I becomes unwieldy. Our results will confirm that FRCG models provide a new universality class of random matrix ensembles. In paper II, we will also study the applications of FRCG models. In this context, we will discuss banded random matrices (BRM) and quantum kicked rotors (QKR) in detail. We will also demonstrate that FRCG models are appropriate models for BRM and QKR.

Appendix

Eigenvalues and Eigenvectors for $d=2$, and $\beta=1,2,4$

[TABLE]

[TABLE]

[TABLE]