An Overview of Emergent Order in Far-from-equilibrium Driven Systems: From Kuramoto Oscillators to Rayleigh-B\'enard Convection

Soft-matter systems when driven out-of-equilibrium often give rise to structures that usually lie in-between the macroscopic scale of the material and microscopic scale of its constituents. In this paper we review three such systems, the two-dimensional square-lattice Ising model, the Kuramoto model and the Rayeligh-B\'enard convection system which when driven out-of-equilibrium give rise to emergent spatio-temporal order through self-organization. A common feature of these systems is that the entities that self-assemble are coupled to one another in some way, either through local interactions or through a continuous media. Therefore, the general nature of non-equilibrium fluctuations of the intrinsic variables in these systems are found to follow similar trends as order emerges. Through this paper, we attempt to look for connections between among these systems and systems in general which give rise to emergent order when driven out-of-equilibrium.


I. INTRODUCTION
A system at equilibrium is indistinguishable from its surroundings. The same system when driven out-of-equilibrium gives rise to flows that forces the system to relax back into its equilibrium state. The rate of relaxation is governed by how far the system has been driven out-of-equilibrium [1][2][3]. Soft-matter systems in this respect are especially fascinating as they give often rise to order as long as the driving field maintains it out-of-equilibrium [1,4,5].
Some prominent examples where self-organization or self-assembly gives rise to emergent order include, liquid crystals, granular material, polymers, gels, and a wide spectrum of biological phenomena/materials [6][7][8][9][10][11][12][13]. This emergent order can vary across several length and time scales, and since they are very sensitive to fluctuations (thermal) they are usually difficult to predict.
In this paper, we discuss how coupling can play a role in driven systems as they selfassemble to give rise to emergent order. We start with the simplest statistical model that shows phase-transition -the two-dimensional square lattice Ising model. Following which, we model the phenomena of synchronization of a large set of coupled oscillators or the Kuramoto model. Both of these systems are numerically modelled and the effect of coupling on the emergence of order is discussed. The results are then compared with an experimental system, the Rayleigh-Bénard convection for different fluid mixtures at varying Rayleigh numbers (10 3 − 10 6 ). Our study indicates that the nature of emergent order, and the second statistical moment of the respective intrinsic variables follow similar trends across the three systems [14][15][16].

A. Ising Model
One of the earliest physical models that studied the emergence of order as a consequence of interaction between the agents is the Ising model [17,18]. Traditionally the Ising system was used to model ferromagnetism in statistical physics, where magnetic dipoles could either have a spin 'up' or 'down'. Since then it has become a prototype for many two-state model systems examples include, protein folding, ligand-receptor interactions, spin glasses, firing of neurons etc [19][20][21][22]. The Hamiltonian for an Ising system in the presence of an external field 'h' is given by, The spins are denoted by σ and the indices represent neighboring lattice sites. The signature of J ij tells us the nature of interaction between the pair (i, j). While the simplest case of the Ising system is the one-dimensional case, interesting features emerge when it is studied on a two-dimensional square lattice. The two-dimensional square-lattice Ising model is one of the simplest statistical models that allows for phase-transition [17,18,23]. In order to numerically solve the problem a two-dimensional partition function is defined, Here, σ assigns a value of either +1 (up) or −1 (down) for each lattice site and the variables 'm' and 'n' denote the rows and columns of the lattice (the special case being m = n = N) with periodic boundary conditions. For the case of isotropic coupling one achieves a phase transition when, The dynamics of the model was simulated using a Monte-Carlo algorithm to randomly assign spin values at every lattice site. One can encounter a practical problem if the spins are to randomly flip at every lattice site with each simulation step, and end up eventually with a checkerboard pattern. Therefore, one needs to create a Markovian decision model where the spin state at a site is the most probable outcome based on spin probabilities in a set of randomly chosen sites within the lattice. While better prediction based on larger regions for decision making makes the simulation faster, this was not really the aim of the model. Macroscopically, the system's dynamics is 'equilibrium-like', but microscopically spin outcomes at each lattice site is inherently stochastic. The emergence of order was further analyzed when the system was externally perturbed under conditions: h = constant and h(t) = A sin ωt (refer Equation 1).

B. Kuramoto Model
Similar to the Ising system one can model collective synchronization in a large population of oscillating elements. The Kuramoto model is a mathematical model that treats a system as an ensemble of limit-cycle oscillators described only by their phases [24][25][26]. In the simplest version of the model, each oscillator in the Kuramoto system has its own intrinsic natural frequency ω i and is coupled to every other oscillator in the system. The intrinsic natural frequencies of the oscillators are drawn from a predefined distribution, usually a normal distribution with well-defined mean and standard deviation. As the system collectively synchronizes, the different frequencies spontaneously locks to a common frequency, Ω.
The Kuramoto model has found several successful applications in condensed matter physics specially in the study of biological phenomena and active matter [25,27,28]. The governing equation for the system is given by, Here, the phase of an oscillator is given by θ i and the coupling strength by κ. Through the following transformation one can solve this nonlinear differential equation for the mean-field case, N → ∞: With R as the order parameter and φ the average phase, one can transform Equation 4 and rewrite it as, Since the oscillators are randomly oriented, the sum over all oscillator phases average to zero. Hence, Equation 6 becomes, For sufficiently strong coupling one achieves a fully synchronized state (R → 1). At a fully synchronized state all the oscillators share a common frequency while their phases may differ. Under steady-state condition (dθ i /dt = 0), the fully synchronized solution for Equation 7 reduces to ω i → Ω = κ sin(θ i ) where Ω is the common frequency of the oscillators.
In this study, the mean-field Kuramoto model was simulated on a two-dimensional square lattice. The effect of coupling strength was observed on the time evolution of the order parameter and simultaneously on the second statistical moment of the angular frequencies of the oscillators. The effect of several other types of coupling mechanisms were also studied (like distance-dependent inverse square), but are not presented in this paper. Interested readers are referred to [29].

C. Rayleigh-Bénard Convection
Finally, we discuss one of the simplest experimental setup to study pattern formation and self-organization. As a thin layer of viscous fluid is heated and convection sets in, one can observe thermal gradients on the surface of the fluid film which are stable in time. The regular pattern of convection cells are known as Bénard cells and the phenomena, Rayleigh-Bénard convection [4,[30][31][32]. To date, it remains one of the most actively and extensively studied physical system. Due to its conceptual richness, the dynamics of a Rayleigh-Bénard convection phenomena connects fundamental ideas from both non-equilibrium thermodynamics and fluid mechanics [33][34][35]. The beauty of this system lies in its simplicity, wherein a dimensionless quantity, the Rayleigh number (Ra), determines the onset of convection which is defined by, Here, the physical quantities g stands for acceleration due to gravity, β for thermal expansion coefficient, ∆T for the thermal gradient across the system, l for fluid film thickness, thermal images that can be converted into a temperature matrix. The temperature of the top layer of the fluid film is obtained by averaging over the entire exposed area. With this setup in place, two types of studies are performed: spatial and temporal. In the temporal study, thermal statistics are recorded from a room temperature equilibrium to a non-equilibrium steady-state as the system is thermally driven by regulating the power input through the heater. Whereas, in the spatial study thermal statistics are obtained once the system has reached a non-equilibrium steady-state. While the temporal study allows us to envision the evolution of order in the system, the spatial study lets us visualize how order is spatially distributed through emergent length-scales. One can find more details on the study: the experiments and the analyses here [15,16].  In this section we discuss the results from our numerical simulations and the experimental study. We are essentially looking at the possibility that these systems which are distinctly different from one another exhibit characteristics that are similar when driven outof-equilibrium. In Figure 1, we report our results from the simulation of the two-dimensional square lattice Ising model. Figure 1a plots magnetization as a function of inverse temperature (β = 1/k B T ). The magnetization (S/S) acts as the order parameter for the system.
The ferromagnetic transition happens around β ≈ 0.44 which corresponds to the Curie temperature (see Equation 3) [18]. The Ising system is then perturbed by an external field, and the evolution of the order parameter is plotted as a function of time in Figure 1b and as a function of inverse temperature in Figure 1c. In the case of no external field one can see that the spins eventually lock into a mixed state with some of the lattice sites with, say 'up' spin and the remaining with 'down' spin. Therefore, the system does not reach a fully spin 'up' or spin'down' state as seen from Figure 1b. This however is not the case when there is an external perturbation as the system eventually directs itself to the direction of the external perturbation as that is energetically more favorable. In case of a sinusoidal time varying field, some oscillations are observed because of periodic aligning and re-aligning. The smaller the temperature, the larger the β and as theory predicts we see phase transitions at sufficiently low enough temperature in Figure 1c. Therefore, for no external perturbation, the transition temperature seems to be the lowest at β ≈ 0.4. In Figure 1d, we plot the standard deviation of the spin as a function of time for the case of no external perturbation at β = 0.2. The standard deviation being a measure of fluctuation steadily decreases as order emerges in the system. This observation although being intuitive yet non-trivial makes us think whether this is a general feature in systems where order emerges as they are driven out-of-equilibrium.
In Figure 2, we plot the scaled standard deviation of the intrinsic variable and emergent order as a function of time for the Kuramoto system and the Rayleigh-Bénard convection.
The intrinsic variable in the Kuramoto system is the angular frequency of the oscillators (ω i ) which collapses to a common frequency (Ω) as the system achieves synchronization.
Similarly, in the Rayleigh-Bénard system spatially-averaged temperature ( T (t) ) plays the same role. As the system approaches a steady-state, T (t) → T ∞ , where T ∞ is the spatially-averaged steady-state temperature of the system. In Figure 2a  At time-step, t = 100 one can observe that atleast more than half of the oscillators present in the system are synchronized (from Figure 2e) and therefore one observes a sharp decline in the scaled standard deviation plot in Figure 2a. Later one can notice that as t ≥ 110 there is a sudden spike in the standard deviation as order increases further. The reason for this could be attributed to a mixture of synchronized and unsynchronized oscillators as R < 1.
As time progresses, the natural frequencies of all the oscillators approach closer to mean-field common frequency. However, due to their equally random phase orientations, some of the oscillators reach the common frequency and lock themselves in that state earlier than the other. A situation like this although reduces the standard deviation when compared to the randomized initial state it however increases the standard deviation at an instant when these two groups of oscillators start oscillating simultaneously, one with low fluctuations and the other with higher fluctuations. As one would expect, this scenario appears to last longer in the case of lower coupling strength among the oscillators because of more unsynchronized oscillators than synchronized ones at any given instant in time. Following our results from There is no well-defined order parameter in a Rayleigh-Bénard system, therefore we define one based on the thermal profile at steady-state as, Here, T (t) represents spatially-averaged temperature at any instant in time, T 0 represents spatially-averaged temperature at initial equilibrium state (room temperature) and T ∞ represents spatially-averaged temperature at a non-equilibrium steady-state. For each of the fluid samples we look at the scaled standard deviation plots as order emerges. In the case of silicone oil sample, the fluid being more stable due to its high viscosity, ν = 150 cSt and low Rayleigh numbers we see a decline in the fluctuations in Figure 2b. This decline can be mapped to the first instance when convection cells start to appear in the system.

The fluctuation reaches a minima when a number of cells have fully formed and nucleated
at the center of the copper pan. As the system has not yet reached a steady-state for atleast another ∼ 10 3 time-steps (see Figure 2f) the temperature keeps on rising and hence, the standard-deviation. With Rayleigh numbers being almost in a similar range, we see a different characteristic with glycerol as our working fluid. Glycerol being a much lighter fluid with viscosity atleast a magnitude lower than silicone oil first nucleates into convection cells which remain stable for sometime, but quickly divides into smaller cells. This twostep nucleation results into two regions of decline in the standard deviation plot as shown in Figure 2c. In Figure 2d, we look at thermal fluctuations in glycerol-water mixtures.
The standard-deviation appears to decline much faster and earlier than the earlier plots (at around time-step, t = 10 2 ), but it lasts for a much shorter duration. The reason for this appears to be lower viscosities (∼ 10 −2 cSt) and higher Rayleigh numbers for the glycerol-water mixtures. Therefore, nucleation not only happens early but also spreads at a faster rate throughout the pan. Following which, they break down into smaller and smaller domains which dissipate heat rather chaotically as the system enters a turbulent regime. One can observe this from the amount of noise in the standard deviation plots, the magnitude of which keeps on increasing with time. Moreover, the system also does not reach a steady-state (temperature graph not shown) as one can see in Figure 2h, where the order as a function of time seems to be monotonic near R = 1 rather than being asymptotic. The similarities between the Kuramoto model and the Rayleigh-Bénard convection are striking. If the extent of synchronization in the Kuramoto model is considered as a measure of order then the Rayleigh-Bénard convection also shows similar trends as it reaches a non-equilibrium steady-state. Since, matter does not leave the system, the continuity equation is preserved.
At room temperature equilibrium state, the velocity vectors are randomly oriented, therefore the net directional component of the velocity field cancels itself out. While, a steady-state leads to a well-defined (and directed) velocity field which transports heat from the bottom of the copper pan (hot) to the top layer of the fluid film (cold). Therefore, emergent order in the Rayleigh-Bénard system corresponds to synchronization of the frequencies of individual convection cells as the system reaches a steady-state temperature. Thus, ω = dθ/dt = 2u ∞ /l where l/2 is the half thickness of the fluid film and the steady-state velocity, u ∞ ∝ ∇T where ∇T is the thermal gradient across the fluid film thickness.
In Figure 3, Note that in the above equation, x = δω ⋆ . A realistic approximation to such a distribution when there are tails in the data is a Lorentzian function, top cold surface. In Figure 3c, we plot the kernel density estimates to determine the shape of the probability density function for the two experimental trials with different Rayleigh numbers. In Figure 3d, we proceed to fit the data piece-wise. We choose individual tails and fit them with a pair of Gaussian fit functions (in black) and then with a pair of Lorentzian fit functions (in red). As we can see from our plots in Figure 3d, both Gaussian and Lorentzian fits superimpose over one another. The difference between the center of the two peaks is about 0.04 units with one peaking in the positive domain and the other in the negative.
Therefore, one peak signifies the contribution of the upward plumes and the other of the the downward plumes. We are still unsure of the fact that how the fit functions from the two tails merge into one another. In some of our recent works we discuss the presence of a mixture of local equilibrium regions in the Rayleigh-Bénard convection which describes the bimodal nature of the thermal fluctuations [14][15][16]. To conclude, in the mean-field Kuramoto model the final synchronous state being unique allows for the existence of a sharply peaked delta-type distribution, which in reality is best illustrated by a Lorentizian fit. In the case of a non-turbulent Rayleigh-Bénard convection at steady-state we find that there exist two possible states due to the existence of spatial thermal gradients which are stable in time.
These stable spatial gradients lead to the emergence of two local equilibrium-like regions, fluctuations within which can be best represented by respective Gaussian distributions [16].
In Figure 4, we plot the lattice entropy as a function of time for a two-dimensional Kuramoto model with high coupling strength. We have previously seen that evolution of order in the system is inversely related to the fluctuations of the intrinsic variables. By calculating the Shannon entropy summed over every lattice site at every instant in time, we look at the relationship between entropy reduction and fluctuation decay as order emerges in the system.
Therefore, to make sure that the Kuramoto system reaches a fully synchronized state, a high enough coupling between the oscillators is chosen and the simulation is run for a very long duration. It is not surprising that at κ = 5, the model achieves complete synchronization after just 500 time-steps. The Shannon entropy is calculated from Equation 12 and is scaled by its maximum value, such that 0 < S/S ≤ 1 [21,36]. As order emerges, one can clearly see from Figure 4 that a reduction in system's entropy is accompanied by a reduction in the fluctuation of the intrinsic variable.

III. CONCLUSIONS
In this paper, we consider three systems that can be externally perturbed and driven out-of-equilibrium. While the Ising model and the Kuramoto oscillators are numerically solved, the Rayleigh-Bénard convection on the other hand was experimentally probed. The common feature of all the three systems is the emergence of order as they are driven outof-equilibrium. The Ising and the Kuramoto models self-organize by synchronizing their spins and their natural frequencies, respectively. On the other hand, spatio-temporal order that emerges in the case of a Rayleigh-Bénard convection is a result of the competing forces between viscosity and buoyancy which give rise to convective instabilities. A common observation across the three systems was that fluctuations of the intensive variables decay as order emerges, which may seem counter intuitive given that the system is out-of-equilibrium.
By virtue of being in a non-equilibrium state the fluctuations in the system should dominate and drive the system further away from equilibrium. However, that is not observed across the three systems that we study in this paper. This brings us to a more pertinent question as to how far away are these systems from equilibrium? Although we do not yet have a metric to define that but we can anticipate that these systems even when are driven out-of-equilibrium are in a quasi-equilibrium state, where the local equilibrium hypothesis is held true [37,38].
We believe this is an important observation for any physical system to show stable emergent patterns. We saw in our experimental results that for high Rayleigh numbers, we enter the turbulent regime which though give rise to structures but at the same those structures are found to be not stable in time. We imagine that a system can only give rise to stable patterns when it is not driven too far away. Therefore, existing statistical physics models for synchronization can probably hold the key to explain more complicated systems that give rise to stable structures, like the Rayleigh-Bénard system. Therefore, through this paper we attempt to draw similarities by exploring different systems that show emergent structures through detailed quantitative analysis while asking the reader to find ways to understand complex emergent phenomena through simpler models from statistical physics.