Evolutionary phase space in driven elliptical billiards

We perform the first long-time exploration of the classical dynamics of a driven billiard with a four dimensional phase space. With increasing velocity of the ensemble we observe an evolution from a large chaotic sea with stickiness due to regular islands to thin chaotic channels with diffusive motion leading to Fermi acceleration. As a surprising consequence, we encounter a crossover, which is not parameter induced but rather occurs dynamically, from amplitude dependent tunable subdiffusion to universal normal diffusion in momentum space. In the high velocity case we observe particle focusing in phase space.

Billiards are a widespread paradigm to study complex dynamics in many areas of physics. Two dimensional systems are especially appealing, since unlike one dimensional ones they possess already a rich phase space structure, which can range from integrable over mixed to fully chaotic [1]. Billiard systems can be realized experimentally in different ways, such as using semi-conductor heterostuctures [2], quantum dots [3], atom optical [4] and (superconducting) microwave setups [5,6]. Very recent applications of their dynamics are the investigation of directed emission from laser-microcavities [7,8,9] and the design of improved thermoelectric efficiency, where billiards are used to tailor the desired microscopic properties [10].
A natural generalization are time-dependent billiards, where the boundary is driven, for example harmonically. This allows the study of non-equilibrium processes such as Fermi acceleration (FA) [11,12,13]. There are only few studies of such two dimensional (2D) driven billiards [14,15,16,17,18,19,20], since the long-time numerical iteration of the corresponding four dimensional (4D) implicit mapping is computationally challenging and in particular the resulting phase space is hard to analyze due to its high dimensionality and unboundedness. So far this prevented any investigation addressing the long-time behavior of corresponding dynamical quantities. Key questions are whether there will be FA [11], i.e., unbounded energy gain of an ensemble exposed to time-dependent external potentials and more specifically what laws this FA will follow depending on the geometry and the driving of the billiard. In one dimension, the existence of FA depends exclusively on the driving law [21], in 2D the situation is much more complicated, since the dynamics of the corresponding static systems can be integrable, mixed or chaotic. Recently the existence of FA in the harmonically driven elliptical billiard was demonstrated [20], augmenting the "LRA-conjecture" which says that only time-dependent systems with a mixed or chaotic static counterpart exhibit FA [15]. However, none of the mentioned works analyzed the full 4D phase space and investigated how its structure determines the diffusion in momentum space, especially with respect to the long time behavior of the diffusion and acceleration process.
In this Letter, we present, for the first time, a detailed exploration of the long-term dynamics of a driven billiard with a 4D phase space. It is shown that the phase space geometry evolves, with increasing velocity of the propagating ensemble, from a large chaotic sea with rich stickiness structures due to regular islands, to large dominating regular domains with thin chaotic acceleration channels. This change in the phase space structure leads to a crossover from amplitude dependent tunable subdiffusion to universal normal diffusion in momentum space. Due to its acceleration with increasing number of collisions the particle ensemble propagates effectively in an evolutionary phase space which we study in depth.
The geometry of the time-dependent elliptical billiard is given by where (x(t), y(t)) is a point on the boundary, φ is the corresponding elliptic polar angle, C > 0 is the driving amplitude, ω is the frequency of oscillation and A 0 , B 0 are the equilibrium values of the long and the short half-diameter, respectively. We fix ω = 1, A 0 = 2 and B 0 = 1 without loss of generality (arbitrary units). Due to ballistic motion in between collisions, the dynamics can be described by an implicit 4D map [19] using the variables (ξ n , φ n , α n , v n ), where v n = |v n | is the modulus of the particles velocity, ξ n = ωt mod 2π is the phase of the boundary oscillation and α n is the angle between v n and the tangent on the boundary at the n-th collision. Note that the corresponding static ellipse is integrable [1], since besides the energy the product of the angular momenta around the two foci where ε is the eccentricity, is conserved. The static phase space is globally divided by the separatrix (F = 0), connecting the two hyperbolic fixed points, into rotators (F > 0) and librators (F < 0). Two elliptic fixed points are located at the minima F min = −ε 2 /(1 − ε 2 ). Periodic orbits with short periods are located around the separatrix, whereas the existence of period orbits in the librator region close to the elliptic fixed points requires comparably long periods [14]. This will later be important for the understanding of the distribution of periodic orbits in the driven case. In ref. [20] we showed that the driven ellipse exhibits Fermi acceleration (FA). Here we present a comprehensive analysis of the 4D phase space, specifically we explore how its composition changes with increasing v and demonstrate how this affects the long-time evolution of FA in the system. We first analyze the periodic orbits and their stability and consider subsequently the phase space density ρ of a typical ensemble in detail, by exploring its dependence on the number of collisions n as well as the velocity v. Due to the occurrence of FA, the phase space is open with respect to v. Unlike in the static system, α is not restricted to [0, π], but can now range from 0 to 2π [14]. Nevertheless, the total phase space is not just spanned by ξ, φ, α and v, but possesses a more complicated topology, since e.g. for ξ ∈ [π/2, 3π/2] (the ellipse is contracting) α is restricted to α ∈ [0, π] [22].
Finding (unstable) periodic orbits in such a high dimensional phase space is an intricate task and becomes more difficult with increasing velocity of the orbits, since the typical period grows linearly with v. We use a variant of the method developed in [23]. Since the mapping of the driven ellipse is 4D and area-preserving, there are three possible types of eigenvalues of the corresponding real monodromy matrix: (i) Four complex eigenvalues (elliptic), (ii) two real and two complex eigenvalues (loxodromic) and (iii) four real eigenvalues (hyperbolic). The first case corresponds to stable [25], whereas the latter two represent unstable fixed points [21]. Performing extensive numerical computations we detected the periodic orbits up to period 100. The density ρ(F ) of the periodic orbits as a function of F (φ, α) is shown in Fig. 1. By using F , we can effectively map the φ × α plane onto a single dimension (see the insets of Fig. 1 where the F = const. isolines are shown in black). Fig. 1 shows Fig.  1a), the periodic orbits up to period 100 cover the whole available F -range, and thus the whole φ × α plane, except for a narrow region F ≈ 0 around the separatrix of the static system. This has to be expected, since the col-   2: (Color online) Collision resolved phase space density ρn 1 ,n 2 (F, ξ) (log colormap), n1 = 10, n2 = 10 2 (a) and n1 = 4 · 10 6 , n2 = 10 7 (b). ρn 1 ,n 2 (F, ξ) is largest around F ≈ −2.5 for small n (a), whereas the F (φ, α) ≈ 0 region is predominantly populated for high collision numbers (b). In the insets, the corresponding ρn 1 ,n 2 (φ, α) are shown. the billiard: The particles are fast enough to experience many collisions with the oscillating walls within a narrow interval of the phase of oscillation. The relative momentum change obeys δv/v ≪ 1 at each collision and the angle of incidence is almost the angle of reflection. As a result they trace orbits of the static system. As can be seen in Fig. 1b, the region F ≈ 0 possesses the highest density of period orbits, in particular unstable ones (hyperbolic and loxodromic), whereas the librator region around the two elliptic fixed points of the static system is depleted.
The consequences of the above described changes in the distributions of the periodic orbits with increasing v become immediately clear when studying the evolution (in terms of the number of collisions) of the phase space den- where the indices stand for the ith particle at collision number n, respectively. An ensemble of N p = 10 3 particles, with v 0 = 0.1 and ξ 0 , φ 0 , α 0 distributed uniformly and randomly (ξ 0 , φ 0 ∈ [0, 2π], α 0 ∈ [0, π]), is iterated for 10 7 collisions. The visiting probability is shown in a log scale colormap in Fig. 2. In the beginning (n 1 = 10, n 2 = 100) of the evolution, the region around the elliptic fixed points of the static system possesses the highest visiting probability, see Fig. 2a and especially the inset, where the corresponding ρ n1,n2 (φ, α) is shown. For large collisions numbers (n 1 = 4 · 10 6 , n 2 = 10 7 ), see Fig. 2b and the inset, the ensemble is located predominantly in a region around the separatrix (F (φ, α) = 0) of the corresponding static system. In the insets, the half-space α > π is not shown, since ρ(φ, α) ≈ 0 there. The effect of the population of the region around the separatrix becomes even more obvious when considering the phase space visiting probability velocity resolved rather than collision resolved, ρ v1,v2 (F, ξ) = 1 n2 v2 v1 n ρ(F, ξ, v, n)dv. All collisions between n 1 = 1 and n 2 = 10 7 with v 1 < v < v 2 are projected onto the F × ξ space. For low velocities, v 1 = 0, v 2 = 1 (see Fig. 3a), the complete F × ξ space and thus the complete φ × α plane (see the inset where the corresponding ρ v1,v2 (φ, α) is shown) is populated. For higher velocities, v 1 = 30, v 2 = 40 (see Fig. 3b and the inset), a narrow, sharply bounded region around F (φ, α) = 0 is now exclusively populated, i.e., outside this region no collisional events occur. This kind of population of phase space is generic and does not depend on the initial ensemble, as long as the initial velocity v 0 of the ensemble is sufficiently small. Unlike in dissipative systems the above squeezing of the ensemble is not caused by the existence of attractors but is due to a mechanism which will be explored in the following.
For low v we encounter a mixed phase space, where an infinite hierarchy of regular islands is embedded into a large chaotic sea. With increasing v, large regular regions emanate (we confirmed the existence of these regions among others by calculating the corresponding Lyapunov exponents and determined their extension) from the position of the elliptic fixed points of the static system. Their extension in the φ× α subspace grows rapidly with increasing v until v ≈ 10. From there on the regular regions capture a substantial amount of the available phase space volume, except a small region around F (φ, α) ≈ 0. If we now increase v further, the growth of these regular regions is reduced significantly, leaving the channel F (φ, α) ≈ 0 open for diffusive processes. The precondition to use sufficiently low v 0 in order to observe FA becomes obvious now: It ensures that particles start inside the chaotic sea which is connected to the channels F (φ, α) ≈ 0 that eventually allow for the acceleration process.
The evolution of the phase space density can then be understood as follows. The pronounced localization of the ensemble with respect to the F × ξ and correspondingly to the φ × α space in the v resolved picture (Fig. 3) is smeared out when considering the collision resolved density (Fig. 2). This happens because the increase of the ensemble average v with increasing n comes along with a certain distribution ρ(v) (a Maxwell-Boltzmann like distribution [22]), i.e., the n resolved distribution is a weighted superposition of the v resolved ρ v1,v2 (F, ξ), where the weights are given by ρ(v).
The evolution of the composition of phase space with increasing v has amazing consequences on the mean velocity |v| of the ensemble, in particular on its long time behavior. The rich phase space structure, associated with many spacious regular islands embedded into the low v chaotic sea, causes stickiness of the trajectories. As a result the dynamics is intermittent, exhibiting both laminar phases of oscillations as well as chaotic phases. The chaotic phases of the trajectories lead to stochastic fluctuations of v and thus to a net increase of the mean velocity |v| of the ensemble. This intermittent interplay of laminar and turbulent phases can be effectively described through continuous time random walks [24] predicting the appearance of subdiffusion, which is exactly observed here. With increasing |v| , the ensemble is more and more squeezed into the acceleration channels around F = 0. A detailed analysis based on the periodic orbits shows that there are much fewer and at the same time smaller regular islands in the acceleration channels for higher compared to low velocities, see also Fig. 1. This leads to a significantly reduced stickiness, i.e., fewer and shorter laminar phases once the ensemble is predominantly located in the narrow acceleration channels. Consequently, the diffusion in v space becomes then normal. The crossover at around n c ≈ 10 8 from subdiffusion to normal diffusion can be seen in Fig. 4 for two different values of the amplitude (C = 0.1 and C = 0.5). Such a crossover occurs in other driving modes of the ellipse as well, e.g. in the quadrupole mode, where a(t) = a 0 + C sin ωt and b(t) = b 0 − C sin ωt [22]. In Fig. 4, the average velocity grows according to a power law |v| (n) ∼ n β(C) between n ≈ 5 · 10 4 and n ≈ 10 8 . In this region, the diffusion exponent β(C) is amplitude dependent (β(0.1) = 0.28, β(0.5) = 0.40) and always smaller than 0.5 (subdiffusion). The dependence of β(C) in this regime can also be traced back to the phase space properties, e.g. there are more islands causing stickiness for smaller than for larger amplitudes [22]. For n > 10 8 |v| (n) ∼ n 1/2 , i.e., the exponent obeys β = 1/2 independent of the amplitude. Note however that an ensemble of particles has to be propagated for n ≫ n c ≈ 10 8 (until n = 10 11 in our case), which is extremely demanding and requires large scale computations, to detect this crossover. We checked the convergence of the ensemble with 160 particles shown in Fig. 4 with a reference ensemble consisting of 10 4 particles and shorter iteration times and found perfect agreement.
Our analysis of the 4D phase space of the driven elliptical billiard shows how its composition dynamically changes with increasing velocity v. For high v, thin acceleration channels of diffusive motion evolve from a spacious chaotic sea located at low v. This induces a crossover from sub to normal diffusion in momentum space. The acceleration channels are present in particular in the high v quasistatic limit. We expect this to be a generic feature of many driven dynamical systems which show Fermi acceleration. With increasing v, the ensemble will be focused on phase space regions where the corresponding static system possesses the highest degree of instability, or even chaoticity in the case of non-integrable billiards. In this sense driven billiards, as shown in the case of elliptical geometry, may operate not only as an accelerator but also as a phase space collimator.
Financial support by the Deutsche Forschungsgemeinschaft is gratefully acknowledged. F.L. acknowledges support from the Landesgraduiertenförderung Baden-Württemberg in the framework of the Heidelberg Graduate School of Fundamental Physics (HGSFP). F.K.D. likes to thank the HGSFP for financial support.