Level Compressibility of Certain Random Unitary Matrices

Abstract

The value of spectral form factor at the origin, called level compressibility, is an important characteristic of random spectra. The paper is devoted to analytical calculations of this quantity for different random unitary matrices describing models with intermediate spectral statistics. The computations are based on the approach developed by G. Tanner for chaotic systems. The main ingredient of the method is the determination of eigenvalues of a transition matrix whose matrix elements equal the squared moduli of matrix elements of the initial unitary matrix. The principal result of the paper is the proof that the level compressibility of random unitary matrices derived from the exact quantisation of barrier billiards and consequently of barrier billiards themselves is equal to $1/2$ irrespective of the height and the position of the barrier.

1. Introduction

The leading idea behind statistical descriptions of complex deterministic quantum problems is that quantum characteristics (e.g., eigenenergies) of a large variety of such problems are so erratic and irregular that their precise values are irrelevant (such as the position of a molecule in the air), and only their statistical properties are of importance. As matrices are inherent in quantum mechanics, random matrices occupy a predominant place in the application of statistics to quantum problems [1]. In a typical setting, one tries to find a random matrix ensemble whose eigenvalues have the same statistical distributions as (high-excited) eigenenergies of a given deterministic quantum problem. Until now, this query has been figured out only for two limiting classes of quantum problems: (i) models whose classical limit is integrable [2] and (ii) models whose classical limit is chaotic [3]. For generic integrable models, quantum eigenenergies are distributed as eigenvalues of diagonal matrices with independent identically distributed (i.i.d.) elements which means that their correlation functions after unfolding coincide with the ones of the Poisson distribution [2]. For generic chaotic systems, it was conjectured in [3] that their eigenenergies are distributed as eigenvalues of standard random matrix ensembles (GOE, GUE, GSE) depending only on system symmetries whose correlation functions are known explicitly [4]. The difference between these two cases is clearly seen from the limiting behaviour of their nearest-neighbour distribution $P_0\left(s\right)$, which is the probability density that two eigenvalues are separated by a distance s and there are no other eigenvalues in between. For the Poisson statistics, there is no level repulsion that means that ${lim}_{s\to 0}{P}_0\left(s\right)\ne 0$ and for large argument $P_0\left(s\right)$ decreases exponentially with s. Standard random matrix ensembles levels repel each other, ${lim}_{s\to 0}{P}_0\left(s\right)=0$, and $P_0\left(s\right)\sim exp(-a{s}^2)$ when $s\to \infty$.

These two big conjectures form a cornerstone of quantum chaos and have been successfully applied to various problems from nuclear physics to number theory. Nevertheless, they do not cover all possible types of models. Especially intriguing is the class of pseudo-integrable billiards (see, e.g., [5]) which are two-dimensional polygonal billiards whose angles ${\theta }_j$ are rational multiples of $\pi$

$${\theta }_j=\frac{m_j}{n_j}\pi$$

with co-prime integers $m_j$ and $n_j$. A peculiarity of such billiards is seen in the fact that their classical trajectories belong to a two-dimensional surface of genus g related with angles as follows [6]

$$g=1+\frac{N_n}{2}\sum _j\frac{m_j-1}{n_j}$$

where $N_n$ is the least common multiple of all denominators $n_j$. Consequently, any such model with at least one numerator $m_j>1$ is neither integrable (which would imply that trajectories belong to a two-dimensional torus with $g=1$) nor fully chaotic (in which case trajectories should cover a three-dimensional surface of constant energy), and the aforementioned conjectures cannot be applied to such systems. Numerical calculations show that for many pseudo-integrable billiards, the spectral statistical properties of corresponding quantum problems differ from both the Poisson statistics and the random matrix statistics mentioned above (see, e.g., [7,8] and references therein). In particular, for these models, (i) ${lim}_{s\to 0}{P}_0\left(s\right)=0$ as for standard random matrix ensembles but (ii) $P_0\left(s\right)\sim exp(-bs)$ for large s as for the Poisson statistics. Such hybrid statistics, labeled intermediate statistics, had been first observed in the Anderson model at the metal–insulate transition [9,10], and they constitute a special, interesting but poorly investigated class of spectral statistics.

Probably the simplest example of pseudo-integrable systems is the so-called barrier billiard, which is a rectangular billiard with a barrier inside sketched in Figure 1a. The quantum problem for this model consists of solving the Helmholtz equation

$$(\Delta +E_{\alpha }){\mathsf{\Psi }}_{\alpha }(x,y)=0$$

imposing that eigenfunction ${\mathsf{\Psi }}_{\alpha }(x,y)$ obeys (e.g.,) the Dirichlet conditions on the boundary of the rectangle as well as on the barrier

$${\mathsf{\Psi }}_{\alpha }(x,y){{|}_{\mathrm{boundary}}=0,{\mathsf{\Psi }}_{\alpha }(x,y)|}_{\mathrm{barrier}}=0.$$

Calculating the exact S-matrix for the scattering inside the infinite slab with a barrier depicted in Figure 1b, it has been demonstrated in [7,8] that spectral statistics of this model is the same as the statistics of eigenvalues of the following $N\times N$ random unitary matrix

$$B_{m,n}=e^{i{\mathsf{\Phi }}_m}\frac{L_m{L}_n}{x_m+x_n},m,n=1,\dots ,N$$

(1)

where ${\mathsf{\Phi }}_m$ are i.i.d. random variables uniformly distributed on interval $[0,2\pi )$ and

$$L_m=\sqrt{2x_m\prod _{k\ne m}\frac{x_m+x_k}{x_m-x_k}}$$

(2)

where coordinates $x_m$ depend on the position of the barrier.

Define the following quantities (momenta) of propagating modes in each of three channels indicated in Figure 1b.

$$p_m^{\left(1\right)}=\sqrt{k^2-\frac{{\pi }^2{m}^2}{b^2}},p_m^{\left(2\right)}=\sqrt{k^2-\frac{{\pi }^2{m}^2}{{(b-h)}^2}},p_m^{\left(3\right)}=\sqrt{k^2-\frac{{\pi }^2{m}^2}{h^2}}.$$

If the ratio $h/b$ is an irrational number, coordinates $x_m$ have the following form

$$\overrightarrow{x}=b\left(\underset{N_1}{\underbrace{p_1^{\left(1\right)},\dots ,p_{N_1}^{\left(1\right)}}},\underset{N_2}{\underbrace{-p_1^{\left(2\right)},\dots ,-p_{N_2}^{\left(2\right)}}},\underset{N_3}{\underbrace{-p_1^{\left(3\right)},\dots ,-p_{N_3}^{\left(3\right)}}}\right).$$

(3)

Here, $N_j$ with $j=1,2,3$ are the numbers of propagating modes in each channel

$$N_1=\left[\frac{kb}{\pi }\right],N_2=\left[\frac{k(b-h)}{\pi }\right],N_3=\left[\frac{kh}{\pi }\right]$$

(4)

where $\left[x\right]$ is the largest integer $\le x$ and the total dimension of the B-matrix is $N=N_1+N_2+N_3$.

When the ratio $h/b$ is a rational number, $h/b=p/q$ with co-prime integers p and q ($p<q$), there exist exact plane wave solutions of barrier billiard equal to zero at the whole line passing through the barrier. It is natural to disregard them and take into account only non-trivial eigenvalues. In such case, coordinates $x_m$ have to be chosen as indicated below

$$\overrightarrow{x}=b\left(\underset{k\ne 0modq}{\underbrace{p_1^{\left(1\right)},\dots ,p_k^{\left(1\right)},\dots ,p_{N_1}^{\left(1\right)}}},\underset{N_2}{\underbrace{-p_1^{\left(2\right)},\dots ,-p_{N_2}^{\left(2\right)}}},\underset{k\ne 0modp}{\underbrace{-p_1^{\left(3\right)},\dots ,-p_k^{\left(3\right)},\dots ,-p_{N_3}^{\left(3\right)}}}\right).$$

(5)

The dimension of this vector is

$$N^{\left(\mathrm{r}\right)}=N_1+N_2+N_3-2N_0,N_0=\left[\frac{kb}{\pi q}\right].$$

(6)

The matrix B can be generalised for arbitrary vector $x_m$ provided the following interlacing conditions are fulfilled

$$|x_1|>|x_2|>\dots >|x_N|,x_m={(-1)}^{m+1}\left|x_m\right|.$$

Exact correlation functions for the B-matrix are unknown at present. In [7,8], it was argued that an approximate Wigner-type surmise for this matrix corresponds to the so-called semi-Poisson distribution [11]. In particular, it implies that the probability density $P_n\left(s\right)$ that two levels are separated by a distance s and there are exactly n levels in between (after the standard unfolding) is given by the following expression

$$P_n\left(s\right)=\frac{2^{2n+2}}{(2n+1)!}{s}^{2n+1}{e}^{-2s},n\ge 0.$$

Numerical calculations presented in [7,8] agree with these simple formulas.

Despite the simplicity of the semi-Poisson distribution, it has been observed (mostly numerically) in various models ranging from certain pseudo-integrable models and quantum maps (see references in [7,8]) to the entanglement spectra of two-bits random many body quantum circuits [12].

This paper is devoted to the calculation of another important characteristic of spectral statistics, namely the level (or spectral) compressibility, which is a long-range characteristic of the spectral two-point correlation function. In particular, this quantity determines the limiting behaviour of the variance of the number of levels inside a given interval. More precisely, let $N\left(L\right)$ be the number of eigenvalues in an interval L unfolded to the unit mean density, which means that the mean number of levels in interval L equals L, $\langle N\left(L\right)\rangle =L$. By definition, the number variance is ${\Sigma }^{\left(2\right)}\left(L\right)\equiv \langle {(N\left(L\right)-L)}^2\rangle$. If for large L

$${\Sigma }^{\left(2\right)}\left(L\right)\underset{L\to \infty }{⟶}\chi L$$

(7)

constant $\chi$ is called the level compressibility. The importance of this quantity follows from the fact that for integrable systems with the Poisson statistics, $\chi =1$, but for standard random matrix statistics typical for chaotic models, $\chi =0$. For all examples of intermediate statistics, it was observed that $0<\chi <1$.

The conventional way of determining the level compressibility for dynamical systems is the summation over all periodic orbits in the diagonal approximation initiated in [13]. For the symmetric barrier billiard with $h/b=1/2$ and $d/a=1/2$, it has been demonstrated in [14] that $\chi =1/2$. The same value had been obtained in [15] for the case $h/b=1/2$ and arbitrary barrier height $d/a$. Finally, for $h/b=p/q$ with co-prime integers p, q and irrational values of $d/a$, it has been proven in [16] that

$$\chi =\frac{1}{2}+\frac{1}{q}$$

but in the calculations, exact eigenvalues whose eigenfunctions are zero on the whole line $y=h$ (see Figure 1a) have not been excluded. When these trivial eigenvalues are removed, the answer is $\chi =1/2$ [17]. Therefore, direct (and quite tedious) calculations suggest that for (almost) all positions and heights of the barrier, the level compressibility is the same as for the semi-Poisson statistics [11]

$$\chi =\frac{1}{2}$$

(8)

but the reason of this universality remains obscure.

The purpose of the paper is to find analytically the spectral compressibility for barrier billiards and for a few other models directly from the corresponding random unitary matrices. To achieve the goal, it is convenient to slightly generalise the method developed in [18] for random unitary matrices that appeared in the quantisation of quantum graphs [19]. The method is briefly explained in Section 2. In Section 3, this method is applied to random unitary matrices derived in [20] by the quantisation of a simple interval-exchange map. In this case, the transition matrix is a circulant matrix whose eigenvalues are known explicitly. The results coincide with the exact level compressibility for these models obtained in [21,22,23,24]. This example gives credit to the method and permits explaining its main features without unnecessary complications. In particular, it clarifies the situation (not covered by Refs. [21,22,23,24]) when a parameter that entered the matrix takes an irrational value. Numerically, it has been observed [22] that in such case, the spectral statistics are well described by the ones of chaotic systems (GOE or GUE), although the Lyapunov exponent of the underlying classical map is always zero.

The main part of the paper is devoted to the derivation of the level compressibility for random matrices associated with barrier billiards (1). The calculations are more complicated, as no eigenvalues (except one) are known analytically. The simplest case of the symmetric barrier billiard with $h/b=1/2$ is investigated in Section 4. To obtain tractable expressions, a kind of paraxial approximation is developed, which permits controlling the largest terms. By using such approximation, the transition matrix is transformed into a Toeplitz matrix with a quickly decreasing symbol, which allows us to find its eigenvalues for large matrix dimension. The result of this section is that the level compressibility equals $1/2$ in accordance with periodic orbit calculations in [14,15].

In Section 5, random unitary matrices corresponding to barrier billiards with irrational ratio $h/b$ (with coordinates given by (3)) are considered. In this case, the transition matrix in paraxial approximation contains quickly oscillating terms and, consequently, has forbidden zones in the spectrum. In spite of that, one can argue that largest moduli eigenvalues are insensitive to fast oscillations and are determined solely by a matrix averaged over such oscillations. The final matrix is also a Toeplitz type, which allows analytical calculations, proving that the level compressibility is again equal to $1/2$.

Section 6 addresses the calculation of level compressibility in the most complicated case of barrier billiards with rational ratio $h/b=p/q\ne 1/2$. The computation is cumbersome, but in the end, one comes to the conclusion that the level compressibility remains equal $1/2$. In other words, the level compressibility of barrier billiards is universal (i.e., independent of the barrier position and its height) and coincides with the semi-Poisson prediction [11].

Section 7 gives a brief summary of the results. A few technical points are discussed in Appendix A, Appendix B and Appendix C.

2. Generalities

It is well known that the level compressibility (7) is related with the two-point correlation form factor $K\left(\tau \right)$ as follows

$$\chi =\underset{\tau \to 0}{lim}K\left(\tau \right).$$

(9)

For $N\times N$ random unitary matrices U, the form factor can, conveniently, be written in the following concise form (see, e.g., [18])

$$K\left(\tau \right)=〈\frac{1}{N}{\left|\mathrm{Tr}{U}^n\right|}^2〉,\tau =\frac{n}{N}$$

(10)

where the average is taken either on different realisations of random parameters, or over a small window of $\tau$, or both.

Unitary matrices considered in the paper all have the product form

$$U_{j,k}=e^{i{\mathsf{\Phi }}_j}{w}_{j,k}$$

(11)

where ${\mathsf{\Phi }}_j$ are i.i.d. random phases uniformly distributed between 0 and $2\pi$, and matrix $w_{j,k}$ is a fixed unitary matrix.

In [18] (see [19] for more detailed discussion), it was shown that for such unitary matrices, the averaging over random phases leads in the diagonal approximation to the following formula

$$K^{\left(\mathrm{diag}\right)}\left(\frac{n}{N}\right)=g\frac{n}{N}\mathrm{Tr}\left(T^n\right)$$

(12)

where matrix elements of matrix T, called below the transition matrix, are squared moduli of matrix elements of matrix U

$$T_{j,k}=|U_{j,k}{|}^2={\left|w_{j,k}\right|}^2.$$

(13)

For systems without time-reversal invariance, $g=1$, and for models with time-reversal invariance, $g=2$ (here, only cases of $g=1$ are considered). Due to the unitarily of matrix U, the T-matrix is a double stochastic matrix, ${\sum }_j{T}_{j,k}=1$ and ${\sum }_k{T}_{j,k}=1$, thus having the meaning of classical transition matrix.

Let ${\mathsf{\Lambda }}_{\beta }$ be eigenvalues of the T-matrix

$$T_{j,k}{u}_k\left(\beta \right)={\mathsf{\Lambda }}_{\beta }{u}_j\left(\beta \right),\beta =0,\dots ,N-1.$$

From (12), it follows that

$$K^{\left(\mathrm{diag}\right)}\left(\frac{n}{N}\right)=\frac{n}{N}\sum _{\beta =0}^{N-1}{\mathsf{\Lambda }}_{\beta }^n.$$

(14)

The unitarity imposes that one eigenvalue ${\mathsf{\Lambda }}_0=1$. Perron–Frobenius theorem states that all other eigenvalues $|{\mathsf{\Lambda }}_{\beta }|\le 1$.

Let the set of eigenvalues be ordered $1={\mathsf{\Lambda }}_0\ge {\mathsf{\Lambda }}_1\ge \dots \ge {\mathsf{\Lambda }}_{N-1}$. One has

$$K^{\left(\mathrm{diag}\right)}\left(\tau \right)=\tau \left(\mathrm{Tr}\left(T^n\right)\right)=\tau +\tau \sum _{\beta =1}^{N-1}{\mathsf{\Lambda }}_{\beta }^{N\tau },\sum _{\beta =1}^{N-1}{\mathsf{\Lambda }}_{\beta }^{N\tau }\le (N-1){\left|{\mathsf{\Lambda }}_1\right|}^{N\tau }.$$

(15)

It has been noted in [18] that if

$$\underset{N\to \infty }{lim}{\left|{\mathsf{\Lambda }}_1\right|}^N=0$$

(16)

then the second term in (15) goes to zero for all finite $\tau$ and $K^{\left(\mathrm{diag}\right)}\left(\tau \right)=\tau$ for small $\tau$. Consequently, it is reasonable to conjecture (as it has been done in [18]) that the whole spectral statistics of such matrices will be well described by standard random matrix formulas.

For the matrices discussed below, criterium (16) is not fulfilled. Instead, in all considered cases, the largest moduli eigenvalues of transition matrices have the form $1-\mathcal{O}(1/N)$. To calculate the level compressibility from (14), the summation over all such eigenvalues is performed analytically, and then, the limit $\tau \to 0$ is taken.

It is well known that the form factor is not a self-average quantity. It has strong fluctuations and necessarily requires a smoothing. There exist two different sources of fluctuations. The first is related with random phases in matrices (11). Equation (12) corresponds to the averaging over these random phases in the diagonal approximation. The second has its roots in non-smoothness of $K^{\left(\mathrm{diag}\right)}(n/N)$ for different n and could be removed by a smoothing over a small interval of $\tau$. (A trivial example is ${(-1)}^n$).

3. Interval-Exchange Matrices

This section is devoted to the calculation of the level compressibility for special unitary matrices derived in [20] by quantisation of a simple two-dimensional parabolic map. Slightly generalising their result [21], one can write these matrices in the following form

$$M_{n,m}=e^{i{\mathsf{\Phi }}_n}\frac{1-e^{2\pi i\alpha N}}{N\left(1-e^{2\pi i(n-m+\alpha N)/N}\right)},n,m=1,\dots ,N.$$

(17)

Here, $\alpha$ is a real parameter and ${\mathsf{\Phi }}_n$ are i.i.d. random variables uniformly distributed between 0 and $2\pi$ (the case with ‘time-reversal symmetry’ when ${\mathsf{\Phi }}_{N-n+1}={\mathsf{\Phi }}_n$ requires only multiplication of the formulas below by factor $g=2$, as indicated in (12)).

When $\alpha$ is a rational number $\alpha =p/q$ with co-prime integers p and q, the original classical map is an interval-exchange map and, as it was shown in detail in [21,22,23,24], spectral statistics of matrices (17) in the limit of $N\to \infty$ are unusual and peculiar. It appears that the limiting results depend on the residue of $pNmodq$ (when $pN\equiv 0modq$ matrix (17) have explicit eigenvalues not interesting for our purposes). For example, if $\alpha =1/5$, there are two possibilities: $N\equiv \pm 1$ mod 5 and $N\equiv \pm 2mod5$. In the first case, the nearest-neighbour distribution is

$$P_0\left(s\right)=\frac{5^5}{4!}{s}^4{e}^{-5s}$$

but for the second one, the exact result is different

$$P_0\left(s\right)=(a_2{s}^2+a_3{s}^3+a_4{s}^4+a_5{s}^5+a_6{s}^6)e^{-5s}$$

where $a_2=625/2-275\sqrt{5}/2$, $a_3=3125/2-1375\sqrt{5}/2$, $a_4=71875/48+33125\sqrt{5}/48$, $a_5=-15625/3+9375\sqrt{5}/4$, and $a_6=1015625/288-453125\sqrt{5}/288$.

Although the spectral correlation functions for $\alpha =p/q$ are different for different residues $pNmodq$, the calculations show that the spectral compressibility for all residues remains the same [21,22]

$$\chi =\frac{1}{q}.$$

(18)

Matrices (17) with irrational $\alpha$ were investigated numerically in [22], and it was observed that their spectral statistics are well described by standard random matrix ensembles (GOE or GUE). In particular, it implies that in such case

$$\chi =0.$$

(19)

Below, it is shown that values (18) and (19) can easily be recovered by the discussed method.

The transition matrix (13) for the discussed case has the form

$$T_{n,m}\equiv {\left|M_{n,m}\right|}^2=\frac{{sin}^2\left(\pi \alpha N\right)}{N^2{sin}^2\left(\pi (n-m+\alpha N)/N\right)}.$$

(20)

This is a circulant matrix, and its eigenvalues are simply the Fourier transforms of its matrix elements

$${\mathsf{\Lambda }}_{\beta }=\frac{{sin}^2\left(\pi \alpha N\right)}{N^2}\sum _{s=0}^{N-1}\frac{e^{-2\pi i\beta s/N}}{{sin}^2\left(\pi (s+\alpha N)/N\right)}.$$

(21)

Differentiating both sides of the identity with integer s on z

$$\sum _{m=0}^{N-1}{e}^{2\pi im(s+z)/N}=\frac{1-e^{2\pi iz}}{1-e^{2\pi i(s+z)/N}},$$

after straightforward transformations, one proves that

$$\frac{{sin}^2\left(\pi z\right)}{{sin}^2\left(\pi (s+z)/N\right)}=-2isin\left(\pi z\right)e^{-\pi iz}\sum _{m=0}^{N-1}\left(m-N\frac{e^{\pi iz}}{2isin\left(\pi z\right)}\right)e^{2\pi im(s+z)/N}$$

from which it follows that eigenvalues (21) are

$${\mathsf{\Lambda }}_{\beta }=\left(1-\frac{\beta }{N}\left(1-e^{-2\pi i\alpha N}\right)\right)e^{2\pi i\beta \alpha },\beta =0,1\dots ,N-1.$$

(22)

Notice that ${\mathsf{\Lambda }}_0=1$, ${\mathsf{\Lambda }}_{N-\beta }={\mathsf{\Lambda }}_{\beta }^{*}$ for $\beta =1,\dots ,N-1$.

Consider first the case of rational $\alpha =p/q$. Assume that $N\ne 0$ mod q and calculate the form factor from (14) separately for $n\equiv r$ mod q with $r=0,1,\dots ,q-1$ (i.e., $n=qt+r$ with integer t)

$$K\left(\frac{qt+r}{N}\right)=\frac{qt+r}{N}\left[1+2\mathrm{Re}\sum _{\beta =1}^{N/2}{\left(1-\frac{\beta }{N}\left(1-e^{-2\pi i\alpha N}\right)\right)}^{qt+r}{e}^{2\pi i\beta pr/q}\right].$$

Here, we do not take into account that when N is even the term with $\beta =N/2$ is real and has an additional factor $1/2$.

Notice that the phase depends only on r. For large N, one can put the first factor in the exponent and sum the geometric progression from 1 to infinity. The answer is

$$K_r\left(\tau \right)=\tau +2\tau \mathrm{Re}\left(\frac{e^{-\xi \tau +2\pi ipr/q}}{1-e^{-\xi \tau +2\pi ipr/q}}\right),\tau =\frac{qt+r}{N}$$

(23)

where $\xi =1-e^{-2\pi ipN/q}=1-e^{-2\pi ik/q}$ with $k=pNmodq$.

When $\tau \to 0$, all terms except the one with $r=0$ tend to zero, as they do not have a pole at $\tau =0$. The remaining term equals

$$\underset{\tau \to 0}{lim}{K}_0\left(\tau \right)=2\mathrm{Re}\left(\frac{1}{\xi }\right)=\frac{2\mathrm{Re}(1-e^{2\pi iz})}{2-2cos\left(2\pi z\right)}=1.$$

It means that after the averaging over random phases, the limiting value of the form factor $K_r\left(\tau \right)$ strongly depends on small changes of $\tau =(qt+r)/N$. For $r=0$, $K_0\left(\tau \right)$ tends to 1 for small $\tau \ll 1$ but for all other terms with $r=1,\dots ,q-1$$K_r\left(\tau \right)$ tends to 0. The difference between these different values of $\tau$ is very small, of the order of $1/N$. Therefore, after the averaging over any small (but finite) interval of $\tau$, one gets

$$\underset{\tau \to 0}{lim}\langle K\left(\tau \right)\rangle =\frac{1}{q}$$

which agrees with (18) obtained in [22] by a different method.

For illustration, the results of the direct calculation of the form factor for $\alpha =1/5$ and $N=399$ and $N=398$ are presented in Figure 2. First, eigenvalues of the matrix (17) were calculated numerically, and then using (10), the form factor for different n values was computed. The result is averaged over 1000 realisations of random phases. It is clear that, indeed, for different residues of n modulo 5, the results are different, and when $n\equiv 0$ mod 5, the form factor at small argument is close to 1, but for all other residues, it starts at 0. The average over all 5 residues begins at $1/5$, as expected.

Such a clear picture appears when the form factor is calculated at special values of $\tau$, $\tau =n/N$ with integer n. Computing it at arbitrary arguments leads to an irregular plot but, of course, the average curve remains unchanged.

Exactly the same formulas can be applied for an irrational value of parameter $\alpha$. In this case, one has

$$K\left(\tau \right)=\tau +2\tau \mathrm{Re}\left(\frac{e^{-\zeta \tau }}{1-e^{-\zeta \tau }}\right),\zeta =1-e^{-2\pi i\alpha N}-2\pi i\alpha N.$$

(24)

The exponent $\zeta =2{sin}^2\left(\pi \alpha N\right)+i\left[sin\left(2\pi \alpha N\right)-2\pi \alpha N\right]$ has a large imaginary part when $N\to \infty$. It means that the above expression is a strongly oscillated function of $\tau$. When averaged over a small interval of $\tau$, one obtains $K\left(\tau \right)=\tau$ as it should be for the ensemble of usual random matrices (GUE). This result follows without calculations from the fact that the average of all eigenvalues ${\mathsf{\Lambda }}_{\beta }^{\tau N}$ except $\beta =0$ equals zero as a consequence of rapidly changing phases. (For even N, the term with $\beta =N/2$ is real, but as it tends to zero at large N, its contribution is negligible).

Notice that criterion (16) for matrix (17) with irrational $\alpha$ is not fulfilled. Nevertheless, the spectral statistics of such a matrix is close to GUE statistics. This example illustrates a new mechanism for the appearance of random matrix statistics. The contribution of higher eigenvalues of the transition matrix (13) decreases not because of a gap between the first and the second eigenvalues as has been proposed in (16) but due to rapid oscillations for large matrix dimensions.

4. Symmetric Barrier Billiard with h/b = 1/2

The central problem of the paper is the determination of level compressibility for the B-matrices given by (1) and (2) by employing the method proposed in [18] and used in the previous section for matrices derived from the quantisation of an interval-exchange map. The simplicity of treatment of interval-exchange matrices comes from the fact that their transition matrices are circulant matrices whose eigenvalues are known exactly. For the B-matrices, calculations are more complicated, as there are no explicit formulas for eigenvalues of the corresponding transition matrix.

$$T_{m,n}={\left|B_{m,n}\right|}^2=\frac{L_m^2{L}_n^2}{{(x_m+x_n)}^2},m,n=1,\dots ,N.$$

(25)

This section is devoted to the investigation of the B-matrix corresponding to the symmetric barrier billiard with ratio $h/b=1/2$. In this case, $q=2$, $b-h=h$, and the second part of the vector $\overrightarrow{x}$ in (5) coincides with the third one. Now, trivial eigenfunctions can be removed by considering a desymmetrised rectangular billiard with height $h=b/2$ and imposing the Neumann boundary conditions for negative x and $y=h$. It is the equivalent of dropping the second part of vector (5) and taking coordinates $x_m$ as follows [7]

$$x_m={(-1)}^{m+1}b\sqrt{k^2-\frac{{\pi }^2{m}^2}{b^2}},m=1,\dots ,N,N=\left[\frac{kb}{\pi }\right].$$

(26)

Odd (resp., even) indices describe the first (resp., the third) part of vector (5).

The numerically calculated spectrum of the transition matrix in this case is presented in Figure 3a.

To calculate this spectrum (or, at least, the behaviour of largest-moduli eigenvalues) analytically, a kind of paraxial approximation has been developed. It is based on the fact that the main ingredient of matrices with intermediate statistics is a linear fall-off of matrix elements from the diagonal [25,26]. In the simplest setting, it means that

$$M_{m,n}\sim \frac{R_{m,n}}{m-n+\mathrm{const}.},m,n\gg 1,m-n=\mathcal{O}\left(1\right).$$

Therefore, it is natural to assume that the most important contributions come from the pole terms with $R_{m,n}\approx R_{m,m}$. This type of approximation can be done directly from the definition (2), as it is demonstrated in Appendix A. According to these results, the T-matrix in the paraxial approximation is a block matrix

$$T=\left(\begin{array}{cc}{0}_{\mathrm{o},\mathrm{o}}& t_{\mathrm{o},\mathrm{e}}\\ t_{\mathrm{o},\mathrm{e}}^T& 0_{\mathrm{e},\mathrm{e}}\end{array}\right),t_{\mathrm{o},\mathrm{e}}\equiv t_{2m-1,2n}=\frac{1}{{\pi }^2{(n-m+1/2)}^2}.$$

(27)

Here, subscripts ‘o’ and ‘e’ indicate odd and even indices, respectively.

It is instructive to obtain this answer without the knowledge of the exact B-matrix. One can achieve it by using the instantaneous approximation used in quantum mechanics when the interaction changes suddenly. In optics, such an approximation is analogous to the Fraunhofer diffraction. In the barrier billiard, it corresponds to the situation when a wave with large momentum quickly moving in a channel enters into another channel (cf., Figure 1b). In the instantaneous approximation, eigenfunctions in the new channel are just a re-expansion of the initial eigenfunctions into a complete set of eigenfunctions with correct boundary conditions inside the final channel.

Consider a normalised wave with the Neumann boundary conditions at $y=h=b/2$ and the Dirichlet ones at $y=0$

$${\psi }_{2m-1}^{\left(1\right)}(x,y)=\frac{2}{\sqrt{b}}sin\left(\frac{(2m-1)\pi }{b}y\right)exp\left(i{p}_{2m-1}^{\left(1\right)}x\right),x<0,0\le y\le b/2$$

propagating in the desymmetrised barrier billiard at negative x. When it penetrates into the region of positive x, it has to be expanded into correct waves propagating inside that region

$${\psi }_{2m-1}^{\left(1\right)}(x,y)=\sum _{n=1}^{\infty }{S}_{2m-1,2n}{\psi }_{2n}^{\left(3\right)}(x,y)$$

where ${\psi }_{2n}^{\left(3\right)}(x,y)$ are waves obeying the Dirichlet boundary conditions at $y=0$ and $y=h=b/2$

$${\psi }_{2n}^{\left(3\right)}(x,y)=\frac{2}{\sqrt{b}}sin\left(\frac{2\pi n}{b}y\right)exp\left(i{p}_n^{\left(3\right)}x\right),x>0,0\le y\le b/2.$$

Coefficients $S_{2m-1,2n}$ are the S-matrix for this process. In the paraxial approximation, they are calculated as follows (notice that in the paraxial approximation, $p_{2m-1}^{\left(1\right)}\approx p_n^{\left(3\right)}$)

$$\begin{array}{ccc}\hfill S_{2m-1,2n}^{(\mathrm{p}.\mathrm{a}.)}& =& \frac{4}{b}{\int }_0^{b/2}sin\left(\frac{\pi (2m-1)}{b}\right)sin\left(\frac{2\pi n}{b}y\right)dy=\hfill \\ & =& \frac{{(-1)}^{m+n}}{{\pi }^2(n-m+1/2)}+\frac{{(-1)}^{m+n}}{{\pi }^2(n+m-1/2)}.\hfill \end{array}$$

Taking into account only the pole term (and symmetry of the S-matrix), one obtains for the T-matrix exactly the same expression as (A1).

Thus, the transition matrix (27) is a block Toeplitz matrix. It is plain that its eigenvalues $\mathsf{\Lambda }=\pm \sqrt{{\lambda }_{\beta }}$ where ${\lambda }_{\beta }$ are eigenvalues of a $N_1\times N_1$ matrix (with $N_1=[N/2]$).

$$P_{m,n}=\sum _{k=1}^{N_1}{t}_{m,k}{t}_{n,k},t_{m,k}=\frac{1}{{\pi }^2{(k-m+1/2)}^2}.$$

Dominant contributions to the sum come from regions $k\sim m$ and $k\sim n$. Due to a quick decrease of the summands, the finite summation over k can safely be substituted in the limit $N_1\to \infty$ by the sum over all integer k

$$P_{m,n}\approx \sum _{k=-\infty }^{\infty }{t}_{m,k}{t}_{n,k}.$$

Using (A6), the necessary sum is easily calculated, and the result is

$$P_{m,n}=t_{m-n},t_0=\frac{1}{3},t_{r\ne 0}=\frac{2}{{\pi }^2{r}^2}.$$

(28)

This formula is valid when $m,n\gg 1$ and $m-n=\mathcal{O}\left(1\right)$.

Matrix (28) is a Toeplitz matrix with quickly decreasing matrix elements. It is well known that eigenvalues of the $N\times N$ Toeplitz matrix can be asymptotically calculated as follows (see, e.g., [27,28,29,30,31] and references therein)

$${\lambda }_{\beta }=f\left(\frac{\beta }{2N}+\mathcal{O}\left(\frac{1}{N}\right)\right),\beta =1,\dots ,N$$

(29)

where function $f\left(x\right)$ called the symbol is the Fourier series of $t_r$.

$$f\left(x\right)=\sum _{r=-\infty }^{\infty }{t}_r{e}^{2\pi irx}.$$

(More precise formulas can be found in the above references).

Using (A4) and (A5) one finds that the symbol of matrix (28) is

$$f\left(x\right)=\frac{1}{3}+\frac{2}{{\pi }^2}\underset{r=-\infty }{\overset{\infty }{{\sum }^{\prime }}}\frac{e^{2\pi irx}}{r^2}=\frac{1}{3}+4B_2\left(x\right)={(1-2x)}^2.$$

Therefore, eigenvalues of the P-matrix for large $N_1$ are

$${\lambda }_{\beta }\approx {\left(1-\frac{\beta }{N_1}\right)}^2,\beta =1,\dots ,N_1.$$

Eigenvalues of block matrix (27) $\mathsf{\Lambda }=\pm \sqrt{{\lambda }_{\beta }}$. Taking into account that the dimension of matrix (27) is $N\approx 2N_1$, one concludes that approximately its eigenvalues are

$$\mathsf{\Lambda }=\pm \left(1-\frac{2\beta }{N}\right),\beta =1,\dots ,\frac{1}{2}N.$$

(30)

With the corresponding redefinition of index $\beta$, these eigenvalues can be rewritten in the form

$${\mathsf{\Lambda }}_{\beta }\approx \frac{2\beta }{N}-1,\beta =1,\dots ,N$$

which agrees well with numerical calculations (see Figure 3a).

The form factor in the diagonal approximation is related with transition matrix eigenvalues by (14)

$$K^{\left(\mathrm{diag}\right)}\left(\frac{n}{N}\right)=\frac{n}{N}\sum _{\beta }{\mathsf{\Lambda }}_{\beta }^n.$$

As ${\mathsf{\Lambda }}_{N-\beta }=-{\mathsf{\Lambda }}_{\beta }$ for $\beta =1,\dots ,N-1$ (which is a consequence of the block structure of the transition matrix (27)), the form factor $K^{\mathrm{diag}}(n/N)$ with odd n in the diagonal approximation tends to zero when $\tau =n/N\to 0$

$$K^{\left(\mathrm{diag}\right)}\left(\tau \right)\underset{\tau \to 0}{⟶}0.$$

(31)

However, for even n, one obtains a different answer. Equation (30) may not be accurate for extreme eigenvalues with small $\beta$. For $\tau =2n/N$, one can separate the contribution of small $\beta <{\beta }_0$ and the rest for which (30) is a good approximation

$$K^{\left(\mathrm{diag}\right)}\left(\tau \right)=\tau \left(\mathrm{const}+2\sum _{\beta ={\beta }_0}^{N/2}{\left(1-\frac{2\beta }{N}\right)}^{\tau N}\right)\underset{N\to \infty }{⟶}\tau \left(\mathrm{const}+\frac{2e^{-2j_0\tau }}{1-e^{-2\tau }}\right)\underset{\tau \to 0}{⟶}1.$$

(32)

As has been discussed in the previous section, it means that the spectral compressibility of the B-matrix for symmetric barrier billiard coincides with the semi-Poisson value

$$\chi \equiv K\left(0\right)=\frac{1}{2}.$$

(33)

For illustration, the form factor for the symmetric billiard calculated numerically by the direct diagonalisation of $400\times 400$ matrices (1) with coordinates given by (26) and averaged over 1000 realisations is shown in Figure 3b. Two branches corresponding to odd and even n are clearly seen. The average over odd and even values agrees well with the semi-Poisson expression for the form factor [11] and, in fortiori, the level compressibility is $1/2$, as in (33).

5. Barrier Billiard with Irrational Ratio h/b

The transition matrices for a general barrier billiard with an off-centre barrier remain the same as in (25), but coordinates $x_m$ should have the form (3) for irrational ratio $h/b$ and (5) for rational $h/b=m/q$. The direct calculations of eigenvalues of these matrices reveal that they are more complicated that the ones for symmetric billiard with $h/b=1/2$ discussed in the previous section. As an example, in Figure 4, the spectra of the transition matrices with $h/b=1/\sqrt{5}$ and $h/b=2/5$ are presented. It is clearly seen that although eigenvalues with small moduli are quite irregular and have gaps, the largest moduli eigenvalues are well described by a straight line ${\mathsf{\Lambda }}_{\beta }=2\beta /N-1$.

This section is concentrated on the analytical treatment of billiards with irrational ratio $h/b$. As in the previous section, the first step consists of the calculation of a paraxial S-matrix for the scattering inside the slab with a barrier, as shown in Figure 1b. It can easily be completed in the instantaneous approximation exactly as above. In such an approximation, only transitions from channel 1 to 2 and to 3 and their inverse are non-zero. One has

$$\begin{array}{ccc}\hfill S_{m,n}^{1\to 3}& =& \frac{2}{\sqrt{bh}}{\int }_0^{h}sin\left(\frac{\pi m}{b}y\right)sin\left(\frac{\pi n}{h}y\right)dy=\hfill \\ & =& \frac{{(-1)}^{n}sin(\pi hm/b)}{\pi \sqrt{bh}}\left({(m/b-n/h)}^{-1}-{(m/b+n/h)}^{-1}\right).\hfill \end{array}$$

Similarly,

$$\begin{array}{ccc}\hfill S_{m,n}^{1\to 2}& =& \frac{2}{\sqrt{b(b-h)}}{\int }_h^{b}sin\left(\frac{\pi m}{b}y\right)sin\left(\frac{\pi n}{b-h}(b-y)\right)dy=\hfill \\ & =& \frac{{(-1)}^{n}sin(\pi mh/b)}{\pi \sqrt{b(b-h)}}\left({(m/b-n/(b-h))}^{-1}-{(m/b+n/(b-h))}^{-1}\right).\hfill \end{array}$$

The transition matrix $T={\left|S\right|}^2$ also has the same block structure. Retaining only the pole (the first) terms (and slightly changing the notations), one obtains

$$T=\left(\begin{array}{ccc}{0}_{1,1}& t_{n_1,n_2}^{1\to 2}& t_{n_1,n_3}^{1\to 3}\\ t_{n_1,n_2}^{1\to 2}& 0_{2,2}& 0_{2,3}\\ t_{n_1,n_3}^{1\to 3}& 0_{3,2}& 0_{3,3}\end{array}\right)$$

(34)

where $n_i=1,\dots ,N_i$ with $N_i$ given by (4) and

$$t_{n_1,n_2}^{1\to 2}=\frac{z_2{sin}^2(\pi n_{1}h/b)}{{\pi }^2{\left(z_2{n}_1-n_2\right)}^2},t_{n_1,n_3}^{1\to 3}=\frac{z{sin}^2(\pi n_{1}h/b)}{{\pi }^2{\left(z{n}_1-n_3\right)}^2},z_2=1-z,z=\frac{h}{b}.$$

(35)

Due to the block structure of the transition matrix (34), it follows that its eigenvalues $\mathsf{\Lambda }$ are determined by the relation ${\mathsf{\Lambda }}^2={\lambda }_{\beta }$ where ${\lambda }_{\beta }$ are eigenvalues of $N_1\times N_1$ matrix

$$P_{m,n}=\sum _{k=1}^{N_2}{t}_{m,k}^{1\to 2}{t}_{n,k}^{1\to 2}+\sum _{k=1}^{N_3}{t}_{m,k}^{1\to 3}{t}_{n,k}^{1\to 3}.$$

For large matrix dimensions, the summation can be extended over all integer k, and the sums can be calculated explicitly by using (A7) from Appendix C. The results are

$$P_{m,m}=\frac{z^2+{(1-z)}^2}{3}(3-2{sin}^2\left(\pi mz\right)),$$

(36)

and for $m\ne n$

$$P_{m,n}=\frac{2\left({sin}^2\left(\pi zm\right)+{sin}^2\left(\pi nz\right)\right)}{{\pi }^2{(m-n)}^2}-\frac{2(1-2z)sin\left(\pi mz\right)sin\left(\pi nz\right)sin\left(\pi \right(m-n\left)z\right)}{{\pi }^{3}z(1-z){(m-n)}^3}.$$

(37)

This matrix is a combination of Toeplitz terms dependent on the difference $m-n$ and oscillating terms (which explains the existence of forbidden zones in its spectrum; see Figure 4a).

Due to the unitarity of the B-matrix, the exact transition matrix $T={\left|B\right|}^2$ has the largest eigenvalue equal to 1, whose corresponding eigenvector is $(1,1,\dots ,1)$. It is natural (and is confirmed by calculations) that eigenvectors of the P-matrix corresponding to large moduli eigenvalues are slowly varying functions. Consequently, all oscillating terms in (36) and (37) for large m and n could be ignored. These arguments lead to the following recipe of the next step of approximation. Put $m=n+r$ and average all matrix elements of the P-matrix over quickly changing phase $\pi nz$. The calculations are straightforward and

$$\langle P_{m,n}\rangle \equiv \underset{N\to \infty }{lim}\frac{1}{N}\sum _{n=1}^N{P}_{n+r,n}=f_{m-n}$$

(38)

where

$$f_0=\frac{2(z^2+{(1-z)}^2)}{3},f_{r\ne 0}=\frac{2}{{\pi }^2{r}^2}-\frac{(1-2z)sin\left(2\pi rz\right)}{2{\pi }^{3}z(1-z)r^3}.$$

Eigenvalues of such a matrix for large N are calculated by the Fourier transform of this symbol

$${\lambda }_{\beta }=\sum _{r=-\infty }^{\infty }{f}_r{e}^{2\pi irx},x=\frac{\beta }{2N_1}.$$

The necessary sums are expressed through the Bernoulli polynomials (A4), (A5), and the result is

$${\lambda }_{\beta }=\frac{2}{3}(z^2+{(1-z)}^2)+4B_2\left(x\right)-\frac{1-2z}{3z(1-z)}\left(B_3\left(\{x+z\}\right)-B_3\left(\{x-z\}\right)\right).$$

From the beginning, one can assume that $h<h-b$, i.e., $z=h/b<1/2$ (the case $h/b=1/2$ was discussed in Section 4). Then

$${\lambda }_{\beta }=\left\{\begin{array}{cc}1-4x+{\displaystyle \frac{x^2}{z(1-z)}},& 0\le x\le z\\ {\displaystyle \frac{1-2z+2z^2}{1-z}}-{\displaystyle \frac{2x(1-x)}{1-z}},& z\le x\le \frac{1}{2}\end{array}\right.,x=\frac{\beta }{2N_1}.$$

(39)

As eigenvalues of the block matrix (34) $\mathsf{\Lambda }=\pm \sqrt{{\lambda }_{\beta }}$, it follows that close to maximum value (i.e., with small $\beta$),

$$\mathsf{\Lambda }\approx \pm \left(1-\frac{\beta }{N_1}+\mathcal{O}\left(\frac{1}{N_1}\right)\right)\approx \pm \left(1-\frac{2\beta }{N}\right).$$

As in the calculation of the form factor (14) small moduli eigenvalues are irrelevant, one can ignore higher-order terms in the above expression, which gives the same expression as in (30). It means that the level compressibility of barrier billiards with irrational ratio $h/b$ has the same value as in the preceding sections $\chi =\frac{1}{2}$.

In Figure 5a, the above formulas are compared with the results of direct calculations for the P-matrix with $h/b=1/\sqrt{5}$. As has been demonstrated, the approximate expression (39) is tangent to the exact spectrum close to 1. The form factor computed numerically for the same ratio $h/b$ is presented in Figure 5b. The agreement with the above result is clearly seen.

6. Barrier Billiard with Rational Ratio h/b = p/q

The calculation of transition matrix eigenvalues when the ratio $h/b$ is a rational number can be completed by a similar method. An additional difficulty in such case is that one has to select special combinations of states in the second and the third channels to remove trivial eigenvalues equal to zero on the whole line passing through the barrier. It has been discussed in detail in [8] and briefly reviewed in Appendix B. Combining all terms together, one concludes that the transition matrix when $h/b=p/q$ with p and q being co-prime integers has the block form similar to (34) but with one more block

$$T=\left(\begin{array}{cccc}{0}_{1,1}& t_{n_1,n_2}^{1\to 2}& t_{n_1,n_3}^{1\to 3}& t_{n_1,n4}^{1\to 4}\\ t_{n_1,n_2}^{1\to 2}& 0_{2,2}& 0_{2,3}& 0_{2,4}\\ t_{n_1,n_3}^{1\to 3}& 0_{3,2}& 0_{3,3}& 0_{3,4}\\ t_{n_1,n_4}^{1\to 4}& 0_{4,2}& 0_{4,3}& 0_{4,4}\end{array}\right).$$

(40)

Here, indices $n_j=1,\dots N_j^{\left(\mathrm{r}\right)}$ have the following restrictions

$$n_1\not\equiv 0modq,n_2\not\equiv 0modq-p,n_3\not\equiv 0modp$$

and

$$N_1^{\left(\mathrm{r}\right)}=N_1-N_0,N_2^{\left(\mathrm{r}\right)}=N_2-N_0,N_3^{\left(\mathrm{r}\right)}=N_3-N_0,N_4^{\left(\mathrm{r}\right)}=N_0$$

(41)

with $N_1,N_2,N_3$ given by (4) and $N_0$ is determined by (A2) or (6). The total matrix dimension is $N^{\left(\mathrm{r}\right)}={\sum }_{j=1}^4{N}_j^{\left(\mathrm{r}\right)}=N_1+N_2+N_3-2N_0$ as in (6).

Matrices $t^{1\to 2}$ and $t^{1\to 3}$ are the same as in (35) and $t^{1\to 4}$ given by (A3) from Appendix B

$$t_{n_1,n_4}^{1\to 4}=\frac{{sin}^2(\pi mh/b)}{{\pi }^{2}p(q-p){(n_1/q-n_4)}^2}.$$

The eigenvalues of block matrix (40) $\mathsf{\Lambda }$ are $\mathsf{\Lambda }={\lambda }_{\beta }^2$ where ${\lambda }_{\beta }$ with $\beta =1,\dots ,N_1$ are eigenvalues of matrix (superscript (res), indicating that the matrix describes the resonance case $h/b=p/q$)

$$P_{m,n}^{\left(\mathrm{res}\right)}=\sum _{k\ne 0mod(q-p)}{t}_{m,k}^{1\to 2}{t}_{n,k}^{1\to 2}+\sum _{k\ne 0modp}{t}_{m,k}^{1\to 3}{t}_{n,k}^{1\to 3}+\sum _k{t}_{m,k}^{1\to 4}{t}_{n,k}^{1\to 4}.$$

Using an evident relation

$$\sum _{k\ne 0modr}f\left(k\right)=\sum _{k}f\left(k\right)-\sum _{k}f\left(rk\right)$$

and (A7), the above sums can be explicitly calculated.

The results are

$$\begin{array}{ccc}\hfill P_{m,m}^{\left(\mathrm{res}\right)}& =& \frac{1}{3q^2}\left(p^2+{(q-p)}^2\right)\left(3-2{sin}^2(\pi pm/q)\right)\hfill \\ & +& \frac{2{sin}^4(\pi pm/q)\left(3-2{sin}^2(\pi m/q)\right)}{3p(q-p)q^2{sin}^4(\pi m/q)}\hfill \end{array}$$

(42)

and when $m\ne n$

$$P_{m,n}^{\left(\mathrm{res}\right)}=\frac{b_{m,n}}{{(m-n)}^2}-\frac{c_{m,n}sin(\pi p(m-n)/q)}{{(m-n)}^3}-\frac{d_{m,n}sin(\pi (m-n)/q)}{{(m-n)}^3}$$

(43)

where

$$\begin{array}{ccc}\hfill b_{m,n}& =& \frac{2}{{\pi }^2}\left({sin}^2(\pi pm/q)+{sin}^2(\pi pn/q)\right)\hfill \\ & +& \frac{2{sin}^2(\pi pm/q){sin}^2(\pi pn/q)}{p(q-p){\pi }^2}\left(\frac{1}{{sin}^2(\pi m/q)}+\frac{1}{{sin}^2(\pi n/q)}\right),\hfill \\ \hfill c_{m,n}& =& \frac{2q(q-2p)}{{\pi }^{3}p(q-p)}sin(\pi pm/q)sin(\pi pn/q),\hfill \\ \hfill d_{m,n}& =& \frac{4q{sin}^2(\pi pm/q){sin}^2(\pi pn/q)}{{\pi }^{3}p(q-p)sin(\pi m/q)sin(\pi n/q)}.\hfill \end{array}$$

(44)

Although these expressions are indexed by integers m and n, this notation is symbolic. The point is that by construction, these integers cannot be arbitrary but have to be not divisible by q. Let us order such numbers and let $\nu \left(k\right)$ with $k=1,2,\dots ,$ be the $k^{\mathrm{th}}$ integer $\not\equiv 0modq$. Then, indices of matrix $P_{m,n}^{\left(\mathrm{res}\right)}$ have to be considered as follows: $m=\nu \left(j\right)$, $n=\nu \left(k\right)$ with $j,k=1,2,\dots ,N_1^{\left(\mathrm{r}\right)}$ with $N_1^{\left(\mathrm{r}\right)}$ defined in (41). In such notation, matrix $P^{\left(\mathrm{res}\right)}$ is $N_1^{\left(\mathrm{r}\right)}\times N_1^{\left(\mathrm{r}\right)}$ matrix

$$P_{m,n}^{\left(\mathrm{res}\right)}\equiv P_{\nu \left(j\right),\nu \left(k\right)}^{\left(\mathrm{res}\right)},j,k=1,\dots ,N_1^{\left(\mathrm{r}\right)}.$$

The next step, as in the previous section (cf., (38)), consists in the substitution instead of the above exact expressions, their mean values with fixed difference between the indices

$$\langle P_{m,n}^{\left(\mathrm{res}\right)}\rangle \left(r\right)=\underset{N\to \infty }{lim}\frac{1}{N}\sum _{n=1}^N{P}_{n+r,n}^{\left(\mathrm{res}\right)}$$

where both integers n and $n+r$ have to be not divisible by q.

According to (42) and (43), the $P^{\left(\mathrm{res}\right)}$ matrix is a mixture of functions depending explicitly on the differences of indices and certain coefficients depending on indices modulus q. Only the latter requires the explicit averaging. Using (A8)–(A11) from Appendix C, one obtains that

$$\begin{array}{ccc}& & \sum _{m=1}^{q-1}{P}_{m,m}=\frac{(p^2+{(q-p)}^2)(2q-3)}{3q^2}+\frac{2p^2(2q-3p)}{3(q-p)q^2},\hfill \\ & & \underset{m=1}{\overset{q-1}{{\sum }^{\prime }}}{b}_{m,m+r}=\frac{2q(q-p+1-2{sin}^2(\pi pr/q))}{{\pi }^2(q-p)},\hfill \\ & & \underset{m=1}{\overset{q-1}{{\sum }^{\prime }}}{c}_{m,m+r}=\frac{q^2(q-2p)cos(\pi pr/q)}{{\pi }^{3}p(q-p)},\underset{m=1}{\overset{q-1}{{\sum }^{\prime }}}{d}_{m,m+r}=\frac{q^{2}sin(2\pi pr/q)}{{\pi }^{3}p(q-p)sin(\pi r/q)}.\hfill \end{array}$$

Here, it is taken into account that $p/q<1/2$. The superscript ′ in these sums indicates that the term with $m+r\equiv 0modq$ is omitted. The latter condition implies that the number of independent terms equal $q-1$ if $r\equiv 0modq$ or $q-2$ otherwise. Finally, one obtains

$$\langle P_{m,n}^{\left(\mathrm{res}\right)}\rangle ={\alpha }_1{\delta }_{m,n}+f_{m-n}$$

(45)

with

$$f_{m-n}=\left\{\begin{array}{cc}{\displaystyle \frac{{\alpha }_2}{{(m-n)}^2}},& (m-n)\equiv 0modq\\ & \\ {\displaystyle \frac{{\alpha }_3}{{(m-n)}^2}}+{\displaystyle \frac{{\alpha }_{4}cos(2\pi pr/q)}{{(m-n)}^2}}+{\displaystyle \frac{{\alpha }_{5}sin(2\pi pr/q)}{{(m-n)}^3}},& (m-n)\equiv r\ne 0modq\end{array}\right.$$

where constants ${\alpha }_j$ are

$$\begin{array}{ccc}\hfill {\alpha }_1& =& \frac{(p^2+{(q-p)}^2)(2q-3)}{3q^2(q-1)}+\frac{2p^2(2q-3p)}{3(q-p)q^2(q-1)}{\alpha }_2=\frac{2q(q-p+1)}{{\pi }^2(q-p)(q-1)},\hfill \\ \hfill {\alpha }_3& =& \frac{2q}{{\pi }^2(q-2)},{\alpha }_4=\frac{2q}{{\pi }^2(q-p)(q-2)},{\alpha }_5=\frac{q^2(2p-q-2)}{2{\pi }^{3}p(q-p)(q-2))}.\hfill \end{array}$$

Although this matrix depends only on the difference of indices $m-n$, it is not a Toeplitz matrix as m and n are not arbitrary numbers but only integers not divisible by q. Nevertheless, one can argue that the largest eigenvalues for large matrix dimension can be calculated by a formula similar to Toeplitz matrices (which is a kind of variational method)

$${\lambda }_{\beta }=\frac{1}{N_1^{\left(\mathrm{r}\right)}}\sum _{j,k=1}^{N_1^{\left(\mathrm{r}\right)}}\langle P_{\nu \left(j\right),\nu \left(k\right)}^{\left(\mathrm{res}\right)}\rangle e^{2\pi ix(j-k)},x=\frac{\beta }{2N_1^{\left(\mathrm{r}\right)}}.$$

(46)

Here, as above, $\nu \left(k\right)$ is the $k^{\mathrm{th}}$ integer $\not\equiv 0modq$.

In Appendix C (see (A13)), it is shown that such a sum can be written as follows

$$\begin{array}{ccc}\hfill {\lambda }_{\beta }& =& {\alpha }_1+{\alpha }_2\sum _{u=-\infty }^{\infty }\frac{e^{2\pi ix(q-1)u}}{{\left(qu\right)}^2}+2\mathrm{Re}\sum _{t=1}^{q-1}\left(1-\frac{t}{q-1}\right)\times \hfill \\ & \times & \sum _{u=-\infty }^{\infty }{e}^{2\pi ix\left(\right(q-1)u+t)}\left(\frac{{\alpha }_3}{{(uq+t)}^2}+\frac{{\alpha }_{4}cos(2\pi pt/q)}{{(uq+t)}^2}+\frac{{\alpha }_{5}sin(2\pi pr/q)}{{(uq+t)}^3}\right).\hfill \end{array}$$

The first sum is calculated through the Bernoulli polynomial $B_2\left(x\right)$ (see (A4)). The last sums are expressed through two functions

$$G(x,r)=\sum _{k=-\infty }^{\infty }\frac{e^{2\pi ix(kq+r)}}{{(kq+r)}^2},F(x,r)=\sum _{k=-\infty }^{\infty }\frac{e^{2\pi ix(kq+r)}}{{(kq+r)}^3}.$$

The explicit expressions of these function can be obtained as follows.

Define one more function

$$g(x,r)=\sum _{k=-\infty }^{\infty }\frac{e^{2\pi ix(kq+r)}}{kq+r}.$$

By the differentiation over x, one has $G^{\prime }(x,r)=2\pi ig(x,r)$ and $F^{\prime }(x,r)=2\pi iG(x,r)$. As the differentiation of $g(x,r)$ over x gives the sum of $\delta$-function, it is plain that $g(x,r)$ is the piece-wise constant function in interval $[j/q,(j+1)/q]$. Using (A6), one gets

$$g(j/q+x,r)=\frac{\pi }{qsin(\pi r/q)}{e}^{i\pi r(2j+1)/q}.$$

Correspondingly, function $G(x,r)$ is a piece-wise linear function in the same intervals

$$G(j/q+x,r)=\frac{{\pi }^2}{q^2{sin}^2(\pi r/q)}{e}^{2\pi irj/q}+\frac{2{\pi }^{2}ix}{qsin(\pi r/q)}{e}^{i\pi r(2j+1)/q}.$$

(47)

In the same way, one proves that function $F(x,r)$ is a piece-wise quadratic function

$$\begin{array}{ccc}\hfill F(j/q+x,r)& =& \frac{{\pi }^{3}cos(\pi r/q)}{q^3{sin}^3(\pi r/q)}{e}^{2\pi irj/q}+\frac{2{\pi }^{3}ix}{q^2{sin}^2(\pi r/q)}{e}^{2\pi irj/q}-\hfill \\ & -& \frac{2{\pi }^3{x}^2}{qsin(\pi r/q)}{e}^{i\pi r(2j+1)/q}.\hfill \end{array}$$

(48)

In all these formulas, $j=0,\dots ,q-1$ and $0\le x\le 1/q$.

Combining all terms together, one finds

$$\begin{array}{ccc}& & {\lambda }_{\beta }={\alpha }_1+\frac{2{\pi }^2{\alpha }_2}{q^2}{B}_2\left(\left\{x(q-1)\right\}\right)+2\mathrm{Re}\sum _{t=1}^{q-1}\left(1-\frac{t}{q-1}\right)\times \hfill \\ & \times & e^{ixt/q}\left[\left({\alpha }_3+{\alpha }_{4}cos(2\pi pt/q)\right)G\left(\frac{x(q-1)}{q},t\right)+{\alpha }_{5}sin(2\pi pt/q)F\left(\frac{x(q-1)}{q},t\right)\right].\hfill \end{array}$$

The main interest for the calculation of the form factor is the behaviour of the largest eigenvalues for x close to zero. Using (47) and (48), one concludes that

$${\lambda }_{\beta }=C_0+C_1\frac{\beta }{N_1^{\left(\mathrm{r}\right)}},0\le \beta \ll N_1^{\left(\mathrm{r}\right)}.$$

Here

$$\begin{array}{ccc}\hfill C_0& =& {\alpha }_1+\frac{{\pi }^2{\alpha }_2}{3q^2}+\hfill \\ & +& \frac{2{\pi }^2}{q^2}\sum _{t=1}^{q-1}\left(1-\frac{t}{q-1}\right)\left(\frac{{\alpha }_3+{\alpha }_{4}cos(2\pi pt/q)}{{sin}^2(\pi t/q)}+\frac{\pi {\alpha }_{5}sin(2\pi pr/q)cos(\pi t/q)}{q{sin}^3(\pi t/q)}\right)\hfill \end{array}$$

and

$$C_1=-\frac{2{\pi }^2{\alpha }_2(q-1)}{2q^2}-\frac{2{\pi }^2(q-1)}{q^2}\sum _{t=1}^{q-1}\left(1-\frac{t}{q-1}\right)\left({\alpha }_3+{\alpha }_{4}cos(2\pi pt/q)\right).$$

The sum over residues is of the form

$$\sum _{t=1}^{q-1}\left(1-\frac{t}{q-1}\right)h\left(t\right)$$

and (as it is easy to check) in the considered case $h(q-t)=h\left(t\right)$. Therefore

$$\sum _{t=1}^{q-1}th\left(t\right)=\sum _{t=1}^{q-1}(q-t)h(q-t)=\sum _{t=1}^{q-1}(q-t)h\left(t\right)⟶\sum _{t=1}^{q-1}th\left(t\right)=\frac{q}{2}\sum _{t=1}^{q-1}h\left(t\right).$$

Consequently

$$\begin{array}{ccc}\hfill C_0& =& {\alpha }_1+\frac{{\pi }^2{\alpha }_2}{3q^2}+\hfill \\ & +& \frac{2{\pi }^2(q-2)}{q^2(q-1)}\sum _{t=1}^{q-1}\left(\frac{{\alpha }_3+{\alpha }_{4}cos(2\pi pt/q)}{{sin}^2(\pi t/q)}+\frac{\pi {\alpha }_{5}sin(2\pi pr/q)cos(\pi t/q)}{q{sin}^3(\pi t/q)}\right)\hfill \end{array}$$

and

$$C_1=-\frac{2{\pi }^2{\alpha }_2(q-1)}{2q^2}-\frac{2{\pi }^2(q-2)}{q^2}\sum _{t=1}^{q-1}\left({\alpha }_3+{\alpha }_{4}cos(2\pi pt/q)\right).$$

Using sums indicated in Appendix C and collecting all terms in the end, one finds that

$$C_0=1,C_1=-2.$$

This result signifies that the largest moduli eigenvalues of the transition matrix for the barrier billiard with rational ratio $h/b=p/q$ are (i) independent of values of integers p and q and (ii) have the same asymptotic expression as in (30) (taking into account that $N^{\left(\mathrm{r}\right)}\approx 2N_1^{\left(\mathrm{r}\right)}$)

$$\mathsf{\Lambda }=\pm \left(1-\frac{\beta }{N_1^{\left(\mathrm{r}\right)}}+\mathcal{O}\left(\frac{1}{N_1^{\left(\mathrm{r}\right)}}\right)\right)\approx \pm \left(1-\frac{2\beta }{N^{\left(\mathrm{r}\right)}}\right).$$

As it has been explained above, it implies that (i) the form factor for barrier billiards is different for odd and even n and (ii) the spectral compressibility is exactly equal $1/2$ for all positions of the barrier

$$\chi \equiv K\left(0\right)=\frac{1}{2}.$$

The numerical calculations exemplified in Figure 6 confirm well these results.

7. Summary

It is demonstrated that the method of calculation of the level of compressibility proposed by G. Tanner in [18] for chaotic systems can successfully be applied for intermediate statistics models. The criterium discussed in [18] states that if the difference between the dominant eigenvalue of the transition matrix $\mathsf{\Lambda }=1$ and the second in magnitude eigenvalue is big enough, then only the dominant eigenvalues contributes to the form factor, and one obtains the usual value of the form factor corresponding to standard random matrix ensembles. Notably, the level compressibility is zero.

For models considered in the paper, no individual transition matrix eigenvalues dominate, and one has to sum over many of them with moduli close to 1. Two types of random unitary matrices were investigated. The first corresponding to a quantisation of an interval-exchange map [20] has been discussed in detail in [21,22,23,24]. In particular, the values of the level compressibility were derived. The application of the transition matrix approach for this case serves first of all to check the validity of Tanner’s method for intermediate statistics models. It appears that interval-exchange matrices lead to circulant transition matrices whose eigenvalues are explicitly known and all necessary sums are easily estimated. In the end, one obtains the same values of the level compressibility as obtained in [21,22,23,24] but with much simpler and transparent calculations. An example is of a special interest. It corresponds to interval-exchange matrices with an irrational value of a parameter (which strictly speaking describes not an interval-exchange map but only a parabolic one). Numerically, it has observed in [22] that in such case, spectral statistics is of usual random matrix type (GOE or GUE depending on a symmetry) as for chaotic systems, which looks strange as the Lyapunov exponent of any parabolic map is zero. The transition matrix approach clearly indicates that although there is no dominant eigenvalue as has been discussed in [18], all eigenvalues except $\mathsf{\Lambda }=1$ for large matrix dimensions are so quickly oscillating that averaging over a small interval of the argument effectively removes their contributions, producing the standard random matrix result.

The main part of the paper contains the calculation of the level compressibility for random unitary matrices derived from the exact quantisation of barrier billiards in [7,8]. The importance of such matrices comes from the fact that they have the same spectral statistics as high-excited states of barrier billiards, which are the simplest examples of pseudo-integrable models for which very little is known analytically.

The barrier billiard transition matrices are more complicated than the ones for interval-exchange matrices. Their spectra contain forbidden zones, and their exact eigenvalues seems not to be accessible in closed form. Nevertheless, as the level compressibility requires the control only of largest moduli eigenvalues of the transition matrix, it is possible to find such eigenvalues for large matrix dimensions precisely. The main simplification comes from the fact that eigenvectors corresponding to largest moduli eigenvalues are slow oscillating functions. Therefore, quickly oscillating terms in matrix elements will give negligible contributions on these eigenvectors, and one can substitute instead of exact matrix elements their average over fast oscillations. The resulting matrices are simpler and permit finding their large moduli eigenvalues analytically. In the end, one proves that the level compressibility of barrier billiards for all positions and heights of the barrier is the same and equals $1/2$. This result strongly indicates that spectral statistics of the B-matrices associated with barrier billiards is universal (i.e., independent on the barrier position) and well described by the semi-Poisson distribution.