# Equilibration of energy in slow–fast systems

^{a}Department of Electrical Engineering and Computer Science, Indian Institute of Science Education and Research, Bhopal 462066, India;^{b}Department of Mathematics, Imperial College, London SW7 2AZ, United Kingdom;^{c}Lobachevsky University of Nizhny, Novgorod 603950, Russia;^{d}Mathematics Institute, University of Warwick, Coventry CV4 7AL, United Kingdom;^{e}Department of Computer Science and Applied Mathematics, The Weizmann Institute of Science, Rehovot 76100, Israel

See allHide authors and affiliations

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 M is taken to be much larger than the particle mass 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

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. 1*B*). 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. 1*A*) and springy mushroom (Fig. 1*C*)]. Such systems, with variable partition of the fast subspace (VFS), are hereafter called VFS systems.

All of our numerical simulations for nonzero values of

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

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

The “bar–particle” system is Hamiltonian and has

If the frozen billiard is ergodic for each value of **2** becomes

A straightforward computation shows that solutions of Eq. **5** indeed follow the level sets of J. By noting that the adiabatic law given by Eq. **3** takes the form

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. 1*B*), 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 L and height *A*). The particle’s horizontal speed **3** with **2** becomes

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

Hence, we suggest that the bar–particle system is well-approximated by the following

The bar dynamics follow the level lines of

In Fig. 2*A*, we plot the projection of direct simulations of the bar–particle system onto the *A* shows that the bar motion indeed follows the level sets of the corresponding J values and that the trajectory switches between level sets (Fig. 2*A*, *Upper Left Inset*); while the bar energy changes significantly when the particle is above/below the bar (Fig. 2*A*, *Lower Right Inset*), the corresponding J values remain approximately constant for each such interval (Fig. 2*A*, *Lower Left Inset*).

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,

In all of our simulations, we find that *A*). However, for each fixed 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 *Numerical Algorithms*.

We examine how the rates depend on the mass ratio. For sufficiently small m, the rates do not display a significant dependence on m (Fig. 3*A*) 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. 4*A*, we obtain very similar rates; 10 runs of the

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. 1*C*). 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 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 w, and the island vanishes when

In the springy mushroom, the bar located at the bottom of the stem can move vertically. It has mass *C*). 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: **11**). Then, the island size depends on the bar position

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:**4** (with

Expressions for the volumes in Eq. **8** can be found explicitly. Indeed, at the energy level

When the particle is in the chaotic zone, it acts on the bar with average force **8**, we conclude that the averaged equations for the bar motion are of the form

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. 2*B* shows that trajectories of the springy mushroom follow a level set of J (Fig. 2*B*, red circles) and then, by following a level set of the free motion (Fig. 2*B*, green circles), switch to follow another level set of J. Fig. 2*B*, *Lower Left Inset* shows that J is approximately constant during the time intervals in which the particle remains in the chaotic zone, and Fig. 2*B*, *Lower Right Inset* shows the corresponding oscillations of

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

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

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 **10**. If captured, the particle is released at time

Now we are ready to describe the mechanism for the ergodization in the springy mushroom. Choosing a nonmonotonic function

For example, in our numerical experiments, we use the nonmonotonic function

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. 3*B* shows that the kinetic energies of the springy mushroom equilibrate. Defining

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 *B*). Indeed, the rates for 10 runs of

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. 1*B*), which can be obtained from the springy mushroom by removing the collar (i.e., by using **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 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 m, small fluctuations ultimately destroy the adiabatic conservation law. In particular, defining *B*). For small values of m, we observe smaller equilibration rates (but not much smaller) than in the springy mushrooms with comparable parameters (Fig. 4*B*). Using the best fit line to extrapolate the computed rates to

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 m values, we showed that the equilibration rate for the considered EFS system (springy stadium) vanishes in the limit

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,

### 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

Theoretically, for an infinitely large ensemble, *A*). In the case of the springy mushroom and stadium, a higher accuracy is needed to study the behavior of the m-dependent rates in the limit

The computation of a trajectory segment for *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

#### 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 **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 **7**. If the phase is in

#### 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

**10**, where

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

** Algorithm 2—Switching without integration:** We begin by preparing tables, which provide, for each value of the adiabatic invariant 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

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

- ↵
^{1}To 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.

- Copyright © 2017 the Author(s). Published by PNAS.

This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

## References

- ↵
- Landau LD,
- Lifshitz EM

- ↵
- Huang K

- ↵
- Ruelle D

- ↵
- Dorfman JR

- ↵
- Hilbert S,
- Hänggi P,
- Dunkel J

- ↵
- ↵
- Sinai YG,
- Chernov N

- ↵
- Simanyi N,
- Szasz D

- ↵
- ↵
- Rapoport A,
- Rom-Kedar V,
- Turaev D

- ↵
- ↵
- Zaslavsky GM

- ↵
- Lichtenberg AJ,
- Lieberman M

- ↵
- Anosov DV

- ↵
- Kasuga T

- ↵
- Lochak P,
- Meunier C

- ↵
- Lieb EH

- ↵
- Szasz D

- Lebowitz JL,
- Piasecki J,
- Sinai YG

- ↵
- Neishtadt AI,
- Sinai YG

- ↵
- Wright P

- ↵
- Dolgopyat D,
- Liverani C

- ↵
- Grigo A,
- Khanin K,
- Szasz D

- ↵
- Balint P,
- Gilbert T,
- Nándori P,
- Domokos S,
- Toth IP

- ↵
- Chernov N,
- Markarian R

- ↵
- Hannay JH,
- McCraw RJ

- ↵
- Bunimovich LA

- ↵
- ↵
- Schmidt G,
- Chen Q

- ↵
- Lamba H

- ↵
- ↵
- Luo ACJ,
- Guo Y

- ↵
- Lichtenberg AJ,
- Lieberman MA,
- Cohen RH

- ↵
- Jarzynski C

- ↵
- Loskutov A,
- Ryabov AB,
- Akinshin LG

- ↵
- de Carvalho RE,
- Souza FC,
- Leonel ED

- ↵
- ↵
- Dolgopyat D

- ↵
- Shah K

- ↵
- ↵
- ↵
- Gelfreich V,
- Rom-Kedar V,
- Turaev D

- ↵
- Batistic B

- ↵
- Turaev D

- ↵
- ↵
- Shah K,
- Turaev D,
- Rom-Kedar V

- ↵
- Pereira T,
- Turaev D

- ↵
- Khanin K,
- Yarmola T

- ↵
- Mokhtar-Kharroubi M,
- Rudnicki R

- ↵
- Eckmann J,
- Young L

- ↵
- Konishi T,
- Yanagita T

- ↵
- Nitecki Z,
- Robinson C

- Burton R,
- Easton RW

- ↵
- Sivaramakrishnan A

- ↵
- Berger P,
- Carrasco P

- ↵
- Blumenthal A,
- Xue J,
- Young LS

- ↵
- Primack H,
- Smilansky U

- ↵
- ↵
- Vainshtein DL,
- Vasiliev AA,
- Neishtadt AI

- ↵
- Neishtadt AN,
- Chaikovskii DK,
- Chernikov AA

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Applied Mathematics