• Explore more than 100 years of PNAS content!
  • Sign-up for PNAS eTOC Alerts

Equilibration of energy in slow–fast systems

  1. Vered Rom-Kedare,1
  1. aDepartment of Electrical Engineering and Computer Science, Indian Institute of Science Education and Research, Bhopal 462066, India;
  2. bDepartment of Mathematics, Imperial College, London SW7 2AZ, United Kingdom;
  3. cLobachevsky University of Nizhny, Novgorod 603950, Russia;
  4. dMathematics Institute, University of Warwick, Coventry CV4 7AL, United Kingdom;
  5. eDepartment of Computer Science and Applied Mathematics, The Weizmann Institute of Science, Rehovot 76100, Israel
  1. Edited by Jean-Pierre Eckmann, University of Geneva, Geneva, Switzerland, and accepted by Editorial Board Member Herbert Levine October 29, 2017 (received for review April 17, 2017)

Significance

Do partial energies in slow–fast Hamiltonian systems equilibrate? This is a long-standing problem related to the foundation of statistical mechanics. Altering the traditional ergodic assumption, we propose that nonergodicity in the fast subsystem leads to equilibration of the whole system. To show this principle, we introduce a set of mechanical toy models—the springy billiards—and describe stochastic processes corresponding to their adiabatic behavior. We expect that these models and this principle will play an important role in the quest to establish and study the underlying postulates of statistical mechanics, one of the long-standing scientific grails.

Abstract

Ergodicity is a fundamental requirement for a dynamical system to reach a state of statistical equilibrium. However, in systems with several characteristic timescales, the ergodicity of the fast subsystem impedes the equilibration of the whole system because of the presence of an adiabatic invariant. In this paper, we show that violation of ergodicity in the fast dynamics can drive the whole system to equilibrium. To show this principle, we investigate the dynamics of springy billiards, which are mechanical systems composed of a small particle bouncing elastically in a bounded domain, where one of the boundary walls has finite mass and is attached to a linear spring. Numerical simulations show that the springy billiard systems approach equilibrium at an exponential rate. However, in the limit of vanishing particle-to-wall mass ratio, the equilibration rates remain strictly positive only when the fast particle dynamics reveal two or more ergodic components for a range of wall positions. For this case, we show that the slow dynamics of the moving wall can be modeled by a random process. Numerical simulations of the corresponding springy billiards and their random models show equilibration with similar positive rates.

In classical statistical mechanics, one deals with systems with a large number of degrees of freedom, where the exact knowledge of the state of the system at every moment in time is impossible to obtain or is irrelevant. The state variables are, therefore, declared “microscopic,” and an attempt is made to describe the properties of “macroscopic” variables (certain functions of the microscopic state) in an ensemble of many systems similar to the given one. In other words, statistical mechanics examines averaging over the phase space of a dynamical system, typically a Hamiltonian one. There is no a priori way of choosing the probability distribution over which the averaging is performed: given a dynamical system, its evolution is different for different initial conditions, and the distribution of the initial conditions is not encoded in the system and can be arbitrary. However, as the simplest option, one can use the so-called microcanonical ensemble [i.e., assume that the initial conditions are uniformly distributed in the phase space at a fixed value of the system’s energy (1??4)]. After this choice of the ensemble is made, standard results of statistical mechanics are recovered. In particular, one obtains the equipartition theorem, which describes the distribution of energy among different degrees of freedom (5). For multiparticle systems, this theorem implies that the ensemble averages of the kinetic energy per degree of freedom are equal for all particles, and this is used to derive fundamental thermodynamic properties of the system.

The overwhelming success of statistical mechanics indicates that there must be some deeper rationale beyond the assumption of uniform distribution of initial conditions. Since many physical experiments measure time averages, a classical approach is to invoke the ergodic hypothesis, which implies that the time average of an observable coincides with its phase space average over the microcanonical ensemble at the corresponding energy level. Thus, the derivation of thermodynamic laws for deterministic multiparticle dynamics relies on establishing ergodicity of the corresponding Hamiltonian system.

The most critical problem is that Hamiltonian systems are usually not ergodic in a range of energy levels, even if the number of degrees of freedom is large. For example, consider the paradigmatic model of the Boltzmann gas of hard spheres. It is, most probably, ergodic (6?8), but replacing the instantaneous collisions of the spheres with mutual repulsion is expected to destroy the ergodicity, even for an arbitrarily steep repulsing potential (9?11). In general, the dynamics in a smooth potential consist of a “chaotic sea” and “stability islands” (12, 13). A stability island contains a set of quasiperiodic motions confined only to a portion of the energy level because of which the ergodicity is broken. If the islands occupy a noticeable portion of the phase space, use of the microcanonical ensemble for averaging is unfounded. In general, this issue remains unresolved.

In this paper, we propose a mechanism for the onset of apparent ergodicity and mixing in slow–fast Hamiltonian systems. It is not based on either the usual assumption of a large number of degrees of freedom or the inherent instability of dispersing geometries (such as the hard spheres models).

We consider an isolated system, in which certain degrees of freedom evolve slower than the rest. One may think of the slow variables as the system’s parameters that evolve in reaction to the dynamics of the fast variables. The system obtained by freezing the values of these parameters is called the fast subsystem. If it is ergodic for every value of the frozen parameters, the evolution of the slow variables over a long time interval is accurately described by the system averaged over the microcanonical ensemble in the fast-phase space [a rigorous formulation of this fact is given by the Anosov–Kasuga averaging theorem (14?16)]. This averaged system has an additional conserved quantity, the Gibbs volume entropy of the fast subsystem: the fast subsystem responds adiabatically to the slow change in the parameters, and adiabatic processes are known to keep the entropy constant (5). Therefore, the averaged system is not ergodic.

We see that, on the timescale of the adiabatic approximation, the ergodicity of the fast subsystem prevents the whole slow–fast system from behaving ergodically. A similar phenomenon is observed in the “notorious piston problem,” where the adiabatic compression law implies that the system of two ideal gases at different temperatures contained in a finite cylinder and separated by a heavy piston never comes to equilibrium, which seems to defy the second law of thermodynamics (17??20).

Note that the fast subsystem’s entropy is preserved only approximately, and therefore, although the averaged system is not ergodic, the whole original system may still be. Deriving equations for the energy transfer between degrees of freedom in such systems is a challenging task, with results being sensitive to the structure of the interaction terms and the choice of scaling limit (21?23).

In this paper, we consider a more natural situation, in which the fast subsystem is not always ergodic, and the partition of the fast subspace into ergodic components depends on the slow variables. For example, the fast subsystem can be ergodic for some values of the slow variables and have more than one ergodic component for other values. Then, the adiabatic invariants obtained by averaging over the fast dynamics are destroyed. We show in this setting that ensemble averages of slow observables converge exponentially to the vicinity of the values obtained by averaging over the microcanonical ensemble in the whole-phase space. This suggests, somewhat unexpectedly, that the nonergodic behavior of the fast degrees of freedom can lead to effective equilibration of the whole system!

To elucidate this principle, we consider a special class of slow–fast Hamiltonian systems, springy billiards. These models often allow for a detailed explicit analysis, and their direct numerical simulations do not require integration of stiff differential equations. Recall that a billiard is a dynamical system representing a point particle, which moves inertially inside a domain and is reflected elastically from the domain boundary. Particle dynamics depend on the billiard domain shape. Thus, the motion may be integrable, or it may be chaotic: weakly chaotic with a zero Lyapunov exponent, chaotic and ergodic with a positive Lyapunov exponent, or chaotic with stability islands (24). The barred rectangle (25), stadium (26), and mushroom (27) billiard domains have been extensively studied as paradigmatic models for the different types of chaotic behavior. We use these models to create slow–fast systems with the controlled properties of the fast chaotic dynamics (Fig. 1). We allow one of the billiard boundaries, hereafter called the bar, to move. We assume that both the particle and the bar are of finite mass and that the bar is attached to a linear spring (refs. 28??31 have related setups). The bar mass <mml:math><mml:mi>M</mml:mi></mml:math>M is taken to be much larger than the particle mass <mml:math><mml:mi>m</mml:mi></mml:math>m, thus ensuring that the bar typically moves much slower than the particle. The billiard obtained by retaining the bar in a fixed position is called the frozen billiard.

The case of an infinitely heavy bar has been studied under the name of Fermi acceleration (32??????????43). In this case, the motion of the bar is not affected by the particle, while the particle gains or loses speed at every collision with the moving bar. Two types of behavior have been found. If the frozen billiard is ergodic for all possible positions of the bar (ergodic accelerators; e.g., stadium), then the ensemble-averaged kinetic energy may grow at most quadratically with time (34, 40, 44). If the ergodicity of the frozen billiard is broken for a part of the period of the bar’s motion (a multicomponent accelerator; e.g., barred rectangle and mushroom), then the particle kinetic energy grows exponentially with time for almost every initial condition (40?42, 44, 45). Thus, the nonergodicity of the frozen billiard accelerates the energy transfer from the bar to the particle.

Here, we consider a different case by taking the mass of the bar as finite while keeping the mass ratio <mml:math><mml:mrow><mml:mi>m</mml:mi><mml:mo>/</mml:mo><mml:mi>M</mml:mi></mml:mrow></mml:math>m/M small. Then, the question of Fermi acceleration is replaced by the question of equilibration—does the energy transfer between the bar and the particle lead the system to statistical equilibrium? Do ensemble averages converge to averages over the microcanonical ensemble? If so, then by the equipartition theorem, the averaged kinetic energies of the bar and particle should eventually become equal.

Results from studies of Fermi acceleration suggest that the equilibration depends on the geometry of the frozen billiard (i.e., on the dynamics of the fast subsystem). In general, we divide slow–fast Hamiltonian systems into two groups. The first consists of systems with an ergodic fast subsystem (EFS) for almost all values of the slow variables (EFS systems). An example is given by the springy stadium (Fig. 1B). The second is comprised of systems for which the structure of the partition of the fast subspace into ergodic components bifurcates as the slow variables change [examples are the springy barred rectangle (SBR) (Fig. 1A) and springy mushroom (Fig. 1C)]. Such systems, with variable partition of the fast subspace (VFS), are hereafter called VFS systems.

All of our numerical simulations for nonzero values of <mml:math><mml:mrow><mml:mi>m</mml:mi><mml:mo>/</mml:mo><mml:mi>M</mml:mi></mml:mrow></mml:math>m/M show exponential equilibration, suggesting apparent ergodicity and even mixing in both EFS and VFS springy billiards. However, we find that there is a fundamental difference between these two classes: the relaxation time in EFS systems tends to infinity as <mml:math><mml:mrow><mml:mrow><mml:mi>m</mml:mi><mml:mo>/</mml:mo><mml:mpadded width="+1.7pt"><mml:mi>M</mml:mi></mml:mpadded></mml:mrow><mml:mo>→</mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>m/M→?0, whereas the relaxation time stays bounded in VFS systems.

To analyze this phenomenon, we derive a stochastic model for the slow bar motion. In this model, the force exerted on the bar by the particle is found by averaging over one of the fast subsystem’s ergodic components, while the probability of switching between the components is determined by the bar position and velocity. The resulting Markov process of hopping between the ergodic components of the fast subsystem leads to equilibration. Numerics shows that the bar motion is quite accurately represented by this model, and the equilibration rates remain nonzero and approach the Markov process rates as <mml:math><mml:mrow><mml:mrow><mml:mi>m</mml:mi><mml:mo>/</mml:mo><mml:mpadded width="+1.7pt"><mml:mi>M</mml:mi></mml:mpadded></mml:mrow><mml:mo>→</mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>m/M→?0.

The derivation of the proposed stochastic models is not specific to billiards. Similar Markov processes can be constructed for a general VFS system by a procedure analogous to that used to study the Fermi acceleration in a homogeneous potential (46). Under standard irreducibility and aperiodicity conditions, such processes should converge to a unique stationary measure (cf. refs. 22, 23, 47, and 48), and by uniqueness, it must correspond to the Liouville measure of the whole system. Therefore, we believe that the proposed apparent ergodization mechanism should be universally applicable.

Particle in a Springy Billiard

Consider a particle of mass <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>m</mml:mi></mml:mpadded><mml:mo>?</mml:mo><mml:mn>?1</mml:mn></mml:mrow></mml:math>m??1 in a <mml:math><mml:mi>D</mml:mi></mml:math>D-dimensional billiard. One of the billiard walls is a bar of mass <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>M</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mn>?1</mml:mn></mml:mrow></mml:math>M=?1, which oscillates vertically by virtue of being suspended from a spring of spring constant <mml:math><mml:mi>k</mml:mi></mml:math>k. At impact, the bar and the particle undergo an elastic collision, leading to the exchange of momentum and energy between them:<mml:math display="block"><mml:mrow><mml:mrow><mml:mrow><mml:msubsup><mml:mi>v</mml:mi><mml:mi>p</mml:mi><mml:mo>+</mml:mo></mml:msubsup><mml:mo>=</mml:mo><mml:mpadded width="+1.7pt"><mml:mfrac><mml:mrow><mml:mrow><mml:mn>2</mml:mn><mml:msubsup><mml:mi>v</mml:mi><mml:mi>b</mml:mi><mml:mo>?</mml:mo></mml:msubsup></mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:mrow><mml:mo>(</mml:mo><mml:mrow><mml:mn>1</mml:mn><mml:mo>?</mml:mo><mml:mi>m</mml:mi></mml:mrow><mml:mo>)</mml:mo></mml:mrow><mml:msubsup><mml:mi>v</mml:mi><mml:mi>p</mml:mi><mml:mo>?</mml:mo></mml:msubsup></mml:mrow></mml:mrow><mml:mrow><mml:mn>1</mml:mn><mml:mo>+</mml:mo><mml:mi>m</mml:mi></mml:mrow></mml:mfrac></mml:mpadded></mml:mrow><mml:mo>,</mml:mo><mml:mo>?</mml:mo><mml:mrow><mml:msubsup><mml:mi>v</mml:mi><mml:mi>b</mml:mi><mml:mo>+</mml:mo></mml:msubsup><mml:mo>=</mml:mo><mml:mfrac><mml:mrow><mml:mrow><mml:mrow><mml:mo>(</mml:mo><mml:mrow><mml:mn>1</mml:mn><mml:mo>?</mml:mo><mml:mi>m</mml:mi></mml:mrow><mml:mo>)</mml:mo></mml:mrow><mml:msubsup><mml:mi>v</mml:mi><mml:mi>b</mml:mi><mml:mo>?</mml:mo></mml:msubsup></mml:mrow><mml:mo>+</mml:mo><mml:mrow><mml:mn>2</mml:mn><mml:mi>m</mml:mi><mml:msubsup><mml:mi>v</mml:mi><mml:mi>p</mml:mi><mml:mo>?</mml:mo></mml:msubsup></mml:mrow></mml:mrow><mml:mrow><mml:mn>1</mml:mn><mml:mo>+</mml:mo><mml:mi>m</mml:mi></mml:mrow></mml:mfrac></mml:mrow></mml:mrow><mml:mo>,</mml:mo></mml:mrow></mml:math>vp+=2vb??(1?m)vp?1+m,?vb+=(1?m)vb?+2mvp?1+m,[1]where <mml:math><mml:msubsup><mml:mi>v</mml:mi><mml:mi>b</mml:mi><mml:mo>?</mml:mo></mml:msubsup></mml:math>vb?, <mml:math><mml:msubsup><mml:mi>v</mml:mi><mml:mi>p</mml:mi><mml:mo>?</mml:mo></mml:msubsup></mml:math>vp? and <mml:math><mml:msubsup><mml:mi>v</mml:mi><mml:mi>b</mml:mi><mml:mo>+</mml:mo></mml:msubsup></mml:math>vb+, <mml:math><mml:msubsup><mml:mi>v</mml:mi><mml:mi>p</mml:mi><mml:mo>+</mml:mo></mml:msubsup></mml:math>vp+ are the vertical velocities of the bar <mml:math><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>b</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:math>(b) and the particle <mml:math><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>p</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:math>(p) just before and just after the collision, respectively. The total energy of the system, <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>E</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub><mml:mo>+</mml:mo><mml:msub><mml:mi>E</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow></mml:mrow></mml:math>E=Ep+Eb, is preserved, whereas the particle kinetic energy <mml:math><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:math>Ep and the bar energy <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>E</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:msubsup><mml:mi>v</mml:mi><mml:mi>b</mml:mi><mml:mn>2</mml:mn></mml:msubsup><mml:mo>+</mml:mo><mml:mrow><mml:mrow><mml:mi>k</mml:mi><mml:msubsup><mml:mi>y</mml:mi><mml:mi>b</mml:mi><mml:mn>2</mml:mn></mml:msubsup></mml:mrow><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow></mml:mrow></mml:mrow></mml:math>Eb=vb2+kyb2/2 (<mml:math><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:math>yb is the bar position) change at impact. Note that, between collisions, the particle moves inertially, while the bar oscillates harmonically. Therefore, the numerical simulations do not require integration to determine the dynamics of the system between the impacts.

The “bar–particle” system is Hamiltonian and has <mml:math><mml:mrow><mml:mi>D</mml:mi><mml:mo>+</mml:mo><mml:mn>1</mml:mn></mml:mrow></mml:math>D+1 degrees of freedom. It is slow–fast, since the bar speed <mml:math><mml:mrow><mml:mrow><mml:mo stretchy="false">|</mml:mo><mml:msub><mml:mi>v</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">|</mml:mo></mml:mrow><mml:mo>≤</mml:mo><mml:msqrt><mml:mrow><mml:mn>2</mml:mn><mml:mi>E</mml:mi></mml:mrow></mml:msqrt></mml:mrow></mml:math>|vb|≤2E is of order 1, while the particle speed <mml:math><mml:mrow><mml:mrow><mml:mo stretchy="false">|</mml:mo><mml:msub><mml:mi>v</mml:mi><mml:mi>p</mml:mi></mml:msub><mml:mo stretchy="false">|</mml:mo></mml:mrow><mml:mo>=</mml:mo><mml:msqrt><mml:mrow><mml:mrow><mml:mn>2</mml:mn><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mrow><mml:mo>/</mml:mo><mml:mi>m</mml:mi></mml:mrow></mml:msqrt></mml:mrow></mml:math>|vp|=2Ep/m is typically large. As the particle moves fast, many collisions with the bar occur during each short time interval; the averaged motion of the bar is governed by the equation<mml:math display="block"><mml:mrow><mml:mrow><mml:mrow><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>¨</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi></mml:msub><mml:mo>+</mml:mo><mml:mfrac><mml:mrow><mml:mi>d</mml:mi><mml:mi>U</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mrow><mml:mi>d</mml:mi><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow></mml:mfrac></mml:mrow><mml:mo>=</mml:mo><mml:mi>F</mml:mi></mml:mrow><mml:mo>,</mml:mo></mml:mrow></mml:math>y¨b+dU(yb)dyb=F,[2]where <mml:math><mml:mrow><mml:mrow><mml:mi>U</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:mn>?1</mml:mn><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow><mml:mi>k</mml:mi><mml:msubsup><mml:mi>y</mml:mi><mml:mi>b</mml:mi><mml:mn>2</mml:mn></mml:msubsup></mml:mrow></mml:mrow></mml:math>U(yb)=?1/2kyb2 is the spring potential and <mml:math><mml:mi>F</mml:mi></mml:math>F denotes the averaged force exerted on the bar at position <mml:math><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:math>yb by the particle with energy<mml:math><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:math>Ep. The averaging is performed over many collisions at a frozen value of <mml:math><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:math>yb. Since the work done by this force corresponds to the change in particle energy, it follows that <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>F</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:mrow><mml:mrow><mml:mi>d</mml:mi><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mrow><mml:mo>/</mml:mo><mml:mi>d</mml:mi></mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow></mml:mrow></mml:mrow></mml:math>F=?dEp/dyb.

If the frozen billiard is ergodic for each value of <mml:math><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:math>yb, the Anosov–Kasuga theorem (20) implies that the phase–space volume at a given energy level is approximately preserved (for most trajectories, for sufficiently small <mml:math><mml:mi>m</mml:mi></mml:math>m, and any finite interval of the slow time). Hence,<mml:math display="block"><mml:mrow><mml:mrow><mml:mi>J</mml:mi><mml:mo>=</mml:mo><mml:mrow><mml:msubsup><mml:mi>E</mml:mi><mml:mi>p</mml:mi><mml:mfrac><mml:mi>D</mml:mi><mml:mn>2</mml:mn></mml:mfrac></mml:msubsup><mml:mi>V</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>≈</mml:mo><mml:mi>const</mml:mi></mml:mrow><mml:mo>,</mml:mo></mml:mrow></mml:math>J=EpD2V(yb)≈const,[3]where <mml:math><mml:mrow><mml:mi>V</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>V(yb) is the volume of the frozen billiard domain. Thus, in the adiabatic approximation,<mml:math display="block"><mml:mrow><mml:mrow><mml:mfrac><mml:mrow><mml:mi>d</mml:mi><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mrow><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mfrac><mml:mo>=</mml:mo><mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:mfrac><mml:mn>2</mml:mn><mml:mi>D</mml:mi></mml:mfrac><mml:mfrac><mml:mrow><mml:mi>d</mml:mi><mml:mi>V</mml:mi></mml:mrow><mml:mi>V</mml:mi></mml:mfrac></mml:mrow></mml:mrow></mml:mrow><mml:mo>.</mml:mo></mml:mrow></mml:math>dEpEp=?2DdVV.[4]This implies that the force acting on the bar is <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>F</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:mrow><mml:mrow><mml:mrow><mml:mrow><mml:mn>?2</mml:mn><mml:mo>/</mml:mo><mml:mpadded width="+1.7pt"><mml:mi>D</mml:mi></mml:mpadded></mml:mrow><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mrow><mml:mo>/</mml:mo><mml:mi>V</mml:mi></mml:mrow><mml:mi>d</mml:mi><mml:mi>V</mml:mi></mml:mrow><mml:mo>/</mml:mo><mml:mi>d</mml:mi></mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow></mml:mrow></mml:math>F=?2/DEp/VdV/dyb. Since <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:mi>E</mml:mi><mml:mo>?</mml:mo><mml:msub><mml:mi>E</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mi>E</mml:mi><mml:mo>?</mml:mo><mml:mrow><mml:mrow><mml:mn>1</mml:mn><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow><mml:msubsup><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi><mml:mn>2</mml:mn></mml:msubsup></mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:mi>U</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:mrow></mml:math>Ep=E?Eb=E?1/2y˙b2?U(yb), Eq. 2 becomes<mml:math display="block"><mml:mrow><mml:mrow><mml:mrow><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>¨</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi></mml:msub><mml:mo>+</mml:mo><mml:mfrac><mml:mrow><mml:mi>d</mml:mi><mml:mi>U</mml:mi><mml:mrow><mml:mo>(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo>)</mml:mo></mml:mrow></mml:mrow><mml:mrow><mml:mi>d</mml:mi><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow></mml:mfrac></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mfrac><mml:mn>2</mml:mn><mml:mi>D</mml:mi></mml:mfrac><mml:mrow><mml:mo>[</mml:mo><mml:mrow><mml:mi>E</mml:mi><mml:mo>?</mml:mo><mml:mrow><mml:mstyle displaystyle="false"><mml:mfrac><mml:mn>1</mml:mn><mml:mn>2</mml:mn></mml:mfrac></mml:mstyle><mml:msubsup><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi><mml:mn>2</mml:mn></mml:msubsup></mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:mi>U</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow><mml:mo>]</mml:mo></mml:mrow><mml:mfrac><mml:mrow><mml:mi>d</mml:mi><mml:mi>V</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mrow><mml:mi>V</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mi>d</mml:mi><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow></mml:mfrac></mml:mrow></mml:mrow><mml:mo>.</mml:mo></mml:mrow></mml:math>y¨b+dU(yb)dyb=2D[E?12y˙b2?U(yb)]dV(yb)V(yb)dyb.[5]Note that the same formula follows from the ideal gas law.

A straightforward computation shows that solutions of Eq. 5 indeed follow the level sets of <mml:math><mml:mi>J</mml:mi></mml:math>J. By noting that the adiabatic law given by Eq. 3 takes the form <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msup><mml:mrow><mml:mo>(</mml:mo><mml:mrow><mml:mrow><mml:mi>J</mml:mi><mml:mo>/</mml:mo><mml:mi>V</mml:mi></mml:mrow><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>)</mml:mo></mml:mrow><mml:mrow><mml:mn>2</mml:mn><mml:mo>/</mml:mo><mml:mi>D</mml:mi></mml:mrow></mml:msup></mml:mpadded><mml:mo>=</mml:mo><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mrow></mml:math>(J/V(yb))2/D=Ep, we find that the adiabatic bar motion is governed by an effective potential<mml:math display="block"><mml:mrow><mml:mrow><mml:mrow><mml:msub><mml:mi>U</mml:mi><mml:mtext>eff</mml:mtext></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:mi>U</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>+</mml:mo><mml:msup><mml:mrow><mml:mo>(</mml:mo><mml:mfrac><mml:mi>J</mml:mi><mml:mrow><mml:mi>V</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mfrac><mml:mo>)</mml:mo></mml:mrow><mml:mfrac><mml:mn>2</mml:mn><mml:mi>D</mml:mi></mml:mfrac></mml:msup></mml:mrow></mml:mrow><mml:mo>,</mml:mo></mml:mrow></mml:math>Ueff(yb)=U(yb)+(JV(yb))2D,[6]where <mml:math><mml:mi>J</mml:mi></mml:math>J is determined by the initial condition. A similar effective potential was derived in the context of the piston problem (19, 20). It is easy to see that <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mn>0</mml:mn></mml:mpadded><mml:mo>≤</mml:mo><mml:mpadded width="+1.7pt"><mml:mi>J</mml:mi></mml:mpadded><mml:mo>≤</mml:mo><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>J</mml:mi><mml:mi>f</mml:mi></mml:msub></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:msup><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mi>E</mml:mi><mml:mo>?</mml:mo><mml:mrow><mml:mi>U</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>f</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mrow><mml:mn>2</mml:mn><mml:mo>/</mml:mo><mml:mi>D</mml:mi></mml:mrow></mml:msup><mml:mi>V</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>f</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>0≤J≤Jf=(E?U(yf))2/DV(yf). Here, <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mpadded><mml:mo>=</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>f</mml:mi></mml:msub></mml:mrow></mml:math>yb=yf corresponds to the pressure equilibrium, where the spring compensates for the pressure on the bar caused by the collisions with the particle (<mml:math><mml:msub><mml:mi>y</mml:mi><mml:mi>f</mml:mi></mml:msub></mml:math>yf depends on the total energy <mml:math><mml:mi>E</mml:mi></mml:math>E). At the other extreme (i.e., <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>J</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>J=?0), the particle does not move, and all of the energy is in the oscillating bar.

We see that, when the frozen billiard is ergodic for every position of the bar (an EFS system), the adiabatic approximation predicts that the bar motion will be periodic in time. Consequently, the springy billiard does not equilibrate on the timescale for which the adiabatic invariant is preserved. An example of such a system is a springy stadium (Fig. 1B), which is obtained by a modification of the famous Bunimovich stadium. We will discuss the process of equilibration in this billiard later in the paper.

The averaging theory can sometimes be extended to systems with a nonergodic fast subsystem (for example, refs. 41 and 44). Such extension requires detailed knowledge of the partition of the fast subspace into ergodic components. In the general theory of Hamiltonian systems, this partition is not known. However, mathematical billiards provide a rich set of examples, for which the ergodic decomposition can be described in full detail. We use this knowledge to show the accelerated equilibration for two different examples of springy billiards with nonergodic fast subspace (VFS systems). The first is the SBR, and the second is a modification of the Bunimovich mushroom (Fig. 1 A and C).

Springy Barred Rectangle

Consider a particle moving in a rectangle of length <mml:math><mml:mi>L</mml:mi></mml:math>L and height <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>h</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mn>?2</mml:mn></mml:mrow></mml:math>h=?2, which is partially split by a horizontal bar of length <mml:math><mml:mi>λ</mml:mi></mml:math>λ attached to a spring (Fig. 1A). The particle’s horizontal speed <mml:math><mml:mrow><mml:mo stretchy="false">|</mml:mo><mml:msub><mml:mi>u</mml:mi><mml:mi>p</mml:mi></mml:msub><mml:mo stretchy="false">|</mml:mo></mml:mrow></mml:math>|up| is preserved, and thus, the horizontal motion is periodic with period <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>T</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:mn>?2</mml:mn><mml:mi>L</mml:mi></mml:mrow><mml:mo>/</mml:mo><mml:mrow><mml:mo stretchy="false">|</mml:mo><mml:msub><mml:mi>u</mml:mi><mml:mi>p</mml:mi></mml:msub><mml:mo stretchy="false">|</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>Tp=?2L/|up|. The particle moves with a large vertical speed and moderate horizontal speed, and therefore, the fast subsystem corresponds to the particle’s vertical motion, while its horizontal position and the bar position are frozen. The period <mml:math><mml:msub><mml:mi>T</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:math>Tp is divided into two time intervals. In the first, the particle does not hit the bar, and its vertical speed is preserved; the fast subsystem, therefore, has a single ergodic component corresponding to the particle moving in the segment <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>y</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mpadded><mml:mo>∈</mml:mo><mml:mrow><mml:mo stretchy="false">[</mml:mo><mml:mrow><mml:mo>?</mml:mo><mml:mn>1</mml:mn></mml:mrow><mml:mo>,</mml:mo><mml:mn>1</mml:mn><mml:mo stretchy="false">]</mml:mo></mml:mrow></mml:mrow></mml:math>yp∈[?1,1]. In the second interval of length <mml:math><mml:mrow><mml:mi>τ</mml:mi><mml:msub><mml:mi>T</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mrow></mml:math>τTp, where <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>τ</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:mi>λ</mml:mi><mml:mo>/</mml:mo><mml:mi>L</mml:mi></mml:mrow></mml:mrow></mml:math>τ=λ/L, the particle enters a chamber above or below the bar, where it gains or loses vertical speed as it hits the moving bar many times. Here, the fast subsystem has two ergodic components, one corresponding to the motion above the bar, <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>y</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mpadded><mml:mo>∈</mml:mo><mml:mrow><mml:mo stretchy="false">[</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo>,</mml:mo><mml:mn>1</mml:mn><mml:mo stretchy="false">]</mml:mo></mml:mrow></mml:mrow></mml:math>yp∈[yb,1], and the other corresponding to the motion below the bar, <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>y</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mpadded><mml:mo>∈</mml:mo><mml:mrow><mml:mo stretchy="false">[</mml:mo><mml:mrow><mml:mo>?</mml:mo><mml:mn>1</mml:mn></mml:mrow><mml:mo>,</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">]</mml:mo></mml:mrow></mml:mrow></mml:math>yp∈[?1,yb]. During a single passage above or below the bar, the vertical speed <mml:math><mml:mrow><mml:mo stretchy="false">|</mml:mo><mml:msub><mml:mi>v</mml:mi><mml:mi>p</mml:mi></mml:msub><mml:mo stretchy="false">|</mml:mo></mml:mrow></mml:math>|vp| obeys the adiabatic law given by Eq. 3 with <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>D</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mn>?1</mml:mn></mml:mrow></mml:math>D=?1, <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:mn>?1</mml:mn><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow><mml:mi>m</mml:mi><mml:mpadded width="+1.7pt"><mml:msubsup><mml:mi>v</mml:mi><mml:mi>p</mml:mi><mml:mn>2</mml:mn></mml:msubsup></mml:mpadded></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mi>E</mml:mi><mml:mo>?</mml:mo><mml:mrow><mml:mrow><mml:mn>1</mml:mn><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow><mml:msubsup><mml:mi>v</mml:mi><mml:mi>b</mml:mi><mml:mn>2</mml:mn></mml:msubsup></mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:mrow><mml:mn>1</mml:mn><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow><mml:mi>k</mml:mi><mml:msubsup><mml:mi>y</mml:mi><mml:mi>b</mml:mi><mml:mn>2</mml:mn></mml:msubsup></mml:mrow></mml:mrow></mml:mrow></mml:math>Ep=?1/2mvp2=E?1/2vb2?1/2kyb2, and <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>V</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:msup><mml:mi>V</mml:mi><mml:mrow><mml:mi>u</mml:mi><mml:mi>p</mml:mi></mml:mrow></mml:msup><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mn>?1</mml:mn><mml:mo>?</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow></mml:mrow></mml:math>V=Vup(yb)=?1?yb when the particle is above the bar or <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>V</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:msup><mml:mi>V</mml:mi><mml:mrow><mml:mi>d</mml:mi><mml:mi>o</mml:mi><mml:mi>w</mml:mi><mml:mi>n</mml:mi></mml:mrow></mml:msup><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mn>?1</mml:mn><mml:mo>+</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow></mml:mrow></mml:math>V=Vdown(yb)=?1+yb when the particle is below it. Hence, <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>J</mml:mi><mml:mrow><mml:mrow><mml:mrow><mml:mi>u</mml:mi><mml:mi>p</mml:mi></mml:mrow><mml:mo>/</mml:mo><mml:mi>d</mml:mi></mml:mrow><mml:mi>o</mml:mi><mml:mi>w</mml:mi><mml:mi>n</mml:mi></mml:mrow></mml:msub></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:msqrt><mml:mrow><mml:mi>E</mml:mi><mml:mo>?</mml:mo><mml:mrow><mml:mrow><mml:mn>1</mml:mn><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow><mml:msubsup><mml:mi>v</mml:mi><mml:mi>b</mml:mi><mml:mn>2</mml:mn></mml:msubsup></mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:mrow><mml:mn>1</mml:mn><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow><mml:mi>k</mml:mi><mml:msubsup><mml:mi>y</mml:mi><mml:mi>b</mml:mi><mml:mn>2</mml:mn></mml:msubsup></mml:mrow></mml:mrow></mml:msqrt><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mn>1</mml:mn><mml:mo>?</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mi>const</mml:mi></mml:mrow></mml:math>Jup/down=E?1/2vb2?1/2kyb2(1?yb)=const. Since <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:msubsup><mml:mi>J</mml:mi><mml:mrow><mml:mrow><mml:mrow><mml:mi>u</mml:mi><mml:mi>p</mml:mi></mml:mrow><mml:mo>/</mml:mo><mml:mi>d</mml:mi></mml:mrow><mml:mi>o</mml:mi><mml:mi>w</mml:mi><mml:mi>n</mml:mi></mml:mrow><mml:mn>2</mml:mn></mml:msubsup><mml:mo>/</mml:mo><mml:msup><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mn>1</mml:mn><mml:mo>?</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mn>2</mml:mn></mml:msup></mml:mrow></mml:mrow></mml:math>Ep=Jup/down2/(1?yb)2, the right-hand side of Eq. 2 becomes <mml:math><mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:mrow><mml:mn>2</mml:mn><mml:msubsup><mml:mi>J</mml:mi><mml:mrow><mml:mrow><mml:mrow><mml:mi>u</mml:mi><mml:mi>p</mml:mi></mml:mrow><mml:mo>/</mml:mo><mml:mi>d</mml:mi></mml:mrow><mml:mi>o</mml:mi><mml:mi>w</mml:mi><mml:mi>n</mml:mi></mml:mrow><mml:mn>2</mml:mn></mml:msubsup></mml:mrow><mml:mo>/</mml:mo><mml:msup><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mn>1</mml:mn><mml:mo>?</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mn>3</mml:mn></mml:msup></mml:mrow></mml:mrow></mml:math>?2Jup/down2/(1?yb)3 for the particle above/below the bar, respectively. Note that, here, <mml:math><mml:mi>E</mml:mi></mml:math>E and <mml:math><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:math>Ep do not include the particle’s horizontal kinetic energy, which is decoupled from the dynamics.

To construct effective equations for the bar’s average motion, we approximate the deterministic process by a stochastic one, assuming that the probability of entering the chamber above or below the bar is proportional to the width of the corresponding gap between the bar and the boundary of the rectangle. These probabilities are equal to <mml:math><mml:mrow><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mn>1</mml:mn><mml:mo>?</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow></mml:math>(1?yb)/2, where <mml:math><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:math>yb is taken at the moment of entrance into the chamber. The same assumptions were successfully used in ref. 45 for the study of Fermi acceleration in a rectangle with an infinitely heavy oscillating bar.

Hence, we suggest that the bar–particle system is well-approximated by the following <mml:math><mml:msub><mml:mi>T</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:math>Tp-periodic probabilistic hybrid model:<mml:math display="block"><mml:mrow><mml:mrow><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>¨</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi></mml:msub><mml:mo>+</mml:mo><mml:mrow><mml:mi>k</mml:mi><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mo>{</mml:mo><mml:mtable displaystyle="true"><mml:mtr><mml:mtd columnalign="center"><mml:mrow><mml:mrow><mml:msub><mml:mi>F</mml:mi><mml:mi>j</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo>,</mml:mo><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>,</mml:mo></mml:mrow></mml:mtd><mml:mtd columnalign="left"><mml:mrow><mml:mrow><mml:mrow><mml:mtext>for</mml:mtext><mml:mo>?</mml:mo><mml:mi>j</mml:mi><mml:msub><mml:mi>T</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mrow><mml:mo>≤</mml:mo><mml:mi>t</mml:mi><mml:mo><</mml:mo><mml:mrow><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mi>j</mml:mi><mml:mo>+</mml:mo><mml:mi>τ</mml:mi></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:msub><mml:mi>T</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mrow></mml:mrow><mml:mo>,</mml:mo></mml:mrow></mml:mtd></mml:mtr><mml:mtr><mml:mtd columnalign="center"><mml:mrow><mml:mn>0</mml:mn><mml:mo>,</mml:mo></mml:mrow></mml:mtd><mml:mtd columnalign="left"><mml:mrow><mml:mrow><mml:mrow><mml:mtext>for</mml:mtext><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mi>j</mml:mi><mml:mo>+</mml:mo><mml:mi>τ</mml:mi></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:msub><mml:mi>T</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mrow><mml:mo>≤</mml:mo><mml:mi>t</mml:mi><mml:mo><</mml:mo><mml:mrow><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mi>j</mml:mi><mml:mo>+</mml:mo><mml:mn>1</mml:mn></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:msub><mml:mi>T</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mrow></mml:mrow><mml:mo>,</mml:mo></mml:mrow></mml:mtd></mml:mtr></mml:mtable></mml:mrow></mml:mrow></mml:math>y¨b+kyb={Fj(yb,y˙b),for?jTp≤t<(j+τ)Tp,0,for(j+τ)Tp≤t<(j+1)Tp,[7]where <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>j</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:mi>t</mml:mi><mml:mo>/</mml:mo><mml:msub><mml:mi>T</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mrow><mml:mo>?</mml:mo></mml:mrow></mml:mrow></mml:math>j=?t/Tp? is the number of the period. For each period, a choice is made between <mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>F</mml:mi><mml:mi>j</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:mrow><mml:mn>2</mml:mn><mml:msubsup><mml:mi>J</mml:mi><mml:mrow><mml:mi>j</mml:mi><mml:mo>,</mml:mo><mml:mrow><mml:mi>u</mml:mi><mml:mi>p</mml:mi></mml:mrow></mml:mrow><mml:mn>2</mml:mn></mml:msubsup></mml:mrow><mml:mo>/</mml:mo><mml:msup><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mn>1</mml:mn><mml:mo>?</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mn>3</mml:mn></mml:msup></mml:mrow></mml:mrow></mml:mrow></mml:math>Fj(yb)=?2Jj,up2/(1?yb)3 (with probability <mml:math><mml:msub><mml:mi>β</mml:mi><mml:mi>j</mml:mi></mml:msub></mml:math>βj) and <mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>F</mml:mi><mml:mi>j</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:mn>?2</mml:mn><mml:msubsup><mml:mi>J</mml:mi><mml:mrow><mml:mi>j</mml:mi><mml:mo>,</mml:mo><mml:mrow><mml:mi>d</mml:mi><mml:mi>o</mml:mi><mml:mi>w</mml:mi><mml:mi>n</mml:mi></mml:mrow></mml:mrow><mml:mn>2</mml:mn></mml:msubsup></mml:mrow><mml:mo>/</mml:mo><mml:msup><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mn>1</mml:mn><mml:mo>+</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mn>3</mml:mn></mml:msup></mml:mrow></mml:mrow></mml:math>Fj(yb)=?2Jj,down2/(1+yb)3 (with probability <mml:math><mml:mrow><mml:mn>1</mml:mn><mml:mo>?</mml:mo><mml:msub><mml:mi>β</mml:mi><mml:mi>j</mml:mi></mml:msub></mml:mrow></mml:math>1?βj). Here, <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>J</mml:mi><mml:mrow><mml:mi>j</mml:mi><mml:mo>,</mml:mo><mml:mrow><mml:mrow><mml:mrow><mml:mi>u</mml:mi><mml:mi>p</mml:mi></mml:mrow><mml:mo>/</mml:mo><mml:mi>d</mml:mi></mml:mrow><mml:mi>o</mml:mi><mml:mi>w</mml:mi><mml:mi>n</mml:mi></mml:mrow></mml:mrow></mml:msub></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:msqrt><mml:mrow><mml:mi>E</mml:mi><mml:mo>?</mml:mo><mml:msub><mml:mi>E</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow></mml:msqrt><mml:mrow><mml:mo>(</mml:mo><mml:mrow><mml:mn>1</mml:mn><mml:mo>?</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow><mml:mo>)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>Jj,up/down=E?Eb(1?yb), <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>E</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:mo>(</mml:mo><mml:mrow><mml:msubsup><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi><mml:mn>2</mml:mn></mml:msubsup><mml:mo>+</mml:mo><mml:mrow><mml:mi>k</mml:mi><mml:msubsup><mml:mi>y</mml:mi><mml:mi>b</mml:mi><mml:mn>2</mml:mn></mml:msubsup></mml:mrow></mml:mrow><mml:mo>)</mml:mo></mml:mrow><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow></mml:mrow></mml:math>Eb=(y˙b2+kyb2)/2, and <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>β</mml:mi><mml:mi>j</mml:mi></mml:msub></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mn>1</mml:mn><mml:mo>?</mml:mo><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mi>j</mml:mi><mml:msub><mml:mi>T</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow></mml:mrow></mml:math>βj=(1?yb(jTp))/2.

The bar dynamics follow the level lines of <mml:math><mml:msub><mml:mi>J</mml:mi><mml:mrow><mml:mi>u</mml:mi><mml:mi>p</mml:mi></mml:mrow></mml:msub></mml:math>Jup, <mml:math><mml:msub><mml:mi>J</mml:mi><mml:mrow><mml:mi>d</mml:mi><mml:mi>o</mml:mi><mml:mi>w</mml:mi><mml:mi>n</mml:mi></mml:mrow></mml:msub></mml:math>Jdown, or <mml:math><mml:mrow><mml:msub><mml:mi>J</mml:mi><mml:mn>0</mml:mn></mml:msub><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:mn>1</mml:mn><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:msubsup><mml:mi>v</mml:mi><mml:mi>b</mml:mi><mml:mn>2</mml:mn></mml:msubsup><mml:mo>+</mml:mo><mml:mrow><mml:mi>k</mml:mi><mml:msubsup><mml:mi>y</mml:mi><mml:mi>b</mml:mi><mml:mn>2</mml:mn></mml:msubsup></mml:mrow></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>J0=1/2(vb2+kyb2). At <mml:math><mml:mrow><mml:mi>t</mml:mi><mml:mo>=</mml:mo><mml:mrow><mml:mi>j</mml:mi><mml:msub><mml:mi>T</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mrow></mml:mrow></mml:math>t=jTp, the bar height <mml:math><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mi>j</mml:mi><mml:msub><mml:mi>T</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>yb(jTp) determines the probability <mml:math><mml:msub><mml:mi>β</mml:mi><mml:mi>j</mml:mi></mml:msub></mml:math>βj of choosing <mml:math><mml:msub><mml:mi>J</mml:mi><mml:mrow><mml:mi>u</mml:mi><mml:mi>p</mml:mi></mml:mrow></mml:msub></mml:math>Jup vs. <mml:math><mml:msub><mml:mi>J</mml:mi><mml:mrow><mml:mi>d</mml:mi><mml:mi>o</mml:mi><mml:mi>w</mml:mi><mml:mi>n</mml:mi></mml:mrow></mml:msub></mml:math>Jdown, and <mml:math><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mi>j</mml:mi><mml:msub><mml:mi>T</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>,</mml:mo><mml:mrow><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mi>j</mml:mi><mml:msub><mml:mi>T</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:math>(yb(jTp),y˙b(jTp)) determines the particular <mml:math><mml:msub><mml:mi>J</mml:mi><mml:mrow><mml:mrow><mml:mrow><mml:mi>u</mml:mi><mml:mi>p</mml:mi></mml:mrow><mml:mo>/</mml:mo><mml:mi>d</mml:mi></mml:mrow><mml:mi>o</mml:mi><mml:mi>w</mml:mi><mml:mi>n</mml:mi></mml:mrow></mml:msub></mml:math>Jup/down level set along which the motion will continue.

In Fig. 2A, we plot the projection of direct simulations of the bar–particle system onto the <mml:math><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo>,</mml:mo><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:math>(yb,y˙b) plane for a few bar oscillations. Fig. 2A shows that the bar motion indeed follows the level sets of the corresponding <mml:math><mml:mi>J</mml:mi></mml:math>J values and that the trajectory switches between level sets (Fig. 2A, Upper Left Inset); while the bar energy changes significantly when the particle is above/below the bar (Fig. 2A, Lower Right Inset), the corresponding <mml:math><mml:mi>J</mml:mi></mml:math>J values remain approximately constant for each such interval (Fig. 2A, Lower Left Inset).

Fig. 2.

Piecewise adiabatic bar motion for the (A) SBR and (B) springy mushroom. Two trajectories of the springy billiards with particle mass <mml:math><mml:mrow><mml:mi>m</mml:mi><mml:mo>=</mml:mo><mml:msup><mml:mn>10</mml:mn><mml:mrow><mml:mo>?</mml:mo><mml:mn>8</mml:mn></mml:mrow></mml:msup></mml:mrow></mml:math>m=10?8, one with <mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mn>0</mml:mn><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mn>0.9</mml:mn></mml:mrow></mml:math>Ep(0)=0.9 (at the center) and one with <mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mn>0</mml:mn><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mn>0.1</mml:mn></mml:mrow></mml:math>Ep(0)=0.1 (at the periphery), are projected to the bar-phase plane, <mml:math><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mtext>,</mml:mtext><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi></mml:msub></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:math>(yb,y˙b), and plotted (black lines) on top of level lines of the adiabatic invariant. Green indicates the free springy bar motion. Red and blue indicate level lines for the bar motion in the effective potential created when the particle hits the bar from above and below, respectively. Each orbit closely follows the level lines, switching occasionally between them. Lower Left Insets show that the corresponding adiabatic invariants are approximately piecewise constant, and Lower Right Insets show the bar energy oscillations. In A, Upper Left Inset shows a magnification of the transition from a level set of one averaged law of motion to another.

To study the equilibration process, we conduct the following numerical experiment. We generate an ensemble of initial conditions and for each initial condition, simulate a single-particle motion in the springy billiard. At each moment in time, we compute the ensemble-averaged difference between the kinetic energy of the bar and the vertical kinetic energy of the particle, <mml:math><mml:mrow><mml:mrow><mml:mi mathvariant="normal">Δ</mml:mi><mml:mi>K</mml:mi><mml:mpadded width="+1.7pt"><mml:mi>E</mml:mi></mml:mpadded></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:mrow><mml:mi>M</mml:mi><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow><mml:msubsup><mml:mi>v</mml:mi><mml:mi>b</mml:mi><mml:mn>2</mml:mn></mml:msubsup></mml:mrow><mml:mo>?</mml:mo></mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:mrow><mml:mi>m</mml:mi><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow><mml:msubsup><mml:mi>v</mml:mi><mml:mi>p</mml:mi><mml:mn>2</mml:mn></mml:msubsup></mml:mrow><mml:mo>?</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>ΔKE=?M/2vb2???m/2vp2?. We study the time evolution of <mml:math><mml:mrow><mml:mi mathvariant="normal">Δ</mml:mi><mml:mi>K</mml:mi><mml:mi>E</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>t</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>ΔKE(t) for ensembles of <mml:math><mml:msup><mml:mn>10</mml:mn><mml:mn>4</mml:mn></mml:msup></mml:math>104 initial states with a fixed initial bar energy. The initial phase of the bar oscillation is chosen randomly (with the uniform distribution). The position of the particle is also chosen randomly inside the rectangle. Finally, the vertical speed of the particle is determined from the initial particle energy, while the direction of motion (up or down) is chosen randomly. We repeat the experiment for different values of particle mass <mml:math><mml:mi>m</mml:mi></mml:math>m and different initial values of the bar energy, keeping the total energy of the system constant, <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>E</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mn>?1</mml:mn></mml:mrow></mml:math>E=?1.

In all of our simulations, we find that <mml:math><mml:mrow><mml:mrow><mml:mi mathvariant="normal">Δ</mml:mi><mml:mi>K</mml:mi><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>E</mml:mi><mml:mrow><mml:mi>t</mml:mi><mml:mo>→</mml:mo><mml:mrow><mml:mo>+</mml:mo><mml:mi mathvariant="normal">∞</mml:mi></mml:mrow></mml:mrow></mml:msub></mml:mpadded></mml:mrow><mml:mo>→</mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>ΔKEt→+∞→?0 exponentially (Fig. 3A). However, for each fixed <mml:math><mml:mi>m</mml:mi></mml:math>m, the rate of convergence depends on the choice of the initial ensemble and to some extent, the choice of fitting interval; similar phenomena for different setups were reported in refs. 49 and 50. To enable a comparison of the rates of convergence in different numerical experiments, we fix a practical definition of the equilibration rate as the best fitted slope to <mml:math><mml:mrow><mml:mi>log</mml:mi><mml:mrow><mml:mo>|</mml:mo><mml:mrow><mml:mi mathvariant="normal">Δ</mml:mi><mml:mi>K</mml:mi><mml:mi>E</mml:mi></mml:mrow><mml:mo>|</mml:mo></mml:mrow></mml:mrow></mml:math>log|ΔKE| vs. <mml:math><mml:mi>t</mml:mi></mml:math>t over a time interval <mml:math><mml:mrow><mml:mo stretchy="false">[</mml:mo><mml:mn>0</mml:mn><mml:mo>,</mml:mo><mml:msub><mml:mi>T</mml:mi><mml:mn>0</mml:mn></mml:msub><mml:mo stretchy="false">]</mml:mo></mml:mrow></mml:math>[0,T0], where <mml:math><mml:msub><mml:mi>T</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:math>T0 is defined by <mml:math><mml:mrow><mml:mrow><mml:mi>log</mml:mi><mml:mrow><mml:mo stretchy="false">|</mml:mo><mml:mrow><mml:mrow><mml:mrow><mml:mi mathvariant="normal">Δ</mml:mi><mml:mi>K</mml:mi><mml:mi>E</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mn>0</mml:mn><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>/</mml:mo><mml:mi mathvariant="normal">Δ</mml:mi></mml:mrow><mml:mi>K</mml:mi><mml:mi>E</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>T</mml:mi><mml:mn>0</mml:mn></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo stretchy="false">|</mml:mo></mml:mrow></mml:mrow><mml:mo>≈</mml:mo><mml:mn>?1</mml:mn></mml:mrow></mml:math>log|ΔKE(0)/ΔKE(T0)|≈?1. The details of the computations are described in Numerical Algorithms.

Fig. 3.

Kinetic energy equilibration for springy billiards. For two ensembles of <mml:math><mml:msup><mml:mn>10</mml:mn><mml:mn>4</mml:mn></mml:msup></mml:math>104 initial states, one with <mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>E</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mn>0</mml:mn><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mn>0.9</mml:mn></mml:mrow></mml:math>Eb(0)=0.9 and the other with <mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>E</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mn>0</mml:mn><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mn>0.1</mml:mn></mml:mrow></mml:math>Eb(0)=0.1, the ensemble average, <mml:math><mml:mrow><mml:mi mathvariant="normal">Δ</mml:mi><mml:mi>K</mml:mi><mml:mi>E</mml:mi></mml:mrow></mml:math>ΔKE, of the difference between the kinetic energies per degree of freedom for the bar and the particle is plotted as a function of time. (A) The SBR equilibration process is essentially the same for the mass ratios <mml:math><mml:mrow><mml:mrow><mml:mi>m</mml:mi><mml:mo>/</mml:mo><mml:mi>M</mml:mi></mml:mrow><mml:mo>=</mml:mo><mml:msup><mml:mn>10</mml:mn><mml:mrow><mml:mo>?</mml:mo><mml:mn>5</mml:mn></mml:mrow></mml:msup></mml:mrow></mml:math>m/M=10?5, <mml:math><mml:msup><mml:mn>10</mml:mn><mml:mrow><mml:mo>?</mml:mo><mml:mn>7</mml:mn></mml:mrow></mml:msup></mml:math>10?7 and for the random model of Eq. 7. (B) The stadium and mushroom equilibration process at <mml:math><mml:mrow><mml:mrow><mml:mi>m</mml:mi><mml:mo>/</mml:mo><mml:mi>M</mml:mi></mml:mrow><mml:mo>=</mml:mo><mml:msup><mml:mn>10</mml:mn><mml:mrow><mml:mo>?</mml:mo><mml:mn>5</mml:mn></mml:mrow></mml:msup></mml:mrow></mml:math>m/M=10?5. The springy stadium equilibrates more slowly than the springy mushroom. Insets show a semilog plot, which shows the exponential decay of <mml:math><mml:mrow><mml:mi mathvariant="normal">Δ</mml:mi><mml:mi>K</mml:mi><mml:mi>E</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>t</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>ΔKE(t).

We examine how the rates depend on the mass ratio. For sufficiently small <mml:math><mml:mi>m</mml:mi></mml:math>m, the rates do not display a significant dependence on <mml:math><mml:mi>m</mml:mi></mml:math>m (Fig. 3A) and do not approach zero. To validate the result, we compare the equilibration rates of the direct bar–particle simulations with those of the simulations of the theoretical stochastic model described by Eq. 7. As can be seen in Fig. 4A, we obtain very similar rates; 10 runs of the <mml:math><mml:msup><mml:mn>10</mml:mn><mml:mn>4</mml:mn></mml:msup></mml:math>104 particle ensembles in the stochastic model lead to the rates <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mn>0.111</mml:mn></mml:mpadded><mml:mo>±</mml:mo><mml:mn>?0.002</mml:mn></mml:mrow></mml:math>0.111±?0.002 and <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mn>0.126</mml:mn></mml:mpadded><mml:mo>±</mml:mo><mml:mn>?0.004</mml:mn></mml:mrow></mml:math>0.126±?0.004 for <mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>E</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mn>0</mml:mn><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mn>?0.1</mml:mn></mml:mrow></mml:math>Eb(0)=?0.1 and <mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>E</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mn>0</mml:mn><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mn>?0.9</mml:mn></mml:mrow></mml:math>Eb(0)=?0.9, respectively. Similarly, averaging over 100 runs of direct simulations of the SBR, 10 runs for each of the 10 different <mml:math><mml:mi>m</mml:mi></mml:math>m values, produces the rates <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mn>0.112</mml:mn></mml:mpadded><mml:mo>±</mml:mo><mml:mn>?0.002</mml:mn></mml:mrow></mml:math>0.112±?0.002 and <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mn>0.128</mml:mn></mml:mpadded><mml:mo>±</mml:mo><mml:mn>?0.003</mml:mn></mml:mrow></mml:math>0.128±?0.003, respectively. Since the theoretical model has no adjustable parameters, the agreement suggests that the SBR equilibrates, with the relaxation time remaining bounded as <mml:math><mml:mrow><mml:mrow><mml:mi>m</mml:mi><mml:mo>/</mml:mo><mml:mpadded width="+1.7pt"><mml:mi>M</mml:mi></mml:mpadded></mml:mrow><mml:mo>→</mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>m/M→?0.

Fig. 4.

Energy equilibration rate dependence on the mass ratio. Each data point corresponds to the decay rate for a fixed value of <mml:math><mml:msqrt><mml:mrow><mml:mi>m</mml:mi><mml:mo>/</mml:mo><mml:mi>M</mml:mi></mml:mrow></mml:msqrt></mml:math>m/M and a fixed value of <mml:math><mml:mrow><mml:msub><mml:mi>E</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mn>0</mml:mn><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>Eb(0) of the ensemble-averaged difference <mml:math><mml:mrow><mml:mi mathvariant="normal">Δ</mml:mi><mml:mi>K</mml:mi><mml:mi>E</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>t</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>ΔKE(t) between the kinetic energies per degree of freedom for the bar and particle. (A) In the SBR case, there are 10 red [<mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>E</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mn>0</mml:mn><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mn>?0.1</mml:mn></mml:mrow></mml:math>Eb(0)=?0.1] and 10 blue [<mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>E</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mn>0</mml:mn><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mn>?0.9</mml:mn></mml:mrow></mml:math>Eb(0)=?0.9] dots for each value of <mml:math><mml:mrow><mml:mi>m</mml:mi><mml:mo>/</mml:mo><mml:mi>M</mml:mi></mml:mrow></mml:math>m/M. The dots on the y axis are the rates obtained from simulations of the stochastic model (Eq. 7). The solid blue/red lines show the mean rates of the stochastic model; the dashed lines bound the bands of three standard deviations. We see that 98% of the rates for <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>m</mml:mi></mml:mpadded><mml:mo>></mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>m>?0 are within these bands. (B) In the springy stadium case, the rates decrease to zero as <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>m</mml:mi></mml:mpadded><mml:mo>→</mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>m→?0, whereas for the mushroom, they decrease to positive values close to the rates predicted by simulations of the random model of Eq. 9. The open rectangle and open circle on the y axis in Inset represent the equilibration rates <mml:math><mml:mrow><mml:mn>1.42</mml:mn><mml:mo>×</mml:mo><mml:msup><mml:mn>10</mml:mn><mml:mrow><mml:mo>?</mml:mo><mml:mn>4</mml:mn></mml:mrow></mml:msup></mml:mrow></mml:math>1.42×10?4 and <mml:math><mml:mrow><mml:mn>2.71</mml:mn><mml:mo>×</mml:mo><mml:msup><mml:mn>10</mml:mn><mml:mrow><mml:mo>?</mml:mo><mml:mn>4</mml:mn></mml:mrow></mml:msup></mml:mrow></mml:math>2.71×10?4 obtained from the stochastic model for <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>E</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mpadded><mml:mo>=</mml:mo><mml:mn>?0.9</mml:mn></mml:mrow></mml:math>Eb=?0.9 and <mml:math><mml:mn>0.1</mml:mn></mml:math>0.1, respectively.

The SBR can serve as a prototype for equilibration mechanisms in a much more general class of VFS systems. For example, separation of the fast space into two ergodic components can be achieved by splitting a chaotic billiard domain into two separate spatial regions (in the Fermi acceleration case, such systems were studied in refs. 40 and 44). However, the SBR is a rather special model, because the dynamics in each frozen ergodic component are integrable. Moreover, transitions between the fast ergodic components only happen at prescribed moments of time. One could potentially use these two properties to rigorously establish the ergodicity of the SBR using hyperbolic theory techniques, as this problem has features similar to linked twist maps and their generalizations described in refs. 51??54. However, one does not expect a general VFS system to have these properties.

To show that the equilibration observed in the SBR is not a peculiarity of that particular model, in the following, we consider another class of springy billiard systems, which have a different and more general structure of the decomposition of the fast subspace into ergodic components. For this class, we again show that violation of ergodicity in the fast-phase space leads to energy equilibration in the whole slow–fast system.

Springy Mushroom

When there are several fast degrees of freedom, the fast subspace of a slow–fast Hamiltonian system typically includes both regular and chaotic components (12, 32). As the slow variables change, a trajectory can transfer between these components. To study the influence of these transitions on equilibration, we introduce the springy mushroom model.

A mushroom billiard consists of a circular cap and a rectangular stem. It was proposed in ref. 27 to serve as a paradigm for systems with a mixed phase space, as it has the simplest possible mixed-phase space structure: there is a unique chaotic component and one regular (and completely integrable) component—a stability island that corresponds to the particle trapped forever in the cap.

In this paper, we use a variation of the original Bunimovich mushroom (Fig. 1C). Our mushroom billiard is bounded by a semicircular cap, two slanted side walls, and a straight bar at the bottom. This shape is a version of the famous Bunimovich stadium (26); the slanted walls prevent the appearance of horizontal parabolic orbits, which could spoil the statistics (cf. ref. 55). The corresponding billiard has a single chaotic component, which occupies the whole phase space. To create an integrable island, while keeping the billiard area unchanged, we add two horizontal straight segments, which leave a throat of radius <mml:math><mml:mi>w</mml:mi></mml:math>w at the cap diameter. In this case, the trajectories that bounce inside the cap close to the circular boundary remain forever in the cap, and their motion is integrable (the absolute value of the angular momentum is preserved for such orbits). The size of this island depends on <mml:math><mml:mi>w</mml:mi></mml:math>w, and the island vanishes when <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>w</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mi>r</mml:mi></mml:mrow></mml:math>w=r, where <mml:math><mml:mi>r</mml:mi></mml:math>r is the radius of the cap. However, the trajectories that move through the throat between the cap and the stem form a single chaotic component in the phase space (26, 27).

In the springy mushroom, the bar located at the bottom of the stem can move vertically. It has mass <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>M</mml:mi></mml:mpadded><mml:mo>></mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>M>?0 and is attached to a linear spring (Fig. 1C). The bar moves harmonically between collisions with the particle, whereas the phase and amplitude of the oscillations change at every collision.

Finally, we take the radius of the throat change to be a prescribed function of the bar position: <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>w</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:mi>w</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>w=w(yb) (e.g., we use Eq. 11). Then, the island size depends on the bar position <mml:math><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:math>yb. At a moment of time when the particle is in the stem, it belongs to the chaotic component of the frozen mushroom. After some collisions with the walls and the bar, the particle typically enters the cap. If the throat shrinks while the particle is in the cap, the particle may be captured in the cap, and it will stay there until the throat opens wide enough to let the particle back into the stem. During the period of capture, the particle follows the integrable dynamics of a semicircular billiard. Since the particle motion is fast compared with the bar, we can interpret these capture and release events as transitions between the chaotic and regular components of the frozen mushroom.

To construct effective equations for the average bar motion, we need to find the effective pressure exerted by the particle on the bar. Since the frozen billiard is not ergodic, the standard adiabatic law given by Eq. 4 is not applicable. Instead, we use the “leaky adiabatic law” of ref. 41:<mml:math display="block"><mml:mrow><mml:mrow><mml:mfrac><mml:mrow><mml:mi>d</mml:mi><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mrow><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mfrac><mml:mo>=</mml:mo><mml:mrow><mml:mo>?</mml:mo><mml:mfrac><mml:mrow><mml:mi>d</mml:mi><mml:mi>V</mml:mi></mml:mrow><mml:msub><mml:mi>V</mml:mi><mml:mi>c</mml:mi></mml:msub></mml:mfrac></mml:mrow></mml:mrow><mml:mo>,</mml:mo></mml:mrow></mml:math>dEpEp=?dVVc,[8]where <mml:math><mml:mrow><mml:mi>V</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>V(yb) is the phase space volume of the energy level <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:mn>?1</mml:mn><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow></mml:mrow></mml:math>Ep=?1/2 in the phase space of the frozen billiard and <mml:math><mml:mrow><mml:msub><mml:mi>V</mml:mi><mml:mi>c</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>Vc(yb) is the phase space volume of the chaotic component for the same energy level. For billiards, <mml:math><mml:mi>V</mml:mi></mml:math>V is the area of the billiard domain times <mml:math><mml:mrow><mml:mn>2</mml:mn><mml:mi>π</mml:mi></mml:mrow></mml:math> (the velocity circle). The effect of collisions between the particle and the bar is averaged over the chaotic component only. To account for this, <mml:math><mml:mi>V</mml:mi></mml:math>V is replaced by <mml:math><mml:msub><mml:mi>V</mml:mi><mml:mi>c</mml:mi></mml:msub></mml:math>Vc in the denominator of Eq. 4 (with <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>D</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mn>?2</mml:mn></mml:mrow></mml:math>D=?2). The term <mml:math><mml:mrow><mml:mi>d</mml:mi><mml:mi>V</mml:mi></mml:mrow></mml:math>dV is kept in the numerator, because the work is done only by the bar and is proportional to the change in volume of the entire billiard table.

Expressions for the volumes in Eq. 8 can be found explicitly. Indeed, at the energy level <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:mn>?1</mml:mn><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow></mml:mrow></mml:math>Ep=?1/2 for the frozen billiard, we get <mml:math><mml:mrow><mml:mrow><mml:mi>V</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:msup><mml:mi>π</mml:mi><mml:mn>2</mml:mn></mml:msup><mml:mpadded width="+1.7pt"><mml:msup><mml:mi>r</mml:mi><mml:mn>2</mml:mn></mml:msup></mml:mpadded></mml:mrow><mml:mo>+</mml:mo><mml:mrow><mml:mn>?2</mml:mn><mml:mi>π</mml:mi><mml:mrow><mml:mo>(</mml:mo><mml:mrow><mml:mrow><mml:mn>2</mml:mn><mml:mi>r</mml:mi><mml:mi mathvariant="normal">?</mml:mi></mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:msup><mml:mi mathvariant="normal">?</mml:mi><mml:mn>2</mml:mn></mml:msup><mml:mrow><mml:mi>tan</mml:mi><mml:mo>?</mml:mo><mml:mi>θ</mml:mi></mml:mrow></mml:mrow></mml:mrow><mml:mo>)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:mrow></mml:math>V(yb)=π2r2+?2π(2r???2tan?θ), where <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi mathvariant="normal">?</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi mathvariant="normal">?</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:mpadded><mml:mo>?</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow></mml:mrow></mml:math>?=?0?yb is the stem length. The chaotic component is the complement of the integrable one. Thus, <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>V</mml:mi><mml:mi>c</mml:mi></mml:msub></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>V</mml:mi></mml:mpadded><mml:mo>?</mml:mo><mml:msub><mml:mi>V</mml:mi><mml:mtext>isl</mml:mtext></mml:msub></mml:mrow></mml:mrow></mml:math>Vc=V?Visl, where <mml:math><mml:msub><mml:mi>V</mml:mi><mml:mtext>isl</mml:mtext></mml:msub></mml:math>Visl is the volume of the integrable island. It can be found by using the conservation of the angular momentum, <mml:math><mml:mrow><mml:mrow><mml:mrow><mml:msub><mml:mi>V</mml:mi><mml:mtext>isl</mml:mtext></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mn>?2</mml:mn><mml:mi>π</mml:mi><mml:msup><mml:mi>r</mml:mi><mml:mn>2</mml:mn></mml:msup><mml:mrow><mml:mo>(</mml:mo><mml:mrow><mml:mrow><mml:msup><mml:mi>cos</mml:mi><mml:mrow><mml:mo>?</mml:mo><mml:mn>1</mml:mn></mml:mrow></mml:msup><mml:mrow><mml:mi>w</mml:mi><mml:mo>/</mml:mo><mml:mi>r</mml:mi></mml:mrow></mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:mrow><mml:mi>w</mml:mi><mml:mo>/</mml:mo><mml:mi>r</mml:mi></mml:mrow><mml:msqrt><mml:mrow><mml:mn>1</mml:mn><mml:mo>?</mml:mo><mml:msup><mml:mrow><mml:mo>(</mml:mo><mml:mrow><mml:mi>w</mml:mi><mml:mo>/</mml:mo><mml:mi>r</mml:mi></mml:mrow><mml:mo>)</mml:mo></mml:mrow><mml:mn>2</mml:mn></mml:msup></mml:mrow></mml:msqrt></mml:mrow></mml:mrow><mml:mo>)</mml:mo></mml:mrow></mml:mrow></mml:mrow><mml:mo>,</mml:mo></mml:mrow></mml:math>Visl(yb)=?2πr2(cos?1w/r?w/r1?(w/r)2), where <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>w</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:mi>w</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>w=w(yb) is the width of the cap throat. Ref. 41 has a detailed derivation of these formulas.

When the particle is in the chaotic zone, it acts on the bar with average force <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>F</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:mrow><mml:mrow><mml:mi>d</mml:mi><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mrow><mml:mo>/</mml:mo><mml:mi>d</mml:mi></mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow></mml:mrow></mml:mrow></mml:math>F=?dEp/dyb. However, when the particle is in the cap, it does not collide with the bar. Using Eq. 8, we conclude that the averaged equations for the bar motion are of the form<mml:math display="block"><mml:mrow><mml:mrow><mml:mrow><mml:mover accent="true"><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo>¨</mml:mo></mml:mover></mml:mrow><mml:mo>+</mml:mo><mml:mfrac><mml:mrow><mml:mi>d</mml:mi><mml:mi>U</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mrow><mml:mi>d</mml:mi><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow></mml:mfrac></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mo>{</mml:mo><mml:mtable displaystyle="true"><mml:mtr><mml:mtd columnalign="left"><mml:mrow><mml:mpadded width="+1.7pt"><mml:mfrac><mml:mrow><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo>,</mml:mo><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mrow><mml:msub><mml:mi>V</mml:mi><mml:mi>c</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mfrac></mml:mpadded><mml:mfrac><mml:mrow><mml:mi>d</mml:mi><mml:mi>V</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mrow><mml:mi>d</mml:mi><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow></mml:mfrac></mml:mrow></mml:mtd><mml:mtd columnalign="left"><mml:mrow><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mtext>particle in</mml:mtext><mml:msub><mml:mi>V</mml:mi><mml:mi>c</mml:mi></mml:msub></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mo>,</mml:mo></mml:mrow></mml:mtd></mml:mtr><mml:mtr><mml:mtd columnalign="left"><mml:mn>0</mml:mn></mml:mtd><mml:mtd columnalign="left"><mml:mrow><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mtext>particle in</mml:mtext><mml:msub><mml:mi>V</mml:mi><mml:mtext>isl</mml:mtext></mml:msub></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mo>,</mml:mo></mml:mrow></mml:mtd></mml:mtr></mml:mtable></mml:mrow></mml:mrow></mml:math>yb¨+dU(yb)dyb={Ep(yb,y˙b)Vc(yb)dV(yb)dyb(particle inVc),0(particle inVisl),[9]where we use energy conservation to express the particle energy <mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo>,</mml:mo><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>E</mml:mi></mml:mpadded><mml:mo>?</mml:mo><mml:mrow><mml:mrow><mml:mn>?1</mml:mn><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow><mml:mpadded width="+1.7pt"><mml:msubsup><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi><mml:mn>2</mml:mn></mml:msubsup></mml:mpadded></mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:mi>U</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:mrow></mml:math>Ep(yb,y˙b)=E??1/2y˙b2?U(yb) as a function of <mml:math><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:math>yb and <mml:math><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi></mml:msub></mml:math>y˙b. The first of the equations conserves <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>J</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo>,</mml:mo><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mi>G</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>J=Ep(yb,y˙b)G(yb), where <mml:math><mml:mrow><mml:mrow><mml:mi>G</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mi>exp</mml:mi><mml:mrow><mml:mo>(</mml:mo><mml:mrow><mml:msubsup><mml:mo largeop="true" symmetric="true">∫</mml:mo><mml:mn>0</mml:mn><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:msubsup><mml:mrow><mml:mrow><mml:mrow><mml:mrow><mml:mrow><mml:mi>d</mml:mi><mml:mi>V</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>s</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>/</mml:mo><mml:mi>d</mml:mi></mml:mrow><mml:mi>s</mml:mi><mml:mi>d</mml:mi><mml:mi>s</mml:mi></mml:mrow><mml:mo>/</mml:mo><mml:msub><mml:mi>V</mml:mi><mml:mi>c</mml:mi></mml:msub></mml:mrow><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>s</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow><mml:mo>)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>G(yb)=exp(∫0ybdV(s)/dsds/Vc(s)). Thus, the fast collisions of the particle with the bar create, on average, an effective potential: <mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>U</mml:mi><mml:mtext>eff</mml:mtext></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo>;</mml:mo><mml:msub><mml:mi>J</mml:mi><mml:mn>0</mml:mn></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:mi>U</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>+</mml:mo><mml:mrow><mml:mrow><mml:msub><mml:mi>J</mml:mi><mml:mn>0</mml:mn></mml:msub><mml:mo>/</mml:mo><mml:mi>G</mml:mi></mml:mrow><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:mrow></mml:math>Ueff(yb;J0)=U(yb)+J0/G(yb), where <mml:math><mml:msub><mml:mi>J</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:math>J0 is the value of <mml:math><mml:mi>J</mml:mi></mml:math>J at the moment the particle enters the chaotic zone. While the particle is in the chaotic component, the bar motion occurs along the level line <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>J</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:msub><mml:mi>J</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:mrow></mml:math>J=J0, which is also a level line of the Hamiltonian defined by the effective potential <mml:math><mml:mrow><mml:msub><mml:mi>U</mml:mi><mml:mtext>??????</mml:mtext></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo>;</mml:mo><mml:msub><mml:mi>J</mml:mi><mml:mn>0</mml:mn></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>U??????(yb;J0). Similarly, when the particle is in the integrable component, the bar motion occurs along the level lines of the Hamiltonian defined by the free spring potential <mml:math><mml:mrow><mml:mi>U</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>U(yb).

We conclude that, similar to the SBR case, the bar’s averaged motion switches between level sets of different functions. Numerical simulations of the bar–particle dynamics confirmed this property for small mass ratios: Fig. 2B shows that trajectories of the springy mushroom follow a level set of <mml:math><mml:mi>J</mml:mi></mml:math>J (Fig. 2B, red circles) and then, by following a level set of the free motion (Fig. 2B, green circles), switch to follow another level set of <mml:math><mml:mi>J</mml:mi></mml:math>J. Fig. 2B, Lower Left Inset shows that <mml:math><mml:mi>J</mml:mi></mml:math>J is approximately constant during the time intervals in which the particle remains in the chaotic zone, and Fig. 2B, Lower Right Inset shows the corresponding oscillations of <mml:math><mml:msub><mml:mi>E</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:math>Eb.

To complete the derivation of the stochastic model for the motion of the bar, we describe the transition rules between the two possible states of the particle. A transition from the chaotic component to the integrable one can be modeled by a random process as follows. Suppose that the particle is in the chaotic component at time <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>t</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>c</mml:mi></mml:msub></mml:mrow></mml:math>t=tc. Then, the motion is ergodic, and the probability of transfer during a time interval <mml:math><mml:mrow><mml:mo>(</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>c</mml:mi></mml:msub><mml:mo>,</mml:mo><mml:mrow><mml:msub><mml:mi>t</mml:mi><mml:mi>c</mml:mi></mml:msub><mml:mo>+</mml:mo><mml:mrow><mml:mi>d</mml:mi><mml:mi>t</mml:mi></mml:mrow></mml:mrow><mml:mo>)</mml:mo></mml:mrow></mml:math>(tc,tc+dt) coincides with the ratio of the transferred phase space volume to the chaotic zone volume. The transferred volume is positive only when the throat is shrinking (i.e., when <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mrow><mml:mover accent="true"><mml:mi>w</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow></mml:mpadded><mml:mo><</mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>w˙<?0 at <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>t</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>c</mml:mi></mml:msub></mml:mrow></mml:math>t=tc). In this case, the probability is given by<mml:math display="block"><mml:mrow><mml:mrow><mml:mrow><mml:msub><mml:mi>P</mml:mi><mml:mrow><mml:mtext>cha</mml:mtext><mml:mo>→</mml:mo><mml:mtext>isl</mml:mtext></mml:mrow></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>c</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mfrac><mml:mrow><mml:msubsup><mml:mi>V</mml:mi><mml:mtext>isl</mml:mtext><mml:mo>′</mml:mo></mml:msubsup><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi></mml:msub></mml:mrow><mml:mrow><mml:msub><mml:mi>V</mml:mi><mml:mi>c</mml:mi></mml:msub><mml:mrow><mml:mo>(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo>)</mml:mo></mml:mrow></mml:mrow></mml:mfrac><mml:mi>d</mml:mi><mml:mi>t</mml:mi></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mi>d</mml:mi><mml:mrow><mml:mi>ln</mml:mi><mml:mfrac><mml:mrow><mml:mi>G</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mrow><mml:msub><mml:mi>V</mml:mi><mml:mi>c</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mfrac></mml:mrow></mml:mrow></mml:mrow><mml:mo>.</mml:mo></mml:mrow></mml:math>Pcha→isl(tc)=Visl′(yb)y˙bVc(yb)dt=dlnG(yb)Vc(yb).[10]If <mml:math><mml:mrow><mml:mrow><mml:mover accent="true"><mml:mi>w</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow><mml:mo>></mml:mo><mml:mn>0</mml:mn></mml:mrow></mml:math>w˙>0, then <mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>P</mml:mi><mml:mrow><mml:mtext>cha</mml:mtext><mml:mo>→</mml:mo><mml:mtext>isl</mml:mtext></mml:mrow></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>c</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mn>0</mml:mn></mml:mrow></mml:math>Pcha→isl(tc)=0.

The transition from the integrable component back to the chaotic one is described by the following deterministic rule (41). Suppose that the particle is captured in the integrable component at time <mml:math><mml:msub><mml:mi>t</mml:mi><mml:mi>c</mml:mi></mml:msub></mml:math>tc. While the particle is captured, the bar moves harmonically. Therefore, the bar position and the throat size change periodically with time. Let <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>t</mml:mi><mml:mi>r</mml:mi></mml:msub></mml:mpadded><mml:mo>></mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>c</mml:mi></mml:msub></mml:mrow></mml:math>tr>tc be the first instance of time, such that <mml:math><mml:mrow><mml:mrow><mml:mi>w</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>r</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mi>w</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>c</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>w(yb(tr))=w(yb(tc)). Then, <mml:math><mml:mrow><mml:mrow><mml:mi>w</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>t</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>></mml:mo><mml:mrow><mml:mi>w</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>r</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>w(yb(t))>w(yb(tr)) for <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>t</mml:mi><mml:mi>c</mml:mi></mml:msub></mml:mpadded><mml:mo><</mml:mo><mml:mpadded width="+1.7pt"><mml:mi>t</mml:mi></mml:mpadded><mml:mo><</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>r</mml:mi></mml:msub></mml:mrow></mml:math>tc<t<tr. Since typically, the particle trajectory in the integrable zone densely fills a strip between the semicircle of radius <mml:math><mml:mrow><mml:mi>w</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>c</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>w(yb(tc)) and the cap boundary, the particle is released from the cap at a time interval <mml:math><mml:mrow><mml:mo>(</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>r</mml:mi></mml:msub><mml:mo>,</mml:mo><mml:mrow><mml:msub><mml:mi>t</mml:mi><mml:mi>r</mml:mi></mml:msub><mml:mo>+</mml:mo><mml:mrow><mml:mi>d</mml:mi><mml:mi>t</mml:mi></mml:mrow></mml:mrow><mml:mo>)</mml:mo></mml:mrow></mml:math>(tr,tr+dt). In the adiabatic limit, the indeterminacy in the exit moment, <mml:math><mml:mrow><mml:mi>d</mml:mi><mml:mi>t</mml:mi></mml:mrow></mml:math>dt, vanishes.

These arguments provide the transition rules for the stochastic model of Eq. 9: when the throat shrinks, the particle can move from the chaotic component to the integrable one at time <mml:math><mml:msub><mml:mi>t</mml:mi><mml:mi>c</mml:mi></mml:msub></mml:math>tc with probability <mml:math><mml:mrow><mml:msub><mml:mi>P</mml:mi><mml:mrow><mml:mtext>cha</mml:mtext><mml:mo>→</mml:mo><mml:mtext>isl</mml:mtext></mml:mrow></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>c</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>Pcha→isl(tc) given by Eq. 10. If captured, the particle is released at time <mml:math><mml:msub><mml:mi>t</mml:mi><mml:mi>r</mml:mi></mml:msub></mml:math>tr when the throat restores its size [i.e., when <mml:math><mml:mrow><mml:mrow><mml:mi>w</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>r</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mi>w</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>c</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>w(yb(tr))=w(yb(tc)) for the first time].

Now we are ready to describe the mechanism for the ergodization in the springy mushroom. Choosing a nonmonotonic function <mml:math><mml:mrow><mml:mi>w</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>w(yb) allows for the possibility that the bar positions at the moments of capture and release are different: <mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>c</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>≠</mml:mo><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>r</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>yb(tc)≠yb(tr). Then, the monotonicity of the function <mml:math><mml:mrow><mml:mi>G</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>G(yb) implies that the particle returns to the chaotic zone with a different value of the adiabatic invariant <mml:math><mml:mrow><mml:mrow><mml:mi>J</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>r</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>,</mml:mo><mml:mrow><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>r</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>≠</mml:mo><mml:mrow><mml:mi>J</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>c</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>,</mml:mo><mml:mrow><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>c</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>J(yb(tr),y˙b(tr))≠J(yb(tc),y˙b(tc)). Consequently, this capture–release event moves the particle to a new level line of the adiabatic invariant. The destruction of the adiabatic invariant enhances the energy transfer between the fast and slow degrees of freedom, and therefore, repeated capture–release events lead to faster equilibration.

For example, in our numerical experiments, we use the nonmonotonic function<mml:math display="block"><mml:mrow><mml:mrow><mml:mi>w</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mo>{</mml:mo><mml:mtable displaystyle="true"><mml:mtr><mml:mtd columnalign="left"><mml:mrow><mml:mrow><mml:mi>min</mml:mi><mml:mrow><mml:mo>{</mml:mo><mml:mn>1</mml:mn><mml:mo>,</mml:mo><mml:mrow><mml:mn>?0.7</mml:mn><mml:mo>+</mml:mo><mml:mrow><mml:mn>0.6</mml:mn><mml:msup><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo>?</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mn>2</mml:mn></mml:msup></mml:mrow></mml:mrow><mml:mo>}</mml:mo></mml:mrow></mml:mrow><mml:mo>?</mml:mo></mml:mrow></mml:mtd><mml:mtd columnalign="left"><mml:mrow><mml:mrow><mml:mrow><mml:mtext>for</mml:mtext><mml:mo>?</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow><mml:mo><</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:mrow><mml:mo>,</mml:mo></mml:mrow></mml:mtd></mml:mtr><mml:mtr><mml:mtd columnalign="left"><mml:mrow><mml:mrow><mml:mi>min</mml:mi><mml:mrow><mml:mo>{</mml:mo><mml:mn>1</mml:mn><mml:mo>,</mml:mo><mml:mrow><mml:mn>?0.7</mml:mn><mml:mo>+</mml:mo><mml:mrow><mml:mn>6.0</mml:mn><mml:msup><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo>?</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mn>2</mml:mn></mml:msup></mml:mrow></mml:mrow><mml:mo>}</mml:mo></mml:mrow></mml:mrow><mml:mo>?</mml:mo></mml:mrow></mml:mtd><mml:mtd columnalign="left"><mml:mrow><mml:mrow><mml:mrow><mml:mtext>for</mml:mtext><mml:mo>?</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow><mml:mo>≥</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:mrow><mml:mo>,</mml:mo></mml:mrow></mml:mtd></mml:mtr></mml:mtable></mml:mrow></mml:mrow></mml:math>w(yb)={min{1,?0.7+0.6(yb?y0)2}?for?yb<y0,min{1,?0.7+6.0(yb?y0)2}?for?yb≥y0,[11]where <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>y</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:mo>?</mml:mo><mml:mn>0.2436</mml:mn></mml:mrow></mml:mrow></mml:math>y0=?0.2436. This choice of <mml:math><mml:mi>w</mml:mi></mml:math>w eliminates some evident obstacles for ergodicity in the springy mushroom and its random model. In particular, we have <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mn>0</mml:mn></mml:mpadded><mml:mo>≤</mml:mo><mml:mrow><mml:mi>w</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>y</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>≤</mml:mo><mml:mn>?1</mml:mn></mml:mrow></mml:math>0≤w(y)≤?1 for all <mml:math><mml:mi>y</mml:mi></mml:math>y and <mml:math><mml:mrow><mml:mrow><mml:mi>w</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>y</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mn>?1</mml:mn></mml:mrow></mml:math>w(y)=?1 for <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>y</mml:mi></mml:mpadded><mml:mo>≥</mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>y≥?0. The latter property ensures that the integrable component in the fast subspace disappears at <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mpadded><mml:mo>≥</mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>yb≥?0, and therefore, no initial condition can correspond to a particle trapped in the cap forever. In addition, we require the function <mml:math><mml:mi>w</mml:mi></mml:math>w to have a unique minimum at a point <mml:math><mml:msub><mml:mi>y</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:math>y0, such that <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>y</mml:mi><mml:mi>f</mml:mi></mml:msub></mml:mpadded><mml:mo>≤</mml:mo><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>y</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:mpadded><mml:mo><</mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>yf≤y0<?0 (where <mml:math><mml:msub><mml:mi>y</mml:mi><mml:mi>f</mml:mi></mml:msub></mml:math>yf is the pressure equilibrium for the energy level that we are studying). One can check that this inequality ensures that, on every level line of <mml:math><mml:mi>J</mml:mi></mml:math>J, there exist capture moments, such that necessarily <mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>c</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>≠</mml:mo><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>r</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>yb(tc)≠yb(tr), the crucial condition for the destruction of the adiabatic invariant. Finally, taking <mml:math><mml:mrow><mml:mrow><mml:mi>w</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mn>0</mml:mn></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>≥</mml:mo><mml:mn>?0.7</mml:mn></mml:mrow></mml:math>w(y0)≥?0.7 eliminates vertical parabolic orbits in the chaotic component. We expect that the conclusions presented below will apply to other functions fulfilling these requirements.

To test that the proposed mechanism correctly describes the relaxation to the energy equipartition, we conducted a series of numerical experiments with both the springy mushroom and the corresponding stochastic model. Fig. 3B shows that the kinetic energies of the springy mushroom equilibrate. Defining <mml:math><mml:mrow><mml:mrow><mml:mi mathvariant="normal">Δ</mml:mi><mml:mi>K</mml:mi><mml:mi>E</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>t</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:mrow><mml:mi>M</mml:mi><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow><mml:msubsup><mml:mi>v</mml:mi><mml:mi>b</mml:mi><mml:mn>2</mml:mn></mml:msubsup></mml:mrow><mml:mo>?</mml:mo></mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:mrow><mml:mn>1</mml:mn><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow><mml:mrow><mml:mo>?</mml:mo><mml:msub><mml:mi>E</mml:mi><mml:mi>p</mml:mi></mml:msub><mml:mo>?</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:mrow></mml:math>ΔKE(t)=?M/2vb2??1/2?Ep?, we find that, in all our simulations, <mml:math><mml:mrow><mml:mrow><mml:mi mathvariant="normal">Δ</mml:mi><mml:mi>K</mml:mi><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>E</mml:mi><mml:mrow><mml:mi>t</mml:mi><mml:mo>→</mml:mo><mml:mrow><mml:mo>+</mml:mo><mml:mi mathvariant="normal">∞</mml:mi></mml:mrow></mml:mrow></mml:msub></mml:mpadded></mml:mrow><mml:mo>→</mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>ΔKEt→+∞→?0. Here, the ensembles with fixed bar energy are chosen with a random initial phase of the bar motion, a random position of the particle in the stem (we do not start in the cap), and random direction of the particles velocity. Calculating the equilibration rates using the procedure outlined for the SBR model, we again find, as with the SBR model, the sensitivity to the initial energy of the bar and the choice of fitting interval. Nonetheless, for every fixed setting, we find that the rates are well-defined and depend linearly on the square root of the mass ratio, <mml:math><mml:msqrt><mml:mi>m</mml:mi></mml:msqrt></mml:math>m. Extrapolating these rates to <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>m</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>m=?0, we obtain strictly positive limits <mml:math><mml:mrow><mml:mn>1.75</mml:mn><mml:mo>×</mml:mo><mml:msup><mml:mn>10</mml:mn><mml:mrow><mml:mo>?</mml:mo><mml:mn>4</mml:mn></mml:mrow></mml:msup></mml:mrow></mml:math>1.75×10?4 for the ensemble with <mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>E</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mn>0</mml:mn><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mn>?0.9</mml:mn></mml:mrow></mml:math>Eb(0)=?0.9 and <mml:math><mml:mrow><mml:mn>3.08</mml:mn><mml:mo>×</mml:mo><mml:msup><mml:mn>10</mml:mn><mml:mrow><mml:mo>?</mml:mo><mml:mn>4</mml:mn></mml:mrow></mml:msup></mml:mrow></mml:math>3.08×10?4 for the ensemble with <mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>E</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mn>0</mml:mn><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mn>?0.1</mml:mn></mml:mrow></mml:math>Eb(0)=?0.1.

Simulating the stochastic model for the corresponding ensembles of initial conditions (fixed bar energy, random bar phase, motion starting in the chaotic phase) (more details are in Numerical Algorithms), we find that the random model equilibrates in a similar fashion, with positive equilibration rates close to those obtained by the springy billiard simulations when extrapolated to <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>m</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>m=?0 (Fig. 4B). Indeed, the rates for 10 runs of <mml:math><mml:msup><mml:mn>10</mml:mn><mml:mn>4</mml:mn></mml:msup></mml:math>104 particles each are <mml:math><mml:mrow><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mpadded width="+1.7pt"><mml:mn>1.42</mml:mn></mml:mpadded><mml:mo>±</mml:mo><mml:mn>?0.03</mml:mn></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mo>×</mml:mo><mml:msup><mml:mn>?10</mml:mn><mml:mrow><mml:mo>?</mml:mo><mml:mn>4</mml:mn></mml:mrow></mml:msup></mml:mrow></mml:math>(1.42±?0.03)×?10?4 and <mml:math><mml:mrow><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mpadded width="+1.7pt"><mml:mn>2.71</mml:mn></mml:mpadded><mml:mo>±</mml:mo><mml:mn>?0.11</mml:mn></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mo>×</mml:mo><mml:msup><mml:mn>?10</mml:mn><mml:mrow><mml:mo>?</mml:mo><mml:mn>4</mml:mn></mml:mrow></mml:msup></mml:mrow></mml:math>(2.71±?0.11)×?10?4 for <mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>E</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mn>0</mml:mn><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mn>?0.9</mml:mn></mml:mrow></mml:math>Eb(0)=?0.9 and <mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>E</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mn>0</mml:mn><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mn>?0.1</mml:mn></mml:mrow></mml:math>Eb(0)=?0.1, respectively.

We conclude that the stochastic models, with no adaptable parameters, provide a reasonable approximation of the billiard dynamics observed in our numerical experiments. These results support our claim that slow–fast VFS systems are well-approximated by processes with random switching between several different equations for the slow variables.

Springy Stadium

As a control, we examine the equilibration process in a springy billiard, where the frozen billiard is chaotic for all positions of the bar (an EFS system). Specifically, we consider a springy half-stadium (Fig. 1B), which can be obtained from the springy mushroom by removing the collar (i.e., by using <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>w</mml:mi></mml:mpadded><mml:mo>≡</mml:mo><mml:mn>?1</mml:mn></mml:mrow></mml:math>w≡?1 instead of Eq. 11) and keeping the values of all other parameters. Ergodicity of the frozen billiard is proven in ref. 26.

We show that, in such springy billiards, the presence of the adiabatic invariant slows down the relaxation to energy equipartition and that the equilibration rates approach zero as the mass ratio is decreased.

Indeed, direct simulations of the springy stadium show that the adiabatic invariant <mml:math><mml:mi>J</mml:mi></mml:math>J defined by Eq. 3 remains approximately constant and that the bar energy behaves approximately periodically for long time intervals. However, these simulations also show that, for finite <mml:math><mml:mi>m</mml:mi></mml:math>m, small fluctuations ultimately destroy the adiabatic conservation law. In particular, defining <mml:math><mml:mrow><mml:mi mathvariant="normal">Δ</mml:mi><mml:mi>K</mml:mi><mml:mi>E</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>t</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>ΔKE(t) as for the springy mushroom and repeating exactly the same direct numerical simulations for the springy stadium case, we observe again that, for all initial ensembles, <mml:math><mml:mrow><mml:mrow><mml:mi mathvariant="normal">Δ</mml:mi><mml:mi>K</mml:mi><mml:mi>E</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>t</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>→</mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>ΔKE(t)→?0 as <mml:math><mml:mi>t</mml:mi></mml:math>t grows (Fig. 3B). For small values of <mml:math><mml:mi>m</mml:mi></mml:math>m, we observe smaller equilibration rates (but not much smaller) than in the springy mushrooms with comparable parameters (Fig. 4B). Using the best fit line to extrapolate the computed rates to <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>m</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>m=?0, we see that the limit rate vanishes within the accuracy of our computations. Repeating the computations with increased precision, larger ensemble size, and smaller <mml:math><mml:mi>m</mml:mi></mml:math>m, we obtain even smaller limit values for the equilibration rates. We conclude that the limit equilibration rates of the springy stadium are at least two orders of magnitude smaller than the corresponding limit rates for the springy mushroom.

This is our main finding, as it reveals the profound difference between EFS and VFS systems. In agreement with adiabatic theory, the limit equilibration rates vanish in EFS systems. The positive limit rates achieved for VFS systems support our claim of efficient mixing induced by the hopping between ergodic components of the fast subsystem.

Discussion

Springy billiards show an important principle in slow–fast Hamiltonian systems: ergodicity of the fast subsystem impedes ergodization and equilibration of the whole slow–fast system, whereas its violation can lead to equilibration and equipartition of energies. For EFS systems, the averaged lower-dimensional system is known to have an adiabatic invariant, and it is thus not ergodic and does not equilibrate. In contrast, for VFS systems, where the ergodicity of the fast subsystem is broken for a range of values of the slow variables, we showed that the adiabatic behavior is well-approximated by a random process. This process is defined by random switching among several systems that govern the evolution of the slow variables. Importantly, the magnitude and probabilities of the jumps in the stochastic model remain strictly positive in the adiabatic limit. Under mild conditions, such random processes can lead to equilibration of energies with positive rates. We constructed both EFS- and VFS-type springy billiard systems and derived the random models for the VFS springy billiards. We computed the equilibration rates for both the springy billiards and the random models. Extrapolating the rates found at decreasing <mml:math><mml:mi>m</mml:mi></mml:math>m values, we showed that the equilibration rate for the considered EFS system (springy stadium) vanishes in the limit <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>m</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>m=?0, whereas for the VFS systems (SBR and springy mushroom), it remains positive and matches the equilibration rate of the corresponding stochastic model (Fig. 4).

Our equilibration mechanism is also expected to appear in smooth slow–fast systems. The chaotization mechanism that we propose for VFS systems is reminiscent of the phenomena of adiabatic chaos—chaotization of smooth slow—fast systems in which the fast dynamics is integrable, yet the structure of the fast-phase space changes as the slow variables are changed (56?58). We suggest that this mechanism is universal and not restricted to fast subsystems that are integrable. In fact, if some of the ergodic components of the fast system are chaotic and additionally, if there exists a range of slow variables for which the chaotic ergodic component occupies a large portion of the fast subsystem phase space, the equilibration process may be particularly fast.

The proposed framework naturally opens several avenues of research. (i) We saw that the EFS system can equilibrate, despite the presence of an adiabatic invariant. Can a comprehensive theory be built for the behavior beyond the adiabatic timescale? (ii) In the stochastic model for a VFS system, we should typically have convergence to a stationary measure. What if the irreducibility or aperiodicity conditions are violated? For example, the fast subsystem can have a stability island, which persists for all values of slow variables. This could create an invariant domain in the phase space, which has to be excluded from the microcanonical ensemble. (iii) In both EFS and VFS systems, we observe ensemble-dependent rates of convergence to equilibrium. How typical is this phenomenon, and how does it affect the equilibration? (iv) A slow–fast system ceases to be slow–fast near stationary points of the fast subsystem (the fast system freezes for some time). How can one incorporate the freezing incidents into the theory?

Finally, we propose a broader viewpoint on this work. Here, our slow and fast systems were just two mechanical components, and the timescale separation stemmed from the mass ratio. In the broader statistical mechanics context, the fast system governs the motion of many particles (and is, thus, high-dimensional), and the slow macroscopic variables are defined as certain averages over the fast microscopic system. When the structure of the fast system changes, for example, from a gas to a liquid state, a phase transition occurs. Usually, for macroscopic values near the phase transition, the microscopic phase space structure is complex, with long-lasting structures in which the two states coexist. We may say that, in this range of macroscopic variables, the fast subsystem has chaotic (say, gas phase) and ordered (say, liquid state) components. Therefore, we conjecture that phase transitions play a central role in the equilibration process between microscopic and macroscopic variables (e.g., consider the analogue of the notorious piston problem in a multiphase gas).

Numerical Algorithms

Numerical simulations of smooth slow–fast systems are challenging; by considering springy billiards, we avoid the need for numerical integration of stiff equations. Nevertheless, the small mass ratio, <mml:math><mml:mrow><mml:mrow><mml:mi>m</mml:mi><mml:mo>/</mml:mo><mml:mpadded width="+1.7pt"><mml:mi>M</mml:mi></mml:mpadded></mml:mrow><mml:mo>∈</mml:mo><mml:mrow><mml:mo stretchy="false">[</mml:mo><mml:msup><mml:mn>10</mml:mn><mml:mrow><mml:mo>?</mml:mo><mml:mn>7</mml:mn></mml:mrow></mml:msup><mml:mo>,</mml:mo><mml:msup><mml:mn>10</mml:mn><mml:mrow><mml:mo>?</mml:mo><mml:mn>4</mml:mn></mml:mrow></mml:msup><mml:mo stretchy="false">]</mml:mo></mml:mrow></mml:mrow></mml:math>m/M∈[10?7,10?4], leads to a large number of billiard reflections per unit time, and therefore, we have to deal with long-time simulations of chaotic dynamics. Here, we describe the algorithms used for simulations of the springy billiard dynamics and the associated stochastic models. We also include a list of various consistency checks used to verify the accuracy of the computations.

Direct Simulations of Springy Billiards.

We use an explicit formula for the harmonic motion of the bar and the inertial motion of the particle between collisions. The time and place of the collisions are evaluated with the maximal precision allowed by the computer (we use at least double precision). Most of the collisions are with the static walls, and these are computed explicitly. To compute a collision of the particle with the moving bar, we find explicit expressions for the monotonicity intervals of the distance between the particle and the bar as a function of time. Then, we pick out the first monotonicity interval, such that the function changes sign at its ends. This interval contains a single zero of the distance function, which is found using a combination of the bisection and Newton method.

In every series of numerical experiments, we investigate energy equilibration on a fixed energy surface (typically <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>E</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mn>?1</mml:mn></mml:mrow></mml:math>E=?1) in a time interval <mml:math><mml:mrow><mml:mo stretchy="false">[</mml:mo><mml:mn>0</mml:mn><mml:mo>,</mml:mo><mml:msub><mml:mi>T</mml:mi><mml:mn>1</mml:mn></mml:msub><mml:mo stretchy="false">]</mml:mo></mml:mrow></mml:math>[0,T1]. We construct an ensemble of initial conditions that share the same initial value of the bar energy <mml:math><mml:msub><mml:mi>E</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:math>Eb (the details of the ensembles are in the relevant parts of the paper). For each initial condition, we compute the trajectory of the springy billiard for <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mn>0</mml:mn></mml:mpadded><mml:mo>≤</mml:mo><mml:mpadded width="+1.7pt"><mml:mi>t</mml:mi></mml:mpadded><mml:mo>≤</mml:mo><mml:msub><mml:mi>T</mml:mi><mml:mn>1</mml:mn></mml:msub></mml:mrow></mml:math>0≤t≤T1, evaluating the energy components (kinetic and potential energies) at equally spaced moments in time (typically <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>t</mml:mi><mml:mi>k</mml:mi></mml:msub></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:mi>k</mml:mi><mml:mi>h</mml:mi></mml:mrow></mml:mrow></mml:math>tk=kh with <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>h</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:msup><mml:mn>?10</mml:mn><mml:mrow><mml:mo>?</mml:mo><mml:mn>3</mml:mn></mml:mrow></mml:msup></mml:mrow></mml:math>h=?10?3). Then, the ensemble-averaged values of the energy components are computed for every <mml:math><mml:msub><mml:mi>t</mml:mi><mml:mi>k</mml:mi></mml:msub></mml:math>tk and stored to be used for the computation of the equilibration rate. The rate is defined as the best fitted slope to <mml:math><mml:mrow><mml:mi>log</mml:mi><mml:mrow><mml:mo>|</mml:mo><mml:mrow><mml:mi mathvariant="normal">Δ</mml:mi><mml:mi>K</mml:mi><mml:mi>E</mml:mi></mml:mrow><mml:mo>|</mml:mo></mml:mrow></mml:mrow></mml:math>log|ΔKE| vs. <mml:math><mml:mi>t</mml:mi></mml:math>t over a time interval <mml:math><mml:mrow><mml:mo stretchy="false">[</mml:mo><mml:mn>0</mml:mn><mml:mo>,</mml:mo><mml:msub><mml:mi>T</mml:mi><mml:mn>0</mml:mn></mml:msub><mml:mo stretchy="false">]</mml:mo></mml:mrow></mml:math>[0,T0], so that <mml:math><mml:mrow><mml:mrow><mml:mi>log</mml:mi><mml:mrow><mml:mo stretchy="false">|</mml:mo><mml:mrow><mml:mrow><mml:mrow><mml:mi mathvariant="normal">Δ</mml:mi><mml:mi>K</mml:mi><mml:mi>E</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mn>0</mml:mn><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>/</mml:mo><mml:mi mathvariant="normal">Δ</mml:mi></mml:mrow><mml:mi>K</mml:mi><mml:mi>E</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>T</mml:mi><mml:mn>0</mml:mn></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo stretchy="false">|</mml:mo></mml:mrow></mml:mrow><mml:mo>≈</mml:mo><mml:mn>?1</mml:mn></mml:mrow></mml:math>log|ΔKE(0)/ΔKE(T0)|≈?1.

Theoretically, for an infinitely large ensemble, <mml:math><mml:msub><mml:mi>T</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:math>T0 can be precisely defined by setting it to be the first time where <mml:math><mml:mrow><mml:mrow><mml:mi>log</mml:mi><mml:mrow><mml:mo stretchy="false">|</mml:mo><mml:mrow><mml:mrow><mml:mrow><mml:mi mathvariant="normal">Δ</mml:mi><mml:mi>K</mml:mi><mml:mi>E</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mn>0</mml:mn><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>/</mml:mo><mml:mi mathvariant="normal">Δ</mml:mi></mml:mrow><mml:mi>K</mml:mi><mml:mi>E</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>T</mml:mi><mml:mn>0</mml:mn></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo stretchy="false">|</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mn>?1</mml:mn></mml:mrow></mml:math>log|ΔKE(0)/ΔKE(T0)|=?1. For finite ensembles, this definition leads to noisy <mml:math><mml:msub><mml:mi>T</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:math>T0 for different realizations with a small systematic underestimate of <mml:math><mml:msub><mml:mi>T</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:math>T0. In the case of the SBR, for all <mml:math><mml:mi>m</mml:mi></mml:math>m, we find that <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>T</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:mpadded><mml:mo>≈</mml:mo><mml:mn>?10.5</mml:mn></mml:mrow></mml:math>T0≈?10.5, and its variations do not affect the consistency of the results, as different realizations produce similar rates (Fig. 4A). In the case of the springy mushroom and stadium, a higher accuracy is needed to study the behavior of the <mml:math><mml:mi>m</mml:mi></mml:math>m-dependent rates in the limit <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>m</mml:mi></mml:mpadded><mml:mo>→</mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>m→?0. Therefore, to help reduce uncertainty in the definition of the rates, we choose an ensemble-independent <mml:math><mml:msub><mml:mi>T</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:math>T0. To set <mml:math><mml:msub><mml:mi>T</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:math>T0 as a function of <mml:math><mml:mi>m</mml:mi></mml:math>m only, we first find the <mml:math><mml:mi>m</mml:mi></mml:math>m-dependent rates using ensemble-dependent fitting time intervals. Then, we define a linear function <mml:math><mml:mrow><mml:mi>R</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msqrt><mml:mi>m</mml:mi></mml:msqrt><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>R(m) by the best linear fit for these rates as a function of <mml:math><mml:mi>m</mml:mi></mml:math>m, and let <mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>T</mml:mi><mml:mn>0</mml:mn></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>m</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:mn>?1</mml:mn><mml:mo>/</mml:mo><mml:mi>R</mml:mi></mml:mrow><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msqrt><mml:mi>m</mml:mi></mml:msqrt><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>T0(m)=?1/R(m). Finally, we recompute the rates using <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>T</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:msub><mml:mi>T</mml:mi><mml:mn>0</mml:mn></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>m</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>T0=T0(m). This procedure leads to only a small change in the obtained rate values (the difference is about 2% of their first estimate). Nevertheless, the extrapolation to <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>m</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>m=?0 is improved; as an example, for the stadium, the extrapolated rates get closer to zero, which is the theoretical rate.

The computation of a trajectory segment for <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>t</mml:mi></mml:mpadded><mml:mo>∈</mml:mo><mml:mrow><mml:mo stretchy="false">[</mml:mo><mml:mn>0</mml:mn><mml:mo>,</mml:mo><mml:msub><mml:mi>T</mml:mi><mml:mn>1</mml:mn></mml:msub><mml:mo stretchy="false">]</mml:mo></mml:mrow></mml:mrow></mml:math>t∈[0,T1] (<mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:msub><mml:mi>T</mml:mi><mml:mn>1</mml:mn></mml:msub></mml:mpadded><mml:mo>=</mml:mo><mml:mn>21</mml:mn></mml:mrow></mml:math>T1=21 for the SBR; <mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi>T</mml:mi><mml:mn>1</mml:mn></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>m</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>∈</mml:mo><mml:mrow><mml:mo stretchy="false">[</mml:mo><mml:mn>200,3,000</mml:mn><mml:mo stretchy="false">]</mml:mo></mml:mrow></mml:mrow></mml:math>T1(m)∈[200,3,000] for the mushroom and the stadium) includes <mml:math><mml:msup><mml:mn>10</mml:mn><mml:mn>3</mml:mn></mml:msup></mml:math>103<mml:math><mml:msup><mml:mn>10</mml:mn><mml:mn>6</mml:mn></mml:msup></mml:math>106 collisions. To maintain a high level of accuracy, the following checks of the code are carried out. (i) The total energy of the system for the full duration of the simulations is conserved up to at least nine significant digits. (ii) When the code is run for a static billiard (i.e., the position of the bar is fixed), the particle energy is preserved with at least nine significant digits. (iii) For selected parameter values, the code is run several times to estimate the size of the fluctuations caused by the finite size of the ensemble (e.g., the reported SDs for the equilibration rates). (iv) The experiments are repeated with different ensemble sizes.

Stochastic Model Simulations.

In the stochastic models, the motion of the bar is described by several explicitly defined Hamiltonian systems with one degree of freedom (Eqs. 7 and 9) and a probabilistic law, which determines transitions between these systems. The stochastic models correspond to the limit of <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>m</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mn>?0</mml:mn></mml:mrow></mml:math>m=?0 and are independent of the mass ratio.

Stochastic simulations for the SBR.

In this case, the bar motion is described by Eq. 7, which is integrated with the help of the standard Runge–Kutta method of order 4 and a time step size of <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>h</mml:mi></mml:mpadded><mml:mo>≈</mml:mo><mml:msup><mml:mn>?10</mml:mn><mml:mrow><mml:mo>?</mml:mo><mml:mn>4</mml:mn></mml:mrow></mml:msup></mml:mrow></mml:math>h≈?10?4. When <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>t</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:mn>?0</mml:mn><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mo>mod</mml:mo><mml:mi>T</mml:mi></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>t=?0(modT), a random choice of <mml:math><mml:msub><mml:mi>F</mml:mi><mml:mi>j</mml:mi></mml:msub></mml:math>Fj is made with the prescribed probabilities, and at <mml:math><mml:mrow><mml:mpadded width="+1.7pt"><mml:mi>t</mml:mi></mml:mpadded><mml:mo>=</mml:mo><mml:mrow><mml:mi>τ</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mo>mod</mml:mo><mml:mi>T</mml:mi></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>t=τ(modT), the right-hand side of Eq. 5 is reset to zero. The initial conditions for the bar motion are set to reflect the properties of the ensembles used in the simulations of the SBR: the initial energy of the bar takes the same value, while its phase on the <mml:math><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo>,</mml:mo><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:math>(yb,y˙b) plane is random. To mimic the uniform distribution of the particle position in the box, a random phase is chosen from which one starts the integration of Eq. 7. If the phase is in <mml:math><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mn>0</mml:mn><mml:mo>,</mml:mo><mml:mrow><mml:mrow><mml:mi>τ</mml:mi><mml:msub><mml:mi>T</mml:mi><mml:mi>p</mml:mi></mml:msub></mml:mrow><mml:mo>/</mml:mo><mml:mi>T</mml:mi></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:math>(0,τTp/T), a random choice of the up or down system is made with probabilities <mml:math><mml:msub><mml:mi>β</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:math>β0 and <mml:math><mml:mrow><mml:mn>1</mml:mn><mml:mo>?</mml:mo><mml:msub><mml:mi>β</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:mrow></mml:math>1?β0, respectively.

Stochastic simulations for the springy mushroom.

In contrast to the SBR model, the switching time between the two forces in Eq. 9 is stochastic and governed by the probabilistic law of Eq. 10. We wrote two different numerical codes to carry out these simulations. The rates obtained by these two methods coincide within the accuracy determined by the size of the ensembles.

Algorithm 1—Direct integration and switching: We integrate Eq. 9 in the interval <mml:math><mml:mrow><mml:mo stretchy="false">[</mml:mo><mml:mi>t</mml:mi><mml:mo>,</mml:mo><mml:mrow><mml:mi>t</mml:mi><mml:mo>+</mml:mo><mml:mrow><mml:mi>d</mml:mi><mml:mi>t</mml:mi></mml:mrow></mml:mrow><mml:mo stretchy="false">]</mml:mo></mml:mrow></mml:math>[t,t+dt] and find <mml:math><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mi>t</mml:mi><mml:mo>+</mml:mo><mml:mrow><mml:mi>d</mml:mi><mml:mi>t</mml:mi></mml:mrow></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>yb(t+dt) and <mml:math><mml:mrow><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mi>t</mml:mi><mml:mo>+</mml:mo><mml:mrow><mml:mi>d</mml:mi><mml:mi>t</mml:mi></mml:mrow></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>y˙b(t+dt). Then, we switch to the harmonic motion of the bar with probability <mml:math><mml:mrow><mml:mi>P</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>t</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mi>d</mml:mi><mml:mi>t</mml:mi></mml:mrow></mml:math>P(t)dt given by Eq. 10, where <mml:math><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:math>yb is replaced by <mml:math><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mi>t</mml:mi><mml:mo>+</mml:mo><mml:mrow><mml:mi>d</mml:mi><mml:mi>t</mml:mi></mml:mrow></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>yb(t+dt) and <mml:math><mml:mrow><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi></mml:msub><mml:mi>d</mml:mi><mml:mi>t</mml:mi></mml:mrow></mml:math>y˙bdt is replaced by <mml:math><mml:mrow><mml:mrow><mml:mi>d</mml:mi><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mi>t</mml:mi><mml:mo>+</mml:mo><mml:mrow><mml:mi>d</mml:mi><mml:mi>t</mml:mi></mml:mrow></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>t</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:mrow></mml:math>dyb=yb(t+dt)?yb(t). After the particle is captured, the position and velocity of the bar at the time of release are calculated explicitly using the equations of harmonic motion.

The initial ensemble for the bar position is the same as in the direct simulations, and since in the direct simulations, the particle always starts from the stem, the chaotic simulation always starts with the nontrapped particle. To test the accuracy of our code, we checked that the achieved rates converge as the integration step is decreased (we typically used <mml:math><mml:mrow><mml:mrow><mml:mi>d</mml:mi><mml:mpadded width="+1.7pt"><mml:mi>t</mml:mi></mml:mpadded></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:msub><mml:mi>T</mml:mi><mml:mn>1</mml:mn></mml:msub><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow><mml:mo>,</mml:mo><mml:mn>000</mml:mn></mml:mrow></mml:mrow></mml:math>dt=T1/2,000 in the stochastic simulations). We also checked that repeated runs produce statistically similar results.

Algorithm 2—Switching without integration: We begin by preparing tables, which provide, for each value of the adiabatic invariant <mml:math><mml:mi>J</mml:mi></mml:math>J, the period of the bar’s oscillation in the effective potential of Eq. 9, the probability of the particle being captured in the cap during one period, and when the particle is captured, the cumulative distribution function (CDF) for the bar position <mml:math><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:math>yb at the moment of capture. After that, for any initial condition <mml:math><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub><mml:mo>,</mml:mo><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>y</mml:mi><mml:mo>˙</mml:mo></mml:mover></mml:mrow><mml:mi>b</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:math>(yb,y˙b) (assuming that the particle is in the chaotic component), we generate the number of completed rounds before capture and the value of <mml:math><mml:msub><mml:mi>y</mml:mi><mml:mi>b</mml:mi></mml:msub></mml:math>yb at the moment of capture as random variables, with the CDFs obtained by interpolation of the precomputed tables. After the capture point is generated, the time of release and the new position and velocity of the bar are computed from the harmonic motion law. Then, the process is repeated until the time exceeds <mml:math><mml:msub><mml:mi>T</mml:mi><mml:mn>1</mml:mn></mml:msub></mml:math>T1.

The outcome of this stochastic algorithm is a list of capture and release points (time, bar position, and bar velocity). From this list, we find the bar energy at fixed time intervals: it is constant in the capture phase and found by using the adiabatic invariance in the chaotic phase. This algorithm enables fast simulation without integration. Agreement with the results of the simulations run using Algorithm 1 provides additional verification of those results.

Acknowledgments

Part of this work was done while K.S. was at the Indian Institute of Technology Delhi. K.S. was supported by Science and Engineering Research Board, Government of India File SR/FTP/PS-108/2012. K.S., D.T., and V.G. acknowledge the financial support and hospitality of the Weizmann Institute of Science, where part of this work was done. D.T. was supported by Russian Science Foundation Grant 14-41-00044 for this research and Engineering and Physical Sciences Research Council (EPSRC) Grant EP/P026001/1 and thanks the Royal Society. The research of V.G. was supported by EPSRC Grant EP/J003948/1. V.R.-K. acknowledges support from Israel Science Foundation Grant 1208/16.

Footnotes

  • ?1To whom correspondence should be addressed. Email: vered.rom-kedar{at}weizmann.ac.il.
  • Author contributions: K.S., D.T., V.G., and V.R.-K. designed research; K.S., D.T., V.G., and V.R.-K. performed research; K.S., D.T., V.G., and V.R.-K. contributed new reagents/analytic tools; K.S., D.T., V.G., and V.R.-K. analyzed data; K.S., D.T., V.G., and V.R.-K. conceived work ideas by discussions; K.S., V.G., and V.R.-K. performed simulations; and K.S., D.T., V.G., and V.R.-K. wrote the paper.

  • The authors declare no conflict of interest.

  • This article is a PNAS Direct Submission. J.-P.E. is a guest editor invited by the Editorial Board.

References

  1. ?
    .
  2. ?
    .
  3. ?
    .
  4. ?
    .
  5. ?
    .
  6. ?
    .
  7. ?
    .
  8. ?
    .
  9. ?
    .
  10. ?
    .
  11. ?
    .
  12. ?
    .
  13. ?
    .
  14. ?
    .
  15. ?
    .
  16. ?
    .
  17. ?
    .
  18. ?
    .
  19. ?
    .
  20. ?
    .
  21. ?
    .
  22. ?
    .
  23. ?
    .
  24. ?
    .
  25. ?
    .
  26. ?
    .
  27. ?
    .
  28. ?
    .
  29. ?
    .
  30. ?
    .
  31. ?
    .
  32. ?
    .
  33. ?
    .
  34. ?
    .
  35. ?
    .
  36. ?
    .
  37. ?
    .
  38. ?
    .
  39. ?
    .
  40. ?
    .
  41. ?
    .
  42. ?
    .
  43. ?
    .
  44. ?
    .
  45. ?
    .
  46. ?
    .
  47. ?
    .
  48. ?
    .
  49. ?
    .
  50. ?
    .
  51. ?
    .
  52. ?
    .
  53. ?
    .
  54. ?
    .
  55. ?
    .
  56. ?
    .
  57. ?
    .
  58. ?
    .

Online Impact