# Quantization of a billiard model for interacting particles

## Abstract

We consider a billiard model of a self-bound, interacting three-body system in two spatial dimensions. Numerical studies show that the classical dynamics is chaotic. The corresponding quantum system displays spectral fluctuations that exhibit small deviations from random matrix theory predictions. These can be understood in terms of scarring caused by a 1-parameter family of orbits inside the collinear manifold.

Billiards are interesting and useful models to study the quantum mechanics of classically chaotic systems . In particular, the study of the completely chaotic Sinai billiard  and Bunimovich stadium  showed that the quantum spectra and wave functions of classically chaotic system exhibit universal properties (e.g. spectral fluctuations)  as well as deviations (e.g. scars of periodic orbits)  when compared to random matrix theory (RMT) predictions of the Gaussian orthogonal ensemble (GOE) .

While the theory of wave function scarring has reached a mature state in two dimensional systems  there is a richer structure in more than two dimensions. In particular, invariant manifolds in billiards  and systems of identical particles  may lead to an enhancement in the amplitude of wave functions  provided classical motion is not too unstable in their vicinities.

The purpose of this letter is twofold. First we want to quantize a billiard model for three interacting particles and study a new type of wave function scarring. Second, our numerical results strongly suggest that the system under consideration is chaotic and ergodic. This is interesting in view of recent efforts to construct higher-dimensional chaotic billiards .

This letter is organized as follows. First we introduce a billiard model of an interacting three-body system and study its classical dynamics. Second we compute highly excited eigenstates of the corresponding quantum system and compare the results with RMT predictions.

Recently, a self-bound many-body system realized as a billiard has been studied in the framework of nuclear physics . We want to consider the corresponding three-body system with the Hamiltonian

$$H=\underset{i=1}{\overset{3}{}}\frac{\stackrel{}{p}_i^2}{2m}+\underset{i<j}{}V(|\stackrel{}{r}_i\stackrel{}{r}_j|),$$
(1)   
where $`\stackrel{}{r}_i`$ is a two–dimensional position vector of the $`i`$-th particle and $`\stackrel{}{p}_i`$ is its conjugate momentum. The two-body potential is

$`V(r)=\{\begin{array}{cc}0\hfill & \text{for }r<a,\hfill \\ \mathrm{}\hfill & \text{for }ra.\hfill \end{array}`$ (4)   
The particles thus move freely within a convex billiard in six-dimensional configuration space and undergo elastic reflections at the walls. Besides the energy $`E`$, the total momentum $`\stackrel{}{P}`$ and angular momentum $`L`$ are conserved quantities which leaves us with three degrees of freedom. In what follows we consider the case $`\stackrel{}{P}=0,L=0`$.

To study the classical dynamics it is convenient to fix the velocity $`\stackrel{}{v}_1^2+\stackrel{}{v}_2^2+\stackrel{}{v}_3^2=1`$ and perform computations with the full Hamiltonian (1) without transforming to the subspace $`\stackrel{}{P}=0,L=0`$. We want to compute the Lyapunov exponents of several trajectories. To this purpose we draw initial conditions at random and compute the tangent map  while following their time evolution. To ensure good statistics and good convergence we follow an ensemble of $`7\times 10^4`$ trajectories for $`10^5`$ bounces off the boundary. All followed trajectories have positive Lyapunov exponents. The ensemble averaged value of the maximal Lyapunov exponent and its RMS deviation are $`\lambda a=0.3933\pm 0.0015`$, while the second Lyapunov exponent is also always positive. Thus, the system is chaotic for practical purposes. However, we have no general proof that no stable orbits exist. The reliability of the numerical computation was checked by (i) comparing forward with backward evolution, (ii) observing that energy, total momentum and angular momentum are conserved to high accuracy during the evolution and (iii) using an alternative method  to determine the Lyapunov exponent.

The considered billiard possesses two low-dimensional invariant manifolds that correspond to symmetry planes. The first “collinear” manifold is defined by configurations where all three particles move on a line. The dynamics inside this manifold is governed by the one-dimensional analogon of Hamiltonian (1). After separation of the center-of-mass motion one obtains a two-dimensional billiard with the shape of a regular hexagon. This system is known to be pseudo-integrable . To study the motion in the vicinity of the collinear manifold we compute the full phase space stability matrix for several periodic orbits inside the collinear manifold which come in 1-dim families and can be systematically enumerated using the tiling property of the hexagon. All considered types of orbits except two are unstable in the transverse direction: (i) The family of bouncing ball orbits (i.e. two particles bouncing, the third one at rest in between) is marginally stable (parabolic) in full phase space. (ii) The family of orbits depicted in Fig. 1 is stable (elliptic) in two transversal directions and marginally stable (parabolic) in the other 10 directions of 12-dim phase space. Though this behavior does not spoil the ergodicity of the billiard one may expect that it causes the quantum system to display deviations from RMT predictions. Note that this family of periodic orbits differs from the bouncing ball orbits which have been extensively studied in two- and three-dimensional billiards  since (i) it is restricted to a lower dimensional invariant manifold, and (ii) it is elliptic (complex unimodular pair of eigenvalues) in one conjugate pair of directions.

The second invariant manifold is defined by those configurations where two particles are mirror images of each other while the third particle is restricted to the motion on the (arbitrarily chosen) symmetry line. Inside this manifold one finds mixed (i.e. partly regular and partly chaotic) dynamics. However, the motion is infinitely unstable in the transverse direction due to non-regularizable three-body collisions.

The quantum mechanics is done using the coordinates

$`\stackrel{}{x}`$ $`=`$ $`\left(\stackrel{}{r_1}+\stackrel{}{r_2}+\stackrel{}{r_3}\right)/\sqrt{3},`$ (5)   
$`\rho \mathrm{cos}{\displaystyle \frac{\theta ^{}}{2}}\left(\begin{array}{c}\mathrm{cos}(\varphi \phi ^{}/2)\\ \mathrm{sin}(\varphi \phi ^{}/2)\end{array}\right)`$ $`=`$ $`\left(\stackrel{}{r_1}\stackrel{}{r_2}\right)/\sqrt{2},`$ (8)   
$`\rho \mathrm{sin}{\displaystyle \frac{\theta ^{}}{2}}\left(\begin{array}{c}\mathrm{cos}(\varphi +\phi ^{}/2)\\ \mathrm{sin}(\varphi +\phi ^{}/2)\end{array}\right)`$ $`=`$ $`\left(\stackrel{}{r_1}+\stackrel{}{r_2}2\stackrel{}{r_3}\right)/\sqrt{6}.`$ (11)   
Here $`\rho ,\theta ^{}`$ and $`\phi ^{}`$ describe the intrinsic motion of the three-body system while $`\stackrel{}{x}`$ and $`\varphi `$ are the center of mass and the global orientation, respectively. In a second transformation we apply a rotation of $`\pi /2`$ around the abscissa corresponding to spherical coordinates $`(\rho ,\theta ^{},\phi ^{})`$, namely $`\mathrm{tan}\phi =\mathrm{cot}\theta ^{}/\mathrm{cos}\phi ^{}`$ and $`\mathrm{cos}\theta =\mathrm{sin}\theta ^{}\mathrm{sin}\phi ^{}`$ and obtain for the Laplacian in the subspace $`\stackrel{}{P}=0,L=0`$

$$\mathrm{\Delta }=\frac{^2}{\rho ^2}+\frac{3}{\rho }\frac{}{\rho }+\frac{4}{\rho ^2}\left(\frac{^2}{\theta ^2}+\mathrm{cot}\theta \frac{}{\theta }+\frac{1}{\mathrm{sin}^2\theta }\frac{^2}{\phi ^2}\right).$$
(12)   
Products of Bessel functions and spherical harmonics

$$\psi _{k,l,l_z}(\rho ,\theta ,\phi )=(k\rho )^1J_{2l+1}(k\rho )Y_l^{l_z}(\theta ,\phi )$$
(13)   
are eigenfunctions of the Laplacian (12)

$$\mathrm{\Delta }\psi _{k,l,l_z}(\rho ,\theta ,\phi )=k^2\psi _{k,l,l_z}(\rho ,\theta ,\phi ),$$
(14)   
with the usual relation between wavevector and energy $`k=\mathrm{}^1(2mE)^{1/2}`$. Fig. 2 shows a picture of the billiard taking $`(\rho ,\theta ,\phi )`$ as spherical coordinates. The billiard possesses a $`D_{3h}`$ symmetry. In the fundamental domain $`(\theta ,\phi )(0,\pi /2)\times (\pi /6,\pi /6)`$ the boundary is given by

$$\rho _B(\theta ,\phi )=\frac{a}{\sqrt{1+\mathrm{sin}\theta \mathrm{sin}(\phi +\pi /3)}}.$$
(15)   
The collinear manifold is the equatorial plane $`\theta =\pi /2`$. The second invariant manifold is given by the vertical symmetry planes $`\varphi =\pm \pi /6`$. Note that in this representation, classical geodesics of the billiard between two successive collisions are not straight lines since the centrifugal potential is stronger than in Euclidean case. In what follows we restrict ourselves to the fundamental domain and choose basis functions that fulfill Dirichlet boundary conditions. These are bosonic states.

We are interested in highly excited eigenstates. These may be accurately computed numerically by using the scaling method developed in ref.  and applied to a three-dimensional billiard by one of the authors . This method works efficiently only when a suitable positive weight function is introduced in a boundary integral. To this purpose we note that the radial part of (12) looks like a 4-dim Laplacian. Extending the results of refs.  to four dimensions yields the appropriate weight function, which has a remarkably simple form in our coordinates, namely we minimize the following functional

$$f[\mathrm{\Psi }_k]=_0^1d\mathrm{cos}\theta _{\pi /6}^{\pi /6}𝑑\phi \rho _B^4(\theta ,\phi )|\mathrm{\Psi }_k(\rho _B(\theta ,\phi ),\theta ,\phi )|^2$$

where the wave-function is expressed in terms scaling functions (13), $`\mathrm{\Psi }_k=_lc_{l,l_z}\psi _{k,l,l_z}`$. Due to our particular choice of boundary conditions we consider only the terms for which $`l+l_z`$ is odd and $`l_z=3m`$, and truncate at $`l=l_{\mathrm{max}}=ka/2+\mathrm{\Delta }lka/2`$.

We have computed three stretches of highly excited states. They consist of $`7430`$, $`1813`$ and $`2362`$ consecutive eigenstates with $`120<ka<235`$, $`290<ka<300`$ and $`393<ka<400`$, respectively. The last two stretches comprise levels with sequential quantum numbers around $`20000`$ and $`45000`$, respectively. The completeness of the series was checked by comparing the number of obtained eigenstates with the leading order prediction from the Weyl-formula $`\overline{d}(k)=c(24\pi ^2)^1(ak)^2,c0.51349`$.

Fig. 3 shows that the nearest neighbor spacing distribution agrees very well with RMT predictions already for the lower energy spectral stretch $`120ka235`$. The other series show well agreement, too. As for the long-range spectral correlations, the number variance $`\mathrm{\Sigma }^2(L)`$ deviates from RMT predictions for interval length of more than ten mean level spacings which we believe is due to the parabolic-elliptic family of periodic orbits in the collinear manifold (Fig. 1). The deviation from RMT decreases with increasing $`k`$. For the highest spectral stretch ($`ka400`$) the number variance increases linearly, $`\mathrm{\Sigma }^2(L)\mathrm{\Sigma }_{\mathrm{goe}}^2(L)+\epsilon L`$ up to $`L250`$, with $`\epsilon 0.04`$. This finding is consistent with the model of a statistically independent fraction $`\epsilon `$ of strongly scarred states . $`\mathrm{\Sigma }^2(L)`$ reaches its maximum and begins to oscillate at the saturation length $`L^{}`$ which scales as $`L^{}k^2`$ in agreement with the prediction of ref. .

The length spectrum $`D(r)`$, i.e. the cosine transform of the oscillatory part of the spectral density $`d_{\mathrm{osc}}(k)=_n\delta (kk_n)\overline{d}(k)`$, gives further information about long-range spectral fluctuations. For finite stretches of consecutive levels in the interval $`[k_1,k_2]`$ one uses a Welsh window function $`w(k;k_1,k_2)=(k_2k)(kk_1)/(6(k_2k_1)^3)`$ in the actual computation and obtains $`D(r)=_{k_1}^{k_2}𝑑kw(k;k_1,k_2)\mathrm{cos}(kr)d_{\mathrm{osc}}(k)`$ (see e.g. ). Fig. 4 shows that orbits of length $`r=\sqrt{2}a`$ and its integer multiples cause dominant peaks in the length spectrum.

To investigate the observed deviations from RMT predictions in more detail it is useful to compute the inverse participation ratio (IPR) of the wave functions in some basis . In the case of billiards the use of the angular momentum basis (13) is particularly convenient and suitable since periodic orbits correspond to sets of isolated points within this representation. Let $`c_{l,l_z}^{(n)}`$ denote the expansion coefficients of $`n`$th eigenstate $`\mathrm{\Psi }_{k_n}`$. We compute the IPR over a set of $`N`$ consecutive eigenstates as $`\mathrm{IPR}(l,l_z)=N_n|c_{l,l_z}^{(n)}|^4/(_n|c_{l,l_z}^{(n)}|^2)^2`$. The predicted RMT value for ideally quantum ergodic states is $`\mathrm{IPR}_{\mathrm{goe}}=3`$. Fig. 5 shows the IPR for the two sets of eigenstates with $`170<k<200`$ and $`290<k<300`$, respectively. The agreement with RMT predictions is rather good in both cases. This confirms that the billiard under consideration is dominantly chaotic and ergodic. However, the IPR is slightly enhanced in the region around $`l=l_zka/3`$. This is a robust phenomenon (present at all energy ranges), although the region of enhancement shrinks with increasing $`k`$. This finding is compatible with the expectation of uniform quantum ergodicity in the semi-classical limit. Note that the region $`ll_z`$ corresponds to the vicinity of the collinear manifold. Note further that the orbits belonging to the parabolic-elliptic family depicted in Fig. 1 have length $`\sqrt{2}a`$ and angular momenta in the region $`l/ka=l_z/ka(1/2\sqrt{6},1/\sqrt{6})`$. This is precisely the region where the IPR exhibits its enhancement while the orbits’ lengths coincide with the prominent peaks of the length spectrum in Fig. 4.

Thus, the deviations from RMT predictions observed for the spectrum and for the wave functions are associated with the family of parabolic-elliptic periodic orbits inside the collinear manifold. The special stability properties of this family lead to scars in the wave functions of the quantum system. This is an exciting new type of scars of invariant manifolds and complements results previously found in a three-dimensional billiard  and in interacting few-body systems . Note that the family of parabolic bouncing ball orbits inside the collinear manifold does not cause statistically detectable scarring. The orbits of this family correspond to points in angular momentum space with $`l=l_z<ka/2\sqrt{6}`$ and do not exhibit an enhancement in the IPR since the classical motion is too unstable in their vicinity.

In summary we have investigated an interacting three-body system realized as a billiard. Numerical results show that the classical dynamics is dominantly chaotic and no deviation from ergodic behavior is found. The spectral fluctuations of the quantum system agree well with random matrix theory predictions on energy scales of a few mean level spacings. However, wave function intensities and long ranged spectral fluctuations display deviations. These can be explained in terms of scars of a family of periodic orbits inside the collinear manifold.

We thank L. Kaplan for stimulating discussions. The hospitality of the Centro Internacional de Ciencias, Cuernavaca, Mexico, is gratefully acknowledged.