Resonant Population Transfer in the TimeDependent Quantum Elliptical Billiard
Abstract
We analyze the quantum dynamics of the timedependent elliptical billiard using the example of a certain breathing mode. A numerical method for the timepropagation of an arbitrary initial state is developed, based on a series of transformations thereby removing the timedependence of the boundary conditions. The timeevolution of the energies of different initial states is studied. The maximal and minimal energy that is reached during the timeevolution shows a series of resonances as a function of the applied driving frequency. At these resonances, higher (or lower) lying states are periodically populated, leading to the observed change in energy. The resonances occur when the driving frequency or a multiple of it matches exactly the mean energetic difference between the two involved states. This picture is confirmed by a fewlevel Rabilike model with periodic couplings, reproducing the key results of our numerical study.
pacs:
03.65.Ge, 03.65.Ge, 05.45.Mt, , ,
1 Introduction
Recently, the dynamical properties of classical timedependent twodimensional billiards have attracted major attention, especially in the context of Fermi acceleration [1, 2, 3, 4, 5, 6, 7, 8], which is defined as the unbounded energy gain of an ensemble of particles exposed to some external driving force [9]. Concerning the quantum dynamics of timedependent billiards, there exist several studies investigating the quantum version of the onedimensional FermiUlam model (or variants of it) [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21], but, to our knowledge, only two studies [22, 23] investigate timedependent billiards in higher dimensions.
Since already the onedimensional FermiUlam model involves timedependent Dirichlet boundary conditions which are difficult to treat, both analytically and numerically, most of the works analyze which (nonperiodic) movements of the wall allow for exact solutions: Linearly expanding or contracting wall motion is considered in this context in Refs. [10] and [14]. The authors of Refs. [12, 13] find that exact solutions are not only possible for a linear wall motion, but also for a timelaw of the form (with appropriate real constants and ). The same timelaw, that allows for analytic solutions, is found from a different perspective in Refs. [15, 18] and by means of a supersymmetry formalism in Ref. [20]. If additionally specific electromagnetic fields are superimposed, exact solutions can be obtained for arbitrary timedependencies of the wall motion [19]. Periodic driving laws are considered in Refs. [11] and [21]. The Floquet or quasienergy spectrum is found to be pure point like for most laws of the wall oscillation, resulting in a recurrent dynamics and thus to bounded energy growth [11]. Only certain (nonsmooth) driving laws yield a continuous quasienergy spectrum and thus allow for an unbounded acceleration, for example via resonance excitations, similarly to the quantum kicked rotator [24]. Recently, the authors of Ref. [21] studied the quantum FUM numerically by expanding the wave function in terms of the instantaneous eigenstates of the corresponding static system. They find that the dynamics of the timedependent expansion coefficients is irregular in an intermediate frequency regime, whereas for very low and very high frequencies a periodic behavior of the expansion coefficients is obtained.
Concerning higher dimensional billiards, the onepulse response of a twodimensional stadium billiard to a deformation of the boundary has been studied [23] by analyzing the evolving energy distribution. For small wall velocities, the spreading mechanism of this distribution is dominated by transitions between neighboring levels, while this is not the case for nonadiabatic wall velocities. In Ref. [22], the radially vibrating threedimensional spherical billiard is investigated. The authors claim ^{1}^{1}1Among others, their arguments are based on the correspondence between the classical and quantum 3D driven spherical billiard. By doing so, they state that in the classical system, the angular momentum, which is a constant of the motion in the static case, gets destroyed by the driving. This is not correct, since it is known, see for example Ref. [2], that in the radially oscillating twodimensional circular billiard the angular momentum is still conserved and the arguments can be generalized straightforwardly to the 3D spherical billiard. that only superpositions of two or more states that share the same common rotational symmetry yield quantum chaos, since the orthogonality relations of the instantaneous eigenstates allow in any other case a reduction to a timeindependent one degree of freedom Hamiltonian which cannot be chaotic.
In the present work, we investigate the quantum version of one of the most studied classical twodimensional timedependent billiards, the elliptical billiard with harmonically oscillating boundaries. Already the classical system shows several astonishing results: the driven elliptical billiard exhibits (tunable) Fermi acceleration [6] and the corresponding diffusion in momentum space experiences a dynamical crossover from sub to normal diffusion [25, 26]. A more technical appealing aspect of using the elliptical billiard is that the quantum version of the static billiard can be solved exactly in terms of the Mathieu functions, which allows an intuitive analysis in terms of the instantaneous eigenstates when driving the system. In particular, we will study the timeevolution of the energy in case of an initial eigenstate of the corresponding static billiard. The maximal and minimal energy that is reached as a function of the driving frequency show a series of resonances which we explain in detail by means of a detailed population analysis. We develop a Rabilike model with periodic couplings which nicely reproduces the results of our full numerical simulations, confirming the intuitive picture of the level dynamics of the driven elliptical billiard.
The work is structured as follows. In section 2 we summarize how the eigenstates of the static elliptical billiard can be obtained by using the Mathieu functions and discuss the symmetry properties of the eigenstates. Subsequently, we tackle the problem of how to solve the driven elliptical billiard numerically in section 3. The results of the simulations are presented in section 4 and interpreted, explained and modeled in section 5. Finally, a short summary and outlook is provided in section 6.
2 Static elliptical billiard
As a precursor for the analysis of the timedependent elliptical billiard, we briefly study the static elliptical billiard. In particular, we will use in the analysis of the timedependent system projections on eigenstates of the static ellipse and the symmetry properties of these eigenstates. To find the quantum mechanical eigenstates of a static elliptical billiard, we have to solve the corresponding stationary Schrödinger equation
(1) 
where is the eigenenergy belonging to and the Hamiltonian is given by
(2) 
where denotes the mass. The potential is zero inside and infinite outside the static elliptical billiard
(3) 
This means the eigenvalue equation (1) has to be solved under Dirichlet boundary conditions
(4) 
where and are the semimajor and semiminor axis of the ellipse, respectively.
Due to the boundary condition of Eq. (4), solving the eigenvalue equation (1) is actually not trivial and involves some subtleties [27, 28, 29, 30, 31, 32, 33]. Here, we follow the discussion provided in Ref. [33] and additionally use results presented in Refs. [32] and [31]. We start by transforming from Cartesian to elliptic coordinates:
(5a)  
(5b) 
where is the semifocal distance of the ellipse and . The stationary Schrödinger equation reads
(5f) 
where is an eigenfunction with energy and is the Hamiltonian:
(5g) 
The potential energy in elliptic coordinates is now
As a consequence, the eigenfunctions have to satisfy the Dirichlet boundary condition and also the periodicity condition . The Laplacian operator in elliptic coordinates is given by
(5h) 
The twodimensional (2D) Schrödinger equation inside the elliptical billiard (i.e. ) then simplifies to
(5i) 
with . Obviously, can be separated with the standard ansatz
(5j) 
into a radial () and an angular () part. Plugging this separation ansatz in the Schrödinger equation (5i), we obtain two ordinary differential equations:
(5ka)  
(5kb) 
where is the separation constant and is the dimensionless, rescaled energy
(5kl) 
Equations (5ka), (5kb) are the standard form of the Mathieu equations, Eq. (5kb) is the ordinary Mathieu equation (OME) and Eq. (5ka) is the modified Mathieu equation (MME). The corresponding solutions are the ordinary and modified Mathieu functions (OMF, MMF), respectively. Note that the change of variables transforms the OME into the MME. Even though the ansatz (5j) decouples the Schrödinger equation (5i), yielding two ordinary differential equations, the separation constants and do not decouple [34], as we will see in the following.
The physical solutions of the OME have to be periodic, i.e. they can be expanded in a Fourier series (actually, Floquet theory guarantees the existence of periodic solutions [31])
(5km) 
For a fixed , it is known that only certain values of allow for periodic solutions. These eigenvalues are called the characteristic values and , where is the order ( for and for ). Since the Mathieu equation is of SturmLiouville type, all eigenvalues are real, positive and can be ordered . In the plane, the functions and are curves that do not intersect. The and with even order have a period of , the ones with odd order have a period of . The order specifies the number of zeros of and in the interval . The expansion coefficients and are determined by the recurrence relations for the Mathieu equations [35].
The solutions of the MME Eq. (5ka), i.e. the radial part, can be obtained from the solutions of the OME by the mentioned change of variable , yielding
(5kn) 
Note that the expansion coefficients and are the same as the ones in the solutions (5km) of the OME. The solutions of the MME can also be viewed as an independent eigenvalue problem for the ’s (with fixed and ), now has to the satisfy the Dirichlet boundary condition .
The eigenvalues associated to and are denoted by and , where the index refers to the even and odd MMF and is the th zero of the MMF, not counting a possible root at . In the plane, the functions and are curves that again do not intersect. The solution of the Schrödinger equation (5i) is the product (5j) of the OMF and MMF. Both eigenvalue problems have to be solved simultaneously, since (or ) and (or ) in the OMF and MMF of the solution cannot be chosen independently. The pairs (or ) that satisfy both eigenvalue problems simultaneously are of course the crossing points of the family of curves and with the and . Thus, there are two types of crossing points that correspond to solutions of the simultaneous eigenvalues problem:

with (even modes)

with (odd modes)
These two types of crossing points are related to the symmetries of the solution . This yields the eigenfunctions
(5koa)  
(5kob) 
where and correspondingly . The indices and are the angular and radial quantum numbers of the even and odd modes, respectively. is the order of the OMF and specifies the number of zeros of the MMF in the interval , not counting the zeros at the origin i.e. is the th zero of the MMF of order . Since the MMF has to be zero at , it follows that is greater or equal to one. According to Eq. (5kl), the separation constant is directly proportional to the energy and the energy eigenvalues corresponding to the eigenfunctions can be written as
(5kop) 
where is the eccentricity of the ellipse.
Symmetries
The symmetries of the eigenstates of the elliptical billiard are determined by the OMF (5km) [31]. For an even OMF, i.e. for , the wavefunction is symmetric with respect to reflections at the axis (), this symmetry is denoted by . Odd OMF imply an antisymmetric eigenstate with respect to the axis (), denoted by . The symmetry of with respect to the axis is fixed by the angular quantum number , which is the order of the OMF. For even Mathieu functions, is symmetric with respect to reflections at the axis () if is even and antisymmetric () if is odd. Naturally, this symmetry is denoted by . For odd Mathieu functions, is symmetric to the axis if is odd and antisymmetric if is even. Overall, this yields four possible combinations concerning the parity of the eigenstates , which are summarized in the following list:

even solutions with even

even solutions with odd

odd solutions with odd

odd solutions with even
The four parity combinations correspond to the characterization of the symmetry reduced quarter elliptic billiard. However, the boundary conditions along the coordinate axes of the quarter billiard have to be adjusted according to the different parities: Dirichlet boundary conditions are required along the axis (axis) for () and Neumann boundary conditions along the axis (axis) for ().
Label  even/odd  

1  4,267  0  1  e  
2  9,058  1  1  e  
3  12,577  1  1  o  
4  15,993  2  1  e  
5  19,358  2  1  o  
6  25,061  3  1  e  
7  25,895  0  2  e  
8  27,998  3  1  o  
9  35,156  1  2  e  
10  36,178  4  1  e  
11  38,511  4  1  o  
12  44,040  1  2  o  
13  46,406  2  2  e  
14  49,199  5  1  e 
3 Driven elliptical billiard
To study the time evolution of an initial state in the harmonically oscillating elliptical billiard, we have to solve the timedependent Schrödinger equation
(5koq) 
where the Hamiltonian is given by
(5kor) 
and the potential is zero inside and infinity outside the timedependent elliptical billiard:
(5kos) 
This means the Schrödinger equation (5koq) has to be solved under Dirichlet boundary conditions
(5kot) 
where and are determined by the driving law of the elliptical billiard.
Transformations
The explicit timedependence of the boundary conditions is very difficult to be treated by standard numerical techniques. Therefore, to remove the explicit timedependence of the boundary conditions, we apply the following timedependent coordinate transformation, whose 1D variant has been successfully used to remove timedependent boundary conditions [10, 12, 13, 14, 18, 19, 20, 21] in 1D systems
(5kou) 
yielding the timedependent Schrödinger equation
(5kov) 
The advantage of applying this timedependent coordinate transformation is that the boundary condition in the new coordinates is extremely simple: the wave function has to vanish on the circle given by . The prize we have to pay is that now the coordinates themselves are explicitly timedependent, and . This has to be taken into account when applying the differential operator on in Eq. (5kov). The additional terms resulting from the left hand side of Eq. (5kov) can be put into an effective Hamiltonian, yielding
(5kow) 
and the corresponding timedependent Schrödinger equation
(5kox) 
To remove the terms proportional to and , we apply the following timedependent unitary transformation, yielding the new wave function
(5koy)  
(5koz) 
With this transformation we obtain the Schrödinger equation
(5koaa) 
with the effective Hamiltonian
(5koab) 
and the timeindependent boundary condition for . This effective Hamiltonian can be interpreted as a twodimensional anisotropic harmonic oscillator with timedependent frequencies and timedependent masses and the above boundary condition.
We expand in terms of the eigenfunctions of the (unit) static circular billiard, since this ansatz automatically fulfills the boundary condition. The eigenfunctions of the static circular billiard are best described in polar coordinates, so we transform and to polar coordinates:
(5koac) 
The Hamiltonian (3) in polar coordinates reads
(5koad) 
and the corresponding Schrödinger equation is
(5koae) 
Ansatz
As already written, we expand the wave function in terms of the eigenfunctions of the static circular billiard. To carry out numerical simulations, this expansion has to be truncated, yielding
(5koaf) 
with timedependent expansion coefficients , where is the radial and the azimuthal quantum number. The timeindependent eigenfunctions of the static circular billiard factorize
(5koag) 
The normalized radial and azimuthal (angular) functions are given by [36, 37]
(5koah)  
(5koai) 
where is the Bessel function of the first kind of order and is the th zero of .
To determine the timedependent expansion coefficients , we insert the ansatz (5koaf) into the timedependent Schrödinger equation (5koae), project with the bra
(5koaj) 
and thus obtain a set of coupled ordinary differential equations (ODE) for the . This coupled ODE system can be written in the following form, for details of the derivation see Ref. [38]
(5koak) 
The can be written in the following form (exemplarily shown for )
where the are simple timedependent functions of and , whereas the are timeindependent. For a definition of the , the and the , see A. More specific, the are composed of integrals over products of Bessel functions, which in general cannot be evaluated analytically. However, since the corresponding integrals are timeindependent, they have to be evaluated only once, which can be done numerically with high precision. Thus, by introducing a linear index the coupled ODE system (3) is in the canonical form and can be solved numerically by standard techniques.
Observables
The timeevolution of the expectation value of the energy of a certain initial state propagating in the driven elliptical is given by
(5koal) 
To calculate at a given time of the original laboratory system, we cannot directly use the wavefunction (5koaf) together with the Hamiltonian (3). Instead we have to apply the following steps: First, the inverse of the transformations (5kou), (5koy) and (5koac) operate on and we obtain
(5koam) 
Secondl, we apply the coordinate transformation as given in Eq. (5kou), but now as a timeindependent coordinate transformation, which just depends parametrically on the time :
(5koan) 
The wavefunction and the Hamiltonian are then given by
(5koao) 
(5koap) 
Finally, changing again to polar coordinates yields
(5koaq) 
and the corresponding Hamiltonian is given by Eq. (3) if we set , i.e.
(5koar) 
Another quantity of interest is the spectral composition of the energy in terms of the population of the instantaneous eigenstates of the elliptical billiard. By ‘instantaneous eigenstates’ we mean the following: Consider the boundary of the elliptical billiard for a certain time instant. If we now take this particular billiard configuration, treat it as a static ellipse and calculate the corresponding eigenstates assuming Dirichlet boundary conditions, we obtain the instantaneous eigenstates. If we denote the eigenstates by ( labels the degree of energetic excitation), the population of an instantaneous eigenstate is then of course simply the projection of the wave function onto this state, i.e. .
It is illustrative to analyze how the energy of the instantaneous eigenstates changes during one oscillation of the elliptical billiard. From Eq. (5kop) it can be seen that the energy of the eigenstates of the elliptical billiard depends (for a given semimajor axis ) on the eccentricity , so . Since the periodic (with period ) driving law , changes the eccentricity, we may write for the energy of the instantaneous eigenstates , where is the phase of the oscillation, measured in terms of , i.e. . for energetically lowest 14 states is shown in figure 1. There are several crossings, note that these are real crossings and not avoided crossings, see the inset. Such crossings are possible since the static elliptical billiard is an integrable system, thus the corresponding eigenstates posses ‘good’ quantum numbers, in agreement with the theorem of Wigner and von Neumann [39].
4 Results
In the following, we will take as an initial state an eigenstate of the static elliptical billiard, i.e. , and let it evolve in the driven ellipse. We choose for the driving law of the elliptical billiard and . For the equilibrium positions of the semimajor and semiminor axes we use , , resulting in an eccentricity of . The driving amplitude is set to , whereas the frequency will be varied.
In figure 2, the evolution of the energy is shown for a driving frequency of and two different initial states, namely the ground state (first state) and the third excited state (fourth state), i.e. and , respectively. These are the first two states that have the same symmetry, namely , see table 1. When starting in the instantaneous ground state (at ), the wave function follows the motion of the elliptic boundary adiabatically. In other words, the projection of on the instantaneous ground state is approximately one, . The periodic oscillations of the energy are then simply present due to the fact that the energy of the instantaneous ground state changes periodically (approximately sinusoidal), cf. figure 1.
The situation is different when starting in the fourth state (). Even though performs sinusoidallike oscillations corresponding to the energy oscillations of the fourth instantaneous eigenstate, there is a small modulation (beating) on top, with a period much larger then the period of the applied driving law. This largeperiodic modulation is better visible in the timeevolution of the population coefficients , which are shown in figure 3. Since , and all the other are zero. The fourth instantaneous eigenstate stays populated () throughout the whole timeevolution, however  besides small, fast oscillations with the period of the applied driving law  there is now a beating with a much larger period . Correspondingly, other states get noticeably populated, primarily and are now greater than zero, with having the same large beating period . Note that the instantaneous eigenstates corresponding to , and share the same symmetry properties , cf. Table 1.
When choosing other low energy eigenfunctions as an initial state, e.g. , , or , the result is very much like the one when starting with , in the sense that the energy performs periodic oscillations matching the sinusoidallike oscillations of the corresponding instantaneous eigenstate. This behavior is also observed for smaller driving frequencies (), now also when starting in the fourth state, i.e. .
For larger driving frequencies, the beating phenomenon becomes even more pronounced. In figure 4, a driving frequency of is applied to and . For the resulting shows besides the large amplitude beating some irregular fluctuations, however stays tightly within a welldefined band. Now for the beating is much more pronounced. Again, the period of the beating is much larger than the driving period , note however that it is not the same beating period we obtained in the case of , cf. figure 3. Like in the case , the timeevolution of the population coefficients, see figure 6, reflects the beating phenomenon. Again (for ), the coefficient is the one right behind that gets populated the most, note that at some times even . Other states with noticeable population are , and , again these are all states with symmetries , cf. Table 1.
By increasing the driving frequency further to , the beating phenomenon is now observed equally when starting in the instantaneous ground state () and when starting in the fourth state (), see figure 5. Unlike in the case , the largeperiodic modulation (in the case ) considerably lowers rather than rises the energy . Consequently, the state besides that gets substantially populated is now and not the higher energetic , see figure 7. Even more, reaches almost one, i.e. almost the whole population is periodically transferred between to .
The beating phenomenon gets destroyed for sufficiently high driving frequencies, see figure 8, where for with () and () is shown. Now a considerable number of states with the same symmetry as the initial state get noticeably populated. Consequently, the energy does not stay in a narrow interval as in the case of small driving frequencies, but rather fluctuates irregularly over a wide range. More specific, increases rapidly within the first few oscillations of the elliptical billiard, then slowly grows further until and saturates for . This means that the energy stays bounded, unlike in the corresponding classical system [25, 26], there is no Fermi acceleration in the quantum version. This is due to dynamical localization: the wave function gets exponentially localized in momentum space, thereby stopping the energy diffusion.
In general, when starting at with an instantaneous eigenstate with energy , the maximal and minimal energy of is a good indicator whether the timeevolution of is adiabatic or not. If the evolution of is adiabatic, are simply given by the corresponding eigenvalues of the instantaneous eigenstate, more specific, in this case () and , cf. figure 1. As soon as the timeevolution is not purely adiabatic, the beating phenomenon sets in and either is raised, or is lowered, or both. This can be seen in figure 9, where and are shown as a function of the driving frequency when starting with the fourth state (). For small driving frequencies (), and are given by the energies of the fourth instantaneous eigenstate, i.e. and , respectively, which means that the timeevolution is adiabatic. For higher driving frequencies, a series of resonances appear: has sharp resonances at and broad resonances with double peaks at and . On the other hand, shows considerably fewer resonances, sharp ones at and a broad one at . The fact that has considerably fewer resonances than is immediately clear, since the only way that the energy, when starting in the fourth state (), can be lowered is to populate , whereas many different energetically higher lying instantaneous states can be populated, leading to an increase of .
To gain an understanding of the above presented results, in particular of the resonances of and , we develop a Rabilike model in the following section.
10.1833  5.0917  3.3944  2.5458  
10.3364  5.1682  3.4455  2.5841  
10.7058  5.3529  3.5686  2.6764  
11.9227  5.9613  3.9742  2.9807  
17.5667  8.7833  5.8556  4.3917  
20.5198  10.2599  6.8399  5.1299  
20.8891  10.4446  6.9630  5.2223  
22.2591  11.1296  7.4197  5.5648  
28.2724  14.1362  9.4241  7.0681  
31.2255  15.6128  10.4085  7.8064  
32.4425  16.2212  10.8142  8.1106  
38.4558  19.2279  12.8186  9.6139  
43.1482  21.5741  14.3827  10.7871  
48.7922  24.3961  16.2641  12.1981  
60.7149  30.3575  20.2383  15.1787 