Singular statistics revised

In the paper we analyze the"singular statistics"of pseudointegrable Seba billiards and show that taking into account growing number of resonances one observes the transition from"semi-Poissonian"-like statistics to Poissonian. This observation is in agreement with an argument that a classical particle does not feel a point perturbation. However, our findings contradict results reported earlier (P. Seba, Phys. Rev. Lett. 64, 1855 (1990)).

(a) (b) Figure 1. Physical scatterers of the same scattering length but different radii introduced into the rectangular billiard. Black circle shows the scatterer, the radius of the dashed circle is equal to the scattering length.

Introduction
The singular perturbed square billiard, also calledŠeba billiard [1] by several authors, is one of the key models for quantum chaotic systems. While the unperturbed square billiard with side ratio 1/( √ 5 − 1) shows Poissonian level-spacings statistics [2] and the Sinai billiard, being proved to be a fully chaotic system [3], exhibits GOE-statistics [4], the "intermediate" case of a singular perturbed billiard is expected to demonstrate some transient behavior [1,5,6,7]. However it was reported, that the billiard with a point perturbation can exhibit "fully developed quantum chaos" [1,6]. It may seem strange, since a point perturbation has almost no influence on the classical phase space of the billiard. Thus it is natural to assume that in the semiclassical limit the billiard with a point perturbation shows a similar statistics to an unperturbed one.
The proposed explanation of the given paradox was based on the argument that for zero-range perturbation for any wavelength one can never reach the limit of the classical billiard with a point-perturbation, since the wavelength is finite while the perturbation radius is zero. Therefore the quantum system becomes chaotic while its classical analog is almost integrable. This argument seemingly was justified experimentally [6].
In the presented report we show that the level-spacings statistics ofŠeba billiards actually tends to Poissonian when the number of taken eigenvalues tends to infinity. These findings are in accordance with the intuitive "classical" argument given above, but are in contradiction with previous theoretical [1,7] and experimental [6] results. For a narrow window of eigenvalues some conclusions of Refs. [1,7] remain valid, however one cannot directly apply the results to a wide eigenvalue range. This discrepancy traces back to the procedure applied by the authors to take care of the singularity of the Green function, by replacing the "bare" coupling constant by a renormalized one absorbing the infinity. Although this renormalization technique is standard in quantum electrodynamics, it is not appropriate to compute the spectrum of aŠeba billiard. In this paper we present the suitable renormalization procedure.
Let us turn to the experimental microwaveŠeba billiard [6] and explain why the interpretation of the obtained distributions was erroneous. Theoretically in all cases the scatterer was treated as a point scatterer, which means physically that the characteristic (a) (b) (c) Figure 2. Sketch of the level-spacings statistics evolution computed for a given fixed number of eigenvalues of a billiard with a small perturbation whose radius a is significantly smaller than its scattering wavelength β. The direction from left to right corresponds to increasing eigenvalues. The "semi-Poissonian" statistics (a) corresponds to the range a ≪ β λ, the "Poissonian" statistics (b) corresponds to the range a ≪ λ ≪ β and the "GOE" statistics (c) corresponds to the range a ∼ λ ≪ β. Here λ is the characteristic wavelength. Here λ is the characteristic wavelength. The names of distributions are written in parenthesis since the corresponding plotted curves keep the essentials of these statistics but may differ from real distributions.
wavelength λ of the field inside the cavity is much larger than the radius a of the scatterer. At the same time, the scattering length β of the given scatterer is, generally speaking, a free parameter depending on the internal structure of the scatterer (e. g. given material of the coating, radius of the metallic core etc). We show below that the influence of the point scatterer is significant when λ β and vanishes when λ ≪ β in accord with a classical limit. To treat the scatterer as a point perturbation we obligatory should require a ≪ λ. Combining the last two estimations we conclude that to cover experimentally the classical limit of aŠeba billiard one needs to create a scatterer with a ≪ β (see figure 1 (b)).
However, experimentally one often has the situation shown in figure 1 (a), i. e. a ∼ β. This means that in the regime λ ∼ β the corresponding billiard should be treated as a quantum Sinai billiard, but not aŠeba billiard. In this case the "classical" limit with a point perturbation can not be achieved.
The expected "experimental" evolution of the level-spacings statistics computed for a given number of resonances taken at different frequencies is plotted in Fig. 2. It has been assumed that the radius of a small scatterer a ≪ β, but remains finite. Figures 2 (a), 2 (b) correspond to theŠeba billiard approximation when the radius of the perturbation can be neglected. Figure 2 (c) shows the "GOE" statistics of the Sinai billiard in the regime where the wavelength of the field is comparable with the radius of the perturbation.
In what follows we take the limit a → 0 which means that figure 2 (c) can not be reproduced within the framework of the considered approach. Thus we restrict ourself to the mathematical model of a point perturbation as it was done by the authors of Refs. [1,7].

Point perturbation of the billiard
Let us now come to the theoretical and numerical study of theŠeba billiard and show why the previous treatments have given incorrect results. Following [8,9,1] let us first introduce a point perturbation of the billiard at point R. Basically we will construct the "self-adjoint extension" of the unperturbed "Hamiltonian". This approach was already previously used in [10].
For the unperturbed billiard the eigenfunctions ψ and eigenvalues k 2 = k 2 n obey the equation If ψ(R) = 0, the corresponding states do not feel the perturbation, thus these eigenfunctions and the corresponding eigenvalues are identical for the unperturbed and perturbed billiards. Next we assume that the perturbed eigenfunctions G obey the equation outside of the scatterer of the radius a. In what follows we assume a → 0. We shall see below that G is nothing but the Green function of the unperturbed system. To recover an appropriate boundary condition at the perturbation point let us consider the asymptotics of the function G(r, R; k) outside of the scatterer when r tends to R. Rewriting (2) in cylindrical coordinates we obtain where ρ = |r − R| and ϕ is the angle between the vector r − R and x-axis going along a side of the rectangle. We require G(r, R; k) to be cylindrically symmetric in the vicinity of the point r = R, which assumes that the scatterer is cylindrically symmetric. Thus we obtain 1 ρ The solution of the last equation is where J 0 and Y 0 are Bessel functions of the first and the second kind respectively and c 1 , c 2 are constants parametrically depending on R and k. The last equality should be understood in asymptotic sense only, since J 0 and Y 0 do not obey the proper conditions at the outer boundary of the billiard.
Using the asymptotic form of Y 0 (z → 0) [11] Y 0 (z → 0) = 2 π ln where γ is the Euler constant, and the "identity" we find that in the limit a → 0 the function G(r, R; k) obeys the equation if we assume c 2 (R; k) = 1/4. Another choice of the constant c 2 (R; k) would only lead to a different normalization. Taking into account the boundary conditions for the function G(r, R; k) at the outer boundary of the billiard and (8) we conclude that the perturbed eigenfunctions are the Green functions of the unperturbed billiard.
We are now going to derive the proper boundary condition at the perturbation point. First we separate the Green function into its regular and singular part, respectively.
Combining (5), (6) we obtain where z = kρ, b is some arbitrary length, and Here ξ b (R; k) is the renormalized Green function. In figure 3 we illustrate the singularity of the Green function near the perturbation point.
Let us now consider the rectangle with a small pricked circle of radius a whose center is situated at the point R. We denote it Ω a . Then we consider the linear space of functions consisting of two subspaces: (1) the subspace of functions f (1) (r, R) vanishing at the outer boundary and possessing the asymptotics where B and ξ are some constants, and (2) the subspace of regular functions f (2) (r) vanishing at the outer boundary of the billiard such that f (2) (R) = 0. G(r, R; k) belongs to the subspace (1) because of its asymptotic behavior (9). The unperturbed eigenfunctions ψ(r) belong to the subspace (2). Now we study the action of the the operator −∆ = −∇ 2 on the space of functions defined above. The requirement of hermicity gives where f (i) , f (j) with i, j = 1, 2 are arbitrary functions taken from the subspaces (1) and (2), respectively, and By means of the Green's theorem we obtain Equation (14) shows that (12) holds automatically if f (i) and f (j) both belong to the subspace (2). Assume now that i = 1 and j = 2 or vice versa. Using the asymptotics (9) we see that (12) again holds automatically for any f (1) and f (2) . Thus the case i, j = 1 implies the only nontrivial condition superimposed by the hermicity requirement of the constructed operator. Let us take two functions f Substituting (15) into (14) we find The equality (16) must hold for any values of B i , ξ i . This leads to the conclusion that for all functions from the subspace (2) the constant ξ in (11) is real and the same. Let us chose a certain value ξ = −D. Then the boundary condition at the perturbation point reads: provided that the length b is fixed. Comparing (9) with (11) we find that for the perturbed eigenfunctions G(r, R; k) the constant ξ is equal to ξ b (R; k). Then (17) gives Substituting (10) into (18) we obtain The proper boundary condition cannot depend on the arbitrary length b, but rather should depend on a parameter, characterizing the inner structure of the perturbation. Therefore the length b should be canceled in (19) by a proper choice of D. This can be achieved by the following choice: where β is the scattering length of the perturbation. The value of the length β can not be obtained from the consideration above since we cut out the area containing the perturbation, thereby loosing the information on it. Thus we draw the conclusion that the scattering length is the only parameter describing the perturbation in the limit ka ≪ 1. Substituting (20) in (18) and using (10), we find The perturbed part of the spectrum may now be obtained from the solutions k 2 = k 2 n of (21). Replacing b by β, k by k n in (9) and using the equality ξ β (R, k n ) = 0 we find the asymptotic expansion of the perturbed eigenfunction corresponding to the eigenvalue k 2 n : The leading term of the asymptotics (22) becomes zero when ρ = β. This fact can be used to determine experimentally a scattering length of a given perturbation. We note that the proper definition of the scattering length was missing in Refs. [1,12]. This has led to the deficiencies discussed above. However in the monograph [8], devoted to point perturbations, the scattering length in two-dimensional problems was properly introduced.
Several conclusions on the level-spacings distribution can be drawn already from (21). Indeed, from (9) and (10) we conclude that c 1 (R; k) has the same poles as the Green function of the unperturbed billiard. From (21) we obtain that when k tends to infinity, the eigenvalues of the perturbed billiard approach the eigenvalues of the unperturbed one. Indeed, close to the eigenvalue k n the function c 1 (R, k) may be approximated by const/(k 2 − k 2 n ) whence follows: const Then k 2 − k 2 n = −2π const/ ln(kβ/2). In the limit k → ∞ we recover the original spectrum of the billiard! Since the statistics of inter-level spacings for the unperturbed rectangular billiard with chosen side ratio is Poissonian [2], we conclude that the same statistics for high-lying eigenvalues of theŠeba billiard is also Poissonian. This is the most important conclusion of the paper, which contradicts the prediction given in [1].

Ewald's representation of the renormalized Green function
For an explicit calculation of the spectrum of the perturbed billiard from (21) an expression for ξ β (R; k) is needed, which for the general cases by no means is a trivial task. For the rectangle it can be obtained by an application of Ewald's method. The derivation is technical, and anybody not interested in the details may proceed directly to (55). In the paper [7] in (21) the logarithmic dependence of the dimensionless scattering strength η β (k) was missed. Therefore the proper spectral statistics of theŠeba billiard still has to be computed. In this section we explain the numerical procedure, based on Ewald's method [13,14,15,16,17,18], which allows us to compute the renormalized Green function and then, from (21), the spectral statistics.
We start from the eigenfunction representation of the Green function for the free billiard where R = (x ′ , y ′ ), d x and d y are the two sides of the rectangle. When x → x ′ and y → y ′ then the series (24) diverges logarithmically. This is just another manifestation of the well-known singularity of the Green function for r → R, see (9), which is a local feature and does not depend on outer boundary conditions. This suggests that the eigenfunctions representation is not the appropriate choice to compute the renormalized Green function, but that the images representation [16] might be preferable. The images representation of the Green function reads where and G f (r, R; k) is the Green function for the two-dimensional plane (see figure 4). The images representation is much better suited to compute the renormalized Green function, since when r tends to R, the term n = m = 0, s 1 = s 2 = 0 is the only one in the representation logarithmically tending to infinity. Now the divergency can be subtracted analytically. However the image representation does not solve the problem yet since it converges absolutely only if Im k > 0. To overcome this obstacle Ewald [13] proposed the dual representation keeping features of the images as well as the eigenmodes representation.
Below we follow the works [14,15,16]. Let us first find a convenient representation for G f . To this end we consider the following initial-value problem for the function g(t; r, R, k) Then G f can be written as The contour C should start at t = 0 and tend to infinity in such a way that g tends to zero. Obviously where the heat kernel K(t; r, R) can be found by a separation of variables Finally we obtain The simplest contour of the integration is the imaginary half-line going from zero to i∞ (contour C 1 in figure 5 (a) ). Using this contour we conclude ( [11], Entry B.187 (2)) that where H (1) 0 is the Hankel function of the first kind. However, the asymptotic behavior (9) becomes hidden in this representation. To recover the asymptotic behavior we will use another contour of the integration. We see that for small values of |t| the best convergency provides an interval lying on the real axis from zero to some positive value t Ew (see figure 5 (a) ). We will call this value the Ewald parameter. From the other side the best convergency for large values of |t| would provide the half-line going from some negative value (we choose it to be equal to −t Ew ) to −∞. What remains is to connect these parts to make a contour. We connect them by a half-circle C 2 . Finally the constructed contour is equivalent to C 1 since the integrals along the quarter-circles C 3 , C 4 (plotted by dashed lines in figure 5 (a) ) tend to zero when the radius of C 3 tends to infinity and the radius of C 4 tends to zero. Using the constructed contour we can obtain the asymptotics of the type (9) from the representation (33). The developed technic will be used further to compute the renormalized Green function. We write where Introducing the notations u = (r − R) 2 /(4t Ew ), v = k 2 t Ew , we write G (1) f as follows When r tends to R, then u tends to zero. To compute the asymptotics of (37) for u → 0 we make the following transformations: where (see [11] (Entry 8.367.12)). Since g 0 (v) has no singularity at t = 0 it can be written as where C 0 is a half-circle of the unit radius. Since G (2) f has no singularity when r → R, the leading term of its asymptotics is Using (38), (40), (41) we find The last equality does not depend on the choice of t Ew . To prove it we can take v as the independent parameter, then u = k 2 (r − R) 2 /(4v). Differentiation of the last equality with respect to v gives zero. The reasonable choice of v should not lead to exponentially large values of g 0 (v). Indeed due to (42) such a large contribution is somehow artificial since it is annihilated by G (2) f . Thus it is natural to take v = 1. Then (42) gives This calculation has demonstrated that it is important to divide the free Green function in two parts: G (1) f and G (2) f . The first part describes the space singularity, and the second part makes the contribution into the regular part of the asymptotics.
Let us turn now to the image representation (27) of the Green function of the rectangular cavity. We again write the free Green function G f in the form where G (1) f is defined by (37) and The best choice of the contour C 5 will be discussed later. We shall see that the choice used to compute the asymptotics (43) does not fit. Substituting (44) into (27) we obtain where To improve the convergency of series for G (2) we use the identity which can be proved by applying the Poisson sum rule [16] to the function Performing the resummation we obtain Now we see that the integral over C 5 should converge for positive as well as for negative values of the real part of k 2 − (πn/d x ) 2 − (πm/d y ) 2 provided that Im k > 0. Therefore we have to assume t → i∞ along the contour C 5 . This leads to the choice of the contour shown in figure 5 (b).
Performing the integration we obtain Summarizing over s 1 , s 2 we get Now we can perform summations over n and m: Formulas (37), (47) and (54) give the Ewald representation of the Green function for the rectangular billiard. The integral in (37) has to be computed numerically. Now we can recapitulate the advantages of the Ewald's method. First of all, both series G (1) and G (2) are exponentially convergent. Thus we can take the analytic continuation and choose real k. Second, we have separated the part G (1) responsible for the space singularity from the part G (2) responsible for the poles information. Indeed, G (2) exponentially converges even when r = R. Third, only in G (1) there is the term corresponding to s 1 = s 2 = 0, n = m = 0 which asymptotically tends to infinity when r → R. The rest of the series is exponentially convergent. The last observation allows to compute the renormalized Green function.
Though the Poisson resummation is a common tool used to get the Ewald's representation of the Green function [13,14,15,16,18], one can avoid it and obtain formulas (47) and (54) easier (see the Appendix for details).
Using the asymptotic expansion (38) we obtain the exact Ewald representation for the renormalized Green function: with G (1) f to be computed numerically from the integral (37), and G (2) from the sum (54). Now we can compute the perturbed part of the spectrum from the condition ξ β (R, k n ) = 0, see (21), using the representation (55). This final equation in contrast to (19) lost the clearness since it depends on the as yet not defined Ewald parameter t Ew . To define it we first consider large values of k. Then to avoid exponentially large values of the function g 0 (k 2 t Ew ) as well as exponentially large amplitudes of terms with small numbers n, m in the expansion of G (2) we put t Ew = 1/k 2 . Obviously this choice is inappropriate for k → 0, since this would mean to compute a huge number of terms in G (1) . So, finally the Ewald parameter can be chosen as where k 2 0 = (π/d x ) 2 + (π/d y ) 2 is the lowest eigenvalue of the unperturbed system. To investigate the spectral statistics we can assume t Ew = 1/k 2 . Then (21) reads 1 4π [g 0 (1) Now (57) resembles (19), so the main conclusions made above could be repeated. The equation (57) is alike (3) in [7], apart from the fact that G (2) (R, R; k) is not a finite sum and the rest in (57) is not a polynomial as a function of k 2 . Equation (57) is exact and especially fits for the numerical study, since it contains exponentially convergent series. For large k the double sum in (57) can be neglected and the rest looks very similar to the "N -poles" approximation [10]. Indeed in this regime the spectrum of the billiard can be found from the equation where only a finite number of terms in the expansion of G (2) can be taken into account due to the exponential convergency. Figure 6 shows a graphical interpretation of (57). In our calculations we found that approximation (58) works perfectly above the first resonance already.

Integrated density of states
In what follows we are interested in level-spacing statistics for the subset of perturbed eigenvalues of theŠeba billiard. There are several reasons to restrict ourselves to the statistics of the subspectrum. First of all the influence of the perturbation is more pronounced if one considers only the perturbed part of the spectrum. This is probably the reason why in the pioneering work [1] only the statistics of the subspectrum is considered. Another reason to consider subspectrum's statistics is (21), which determines only the perturbed subspectrum. The graphical interpretation of (57) gives already an idea on the structure of the perturbed subspectrum (see figure 6), while considering the unperturbed subspectrum as well we loose the clearness. The last reason to consider the statistics of the perturbed subspectrum only is the direct correspondence of the perturbed subspectrum to the spectrum obtained from the reflection measurement with a single antenna introduced at the point of the perturbation. In such an experiment the unperturbed subspectrum is not seen at all since the corresponding eigenstates, vanishing at the perturbation point, can not be excited. If one considers only the perturbed subspectrum it makes a difference whether the ratios x ′ /d x and y ′ /d y are rational or irrational numbers (see figure 7). The difference arises from the fact that for irrational numbers all eigenfunctions are perturbed while for rational ones a part of eigenfunctions remains unperturbed.
To compute the statistics and compare it with GOE, semi-Poissonian and Poissonian predictions one should first unfold the spectrum to a mean level spacing of one. This can be achieved by the following definition of the scaled eigenvalues: where N(z) is a smoothed function counting a total number of eigenvalues k 2 n less then z, i. e. the integrated density of states. If a spectrum of a system is known, the function N(z) can be obtained from a numerical fit. For a conventional unperturbed two-dimensional billiard one can use the Weyl estimation of the integrated density of states (see e. g. [19]) where A 1 is the area of the billiard, A 2 is its circumference and A W is a constant. For the unperturbed rectangular billiard with the sides d x , d y we obtain A 1 = d x d y and A 2 = 2(d x + d y ). While the Weyl estimation holds for the whole spectrum of the unperturbed billiard it can not be directly applied to its subspectrum as well as to the perturbed subspectrum of theŠeba billiard. However to fit numerically the integrated density of states one can still assume that the function to be found has the form (60) with some unknown coefficients A 1 , A 2 , A W . The scaled level spacing corresponding to the nearest eigenvalues k 2 n and k 2 n+1 is Thus the constant term A W in (60) does not influence the statistics. The mean level spacing when M → ∞ in accordance with the rescaling requirement. Though the integrated density of states corresponding to the perturbed subspectrum can be fitted numerically, it is possible to estimate it a priori. Indeed in figure 6 one sees that between two successive eigenvalues of the unperturbed spectrum corresponding to poles of the function G (2) there always exists an eigenvalue of the perturbed billiard. Thus the number of perturbed eigenvalues of the billiard below z should coincide (up to a single eigenvalue) with a number of eigenvalues of the unperturbed billiard below z corresponding to nonvanishing eigenfunctions at the point of the perturbation. Obviously this result does not depend on the value of the scattering length.
Following the argument given above we can compute the expected function N e (z) right from the unperturbed billiard, where N e (z) is equal to the integrated density of those states whose eigenfunctions do not vanish at the perturbation point (x ′ , y ′ ). Let us assume that x ′ = d x p 1 /q 1 , y ′ = d y p 2 /q 2 , where p 1 /q 1 and p 2 /q 2 are irreducible fractions. From figure 7 one can draw the conclusion that eigenfunctions of the small hatched "subbilliards" with Dirichlet conditions at all boundaries are eigenfunctions of the initial billiard and vanish at the point (x ′ , y ′ ). According to the Weyl formula the number of eigenvalues below z can be estimated as for the vertical and horizontal hatched billiards respectively. Here A v and A h are some constants. We have computed twice the eigenvalues of the billiard obtained as an intersection of these subbilliards. Its number of eigenvalues can be estimated as Finally, the number of eigenvalues corresponding to vanishing eigenfunctions is Subtracting the last estimation from the total number of eigenvalues below z we obtain the following estimation for the number of eigenvalues corresponding to nonvanishing eigenfunctions: where . Surprisingly the surface contribution ∼ √ z vanishes.
In figure 8 the estimation (68) and the numerical fit of the form (60) for the integrated density of states are shown for the comparison. One can see that the estimation (68) works very well.
From (58), (68) one can compute a shift of a resonance induced by a point perturbation as compared to the mean level spacing. Indeed when k → E nm the function G (2) (R, R; k) (54) tends to the following expression: From (68) we find the mean level-spacing ∆E : Substituting (69) in (58) and using (70) we find for the relative shift of the resonance: Let us estimate the number of resonances needed to show the transition to the Poissonian level-spacing statistics. Then the relative shift should be very small for all sufficiently large numbers n and m. The sufficient condition is Depending on the values of Q and β the value of E nm can be very large.

Level-spacing statistics
In this section we are presenting a number of numerical results. Figures 9-11 show levelspacings distributions for the perturbed part of the spectrum for various situations to be discussed in detail below. For comparison the curves corresponding to Poissonian, semi-Poissonian and GOE distributions are plotted by solid, dashed, and solid lines, respectively. In each of the figures subfigure (a) shows the level-spacings distribution for the unperturbed system to make sure that the distribution is really Poissonian, since it is well-known that there may be deviations for small distances depending on the side ratio of the rectangle. Subfigures (b)-(d) show level-spacings distributions for β = 1, 0.1, 0.02 respectively. Note that a decreasing value of β means an increase of the perturbation.
In figure 9 the scatterer is placed in the center whereas in figure 10 it is at the point (0.55d x , 0.65d y ). In both cases about 1500 lowest perturbed eigenvalues have been considered. For β = 1, i.e. for a weak perturbation, the distribution shows a linear repulsion for small distances and an exponential tail, but not a semi-Poissonian behavior in the strict sense (dashed line). With β = 0.1, 0.02 there is a gradual transition to a broader distribution, resembling GOE one for the scatterer in the center (figure 9(d) ) and a semi-Poissonian distribution for the scatterer in the off-center position ( figure  10(d) ). This observation would deserve more quantitative treatment, but this goes beyond the scope of this paper. Qualitatively it may be understood from the fact that in the center of the billiard all perturbed eigenfunctions have the same value, while for an off-center position there is a distribution of eigenfunctions amplitudes giving rise to a corresponding distribution of resonances shifts. Figure 11 finally shows the level-spacings distribution again with the perturbation in the center but now for the numbers of perturbed eigenvalues from 25000 to 27000. Comparison of figures 9 and 11 shows a pronounced change of the distribution towards Poissonian with increasing eigenvalues numbers. This is particularly evident for the weaker perturbations β = 1, 0.1 ( figure 11(b)-(c) and demonstrates the main result of this paper: with increasing eigenvalues numbers eventually the level-spacings distribution of the unperturbed system is recovered.

Conclusions
Let us now conclude. First of all we have presented in the paper the complete solution of the spectral problem for the rectangular billiard with a single point perturbation. In contrast to previous studies [1,7] we have shown that the statistics of theŠeba billiard tends to a Poissonian when the number of levels taken into account tends to infinity. The estimation given at the end of Section 5 showed, however, that the transition to Poissonian statistics appears, depending on the scattering length, only at exponentially large quantum numbers. The solution is based on the Ewald representation of the renormalized Green function (55). This representation contains exponentially rapidly convergent series. Together with the Ewald representation of the usual Green function (37), (47), (54) the presented approach is a powerful tool to analyze various experiments Let us consider the initial problem (29) in the rectangular billiard with proper boundary conditions. Then the solution can be written in two equivalent forms: in the form of the images representation g i (t; r, R, k) and in the form of eigenmodes representation g e (t; r, R, k): g i (t; r, R, k) =