Scarring in Rough Rectangular Billiards

We study the mechanism of scarring of eigenstates in rectangular billiards with slightly corrugated surfaces and show that it is very different from that known in Sinai and Bunimovich billiards. We demonstrate that there are two sets of scar states. One set is related to the bouncing ball trajectories in the configuration space of the corresponding classical billiard. A second set of scar-like states emerges in the momentum space, which originated from the plane-wave states of the unperturbed flat billiard. In the case of billiards with one rough surface, the numerical data demonstrate the repulsion of eigenstates from this surface. When two horizontal rough surfaces are considered, the repulsion effect is either enhanced or canceled depending on whether the rough profiles are symmetric or antisymmetric. The effect of repulsion is quite strong and influences the structure of all eigenstates, indicating that the symmetric properties of the rough profiles are important for the problem of scattering of electromagnetic (or electron) waves through quasi-one-dimensional waveguides. Our approach is based on the reduction of the model of one particle in the billiard with corrugated surfaces to a model of two artificial particles in the billiard with flat surfaces, however, with an effective interaction between these particles. As a result, the analysis is conducted in terms of a two-particle basis, and the roughness of the billiard boundaries is absorbed by a quite complicated potential.


Introduction
Recently, there has been a sharp increase in articles discussing the occurrence of ergodicity in physical systems, both classical and quantum [1][2][3][4][5][6][7][8]. This interest is due to the practical problem of statistical description of the behavior of a closed many-particle system whose equations of motion are strictly deterministic, that is, do not contain any random parameters. For a long time, this problem was considered fundamental for the justification of statistical mechanics, and it was assumed that its solution is associated with the so-called ergodicity of the motion of the system under consideration. The classical trajectory of ergodic systems densely and uniformly covers the entire phase space of the system, as a result of which the time average value of any observable is equal to the phase space average as time tends to infinity.
As a result of numerous attempts to prove this property in general terms for real physical systems, it was realized that the proof is possible only for special systems. The most well-known of such systems are the so-called scattering billiards, in particular, Sinai and Bunimovich billiards, in which the trajectory of a particle due to multiple elastic reflections from the walls is unstable for any initial condition except for special values whose measure is equal to zero. Thus, the behavior of the system is determined by a single trajectory, which makes it possible to speak of an ergodic type of motion.
It should be noted that the instability of motion in such billiards is exponentially strong, at which the distance between two close trajectories grows, on average, exponentially fast with time. When trajectories are reflected from the walls of the billiard, the so-called mixing occurs, as a result of which the motion becomes chaotic and indistinguishable from random. Indeed, as already discussed by the Soviet physicist N. S. Krylov in his book [9], the mixing is the main mechanism for the emergence of the most characteristic property of the statistical behavior, namely, the relaxation of a system to statistical equilibrium.
Conversely, examples can be given when a given system is ergodic but its movement is regular, meaning that the relaxation of the system occurs at infinite times, which does not correspond to physical reality.
For quantum systems, the situation is much more complicated, since there is no rigorous and generally accepted definition of quantum ergodicity. The only reliable situation when one can speak of quantum ergodicity is for the systems that are fully ergodic in the classical limit. Needless to say, there are just a few such systems. As an example of the "true" quantum ergodicity, one can consider completely random matrices of large size for which the theory is well developed. For such matrices, the components of the eigenfunctions are Gaussian distributed, which is a property that is a consequence of the uniform distribution of the eigenfunction vector over the surface of the N-dimensional sphere with N 1. It can be seen that in this case, the eigenfunctions are not only ergodic but also maximally chaotic (random), see e.g., [10].
However, for realistic physical systems, many-body matrix elements in a physically chosen basis, as a rule, are not completely random and fill only a band of finite width (determined by a finite radius of interaction in the energy space). Therefore, one can speak of ergodic eigenfunctions occupying only some energy shell of finite size. Clearly, determining the shape and size of this shell is extremely difficult [11]. In any case, at the moment, there is no unified theory of quantum ergodicity of isolated systems which takes into account the influence of the finiteness of the energy shell.
Returning to systems with a small number of particles (small number of degrees of freedom), for example chaotic billiards for which classical ergodicity is rigorously proven, after several years of their study, it was found that many of the eigenfunctions are nonergodic even for very high energies [12,13]. In principle, such non-ergodic states do not contradict Shnirel'man's theorem [14], according to which their number must decrease with increasing energy. Numerical experiments with ergodic billiards have shown that in many cases, these "scarring" states in the configuration space can be associated with unstable periodic orbits, the number of which increases with energy. However, since these trajectories are isolated, when averaged over infinite time, they do not contribute to the mean values of the observables.
To date, there have been many papers that investigate the properties of many-particle scarring states in various physical systems (see, for example, [15][16][17][18][19][20][21] and references therein). The main interest is related to the mechanism for the emergence of such states immersed in a set of ergodic states. In the case of quantum systems that do not have a corresponding classical analogue, this issue becomes extremely complicated. It was found that the main characteristic of such scars in the quantum description is a small number of components of the corresponding eigenstates treated on some physically justified basis in comparison with the large number of states that can be considered as ergodic. Thus, scars can be defined as the localized states, associated with unstable periodic orbits when applicable, and embedded into the set of the ergodic states. As a result of numerous studies of quantum many-body systems, it became clear that the typical mechanism for scars is the presence of local symmetries of the system, which can correspond to the presence of local integrals of motion [20,21].
Our interest in this article is to study the mechanism for the emergence of scar states in rectangular billiards with rough horizontal surfaces. In these billiards, it is relatively easy to trace their emergence since such systems can be considered close to integrable, despite the fact that they are seem to be rigorously ergodic even for any weak boundary roughness. Such billiards of finite size were under close investigation in view of surface scattering in quasi-one-dimensional waveguides. The conventional theory of scattering in the waveguides with a weak roughness was developed long ago; however, recently, researchers found a new mechanism of surface scattering which is due to the so-called square-gradient roughness [22][23][24][25]. This mechanism was neglected in previous studies; however, it should be taken into account when the correlation length computed along the scattering profiles is very small. In particular, it was analytically shown that the scattering properties strongly depend on the correlations between upper and lower horizontal profiles. Specifically, the scattering is very different whether the two random profiles are symmetric or antisymmetric. Our study also shows a quite strong dependence of the degree of localization of the scar states on the type of correlations between the scattering profiles.
In Section 2, we describe the model under study, giving exact expressions for the matrix elements of an effective Hamiltonian in close correspondence with the model. In Section 3, we explain how we measure the degree of localization of the eigenstates in the energy basis corresponding to non-interacting artificial particles and therefore to the billiard with flat horizontal surfaces. In this representation, one easily identifies one set of scar states related to classical trajectories which are parallel to the rough boundaries of the billiard.
We also show what these eigenstates look like in the configuration representation for one rough surface. In addition, we demonstrate that there is another set of scar states related to bouncing ball classical trajectories which are perpendicular to the rough boundaries of the billiard. These states consist of many "unperturbed eigenstates", which is in contrast with those from the first set. In Section 5, we analytically show how the number of strongly localized states from the first set, which corresponds to classical plane waves with the longest perpendicular wavelength, decreases with energy. In Section 6, we study billiards with two horizontal rough profiles and show how the structure of localized eigenstates depends on whether the profiles are symmetric or antisymmetric. The main attention here is paid to the effect of the repulsion of eigenstates from the corrugated surfaces. In Section 7, we draw our conclusions.

Model of Rough Billiards
We consider billiards that are periodic in the longitudinal coordinate x, with Dirichlet boundary conditions on the upper f 1 and lower f 2 surfaces, Here, L y is the average width of the billiard and ξ 1,2 (x + L x ) = ξ 1,2 (x) with ξ 1,2 (x) = 0. The angular brackets stand for the average over one period L x or, in the case of a random profile, over different realizations of ξ 1,2 (x). This model has been thoroughly studied in refs. [26][27][28][29][30][31][32][33][34][35] for the specific case ξ 1 (x) = cos(2πx/L x ) and ξ 2 (x) = 0, the so-called "cosine" or "ripple" billiard [36]. The main interest was in the properties of the energy spectrum [26,28] and in the quantumclassical correspondence for the shape of eigenstates (SE) and local density of states (LDOS) [29][30][31][32]. It was shown that for highly excited states, the global properties of the SES and LDOS in the quantum regime are similar to those described by the equations of motion for a classical particle moving inside the billiard. On the other hand, quite strong quantum effects have been revealed for individual eigenstates in a deep semiclassical region [31,32].
In this paper, we address the case of billiards with rough surfaces, focusing on the structure of eigenstates, see also [22][23][24][25][37][38][39]. The surfaces are modeled by a large sum of harmonics A k drawn from a flat random distribution defined in the interval [−A, A] (with A such that |ξ 1,2 (x)| ≤ 1). With the increase of N T , the degree of complexity of ξ 1,2 (x) also increases, and for N T 1, the surfaces can be treated as random: see an example in Figure 1. In what follows, we focus on rough billiards with N T = 100. For simplicity, we choose ξ 1 (x) = ξ 2 (x). In the first part of our paper, we review the case of billiards with one flat boundary; that is W 2 = 0, see, e.g., Figure 1. Below, in Section 6, we consider the billiards with two antisymmetric, W 1 = W 2 , and symmetric surfaces, Originally, the model is described by the Hamiltonian for a free particle inside the billiard with rough boundaries (such that the corresponding wave function Ψ(x, y) obeys Dirichlet boundary conditions, Ψ(x, y) = 0 at y = f 1 (x) and y = f 2 (x)). However, in order to solve it numerically, it is useful to make a canonical transformation to new variables in which the new Hamiltonian incorporates surfacescattering effects into an effective interaction potential between artificial particles identified with the two new degrees of freedom. This can be achieved by the transformation to new canonical coordinates, As a result, the boundary conditions for the new wave function are trivial: Ψ(u, v) = 0 at v = 0 and v = 1. The Schrödinger equation in new coordinates can be obtained from the covariant expression for a particle moving (in the absence of potentials) in a Riemannian curved space [40], where the quantum Hamiltonian in covariant form [41] iŝ and the covariant momenta arê Here, α, β = u, v, the operator ∆ cov is the Laplace-Beltrami operator, g is the metric, and g αβ is the metric tensor. The wave functions Ψ(u, v) are normalized as Even though (5) is still the kinetic energy, the resulting differential equation takes a much more complicated form than the ordinal Laplacian. This is the price to pay when transferring the effect of the boundaries onto the operator. Then, by substituting in (6) the explicit expressions for the metric tensor where f 1,2 ≡ ∂ f 1,2 (u)/∂u, and the metric we obtainĤ Note that the rough boundaries f 1,2 and their derivatives f 1,2 are fully incorporated into the Hamiltonian operator. We remark that this representation allows us to treat the original model of one free particle in the rough billiard as a model of two interacting "particles" identified with the two degrees of freedom u and v, where the Hamiltonian (10) is separated aŝ Here,P u = −ih[∂ u + (1/4)∂ u ln(g)],P v = −ih∂ v , andV stands for an effective interaction potential between the "particles" u and v. In the following, we treatĤ 0 as the Hamiltonian of two non-interacting particles. The eigenstates ofĤ 0 define the unperturbed basis in which the eigenstates of the total HamiltonianĤ(u, v,P u ,P v ) are expanded. Such a representation turns out to be convenient for the study of the chaotic properties of the model, since one can use the tools and concepts developed in the theory of interacting particles (see, for example, [42,43]). Now, we expand the αth eigenstate of energy E α as where φ mn (u, v) = u, v | m, n are the eigenstates of the unperturbed HamiltonianĤ 0 (u, v), and the indexes n and m are the quantum numbers corresponding to the u and v coordinates, respectively. The function φ mn (u, v) has to satisfy the boundary conditions of the problem in order to form a Galerkin basis [44]. The eigenstates ofĤ 0 are with eigenvalues Here, formally, −∞ < n < ∞ and 1 ≤ m < ∞, but in practice, for the numerical computations, the unperturbed basis is truncated such that −N max < n < N max and 1 ≤ m < M max . The factor π −1/2 g −1/4 in (14) arises from the orthonormality condition (8). Then, we can solve the eigenvalue problem by diagonalizing the Hamiltonian (10) on the basis φ mn (u, v): where the matrix elements are Since the Hamiltonian (10) is periodic in u, its eigenstates are Bloch states. This allows us to write the solution of the Schrödinger equation in the form The statistical properties of eigenstates do not depend on a specific value of the Bloch index k inside the band, except at k = 0 and k = ±π/L x , so we avoid these values of k in our numerical calculations.
Finally after some algebra, we obtain where Notice that depending on the billiard geometry, some of the integrals above may vanish. For example, in the case W 2 = 0 (i.e., for the billiard of Figure 1), we have J 1 = J 4 = J 7 = 0.

Structure of the Hamiltonian Matrix
In order to study the structure of eigenstates of the total HamiltonianĤ(u, v) one needs, first, to choose a way of ordering the unperturbed basis in which we represent the Hamiltonian matrix H l,l (k) = l |Ĥ | l k . Specifically, we have to assign an index l, labeling the basis state | l k ≡| m, n k to each pair of indecies (m, n) (note that although the energy spectra are independent of the assignment (m, n) → l, the structure of the eigenstates is not). In addition, note that the size of the Hamiltonian matrix is determined by the maximum values of n and m: −N max ≤ n ≤ N max and 1 ≤ m ≤ M max .
A natural assignment is the following one [30]. Let us fix the lowest value of n (−N max ) and sweep all values of m (1 ≤ m ≤ M max ). This gives l = 1, 2, . . . , M max . Then, do the same for n = −N max + 1, which gives M max + 1 ≤ l ≤ 2M max , and so on, until finally, we have 1 ≤ l ≤ L max , where L max = (2N max + 1)M max defines the Hamiltonian matrix size, L max × L max . This rule results in a block structure of the Hamiltonian matrix with a block size equal to M max . Figure 2 (left) shows the lower part of a 4030 × 4030 matrix with (N max , M max ) = (32, 62). Here, we can see a number of blocks of size 62 × 62 corresponding to n, n = −32, . . . , −17. In this representation, the matrix is clearly a block matrix.
The above way of ordering the unperturbed basis corresponds to the "channel representation" (or momentum representation), since the index m labels a specific transverse channel for the propagation of the wave through the billiard; see Equation (14). However, for our purposes, it is essential to use the "energy representation", according to which the unperturbed basis is ordered in increasing energy, E 0 . This defines a new rule l → l new = l new (n, m). In Figure 2 (right), we show the lower part of the Hamiltonian matrix of Figure 2 (left) but now in the "energy representation". Notice that in contrast to the block structure of the Hamiltonian matrix in the "channel representation", in the "energy representation", the Hamiltonian matrix shows a band-like structure. In order to analyze the structure of the eigenstates ofĤ(u, v) in detail, we diagonalize Hamiltonian matrices in the "energy representation" and construct the "state matrices" |C α l new | 2 . In Figure 3, we show the lower part of the state matrix corresponding to the Hamiltonian matrix of Figure 2 (right). Here, C α l new are the amplitudes of the eigenstates in the (energy-ordered) basis representation given by the index l new . Namely, the index l new refers to unperturbed basis states that correspond to the unperturbed HamiltonianĤ 0 . The index α refers to a specific exact eigenstate. All eigenstates are ordered in increasing energy, with α = 1 the ground state. Therefore, to understand how strongly localized/extended the exact eigenstates in the unperturbed basis are, one should fix the value of α and explore the dependence of |C α l new | 2 on l new . A crucial point in our study is that the eigenstates of the Hamiltonian in the "energy representation" have a very convenient form for the analysis. The advantage of the "energy representation" over the "channel representation" (i.e., when the Hamiltonian matrix has a block structure) is clearly seen in Figure 4 where an arbitrarily chosen eigenstate is given in the two representations. One can see that in the "channel representation", the eigenstate has a kind of regular and extended structure, while in the "energy representation", the eigenstate is compressed. In the latter case, one may use a statistical approach to describe the global properties of such eigenstates; see refs. [43,45]. Specifically, these eigenstates can be characterized by introducing an envelope around which the components are expected to fluctuate in a pseudo-random way. We stress that by using this energy ordering, it is possible to relate the global form of eigenstates in the energy representation with its classical counterpart; see, e.g., refs. [29][30][31][32].

Eigenstates in Energy Representation
In Figure 5, we present two typical pairs of consecutive eigenstates (α = 390 and 391 and α = 407 and 408). The difference between the eigenstates on the left panels (α = 390 and 407) and the other two is clearly qualitative. More specifically, while the states α = 390 and 407 are extended (in energy) eigenstates, constituted by practically all basis states within the shown energy range, the eigenstates α = 391 and 408 are mostly unperturbed: they are extremely localized in energy. Indeed, by neglecting all small amplitude components surrounding the main component (see left panels of Figure 5), we can determine the basis state l new , defined by the pair (m, n), that most closely resembles the exact eigenstates. We find that this always corresponds to the lowest values of the transversal mode m. This fact can be understood by the following physical argument. Consider an eigenstate of the flat billiard φ m,n (X, Y, k = 0) ∝ sin(mπY/L y ) exp(iK x X) with energy E 0 = (h 2 /2m e )(K 2 x + K 2 y ), where K y = mπ/L y = 2π/Λ y . Turning on the perturbation (flat to rough billiard) will affect the high energy unperturbed states differently depending mainly on the value of Λ y . For example, for m = 1, the ratio Λ y /W 1 is 2L y /W 1 ≈ 33 (with W 1 = 0.06L y , the value of W 1 we use throughout this work), which is so large that the state cannot "see" the roughness and thus will remain essentially unperturbed. In contrast, for unperturbed states with the same (or about the same) energy but with large values of m (say, m = 62 = M max and correspondingly small K x ), their Λ y is sufficiently small compared to the amplitude of the roughness (Λ y /W 1 ≈ 0.5), so that the rough boundary produces a strong mixing of unperturbed levels. The resulting exact eigenstate will consist of many components extended over the energy. Note that we may treat strongly localized states in the channel representation as a kind of scar state. The fact that the extremely localized (in energy) eigenstates can be identified with the plane waves φ m,n (x, y) with small m, proper of the flat billiard, makes us expect that they will also be similar to plane waves when presented in configuration representation |Ψ α (x, y)| 2 . Then, in Figure 6 (right), the two localized (in energy) eigenstates of Figure 5 (right) are shown in configuration representation. It is quite unexpected that these eigenstates are very different from the unperturbed ones, even though they are similar in energy representation. Figure 6 (right) shows that the rough boundary "pushes" the probability |Ψ α (x, y)| 2 away from it. Thus, the eigenstate α = 391, whose main component in energy representation is identified with m = 1, differs importantly from the unperturbed mode with m = 1 whose maximum is at the billiard center y = L y /2 (as in the case of the lowest eigenmode of a box of width L y with hard walls). Similar repulsion occurs for the eigenstate α = 408 identified with m = 2. We stress that this repulsion effect occurs only for rough billiards and is stronger the more rough the billiard boundary is, see, e.g., ref. [39]. A detailed analysis of this repulsion effect will be performed in Section 6. Finally, for comparison purposes, in Figure 6 (left), we also present the two extended (in energy) eigenstates of Figure 5 (left).
In order to characterize quantitatively the eigenstates, we compute various localization measures. The first one is the so-called entropy localization length l H , Here, H stands for the Shannon entropy of an eigenstate in a given basis, and H GOE is the entropy of a completely chaotic state which is characterized by Gaussian fluctuations (for N → ∞) of all components C α l new with the same variance w α l new = 1/N, where w α l new = |C α l new | 2 . The latter property occurs for completely random matrices belong-ing to a Gaussian Orthogonal Ensemble (GOE). Defined in this way, the quantity l H gives the measure of the effective number of components in an eigenstate. For example, the eigenstates of Figure   The second quantity, l ipr , which gives another measure of the effective number of components in an eigenstate, is expressed via the inverse participation ratio P, where P GOE ≈ 3 is chosen in order to obtain l ipr = N in the GOE limit case. Correspondingly, the eigenstates of Figure 5 have l ipr (α = 390) = 1061.2, l ipr (α = 391) = 6.32, l ipr (α = 407) = 1464.4, and l ipr (α = 408) = 14.99, where a high correlation with l H can be seen. The above two definitions of localization lengths are the most frequently used when describing the global structure of eigenstates. One should note that these quantities provide an estimate of the effective number of large components, independently on the location of these components, in the unperturbed basis.
To obtain a complete panorama, in Figure 7a,b we plot these two measures, l H and l ipr , for the eigenstates |α of the rough billiard. The strong fluctuations of the localization measures are evident in these figures. We can see that neighboring high-energy eigenstates may have drastically different localization measures, which is in agreement with the discussion above about the existence of localized and extended eigenstates. Moreover, these figures give us information about the relative number of each type (localized and extended eigenstates) to be found in a given energy range.  Additional information about the structure of eigenstates can be obtained from the "width" or mean square root l σ of an eigenstate, which is computed as where n c = ∑ l new l new w α l new determines the centroid of an eigenstate in the unperturbed basis.
Comparison of the width l σ with l H and l ipr gives the possibility of detecting the so-called sparsity of eigenstates. Indeed, small values of the ratio l H /l σ (or l ipr /l σ ) indicate that there are "holes" in the structure of the eigenstates; therefore, such eigenstates are sparse; see Figure 7c. A detailed analysis shows the existence of sparsed eigenstates. As for the centroid, n c > α observed in Figure 7d indicates that the interaction strengthV is relatively strong compared to the unperturbed partĤ 0 .
The data of Figure 7a,b show the existence of a wide range of values of both localization measures even at high energies; moreover, there are eigenstates with small localization lengths (visible as clusters of points in the lower part of the plots) along the entire energy range. This situation is also present for the same quantities computed for the local density of states (LDOS): that is, for the basis eigenstates | l new expanded in the exact basis | α ; see Figure 8. Note that to compute the localizations lengths, as well as the centroid, for the LDOS, one must interchange α ↔ l new in expressions (19), (21) and (22). It is relevant to stress that the localization lengths for the LDOS show well-defined patterns, see Figure 8, which are not present for the exact eigenstates, see Figure 7. Inspection of  This structure is expected to prevail at all energy ranges, since in any sufficiently large range of energies, there are unperturbed states with all values of m ∈ [1, M max ]. Even deep in the semiclassical regime, extremely localized and sparse states will appear but less and less frequently since the energy differences between states of the same type increases with energy.

Scar-Like States in Energy Space
In this section, we discuss the origin of the scar-like states emerging in the energy representation; they seem to be generic for plane billiards with rough surfaces. The starting point is the observation that for a weak roughness, these states can be associated with those unperturbed states that have a small value of the index m, especially for m = 1. Since the localized eigenstates identified above with m = 1 exist at all energies, the important question is about the fraction of these states relative to the total number of states N(E) ≡ ∑ i {i|E i ≤ E} as a function of energy. This question is important in view of Shnirel'man's theorem [14] stating that the eigenfunctions of a classically ergodic system are equidistributed over the energy shell in the classical limit. In other words, we expect that where N m=1 (E) is the number of m = 1 eigenstates up to some energy E: Neither Shnirel'man's theorem nor Equation (23) says anything about how the limit is approached. Moreover, Shnirel'man's theorem does not consider the possible existence of parabolic fixed points. These occur in our billiard; they are the bouncing ball orbits or the continuous set of all horizontal trajectories (P y = 0) which do not hit the two boundaries. Quantum mechanically, there are no P y = 0 states, but the eigenstates of the flat channel with a minimum value of P y are the m = 1 states. Thus, we have referred to the m = 1 eigenstates as "bouncing ball states" in [38]. For W 1 /L y 1, one can use the unperturbed spectra given by Equation (15) to obtain a good estimate for N m=1 (E) and N(E). For large N max and M max , Equation (15) represents, ignoring the k 2 term, half of the ellipse n,m (0). In Figure 9, we show the n − m plane for the case of N max = 32 and M max = 62; there, the dots represent basis states and the red ellipse is Equation (25) with E = 400(h 2 /2m e ) and L x = Ly = 2π. Then, for large E, the number N(E) equals the area of the half ellipse, and the number N m=1 (E) of localized eigenstates equals twice the size or its minor axis: therefore, where σ = L x L y is the area of one period of the billiard. We see that Equation (26) is precisely the first term in the Weyl series for the integrated density of states (see for example [46], Section 7). Equation (28) indicates that (i) for fixed E and σ, the portion of m = 1 states, relative to the total number of states, is larger the narrower the billiard is and (ii) the convergence of N m=1 (E)/N(E) to zero is rather slow; it decreases as a power-law of E.
It is important to stress that in the case of a rough billiard having a highly modulated boundary composed of many harmonics, the expression for N(E) is not simply given by (26), since now, the perimeter γ of the boundary, as well as its curvature, contribute significantly to N(E) as [47] The coefficients C 2r+1 , r ≥ 1, in the sum of the Weyl series (29), depend on the curvature (and its derivative) of the modulated boundary (see Table 3 of [47]). As is shown in [47], the higher the index r of the Weyl coefficient C 2r+1 , the more complex the expansions for C r are, involving higher and higher powers of the curvature and its derivatives. Clearly, the larger the number of harmonics N T , the more important the perimeter and curvature terms become.
In addition, it is also instructive to look at the location of the eigenstates on the n − m plane. Thus, in Figure 10, we present the eigenstates of Figure 5 also in energy representation but now on the n − m plane. From this figure, we can observe that: (i) the components of the energy-extended eigenstates of Figure 5 (left) are not equidistributed on the n − m plane; instead, their main components are concentrated around the region with n ∼ 0 and large m. That the components are concentrated around n ∼ 0 supports the fact that these eigenstates are not extended over the x-axis when plotted in configuration representation; see Figure 6 (left). (ii) The components of the energy-localized eigenstates of Figure 5 (right) are also localized on the n − m plane; they are indeed characterized by a single value of n and several but small values of m (a detailed decoding of this m dependence will be performed in the next section).

Scars and Billiard Symmetries; Repulsion Effect
In this section, we study in detail the repulsion, in configuration representation, that suffer the energy-localized eigenstates of rough billiards from the rough boundaries. This effect, already shown in Figure 6right panels, was first reported in [38,39] but also observed in [24,48,49]. Moreover, in addition to the one-flat-boundary billiard of Figure 11a, we now consider two additional billiard geometries: the antisymmetric billiard, Figure 11b, and the symmetric billiard, Figure 11c; see also [24,25].
We start by presenting in Figure 12 few typical extremely energy-localized eigenstates for the billiards of Figure 11 in three different energy regions. The main component of all those eigenstates can be identified with an unperturbed state with m = 1 (see also Figure 13) where we present again the eigenstates of Figure 12 but now in configuration representation). Close inspection of Figure 12 shows that for the one-flat-boundary billiard and the symmetric billiard, panels (a) and (c) of Figure 12, respectively, the localized eigenstates have several components of appreciable magnitude. This is in contrast with the eigenstates of the antisymmetric billiard, see panel (b) of Figure 12, which has only one effective component with magnitude very close to unity.
The decoding of the main components of the energy-localized eigenstates of the oneflat-boundary billiard, corresponding to panels (a) of Figure 12, shows that the second most important component of each of the localized states corresponds to m = 2, the next one to m = 3 and so on, all with the same value of n. That is, there is no mixing between different values of n. Moreover, the amplitudes of these eigenstate components decay exponentially as a function of m; see Figure 14 (left).
On the other hand, by decoding the main components of the energy-localized eigenstates of the symmetric billiard, corresponding to panels (c) of Figure 12, we see that the second most important component of each of the localized states corresponds to m = 3, the next one to m = 5 and so on, all with the same value of n but with alternating signs. Again, there is no mixing between different values of n. In addition, the amplitudes of these eigenstate components show an exponential decay as a function of m, which is in this case modulated by a sinus function; see Figure 14 (right).  Figure 11. The rough billiards analyzed below: (a) a billiard with one flat boundary, W 2 = 0; (b) an antisymmetric billiard, W 1 = W 2 ; and (c) a symmetric billiard W 1 = −W 2 . N T = 100, W 1 /L y = 0.06, L x /L y = 1, L x = 2π, and ξ 1 (x) = ξ 2 (x).
The parameters S a = 0.8, β a = 0.53, S c = 0.9, and β c = 0.426 are obtained by fitting the data to the dependence given by the expressions for C α mn above; see the fittings in Figure 14. It is important to remark that for a fixed N T , all energy-localized states are characterized by the same values of S a,b and β a,b . Substituting the dependence of C α mn into Equation (13) gives where C a,b,c arise to satisfy the orthonormality condition in curvilinear coordinates, Then, it follows that These latter expressions give the average shape of the localized eigenstates in the configuration representation, projected onto the y− | Ψ α | 2 plane. In addition, note that Ψ α loc (y) does not depend on the parameters S a,c , so the relevant parameters are β a,b . To obtain | Ψ α loc (y) | 2 above, we have made the approximation v ≈ y and g ≈ 1 since W 1 /L y = 0.06 1. Finally, after some algebra, we obtain Note that when β a,b → ∞, i.e., when the energy-localized eigenstates have only one single component, | Ψ α loc (y) | 2 [a] =| Ψ α loc (y) | 2 [c] = (1/π) sin πy/L y , as required. In Figure 15, we plot Equations (34)-(36) together with the numerical data from the eigenstates of Figure 13. From this figure, it is clear that (i) in the one-flat-boundary billiard, the energy-localized eigenstates are repelled from the rough boundary toward the flat boundary; see the curves labeled with (a), as already observed in Figure 6; (ii) in the symmetric billiard, the energy-localized eigenstates are repelled from both rough boundaries toward the billiard center; see the curves labeled with (c); while (iii) in the antisymmetric billiard, the repulsion effect is completely absent, so the eigenstates correspond to those of the flat billiard; see the curves labeled with (b). A detailed study [24] revealed that the absence of repulsion in the antisymmetric billiard is due to the fact that in this billiard, there is no square-gradient scattering between different modes or channels, in contrast to the one-flat-boundary billiard and the symmetric billiard. Thus, one can conclude that the main contribution to the repulsion is due to the intermode square-gradient scattering terms of the Hamiltonian matrix (i.e., the integrals J 3 , J 6 , and J 7 in Equation (17), which vanish for the antisymmetric billiard). This conclusion supports the observation made above according to which the effect of repulsion can be explained as due to a strong localization in the channel space. This localization occurs due to a relatively strong interaction between different conducting channels (or billiard modes).  Note that since the parameters S a,b and β a,b are the same for all m = 1 eigenstates for a fixed N T , one can find S a,b and β a,b for some low-energy eigenstate (so that the matrix to be diagonalized is small) and with this one can infer the global shape of all m = 1 eigenstates. Conversely, knowing the shape of an m = 1 eigenstate within any energy range, one may assess by looking at the amount of repulsion the degree of complexity, presumably unknown, of the boundaries.
The excellent agreement between Equations (34)- (36) and the numerical data from the eigenstates of Figure 13 indicates that for extremely localized eigenstates with m = 1, there is no need to diagonalize the whole matrix since, there is no mixing between different values of n. In order to obtain a good approximation of these eigenstates, one only needs to diagonalize the block of size M max × M max corresponding to the quantum number n α , which is a good quantum number in the presence of perturbation.
Within this approach, it is also possible to study energy-localized eigenstates with m > 1 using blocks of the Hamiltonian (17) of size M max × M max . So, in Figure 16, we plot the 2nd, 3rd, and 4th lowest eigenstates in configuration representation of a block of the Hamiltonian (17) of size 100 × 100 for the one-flat-boundary billiard (left panels), an antisymmetric billiard (central panels), and the symmetric billiard (right panels). They correspond to eigenstates characterized by m = 2, 3, and 4, respectively. As expected, they have 2, 3, and 4 maxima in the transverse direction, respectively. As well as for the eigenstates characterized by m = 1, see Figures 13 and 15, the eigenstates of Figure 16 suffer repulsion from the rough boundaries in the one-flat-boundary billiard and the symmetric billiard, while in the antisymmetric billiard, the repulsion is absent.  Figure 11b,c, respectively. N T = 100, W 1 /L y = 0.06, L x /L y = 1, L x = 2π, ξ 1 (x) = ξ 2 (x), M max = 100, n α = 0, and k = 0 were used.

Conclusions
In this paper, we have studied the properties of rectangular billiards with one and two rough boundaries. In the classical description, such billiards can be considered as completely chaotic, although this has not been rigorously proven. Our interest was to understand the typical structure of the eigenfunctions in the quantum description and to discover the mechanism of occurrence of scar states. For numerical simulation, we have chosen the method of reducing the original model with one particle moving inside the billiard to a model of two artificial particles in a billiard with flat boundaries, however with an effective interaction between the particles. Thus, the complexity of the rough boundary is embedded into a quite complicated interaction potential. In this representation, the Hamiltonian of the system can naturally be represented as the sum of two terms, one of which corresponds to the motion of two particles in a billiard with flat walls, and the second term describes the interaction between particles. As can be seen, such a scheme is similar to that which is often used in physics when describing isolated systems with many interacting particles. Thus, the question arises about the correspondence between the properties of eigenfunctions in a many-particle basis (in our case, a two-particle basis) and the properties of eigenfunctions in the configuration space. As was found, this correspondence is not simple, especially if we are interested in the degree of localization of the eigenfunctions and, in particular, in the occurrence of scar states.
As a result of our study, we came to the conclusion that one can roughly speak of three different types of eigenstates. The first type includes states whose density in the configuration basis is concentrated in the vicinity of bouncing ball trajectories that are perpendicular to two horizontal surfaces. Such highly localized functions (with some degree of chaos) are analogous to well-studied scar states in stadium billiards. The second type of eigenfunctions, which are strongly localized in the momentum space, can be associated with the plane waves in the billiard with flat boundaries. To some extent, they can also be termed as scar states because of their strong correspondence to unperturbed eigenstates. There is a third type of eigenfunctions that consists of many components in both the energy space of the Hamiltonian and in the configuration space. These eigenstates are quite complicated and can be treated as partially chaotic. Thus, it can be assumed that these eigenstates quickly become completely chaotic with increasing energy in comparison with those functions that are strongly localized either in energy or in configuration space.
We have considered in detail how the number of eigenstates strongly localized in the energy space decreases with increasing energy. This question is related to Shnirel'man's theorem, according to which, in the classical limit, the measure (number) of scar functions decreases to zero. It is worth noting that this theorem does not predict at what rate (in energy) all eigenfunctions become ergodic. In the general case of ergodic billiards, this question remains open, but in our case, the estimate indicates that the number of scar states (in the energy space) decreases rather slowly, namely, inversely proportional to the square of the energy.
We were also interested in the "repulsion" effect, discovered in [22], according to which for a billiard with one corrugated surface, some of the eigenstates suffer a shift of the maximal density toward the surface, which is flat. Our numerical data showed that this effect is highly pronounced for those eigenfunctions that are close to plane waves traveling in the horizontal direction. For a subset of such eigenstates, we have derived an approximate expression that clearly indicates an exponential localization in the momentum space. Even more interesting was to compare the structure of this type of eigenstate for the billiard with two corrugated surfaces. We have found that for symmetric rough surfaces, the repulsion effect is strongly enhanced. On the other hand, for antisymmetric profiles, the repulsion is absent. The origin of this phenomenon can be explained by analyzing the structure of the Hamiltonian matrix, which is strongly influenced by the form of the off-diagonal matrix elements that depend on the amplitude, first derivative, and the second derivative of the billiard profiles. Apart from that, there are terms that describe the intercorrelations between profiles. Our analysis explains both the repulsion enhancement and the disappearance of repulsion. Due to this analysis, one can predict the impact of the type of symmetry between the profiles on the scattering properties of quasi-one-dimensional waveguides with corrugated surfaces.
This article is dedicated to Professor Giulio Casati on the occasion of his 80th birthday. We wish him to be healthy and scientifically productive for the next 20 years.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.