Controlling a complex system near its critical point via temporal correlations

Many complex systems exhibit large fluctuations both across space and over time. These fluctuations have often been linked to the presence of some kind of critical phenomena, where it is well known that the emerging correlation functions in space and time are closely related to each other. Here we test whether the time correlation properties allow systems exhibiting a phase transition to self-tune to their critical point. We describe results in three models: the 2D Ising ferromagnetic model, the 3D Vicsek flocking model and a small-world neuronal network model. We demonstrate that feedback from the autocorrelation function of the order parameter fluctuations shifts the system towards its critical point. Our results rely on universal properties of critical systems and are expected to be relevant to a variety of other settings.

Scientific RepoRtS | (2020) 10:12145 | https://doi.org/10.1038/s41598-020-69154-0 www.nature.com/scientificreports/ We show that a system can be tuned to the vicinity of its finite-size "critical" (i.e. maximum susceptibility) point using the first autocorrelation coefficient AC(1) of the order parameter fluctuations (we define AC(1) as in time-series analysis as the time correlation of the order parameter at t = 1 ). Because AC (1) peaks at the same point as the susceptibility, yet does so more smoothly than the susceptibility, a control feedback seems straightforward. This behavior of AC (1) near criticality responds to the notion of critical slowing down in both equilibrium and non-equilibrium critical dynamics by which perturbations take longer to dissipate near criticality.
The existence of the maximum of AC (1) can be understood from dynamic scaling. The dynamic scaling form of the time correlation is 15 where k is the observation wavevector, ξ is the correlation length, the function g is such that g(t = 0) = 1 , i.e. C 0 (k) is the static correlation function, and the characteristic time obeys where g, Ŵ and are unspecified scaling functions and z is the dynamic scaling exponent. Now For a global quantity, k = 0 (see Suppl. Material) so that where A > 0 is a time dependent constant. Hence, the normalized time correlation has a maximum at T = T c , for fixed δt.
The main idea is demonstrated here by applying it to three well understood systems, namely the ferromagnetic Ising model, the Vicsek model of flocking and a typical neuronal small-world network. We remark that the results are general enough to be also expected in many other systems.
Ising model Figure 1 illustrates the typical behaviour of the 2D ferromagnetic Ising model at increasing temperatures. The system undergoes a second order phase transition at a critical temperature T c , reflected in a steep change in magnetization as well as a sharp peak in susceptibility (Fig. 1A). Equally distinct changes are also demonstrated for the correlation properties of the model computed from appropriate system variables (Fig. 1B). A sharp increase in the average pairwise correlations is observed as the system approaches T c , where the correlation length matches the size of the system. The relatively sharp changes in the spatial correlations contrast with the relatively smoother changes in the temporal correlations, as reflected by the first auto-correlation coefficient AC(1) of the magnetization fluctuations around the mean, which at T c approaches unity. Now we asses how to control the Ising model to stay at the vicinity of the susceptibility peak. According with the discussion in the introduction, we must restrict ourselves to do it using only either local or global information. In that sense, the time correlation evaluated by AC(1) meets such conditions, because it can be computed from a temporally delayed version of a global average of magnetization. In turn, magnetization can be assessed simply by averaging samples of a relatively large number of sites.
To demonstrate control we proceed by choosing an initial random temperature and simulate the dynamics for some large number of Montecarlo (MC) steps, which we denote as an "adaptation iteration step" indexed by i. We proceed by estimating the AC(1) of the fluctuations around the mean magnetization during the lapse of time corresponding to the adaptation iteration step i and monitor the change of AC(1) between two consecutive steps i, defining so that d changes sign when a decrease in AC(1) is detected. We then use the gradient to its maximum value to change the future temperature T(i + 1) (i.e., the control parameter) according to where κ is a constant that determines how slowly the temperature is adjusted. Its exact value is not crucial for the present results. Successive iterations of Eqs. 5-7 demonstrate convergence of the temperature to the expected value at equilibrium T c ∼ 2.3 . Fig. 2 illustrates typical results for various initial temperatures, which in all cases converge to the vicinity of T c . We note that the successive values of the parameters (order, control and AC(1)) obtained during the adaptive simulations over-imposes well (i.e., matches) those obtained from equilibrium simulations.
Vicsek model We were also able to use the AC(1) function to control the Vicsek model 16 , the archetypal model for flocking behavior, towards its critical point. In this model, N self-propelled particles endowed with ; kξ , where S i is a sphere of radius r c centered at r i (t) . The operator R η normalizes its argument and rotates it randomly within a spherical cone centered at it and spanning a solid angle η� d , where d is the area of the unit sphere in The order parameter, which measures the degree of flocking, is the normalized modulus of the average velocity 16,17 , in the the ordered phase. We choose t = r c = 1 , noise amplitude η = 0.5 , the speed v 0 and the number density ρ = N/V , where V = L d is the volume of the (periodic) box. Here, we choose N as the control parameter noting that similar results are obtained using noise amplitude η as control parameter for fixed N (see Suppl. Material). We apply Eqs. 5-7 to this model, using η as control parameter and keeping the density fixed. For comparison we over-plotted results from equilibrium runs with values taken during adaptive control of the simulations (Fig. 3). The close match demonstrates that the technique is able to control the flock model near its critical size N ∼ 560.
Neuronal network model Successful control was further demonstrated for a previously described neural network model 18 consisting of a network of interconnected nodes together with a dynamical rule. The model exhibits a second order phase transition 20 on a region of parameters. The model matrix of interactions follows a small-world topology and each node exhibits discrete state excitable dynamics, following the Greenberg-Hastings www.nature.com/scientificreports/ model 19 . Briefly, each node is assigned one of three states: quiescent Q, excited E, or refractory R, and the transition rules are: (1) Q → E with a small probability r 1 ( ∼ 10 −3 ), or if the sum of the connection weights w ij with the active neighbors (j) is higher than a threshold T h , i.e., w ij > T h and Q → Q otherwise; (2) E → R always; (3) R → Q with a small probability r 2 ( ∼ 10 −1 ) delaying the transition from the R to the Q state for some time steps. Parameters r 1 and r 2 , which determine the time scales of self-excitation and of recovery from the excited state, respectively, were kept fixed and T h was updated according to control Eqs. (5-7). The density of excited nodes, (i.e., in state E) in each time step was taken as the order parameter. As shown for the previous models, the feedback of the AC(1) of that order parameter was able to move and maintain the system near its critical point (here T h ∼ 0.16 ) (Fig. 4). It is important to remark, that the neurons excitability is adjusted not by the network rate of activity, but by the temporal correlations (i.e., AC(1)) of such activity fluctuations.
We have verified that the present results apply, with some small differences, to systems undergoing either 1st or 2nd order phase transitions. We note also that low dimensional dynamical systems exhibiting continuous or discontinuous bifurcations from fixed points to limit cycles which can be controlled, using the same idea, near the bifurcation point. www.nature.com/scientificreports/ In conclusion we have demonstrated, in three paradigmatic cases, how to shift the system towards its critical point using a feedback loop between the control parameter and the first autocorrelation coefficient of the order parameter fluctuations. Our results build at least on two previous lines of work which come close to describe this control strategy. One is the view of self-organized criticality 1 as a feedback between order and control parameters 21 . The other line relates to forecasting of an upcoming tipping point via the generic slowing down present at criticality 22,23 . The current results go beyond these previous approaches by demonstrating an alternative mechanism for the presence of criticality in some systems. Furthermore it provides a strategy of control amenable of practical implementations in different areas. For instance, in neuroscience, this approach could be implemented in conjunction with optogenetical targeting 24 in behaving rodents to clamp cortical networks to any desired dynamical state, helping to predict its influence on perceptual performance during a given task.