Butterfly Effect and Spatial Structure of Information Spreading in a Chaotic Cellular Automaton

Inspired by recent developments in the study of chaos in many-body systems, we construct a measure of local information spreading for a stochastic Cellular Automaton in the form of a spatiotemporally resolved Hamming distance. This decorrelator is a classical version of an Out-of-Time-Order Correlator studied in the context of quantum many-body systems. Focusing on the one-dimensional Kauffman Cellular Automaton, we extract the scaling form of our decorrelator with an associated butterfly velocity $v_b$ and a velocity-dependent Lyapunov exponent $\lambda(v)$. The existence of the latter is not a given in a discrete classical system. Second, we account for the behaviour of the decorrelator in a framework based solely on the boundary of the information spreading, including an effective boundary random walk model yielding the full functional form of the decorrelator. In particular, we obtain analytic results for $v_b$ and the exponent $\beta$ in the scaling ansatz $\lambda(v) \sim \mu (v - v_b)^\beta$, which is usually only obtained numerically. Finally, a full scaling collapse establishes the decorrelator as a unifying diagnostic of information spreading.


Introduction.
A central hallmark of chaotic systems [1] is their sensitivity to perturbations: even a small change in initial conditions leads to entirely unpredictable, and large, differences in the state of the system at later times. This is popularly captured by the "butterfly effect" [2], which in manybody systems encodes two distinct notions: first, the exponential growth of the perturbation giving rise to the notion of the Lyapunov exponent, λ, characterising the growth with time [3]; and second, information spreading in space, whereby the "effect" of the butterfly's wingbeat is felt at a distant location only with a time delay given by the "ballistic" propagation speed known as butterfly velocity, v b [4,5].
A prominent recent strand of investigation of chaos in quantum many-body systems revolves around the study of information scrambling, where Out-of-Time-Order Correlators (OTOCs) have been employed to measure the propagation of quantum chaos [6]. OTOCs can be understood as two-time correlation functions where operators are not chronologically ordered and are a simple measure of the "footprint" of an operator that spreads in space [7], thus, naturally measuring the spread of information. In chaotic systems, this quantity may grow exponentially in time, governed by a Lyapunov exponent λ [8]. Recently, OTOCs have been studied extensively as an early-to-intermediate time diagnostic of quantum chaos or information spreading in a number of different quantum models whose time evolution can be generated by different dynamics such as Floquet dynamics [9,10], random unitary circuits [11][12][13][14], and time-independent Hamiltonians such as integrable spin chains [15,16], generalised SYK models [17], diffusive metals [18] and Luttinger-liquids [19].
In that context, analogues of OTOCs have been recently developed for classical systems [20][21][22][23][24][25]. For example, for spin chains the decorrelator D(x, t) = 1 − S A (x, t) · S B (x, t) between two copies of spin configurations S A/B , which at t = 0 only differ locally by a small spin rotation, is a semiclassical version of an OTOC [20]. For Heisenberg magnets, this exhibits ballistic propagation with a light-cone structure governed by a butterfly velocity even in the high-temperature regime without long-range magnetic order but with spin diffusion [20].
A common feature observed in classical and quantum models is the exponential growth (or decay) of the OTOCs (analogues) along rays of constant velocity v = dx dt quantified by velocity-dependent Lyanupov exponents (VDLEs) [8]. For intermediate/late times it takes the scaling form Intriguingly, a unifying framework of VDLEs with λ(v) = µ(v − v b ) β captures the spatio-temporal structure of information spreading in many-body quantum, semiclassical and classical chaotic systems [1,3,5,20,26].
In this work, we provide a basic description of such information-spreading in a minimal chaotic setting: the Kauffmann Cellular Automaton (KCA). Our system choice is motivated by the fact that KCA are minimal chaotic manybody models [27] which display a rich phenomenology, including a phase transition as a function of a tuning parameter, a probability p. They exhibit universal scaling, e.g. of the directed percolation universality class [28], and lend themselves to analytical insight [29]. Additionally, KCA have a wide range of applicability [30][31][32][33][34][35][36]: initially introduced to study fitness landscapes of biological systems and gene expression [37], they now appear also in optimisation problems [38], random mapping models [39], and most pertinently in the emergence of chaos [40].
Surprisingly, despite decades of research on chaotic CA, see Ref. [41] for a review, the dynamics of chaos has only been investigated in terms of the global Hamming distance but a local diagnostic has thus far been missing. Here, we construct an OTOC analogue for KCA and explore the VDLE phenomenology. This analogue enables us to uncover the ballistic spatio-temporal structure of perturbation spreading in the chaotic phase, Fig. 1(b), in contrast to the decay of such "damage" spreading in the frozen phase, Fig. 1(a). We develop a full microscopic theory of the VDLE, recovering the functional form Eq. (1) including an analytical calculation of the exponent β. Thus, we provide the tools for describing the sensitivity of chaotic many-body systems to perturbations through the framework of VDLEs.
KCA Model and Classical OTOC Analogue. We focus on a generic dissipative dynamical system called an N K model. Concretely, a local KCA is a system of N Boolean elements σ(x, t) = ±1 which evolve in discrete time steps through rules which depend upon 2K nearest neighbours of each site in 1D [41]. Our KCA system evolves under a set of (annealed) local rules {f x,t }: (2) which are random with probability p in space and time: In a pioneering work, Derrida and Stauffer [42] showed that KCA display a chaotic-to-frozen phase transition controlled by the parameter p, see Fig. 1 inset. The two phases are distinguished by the decay or spread of localised perturbations diagnosed with the global Hamming distance, between two copies of the system σ A/B (x, t) which differ by a single inverted site in the initial state at t = 0. This measure is then ensemble-averaged over realisations with the same probability p. The distance grows linearly in the chaotic phase and decays to zero in the frozen phase, see Fig. 1.
Analogous to the classical Heisenberg chain [20], we take the classical OTOC as the local distance between two copies of the same system which only differ by a local perturbation of the initial conditions, thus leading to the decorrelator It is nothing but a local Hamming distance, which is related to the global distance by H(t) = 1 N x D(x, t). Numerical Results. In Fig. 1 we show the spatio-temporal evolution of the decorrelator D(x, t) for representative values of p below and above p c . First, in the frozen phase D(x, t) initially spreads but then decays in time and space to zero. This attenuation of the local decorrelator reflects the vanishing of the long time value of the Hamming distance in the frozen phase. Second, in the chaotic phase D(x, t) spreads with an apparent light-cone structure. Again, from the sum rule the long-time value D 0 = D(x = 0, t → ∞) within the lightcone can be identified as the order parameter associated with the phase transition, see Fig. 1 inset.
Because of the locality of KCA rules, the speed of the damage spreading measured with D(x, t) is necessarily bounded by the maximum velocity v max = K but the actual spread is slower, with the butterfly velocity v b < K. To demonstrate the presence of v b , we plot the behavior of the decorrelator along rays of constant velocity. In Fig. 2 we see that, after We compare the full D(x, t) numerical data (solid lines) and the prediction (dashed lines) of the boundary random walk model. a transient effect, there is a ray (blue lines, corresponding to v = 2.9) along which D(x = v b t, t) is constant for t > 100. This establishes the presence of a butterfly velocity v b < v max and is an efficient way of extracting v b as a function of p, see blue dots in the inset of Fig. 4.
Next, we investigate whether such behaviour follows the general VDLE phenomenology of Eq. (1). We note that in contrast to previous work on Heisenberg spin chains, the presence of Lyapunov exponents in KCA is far from obvious. In fact, for small distances one does not expect an exponential pickup as a function of time because the local perturbation between the two copies σ A/B is necessarily big because of the discrete and bounded nature of the inverted site. However, we find that for distances x far away from the perturbation and after disorder-averaging one may still observe an exponential pickup in time of the OTOC analogue. Because of the discrete nature of the variables the scaling behavior only appears in a small window around v b as the decorrelator quickly saturates to D = D 0 for velocities v < v b and quickly decays to D = 0 Again, a clear picture emerges by studying D(x, t) along "rays" of constant velocity v = x/t. As shown in Fig. 2, we observe an exponential decay of the OTOC in time, with a VDLE λ(v). For rays with v > v b , the linear decay of ln D(x = vt, t) = −λ(v)t in the long-time limit demonstrates that the decorrelator indeed decays exponentially with a VDLE λ(v).
Overall, the numerical results establish the presence of a butterfly velocity v B and validate the VDLE framework in KCA systems at long times. This is a surprising feature given the discrete, and large, nature of "perturbations" in the Boolean network, as such framework is usually observed for continuous and infinitesimal perturbations [8].
To quantitatively analyse the small active region around the wavefront set by v b where the damage actively spreads, we define for each sample the boundary of a spreading perturbation as the furthest point from the centre which differs from the unperturbed system. As shown in Fig. 3, the probability density of boundary positions P (x) approaches a Gaussian profile in the long-time limit with a width σ x ∝ √ t. We find that its mean is in quantitative agreement with the probabilitydependent butterfly velocity v b (p) for all p, see inset of Fig. 4. The Gaussian behavior of the boundary, therefore, motivates a random-walk-like description of the active region.
Boundary Random Walk Model. We can now develop a microscopic statistical model of the boundary dynamics starting from the microscopic KCA rules. We only sketch the main idea and relegate details to the Supplementary Material. The basic ingredient for our random-walk model is the expectation value of the outwards move of a damage site in each time-step: The furthest the boundary could move outwards in one timestep is K, with probability p(K) = p d = 2p(1 − p); and the probability of moving x < K steps is p(x) = p K−x s p d , where p s = p 2 +(1−p) 2 . From this we obtain the mean and variance of the boundary steps and the Central Limit Theorem (CLT) allows us to obtain the full distribution (valid in the long time and p p c limit) as a Gaussian.
This model correctly predicts various aspects of the decorrelator. First, we obtain the long-time value of D(x, t) inside the light-cone as given by D 0 = 2p(1 − p) which corresponds to the spins pointing up and down with random probability p and 1 − p, see the dashed line in the inset of Fig.1(a). Second, the model predicts the butterfly velocity, given by the following closed-form expression in p which can be read- ily expanded into a power series It agrees with the full model's butterfly velocity at large p away from p c as shown in the inset of Fig. 4. Third, the standard deviation of the boundary distribution is which confirms the Gaussian form of the boundary random walk with a variance that scales linearly with time. As detailed in the Supplementary Material, the cumulative distribution of the boundary then allows us to obtain the complete functional form of the OTOC as a function of v, governed by its inverse width µ(p). Figure.2 shows the predicted analytic form of the decorrelator (dashed lines) compared to the numerical data (solid lines). This quantitative agreement in the long-time limit confirms that the boundary controls the dynamics of chaos spreading of the full KCA model. Most notably, in the long-time limit the analytical model recovers the linear decay of the decorrelator in time, described by a VDLE for v close to v b : with an exponent β = 2. Therefore, in this regime we expect a data collapse around v b by plotting ln(D)/t against . Indeed, as shown in Fig. 4, the data of these two variables for a range of probabilities fall onto the single curve of the analytical prediction (black dashed line).
Discussion. Given the discrete nature of the dynamics, the quantitative agreement between our model and the full numerical simulation is somewhat surprising but highlights the universal features of chaos spreading. Crucially, our analysis is only valid after sample averaging and in the long time limit where the effects of the discrete perturbation and dynamics have been smoothed out. In particular, when approaching the critical point p c from above we see systematic deviations due to fluctuations.
Our work on a minimal classical model is inspired by recent developments in the study of chaos in quantum many-body systems where OTOCs have become a powerful quantitative tool. One important prediction in that context is that the butterfly velocity is bounded at high temperatures v b (T → ∞) ∼ v LR [43] where v LR is the Lieb-Robinson velocity which is the upper limit of information propagation in short-range interacting non-relativistic quantum systems [26,44]. For KCA one may connect the model parameter p to a temperature T via p = e −1/T /(e 1/T + e −1/T ) [45]. Then from Eq.6 we can obtain the full temperature dependence of the butterfly velocity, and the high temperature limit is It confirms that also in our classical many-body system the maximum velocity of information spreading v b (T → ∞) = K − 1, is always less than the maximum v max = K allowed by the local dynamics. Conclusion and Outlook. We have constructed a local diagnostic of information spreading for a one-dimensional random CA in analogy with recent semi-classical versions of OTOCs. We demonstrated that it displays ballistic propagation characterised by a butterfly velocity and exponential growth in time captured by a VDLE. We developed a random walk model of the boundary of information spreading which permits the calculation of the full functional form of the classical OTOCs, including the exponent β of the VDLE.
An obvious extension of our work is to consider the 2dimensional KCA where an even richer phenomenology is expected. Even though our boundary random walk model proposed here is simple, it has many physical and mathematical aspects that remain unexplored. For example, if the random walk is approximated to have continuous walk distance which still follows the same probability distribution, can it have a closed form solution for early times without having to invoke the CLT? Alternatively, it would be worthwhile to explore a Langevin like description of the boundary which is particularly useful while considering KCA in higher dimensions. For example, in 2D KCA, we expect a diffusion-like process for the perturbation which could potentially percolate across the lattice. In particular, it is interesting to investigate if this falls into any percolation universality class where established critical exponents could be linked to the exponent of the classical OTOCs.
Generally we expect our decorrelator and theoretical tools to be applicable to other stochastic models with discrete variables, which are widely used to describe the dynamics of both quantum and classical many-body systems [46]. We expect that different variants of the boundary theory developed here should be able to predict the scaling forms in these systems as well. In particular, it would be interesting to see them applied to models with charge or dipole conservation rules [47,48] where conjectured critical exponents could potentially be derived analytically.