# Influence of phase space localization on the energy diffusion in a quantum chaotic billiard

## I Introduction

The quantum dynamics of a classically chaotic system is in the present a topic of very active interest. Specially, with the great development in artificially fabricated small devices as quantum dots , where quantum manifestation of classical chaos  plays an important role.

Systems governed by a parameter dependent Hamiltonian $`H(\mathrm{})`$ are excellent models to study quantum manifestations of classical chaos. Dynamical localization  and phase space localization (scarring phenomena)  have been investigated in those systems. Recently, Wilkinson  studied the dispersion of the energy in a generic non integrable system when the parameter is time dependent. He observed that it spreads diffusively and found asymptotic theoretical expresions (using random matrix theory) for the diffusion constant in the limits of large and small velocities of the parameter. Numerical experiments performed in a random matrix model were in good agreement with those predictions. The physical motivation of these studies is the possibility of modelling the quantum dissipation of a finite-size system of non-interacting fermions, based on that microscopic diffusive behaviour of the energy dispersion.

In this paper, we treat the dispersion of the energy in a more realistic system: a two dimensional chaotic billiard with moving boundary. It is not a simple numerical task to solve the time dependent Schöedinger equation for a 2-D billiard with moving boundary , specially, it is very lengthy from the computationally point of view. We found a simple shape parameter Hamiltonian expansion which enables us to obtain the spectrum (up to second order) and wave functions of the deformed billiard for deformations so large as the characteristic wave length in the working energy region. Using that expansion and defining a shape oscillatory motion, the quantum dynamic of a particle inside the billiard reduces to a system of coupled linear differential equations which can be integrated directly. On the other hand, this expansion can be used to obtain efficiently static properties of the spectrum like curvature distributions or avoided crossing distributions.

We studied the quantum dynamics of the Bunimovich stadium billiard with oscillating walls and calculated the diffusion constant $`D`$ as a function of the boundary velocity. Differences with theoretical predictions were found. For small velocities of the boundary, the system diffuses more than predictions. In this regime, the existence of bouncing ball states enables coherent transport of probability, so the diffusion is enhanced. On the other hand, for large velocities, the system diffuses less than predictions. We propose phase space localization as the possible mechanism to reduce diffusion in this regime. These assertions are supported with numerical simulations on the stadium billiard where the most localized states in the region under study were extracted. In this case, the system behaves in good agreement with the theory.

The outline of the paper is as follows. To make it self-contained, in Sec. II we give a short introduction to the main results obtained previously for the diffusion of the energy in generic chaotic systems. In Sec. III we show how to expand the Hamiltonian of a 2-D planar billiard in powers of the shape parameter. In Sec. IV we give a detailed numerical study of the diffusion constant $`D`$ as a function of the velocity of the boundary for the stadium billiard. Section IV is devoted to final remarks. We include an appendix with the expresions required to obtain the parameter Hamiltonian expansion in the stadium billiard.

## II Quantum energy diffusion in chaotic systems

Let $`H(\mathrm{})`$ be a parametric Hamiltonian related to a generic chaotic system for all $`\mathrm{}`$. The parameter $`\mathrm{}`$ is a time dependent function $`\mathrm{}(t)`$. We restrict to systems with time reversal symmetry where the statistical properties of the spectrum are well described by the eigenvalues of ensambles of gaussian orthogonal random matrices (GOE) . At any $`\mathrm{}(t)`$, the system admits an energy spectrum given by the eigenvalue problem for the instantaneous Hamiltonian

$$\widehat{H}(\mathrm{}(t))\psi _\mu (\mathrm{}(t),𝐫)=\frac{k_\mu ^2(\mathrm{}(t))\mathrm{}^2}{2m}\psi _\mu (\mathrm{}(t),𝐫).$$
(1)   
The system is prepared at $`t=0`$ in a highly excited eigenstate $`\psi _\nu (\mathrm{}(0),𝐫)`$ of the Hamiltonian $`\widehat{H}(\mathrm{}(0))`$. The state for time $`t`$ can be expressed in the basis of eigenstates of the instantaneus Hamiltonian (eq. 1), the so called adiabatic basis,

$`\mathrm{\Psi }(𝐫,t)={\displaystyle \underset{\mu }{}}a_\mu (t)\psi _\mu (\mathrm{}(t),𝐫).`$

The dispersion of the energy for the considered state is:

$$\mathrm{\Delta }^2E(t)=\underset{\mu }{}|a_\mu (t)|^2(E_\mu (\mathrm{}(t))E(t))^2,$$
(2)   
where $`E(t)\mathrm{\Psi }|E|\mathrm{\Psi }=_\mu |a_\mu (t)|^2E_\mu (\mathrm{}(t))`$ is the expectation value of the energy as a function of $`t`$.

Wilkinson  observed that the dispersion of the energy exhibit a diffusive growth when it is averaged over many states, that is

$`\mathrm{\Delta }^2E(t)=2Dt.`$

Using random matrix theory for the statistical properties of the spectrum, he predicted that the diffusion constant $`D`$ is an universal function of $`\rho `$ (the mean energy density), $`\mathrm{}`$, $`\dot{\mathrm{}}`$ and the mean value of the off diagonal elements of $`H/\mathrm{}`$ that he called $`\sigma `$. The dimensionless parameter $`\kappa =\rho ^2\mathrm{}\dot{\mathrm{}}\sigma `$ is a measure of the adiabaticity of the variation in the Hamiltonian. He found the following asymptotic form of the diffusion constant for small and large $`\kappa `$:

$`D=\{\begin{array}{cc}2^{5/4}\pi \mathrm{\Gamma }(3/4)\frac{1}{\rho ^3\mathrm{}}\kappa ^{3/2}\hfill & \text{if }\kappa 1,\hfill \\ \pi \frac{1}{\rho ^3\mathrm{}}\kappa ^2\hfill & \text{if }\kappa 1\text{ }.\hfill \end{array}`$ (5)   
## III Hamiltonian expansion for deformed billiards

This section is devoted to obtain a Hamiltonian expansion for a particle of mass $`m`$ inside a shape parameter dependent planar billiard. For variations $`\delta \mathrm{}\mathrm{}_0`$ of the parameter around $`\mathrm{}_0`$, we propose the following expansion in the basis of eigenfunctions at $`\mathrm{}_0`$,

$$H_{\mu \nu }(\mathrm{}_0+\delta )(H_{\mu \mu }+\delta ^2H_{\mu \mu }^{\prime \prime }/2)\delta _{\mu \nu }+\delta H_{\mu \nu }^{},$$
(6)   
with $`\delta _{\mu \mu }=1`$, and $`\delta _{\mu \nu }=0`$ for $`\mu \nu `$. This allows us to express the eigenfunctions of the deformed billiard in terms of the eigenfunctions at $`\mathrm{}_0`$. Then, when the parameter changes as a function of time, the dynamical evolution of the particle can be integrated easily using a standard Runge-Kutta method (see Sec. IV).

Let $`\zeta `$ be a smooth close curve defining a planar billiard. We use a curvilinear coordinate system around the boundary with $`s`$ along $`\zeta `$ and $`z`$ perpendicular to it at $`s`$ ($`z=0`$ on $`\zeta `$). Consider now that the boundary is deformed and the changes are parametrized by

$$𝐫(s,\delta )=𝐫_\mathrm{𝟎}(s)+z(s,\delta )𝐧,$$
(7)   
with $`𝐫_\mathrm{𝟎}(s)`$ the parametric equation for $`\zeta `$, and $`𝐧`$ the outward normal unit vector to $`\zeta `$ at $`𝐫_0(s)`$ (see Fig. 1). In the appendix we give $`z(s,\delta )`$ for the stadium billiard.

Let $`\{\varphi _\mu (𝐫);\mu =1,2,\mathrm{}N\}`$ be the eigenfunctions of the billiard defined by $`\zeta `$ with eigenwave numbers $`k_\mu `$ around $`k_0`$ ($`|k_\mu k_0|\mathrm{\Delta }k\mathrm{Perimeter}/\mathrm{Area}`$). Taking $`\mathrm{}^2/2m=1`$, we immediately see from Eq. (6) that $`H_{\mu \mu }=k_\mu ^2`$. The eigenfunctions vanish on $`\zeta `$ and can be extended outside the billiard in a smooth way; up to second order as follows:

$$\varphi _\mu (𝐫_\mathrm{𝟎}(s)+z𝐧)=z\frac{\varphi _\mu }{𝐧}(𝐫_\mathrm{𝟎}(s))+𝒪(z^3).$$
(8)   
Moreover, to each function $`\varphi _\mu (𝐫)`$ we associate the scaling function $`\varphi _\mu (k𝐫/k_\mu )`$. This family of functions depending on the scaling parameter $`k`$ verifies Helmholtz equation with wave number $`k`$.

For small variations of the boundary ($`\delta 1`$), it is valid to use perturbation theory to obtain the eigenfunctions and eigenvalues of the deformed billiard. With those aproximated solutions we will obtain a Hamiltonian expansion in powers of $`\delta `$ which allows us to extend the range of deformations to the order of the wave lenght ($`2\pi /k_0`$).

Let $`\psi _\mu (𝐫,\delta )`$ be the eigenfunction of the deformed billiard obtained from $`\varphi _\mu (𝐫)`$ by a continuos variation of the parameter. As we go to obtain the spectrum up to second order in $`\delta `$, we express $`\psi _\mu `$ in terms of $`\varphi _\nu `$ up to first order

$`\psi _\mu (𝐫,\delta )=\varphi _\mu \left({\displaystyle \frac{k_\mu (\delta )}{k_\mu }}𝐫\right)+\delta {\displaystyle \underset{\nu (\mu )}{}}c_{\mu \nu }\varphi _\nu \left({\displaystyle \frac{k_\mu (\delta )}{k_\nu }}𝐫\right),`$

with $`k_\mu (\delta )`$ the wave number for which $`\psi _\mu (𝐫,\delta )`$ vanishes on the deformed boundary (Eq. (7)). Then, by expanding around $`\mathrm{}_0`$ we have

$`0={\displaystyle \frac{d\psi _\mu }{d\mathrm{}}}(𝐫_\mathrm{𝟎}(s),0)=`$ $`[k_\mu ^{}r_n/k_\mu +z^{}(s)]{\displaystyle \frac{\varphi _\mu }{𝐧}}(𝐫_\mathrm{𝟎}(s))+`$ (9)   
$`{\displaystyle \underset{\nu (\mu )}{}}c_{\mu \nu }\varphi _\nu (k_\mu 𝐫_\mathrm{𝟎}(s)/k_\nu ),`$ (10)   
where primes indicate derivation with respect to $`\mathrm{}`$ at $`\mathrm{}_0`$ and $`r_n𝐫.𝐧`$. In what follows, we omit the argument $`(𝐫_\mathrm{𝟎}(s))`$ of the normal derivatives and drop terms including derivatives of second order or more. In particular, we take $`\varphi _\nu (k_\mu 𝐫_\mathrm{𝟎}/k_\nu )(k_\mu /k_\nu 1)\varphi _\nu /𝐧`$. Then, multipying Eq. (10) by $`\varphi _\mu /𝐧`$ and integrating on $`\zeta `$, we have

$$H_{\mu \mu }^{}(k_\mu ^2)^{}=_\zeta z^{}(s)(\frac{\varphi _\mu }{𝐧})^2𝑑s.$$
(11)   
To obtain the last expresion it is necessary to use the following quasiorthogonality relation 

$`{\displaystyle _\zeta }{\displaystyle \frac{\varphi _\mu }{𝐧}}{\displaystyle \frac{\varphi _\nu }{𝐧}}{\displaystyle \frac{r_nds}{2k_\mu k_\nu }}=\delta _{\mu \nu }+{\displaystyle \frac{(k_\mu k_\nu )}{(k_\mu +k_\nu )}}𝒪(1).`$

Now, multiplying Eq. (10) by $`\varphi _\nu /𝐧`$ and integrating on $`\zeta `$, we obtain

$`c_{\mu \nu }(k_\mu ^2k_\nu ^2)={\displaystyle \frac{(k_\mu +k_\nu )}{2k_\nu }}{\displaystyle _\zeta }z^{}(s){\displaystyle \frac{\varphi _\mu }{𝐧}}{\displaystyle \frac{\varphi _\nu }{𝐧}}𝑑s.`$

By perturbation theory the last expresion would be $`H_{\mu \nu }^{}`$; however, it is not symmetrical with respect to $`\mu `$ and $`\nu `$. Working to first order in $`\delta `$ it is sufficient to replace the factor $`(k_\mu +k_\nu )/2k_\nu `$ by one (we assume that $`|k_\mu k_\nu |k_0`$); but to second order it is not the case. Consider two states at the same distance from $`k_\mu `$; that is $`k_\pm =k_\mu \pm ϵ`$. Then, the factor in question (to first order in $`ϵ`$) is $`1ϵ/2k_\mu `$. To solve the problem, we define the symmetrical factors

$`A_{\mu \nu }=2(k_\mu +k_\nu )/2k_0,`$

which give the same asymmetrical contribution to the state $`\mu `$ from states $`\pm `$. On the other hand, we stress that for generic non-integrable billiards and for any non trivial deformation (a dilation is trivial) the interaction is of long range. So, to have into acount the finite dimension of the basis we suggest to multiply the elements of Eq. (11) by the following cut-off:

$`\mathrm{Cf}_{\mu \nu }=\mathrm{exp}\left[2(k_\mu ^2k_\nu ^2)^2/(k_0\mathrm{\Delta }k)^2\right].`$

Finally, non diagonal Hamiltonian elements acquire the expresion

$`H_{\mu \nu }^{}=\mathrm{Cf}_{\mu \nu }A_{\mu \nu }{\displaystyle _\zeta }z^{}(s){\displaystyle \frac{\varphi _\mu }{𝐧}}{\displaystyle \frac{\varphi _\nu }{𝐧}}𝑑s.`$

Working as before, to obtain $`H_{\mu \mu }^{\prime \prime }`$ we multilply equation $`\mathrm{\hspace{0.33em}0}=d^2\psi _\mu /d\mathrm{}^2(𝐫_\mathrm{𝟎}(s),0)`$ by $`\varphi _\mu /𝐧`$ and integrate it on $`\zeta `$. The result is

$`H_{\mu \mu }^{\prime \prime }={\displaystyle \frac{3(H_{\mu \mu }^{})^2}{2H_{\mu \mu }}}{\displaystyle _\zeta }z^{\prime \prime }(s)({\displaystyle \frac{\varphi _\mu }{𝐧}})^2𝑑s.`$

Figure 2 compares the approximated spectrum obtained from eq. (6) with the exact one for the Bunimovich stadium billiard; the agreement is excellent.

## IV Numerical results

We have consider the desymmetrized stadium billiard with radius $`r`$ and straight line of lenght $`a`$. The boundary only depends on the shape parameter $`\mathrm{}=a/r`$ (the area is fixed to the value $`1+\pi /4`$). The parameter $`\mathrm{}`$ oscillates harmonically around $`\mathrm{}=1`$ with frequency $`\omega `$ and amplitude $`\alpha `$, that is,

$`\mathrm{}(t)=1+\alpha \mathrm{sin}(\omega t).`$

Sec. II shows universal expresions for $`D`$ in the asymptotic limits of $`\kappa `$; therefore, $`\kappa `$ is a good quantity to characterize the quantum dynamics of a specific chaotic system. In order to calculate the coefficient $`D`$, we solved the time dependent Schroedinger equation for different initial conditions

$$i\mathrm{}\frac{\mathrm{\Psi }(t,𝐫)}{t}=\widehat{H}(\mathrm{}(t))\mathrm{\Psi }(t,𝐫).$$
(12)   
A very efficient way (computationally) to solve 12 is using the parameter expansion of the Hamiltonian obtained in Sec. III. The solution $`\mathrm{\Psi }(𝐫,t)`$ is expanded in the basis of eigenfunctions at $`\mathrm{}(0)=1`$,

$$\mathrm{\Psi }(𝐫,t)=\underset{\mu }{}b_\mu (t)\varphi _\mu (𝐫).$$
(13)   
Although $`\varphi _\mu (𝐫)`$ satifies boundary conditions for $`\mathrm{}(0)=1`$, it does not cancel outside the boundary (see Eq. 8). Of course, the amplitude of the oscillations is limited to values of the order of the wave lenght. After replacing 13 in 12, we obtain the system of differential equations

$`\dot{b}_\mu (t)={\displaystyle \frac{i}{\mathrm{}}}{\displaystyle \underset{\nu }{}}H_{\mu \nu }(\mathrm{}(t))b_\nu .`$

This system was integrated using a standard fourth-order Runge-Kutta method, and the dispersion of the energy is evaluated using (2) at times $`t=2n\pi /\omega `$ ($`n=0,1,2,\mathrm{}`$). We take $`\mathrm{}=1`$ and $`m=1/2`$. Then, we used a basis of 103 odd-odd eigenfunctions around $`k_0=47.3`$ at $`\mathrm{}=1`$ (see Fig. 2). Time evolutions were stopped before the wave functions had spread into the tails of the spectrum.

As an example, Fig. 3 shows the dispersion of the energy $`<\mathrm{\Delta }^2E>`$ for $`\omega =0.4`$ and $`\alpha =0.05`$. The average is over ten initial states. They were chosen near the center of the spectrum. It is observed clearly a linear spreading (diffusive behavior) of the dispersion of the energy for the first oscillations of the boundary (the slope being equal to $`2D`$); later, the spread goes slowly and eventually saturates.

Figure 4 shows a log-log plot of the difussion constant $`D`$ as a function of $`\kappa `$ which is proporsional to the mean value of the boundary velocity $`<|\dot{\mathrm{}}|>=2\alpha \omega /\pi `$. Simulations were performed for $`\omega =0.05,\mathrm{\hspace{0.33em}0.1},\mathrm{\hspace{0.33em}0.2},\mathrm{\hspace{0.33em}0.4},\mathrm{\hspace{0.33em}1},\mathrm{\hspace{0.33em}2},\mathrm{\hspace{0.33em}5}`$, $`10`$, $`20`$ and $`40`$, with $`\alpha =0.05`$ (filled circles). In all the cases the average was over ten initial conditions. Errorbars were calculated using the standard deviation of the average in the ensamble of initial states. The figure shows theoretical predictions in solid straight lines for a generic chaotic system with time reversal symmetry (Eq. 5), and we can observe significant differences with the numerical data.

In order to directly see the influence of localized eigenfunctions on the quantum dynamics of the system, we extracted from the spectrum the most localized states. This procedure is very simple by using the Hamiltonian expansion presented in Sec. III. In order to extract state $`j`$, we remove simply row and column $`j`$ from the Hamiltonian matrix. If at $`\mathrm{}=1`$ this state collide in an avoided crossing it is necessary to do first the transformation explained in reference . We have removed bouncing ball states with $`k=46.4589`$, $`47.1943`$, $`47.4220`$, and $`47.7548`$, all of them at $`\mathrm{}=1`$ (see Fig. 2). The resulting spectrum is showed in Fig. 5. In the same way as before, we have calculated the difusion constant for $`\omega =0.2,\mathrm{\hspace{0.33em}0.4},\mathrm{\hspace{0.33em}1},\mathrm{\hspace{0.33em}5},\mathrm{\hspace{0.33em}10}`$ and $`20`$ with $`\alpha =0.05`$ , and for $`\omega =0.2,\mathrm{\hspace{0.33em}0.4},\mathrm{\hspace{0.33em}10},\mathrm{\hspace{0.33em}20},\mathrm{\hspace{0.33em}40}`$ and $`60`$ with $`\alpha =0.09`$. In this case, the results are very close to theoretical predictions except for the highest values of $`\omega `$ (see Fig. 4).

## V Final Remarks

We have studied the spreading of the energy dispersion in the Bunimovich stadium billiard with oscillating walls. To the best of our knowledge this is the first computation of diffusion on a quantum chaotic system. We arrived to that possibility by using a Hamiltonian expansion which we develop in this work.

Numerical simulations show that the energy spreads diffusively for a number of oscillations (linear behaviour) and then saturates. This saturation phenomenon is associated to the fact that the eigenfunctions of the evolution operator for one period (Floquet states), which are localized in energy, inhibit the spread of energy for long times . We calculated the diffusion coefficient (the slope of the initial linear behaviour) for different velocities of the boundary, observing significant differences with theoretical predictions based on random matrix theory.

In the adiabatic regime (slow boundary velocities), the system diffuse more than predictions. We associate these differences to the existence of states strongly localized in phase space (their Wigner or Hussimi distributions are strongly localized). In particular, we refer to bouncing ball states which are strongly localized in momentum space. The interaction between these states and generic chaotic ones is much smaller than the typical interaction between generic states. This fact is reflected on the occurrence in the spectrum of structures surviving parametric variation; that is, straight lines interrupted by very small avoided crossings (see Fig. 2). Then, having in mind that in this regime Landau-Zener transitions  are the mechanism of diffusion, those straight lines behave like channels allowing a coherent transport of probability to a wide range of energy in each oscillation. In order to confirm numerically the above explanation we have extracted from the spectrum those bouncing ball states living in the region under study. In this case, the system diffuse according to the theory. Certainly, classical billiards with a bouncing ball continuous family of orbits display anomalous diffusion too .

As the velocity increases, Landau-Zener probability transitions across small avoided crossings approach exponentially one. Then, if the system starts in a bouncing ball state, it moves essentially up and down by a straight line without mixing. And this situation mimics the evolution of a regular system where diffusion is known to be very small . In the spectrum without bouncing balls states, observed diffusion agrees with theory (in Fig. 4 the solid line with highest slope) reasonably well. However, for high velocities simulations and theory fall apart again showing the beginning of saturation for the diffusion process.

Recently we have proposed  the elimination of avoided crossings as the natural mechanism to uncover localized structures embedded in the eigenfunctions of chaotic Hamiltonian systems. Precisely, we have shown that many scars of short periodic orbits are uncover in the stadium billiard transforming the parametric spectrum in a set of smooth curves which cross among them. As we mention above, the elimination of avoided crossings can be carried out dynamically by increasing the boundary velocity; that is, when the probability transition is practically one. In the light of this, saturation of diffusion is expected in chaotic systems for high velocities as a consequence of these localized structures associated to short periodic orbits.

In conclusion, we have verified that asymptotic expresions (5) for the energy diffusion in a chaotic system work very well in a wide range of boundary velocities except for high velocities where diffusion saturates. Moreover, if the system includes a fraction of localized states like bouncing ball (or states living in a classically regular region which are connected to chaotic states by tunnelling), diffusion will be anomalous in the adiabatic regime.

## ACKNOWLEDGMENTS

This work was partially supported by UBACYT (TW35), PICT97 03-00050-01015 and SECYT-ECOS. We would like to thank A. Fendrik, C. Lewenkopf and M. Saraceno for useful discussions.

## A

In this appendix we show the parametrization used to describe the deformation of the desymmetrized stadium billiard. The area is fixed to the value $`A=1+\pi /4`$, so the boundary only depends on the shape parameter $`\mathrm{}=a/r`$. A point $`𝐫_\mathrm{𝟎}=(x,y)`$ on the boundary at $`\mathrm{}_0=1`$ is given in terms of the curvilinear coordinate system defined in Sec. II (see Fig. 1)

$`x(s)=\{\begin{array}{cc}s\hfill & \text{if }s1,\hfill \\ 1+\mathrm{sin}(s1)\hfill & \text{if }s>1\text{ }.\hfill \end{array}`$

$`y(s)=\{\begin{array}{cc}1\hfill & \text{if }s1,\hfill \\ \mathrm{cos}(s1)\hfill & \text{if }s>1\text{ }.\hfill \end{array}`$

Then, if the deformed boundaries are described by Eq. (7), it is a geometrical problem to show that $`z(s,\delta )\delta z^{}(s)+\delta ^2z^{\prime \prime }(s)/2`$, with

$`z^{}(s)=\{\begin{array}{cc}\frac{1}{2A}\hfill & \text{if }s1,\hfill \\ & \\ (1\frac{1}{2A})\mathrm{sin}(s1)\frac{1}{2A}\hfill & \text{if }s>1\text{ }.\hfill \end{array}`$

and

$`z^{\prime \prime }(s)=\{\begin{array}{cc}\frac{3}{4A^2}\hfill & \text{if }s1,\hfill \\ & \\ & \\ \frac{1}{A}1+(1\frac{1}{2A})^2\mathrm{sin}^2(s1)+\hfill & \\ & \\ \frac{1}{2A^2}\frac{1}{A}(1\frac{3}{4A})\mathrm{sin}(s1)\hfill & \text{if }s>1\text{ }.\hfill \end{array}`$