Resonant Population Transfer in the Time-Dependent Quantum Elliptical Billiard

We analyze the quantum dynamics of the time-dependent elliptical billiard using the example of a certain breathing mode. A numerical method for the time-propagation of an arbitrary initial state is developed, based on a series of transformations thereby removing the time-dependence of the boundary conditions. The time-evolution of the energies of different initial states is studied. The maximal and minimal energy that is reached during the time-evolution 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 few-level Rabi-like model with periodic couplings, reproducing the key results of our numerical study.


Introduction
Recently, the dynamical properties of classical time-dependent two-dimensional 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 time-dependent billiards, there exist several studies investigating the quantum version of the one-dimensional Fermi-Ulam 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 time-dependent billiards in higher dimensions.
Since already the one-dimensional Fermi-Ulam model involves time-dependent Dirichlet boundary conditions which are difficult to treat, both analytically and numerically, most of the works analyze which (non-periodic) 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 time-law of the form l(t) = √ at 2 + 2bt + c (with appropriate real constants a, b and c). The same time-law, 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 electro-magnetic fields are superimposed, exact solutions can be obtained for arbitrary time-dependencies 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 (non-smooth) 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 time-dependent 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 one-pulse 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 non-adiabatic wall velocities. In Ref. [22], the radially vibrating three-dimensional spherical billiard is investigated. The authors claim ‡ that ‡ Among 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 two-dimensional circular billiard the angular momentum is still conserved and the arguments can be generalized straightforwardly to the 3D spherical billiard. 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 time-independent 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 two-dimensional time-dependent 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 time-evolution 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 Rabi-like 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.

Static elliptical billiard
As a precursor for the analysis of the time-dependent 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 Ψ n (x, y) of a static elliptical billiard, we have to solve the corresponding stationary Schrödinger equation where E is the eigenenergy belonging to Ψ(x, y) and the Hamiltonian H(x, y) is given by where µ denotes the mass. The potential V (x, y) is zero inside and infinite outside the static elliptical billiard This means the eigenvalue equation (1) has to be solved under Dirichlet boundary conditions where a and b are the semi-major and semi-minor 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: where f = √ a 2 − b 2 is the semi-focal distance of the ellipse and ξ 0 = arctanh(b/a). The stationary Schrödinger equation reads where ϕ(ξ, η) is an eigenfunction with energy E and H(ξ, η) is the Hamiltonian: The potential energy V (ξ, η) in elliptic coordinates is now As a consequence, the eigenfunctions ϕ(ξ, η) have to satisfy the Dirichlet boundary condition ϕ(ξ 0 , η) = 0 and also the periodicity condition ϕ(ξ, η) = ϕ(ξ, η + 2π). The Laplacian operator in elliptic coordinates is given by The two-dimensional (2D) Schrödinger equation inside the elliptical billiard (i.e. ξ ≤ ξ 0 ) then simplifies to with k 2 = 2µE/ 2 . Obviously, ϕ(ξ, η) can be separated with the standard ansatz into a radial (R(ξ)) and an angular (Θ(η)) part. Plugging this separation ansatz in the Schrödinger equation (9), we obtain two ordinary differential equations: where α is the separation constant and q is the dimensionless, rescaled energy Equations (11a), (11b) are the standard form of the Mathieu equations, Eq. (11b) is the ordinary Mathieu equation (OME) and Eq. (11a) 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 η = iξ transforms the OME into the MME. Even though the ansatz (10) decouples the Schrödinger equation (9), yielding two ordinary differential equations, the separation constants α and q 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]) For a fixed q, it is known that only certain values of α allow for periodic solutions. These eigenvalues are called the characteristic values α l and β l , where l is the order (l ≥ 0 for ce l (η) and l ≥ 1 for se l (η)). Since the Mathieu equation is of Sturm-Liouville type, all eigenvalues are real, positive and can be ordered α 0 < β 1 < α 1 < β 2 . . .. In the q × α plane, the functions α l (q) and β l (q) are curves that do not intersect. The ce l (η, q) and se l (η, q) with even order have a period of π, the ones with odd order have a period of 2π. The order l specifies the number of zeros of ce l (η, q) and se l (η, q) in the interval η ∈ [0, π]. The expansion coefficients A k (l, q) and B k (l, q) are determined by the recurrence relations for the Mathieu equations [35]. The solutions of the MME Eq. (11a), i.e. the radial part, can be obtained from the solutions of the OME by the mentioned change of variable η = iξ, yielding Note that the expansion coefficients A k (l, q) and B k (l, q) are the same as the ones in the solutions (13) of the OME. The solutions of the MME can also be viewed as an independent eigenvalue problem for the q's (with fixed α and β), now R l (ξ) has to the satisfy the Dirichlet boundary condition R l (ξ 0 ) = 0. The eigenvalues associated to α and β are denoted by q e r (α) and q o r (β), where the index o, e refers to the even and odd MMF and r ≥ 1 is the rth zero of the MMF, not counting a possible root at ξ = 0. In the q × α plane, the functions q e r (α) and q o r (β) are curves that again do not intersect. The solution ϕ(ξ, η) of the Schrödinger equation (9) is the product (10) of the OMF and MMF. Both eigenvalue problems have to be solved simultaneously, since α (or β) and q e (or q o ) in the OMF and MMF of the solution ϕ(ξ, η) cannot be chosen independently. The pairs (α, q) (or (β, q)) that satisfy both eigenvalue problems simultaneously are of course the crossing points of the family of curves α l (q) and β l (q) with the q e r (α) and q o r (β). Thus, there are two types of crossing points that correspond to solutions of the simultaneous eigenvalues problem: (i) α r (q) with q e m (α) (even modes) (ii) β r (q) with q o m (β) (odd modes) These two types of crossing points are related to the symmetries of the solution ϕ(ξ, η). This yields the eigenfunctions ϕ l,r (ξ, η) ϕ e l,r (ξ, η) = Ce l (ξ, q e l,r )ce l (η, q e l,r ) (15a) where q e l,r = q e r (α l ) and correspondingly q o l,r = q o r (β l ). The indices l and r are the angular and radial quantum numbers of the even and odd modes, respectively. l is the order of the OMF and r specifies the number of zeros of the MMF in the interval ξ ∈ [0, ξ 0 ], not counting the zeros at the origin i.e. q e,o l,r is the r-th zero of the MMF of order l. Since the MMF has to be zero at ξ 0 , it follows that r is greater or equal to one. According to Eq. (12), the separation constant q is directly proportional to the energy and the energy eigenvalues E e,o l,r corresponding to the eigenfunctions ϕ e,o l,r (ξ, η) can be written as where ǫ is the eccentricity of the ellipse.

Symmetries
The symmetries of the eigenstates Ψ e,o l,r (x, y) of the elliptical billiard are determined by the OMF Θ l (η) (13) [31]. For an even OMF, i.e. for ce l (η), the wavefunction Ψ e l,r (x, y) is symmetric with respect to reflections at the x-axis (Ψ e l,r (x, y) = Ψ e l,r (x, −y)), this symmetry is denoted by π y = +1. Odd OMF se l+1 (η) imply an antisymmetric eigenstate Ψ o l,r (x, y) with respect to the x-axis (Ψ e l,r (x, y) = −Ψ e l,r (x, −y)), denoted by π y = −1. The symmetry of Ψ e,o l,r (x, y) with respect to the y-axis is fixed by the angular quantum number l, which is the order of the OMF. For even Mathieu functions, Ψ e l,r (x, y) is symmetric with respect to reflections at the y-axis (Ψ o l,r (x, y) = Ψ o l,r (−x, y)) if l is even and antisymmetric (Ψ o l,r (x, y) = −Ψ o l,r (−x, y)) if l is odd. Naturally, this symmetry is denoted by π x . For odd Mathieu functions, Ψ o l,r (x, y) is symmetric to the y-axis if l is odd and antisymmetric if l is even. Overall, this yields four possible combinations concerning the parity of the eigenstates Ψ e,o l,r (x, y), which are summarized in the following list: (i) π y = +1, π x = +1 ⇒ Ψ e l,r (x, y) = Ψ e l,r (−x, −y) even solutions with even l Table 1. All eigenstates Ψ e,o l,r (x, y) of the elliptical billiard (a = 1, b = √ 0.51 ⇒ eccentricity ǫ = 0.7) with E < 50 (in units of 2 /µ). The eigenstates are ordered according to the energy. π y (π x ) denotes the symmetry of Ψ e,o l,r (x, y) with respect to reflections at the x-axis (y-axis).
Label E[ 2 /µ] l r even/odd π y π x 1 4,267 0 odd solutions with even l 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 x-axis (y-axis) for π y = −1 (π x = −1) and Neumann boundary conditions along the x-axis (y-axis) for π y = +1 (π x = +1).

Driven elliptical billiard
To study the time evolution of an initial state Ψ 0 (x, y) := Ψ(x, y, t = 0) in the harmonically oscillating elliptical billiard, we have to solve the time-dependent Schrödinger equation where the Hamiltonian H(x, y, t) is given by and the potential V (x, y, t) is zero inside and infinity outside the time-dependent elliptical billiard: This means the Schrödinger equation (17) has to be solved under Dirichlet boundary conditions where a(t) and b(t) are determined by the driving law of the elliptical billiard.

Transformations
The explicit time-dependence of the boundary conditions is very difficult to be treated by standard numerical techniques. Therefore, to remove the explicit time-dependence of the boundary conditions, we apply the following time-dependent coordinate transformation, whose 1D variant has been successfully used to remove time-dependent boundary conditions [10,12,13,14,18,19,20,21] in 1D systems yielding the time-dependent Schrödinger equation The advantage of applying this time-dependent coordinate transformation is that the boundary condition in the new coordinates η, ξ is extremely simple: the wave function Ψ(η, ξ, t) has to vanish on the circle given by η 2 + ξ 2 = 1. The prize we have to pay is that now the coordinates themselves are explicitly time-dependent, η = η(t) and ξ = ξ(t). This has to be taken into account when applying the differential operator ∂/∂t on Ψ(η(t), ξ(t), t) in Eq. (22). The additional terms resulting from the left hand side of Eq. (22) can be put into an effective Hamiltonian, yielding and the corresponding time-dependent Schrödinger equation To remove the terms proportional to i η ∂ ∂η and i ξ ∂ ∂ξ , we apply the following timedependent unitary transformation, yielding the new wave function Λ(η, ξ, t) With this transformation we obtain the Schrödinger equation with the effective Hamiltonian and the time-independent boundary condition Λ(η, ξ, t) = 0 for η 2 + ξ 2 = 1. This effective Hamiltonian can be interpreted as a two-dimensional anisotropic harmonic oscillator with time-dependent frequencies and time-dependent masses and the above boundary condition. We expand Λ(η, ξ, t) in terms of the eigenfunctions Φ n,m of the (unit) static circular billiard, since this ansatz automatically fulfills the boundary condition. The eigenfunctions Φ n,m of the static circular billiard are best described in polar coordinates, so we transform η and ξ to polar coordinates: The Hamiltonian (28) in polar coordinates reads and the corresponding Schrödinger equation is

Ansatz
As already written, we expand the wave function in terms of the eigenfunctions Φ n,m (r, φ) of the static circular billiard. To carry out numerical simulations, this expansion has to be truncated, yielding with time-dependent expansion coefficients c n,m (t), where n is the radial and m the azimuthal quantum number. The time-independent eigenfunctions Φ n,m (r, φ) of the static circular billiard factorize The normalized radial and azimuthal (angular) functions are given by [36,37] where J m is the Bessel function of the first kind of order m and k m,n is the nth zero of J m .
To determine the time-dependent expansion coefficients c n,m (t), we insert the ansatz (32) into the time-dependent Schrödinger equation (31), project with the bra Λ| and thus obtain a set of coupled ordinary differential equations (ODE) for the c n,m (t). This coupled ODE system can be written in the following form, for details of the derivation see Ref. [38] i ċ n,m (t) = The d nmn ′ (t) can be written in the following form (exemplarily shown for d nmn ′ are composed of integrals over products of Bessel functions, which in general cannot be evaluated analytically. However, since the corresponding integrals are time-independent, they have to be evaluated only once, which can be done numerically with high precision. Thus, by introducing a linear index (n, m) ↔ i the coupled ODE system (37) is in the canonical form and can be solved numerically by standard techniques.

Observables
The time-evolution of the expectation value of the energy of a certain initial state propagating in the driven elliptical is given by To calculate E(t) at a given time t 0 of the original laboratory system, we cannot directly use the wavefunction Λ(r, φ, t 0 ) (32) together with the Hamiltonian H e (r, φ, t 0 ) (30).
Instead we have to apply the following steps: First, the inverse of the transformations (21), (25) and (29) operate on Λ(r, φ, t 0 ) and we obtain Secondl, we apply the coordinate transformation as given in Eq. (21), but now as a time-independent coordinate transformation, which just depends parametrically on the time t 0 : The wavefunction and the Hamiltonian are then given by Finally, changing again to polar coordinates yields Ψ(r, φ, t 0 ) = e 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 Ψ i (t)| (i labels the degree of energetic excitation), the population p i of an instantaneous eigenstate is then of course simply the projection of the wave function Ψ(x, y, t) onto this state, i.e. p i (t) = | Ψ i (t)|Ψ(t) |.
It is illustrative to analyze how the energy of the instantaneous eigenstates changes during one oscillation of the elliptical billiard. From Eq. (16) it can be seen that the energy of the eigenstates of the elliptical billiard depends (for a given semi-major axis a) on the eccentricity ǫ, so E i = E i (ǫ). Since the periodic (with period T = 2π/ω) driving law a(t), b(t) changes the eccentricity, we may write for the energy of the instantaneous eigenstates E i (ζ), where ζ is the phase of the oscillation, measured in terms of T , i.e.    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].

Results
In the following, we will take as an initial state Ψ 0 = Ψ(t = 0) an eigenstate of the static elliptical billiard, i.e. Ψ 0 = Ψ e,o l,r (x, y), and let it evolve in the driven ellipse. We choose for the driving law of the elliptical billiard a(t) = a 0 + c sin(ωt) and b(t) = b 0 + c sin(ωt). For the equilibrium positions of the semi-major and semi-minor axes we use a 0 = 1, b 0 = √ 0.51, resulting in an eccentricity of ǫ 0 = 0.7. The driving amplitude c is set to 0.1, whereas the frequency ω will be varied.
In figure 2, the evolution of the energy E(t) is shown for a driving frequency of ω = 5 and two different initial states, namely the ground state (first state) and the third excited state (fourth state), i.e. Ψ 0 = Ψ e 0,1 (x, y) and Ψ 0 = Ψ e 2,1 (x, y), respectively. These are the first two states that have the same symmetry, namely π x = π y = +1, see table 1. When starting in the instantaneous ground state (at t = 0), the wave function Ψ(t) follows the motion of the elliptic boundary adiabatically. In other words, the projection p 1 (t) of Ψ(t) on the instantaneous ground state is approximately one, p 1 (t) = | Ψ e 0,1 (t)|Ψ(t) | ≈ 1. The periodic oscillations of the energy E(t) 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 (Ψ 0 = Ψ e 2,1 (x, y)). Even though E(t) performs sinusoidal-like oscillations corresponding to the energy oscillations of the fourth instantaneous eigenstate, there is a small modulation (beating) on top, with a period T b > T = 2π/ω much larger then the period of the applied driving law. This large-periodic modulation is better visible in the time-evolution of the population coefficients p i (t), which are shown in figure 3. Since Ψ 0 = Ψ e 2,1 (x, y), p 4 (0) = 1 and   all the other p i (0) are zero. The fourth instantaneous eigenstate stays populated (p 4 (t) > 0.98, ∀t) throughout the whole time-evolution, however -besides small, fast oscillations with the period T of the applied driving law -there is now a beating with a much larger period T b . Correspondingly, other states get noticeably populated, primarily p 7 (t) and p 1 (t) are now greater than zero, with p 7 having the same large beating period T b . Note that the instantaneous eigenstates corresponding to p 1 , p 4 and p 7 share the same symmetry properties π x = π y = +1, cf. Table 1. When choosing other low energy eigenfunctions as an initial state, e.g. Ψ 0 = Ψ e 1,1 (x, y), Ψ o 1,1 (x, y), Ψ o 2,1 (x, y) or Ψ 0 = Ψ e 3,1 (x, y), the result is very much like the one when starting with Ψ 0 = Ψ e 0,1 (x, y), in the sense that the energy E(t) performs periodic oscillations matching the sinusoidal-like oscillations of the corresponding instantaneous eigenstate. This behavior is also observed for smaller driving frequencies (ω < 4), now also when starting in the fourth state, i.e. Ψ 0 = Ψ e 2,1 (x, y). For larger driving frequencies, the beating phenomenon becomes even more pronounced. In figure 4, a driving frequency of ω = 10 is applied to Ψ 0 = Ψ e 0,1 (x, y) and Ψ 0 = Ψ e 2,1 (x, y). For Ψ 0 = Ψ e 0,1 (x, y) the resulting E(t) shows besides the large amplitude beating some irregular fluctuations, however E(t) stays tightly within a welldefined band. Now for Ψ 0 = Ψ e 2,1 (x, y) the beating is much more pronounced. Again, the period T b of the beating is much larger than the driving period T , note however that it is not the same beating period we obtained in the case of ω = 5, cf. figure 3. Like in the case ω = 5, the time-evolution of the population coefficients, see figure 6, reflects the beating phenomenon. Again (for Ψ 0 = Ψ e 2,1 (x, y)), the coefficient p 7 is the one right behind p 4 that gets populated the most, note that at some times even p 7 > p 4 . Other states with noticeable population are p 10 , p 1 and p 13 , again these are all states with symmetries π x = π y = +1, cf. Table 1.
By increasing the driving frequency further to ω = 12, the beating phenomenon is now observed equally when starting in the instantaneous ground state (Ψ 0 = Ψ e 0,1 (x, y)) and when starting in the fourth state (Ψ 0 = Ψ e 2,1 (x, y)), see figure 5. Unlike in the case ω = 10, the large-periodic modulation (in the case Ψ 0 = Ψ e 2,1 (x, y)) considerably lowers rather than rises the energy E(t). Consequently, the state besides p 4 that gets substantially populated is now p 1 and not the higher energetic p 7 , see figure 7. Even more, p 1 reaches almost one, i.e. almost the whole population is periodically transferred between p 4 to p 1 .
The beating phenomenon gets destroyed for sufficiently high driving frequencies, see figure 8, where E(t) for ω = 100 with Ψ 0 = Ψ e 0,1 (x, y) (p 1 (0) = 1) and Ψ 0 = Ψ e 2,1 (x, y) (p 4 (0) = 1) 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, E(t) increases rapidly within the first few oscillations of the elliptical billiard, then slowly grows further until t ≈ 60 T and saturates for t > 60 T . 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 t = 0 with an instantaneous eigenstate with energy E i , the maximal E max = max(E(t)) and minimal energy E min = min(E(t)) of Ψ(t) is a good indicator whether the time-evolution of Ψ(t) is adiabatic or not. If the evolution of Ψ(t) is adiabatic, E max,min are simply given by the corresponding eigenvalues of the instantaneous eigenstate, more specific, in this case (ω = 1) E max = E i (ζ = 3/4) and E min = E i (ζ = 1/4), cf. figure 1. As soon as the time-evolution is not purely adiabatic, the beating phenomenon sets in and either E max is raised, or E min is lowered, or both. This can be seen in figure 9, where E max and E min are shown as a function of the driving frequency ω when starting with the fourth state (Ψ 0 = Ψ e 0 ,1 (x, y)). For small    E max is immediately clear, since the only way that the energy, when starting in the fourth state (p 4 (0) = 1), can be lowered is to populate p 1 , whereas many different energetically higher lying instantaneous states can be populated, leading to an increase of E max .
To gain an understanding of the above presented results, in particular of the resonances of E max and E min , we develop a Rabi-like model in the following section. Table 2. (color online) Mean (averaged of the phase ζ) energy differences E i↔j together with E i↔j /2, E i↔j /3 and E i↔j /4 for the first six states that share the same symmetry π x = π y = +1. The energy differences which are related to the resonances of figure 9 between just two states are shaded in gray, the ones related to the large double peak at ω ≈ 10.5 in yellow and the ones of the second double peak at ω ≈ 15 . . . ω ≈ 16 in green.

Interpretation and Rabi-like model
The aim is now to explain the resonances in figure 9 firstly by means of a population analysis of the instantaneous eigenstates and secondly via a Rabi-like model with timeperiodic coupling.
Since the instantaneous eigenstates and the corresponding energies in the elliptical billiard are quantized, only certain states can be populated when driving with a specific frequency ω. In particular, we expect an efficient transfer between two states, when their energy difference exactly matches the driving frequency, i.e. when E j − E i = ω. As we saw earlier, the energy E i of an instantaneous eigenstate depends on the phase ζ, it is thus useful to define the mean energy difference E i↔j between two instantaneous eigenstates as These mean differences together with E i↔j /2, E i↔j /3 and E i↔j /4 are shown in table 2 for the first six states that share the same symmetry π x = π y = +1. Now if we start in the fourth state (p 4 (0) = 1), cf. figure 9, the only way to lower the energy is to populate the ground state (p 1 ). Thus, when E 1↔4 ≈ ω, we expect a significant population of the ground state and consequently a decrease of E min . From table 2 E 1↔4 = 11.92 and in figure 9, there is a broad resonance around ω = 11.89. It is not surprising that the resonance is rather broad, since E 1↔4 is just the mean energy difference, actually the difference between the first and the forth state varies between max(E 4 (ζ)−E 1 (ζ)) = 14.6 and min(E 4 (ζ) − E 1 (ζ)) = 9.6, so transitions between the two states should be possible over a broad range around ω = 11.89 and this is exactly what we observe. The transition between the first and the fourth state is a one photon transition in the sense that ω = E 1↔4 , so one time the driving frequency matches the energy difference, yielding a broad resonance. On the other hand, there are also multi-photon transitions possible, such that n · ω = E 1↔4 , n = 2, 3, . . .. Since the efficiency of multi-photon processes drastically reduces with the number of photons involved, these higher-order resonances should be much sharper and are naturally harder to detect. The sharp resonances of E min at ω = 5.95 and ω = 3.97 are the two and three photon resonances of E 1↔4 , respectively, and the numbers are in excellent agreement with the corresponding values of table 2.
In a similar fashion, the sharp resonances of E max can be explained: ω = 2.58, 3.45, 5.17 correspond to the four, three, and two photon transitions of E 4↔7 , respectively and ω = 6.83 to the three photon transition of E 4↔10 , again in excellent agreement with the results of table 2. Note that a population analysis of all the resonances confirms the so far presented considerations, e.g. at ω = 6.83 besides p 4 , p 10 gets populated, whereas at ω = 5.17, besides p 4 now this is the case for p 7 .
Applied to the broad double peaks of E max , a population analysis yields the following: at ω = 10.19 a mixture of p 4 , p 7 and p 10 is present (also p 13 , p 1 and p 20 are noticeably, however less than the first three ones populated), whereas at ω = 10.51 mostly the states p 7 , p 4 , p 20 , p 13 , p 10 and p 1 contribute. The yellow shielded fields of table 2 indicate that there is a multitude of one-to four-photon transitions in this frequency regime between exactly these states that according to the population analysis contribute to the resonance. Note that there are now also processes of the kind p 4 → p 7 → p 10 possible, since the corresponding energy differences E 4↔7 and E 7↔10 are comparable. In contrast to the just described resonance, the broad double peak at ω = 15.6 and ω = 16.4 is mainly due to transitions between just two states. The corresponding population analysis yields p 4 and p 13 for ω = 15.6 and p 4 and p 18 for ω = 16.4, again in excellent agreement with the two-photon transition of E 4↔13 and the three-photon transition of E 4↔18 , see the greenish fields of table 2.
For all the resonances, the population analysis additionally shows that the beating period T b increases dramatically when approaching a resonance, i.e. the population transfer from one instantaneous eigenstate to another becomes slower and slower the closer the driving frequency gets to the corresponding resonance.
In the discussion of the resonances, we already used the projection onto the instantaneous eigenstates to analyze the results. With this picture in mind (see e.g. figure 3, where just two states are considerably populated), a very crude simplification is to interpret the driven elliptical billiard as a two-level system (two instantaneous eigenstates), where the two levels are coupled by a periodic driving with an unknown coupling strength a. Without the driving, the Hamiltonian of the unperturbed system in a suitable basis is then simply given by where E 1 and E 2 are the energies of the two eigenstates Ψ 0 1 and Ψ 0 2 , with time-evolution Ψ 0 n (t) = Ψ 0 n e −iEnt/ , n = 1, 2. Next, we take into account the driving of the system, which leads to a periodic coupling of the two states. In the Hamiltonian, this leads to off-diagonal elements of the form a sin(ωt). However, as can be seen in figure 1, the energy itself of the instantaneous eigenstates depends on time. To this end, we introduce time-dependent energy shifts △E i (t) on the diagonal. In a good approximation, the △E i (t) are harmonic functions, which we obtain by fitting sine functions to the curves of figure 1. Overall, the model Hamiltonian is then where H 1 (t) = △E 1 (t) a sin(ωt) a sin(ωt) △E 2 (t) .
The standard ansatz for the time-dependent wave function Ψ(t) is a superposition with time-dependent coefficients of the unperturbed wave functions Ψ 0 i (t), i.e.
Plugging this ansatz into the time-dependent Schrödinger equation leads to a system of ordinary differential equations for the coefficientṡ Note that at this point, usually (e.g. in atomic systems) the rotating wave approximation is applied, however, in our case this would lead to incorrect results, so we have to solve the full system given by Eqs. (5), which is only possible numerically. Not surprisingly, this simple 2-level model does not lead to correct results, e.g. the 2-level model cannot reproduce all the resonances of figure 9. The obvious reason is, as the population analysis shows, that already for driving frequencies in the range shown in figure 9 (ω < 17) more than just two states are involved. The straightforward generalization to n-levels yields: , or in other words H 1 ii (t) = △E i (t) and for the off-diagonal elements we have H 1 ij (t) = a sin(ωt). It is reasonable to assume that the coupling strength a is equal between all the levels, leaving us with just one free parameter in the model. Again, the time-evolution of the unperturbed states is given by H 0 Ψ 0 n = E n Ψ 0 n ⇒ Ψ 0 n (t) = Ψ 0 n e −iEnt/ and the superposition ansatz reads yielding a set of coupled linear differential equations for the coefficients, which has to be solved numerically: With this Rabi-like model, we calculate the time evolution of Ψ(t), with the same initial conditions as in the discussion of the resonances, i.e. we start in the fourth state, c 4 (0) = 1. We take into account the first six states that can couple, i.e. the states with symmetries π x = π y = +1, namely 1, 4, 7, 10, 13 and 18, so we utilize a six-level Rabilike model. Just like for the full numerical calculations, we compute E max and E min as a function of the applied driving frequency ω. The result is shown in figure 10. Overall, there is very good agreement with the results of the full simulations, cf. figure 9. In particular, all resonances of figure 9 are reproduced. The additional resonances that appear in the Rabi-like model are four-, five-and even six-photon processes: ω = 2.07 corresponds to the five-photon transition of E 4↔7 , ω = 7.81 to the four-photon transition of E 4↔13 and ω = 12.21 to the four-photon transition of E 4↔18 . The resonance at ω = 5.17 is superimposed with a very sharp resonance at ω = 5.20. This is the sixphoton transition of E 4↔13 , additionally, we suppose that this resonance is enhanced by the combination of the two-photon transition E 4↔7 with the four-photon transition E 7↔13 , again cf. table 2.
The additional resonances are extremely sharp and this is one reason why they do not appear in figure 9. A much finer frequency scan would be necessary in the case of the full numerical simulations. Secondly, as already mentioned, the population transfer between two states becomes extremely slow when approaching a resonance (the beating period T b increases), i.e. the simulations have to be carried out over longer times in order to detect the very sharp resonances.
The double peak at ω = 10.19, 10.51 of figure 9 appears in the Rabi-like model as a single, broad resonance, see figure 10. The reason is that the time-dependent energy shifts △E i (t) of the model are obtained by fitting sine functions to the energy levels E i (ζ) of figure 1. This fitting procedure does preserve the mean distances between the energy levels, in particular it reproduces the correct E i↔j , however, the minima of △E i (t) are lowered and the maxima of △E i (t) are raised compared to E i (ζ). As a consequence, the resonances are broadened in the Rabi-like model, such that the mentioned double-peak merges, so they cannot be resolved within the model. Note that when fitting the E i (ζ) in a different fashion, such that the minima and maxima are preserved, the double peak does appear in the model, however this alternative fitting procedure changes the E i↔j , so all the resonances are shifted.
The fact that the beating period T b increases drastically when approaching a resonance is also reproduced by the Rabi-like model, however, unlike for the appearance of the resonances itself, we do not have an intuitive explanation for these beatings. Nevertheless, the success of the Rabi-like model shows that we have identified the essential ingredients to explain the time-evolution of a wave function, when starting in an eigenstate, in the driven elliptical billiard.

Conclusion and Outlook
In this work, we have developed a numerical method to propagate an arbitrary initial state in the time-dependent elliptical billiard characterized by time-dependent boundary conditions. The method is based on a series of transformations of the original Hamiltonian, removing the time-dependent boundary conditions. Subsequently expanding the wave function in a suitable basis, we obtain a large system of coupled ordinary differential equations for the time-dependent expansion coefficients which can be solved by standard techniques. We then propagate eigenstates of the static system and investigate the time-evolution of their energy in the driven billiard. Unlike in the corresponding classical system, the energy does not grow unboundedly but rather saturates, even for large driving frequencies, i.e. there is no Fermi acceleration in the quantum billiard. The maximal and minimal energy that is reached as a function of the applied driving frequency show a series of resonances. We analyze these resonances by projecting the wave function onto the instantaneous eigenstates, yielding the spectral composition of the states. At these resonances, population is periodically, with a period T b , transferred between the initial state and another one, either increasing or decreasing the energy of the wave function. This analysis shows that the resonances appear exactly when the driving frequency ω matches the mean energetic difference between the involved instantaneous states (one photon transition) or multiples of the driving frequency (n-photon transitions, which result in very sharp resonances). The beating period T b of the population transfer is much larger than the period T = 2π/ω of the driving and increases dramatically when approaching a resonance. To confirm our intuitive picture describing the appearance of the resonances, we develop a fewlevel Rabi-like model with time-periodic couplings, and time-dependent energy shifts on the diagonal, resembling the phase dependence of the energies of the instantaneous eigenstates. This model reproduces nicely the essential results of our full numerical simulations, in particular the mentioned resonances and the occurrence of the beating period T b .
So far, we investigated the breathing mode of the elliptical billiard, i.e. the timedependence of the semi-major and semi-minor axes is given by a(t) = a 0 + c sin(ωt) and b(t) = b 0 + c sin(ωt), respectively. A promising alternative would be to use a timelaw obeying a(t)b(t) = const., since this preserves the area of the ellipse. As a result, the instantaneous energy levels exhibit a multitude of crossings, naturally posing the question whether a(t) could be adjusted in such a way that the momentum diffusion does not stop, yielding quantum Fermi acceleration in the driven elliptical billiard.