text:
# Periodic orbits contribution to the 2-point correlation form factor for pseudo-integrable systems

## 1 Introduction

The statistical properties of quantum systems attracted wide attention in the last years (see e.g. ). The investigation of many different models had led to a few accepted conjectures which relate statistical distribution of quantum energy levels with general properties of corresponding classical motion. For generic systems these conjectures are the following: for chaotic systems the level spacing distribution follows the Random Matrix statistics , ; for integrable systems it follows the Poisson statistics . Both conjectures are supported by a lot of numerical evidences and by some analytical arguments -.

These well-established conjectures are applicable only to completely chaotic or integrable models. But there are systems which are neither chaotic nor integrable. Noticeable examples of such systems are plane polygonal billiards with all angles, $`\alpha _i`$, commensurable with $`\pi `$

$$\alpha _i=\pi \frac{m_i}{n_i},$$
(1)
where $`m_i`$, $`n_i`$ are co-prime integers. In such systems all trajectories belong to a 2-dimensional surface of genus

$$g=1+\frac{N}{2}\underset{i}{}\frac{m_i1}{n_i},$$
(2)
where $`N`$ is the least common factor of the $`n_i`$ (see e.g. ). The case where all $`m_i=1`$ corresponds to $`g=1`$ (i.e. to a torus) which is integrable. If some $`m_i>1`$ trajectories belong to a higher genus surface and, consequently, the system is not integrable (at least in the usual sense) but it is not chaotic either since all trajectories belong to a 2-dimensional surface and cannot cover a 3-dimensional energy surface ergodically as is required for chaotic systems. For such reasons these systems are called pseudo-integrable.

A natural question appears: what is the spectral statistics of pseudo-integrable systems? Numerical calculations - clearly demonstrated that statistical properties of such systems differ from standard examples but have many points in common with the statistics of the 3-dimensional Anderson model at the metal-insulator transition point . The full analytical approach to this question meets with difficulties related mostly with the existence of quickly growing terms in the trace formula  which do not permit to use standard methods.

The main purpose of this paper is to compute analytically the value of the 2-point correlation form factor, $`K_2(\tau )`$, in the limit $`\tau 0`$ for two examples of pseudo-integrable systems, (i) a plain billiard in the shape of right triangle with one angle equal $`\pi /n`$ and (ii) a rectangular billiards with a Bohm-Aharonov flux line inside. We argue that in the small-$`\tau `$ limit the diagonal approximation  is valid and the problem reduces to the calculation of distributions of periodic orbit lengths and areas occupied by periodic orbit families. Though for general pseudo-integrable systems very little is known on this subject, triangular billiards in the shape of right triangles with angle $`\pi /n`$ belong to the so-called Veech polygons ,  and have a hidden group structure which make possible explicit calculation of necessary quantities. After the calculations we found a finite value of the 2-point correlation form factor at the origin, $`0<K_2(0)<1`$, which is different from both the Poisson distribution (for which $`K_2(0)=1`$) and the random matrix results (where $`K_2(0)=0`$). Analogous result has also been obtained for rectangular billiards with a Bohm-Aharonov flux line. Non-zero values of the 2-point correlation form factor at the origin confirm peculiar properties of spectral statistics for pseudo-integrable systems. We also discuss the comparison of theoretical predictions with the results of extensive numerical calculations.

The plan of the paper is the following. In Section 2 the discussions of the trace formula and the diagonal approximation for the 2-point correlation form factor are presented. A brief introduction to the Veech structure of certain pseudo-integrable billiards is given in Section 3. For clarity we start in Section 3.1 with a simple example of square billiards where ideas and methods can easily be illustrated. Needed properties of the modular group and the Eisenstein series are shortly revised in Sections 3.2 and 3.3. In Section 3.4 the Veech group for $`\pi /n`$ right triangle is derived and in Section 3.4.2 the density of periodic orbits for this triangle is computed. In Section 3.5 the calculation of 2-point form factor at the origin is performed and the comparison with the results of numerical calculations is discussed. Section 4 is dwelt on the calculation of the 2-point form factor for a rectangular billiard with a flux line. As in the previous Sections the main point is the calculation of areas swept by periodic orbits around the flux line. The result depends on arithmetical properties of ratios of coordinates of the flux line to the corresponding rectangular sides. In Section 5 concluding remarks are presented.

## 2 The form factor in the diagonal approximation

### 2.1 The density of states

The modern semiclassical approximation of multi-dimensional quantum systems is based on various types of trace formulas which express quantum density of states (and other quantities as well) through quantities computed in pure classical mechanics , , . The main step in deriving trace formulas is the semiclassical approximation for the (advanced) Green function

$$G_+(\stackrel{}{x},\stackrel{}{y})=\underset{n}{}\frac{\overline{\mathrm{\Psi }}_n(\stackrel{}{x})\mathrm{\Psi }_n(\stackrel{}{y})}{EE_n+iϵ},$$
(3)
where $`E_n`$ is the set of energy levels and $`\mathrm{\Psi }_n`$ the eigenfunctions as a sum over classical trajectories with energy $`E`$ connecting initial point $`\stackrel{}{x}`$ and final point $`\stackrel{}{y}`$ ,

$$G_+(\stackrel{}{x},\stackrel{}{y})=\underset{tr}{}A_{tr}\mathrm{exp}(\frac{i}{\mathrm{}}S_{cl}i\frac{\pi }{2}\nu ).$$
(4)
$`S_{cl}`$ is the classical action computed along a trajectory, $`A_{tr}`$ is a pre-factor depending on the system considered, and $`\nu `$ is a phase (the Maslov index) which, roughly speaking, counts points where simple semiclassical approximation breaks down.

For 2-dimensional free motion (and for 2-dimensional polygonal billiards) the semi-classical approximation for $`G`$ reads (see e.g. )

$$G_+(\stackrel{}{x},\stackrel{}{y})=\underset{p}{}\frac{e^{ikl_pi\frac{\pi }{2}\nu _pi\frac{3\pi }{4}}}{\sqrt{8\pi kl_p}},$$
(5)
where $`l_p`$ is the geometrical length of the orbit and $`k=\sqrt{E}`$ is the wave vector (in the units $`\mathrm{}=1`$ and $`m=1/2`$).

The knowledge of the Green function permits to find other quantum quantities as well. In particular the quantum density of states

$$d(E)=\underset{n}{}\delta (EE_n)$$
(6)
may be written by the means of the advanced Green function as

$$d(E)=\frac{1}{\pi }m𝑑\stackrel{}{x}G_+(\stackrel{}{x},\stackrel{}{x}).$$
(7)
The contribution from very short trajectories gives the mean level density, $`\overline{d}`$, and the integration over the space selects periodic orbit contributions ,  and determines an oscillating part of level density, $`d^{(osc)}(E)`$. For example, the density of states of an integrable rectangular billiard with sides $`a`$ and $`b`$ is

$$d(E)=\overline{d}+d^{(osc)}(E).$$
(8)
Here the smooth part is

$$\overline{d}=\frac{𝒜}{4\pi },$$
(9)
where $`𝒜`$ is the area of the rectangle (this formula is valid for all 2-dimensional billiards) and the oscillating part is

$$d_{p.o.}(E)=\underset{p.p.o.}{}\underset{n=1}{\overset{\mathrm{}}{}}\frac{𝒜_p}{4\pi }\frac{1}{\sqrt{2\pi knl_p}}e^{iknl_pi\frac{\pi }{2}\nu _pi\frac{\pi }{4}}+c.c.$$
(10)
where

$$l_p=\sqrt{(2Ma)^2+(2Nb)^2}.$$
(11)
In the rectangular billiard, the lengths of periodic orbits are 4 times degenerate in the sum $`(\text{5})`$ because $`(\pm M,\pm N)`$ give the same length. When the integral $`(\text{7})`$ is performed, orbits $`(M,N)`$ and $`(M,N)`$ are absorbed in the same $`𝒜_p`$. The summation in $`(\text{10})`$ is therefore performed over all primitive periodic orbits of length $`l_p`$ with $`M0`$ repeated $`n`$ times (an orbit $`(M,N)`$ and its time-reverse companion $`(M,N)`$ are counted as two different orbits). In all integrable billiards periodic orbits are not isolated but belong to families. $`𝒜_p`$ is the area of the pencil of periodic orbits of length $`l_p`$. For the rectangular billiard $`𝒜_p=2𝒜`$.

Pseudo-integrable systems considered in the paper belong to the class of diffractive systems whose characteristic property is the existence of singularities which make the classical motion undetermined. Each time a classical trajectory hits a singularity there is no unique way to continue it. Quantum mechanics smoothes out these singularities and associates with each (not too strong) singularity a diffraction coefficient, $`D(\stackrel{}{n},\stackrel{}{n}^{})`$, (or scattering amplitude) which defines an amplitude of scattering on this singularity from the initial direction $`\stackrel{}{n}`$ to the final direction $`\stackrel{}{n}^{}`$.

Correspondingly, the semiclassical approximation of the Green function in the presence of a singularity at point $`\stackrel{}{x}_0`$ takes the form

$$G(\stackrel{}{x},\stackrel{}{y})=G_0(\stackrel{}{x},\stackrel{}{y})+\underset{\stackrel{}{n},\stackrel{}{n}^{}}{}G_0(\stackrel{}{x},(\stackrel{}{x}_0,\stackrel{}{n}))D(\stackrel{}{n},\stackrel{}{n}^{})G_0((\stackrel{}{x}_0,\stackrel{}{n}^{}),\stackrel{}{y}),$$
(12)
where $`G_0(\stackrel{}{x},\stackrel{}{y})`$ is the Green function without singularity and $`G_0(\stackrel{}{x},(\stackrel{}{x}_0,\stackrel{}{n}))`$ is a contribution to the Green function from a classical trajectory starting at point $`\stackrel{}{x}`$ and ending at the singularity $`\stackrel{}{x}_0`$ with momentum in the direction $`\stackrel{}{n}`$. $`G_0((\stackrel{}{x}_0,\stackrel{}{n}^{}),\stackrel{}{x}^{})`$ is a contribution to the Green function from a classical trajectory starting at point $`\stackrel{}{x}_0`$ with momentum in the direction $`\stackrel{}{n}^{}`$ and ending at point $`\stackrel{}{x}^{}`$.

This modification of the Green function changes the trace formula. For diffractive systems the density of states can now be written as the sum of three terms , ,

$$d(E)=\overline{d}+d_{p.o.}(E)+d_{d.o.}(E),$$
(13)
where $`\overline{d}`$ is the mean level density, $`d_{p.o.}`$ is the contribution of periodic orbits without singularity, and the third term, $`d_{d.o.}(E)`$, is a contribution from all classical orbits starting and ending at the singularity (with, in general, different momenta). These trajectories are called diffractive orbits and $`d_{d.o.}(E)`$ is a sum over all possible combinations of them

$$d_{d.o.}(E)=\underset{m=1}{\overset{\mathrm{}}{}}\frac{1}{\pi m}\frac{}{E}G(\stackrel{}{n}_1,\stackrel{}{n}_1^{})D(\stackrel{}{n}_1^{},\stackrel{}{n}_2)\mathrm{}G(\stackrel{}{n}_{m1},\stackrel{}{n}_m^{})D(\stackrel{}{n}_m^{},\stackrel{}{n}_1),$$
(14)
where $`G(\stackrel{}{n},\stackrel{}{n}^{})`$ is the contribution to the Green function from a classical trajectory starting at the singular point with initial momenta in direction $`\stackrel{}{n}`$ and ending at it with final momentum in direction $`\stackrel{}{n}^{}`$.

For polygonal billiards the vertices with $`m_i1`$ play the role of singular point . In the case of scattering on the angle $`\alpha `$ the diffraction coefficient can be derived from Sommerfeld’s exact solution

$$D(\theta _f,\theta _i)=\frac{2}{\gamma }\mathrm{sin}\frac{\pi }{\gamma }\left[\frac{1}{\mathrm{cos}\pi /\gamma \mathrm{cos}(\theta _f+\theta _i)/\gamma }\frac{1}{\mathrm{cos}\pi /\gamma \mathrm{cos}(\theta _f\theta _i)/\gamma }\right],$$
(15)
where $`\gamma =\alpha /\pi `$ and $`\theta _f`$ (resp. $`\theta _i`$) is the final (resp. initial) scattering angle.

For rectangular billiards with Aharonov-Bohm flux lines the flux lines themselves are singular points and the exact solution for an infinite plane with a flux line carrying a flux $`\alpha `$  gives

$$D(\theta _f,\theta _i)=\frac{2\mathrm{sin}\pi \alpha }{\mathrm{cos}\frac{\theta _f\theta _i}{2}}e^{i(\theta _f\theta _i)/2}.$$
(16)
The main difference between pseudo-integrable models discussed in this paper and usual diffractive models is the divergence of diffraction coefficients (15) and (16) at certain directions (called optical boundaries because in the simplest case they separate illuminated regions from dark ones). Of course, exact solutions do not diverge even in vicinity of optical boundaries. The divergence comes from artificial separation of exact waves into geometrical and diffraction parts. Nevertheless, this formal divergence has profound effects on the structure of the trace formula. First, multiple diffraction along optical boundaries need a special treatment. Using a kind of uniform approximation in  it was demonstrated that for polygonal billiards such multiple diffraction produces terms proportional up to a numerical factor to $`l/k`$ where $`l`$ is the total length of the diffractive orbit. When $`l`$ is fixed and $`k\mathrm{}`$ (as in the usual approach to trace formulas) these terms are smaller than periodic orbit terms (10) but bigger than diffractive terms (14). But to compute spectral correlation functions one needs to consider a limit when $`k`$ is fixed and $`l\mathrm{}`$. In this limit multiple diffraction terms are bigger than both periodic orbit and diffraction terms. Another difficulty is related with the existence of terms corresponding to diffraction not exactly on optical boundaries but sufficiently close to them so their contributions are also large. Without exact summation of these quickly growing terms it is not possible to find spectral statistics of the systems considered.

It the next Section we argue that, nevertheless, these terms give negligible contribution to the value of the 2-point correlation form factor at the origin and only diagonal contributions of periodic orbits will be important for this quantity.

### 2.2 The 2-point correlation form factor

The 2-point correlation function is related with the level density by the formal expression

$$R_2(ϵ)=d(E+\frac{ϵ}{2})d(E\frac{ϵ}{2}),$$
(17)
where the brackets denote an energy averaging around $`E`$ on an energy window much larger than the mean level spacing $`1/\overline{d}`$, and much smaller than energy $`E`$.

The two-point correlation form factor is the Fourier transform of $`R_2(ϵ)`$ :

$$K_2(\tau )=_{\mathrm{}}^{\mathrm{}}\frac{dϵ}{\overline{d}}d(E+\frac{ϵ}{2})d(E\frac{ϵ}{2})e^{2i\pi \overline{d}ϵ\tau },$$
(18)
(the factors are chosen so that $`\tau `$ and $`K_2`$ are dimensionless).

Trace formulas, roughly speaking, state that the density of states can be represented as a sum over classical orbits (both periodic and diffractive)

$$d^{(osc)}(E)=\underset{p}{}C_pe^{iS_p(E)/\mathrm{}}+c.c.$$
(19)
Substituting this formal expansion into (17) and using the expansion

$$S(E+ϵ)S(E)+T(E)ϵ$$

where $`T(E)`$ is the period of classical motion one obtains

$`d(E+{\displaystyle \frac{ϵ}{2}})d(E{\displaystyle \frac{ϵ}{2}})=`$ (20)
$`{\displaystyle \underset{p_1,p_2}{}}C_{p_1}C_{p_2}^{}<\mathrm{exp}{\displaystyle \frac{i}{\mathrm{}}}(S_{p_1}(E)S_{p_2}(E))>e^{i(T_{p_1}+T_{p_2})ϵ/(2\mathrm{})}.`$
Here the terms corresponding to the sum of actions are omitted as it is assumed that they are washed out by the smoothing procedure.

The corresponding expression for the 2-point correlation form factor is the following:

$$K_2(\tau )=\underset{p_1,p_2}{}\frac{2\pi \mathrm{}}{\overline{d}}C_{p_1}C_{p_2}^{}<e^{i(S_{p_1}(E)S_{p_2}(E))/\mathrm{}}>\delta (\frac{T_{p_1}+T_{p_2}}{2}2\pi \mathrm{}\overline{d}\tau ).$$
(21)
The main difficulty in such an approach is the computation of the mean value of terms with action differences

$$F(E)=<e^{i(S_{p_1}(E)S_{p_2}(E))/\mathrm{}}>.$$
(22)
The best developed approximation (called the diagonal approximation) consists in taking into account only terms with exactly the same actions  i.e.

$$F(E)=\{\begin{array}{cc}1,\text{if}S_{p_1}(E)=S_{p_2}(E)& \\ 0,\text{if}S_{p_1}(E)S_{p_2}(E)& \end{array},$$
(23)
since terms with $`S_{p_1}(E)S_{p_2}(E)`$ will vanish by smoothing over $`E`$. In this approximation (assuming that for orbits with equal actions pre-factors are also equal (which is not always the case)) the 2-point correlation form factor takes the form

$$K_2^{(diag)}(\tau )=\underset{p}{}\frac{2\pi \mathrm{}}{\overline{d}}g_p^2|C_p|^2\delta (T_p2\pi \mathrm{}\overline{d}\tau ),$$
(24)
where $`g_p`$ is the multiplicity of a given periodic orbit (i.e. the number of orbits with exactly the same action) and the summation is performed over orbits with different actions. In particular for integrable and pseudo-integrable systems from Eq. (10) one gets

$$K_2^{(diag)}(\tau )=\frac{1}{8\pi ^2\overline{d}}\underset{p}{}\frac{|𝒜_p|^2}{ł_p}g_p^2\delta (l_p4\pi k\overline{d}\tau ),$$
(25)
where as before $`l_p`$ is the length of a periodic orbit and $`𝒜_p`$ is the surface occupied by a periodic orbit family.

It is instructive to perform the calculation for the simplest example of the rectangular billiard with sizes $`a`$ and $`b`$. A periodic orbit in this billiard is defined by 2 integers $`m`$, $`n`$ and its length is

$$l_p=\sqrt{(2ma)^2+(2nb)^2}.$$
(26)
As pairs $`(m,n)`$ and $`(m,n)`$ belong to the same family (or torus) the degeneracy is $`g_p=2`$ (we remind the reader that in the rectangular billiard the terms corresponding to $`m<0`$ are already taken into account in $`𝒜_p`$), and it is sufficient to compute the density of periodic orbits with positive $`m,n`$

$$\rho (l)=\underset{m,n0}{}\delta (ll_p).$$
(27)
Changing the summation over integers $`(m,n)`$ to the integration and using the substitution $`m=r\mathrm{cos}\varphi /(2a)`$ and $`n=r\mathrm{sin}\varphi /(2b)`$, one obtains by integrating over $`\varphi `$ from $`0`$ to $`\pi /2`$

$$\rho (l)=\frac{\pi l}{8𝒜}.$$
(28)
Since all families of periodic orbits in the rectangle cover the same area $`𝒜_p=2𝒜`$ and the length multiplicity is $`g_p=2`$, the 2-point correlation form factor for the rectangular billiard in the diagonal approximation is

$$K_2(\tau )=\frac{2𝒜^2}{\pi ^2\overline{d}}_0^{\mathrm{}}\frac{1}{l}\delta \left(l4\pi k\overline{d}\tau \right)\rho (l)𝑑l=1,$$
(29)
which is the expected value for the form factor of integrable systems . The diagonal approximation (23) is known (with physical accuracy) to be valid for generic integrable systems  and can be modified  to compute mean values of more than 2 actions in the exponent of (23).

For general systems the validity of the diagonal approximation is restricted only to small values of $`\tau `$ ,  and it is usually used to compute the first non-zero term of the expansion of the 2-point correlation form factor in powers of $`\tau `$.

For diffractive systems with finite diffraction coefficient one can use the diagonal approximation for both periodic orbit terms and diffractive terms. But when the diffraction coefficient diverges in certain directions these calculations lead to difficulties. For example, multiple diffraction on optical boundaries corresponding to $`n`$ repetitions of a primitive periodic orbit in pseudo-integrable billiard gives the following terms

$$d_{mult.diff.}(E)=\underset{l,n}{}\frac{l}{k}c_n\mathrm{cos}(knl),$$
(30)
where $`c_n`$ are certain numerical coefficients. The attempt to use the diagonal approximation for these terms leads to the following result

$$K_2^{(mult.diff.)}k^2\tau ^3,$$
(31)
if we take into account that the density of primitive periodic orbits in pseudo-integrable systems (at least for Veech systems (see the next Section)) differs only by a numerical factor from Eq. (28). But this expression contains powers of momentum $`k`$ and when $`k\mathrm{}`$ it cannot be correct. All terms corresponding to diffraction on or close to optical boundaries give similar quickly growing terms which cannot be treated separately. Without a ressumation of these terms the determination of spectral statistics of such models seems not possible. These arguments suggest the following scenario. The 2-point form factor is a sum of two terms

$$K_2(\tau )=f_1(k^\alpha \tau )+f_2(\tau ),$$
(32)
where $`\alpha `$ is a certain positive quantity. The first function, $`f_1(x)`$, describes a result of resummation of quickly growing terms connected with divergence of the diffraction coefficient and when $`x\mathrm{}`$ $`f_1(x)`$ should go quickly to zero. The second function, $`f_2(x)`$, is a contribution of diffraction far from optical boundaries and can be computed similarly to ordinary diffraction  in perturbation series of $`\tau `$. Of course, this is only a plausible conjecture and more detailed investigation should be done to give credit to it.

Though the divergence of the diffraction coefficient prevents the calculation of the 2-point correlation form factor in the full range, one can still use the trace formula (14) to find its behavior at the origin, $`\tau =0`$. The main point is that, even when the diffraction coefficient formally diverges, the exact waves remain finite and using a uniform approximation  one can demonstrate that the ratio

$$\frac{D(\stackrel{}{n},\stackrel{}{n}^{})}{\sqrt{kl}}$$
(33)
is bounded for all angles, lengths and momenta. Each term in the diffractive trace formula (14) is a product of certain number of these ratios and the total period of the corresponding composite orbit which appears due to the derivative over energy. Therefore it is of order of $`\tau `$ multiplied by a constant and in the limit $`\tau 0`$ all diffractive terms disappear. Only the periodic orbit contribution (10) remains important at small values of $`\tau `$. From Eq. (25) one concludes that for pseudo-integrable systems

$$K_2(0)=\underset{\tau 0}{lim}\frac{1}{8\pi ^2\overline{d}}\underset{p}{}\frac{|𝒜_p|^2}{l_p}g_p^2\delta (l_p4\pi k\overline{d}\tau ).$$
(34)
The main problem now is how to compute the density of periodic orbits and the distribution of the areas of periodic orbit families. For generic pseudo-integrable systems very little is known and no reliable calculations can be done. E.g. for general plane polygonal billiards with angles commensurable with $`\pi `$ it has only been proved $`\text{[27]},\text{[28]}`$ that the number of periodic orbits with length less than $`l`$, $`𝒩(l_pl)`$ obeys inequalities

$$c_1l^2<𝒩(l_pl)<c_2l^2$$
(35)
for certain constants $`c_1`$ and $`c_2`$ (depending on the polygon). But even the existence of an asymptotic law for $`𝒩(l_pl)`$ was not proved.

Fortunately, there is a sub-class of pseudo-integrable billiards for which all necessary quantities can be computed due to the existence of a hidden group structure, and the triangular billiard in the shape of the right triangle with one angle equal to $`\pi /n`$ belongs to this class. In the following Section we focus on those polygons.

## 3 Veech structures for polygonal billiards

We start the discussion of a hidden group structure of certain polygonal billiards with the simple example of the square billiard where the necessary ideas and methods can be illustrated clearly without technical difficulties.

### 3.1 A simple case: the square billiard

How can one evaluate the number of periodic orbits with length less than $`l`$ in a square billiard of size $`2a`$ with periodic boundary conditions? The exact expression for the length of the periodic orbits in the such billiard is, of course,

$$l_p=\sqrt{(2ma)^2+(2na)^2}$$
(36)
with $`m`$ and $`n`$ (see (11)). The number of periodic orbits with length less than $`l`$, $`𝒩(l_pl)`$, reads

$$𝒩(l_pl)=\underset{m,n}{}\mathrm{\Theta }\left(l2a\sqrt{m^2+n^2}\right)$$
(37)
and asymptotically when $`l\mathrm{}`$

$`𝒩(l_pl)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑m{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑n\mathrm{\Theta }\left(l2a\sqrt{m^2+n^2}\right)`$ (38)
$`=`$ $`{\displaystyle \frac{\pi l^2}{8a^2}}`$
if one sets $`m=(r\mathrm{cos}\phi )/2a`$ and $`n=(r\mathrm{sin}\phi )/2a`$. This is the number of all periodic orbits. More interesting questions and rich mathematical structure appear when one is interested in the calculation of the number of primitive periodic orbits $`𝒩_{pp}(l_pl)`$ (that is, orbits with $`m`$ and $`n`$ coprime).

The number of such orbits for a square billiard can easily be computed by using the inclusion-exclusion principle. The number of primitive periodic orbits with length less than $`l`$ is the total number of periodic orbits with length less than $`l`$ minus the number of orbits repeated $`p`$ times with prime $`p`$, to which we add orbits repeated $`p_1p_2`$ times, which had been subtracted twice, etc. Finally one concludes that

$`𝒩_{pp}(l_pl)`$ $`=`$ $`𝒩(l_pl){\displaystyle \underset{p}{}}𝒩(l_p{\displaystyle \frac{l}{p}})+{\displaystyle \underset{p_1,p_2}{}}𝒩(l_p{\displaystyle \frac{l}{p_1p_2}})`$ (39)
$``$ $`{\displaystyle \underset{p_1,p_2,p_3}{}}𝒩(l_p{\displaystyle \frac{l}{p_1p_2p_3}})\mathrm{}`$
Using the $`l^2`$ dependence of $`𝒩`$ in $`(\text{38})`$, we have

$`𝒩_{pp}(l_pl)`$ $`=`$ $`𝒩(l_pl)(1{\displaystyle \underset{p}{}}{\displaystyle \frac{1}{p^2}}+{\displaystyle \underset{p_1,p_2}{}}{\displaystyle \frac{1}{(p_1p_2)^2}}{\displaystyle \underset{p_1,p_2,p_3}{}}{\displaystyle \frac{1}{(p_1p_2p_3)^2}}\mathrm{})`$ (40)
$`=`$ $`𝒩{\displaystyle \underset{p}{}}(1{\displaystyle \frac{1}{p^2}})=𝒩{\displaystyle \frac{1}{\zeta (2)}}={\displaystyle \frac{6}{\pi ^2}}𝒩,`$
where

$$\zeta (s)=\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{n^s}=\underset{p}{}\frac{1}{1p^s}$$
(41)
is the Riemann zeta function.

From $`(\text{38})`$ one gets

$$𝒩_{pp}(l_pl)=\frac{3l^2}{4\pi a^2}.$$
(42)
Our aim is to generalize the previous calculation of $`𝒩_{pp}(l_pl)`$ to certain triangular billiards. This generalization naturally appears  when one considers carefully the usual geometrical picture of the free motion inside the square billiard. It is well known that any trajectory of such a motion can be unfolded to a straight line when instead of the square billiard one considers the motion on the covering space which for square billiard is a plane with infinite square lattice of the side $`2a`$. The vertices of this lattice (which are the images of the vertices of the initial square) have coordinates

$$x=2am,y=2an$$
(43)
with integers $`m`$ and $`n`$ and can be considered as the result of the application of a $`2\times 2`$ matrix with integer coefficient to a horizontal vector $`(2a,0)`$

$$\left(\begin{array}{cc}m& k\\ n& l\end{array}\right)\left(\begin{array}{c}2a\\ 0\end{array}\right)=\left(\begin{array}{c}2am\\ 2an\end{array}\right).$$
(44)
Thus, the periodic orbit lengths (36) are the distances between these vertices and the initial point $`(0,0)`$.

The problem of finding the number $`𝒩_{pp}`$ of primitive periodic orbits with length less than $`l`$ is therefore equivalent to the problem of finding out how many $`2\times 2`$ matrices with integer coefficients and determinant equal to $`1`$ (since $`m`$ and $`n`$ are coprime one can impose $`mlnk=1`$) exist with $`m^2+n^2x^2`$ for a given $`x`$ (or, which is equivalent, with $`n^2+l^2x^2`$). The group of $`2\times 2`$ matrices with integer coefficients and determinant equal to $`1`$ form a group $`SL(2,)`$ and in the next two Sections we shall discuss its main properties. Though this material is well known we find it useful to remind it informally.

### 3.2 The modular group

The subgroup of $`SL(2,)`$ containing all $`2\times 2`$ matrices with integer coefficients and determinant equal to $`1`$ is called the modular group $`SL(2,)`$. The standard representation of this group (see e.g. $`\text{[26]})`$ is the Poincaré half-plane $``$ with measure

$$ds^2=\frac{1}{y^2}(dx^2+dy^2);$$
(45)
A matrix $`gSL(2,)`$ is represented by the isometry

$$\begin{array}{cccc}\hfill g:& \hfill & & \hfill \\ & \hfill z& & \frac{mz+k}{nz+l}\hfill \end{array}$$
(46)
The modular group is generated by the translation $`T:zz+1`$ and the inversion $`S:z1/z`$, which correspond respectively to the matrices

$$\left(\begin{array}{cc}1& \alpha \\ 0& 1\end{array}\right)$$
(47)
(with $`\alpha =1`$ for the modular group) and

$$\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right).$$
(48)
Since the modular group is a discrete group, we can define its fundamental domain $`𝒟`$ (shown in fig. 1), that is the domain of the Poincaré half-plane $``$ that covers $``$ under the action of the representation $`(\text{46})`$ of the group.

In order to compute the number of matrices $`g=\left(\begin{array}{cc}m& k\\ n& l\end{array}\right)SL(2,)`$ verifying $`n^2+l^2x^2`$, we have to evaluate

$$𝒩(x)=\underset{\begin{array}{c}g\mathrm{\Gamma }_{\mathrm{}}\backslash G\\ n^2+l^2x^2\end{array}}{}1$$
(49)
where $`G=SL(2,)`$ and $`\mathrm{\Gamma }_{\mathrm{}}`$ is the subgroup of $`G`$ generated by the translations ($`\mathrm{\Gamma }_{\mathrm{}}=\{T^n,n\})`$: since the left multiplication by matrices of the form $`T^p`$

$$\left(\begin{array}{cc}1& p\\ 0& 1\end{array}\right)\left(\begin{array}{cc}m& k\\ n& l\end{array}\right)=\left(\begin{array}{cc}m+pn& k+pl\\ n& l\end{array}\right)$$
(50)
does not change $`n`$ and $`l`$ it is necessary to consider the quotient $`\mathrm{\Gamma }_{\mathrm{}}\backslash G`$ (i.e. 2 matrices which differ by $`T^p`$ are considered only once), so that the sum is convergent . If we assume that in the limit $`x\mathrm{}`$ the sums can be written as integrals over $`n`$ and $`l`$ with uniform measure (see later) in the form $`(B/\pi )dndl`$, we get

$$𝒩(x)=_{n^2+l^2x^2}\frac{B}{\pi }𝑑n𝑑l=\frac{B}{2}x^2.$$
(51)
### 3.3 Eisenstein series

In order to compute the coefficient $`B`$, let us introduce the Eisenstein series

$$E(z,s)=\underset{g\mathrm{\Gamma }_{\mathrm{}}\backslash G}{}\left(mg(z)\right)^s$$
(52)
for $`s>1`$. From expression $`(\text{46})`$ we get, since $`mlnk=1`$,

$$mg(z)=\frac{y}{|nz+l|^2}$$
(53)
where $`y=mz`$. Since $`mg^{}g(z)=mg(z)`$ for $`g^{}\mathrm{\Gamma }_{\mathrm{}}`$, the sum over $`\mathrm{\Gamma }_{\mathrm{}}\backslash G`$ is well defined. Let us first compute the asymptotic behavior of $`E(z,s)`$ when $`s1`$. For a given $`R`$, we can rewrite the sum $`(\text{52})`$ as a finite sum over elements of $`G`$ for which $`n^2+l^2<R^2`$ and a sum over elements for which $`n^2+l^2>R^2`$ which diverges as $`s1`$. The divergent part, $`n^2+l^2>R^2`$, is

$`E^{\text{div}}(z,s)`$ $`=`$ $`{\displaystyle \frac{y^s}{\pi }}{\displaystyle _{n^2+l^2>R^2}}{\displaystyle \frac{Bdndl}{|nz+l|^{2s}}}`$ (54)
$`=`$ $`{\displaystyle \frac{By^s}{\pi }}{\displaystyle _R^{\mathrm{}}}{\displaystyle _0^\pi }{\displaystyle \frac{r^{12s}drd\varphi }{\left[(x\mathrm{sin}\varphi +\mathrm{cos}\varphi )^2+(y\mathrm{sin}\varphi )^2\right]^s}}`$ (55)
Since

$$_R^{\mathrm{}}r^{12s}𝑑r=\frac{R^{2(1s)}}{2(s1)}\underset{s1}{}\frac{1}{2(s1)},$$
(56)
and the finite part of the Eisenstein series can be neglected as compared with the divergent part, we have

$$E(z,s)\underset{s1}{}\frac{yB}{2\pi (s1)}_0^\pi \frac{d\varphi }{(x\mathrm{sin}\varphi +\mathrm{cos}\varphi )^2+(y\mathrm{sin}\varphi )^2}.$$
(57)
The computation of the integral can be performed the following way : setting $`A=1x^2y^2`$, $`B=2x`$ and $`C=1+x^2+y^2`$, we get

$$_0^\pi \frac{d\varphi }{(x^2+y^2)\mathrm{sin}^2\varphi +\mathrm{cos}^2\varphi +2x\mathrm{sin}\varphi \mathrm{cos}\varphi }=2_0^\pi \frac{d\varphi }{A\mathrm{cos}2\varphi +B\mathrm{sin}2\varphi +C}$$
(58)
and since $`C^2A^2B^2=4y^2`$, the integral is equal to

$$_0^{2\pi }\frac{d\mathrm{\Psi }}{\sqrt{A^2+B^2}\mathrm{cos}\mathrm{\Psi }+C}=\frac{\pi }{y}$$
(59)
Finally

$$E(z,s)\underset{s1}{}\frac{B}{2(s1)}$$
(60)
and this limit does not depend on $`z`$.

Now let us integrate the series in $`(\text{52})`$ with the invariant measure $`d\mu (z)=dxdy/y^2`$ over a part, $`𝒟_Y`$, of the fundamental domain $`𝒟`$ corresponding to a restriction $`yY`$. If $`d`$ is the width of the fundamental domain,

$`{\displaystyle _{𝒟_Y}}E(z,s)𝑑\mu (z)`$ $`=`$ $`{\displaystyle \underset{g\mathrm{\Gamma }_{\mathrm{}}\backslash G}{}}{\displaystyle _{𝒟_Y}}\left(mg(z)\right)^s𝑑\mu (z)`$ (61)
$`=`$ $`{\displaystyle \underset{g\mathrm{\Gamma }_{\mathrm{}}\backslash G}{}}{\displaystyle _{g𝒟_Y}}\left(mz\right)^s𝑑\mu (z)`$
$``$ $`{\displaystyle _0^d}𝑑x{\displaystyle _0^Y}{\displaystyle \frac{dy}{y^2}}y^s`$ (62)
$`=`$ $`d{\displaystyle \frac{Y^{s1}}{s1}}.`$ (63)
In transformation from Eq. (61) to Eq. (62) we take into account that the image of fundamental region $`𝒟`$ (and $`𝒟_Y`$) under the action of $`G`$ is a certain region on the upper-half plane which under the action of $`\mathrm{\Gamma }_{\mathrm{}}`$ can be moved into a vertical strip of width $`d`$ (which is the fundamental region for the group $`\mathrm{\Gamma }_{\mathrm{}}`$). These images can not intersect and when $`Y\mathrm{}`$ will cover the whole strip with $`yY`$.

The asymptotic behavior of this integral is thus the following

$$\underset{s1}{lim}_{𝒟_Y}E(z,s)𝑑\mu (z)=\frac{d}{s1}.$$
(64)
By comparing this expression with Eq. (60) one concludes that the value of the constant $`B`$ is

$$B=2\frac{d}{\text{Vol }𝒟}$$
(65)
and from Eq. (51), the final answer for the density of primitive ($`n`$ and $`l`$ are coprime) periodic orbits in a square billiard is

$$𝒩_{pp}(n^2+l^2x^2)=\frac{d}{\text{Vol }𝒟}x^2$$
(66)
For the modular group $`d=1`$, $`\text{Vol }𝒟=\pi /3`$, and for a square billiard with side $`2a`$ we have $`x=l/(2a)`$ according to $`(\text{36})`$, so

$$𝒩_{pp}(l_p<l)=\frac{3l^2}{4\pi a^2}$$
(67)
which agree with Eq. (42) obtained by a different method.

### 3.4 Veech group for $`\pi /n`$ right triangles

#### 3.4.1 The symmetry group

The previous calculations were possible because we have found a group –the modular group– that relates periodic orbits in a square to a simple vector (see Eq. (44)). In order to generalize this construction for more complicated polygons it is important to point out that the modular group is the symmetry group of the unfolding of the square billiard (that is, the lattice whose unit cell is a $`1\times 1`$ square). Indeed this square lattice is, evidently, invariant under the following two transformations: the rotation by $`\pi /2`$ around the center of the square which we denote by $`S`$ and the translation of one coordinate (say $`x`$) by $`1`$ which we denote by $`T`$. In Cartesian coordinates these transformations are represented by the following matrices

$$S=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)T=\left(\begin{array}{cc}1& 1\\ 0& 1\end{array}\right).$$
(68)
The vertices of the lattice are also unchanged under the action of the group generated by $`S`$ and $`T`$. But it is well known that this group is exactly the modular group $`SL(2,)`$. Therefore the modular group plays a double role for a square billiard. First, it is the group of invariance of unfolding of the square and, second, it generates periodic orbits starting from a fix vector as in Eq. (44).

It has been proved by Veech  that for certain polygons (called the Veech polygons) there exists a group with similar properties. In particular, a $`\pi /n`$ right triangle (i.e. a triangle with angles $`\pi /2,\pi /n,\pi (n2)/2n`$) belongs to the Veech polygons , .

Let us consider this case in details. The geometrical construction of the unfolding of the classical trajectories in such a billiard is slightly different for $`n`$ even and odd. By reflections with respect to the sides corresponding to the $`\pi /n`$ angle the $`\pi /n`$ triangle can be unfolded to the regular $`n`$-gon. For $`n`$ even the opposite sides of this $`n`$-gon should be identified by translations (see Fig. 2a ). For $`n`$ odd one has to consider 2 regular $`n`$-gons reflected with respect to one side and to identify parallel sides by translations as in Fig. 2b. The resulting surface is the surface of genus $`(n1)/2`$ for $`n`$ odd (see 2) to which all trajectories belong.

From this construction it is clear that if a group of invariance exists it should include the rotation by $`2\pi /n`$ around the center of these $`n`$-gons. In Cartesian coordinates this rotation is defined by the following matrix

$$\sigma _n=\left(\begin{array}{cc}\mathrm{cos}\frac{2\pi }{n}& \mathrm{sin}\frac{2\pi }{n}\\ \mathrm{sin}\frac{2\pi }{n}& \mathrm{cos}\frac{2\pi }{n}\end{array}\right).$$
(69)
To find other transformations which leave this surface invariant it is necessary to consider a few families of periodic orbits.

For $`n`$ even, we define in the $`n`$gon two important elementary families of periodic orbits: the first one is the family of horizontal primitive periodic orbits, the second one is the family of primitive periodic orbits making an angle $`\pi /n`$ with the horizontal (see fig. 3). For $`n`$ odd we only define the first family (see fig. 4).

For $`n`$ even ($`n=4p`$ or $`n=4p+2`$), the first family has orbits with lengths

$$L_j=4\mathrm{cos}\frac{(2j1)\pi }{n}\mathrm{cos}\frac{\pi }{n}$$
(70)
and widths

$$W_j=2\mathrm{cos}\frac{(2j1)\pi }{n}\mathrm{sin}\frac{\pi }{n}$$
(71)
with $`1jp`$. The second family has orbits with lengths

$$l_j^{}=4\mathrm{cos}\frac{(2j2)\pi }{n}\mathrm{cos}\frac{\pi }{n}$$
(72)
and widths

$$w_j^{}=2\mathrm{cos}\frac{(2j2)\pi }{n}\mathrm{sin}\frac{\pi }{n}$$
(73)
with $`2jp`$ if $`n=4p`$ or $`2jp+1`$ if $`n=4p+2`$. The orbit with $`j=1`$ is special: it has a length and a width equal to

$$l_1^{}=2\mathrm{cos}\pi /n,w_1^{}=2\mathrm{sin}\pi /n$$
(74)
For $`n`$ odd $`(n=2p+1)`$, the lengths and widths are the following

$$L_j=4\mathrm{sin}\frac{2j\pi }{n}\mathrm{cos}\frac{\pi }{n}$$
(75)
and

$$W_j=2\mathrm{sin}\frac{2j\pi }{n}\mathrm{sin}\frac{\pi }{n}.$$
(76)
with $`1jp`$.

It is of interest to compute the ratio of the length of each periodic orbit family to its width. From the above formulas it follows that for all families, except the one with $`j=1`$ for even $`n`$, this ratio is the same

$$\frac{l}{w}=2\mathrm{cot}\frac{\pi }{n}.$$
(77)
For the exceptional family (74) this ratio is 2 times smaller.

The unfolding of any family of periodic orbits gives an infinite strip of points of period $`L_p`$ and width $`W_p`$. If there is a group of invariance of the unfolded surface, it should include a transformation which leaves invariant periodic orbit strips. Assume that the strip is oriented horizontally. In this case one sees that the shift of the form

$$\left(\begin{array}{cc}1& \alpha \\ 0& 1\end{array}\right)$$
(78)
leaves points of the strip invariant provided that

$$L_p=n\alpha W_p,$$
(79)
where $`n`$ is an integer.

Because for periodic orbits considered the ratio (77) is constant the invariance group should include the following transformation

$$\tau _n=\left(\begin{array}{cc}1& 2\mathrm{cot}\frac{\pi }{n}\\ 0& 1\end{array}\right).$$
(80)
Veech proved ,  that the invariance group for $`\pi /n`$ triangle is a discrete subgroup $`\mathrm{\Gamma }_n`$ of $`SL(2,)`$ generated by the two elements (69) and (80). Similarly to the relation (44) for the square, periodic orbits in this triangle are generated by the action of $`\mathrm{\Gamma }_n`$ over the elementary families of periodic orbits considered above.

We shall call the members of these families “basis orbits” (and label them by the index $`i`$). We define the corresponding basis vectors $`v_i`$ by $`v_j=(L_j,0)`$ and $`v_j^{}=(l_j^{}\mathrm{cos}(\pi /n),l_j^{}\mathrm{sin}(\pi /n))`$, so that $`\{v_i\}=\{v_j,v_j^{}\}`$ for $`n`$ even and $`\{v_i\}=\{v_j\}`$ for $`n`$ odd.

An element $`g`$ of the symmetry group $`\mathrm{\Gamma }_n`$ has the following matrix representation

$$g=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right).$$
(81)
The result of the action of this element to one of the basis vectors $`v_i=(v_{i1},v_{i2})`$ gives the coordinates of a new primitive periodic orbit (more precisely, a periodic orbit situated on the boundary of periodic orbit pencil)

$$gv_i=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)\left(\begin{array}{c}v_{i1}\\ v_{i2}\end{array}\right),$$
(82)
and the length of this primitive periodic orbit is the length of this vector.

The first family of basis vectors for even $`n`$ and for all basis vectors for odd $`n`$ can be chosen in the form $`v_i=(L_i,0)`$ and the lengths of periodic orbits generated by applying the group $`\mathrm{\Gamma }_n`$ is

$$L_g=\sqrt{a^2+c^2}L_i,$$
(83)
where $`a`$, $`c`$ are matrix elements of $`g`$ (82). The second basis periodic orbits, $`v_i^{}`$, are obtained from horizontal vectors by rotation by $`\pi /n`$. But the matrix corresponding to the inverse of this rotation

$$r=\left(\begin{array}{cc}\mathrm{cos}\frac{\pi }{n}& \mathrm{sin}\frac{\pi }{n}\\ \mathrm{sin}\frac{\pi }{n}& \mathrm{cos}\frac{\pi }{n}\end{array}\right)$$
(84)
does not belong to our group $`\mathrm{\Gamma }_n`$. Nevertheless, this matrix plays the role of a Hecke operator, namely, even if it does not belong to $`\mathrm{\Gamma }_n`$ the conjugation of any matrix from this group does belong to $`\mathrm{\Gamma }_n`$: if $`g\mathrm{\Gamma }_n`$, then $`r^1gr\mathrm{\Gamma }_n`$. To prove it let us note that

$$r^1\sigma _n^pr=\sigma _n^p$$
(85)
where $`\sigma _n`$ is the generator (69) because all rotations commute, and it is easy to check that

$$r^1\tau _nr=\sigma _n\tau _n^1,r^1\tau _n^1=\tau _n\sigma _n^1.$$
(86)
The right-hand sides of these relations belong to $`\mathrm{\Gamma }_n`$ and as all matrices from $`\mathrm{\Gamma }_n`$ can be written as a product of generators we get $`r^1gr\mathrm{\Gamma }_n`$ for $`g\mathrm{\Gamma }_n`$.

Using this conjugation one can rotate the second family of periodic orbits for even $`n`$ by $`\pi /n`$ and the lengths of the orbits generated by the vectors $`v_i^{}`$ will be related to matrix elements of $`\mathrm{\Gamma }_n`$ by the same relation as in (83)

$$L_g^{}=\sqrt{a^2+c^2}l_i^{}.$$
(87)
Therefore, to find the distribution of periodic orbit lengths it is necessary to compute the distribution of $`a^2+b^2`$ for matrices from $`\mathrm{\Gamma }_n`$, which has been done for the modular group in the previous section. Eq. (66) can be derived the same way for $`\mathrm{\Gamma }_n`$. According to the previous section one can compute the density of periodic orbits and other quantities as well by investigation of the fundamental domain of $`\mathrm{\Gamma }_n`$.

The distribution of areas of periodic orbit families is also easy to obtain: as all matrices from $`\mathrm{\Gamma }_n`$ have unit determinant, the area covered by the pencil corresponding to $`gv_i`$ ($`g\mathrm{\Gamma }_n`$) is the same as the area covered by the pencil corresponding to basis vectors $`v_i`$, i.e. it is equal to $`L_iW_i`$. In other words, there is a one to one correspondence between pencils of primitive periodic orbits and vectors $`gv_i`$ for $`g\mathrm{\Gamma }_n`$. The discrete group $`\mathrm{\Gamma }_n`$ is related to periodic orbits in the $`\pi /n`$ right triangle in the same way as the modular group is related to periodic orbits in the square.

#### 3.4.2 The density of periodic orbits

The fundamental domains of the symmetry groups $`\mathrm{\Gamma }_n`$ for $`n`$ even and odd are described in Figs. 5 and 6 respectively. For $`n`$ even, it is the union of two triangles with angles $`2\pi /n,0,0`$ on the Poincaré half-plane : the area of the domain is $`\text{Vol }𝒟=2\pi (n2)/n`$, and its width is $`2\mathrm{cot}(\pi /n)`$.

For $`n`$ odd the two triangles have angles $`\pi /n`$, $`\pi /2`$ and $`0`$, therefore $`\text{Vol }𝒟=\pi (n2)/n`$; the width of the fundamental domain is $`2\mathrm{cot}\pi /n`$. These shapes of fundamental domains can be obtained by taking into account that the group $`\mathrm{\Gamma }_n`$ considered as a group acting on the Poincaré upper-half plane as in Section 3.2 includes (i) the translation by $`2\mathrm{cot}\pi /n`$ and (ii) the rotation around point $`i`$ by $`2\pi /n`$ for $`n`$ odd and by $`4\pi /n`$ for $`n`$ even. This difference between even and odd $`n`$ is related to the fact that the rotation by angle $`\pi `$ corresponds to the transformation $`gg`$, but these 2 matrices are represented by the same function on $``$ (see $`(\text{46})`$). For odd $`n`$ ($`n=2q+1`$), the group generated by the generator (69) contains rotations by angles $`2\pi j/(2q+1)`$ for $`j=0,1,\mathrm{}2q`$. The value $`j=q+1`$ corresponds to the rotation by $`\pi +\pi /n`$. As the rotation by $`\pi `$ is the identity transformation, the rotation by $`\pi /n`$ belongs to $`\mathrm{\Gamma }_n`$ and it is a primitive generator of the subgroup $`\{\sigma _n^p,p\}`$ of $`\mathrm{\Gamma }_n`$.

For even $`n`$ ($`n=2q`$), the rotation by $`2\pi q/n`$ is the identity, therefore the rotation by $`\pi /n`$ does not belong to the group and the primitive generator of the subgroup $`\{\sigma _n^p,p\}`$ is the rotation by $`2\pi /n`$.

Due to $`(\text{66})`$, we now get the number of matrices

$$g=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)\mathrm{\Gamma }_n$$
(88)
verifying $`a^2+c^2x^2`$ :

$$𝒩(a^2+c^2x^2)=\{\begin{array}{c}\frac{n}{\pi (n2)}\mathrm{cot}\frac{\pi }{n}x^2\text{(}n\text{ even)}\\ \\ \frac{2n}{\pi (n2)}\mathrm{cot}\frac{\pi }{n}x^2\text{(}n\text{ odd)}\end{array}$$
(89)
These formulas give the total number of matrices verifying $`a^2+c^2x^2`$. But due to the existence of rotation matrices in the group $`\mathrm{\Gamma }_n`$ each primitive periodic orbit length appears a few times in the above calculations. This multiplicity corresponds to different unfoldings of a given periodic orbit.

For $`n`$ odd the $`2n`$ matrices of the form $`\{\pm \beta _n^kg,0k2n1\}`$, where $`\beta _n`$ is the matrix of rotation by $`\pi /n`$, give rise to one primitive periodic orbit. For $`n`$ even there exist $`n`$ matrices of the form $`\{\pm \sigma _n^kg,0kn1\}`$ (where $`\sigma _n`$ is the matrix of rotation by $`2\pi /n`$ (69)) which describe one periodic orbit.

For $`g`$ given by $`(\text{88})`$, the length of $`gv_i`$ is $`L_i\sqrt{a^2+c^2}`$. So the number of primitive periodic orbits of type $`gv_i`$ less than $`l`$ is

$$𝒩_{i,pp}(L_p<l)=\frac{1}{\pi (n2)}\mathrm{cot}\frac{\pi }{n}\left(\frac{l}{L_i}\right)^2.$$
(90)
The number of all primitive periodic orbits is the sum over all such contributions:

$$𝒩_{pp}(L_p<l)=𝒞\frac{l^2}{𝒜},$$
(91)
where

$$𝒜=\frac{1}{4}\mathrm{sin}\frac{2\pi }{n}$$
(92)
is the area of our triangle and

$$𝒞=\frac{1}{2\pi (n2)}\mathrm{cos}^2\frac{\pi }{n}\underset{i}{}\frac{1}{L_i^2}.$$
(93)
For $`n`$ odd ($`n=2p+1`$) Eq. (75) gives

$$\underset{k=1}{\overset{p}{}}\frac{1}{L_k^2}=\frac{1}{16\mathrm{cos}^2\pi /n}\underset{k=1}{\overset{p}{}}\frac{1}{\mathrm{sin}^2(2\pi k/n)}.$$
(94)
For $`n`$ even ($`n=4p+2ϵ`$, $`epsilon=0,1`$) from Eqs. (70), (72), and (74) it follows

$`{\displaystyle \underset{i=1}{\overset{p}{}}}{\displaystyle \frac{1}{L_i^2}}={\displaystyle \frac{1}{16\mathrm{cos}^2\pi /n}}\left({\displaystyle \underset{j=1}{\overset{p}{}}}{\displaystyle \frac{1}{\mathrm{cos}^2(2j1)\pi /n}}+{\displaystyle \underset{j=2}{\overset{p+ϵ}{}}}{\displaystyle \frac{1}{\mathrm{cos}^2(2j2)\pi /n}}+4\right)`$
$`={\displaystyle \frac{1}{16\mathrm{cos}^2\pi /n}}\left({\displaystyle \underset{k=1}{\overset{(n2)/2}{}}}{\displaystyle \frac{1}{\mathrm{cos}^2k\pi /n}}+4\right).`$ (95)
The last sums can be calculated using the evident formulas (for another method of calculation see )

$$\frac{1}{\mathrm{sin}^2x}=\underset{q=\mathrm{}}{\overset{\mathrm{}}{}}\frac{1}{(xq\pi )^2},$$
(96)
and

$$\frac{1}{\mathrm{cos}^2\pi x/2}=\frac{4}{\pi ^2}\underset{q=1}{\overset{\mathrm{}}{}}(\frac{1}{(2q1x)^2}+\frac{1}{(2q1+x)^2}).$$
(97)
Taking into account that $`\mathrm{sin}^2(k\pi /n)=\mathrm{sin}^2((nk)\pi /n)`$ for odd $`n`$ and performing the following transformations

$$\underset{k=1}{\overset{n1}{}}\underset{q=\mathrm{}}{\overset{\mathrm{}}{}}\frac{1}{(kqn)^2}=\underset{t=\mathrm{}}{\overset{\mathrm{}}{}}\frac{1}{t^2}\underset{q=\mathrm{}}{\overset{\mathrm{}}{}}\frac{1}{n^2q^2}=\frac{\pi ^2}{3}(1\frac{1}{n^2}),$$
(98)
one obtains that for odd $`n`$

$$\underset{k=1}{\overset{(n1)/2}{}}\frac{1}{\mathrm{sin}^2(2\pi k/n)}=\frac{n^21}{6}.$$
(99)
Similarly for even $`n`$

$$\underset{k=1}{\overset{(n2)/2}{}}\frac{1}{\mathrm{cos}^2k\pi /n}=\frac{n^24}{6}.$$
(100)
Therefore the value of constant $`𝒞`$ is

$$𝒞=\frac{1}{192\pi (n2)}\{\begin{array}{cc}(n^21)& \text{for odd}n\\ (n^2+20)& \text{for even}n\end{array}.$$
(101)
This result corresponds to primitive periodic orbits with geometrically different lengths: time-reversal orbits are not included in the summation. They give an additional factor of 2 in Eq. (101). In  the orbits corresponding to different unfoldings of the same periodic orbit have been included in the asymptotic formula, which leads asyptotically to the additional factor $`n`$ for even $`n`$ and $`2n`$ for odd $`n`$ in Eq. (101). Furthermore, if one needs all periodic orbits including repetitions Eq. (101) should be multiplied by $`\pi ^2/6`$ as in Section 3.1.

In Fig. 7 we present numerical results of the cumulative density of primitive periodic orbits in $`\pi /8`$ right triangular billiard (with area $`𝒜=4\pi `$) when all orbits (time-reversal and for different unfoldings) are included. The solid line is the best quadratic fit to these data

$$𝒩_{pp}(L_p<l)=.0294l^2.6055l+.3617.$$
(102)
One sees that this fit can hardy be distinguished from numerical results. The theoretical prediction for the coefficient in front of $`l^2`$ is $`𝒞/𝒜`$ according to Eq. (91); $`𝒞`$ is given by Eq. (101), and has to be multiplied by $`n`$ since the numerical computation has taken into account the repetitions of each orbit, and by $`2`$ as time-reversal orbits are taken into account as well:

$$\frac{𝒞}{𝒜}=\frac{7}{24\pi ^2}.0295,$$
(103)
which is in excellent agreement with numerical calculations.

### 3.5 Explicit calculation of the 2-point form factor for the $`\pi /n`$ right triangle

#### 3.5.1 First case: no degeneracy of the lengths

We assume in this section that there is no degeneracy between the lengths of the periodic orbits (except the ones connected by the time-reversal transformation) or more carefully, that there is no pair of primitive periodic orbits whose lengths are commensurable. Since the lengths of the $`gv_i`$ are proportional to the lengths of the $`v_i`$, the necessary requirement for the validity of this condition is the absence of commensurability relations between the $`L_i`$.

In this case the 2-point correlation form factor in the diagonal approximation is done by Eq. (25). The sum over different periodic orbits can be split into a sum over all types of periodic orbits, then the sum over periodic orbits of each type can be replaced by an integral; since the density $`\rho _i`$ of periodic orbits of type $`i`$ only takes into account once the periodic orbit and its time-reverse, the degeneracy is $`g_p=2`$ and

$$\underset{p}{}g_p^2=4\underset{i}{}_0^{\mathrm{}}𝑑l\rho _i(l).$$
(104)
In $`(\text{25})`$, $`𝒜_p`$ is the area occupied by a pencil of periodic orbits of length $`L_p`$; but this area is the same for all trajectories belonging to the same family $`i`$. So we just have to evaluate the area $`𝒜_i=L_iW_i`$ occupied by an elementary orbit of type $`i`$. The lengths $`(\text{70})`$, $`(\text{72})`$ and $`(\text{75})`$ are

$$L_k=4\mathrm{cos}\frac{k\pi }{n}\mathrm{cos}\frac{\pi }{n},1kp1$$
(105)
and $`L_0=2\mathrm{cos}\pi /n`$ if $`n=2p`$, and

$$L_k=4\mathrm{sin}\frac{2k\pi }{n}\mathrm{cos}\frac{\pi }{n},1kp$$
(106)
if $`n=2p+1`$. The widths are only half the widths $`W_i`$ given by $`(\text{71})`$, $`(\text{73})`$ and $`(\text{76})`$ since each fundamental pencil is symmetric with respect to the line joining two images of the $`\pi /2`$ corner of the triangle (see fig. 3 and 4). So the $`W_k`$ are

$$W_k=\mathrm{cos}\frac{k\pi }{n}\mathrm{sin}\frac{\pi }{n},0kp1$$
(107)
if $`n=2p`$, and

$$W_k=\mathrm{sin}\frac{2k\pi }{n}\mathrm{sin}\frac{\pi }{n},1kp$$
(108)
if $`n=2p+1`$. Then for small $`\tau `$

$$K_2(\tau )=4\underset{i}{}\frac{𝒜_i^2}{8\pi ^2}_0^{\mathrm{}}\frac{1}{l\overline{d}}\delta \left(l4\pi k\overline{d}\tau \right)\rho _i(l)𝑑l,$$
(109)
and replacing the density of orbits of type $`i`$ by its mean value

$$\rho _i=\frac{\pi ^2}{6}\frac{d𝒩_{i,pp}(L_p<l)}{dl}$$
(110)
where $`𝒩_{i,pp}(L_p<l)`$ is given by Eq. (90) we obtain, when performing the integral,

$$K_2(\tau )=\frac{\mathrm{cot}\pi /n}{6\pi (n2)\overline{d}}\underset{k}{}W_k^2;$$
(111)
where the average density of states is $`\overline{d}=𝒜/4\pi `$ with

$$𝒜=\frac{1}{4}\mathrm{sin}\frac{2\pi }{n}$$
(112)
is the area of the triangle.

The sum $`(\text{109})`$ over the widths $`(\text{107})`$ and $`(\text{108})`$ gives

$$\underset{k}{}W_k^2=\{\begin{array}{c}\frac{n+2}{4}\mathrm{sin}^2\frac{\pi }{n},n\text{ even}\\ \\ \frac{n}{4}\mathrm{sin}^2\frac{\pi }{n},n\text{ odd}\end{array}.$$
(113)
So we finally get

$$K_2(\tau )=\frac{n+ϵ(n)}{3(n2)}$$
(114)
with

$`ϵ(n)=0\text{when }n\text{ is odd},`$ (115)
$`ϵ(n)=2\text{when }n\text{ is even}.`$ (116)
#### 3.5.2 Second case: degeneracy of the lengths

We have assumed in the previous Section that the lengths of all primitive periodic orbits were non-commensurable. In the case of the $`\pi /n`$ right triangle, there may exist a commensurability relation between the $`L_k`$ given by $`(\text{105})`$ or $`(\text{106})`$ if there is one between the $`\mathrm{cos}(k\pi /n)`$ ($`0kp1`$) for $`n`$ even, or between the $`\mathrm{sin}(2k\pi /n)`$ ($`1kp`$) for $`n`$ odd. It is shown in  that if $`n`$ is an odd prime, there is no such relation between the $`\mathrm{sin}(2k\pi /n)`$. Ref. $`(\text{[31]})`$ deals with the case $`(k,n)=1`$ and gives the same conclusion. It seems that in the general case, only one relation of that kind exists between the $`\mathrm{cos}(k\pi /n)`$, which is

$$2\mathrm{cos}(\frac{n}{3}\frac{\pi }{n})=\mathrm{cos}(0)$$
(117)
and that no relation exists between terms with sinus. Therefore the only degeneracy occurs in the case where $`n`$ is even and $`3|n`$, that is $`n6`$. In that case, from $`(\text{10})`$ we get

$`d(E+{\displaystyle \frac{ϵ}{2}})d(E{\displaystyle \frac{ϵ}{2}})=K^{diag}+`$
$`4{\displaystyle \underset{p^+p^{}_{}{}^{}+}{}}{\displaystyle \frac{𝒜_p𝒜_p^{}}{16\pi ^2}}{\displaystyle \frac{1}{2\pi k\sqrt{l_pl_p^{}}}}e^{ik(l_pl_p^{})+i\frac{ϵ}{4k}(l_p+l_p^{})}+\text{c.c.,}`$ (118)
where $`K^{diag}`$ is the usual diagonal approximation ($`\text{2.2})`$. $`p^+`$ means that we only count for one orbit in the sum the orbit and its time-reverse, therefore there is a coefficient $`4`$.

If there is a relation $`m_1L_1=m_2L_2`$ (with $`m_1`$ and $`m_2`$ coprime) between two lengths of primitive periodic orbits, we have a contribution $`R_2^{\text{deg}}`$ to the $`2`$-point correlation function $`(\text{17})`$ which comes from orbits of lengths $`qm_1L_1`$ and $`qm_2L_2`$, $`q^{}`$ :

$`R_2^{\text{deg}}`$ $`=`$ $`{\displaystyle \frac{2𝒜_1𝒜_2}{4\pi ^2.2\pi k}}{\displaystyle \underset{q}{}}{\displaystyle \frac{1}{\sqrt{qm_1L_1qm_2L_2}}}e^{i\frac{ϵ}{4k}(qm_1L_1+qm_2L_2)}`$ (119)
$`=`$ $`{\displaystyle \frac{2𝒜_1𝒜_2}{4\pi ^2.2\pi k}}{\displaystyle \underset{q}{}}{\displaystyle \frac{1}{qm_1L_1}}e^{i\frac{ϵ}{2k}(qm_1L_1)}.`$
The sum over all repetition numbers $`q`$ of a function of $`qL_1`$ (where $`L_1`$ is a primitive periodic orbit) can be replaced by an integral :

$$R_2^{\text{deg}}=\frac{𝒜_1𝒜_2}{4\pi ^3k}_0^{\mathrm{}}𝑑l\frac{1}{m_1l}e^{i\frac{ϵ}{2k}m_1l}\rho _1(l)$$
(120)
where the density $`\rho _1(l)`$ is the density of periodic orbits of type $`1`$ (that is, with length proportional to $`L_1`$) with length less than $`l`$, given by $`(\text{110})`$. Performing the Fourier transform $`(\text{18})`$ and the integral over $`l`$ gives

$$K_2(0)=\frac{\mathrm{cot}\pi /n}{6\pi (n2)\overline{d}}\left[\underset{i}{}W_i^2+2\frac{W_1W_2}{m_1m_2}\right].$$
(121)
In our case the degeneracy is given by $`(\text{117})`$ and

$$W_0W_{{\displaystyle \frac{n}{3}}}=\frac{1}{2}\mathrm{sin}^2\frac{\pi }{n}.$$
(122)
We finally obtain

$$K_2(0)=\frac{n+ϵ(n)}{3(n2)}$$
(123)
with

$$ϵ(n)=\{\begin{array}{cc}0\hfill & \text{when }n\text{ is odd}\hfill \\ 2\hfill & \text{when }n\text{ is even and }3n\hfill \\ 6\hfill & \text{when }n\text{ is even and }3n.\hfill \end{array}.$$
(124)
This formula is the main result of our calculations for the triangular billiards. It clearly demonstrated the peculiarities of spectral statistics for pseudo-integrable systems. The non-zero value of the form factor ($`<1`$) at the origin does not correspond to any random matrix ensemble but it is typical for intermediate statistics , .

### 3.6 Comparison with numerical calculations

To compare the prediction (123) with numerical results we have computed 20000 levels for triangular billiards in the shape of a right triangle with one angle $`\pi /n`$ for all $`n=5,7\mathrm{},30`$ (the case of $`n=6`$ is integrable). For each triangle we take levels from 15000 till 20000 and compute numerically the corresponding 2-point correlation form factor. A typical result is presented in Fig. 8.

From data like this, it is quite difficult to find the value of the form factor at the origin because $`\tau 0`$ corresponds, according to Eq. (18), to an infinitely large energy difference in the $`2`$-point correlation function: therefore numerically we always have $`K_2(0)=0`$. We found it convenient first to fit the numerical data to the following simple expression for the form factor,

$$K_2(\tau )=\frac{a^22a+4\pi ^2\tau ^2}{a^2+4\pi ^2\tau ^2}.$$
(125)
and then from it compute $`K_2(0)`$.

$$K_2(0)=1\frac{2}{a}.$$
(126)
The form (125) has been chosen because (i) one wants a simple expression, (ii) when $`\tau \mathrm{}`$ the form factor should go to 1, (iii) to describe the level repulsion it is necessary that

$$_0^{\mathrm{}}\left(1K_2(\tau )\right)𝑑\tau =\frac{1}{2},$$
(127)
and (iv) the expression (125) when $`a=4`$ equals the form factor of the so-called semi-Poisson model ,  which serves as a reference point for intermediate statistics.

We stress that the above expression has no solid theoretical explanations and it is used because it relatively well describes our numerical results. The only fitting parameter is $`K_2(0)`$ related with $`a`$ by Eq. (126). We tried two fitting procedures. First we fit Eq. (125) for the data with all $`\tau `$ or, second, to decrease the influence of very small $`\tau `$, where numerical accuracy is not very good, we did not consider the data with $`0<\tau <0.25`$. In Fig. 8 these two fits are presented. The first one gives $`K_2(0)0.44`$ and the second one $`K_2(0)0.565`$. The expected value $`(\text{123})`$ for $`n=8`$ is $`5/90.56`$.

In Fig. 9 the results of such fitting procedures are given for all triangles. Lower two curves correspond to these fits and the upper curve is the predictions $`(\text{123})`$. (Of course, only points are important. Curves are presented for clarity.)

The numerical results quite well follow theoretical formula $`(\text{123})`$ but there is a small shift which decreases when the region of small $`\tau `$ is ignored. This difference between the curves seems to be a consequence of the fact that the result $`(\text{123})`$ corresponds to asymptotic limit $`k\mathrm{}`$ but numerical calculations have been performed at large but finite energy.

To check this point we present in Fig. 10 the results of the calculation of the mean number variance, $`\mathrm{\Sigma }^{(2)}(L)`$, for the $`\pi /30`$ right triangle, in 10 energy intervals $`[8000k,8000(k+1)]`$, $`0k9`$ (the energy increases from bottom to top). It is well known (see e.g. ) that the behavior of the mean number variance at large distances is related with the value of the form factor at the origin by simple formula

$$\mathrm{\Sigma }^{(2)}(L)K_2(0)L\text{when }L\mathrm{}.$$
(128)
From Fig. 10 it is clearly seen that even for 80000 levels the curve does not stabilize. To find its limiting behavior we extrapolate point by point (with $`L`$ fixed) this ten curves with a fit $`A(L)+B(L)/\sqrt{k}`$. (It means that for each $`L`$ we fit 10 points to find the best $`A(L)`$ and $`B(L)`$.) The limit curve (i.e. $`A(L)`$) is the most upper curve in Fig. 10. It perfectly reproduces the expected features of $`\mathrm{\Sigma }^{(2)}(L)`$: it is a straight line with slope $`K_2(0)=0.38`$ corresponding to the expected value $`(\text{123})`$ for $`n=30`$. In the same way, Prosen and Casati  have computed $`\mathrm{\Sigma }^{(2)}(L)`$ for triangle billiards with angle $`\pi /5`$ for much larger values of the energy, and it seems that such fit works well for their calculations and the result for $`K_2(0)`$ agrees with (123). These (and other) calculations clearly demonstrate that the value of the 2-point correlation form factor at the origin converges slowly to the theoretical result (123) with increasing energy. This behavior may be a consequence of the conjectured existence of two different terms $`(\text{32})`$ in the form factor and, in the final extent, a manifestation of the strong diffraction in vicinity of optical boundaries.

## 4 Rectangular billiard with a flux line

### 4.1 Preliminary calculations

This section is devoted to the study of a rectangular billiard with the Aharonov-Bohm flux line  at a point $`\stackrel{}{r_0}=(x_0,y_0)`$ inside the rectangle. In the polar coordinates, $`r`$, $`\phi `$, around this point the vector potential of the flux line has only $`\phi `$ component

$$A_\phi =\frac{\alpha }{r}$$
(129)
and the 2-dimensional Schrödinger equation for the motion in this potential is (when $`\mathrm{}=c=1`$ and $`m=1/2`$)

$$\left[\frac{^2}{r^2}+\frac{1}{r}\frac{}{r}+\frac{1}{r^2}\left(\frac{}{\phi }i\alpha \right)^2+E_n\right]\mathrm{\Psi }_n(r,\phi )=0.$$
(130)
Similarly to triangular billiards discussed in previous Sections this model belongs to the class of diffractive systems. The diffraction coefficient for the scattering on the flux line (16) diverges in the forward direction but as for pseudo-integrable systems the contribution of diffractive orbits can be neglected when computing the value of the 2-point correlation form factor at the origin.

It is well known that the Aharonov-Bohm potential (129) does not change classical trajectories but gives an additional phase, $`\mathrm{\Delta }\varphi `$, when a trajectory turns $`n`$ times around the flux line

$$\mathrm{\Delta }\varphi =2n\pi \alpha .$$
(131)
Therefore the contribution of periodic orbit to the trace formula (10) will contain an additional phase depending on the winding number of the trajectory around the flux line.

Periodic orbits in the rectangle of sides $`a`$, $`b`$ are determined by two integers $`M`$ and $`N`$ in the usual way and they are characterized by their length

$$l_p=\sqrt{(2Ma)^2+(2Nb)^2},$$
(132)
the area occupied by the periodic orbit family, and the winding number around the flux line. Each pencil of primitive periodic orbits occupies an area $`2ab=2𝒜`$, so its width is $`2𝒜/l_p`$. The images of the flux line in the unfolding of the rectangular billiard are located at the points

$$((\zeta _1+2k)a,(\zeta _2+2k^{})b).$$
(133)
Here $`\zeta _i`$ takes the values $`ϵ_i`$ or $`2ϵ_i`$ ($`i=1`$, $`2`$), where

$$ϵ_1=\frac{x_0}{a}ϵ_2=\frac{y_0}{b}$$
(134)
are the ratios of coordinates of the flux line to the corresponding sides, and $`k,k^{}`$ (see Fig. 11).

Let us define $`[x]`$ as the largest integer less than or equal to $`x`$, so that

$$[x]x<[x]+1,$$
(135)
and $`\{x\}=x[x][0,1[`$.

Each unfolded pencil of a primitive periodic orbit contains two and only two images of the flux line, since it covers twice the area of the rectangle. The periodic orbits from this pencil parallel to the vector $`(M,N)`$ and going through the images of the flux line (which we shall call the saddle connections) split the pencil of primitive periodic orbits parallel to $`(M,N)`$ into three pencils of same length (see Fig. 11). Only the central strip is affected by the presence of the flux line and any trajectory from this strip gets a phase $`2\pi \alpha `$ (according to $`(\text{131})`$ and since the orbit is primitive). So the winding number of a periodic trajectory is nothing but the repetition number of a periodic orbit belonging to the central pencil.

To compute the widths of the central strip as a function of $`M`$ and $`N`$, let us note that the algebraic distance from an image $`((\zeta _1+2k)a,(\zeta _2+2k^{})b)`$ of the flux line to the saddle-connection linking the points $`(0,0)`$ to $`(2Ma,2Nb)`$ is

$$d=\frac{2𝒜}{l_p}\left(\zeta _2M\zeta _1N+2k^{}M2kN\right);$$
(136)
The two images of the flux line which are inside the pencil are the two nearest points among the all images $`((\zeta _1+2k)a,(\zeta _2+2k^{})b)`$ with integer $`k`$ and $`k^{}`$ that have positive distance to the saddle-connection. They correspond to points such that the distance $`(\text{136})`$ is positive and less than $`2𝒜/l_p`$, or less than $`1`$ in units of $`2𝒜/l_p`$.

Let us set

$$Q_\pm =[ϵ_2M\pm ϵ_1N],ϵ_\pm =\{ϵ_2M\pm ϵ_1N\}.$$
(137)
The four possible values of $`(\zeta _1,\zeta _2)`$ ($`(ϵ_1,ϵ_2)`$, $`(ϵ_1,2ϵ_2)`$, $`(2ϵ_1,ϵ_2)`$, and $`(2ϵ_1,2ϵ_2)`$) give four possible families for the distance (136) :

$`d_i`$ $`=`$ $`2k_i\pm (Q_++ϵ_+)`$
$`d_i`$ $`=`$ $`2k_i\pm (Q_{}+ϵ_{})`$ (138)
where $`k_i`$ is a certain integer which depends on $`k,k^{},M`$ and $`N`$. Among these four families, exactly two points correspond to a distance positive and less than $`1`$ : for instance if $`Q_+`$ and $`Q_{}`$ are even, only $`2k_1+Q_++ϵ_+`$ and $`2k_2+Q_{}+ϵ_{}`$ can be positive and less than one. So we must have $`2k_1+Q_+=0`$ and $`2k_2+Q_{}=0`$ and the two images of the flux line that are in the pencil of primitive periodic orbits are at a distance $`ϵ_+`$ and $`ϵ_{}`$ from the saddle-connection $`(0,0)(2Ma,2Nb)`$ ; if both $`Q_+`$ and $`Q_{}`$ are odd, the two distances are $`1ϵ_+`$ and $`1ϵ_{}`$. Since the width of the central strip is the difference between the two distances, it is in both cases $`|ϵ_1ϵ_2|`$. Dealing in the same way with the case where $`Q_+`$ and $`Q_{}`$ have opposite parity, we get that the width of the central strip in units of $`2𝒜/l_p`$ is

$$\eta =\{\begin{array}{c}\left|ϵ_{}ϵ_+\right|,\text{if }Q_+\text{ and }Q_{}\text{ have the same parity}\hfill \\ \\ \left|1ϵ_{}ϵ_+\right|,\text{if }Q_+\text{ and }Q_{}\text{ have opposite parity}\hfill \end{array}.$$
(139)
Both cases can be summed up in the following formula :

$$\eta =\phi (x,y),$$
(140)
with

$$\begin{array}{c}x=ϵ_2M+ϵ_1N\hfill \\ y=ϵ_2Mϵ_1N\hfill \end{array},$$
(141)
where

$$\begin{array}{c}\phi (x,y)=\left|f(x)f(y)\right|,f(x)=(1)^{[x]}\left(\{x\}\frac{1}{2}\right)\hfill \end{array}$$
(142)
$`f(x)`$ is an even function of period $`2`$ ; if we restrict the study of $`\phi `$ to $`[1,1]\times [1,1]`$ we have

$$\phi (x,y)=\{\begin{array}{c}|x+y|\text{if }xy0\hfill \\ |xy|\text{if }xy0\hfill \end{array},$$
(143)
and the Fourier expansion of $`\phi `$ is

$$\phi (x,y)=2\underset{n=1}{\overset{\mathrm{}}{}}\frac{\left(\mathrm{cos}\pi nx\mathrm{cos}\pi ny\right)^2}{\pi ^2n^2}.$$
(144)
Using $`(\text{141})`$ and $`(\text{144})`$, we obtain that the width of the central strip associated to the orbit $`(2M,2N)`$ is the following

$$\eta =\frac{8}{\pi ^2}\underset{n=1}{\overset{\mathrm{}}{}}\frac{\mathrm{sin}^2(\pi nϵ_2M)\mathrm{sin}^2(\pi nϵ_1N)}{n^2}.$$
(145)
### 4.2 Form factor for the billiard with flux line

The value of the two-point correlation form factor in the diagonal approximation given by Eq. (25), when diffractive contributions have been neglected, still holds for billiards with flux line. But $`𝒜_p`$ now includes the phase factor depending on the repetition number of the trajectory. The density of periodic orbits (10) becomes

$$d_{p.o.}(E)=\frac{1}{2}\underset{pp^+,pp^{}}{}\underset{n=1}{\overset{\mathrm{}}{}}\frac{𝒜_{pn}}{2\pi }\frac{1}{\sqrt{2\pi knl_p}}e^{iknl_pi\frac{\pi }{2}\nu _pi\frac{\pi }{4}}+c.c.$$
(146)
Here we distinguish between two types of orbits. The orbits associated with a primitive orbit $`pp^+`$ have a phase $`\mathrm{exp}(2i\pi n\alpha )`$ for the orbit repeated $`n`$ times, and the total coefficient in the trace formula associated with these orbits is

$$𝒜_{p^+n}=𝒜_{p_1}+𝒜_{p_2}e^{2i\pi n\alpha }+𝒜_{p_3},$$
(147)
where $`𝒜_{p_1}`$, $`𝒜_{p_2}`$ and $`𝒜_{p_3}`$ are the areas covered respectively by the three strips in which the pencil of periodic orbits splits. The orbits associated with a primitive orbit $`pp^{}`$ have the complex conjugate phase $`\mathrm{exp}(2i\pi n\alpha )`$ and their contribution to the trace formula is proportional to

$$𝒜_{p^{}n}=𝒜_{p_1}+𝒜_{p_2}e^{2i\pi n\alpha }+𝒜_{p_3}.$$
(148)
When the terms with same length $`l_p=nl_{pp}`$ are gathered together, we get

$$d_{p.o.}(E)=\underset{p^+}{}\frac{𝒜_{np}}{2\pi }\frac{1}{\sqrt{2\pi kl_p}}e^{ikl_pi\frac{\pi }{2}\nu _pi\frac{\pi }{4}}+c.c.$$
(149)
where

$$𝒜_{np}=𝒜_{p_1}+𝒜_{p_2}\mathrm{cos}(2\pi n\alpha )+𝒜_{p_3}$$
(150)
and the sum $`_{p^+}`$ goes over orbits $`M,N0`$. Equation $`(\text{21})`$ now becomes

$$K_2(\tau )=\underset{pp^+}{}\underset{n=1}{\overset{\mathrm{}}{}}\frac{𝒜_{np}^2}{n^2}\frac{1}{2\pi ^2l_{pp}\overline{d}}\delta \left(l_{pp}\frac{4\pi k\overline{d}\tau }{n}\right).$$
(151)
If $`\eta _p`$ is the width of the central strip expressed in units of $`2𝒜/l_p`$, we have

$$𝒜_{np}=2𝒜(1\eta _p+\eta _p\mathrm{cos}2\pi n\alpha )=2𝒜(12\eta _p\mathrm{sin}^2\pi n\alpha ),$$
(152)
and (using the fact that $`\overline{d}=𝒜/4\pi `$) the 2-point correlation form factor at small $`\tau `$ is a sum of 3 terms

$`K_2(\tau )`$ $`=`$ $`{\displaystyle \frac{8𝒜}{\pi }}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n^2}}{\displaystyle \underset{pp}{}}{\displaystyle \frac{1}{l_{pp}}}\delta \left(l_{pp}{\displaystyle \frac{𝒜k\tau }{n}}\right)`$ (153)
$``$ $`{\displaystyle \frac{32𝒜}{\pi }}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mathrm{sin}^2\pi n\alpha }{n^2}}{\displaystyle \underset{pp}{}}{\displaystyle \frac{\eta _{pp}}{l_{pp}}}\delta \left(l_{pp}{\displaystyle \frac{𝒜k\tau }{n}}\right)`$
$`+`$ $`{\displaystyle \frac{32𝒜}{\pi }}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mathrm{sin}^4\pi n\alpha }{n^2}}{\displaystyle \underset{pp}{}}{\displaystyle \frac{\eta _{pp}^2}{l_{pp}}}\delta \left(l_{pp}{\displaystyle \frac{𝒜k\tau }{n}}\right).`$
The summation over primitive periodic orbits can be done by replacing the sum by an integral, taking into account the density of primitive periodic orbits. If $`ϵ_1`$ and $`ϵ_2`$ are rational numbers

$$ϵ_i=\frac{p_i}{q_i},$$
(154)
where $`p_i`$ and $`q_i`$ are co-prime integers, the width of the central strip $`(\text{145})`$ only depends on the remainder $`r_1`$ of $`M`$ modulo $`q_1`$ and $`r_2`$ of $`N`$ modulo $`q_2`$

$$\eta (r_1,r_2)=\frac{8}{\pi ^2}\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{n^2}\mathrm{sin}^2(\pi \frac{p_1}{q_1}nr_1)\mathrm{sin}^2(\pi \frac{p_2}{q_2}nr_2).$$
(155)
There are $`q_1q_2`$ periodic orbit families

$$M=q_1k+r_1,N=q_2k^{}+r_2,$$
(156)
with $`k,k^{}`$.

To compute sums in Eq. (153) one needs to know the mean density of primitive periodic orbits for each family, $`\rho _{pp,r_1,r_2}(l)`$.

Let $`c`$ be the greatest common divisor of $`q_1`$, $`q_2`$: $`c=(q_1,q_2)`$. If $`(r_1,r_2,c)1`$, then $`M`$ and $`N`$ are not coprime and there is no primitive periodic orbit. In the opposite case it is demonstrated in Appendix A that

$$\rho _{pp,r_1,r_2}(l)=\rho _{pp}(l)\alpha (r_1,r_2),$$
(157)
where $`\rho _{pp}(l)`$ is the mean density for all primitive periodic orbits in the rectangle (cf. Eq. $`(\text{42})`$) with $`M,N>0`$

$$\rho _{pp}(l)=\frac{3l}{4\pi 𝒜},$$
(158)
and

$$\alpha (r_1,r_2)=\frac{1}{q_1q_2_{p|\text{lcf}(q_1,q_2)}(11/p^2)}\underset{p|(q_1,r_1),pq_2}{}(1\frac{1}{p})\underset{p|(q_2,r_2),pq_1}{}(1\frac{1}{p}),$$
(159)
and lcf$`(q_1,q_2)`$ is the least common factor of $`q_1`$, $`q_2`$.

The knowledge of the mean density of periodic orbit families permits the computation of mean values of different quantities depending on families. If $`f(r_1,r_2)`$ is such a quantity its mean value is defined as follows

$$<f>=\underset{\genfrac{}{}{0pt}{}{r_imodq_i}{(r_1,r_2,c)=1}}{}f(r_1,r_2)\alpha (r_1,r_2).$$
(160)
In particular

$`{\displaystyle \underset{pp}{}}{\displaystyle \frac{\eta _{pp}^\beta }{l_{pp}}}\delta \left(l_{pp}{\displaystyle \frac{𝒜k\tau }{n}}\right)`$ (161)
$`={\displaystyle \underset{\begin{array}{c}r_imodq_i\\ (r_1,r_2,c)=1\end{array}}{}}\eta ^\beta (r_1,r_2){\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{1}{l}}\rho _{pp,r_1,r_2}(l)\delta \left(l{\displaystyle \frac{𝒜k\tau }{n}}\right)={\displaystyle \frac{3}{4\pi 𝒜}}<\eta ^\beta >.`$ (164)
The sums over $`n`$ that appear in $`(\text{153})`$ can be computed using the standard formula

$$\underset{n=1}{\overset{\mathrm{}}{}}\frac{\mathrm{cos}2\pi nx}{(2\pi n)^2}=x^2x+\frac{1}{6},\text{for }0x1.$$
(165)
It gives

$$\underset{n=1}{\overset{\mathrm{}}{}}\frac{\mathrm{sin}^2\pi n\alpha }{n^2}=\frac{\pi ^2}{2}\overline{\alpha }(1\overline{\alpha })$$
(166)
and

$$\underset{n=1}{\overset{\mathrm{}}{}}\frac{\mathrm{sin}^4\pi n\alpha }{n^2}=\frac{\pi ^2}{4}\overline{\alpha }$$
(167)
where $`\overline{\alpha }`$ is the fractional part $`\{\alpha \}`$ of the flux through the rectangle when $`0\{\alpha \}1/2`$ and $`\overline{\alpha }=1\{\alpha \}`$ when $`1/2\{\alpha \}1`$.

Using $`(\text{161})`$, $`(\text{166})`$ and $`(\text{167})`$ one concludes that the 2-point correlation form factor for $`\tau 0`$ $`(\text{153})`$ is the following

$$K_2(0)=112\overline{\alpha }(1\overline{\alpha })<\eta >+6\overline{\alpha }<\eta ^2>.$$
(168)
To use this formula it is necessary to know the values of $`<\eta >`$ and $`<\eta ^2>`$. In the case where both $`ϵ_1`$ and $`ϵ_2`$ are irrational non-commensurable quantities the fractional parts $`\{nϵ_2M\}`$ and $`\{nϵ_1N\}`$ cover the whole interval $`[1,1]`$ and $`\eta `$ and $`\eta ^2`$ can be computed by integrating expression $`(\text{143})`$ of $`\eta `$ over the square $`[1,1]\times [1,1]`$. Simple calculations show that in this case

$$<\eta >=\frac{1}{3},$$
(169)
$$<\eta ^2>=\frac{1}{6}.$$
(170)
Therefore when the coordinates of the flux line are non-commensurable with the corresponding sides

$$K_2(0)=13\overline{\alpha }+4\overline{\alpha }^2.$$
(171)
In Appendix B, it is shown that for all rational $`ϵ_i`$

$$<\eta >=\frac{1}{3},$$
(172)
like in the irrational case. The average $`<\eta ^2>`$ is more difficult to compute analytically: we have found an analytical expression only when $`q_1`$ divides $`q_2`$ (or similarly $`q_2`$ divides $`q_1`$).

Though the general formula for $`<\eta ^2>`$ is cumbersome, the computation of the $`<\eta ^2>`$ for rational $`ϵ_1`$ and $`ϵ_2`$ can easily be done numerically using Eqs. $`(\text{141})`$ and $`(\text{143})`$. For small denominators the values of $`<\eta ^2>`$ are given in Table 1.

To check the obtained formulas we have computed numerically 1500 first energy levels for the rectangular billiard with sides $`a=4`$ and $`b=\pi `$ and the flux line with coordinates (from low left corner) $`x_0=5a/9`$ and $`y_0=11b/20`$. The typical picture of $`K_2(t)`$ is shown in Fig. 12.

As for the triangular billiards discussed in previous Sections we extrapolated $`K_2(\tau )`$ to small $`\tau `$ using the simple expression (125). The results for different values of the flux are presented in Fig. 13.

We also check the following more suitable fit (which obeys the condition (127) when $`c=(1b)^2`$)

$$K_2(\tau )=\{\begin{array}{c}b+c\tau ,\text{when }\tau <(1b)/c\\ 1,\text{when }\tau >(1b)/c\end{array}.$$
(173)
It gives practically the same results. The existing numerical precision does not permit to distinguish these 2 fits.

In the case where $`x_0=5a/9`$ and $`y_0=11b/20`$, simple calculations show that $`<\eta ^2>=4867/29160`$; Eq. $`(\text{168})`$ gives the expected value of $`K_2(0)`$

$$K_2(0)=1\frac{14573}{4860}\overline{\alpha }+4\overline{\alpha }^2.$$
(174)
which corresponds to the solid curve in fig. 13 (Note that the coefficient of $`\overline{\alpha }`$ equals approximately $`2.99`$ and is practically indistinguishable from the coefficient 3 for irrational $`ϵ_i`$).

Similarly as for triangular billiards there is a small difference between the theoretical and numerical curves. For triangular billiards where more levels are available this difference slowly decreases with energy. We expect the same behavior also for rectangular billiards with a flux line.

## 5 Conclusion

In this paper we have obtained explicit expressions of the 2-point correlation form factor $`K_2(\tau )`$ in the limit $`\tau 0`$ for a few typical examples of pseudo-integrable billiards : triangular billiards in the shape of right triangles with one angle equals $`\pi /n`$, and rectangular billiards with a flux line. The obtained values of $`K_2(0)`$ differ from standard examples of spectral statistics (the Random matrix theory and the Poisson statistics), which confirm analytically the peculiarities of spectral statistics of pseudo-integrable systems. The calculations have been performed by analysing analytically the properties of classical periodic orbits of the systems considered.

In order to elucidate further the special properties of spectral statistics of polygonal billiards, it would be of interest to compute $`K_2(0)`$ for generic triangular billiards without the Veech structure. Moreover, we have taken into account only the diagonal terms and, consequently, were able to obtain only $`K_2(0)`$. The computation of the next terms in the expansion of $`K_2(\tau )`$ in powers of $`\tau `$ should include the exact ressumation of singular contributions, coming from the diffraction close to the optical boundaries. The solutions of these problems require the development of new methods beyond the ones used in this paper.

## Appendix A

Periodic orbits in a rectangle with sides $`a`$, $`b`$ are determined by 2 integers $`M`$ and $`N`$ which count the difference of coordinates of initial, $`(x_i,y_i)`$, and final, $`(x_f,y_f)`$, points

$$x_f=x_i+2aM,y_f=y_i+2bN.$$
(175)
The length of the periodic orbit is the geometrical length of this vector

$$L_p=\sqrt{(2aM)^2+(2bN)^2}.$$
(176)
The mean cumulative density and the corresponding quantity for primitive periodic orbits (when $`M`$, $`N`$ are co-prime integers) can be computed as for the square billiard (see Eqs. (38) and (42)). When $`l\mathrm{}`$ and if only positive $`M`$ are considered

$$N(L_p<l)\frac{\pi l^2}{8ab},$$
(177)
and

$$N_{pp}(L_p<l)\frac{3l^2}{4\pi ab}.$$
(178)
The purpose of this Appendix is the computation of the mean cumulative density of primitive periodic orbits for periodic orbit families when

$$Mr_1modq_1,Nr_2modq_2.$$
(179)
The asymptotics of $`N_{pp}(L_p<l)`$ when $`l\mathrm{}`$ is related with the behavior at small $`x`$ of the $`\mathrm{\Theta }`$-function associated with these periodic orbits

$$\mathrm{\Theta }(x)=\underset{pp}{}e^{xL_p^2}.$$
(180)
If

$$\mathrm{\Theta }(x)\frac{C}{x^\gamma },\text{when }x0,$$
(181)
then

$$N_{pp}(L_p<l)\frac{C}{\gamma \mathrm{\Gamma }(\gamma )}l^{2\gamma },\text{when }l\mathrm{}.$$
(182)
We are interested in the following $`\mathrm{\Theta }`$-function

$$\mathrm{\Theta }(x)=\underset{M,N=\mathrm{}}{\overset{\mathrm{}}{}}e^{x((2aM)^2+(2bN)^2)},$$
(183)
where the summation is performed over all integers $`M`$, $`N`$ with the following constraints

$$(M,N)=1,Mr_1modq_1,Nr_2modq_2.$$
(184)
Note that both positive and negative values of $`M`$, $`N`$ are considered. When only positive $`M`$ are taken into account the formulas below have asymptotically factor $`1/2`$.

To impose the restriction $`Mrmodq`$ it is convenient to introduce the $`\delta `$-function

$$\delta _{t,q}=\{\begin{array}{cc}1,& \text{if }t0modq\\ 0,& \text{otherwise}\end{array}.$$
(185)
Its explicit form may be the following

$$\delta _{t,q}=\frac{1}{q}\underset{k=0}{\overset{q1}{}}e^{2\pi ikt/q}.$$
(186)
As in Section 3.1 the condition $`(M,N)=1`$ can be taken into account by the inclusion-exclusion principle

$$\underset{(M,N)=1}{}f(M,N)=\underset{M,N=\mathrm{}}{\overset{\mathrm{}}{}}\underset{t=1}{\overset{\mathrm{}}{}}f(Mt,Nt)\mu (t),$$
(187)
where $`\mu (t)`$ is the Möbius function equal $`(1)^n`$ if $`t`$ is a product of $`n`$ distinct primes, 0 if $`t`$ contains a squared factor, and $`\mu (1)=1`$.

Combining all the necessary restrictions one finds the final expression for the $`\mathrm{\Theta }`$-function (183)

$`\mathrm{\Theta }(x)`$ $`=`$ $`{\displaystyle \frac{1}{q_1q_2}}{\displaystyle \underset{M,N=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \underset{k_imodq_i}{}}{\displaystyle \underset{t=1}{\overset{\mathrm{}}{}}}\mu (t)e^{4xt^2(M^2a^2+N^2b^2)}`$ (188)
$`\times e^{2\pi ik_1(Mtr_1)/q_1+2\pi ik_2(Ntr_2)/q_2}.`$
Using the Poisson summation formula

$$\underset{M=\mathrm{}}{\overset{\mathrm{}}{}}e^{xM^2+2\pi iyM}=\sqrt{\frac{\pi }{x}}\underset{M=\mathrm{}}{\overset{\mathrm{}}{}}e^{\pi ^2(M+y)^2/x},$$
(189)
one obtains that

$`\mathrm{\Theta }(x)`$ $`=`$ $`{\displaystyle \frac{\pi }{4abxq_1q_2}}{\displaystyle \underset{k_imodq_i}{}}{\displaystyle \underset{t=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mu (t)}{t^2}}e^{2\pi ik_1r_1/q_12\pi ik_2r_2/q_2}`$ (190)
$`\times `$ $`{\displaystyle \underset{M=\mathrm{}}{\overset{\mathrm{}}{}}}e^{\pi ^2(M+k_1t/q_1)^2/(4xt^2a^2)}{\displaystyle \underset{N=\mathrm{}}{\overset{\mathrm{}}{}}}e^{\pi ^2(N+k_2t/q_2)^2/(4xt^2b^2)}.`$
When $`x0`$ the dominant contribution comes from terms with zero exponent, i.e. from terms with

$$M+\frac{k_1t}{q_1}=0,\text{and}N+\frac{k_2t}{q_2}=0,$$
(191)
or $`k_it0modq_i`$. The asymptotics of the $`\mathrm{\Theta }`$-function is therefore the following

$$\mathrm{\Theta }(x)=\frac{\pi }{4abx}\frac{1}{q_1q_2}\underset{k_imodq_i}{}\underset{t=1}{\overset{\mathrm{}}{}}\delta _{k_1t,q_1}\delta _{k_2t,q_2}\frac{\mu (t)}{t^2}e^{2\pi ik_1r_1/q_12\pi ik_2r_2/q_2}.$$
(192)
Using the representation (186) for these $`\delta `$-functions and performing the summation over $`k_i`$ one gets that when $`x0`$

$$\mathrm{\Theta }(x)=\frac{\pi }{4abx}F(r_1,r_2),$$
(193)
where

$$F(r_1,r_2)=\frac{1}{q_1q_2}\underset{l_imodq_i}{}\underset{t=1}{\overset{\mathrm{}}{}}\delta _{l_1tr_1,q_1}\delta _{l_2tr_2,q_2}\frac{\mu (t)}{t^2}.$$
(194)
From Eq. (182) one concludes that the asymptotics of the mean cumulative density of primitive periodic family (184) (with $`M,N>0`$, i.e. with a factor $`1/4`$) is

$$N_{pp}(L_p<l)=\frac{\pi l^2}{16ab}F(r_1,r_2).$$
(195)
To perform the summation over $`t`$ in Eq. (194) it is necessary to know the number of solutions of two equations

$$l_1tr_1modq_1,l_2tr_2modq_2.$$
(196)
It is well known (and can be easily checked) that the number of solutions of the equation $`axbmodq`$ depends on the greatest common divisor of $`a`$ and $`q`$, $`(a,q)=d`$. If $`d=1`$ there is 1 solution, $`xba^1modq`$. If $`d>1`$ and $`db`$ there is no solutions. If $`d|b`$ there is one solution, $`x_0=(b/d)(a/d)^1mod(q/d)`$ and consequently, there are $`d`$ solutions modulo $`q`$: $`x_j=x_0+(q/d)j`$, $`j=0,\mathrm{}d1`$. Therefore

$`F(r_1,r_2)`$ $`=`$ $`{\displaystyle \frac{1}{q_1q_2}}{\displaystyle \underset{t=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mu (t)}{t^2}}(q_1,t)(q_2,t)\delta _{(q_1,t),r_1}\delta _{(q_2,t),r_2}`$ (197)
$`=`$ $`{\displaystyle \frac{1}{q_1q_2}}{\displaystyle \underset{\genfrac{}{}{0pt}{}{d_1|(q_1,r_1)}{d_2|(q_2,r_2)}}{}}{\displaystyle \underset{\genfrac{}{}{0pt}{}{(q_1,t)=d_1}{(q_2,t)=d_2}}{}}{\displaystyle \frac{\mu (t)}{t^2}}d_1d_2.`$
Terms corresponding to $`(d_1,d_2)>1`$ give a $`0`$ contribution to the sum, since in that case $`q_1,q_2,d_1`$ and $`d_2`$ have a common factor, which contradicts condition $`(M,N)=1`$. The sum $`(\text{197})`$ can therefore be restricted to $`(d_1,d_2)=1`$, and the sum over $`t`$ is now a sum over $`t^{}`$ where $`t=d_1d_2t^{}`$. Let us denote by $`P`$ the product of the prime factors of $`q_1`$ that do not divide $`c`$ and by $`P^{}`$ the product of the prime factors of $`q_2`$ that do not divide $`c`$. Now

$$(c,P)=(c,P^{})=(P,P^{})=1.$$
(198)
If $`p`$ is a prime dividing $`d_1`$ and $`c=(q_1,q_2)`$, then $`p`$ divides $`d_1=(q_1,t)`$. Since it divides $`c`$ it also divides $`q_2`$, so $`p|(q_2,t)`$ and $`(d_1,d_2)1`$, which is impossible. So the prime divisors of $`d_1`$ have to be taken among the divisors of $`P`$, and in the same way the prime divisors of $`d_2`$ have to be taken among the divisors of $`P^{}`$. Similarly, we can check that if $`p`$ is a prime divisor of $`c`$ which divides $`t^{}`$, $`p`$ divides $`d_1=(q_1,t)`$, so $`d_1`$ would contain a prime factor of $`c`$, which is impossible; and if $`p`$ is a prime factor of $`P`$ or $`P^{}`$, $`p`$ divides both $`d_1`$ and $`d_2`$. So the sum over $`t^{}`$ must be restricted to $`t^{}`$ which do not contain any prime divisor of $`q_1`$ or $`q_2`$. As $`\mu (ab)=\mu (a)\mu (b)`$ for co-prime $`a`$ and $`b`$, one gets

$$F(r_1,r_2)=\frac{1}{q_1q_2}\underset{d_1,d_2,t^{}}{}d_1d_2\frac{\mu (d_1)\mu (d_2)\mu (t^{})}{d_1^2d_2^2t_{}^{}{}_{}{}^{2}}$$
(199)
where the sum is taken over all $`d_1,d_2,t^{}`$ verifying $`d_1|(P,r_1)`$, $`d_2|(P^{},r_2)`$, $`t^{}q_1`$, $`t^{}q_2`$. Using the identity

$$\underset{p|k}{}(1\frac{1}{p^s})=\underset{\delta |k}{}\frac{\mu (\delta )}{\delta ^s},$$
(200)
we get

$`F(r_1,r_2)`$ $`=`$ $`{\displaystyle \frac{1}{q_1q_2}}{\displaystyle \underset{pq_1,pq_2}{}}(1{\displaystyle \frac{1}{p^2}}){\displaystyle \underset{p|(P,r_1)}{}}(1{\displaystyle \frac{1}{p}}){\displaystyle \underset{p|(P^{},r_2)}{}}(1{\displaystyle \frac{1}{p}})`$ (201)
$`=`$ $`{\displaystyle \underset{\text{all }p}{}}(1{\displaystyle \frac{1}{p^2}})\alpha (r_1,r_2)={\displaystyle \frac{6}{\pi ^2}}\alpha (r_1,r_2),`$
where

$$\alpha (r_1,r_2)=\frac{1}{q_1q_2_{p|\text{lcf}(q_1,q_2)}(11/p^2)}\underset{p|(q_1,r_1),pc}{}(1\frac{1}{p})\underset{p|(q_2,r_2),pc}{}(1\frac{1}{p}),$$
(202)
and lcf$`(q_1,q_2)`$ is the least common factor of $`q_1`$, $`q_2`$ (lcf$`(q_1,q_2)=q_1q_2/c`$).

It is also instructive to check directly that $`\alpha (r_1,r_2)`$ are normalized correctly

$$\underset{\genfrac{}{}{0pt}{}{r_imodq_i}{(r_1,r_2,c)=1}}{}\alpha (r_1,r_2)=1,$$
(203)
where as above $`c=(q_1,q_2)`$. We use once more the inclusion-exclusion principle

$$\underset{\genfrac{}{}{0pt}{}{r_imodq_i}{(r_1,r_2,c)=1}}{}f(r_1,r_2)=\underset{t|c}{}\mu (t)\underset{\genfrac{}{}{0pt}{}{r_imodq_i}{t|r_i}}{}f(r_1,r_2).$$
(204)
If we set

$$D=q_1q_2\underset{p|\text{lcf }(q_1,q_2)}{}(1\frac{1}{p^2}),$$
(205)
we get, with $`(\text{202})`$ and $`(\text{204})`$,

$`{\displaystyle \underset{\genfrac{}{}{0pt}{}{r_imodq_i}{(r_1,r_2,c)=1}}{}}\alpha (r_1,r_2)`$ $`=`$ $`{\displaystyle \frac{1}{D}}{\displaystyle \underset{t|c}{}}\mu (t){\displaystyle \underset{r_1(t|r_1)}{}}{\displaystyle \underset{r_2(t|r_2)}{}}{\displaystyle \underset{\delta _1|(P,r_1)}{}}{\displaystyle \frac{\mu (\delta _1)}{\delta _1}}{\displaystyle \underset{\delta _2|(P^{},r_2)}{}}{\displaystyle \frac{\mu (\delta _2)}{\delta _2}}`$ (206)
$`=`$ $`{\displaystyle \frac{1}{D}}{\displaystyle \underset{t|c}{}}\mu (t){\displaystyle \underset{\delta _1|P}{}}{\displaystyle \underset{\delta _2|P^{}}{}}{\displaystyle \underset{r_1(t\delta _1|r_1)}{}}{\displaystyle \underset{r_2(t\delta _2|r_2)}{}}{\displaystyle \frac{\mu (\delta _1)}{\delta _1}}{\displaystyle \frac{\mu (\delta _2)}{\delta _2}}`$
$`=`$ $`{\displaystyle \frac{1}{D}}q_1q_2{\displaystyle \underset{t|c}{}}{\displaystyle \frac{\mu (t)}{t^2}}{\displaystyle \underset{\delta _1|P}{}}{\displaystyle \frac{\mu (\delta _1)}{\delta _1^2}}{\displaystyle \underset{\delta _2|P^{}}{}}{\displaystyle \frac{\mu (\delta _2)}{\delta _2^2}}.`$
Here we have used the fact that if $`t|c`$ and $`\delta |P`$, $`t`$ and $`\delta `$ have no common factor. In the above sums it is always understood that the summation over $`r_i`$ goes only from $`r_i=0`$ to $`q_i1`$. But the last sum in Eq. (206) exactly equals $`D`$ because $`cPP^{}=\text{lcf}(q_1,q_2)`$ and Eq. (203) holds.

## Appendix B

In the same way one can compute the mean value of $`\eta `$ defined in Eq. (155)

$`<\eta >={\displaystyle \underset{\genfrac{}{}{0pt}{}{r_imodq_i}{(r_1,r_2,c)=1}}{}}\eta (r_1,r_2)\alpha (r_1,r_2)={\displaystyle \frac{8}{\pi ^2D}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n^2}}{\displaystyle \underset{t|c}{}}\mu (t)`$ (207)
$`{\displaystyle \underset{r_1(t|r_1)}{}}{\displaystyle \underset{r_2(t|r_2)}{}}{\displaystyle \underset{\delta _1|(P,r_1)}{}}{\displaystyle \underset{\delta _2|(P^{},r_2)}{}}{\displaystyle \frac{\mu (\delta _1)}{\delta _1}}{\displaystyle \frac{\mu (\delta _2)}{\delta _2}}\mathrm{sin}^2\left(\pi n{\displaystyle \frac{r_1}{q_1}}\right)\mathrm{sin}^2\left(\pi n{\displaystyle \frac{r_2}{q_2}}\right).`$
Since

$$\underset{t\delta |r}{}\mathrm{sin}^2\pi n\frac{r}{q}=\frac{q}{2t\delta }\left(1\delta _{nt\delta ,q}\right),$$
(208)
one obtains

$`<\eta >`$ $`=`$ $`{\displaystyle \frac{2q_1q_2}{\pi ^2D}}{\displaystyle \underset{t|c}{}}{\displaystyle \frac{\mu (t)}{t^2}}{\displaystyle \underset{\delta _1|P}{}}{\displaystyle \frac{\mu (\delta _1)}{\delta _1^2}}{\displaystyle \underset{\delta _2|P^{}}{}}{\displaystyle \frac{\mu (\delta _2)}{\delta _2^2}}`$ (209)
$`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n^2}}\left(1\delta _{nt\delta _1,q_1}\right)\left(1\delta _{nt\delta _2,q_2}\right).`$
The sum over $`n`$ includes 4 terms. The first is the sum over all $`n`$

$$\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{n^2}=\frac{\pi ^2}{6}.$$
(210)
The second sum has the restriction that $`n=(q_1/t\delta _1)m`$ and

$$\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{n^2}\delta _{nt\delta _1,q_1}=\frac{\pi ^2}{6}\left(\frac{t\delta _1}{q_1}\right)^2.$$
(211)
The third sum is the same but with the substitution $`12`$. The fourth sum incorporates two restrictions, $`nt\delta _10modq_1`$ and $`nt\delta _20modq_2`$. Remembering the definition of $`P`$ and $`P^{}`$ (see (198)) one concludes that in this last case the restriction is $`n=(cPP^{}/(t\delta _1\delta _2))m`$ and

$$\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{n^2}\delta _{nt\delta _1,q_1}\delta _{nt\delta _2,q_2}=\frac{\pi ^2}{6}\left(\frac{t\delta _1\delta _2}{cPP^{}}\right)^2.$$
(212)
Performing the summation over $`\delta _i`$ and $`t`$ in Eq. (209) one notes that all three last sums will have as a factor

$$\underset{\delta _1|P}{}\mu (\delta _1)\text{o}r\underset{\delta _2|P}{}\mu (\delta _2).$$
(213)
But for any $`K2`$ we have

$$\underset{\delta |K}{}\mu (\delta )=0.$$
(214)
Since $`q_1`$ and $`q_2`$ are greater than $`1`$, it is impossible that simultaneously $`c=P=1`$, or $`c=P^{}=1`$, the terms $`(\text{213})`$ equal zero. Therefore only the term (210) survives and $`(\text{209})`$ gives

$$<\eta >=\frac{1}{3D}q_1q_2\underset{t|c}{}\frac{\mu (t)}{t^2}\underset{\delta _1|P}{}\frac{\mu (\delta _1)}{\delta _1^2}\underset{\delta _2|P^{}}{}\frac{\mu (\delta _2)}{\delta _2^2}.$$
(215)
These sums are exactly equal to $`D`$ and finally we get

$$<\eta >=\frac{1}{3}.$$
(216)