Emergent order from mixed chaos at low temperature
Pavel Chvykov, Jeremy England

TL;DR
This paper shows how cold environments can lead to orderly behavior in chaotic systems, revealing a new link between thermodynamics and dynamics.
Contribution
The novel connection between thermodynamic and dynamical systems perspectives in the emergence of order at low temperatures is established.
Findings
Numerical evidence supports that mixed chaotic systems exhibit regular behavior when coupled to a cold thermal bath.
Transition temperatures can be predicted using relaxation timescales and a new fluctuation-dissipation relation.
Emergent order is robust and non-trivial, not just reduced motion as in equilibrium.
Abstract
This paper explores a novel connection between a thermodynamic and a dynamical systems perspective on emergent dynamical order. We provide evidence for a conjecture that Hamiltonian systems with mixed chaos spontaneously find regular behavior when minimally coupled to a thermal bath at sufficiently low temperature. Numerical evidence across a diverse set of five dynamical systems supports this conjecture, and allows us to quantify corollaries about the organization timescales and disruption of order at higher temperatures. Balancing the damping-induced phase-space contraction against thermal exploration, we are able to predict the transition temperatures in terms of the relaxation timescales, indicating a novel nonequilibrium fluctuation-dissipation relation, and formally connecting the thermodynamic and dynamical systems views. Our findings suggest that for a wide range of real-world…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsAdvanced Thermodynamics and Statistical Mechanics · Quantum many-body systems · Chaos control and synchronization
Introduction
Emergence of order in space, time, or response properties in dynamical systems generally depends on the selection of a low-entropy or low-phase-space-volume ensemble of microstates, which is stabilized by system dynamics relative to a much wider range of microstates that were available in principle. Thus, understanding general mechanisms by which physical dynamics can lead to the selection of fine-tuned patterns has broad relevance to the study of self-organization and emergent order in nature and applications^1,2^.
Prigogine is known for pointing out that dissipation of energy from an external drive seems to play a key role in various examples of self-organization^3^. In particular, he established the intuition that certain kinds of non-trivial dynamical or spatial order could only emerge once the strength of external driving pushed the system in question beyond the linear-response regime^4^. However, Prigogine’s program ultimately failed to yield a general theory of the so-called dissipative structures he was credited with characterizing^5^.
An interesting hint indicating how Prigogine may have been right comes from the more mathematical work in dynamical systems theory, where there is a common intuition that adding weak damping can lead to emergence of order in mixed-chaotic systems^6–9^. It is, however, unclear how generally or formally this intuition is thought to apply. Additionally, the focus of past work on this phenomenon has been in low-dimensional systems, making the connection to many-body self-organization tenuous. Still, insofar as damping in physical systems is a mechanism of dissipation, there may as a result turn out to be a broad class of dynamical systems that might be said to exhibit Prigogine’s dissipative structure.
The key to putting this connection on firmer physical ground is to realize that damping and temperature are inextricably linked. Subjecting a Hamiltonian system to damping without noise, as was typically done in the aforementioned mixed chaos studies, is equivalent to coupling the system to a zero temperature thermal bath. The Einstein relation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b\,D=k_{B}T$$\end{document} , where D is the diffusion constant, b is damping coefficient, and T is temperature, quantifies this connection, which arises because the bath of particles responsible for the damping drag must also be the source of random thermal fluctuations.
In this sense, the order arising from noiseless damping in mixed-chaotic dynamical systems could potentially be attributed to the impact of being at low temperature. The low-temperature regime is known to bring about spontaneous order in thermal equilibrium, but can it produce a similar effect for driven dynamical systems? In this paper, we study the phenomenon of emergent order brought on by coupling mixed-chaotic dynamical systems to a low-temperature thermal bath. We thus combine the perspectives of dynamical systems theory (via damping) and thermodynamics (via low temperature), and argue for the generality of this mechanism, including in many-body systems.
The clearest thermodynamic argument for emergent order comes from cases involving equilibrium self-assembly, in which a multitude of attractive and repulsive inter-particle forces in a many-body system on the microscale can lead to the formation of useful meso- and macrostructures such as proteins and crystals^10,11^. In typical cases of equilibrium self-assembly, there is a clear explanation for the low-entropy phenomenon that comes from the effect of temperature: when such systems decrease in average energy, they are no longer free to explore a vast diversity of collective states and instead become trapped in rare local energy minima.
Away from equilibrium, for a dynamical system subject to time-varying external forces, we can write a generalization of the Boltzmann distribution for the steady-state of a system with microstates \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_i$$\end{document} at energies \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_i$$\end{document} and temperature \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_B T$$\end{document} :
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} p(x_{i}) \propto \text{ exp }\left[ -\frac{u_i}{k_{B}T}\right] \, \bigg \langle \text{ exp }\left[ -\frac{w_{\rightarrow i}}{k_{B}T}\right] \bigg \rangle ^{-1} \end{aligned}$$\end{document}The path-dependence characteristic of nonequilibrium systems is captured here by the second term, which depends on the work \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w_{\rightarrow i}$$\end{document} done by external driving on the system as it relaxes from an arbitrary initial starting state into the final state i. The averaging \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {\cdot }\rangle$$\end{document} here is done over all stochastic relaxation paths (allowing that the notation for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w_{\rightarrow i}$$\end{document} implies taking the required \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t\rightarrow \infty$$\end{document} limit only after normalizing p). Considering the case of constant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_i$$\end{document} in order to focus on purely dynamical effects, we can see that just as in equilibrium self-assembly, lowering temperature has the potential to concentrate the probability in a rare set of states, resisting the disordering pull of entropy. Here, however, such states are selected for their exceptionally positive-leaning distribution of dissipative work history \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w_{\rightarrow i}$$\end{document} , instead of for their low energy. But what kind of order might these states and associated trajectories reflect? And how is that order stabilized without the attractive interparticle forces that hold things together in the equilibrium case?^12,13^.
The dynamical systems theory perspective on mixed chaotic systems helps to provide an explanation. Mixed chaos is defined as the coexistence of both chaotic and regular (non-chaotic) behaviors in different regions of the dynamical system’s phase space. While this might seem unusual at first, it turns out that this is the “rule rather than the exception” among non-linear Hamiltonian dynamical systems^14^. Intuitively, this happens because it is less likely for a real-world system to be 100% chaotic or regular, than for it to have some of both in its phase space. Formally, this arises as a consequence of KAM (Kolmogorov–Arnold–Moser) theory, which tells us that mixed chaos dynamics will usually emerge when we add small nonlinear interactions to an otherwise linear (or integrable) dynamical system. As the nonlinearity is increased, progressively larger regions of the phase space can become chaotic, while regular behavior can be maintained in some corners of the phase space—leaving islands of regularity in the sea of chaos (Fig.1-iii).
The emergent order in dynamical systems theory is understood in terms of damping, rather than low temperature as in thermodynamics^15–17^ . The key effect of damping on a Hamiltonian system is to break the phase-space volume conservation of Liouville’s theorem. For undamped dynamics \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dot{p} = -\partial H/\partial q$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dot{q}= \partial H/\partial p$$\end{document} , the identity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nabla \cdot [\dot{q},\dot{p}]=\partial ^{2}H/\partial q\partial p-\partial ^{2}H/\partial p\partial q=0$$\end{document} ensures incompressible flow in phase space ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d\rho /dt = 0$$\end{document} ) so that each infinitesimal element dq dp cannot be contracted or expanded over time. Thus, a non-dissipative Hamiltonian system cannot lose memory of initial conditions. Additionally, chaotic and ordered dynamical regions of phase space are forbidden to intermix since the closed regular orbits inside the islands of regularity cannot suddenly start moving chaotically, either forward or backward in time. By contrast, introduction of damping via \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dot{p} = -\partial H/\partial q - bp$$\end{document} , results in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d\rho /dt =-\rho \nabla \cdot [\dot{q},\dot{p}] = b\rho$$\end{document} , which permits phase-space volume contraction into a low-entropy subset of states. Damping is necessary for such concentration of probability, but not sufficient, as chaos in principle could stretch phase-space volumes into space-filling filaments. In mixed chaotic systems, however, initially chaotic trajectories would be expected to explore diffusively until coming arbitrarily close to an ordered trajectory, whereas ordered trajectories stay ordered over time. Due to this asymmetry, the accumulation of probability density in ordered regions becomes favored. Put another way, islands of regularity in Hamiltonian dynamics are regions with vanishing local Lyapunov exponents \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _i=0$$\end{document} , meaning that volumes do not get stretched in any dimension. The addition of weak damping drives exponents below zero and can make these regions attracting^7,18^.
In this paper we combine thermodynamic and dynamical systems perspectives by developing numerical evidence and physical intuition for the following conjecture:
Conjecture 1
Hamiltonian dynamics exhibiting mixed chaos will settle into islands of order upon weakly coupling to a thermal bath at sufficiently low temperature.
In Section 2, we demonstrate this phenomenon in a variety of example systems, including a many-body example suggesting relevance for self-organization. Any formal proofs of this conjecture are left to future work. To better characterize the thermodynamic nature of the phenomenon, we further numerically study the relation between damping strength and relaxation time, discovering a power-law relation (Section 3), as well as the phase diagram that shows how order disappears at higher temperatures (Section 4). Relating these two, we show preliminary evidence for a possible nonequilibrium fluctuation-dissipation relation in these systems.
Steady-state
We begin by showing that our conjecture 1 works consistently across five different dynamical systems in-silico, as shown in the corresponding five rows (a–e) of Fig.1. In each case, we see that regions that used to be islands of regularity under Hamiltonian dynamics (columns (ii) and (iii)), become global attractors as soon as we couple to the thermal bath at low temperature T (column (iv)). Mathematically, by “coupling to a thermal bath” we mean adding damping and noise forces to our system \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta F = -b\,v + \sqrt{2\,b\,T}\;\xi (t)$$\end{document} , where b is the damping coefficient and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi (t)$$\end{document} is white Gaussian noise process (such that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle \xi (t)\rangle =0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle \xi (t),\, \xi (t')\rangle = \delta (t-t')$$\end{document} ). The noise amplitude is determined in terms of T via the Einstein relation. We can thus think of b as the “bath coupling strength,” as it controls both terms, restoring Hamiltonian dynamics in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b \rightarrow 0$$\end{document} limit. This setup clarifies that adding damping without noise – a simplification used in most literature on damped chaos – corresponds to the idiosyncratic \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T=0$$\end{document} thermal bath. Also, throughout this work we assume the bath coupling to be a relatively weak perturbation on top of baseline Hamiltonian dynamics, meaning that it does not significantly change the local phase-space structure as heavy damping would^18^. If damping was large, above a critical value it is known that new strange attractors can appear that have no signature in the undamped Hamiltonian system^6^.Fig. 1. Numerical evidence for conjecture 1 across 5 dynamical systems. Column (ii) shows the values of local Lyapunov exponents for points sampled throughout the system’s phase, quantitatively identifying islands of regularity. These are then also seen in the system trajectories under baseline Hamiltonian dynamics in column (iii), which we then weakly couple to a thermal bath (add damping and noise) to get the steady-state dynamics pictured in column (iv). In (iii) and (iv), we plot stroboscopic system configurations for 256 initial conditions (IC), picked randomly in the shown phase space. Each IC and subsequent trajectory is plotted in its own color (so points of one color correspond to one trajectory). Each IC is evolved under one unique noise realization. The coordinates shown are: (a,b) rotor angle and speed right before kick, (e) same, but only for one of the 16 coupled rotors (the one marked black in the network in (e–i), (c) oscillator position and speed when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega \, t =0$$\end{document} (mod \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\pi$$\end{document} ), (d) plate phase and ball speed immediately after each bounce. (e) is the only many-body example, and is more challenging to visualize: regular behavior is very rare in the 32-dim system phase space, and so none of the randomly sampled IC demonstrate it. The IC shown combine random sampling with sampling from the damped steady-state attractors. Also, since here partial organization is possible, attractors need not be 0-dimensional – so we measure their fractal dimension to quantify degree of organization. In all cases, we see that bath coupling collapses all trajectories to the Lyapunov \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda = 0$$\end{document} islands of regularity, giving emergent order.
To quantitatively characterize the islands of regularity, we measure the local Lyapunov exponents \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _i$$\end{document} throughout each system’s phase space—Fig. 1 column (ii). To do this, for each point we start two trajectories at nearby initial conditions, and run them for 20 steps of our stroboscopic Poincare map (for 50 in (e)), monitoring how their separation grows. This helps to understand the role of local Lyapunov exponent in selecting the dynamically ordered behaviors under thermal bath coupling.
Remarkably, the five examples in Fig. 1 show that our conjecture continues to hold regardless of how weak the bath coupling b is, and regardless of the degree of fine-tuning required (see Fig.2) or of the complexity of the regular motion (e.g., see the attractor shapes in Fig.1(e-iv).
The first system we considered is the kicked rotor (or Chirikov standard map, Fig. 1a), which is the paradigmatic example of mixed chaos, and has been extensively studied in context of dynamical systems theory (as well as quantum chaos and Anderson localization)^19,20^. The standard setup is as a conservative Hamiltonian dynamics, which we then couple to a thermal bath (see Ref.^8^ for a more technical analysis of this system):
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \dot{\theta }&= v \nonumber \\ \dot{v}&= -K\,\sin (\theta )\;\sum _n\delta (t-n) - b\,v + \sqrt{2\,b\,T}\;\xi (t) \end{aligned}$$\end{document}This describes the dynamics \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta (t)$$\end{document} of a particle hinged on a freely-rotating rigid rod, and periodically kicked by a global uniform force field K (where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta$$\end{document} is the Dirac delta and n – any integer). The advantage of this system is that it’s very fast to simulate, as we can integrate the linear dynamics between kicks analytically, effectively turning the system into a discrete stroboscopic map with a single simulation step per kick. The stochastic noise can similarly be added to the linear evolution by scaling its amplitude by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{\Delta t}$$\end{document} – integration time. It has previously been noted that this system regularizes its motion under addition of damping \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b >0$$\end{document} ^8,21^, and we will explore and extend this result.
K is the only parameter controlling the Hamiltonian system dynamics, and thus the phase-space structure. For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K\lesssim 0.97$$\end{document} , the entire phase-space is regular, resembling that of a simple pendulum, while for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K\gtrsim 6.75$$\end{document} it is entirely chaotic. Between these two values, chaotic and regular regions coexist, and trajectories never cross from one to the other. This drastically changes as soon as we add any amount of damping \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b>0$$\end{document} , which now causes all trajectories to collapse to islands of regularity in the steady-state (after sufficient time). Further addition of thermal noise \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T>0$$\end{document} will diffuse the steady-state density around these islands, much like around energy wells at thermodynamic equilibrium (for Fig. 1a, we used \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=2,\, b=0.01,\, T=0.001$$\end{document} ).
A simple way to modify this system is to spring-load the rotor, thus adding a linear restoring force \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-k\,\theta$$\end{document} in Eq.2, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k =$$\end{document} spring constant (Fig. 1b). This system is known as Zaslavsky Web Map, and has been studied in the context of anomalous diffusion and quantum localization, with recent uses in cryptography^22,23^. This map is now characterized by 2 parameters, K and f – the natural frequency of the oscillator (so \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k=(2\pi f)^2$$\end{document} , and we take kick period to be the unit of time). As before, while K controls the transition to chaos, f determines the symmetry of the phase-space structures. So for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f=1/4$$\end{document} or 1/6, we get the entire phase-space tiled with square or hexagonal lattice, respectively, of islands of regularity, making a sort of “web” (hence the name). It is perhaps surprising that islands of regularity are still found for arbitrarily irrational values of f, and take various beautifully unexpected shapes (for sufficiently low K)^19^. This illustrates the robustness of mixed chaos to complex conditions, and that regular dynamics may adapt to become quite complex themselves in such cases – thus leading to rich dynamically ordered behaviors. As before, across all these regimes, our conjecture 1 was seen to hold in this system as well (for Fig. 1b, we used \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=3, f=1/6, b=0.01, T=0.001$$\end{document} ).
To illustrate the conjecture beyond kicked-rotor-like or other discretely driven systems, we show it for the Duffing oscillator (Fig. 1c):
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&\dot{x} = v \nonumber \\&\dot{v} = -x^3 + x + F \,\sin (\omega \, t) - b\, v+ \sqrt{2\,b\,T} \; \xi (t) \end{aligned}$$\end{document}This is similar to a driven harmonic oscillator, except we replace the quadratic potential, with the bistable \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$U(x) = x^4/4 - x^2/2$$\end{document} (Fig. 1c-i). Simulating this system is slower than the kicked systems, and requires full numerical forward-integration using sufficiently small time-steps \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$dt=0.01$$\end{document} . Being a simple non-linear generalization of the harmonic oscillator, while exhibiting a rich phenomenology including mixed chaos, strange attractors (above a critical value of b), stochastic resonance, and hysteresis, this system has been a key bridge between analytically tractable theory and key real-world complex behaviors^9,24,25^. The trajectories shown use the widely studied regime \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F=0.3,\; \omega =1.2$$\end{document} , visualizing the stroboscopic Poincare map at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega \, t = 0$$\end{document} (mod \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\pi$$\end{document} ). When coupling to a thermal bath via \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b=0.01,\; T=0.01$$\end{document} , we again see emergent order through collapse to the regular island. Note that in the Hamiltonian trajectories Fig. 1(c-iii), we see trajectories outside the central chaotic region that seem like they could be regular – however plot (c-ii) reveals that these still have a positive local Lyapunov exponent, unlike the truly integrable island of regularity in the middle.
Another important and interesting example verifying the applicability of our conjecture 1 is the bouncing ball system, well-known as a simple experimental demo of rich chaotic dynamics^26,27^, with recent applications like energy-harvesting^28^. Here an elastic ball is bouncing vertically (1D motion) in a uniform gravitational field on a sinusoidally oscillating plate:
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&Y_{plate}(t) = A\, \sin (\omega \, t) \nonumber \\&\dot{y} = v \nonumber \\&\dot{v} = -g - b \,v + \sqrt{2\,b\,T}\;\xi \nonumber \\&v \rightarrow -v + 2\,\dot{Y}_{plate} \quad \text {at bounce} \end{aligned}$$\end{document}Note that the thermal bath coupling here can be added via air-resistance (as in these equations) or via a coefficient of restitution and noise at each bounce. For small bath couplings, both methods produce the same results. In Fig. 1d, we plot these dynamics stroboscopically at the point of bounce (so, unlike other systems, not at equal time-steps), showing the ball’s post-bounce speed and plate’s phase in the oscillation cycle at that moment. This is done for convenience, as the fastest way to simulate this system is at every bounce to numerically search for the time of the next bounce. It has previously been observed that addition of dissipation (like air resistance or energy loss at bounce) leads to regular dynamics, which are also “quantized” into discrete energy levels^27,29^. Our conjecture explains this phenomenon, showing that such quantization merely reflects the phase-space structure of the non-dissipative Hamiltonian system, just as we saw with the kicked rotor earlier (Fig. 1d – where the plate oscillates at 1 Hz, with 0.1 m amplitude, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g =9.8\, m/s^2$$\end{document} , and bath coupling is introduced via \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b=0.03$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T=0.01$$\end{document} or via coefficient of restitution = 0.98).
An interesting practical application we hypothesize for this result is in understanding and controlling the formation of corrugation, or “washboard effect,” on dirt roads. In the reference frame of a car’s wheel, we suggest that we could model it as a ball bouncing on the oscillating surface of the passing road. Then the spontaneous selection of the bouncing regular states we see in Fig. 1d would lead to a formation of entrenched corrugation. If so, these could therefore be disrupted by driving the system back into chaos, e.g., by reducing dissipation or introducing additional disorder – perhaps in the road surface material, such as by using a gravel mixture with large-scale heterogeneity. Such interventions may not be obvious in the usual way of modeling these corrugations, which is to view them as a linearized dynamical instability or pattern formation^30^. This serves as a good example of a scenario where chaos may be preferable to the dynamical order arising from bath coupling.
Finally, we check that conjecture 1 also holds for systems with many degrees of freedom, pointing to its relevance for many-body self-organization. We consider a natural generalization of the web map which we term “Kicked Harmonic Net” (KHN) (Fig. 1e), which as far as we know has not been studied before. For N kicked rotors, in addition to pinning each one by a spring to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta =0$$\end{document} as in the web map, we also connect them among each other in a harmonic oscillator network, thus giving the dynamics:
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{pmatrix} \vec {\dot{\theta }}\\ \vec {\dot{v}} \end{pmatrix} = - \overbrace{\begin{pmatrix} 0 & -I\\ k\,I + \Lambda & b\, I \end{pmatrix}}^{W}\, \begin{pmatrix} \vec {\theta }\\ {\vec{v}} \end{pmatrix} \;\;-\; \begin{pmatrix} 0\\ K\,\sin (\vec {\theta })\;\end{pmatrix}\delta (t-n) + \sqrt{2\,b\,T}\;\begin{pmatrix}0\\ \vec {\xi } \end{pmatrix} \end{aligned}$$\end{document}Here \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\vec {\theta }$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\vec{v}}$$\end{document} are now vectors of the N rotor coordinates, I is the identity matrix, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda$$\end{document} is the connectivity graph Laplacian. This way, despite the vast complexity and flexibility of this system, it is still very fast to simulate as we have steps of linear evolution given by the matrix W followed by non-linear kicks. Note that kick amplitude and phase are here taken to be identical for all rotors.
The parameter space here is clearly too vast to survey fully, and so for a representative example, in Fig. 1e we took a randomly connected network of 16 rotors, where an edge was included in the network with probability 0.1 for every node pair (Erdös–Rényi model – specific connectivity shown in Fig. 1(e-i), and parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=0.4, \,k=(\pi /5)^2,\, b=0.005,\, T=0.001$$\end{document} , with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda$$\end{document} scaled such that its mean eigenvalue is 2 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Rightarrow \Lambda _{edge}\approx -0.84$$\end{document} ). Since the phase space is 32-dimensional, we plot the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\theta _1,v_1)$$\end{document} phase space for just one of the oscillators. Note that here, the islands of regularity are so highly fine-tuned that no randomly sampled initial condition lands in them. Nonetheless, coupling to the bath reveals the presence of intricate set of ordered or partially ordered attractors—Fig. 1(e-iv). To check that these attractors correspond to islands of regularity of the Hamiltonian system, we measure the local Lyapunov exponents at a set of points sampled from the attractors, as compared to points sampled randomly – plot (e-ii) (note that we used a black background for the plot to highlight that unlike in other plots in column (ii), the phase space here is not sampled exhaustively). Similarly, in (e-iii), we ran Hamiltonian dynamics from set of IC combining random points with points from the attractors – interestingly, only 2 of the 4 attractors in (e-iv) persisted under Hamiltonian evolution.
It’s important to note that unlike the other plots in column (iv), the dynamical order in the attractors here is not necessarily fully deterministic. This might at first seem contradictory to our conjecture 1, for which islands of regularity were seen as integrable, and thus deterministic, dynamics. However, this behavior makes sense in the context of many-body self-organization such as the flocking of birds, where partial order is possible, with some birds forming a cohesive flock, while others are still moving independently and chaotically. In this scenario, it would be unwise to restrict our attention only to global order. Similarly, we can imagine a KHN where the springs pinning each rotor to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta =0$$\end{document} are dominant compared to the ones among different rotors: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k \gg \Lambda$$\end{document} (see Eq. 5). In the limit \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda =0$$\end{document} the rotors are independent, and clearly some can be dynamically ordered, while others are still chaotic, yielding chaotic motion on a lower dimensional submanifold of the full 32-dim. phase space. In this case, to find precisely the degree of partial order we could measure the fractal dimension of the trajectory. As we gradually restore the connectivity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda >0$$\end{document} , the different rotors’ attractors start to interact. For large \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda$$\end{document} as in our example here, the dynamics no longer decompose by individual rotors, but some non-obvious partial order may still emerge – and does emerge, as we see here. To quantify this, we measure the fractal dimension of each attractor (using the “correlation dimension”), as labeled under plot (e-iv).
While the rich phenomenology of KHN in its vast parameter space remains largely unexplored, our conjecture 1 continuous to hold in all the cases we tested. We observed all kinds of different complex attractor structures arise from weak coupling to thermal bath, including complex 1D knots, some with periods of hundreds of kicks before they repeat, 2D curved manifolds with interesting topologies, and higher dimensional attractors more difficult to visualize. This qualitatively continues to hold whether we change values of parameters like K or k, number of rotors in the network, or the network topology (from random, to lattices, to small-world, etc). We invite further exploration of this system and its emergent behaviors in future work.
Another interesting observation about this system is that it maps closely to Recurrent Neural Network architecture, where a weight matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text{ exp }\left[ -W\right]$$\end{document} is being iteratively applied to the state vector \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\vec {\theta },{\vec{v}} )$$\end{document} between nonlinear activation function given by the kicks (see Eq. 5). Providing input via time-varying kick strengths \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\vec {K}(t)$$\end{document} or via the initial conditions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\vec {\theta }_0,{\vec{v}} _0)$$\end{document} , one forward-pass thus corresponds to system evolution, while learning – to changing the connectivity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda$$\end{document} . Emergence of a discrete set of dynamically ordered states implies an inherent tendency for compression and classification of the input signal in this architecture, even before any W learning. This can give more robust outputs and possibly better resistance to adversarial attacks, controlled by the bath coupling b. At the same time, KHN may be a physically-realizable learning architecture, if we appropriately instantiate the dynamics of spring strengths on a slow timescale (see^31^). This builds on the recent interest to leverage self-organization in machine learning^32^.
Transient
While conjecture 1 is powerful in its own right, its practical usefulness depends not only on the steady-state properties, but also on the relaxation timescale. We mentioned that emergence of order should happen for arbitrarily small b, but will take correspondingly longer. We sought to quantify this dependence with numerical experiments, focusing on the simple kicked rotor as it’s the easiest to work with numerically.Fig. 2. Regularization timescale \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} and its dependence on damping b. (a) kicked rotor relaxation dynamics (100 timesteps per frame at 3 points in evolution, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=2, b=0.002$$\end{document} ); (b) corresponding relaxation for the fraction of random initial conditions that have regularized – defining the timescale \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} ; (c) power-law dependence of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} on b, which defines exponent \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} ; (d) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} vs. kick strength K shows no simple pattern, but falls neatly along straight line when plotted against \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon$$\end{document} – defined as the volume fraction occupied by islands of regularity in phase space (small for larger kick strength K). This defines another scaling exponent \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi$$\end{document} , here \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=0.4$$\end{document} .
To measure the regularization timescale we first need a way of detecting when a configuration is chaotic and when it’s regular. For this, we used two methods. The simpler one is simply to count how many unique configuration were visited out of the past 20 steps, and if it is fewer than 16, then we roughly know that this is not chaos. This “counting” method will, however, fail to identify regular motion at non-zero temperature, as well as regular precession, such as in 1(e-iv, red attractor). Therefore, for our second method, we directly measure local Lyapunov exponents by slightly perturbing the configuration, evolving both versions for 10 steps, finding a linear fit to log separation, and thresholding on the resulting slope. This method can cover the failure modes of the first one, but can get confused for weak chaos – such as just outside regular islands. Therefore we use one or the other method depending on the specific context – and verify that the two agree when both are applicable.
In Fig. 2b we see the resulting curve showing the probability of regularization over time at finite temperature. It is characterized by two parameters – the regularization timescale \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} , which we study here, and the steady-state probability of dynamical order p, which we consider in the next section and Fig. 3. To recover the Hamiltonian regime, we expect that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} should diverge as damping \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b \rightarrow 0$$\end{document} – and in Fig. 2c we empirically show that this divergence is a power-law \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau \propto b^{-\eta }$$\end{document} . While the power-law behavior is universal across system parameters and the different systems we tried here, the scaling exponent \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} itself is not. In Fig. 2d, we thus show how \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} varies with the kick strength K. We see that this dependence is complex and non-monotonic. At the same time, K is a parameter specific to the kicked rotor, while to find a general law, we want to use some universally relevant feature of the system. We hypothesize that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} is primarily sensitive to the relative size of the islands of regularity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon$$\end{document} – i.e., the “degree of fine-tuning” required for dynamical order. We define \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon$$\end{document} as the fraction of randomly sampled initial conditions whose evolution is regular (regularity here measured by the “counting” method above). This way, for a set of kick strengths K, we can measure both \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon$$\end{document} , and plot them against one another—Fig. 2e. We see plotting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} against \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon$$\end{document} instead of K neatly collapses all measurements onto one line. Although \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} does not vary over enough decades to see if this dependence is a power-law, at least locally it fits \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta \propto \epsilon ^{-\chi }$$\end{document} with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi = 0.4$$\end{document} .
These results are important for several reasons. First, for regular inertial damped motion \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dot{v}=-b\,v$$\end{document} , the relaxation timescale \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau = 1/b$$\end{document} , and so \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta =1$$\end{document} . As we see in Fig. 2e, this value is approached when the islands of regularity become “macroscopic,” taking up a substantial fraction of the phase space ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon \gtrsim 0.05$$\end{document} ). In this regime, regularization is “trivial” in the sense that it is just a matter of damping out extra velocity to relax to the nearest attractor – since little fine-tuning is needed. Away from this limit, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta >1$$\end{document} indicates non-trivial regularization dynamics. Such anomalous behavior can arise due to the “stickiness” of chaotic trajectories near regular islands, possible anomalous diffusion dynamics in the chaotic regions, or due to some kind of a “search” for the rare fine-tuned regular states.
Second, since the relation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta (\epsilon )$$\end{document} in Fig.2e does not refer to any system-specific quantities, it could be universal. We were able to confirm that it roughly holds for the web map (Fig. 1b), but could not reproduce it for kicked harmonic net (Fig. 1e). There, in addition to the heavier computation required, the high-dimensional phase space makes measuring \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon$$\end{document} variations harder. Preliminary experiments showed that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta \approx 5$$\end{document} seems to be relatively stable across parameters and even different number N of rotors.
Such modification of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta (\epsilon )$$\end{document} relation in many-body systems could be attributed to the possibility of partial order, allowing a more gradual approach to regularity (see^33^). The simplest example to understand this is if instead of seeing Fig. 1b as N different initializations of one kicked rotor, we saw it as describing a KHN with N rotors and vanishing connectivity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda \rightarrow 0$$\end{document} . If we identify dynamical order by thresholding on Lyapunov exponent as we do now, then for this 2N-dimensional system, this threshold would be crossed when enough of the rotors reach regular behavior – i.e., when the regularization probability p in Fig. 1b goes above some threshold. For large N this will happen almost deterministically due to self-averaging, giving an apparently sharp onset of order with timescale \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau \propto b^{-\eta }$$\end{document} same as for a single rotor, independent of N. In contrast, the island area fraction in this 2N-dimensional space will be \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon _{N} = \epsilon _1^{\;N}$$\end{document} . So if we compare such ensembles with different N, we might have varying \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon _N$$\end{document} , but constant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} – similar to our preliminary observation for KHN. While this gives a reason for why \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} might vary differently in a many-body system, it does not explain why \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} should be constant across parameter regimes. This observation might indicate some many-body universality, where dynamics are dominated by collective modes arising from strong chaos and self-averaging, and insensitive to details of individual island structure. These ideas needs further development, and we suggest that rather than being the full picture, our conjecture may be a fundamental building block with which we can build a more general theory of many-body self-organization.
Finally, this relation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta (\epsilon )$$\end{document} is important as it quantifies the usefulness of our main conjecture 1 in practically interesting cases, which may have strong fine-tuning (small \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon$$\end{document} ), such as for many-body self-organization. So if some scaling like \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta \propto \epsilon ^{-0.4}$$\end{document} in Fig. 2e was general, that would indicate that for extreme fine-tuning, this mechanism will have a very large \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} , basically corresponding to random search of globally ordered states, which would predict unrealistically long self-organization timescales unless b is large. This way, the fact that we see \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta (\epsilon )$$\end{document} being modified in our many-body example is promising evidence that the mechanism of conjecture 1 may be on the right track to explain cases of interest.
Noise
We said that conjecture 1 requires a low-temperature bath, hence we now want to quantify how raising the temperature breaks down emergent order. Fig. 3a shows how dynamical order is destroyed at high thermal temperatures, giving a “phase diagram” in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K-T$$\end{document} space (again we focus on kicked rotor here for simplicity). While it is tempting to think of this as an order-disorder phase transition, we must be cautious to remember that thermodynamic phase transitions are defined in terms of order parameters capturing the collective behavior of many degrees of freedom ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\rightarrow \infty$$\end{document} ), while here we are talking about a single kicked rotor dynamics. Nonetheless we can still loosely think of the two regimes as phases if we think of trajectories as a 1D chains along the time-dimension, and so the long-time average of the steady-state becomes our \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\rightarrow \infty$$\end{document} limit. In that case, all the familiar properties of phase transitions can be studied, such as diverging correlation length, scaling exponents, etc.Fig. 3. Effects of thermal noise (on kicked rotor). (a) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K-T$$\end{document} phase diagram at 3 values of b, color showing self-organization probability p at steady-state (see Fig. 2b). (b) phase boundary movement with b: the transition temperature \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T^*$$\end{document} is observed to depend on b as power-law, defining scaling exponent \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} . (c) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} depends on system parameter K, and can be related to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} (recall from Fig. 2d, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} is large when K is). Line shows analytical prediction \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = \eta -1$$\end{document} – suggesting an underlying fluctuation-dissipation relation (in equilibrium systems, Einstein relation gives \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma =0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta =1$$\end{document} ). b is given in units of 1/kick, and T – \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(rad/kick)^2$$\end{document} .
Next, we check how the transition temperature \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T^*$$\end{document} depends on damping b (Fig. 3a, b). Recall that for equilibrium systems in an energy landscape, the steady-state probability is independent of b, which only controls the transients. Mathematically, this is ensured by the Einstein relation, which sets the noise amplitude to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{2\,b\,T}$$\end{document} (see Eq. 2). Thus the observation that here \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T^*$$\end{document} depends on b at all is already interesting, indicating that we cannot simply view islands of regularity as acting effectively like energy wells. Moreover, we find approximately a power-law relation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T^* \propto b^\gamma$$\end{document} , with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} – a system-specific scaling exponent (Fig. 3b).
Interestingly, it seems that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T^*$$\end{document} indeed becomes independent of b ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma \approx 0$$\end{document} ) in the same regime where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau \propto 1/b$$\end{document} ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta =1$$\end{document} ) – i.e., where the relaxation dynamics are diffusion-like. Since Einstein relation is derived precisely by balancing transient relaxation with stochastic fluctuations in the steady-state Boltzmann distribution, and yields \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma =0,\, \eta =1$$\end{document} , we suspect that we can generalize this derivation to relate exponents \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} in a nonequilibrium fluctuation-dissipation relation here.
A simple way to model this relation is to approximate the system as a 2-state Markov process and consider the transition rates R between chaotic c and regular r dynamics. Such a simple model is motivated by realizing that in the undamped Hamiltonian system, islands of regularity r and sea of chaos c are isolated from each other, and so for small bath coupling b, we get perturbatively slow transition rates \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim O(b)$$\end{document} ^7^. This means that we can approximate the dynamics internal to c and r to erase memory of IC faster than the rare transitions among r and c can happen, justifying the Markov assumption. Then, in analogy to equilibrium phase transitions, which arise from competition between the concentrating effects of dissipation and the disordering effects of entropy, here \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{c\rightarrow r}$$\end{document} captures the former and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{r\rightarrow c}$$\end{document} – the latter. Our measurements of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} in Fig. 2 empirically give that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{c\rightarrow r} \propto b^{\eta }$$\end{document} . To get \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{r\rightarrow c}$$\end{document} , we approximate the dynamics inside islands of regularity to be freely diffusing along v-dimension. This gives the exit rate in terms of the first-passage time of a diffusing particle in 1D to leave a region of size a – which is known to be \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\propto a^2/D$$\end{document} . Since our phase space is 2D, we estimate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a^2 \propto \epsilon$$\end{document} , thus giving \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{r\rightarrow c} \propto D/\epsilon$$\end{document} . From Einstein relation of the thermal bath, which drives the diffusion here, we get the diffusion coefficient \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D = b\,T$$\end{document} , and so \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{r\rightarrow c} \propto b\,T/\epsilon$$\end{document} . In the steady-state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_c \, R_{c\rightarrow r} = p_r\, R_{r\rightarrow c}$$\end{document} , from which we get that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_c = p_r$$\end{document} at temperature \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T^* \propto \epsilon \,b^{\eta -1}$$\end{document} , thus predicting that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = \eta -1$$\end{document} . Figure 3c shows that this roughly agrees with the data.
Being a discrete analogue of the derivation of Einstein relation, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b^{\eta }$$\end{document} played the role of motility \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document} , while diffusivity D was given by the usual expression here \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b\,T$$\end{document} , we thus get that
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} D/\mu = b^{1-\eta }\, T \end{aligned}$$\end{document}– a generalized Einstein relation, which reduces to the usual one for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta =1$$\end{document} .
This simple model is clearly not the full story, but has the potential to be general. We were able to verify this relation also for the web map system, but not for the kicked harmonic net. There, the values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta \approx 5,\, \gamma \approx 1$$\end{document} seem to robustly persist across parameter regimes and different numbers of rotors – even though the above derivation would still expect \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = \eta -1$$\end{document} , regardless of the phase space dimensionality. This discrepancy makes sense when we remember that in the many-body case, the attractors we observed were only partially ordered, and therefore not deterministic (see Fig. 1e-iv). This means that the diffusivity inside the nearly-regular islands is not given purely by the thermal bath \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D=b\,T$$\end{document} as before, but can experience “noise amplification” when chaos stretches small thermal fluctuations. This modifies the exit rates \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{r\rightarrow c} \propto D/\epsilon$$\end{document} , and hence the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma \leftrightarrow \eta$$\end{document} relation. To further refine the above model, one could also include a linear restoring “force” \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-b\,v$$\end{document} in the first-passage time calculation, and relax the Markov assumption – an undertaking we leave to future work.
Conclusions
In this work we used five different in silico dynamical systems of varied structure and dimension to test the conjecture that Hamiltonian dynamics exhibiting mixed chaos will settle into islands of regularity upon coupling to a thermal bath at sufficiently low temperature. The practical power of this conjecture comes from showing that if a system is at all capable of regular motion, then cooling it by coupling to a cold thermal bath will tend to induce that regular motion regardless of initial conditions. Note that this means that cooling (or damping) does not generally reduce motion, but rather stabilizes one of the possible ordered dynamical patterns. Finally this conjecture allows making quantitative predictions about stability, relaxation times, and phase diagrams of the dynamically ordered behavior based on the knowledge of the phase space structure of the original Hamiltonian system. This can allow new ways to control dynamical systems by engineering the regular states and modulating their relative noisiness (similar to^1^ for example). While we do not present a formal proof of conjecture 1, leaving it for future work, we can now give some physical motivation for it.
Our conjecture represents a bridge between two very different perspectives on emergence of dynamical order. In the thermodynamic perspective, the emphasis is on energy dissipation and lowering temperature as drivers of entropy reduction and emergent order (e.g., eq.1). In dynamical system theory, the focus is on damping, which allows contracting phase-space volumes, erasing memory of initial conditions, and turning small local Lyapunov exponents negative. But temperature and damping are two aspects of the same molecular dynamics of a thermal bath. Einstein’s relation first established this connection between the thermodynamic and dynamical aspects of fluids. This way, while dissipative drag in undriven systems contracts phase-space volumes, turning all local energy minima into zero-dimensional dynamical attractors, thermal fluctuations add stochastic exploration of the phase space, together realizing the entropy-maximizing Boltzmann distribution. Here we want to generalize this insight for nonequilibrium systems (see eq.6).
The analogy in driven Hamiltonian dynamics exhibiting mixed chaos arises because the dynamical counterpart of energy minima turns out to be the islands of regularity in the undamped driven system’s phase space – which similarly turn into absorbing attractors as soon as we add drag. One way to understand this is in the context of “low rattling” theory^1,34^: Initial conditions in the chaotic region execute a diffusive search of phase space until, at random, they enter a region of regularity. Once in such a region, and if temperature is low, the diffusive search shuts down due to the drop in chaotic motion, and density thereby accumulates around closed, ordered trajectories. This contraction of phase space density into lower sub-volumes is permitted by the frictional drag’s violation of Liouville’s theorem, which originally had to be obeyed by Hamiltonian dynamics in the absence of the heat bath. The key assumption of low rattling required for this argument is that steady-state probability of a state is controlled predominantly by that state’s exit rate – see^34^ for a general study of this assumption.
This way, even for very small damping b we still expect self-organization in the steady-state, though the relaxation time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} can be long: we found that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} scales as power law \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\propto b^{-\eta }$$\end{document} (Fig. 2). The scaling exponent \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} depends on system parameters, and may depend on the degree of fine-tuning required for dynamical order \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon$$\end{document} (defined as volume fraction of regular islands in the original Hamiltonian system) as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta \propto \epsilon ^{-\chi }$$\end{document} . We saw that in many-body systems, where the possibility of partial order makes self-organization more gradual, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} can be less sensitive to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon$$\end{document} . Furthermore, as we vary temperature, our theory predicts a smooth breakdown of regularity, characterized by transition temperature \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T^*$$\end{document} (Fig. 3). We find that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T^*$$\end{document} depends on b, indicating the nonequilibrium nature of the phenomenon, and that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T^* \propto b^\gamma$$\end{document} . At equilibrium, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta =1,\,\gamma =0$$\end{document} . A simple model based on above understanding suggests a more general fluctuation-dissipation-type relation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma =\eta -1$$\end{document} , which we validate experimentally. Interestingly, this relation is modified in systems with more degrees of freedom, indicating that other mechanics may be taking over, such as modification of free diffusivity inside the more complex ordered regions. Because of the broad implications this has for understanding nonequilibrium emergent order, we expect the conjecture on the effects of thermal bath coupling on mixed chaos presented here to stimulate further research across a range of disciplines.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Delande, D. Kicked rotor and anderson localization. Boulder School on Condensed Matter Physics (2013).
- 2Mohammed, A. H., Shibeeb, A. K. & Ahmed, M. H. Image cryptosystem for iot devices using 2-d zaslavsky chaotic map. Int. J. Intell. Eng. Syst. (2022).
- 3Barroso, J. J., Carneiro, M. V. P. & Macau, E. Bouncing ball problem: stability of the periodic modes. Physical review. E, Statistical, nonlinear, and soft matter physics 79 2 Pt 2, 026206 (2009).10.1103/Phys Rev E.79.02620619391819 · doi ↗ · pubmed ↗
- 4Huang, C. et al. Research on dynamics of bouncing ball in triboelectric nanogenerator. Journal of Micromechanics and Microengineering 31 (2021).
- 5Vargas, D. V., Foong, T. Y. & Zhang, H. Dynamical equations with bottom-up self-organizing properties learn accurate dynamical hierarchies without any loss function. ar Xiv preprint 10.48550/ar Xiv.2302.02140 (2023). Ar Xiv:2302.02140 [cs].
