# Synchronization and Maximum Lyapunov Exponents of Cellular Automata

## Abstract

We study the synchronization of totalistic one dimensional cellular automata (CA). The CA with a non zero synchronization threshold exhibit complex non periodic space time patterns and conversely. This synchronization transition is related to directed percolation. We study also the maximum Lyapunov exponent for CA, defined in analogy with continuous dynamical systems as the exponential rate of expansion of the linear map induced by the evolution rule of CA, constructed with the aid of the Boolean derivatives. The synchronization threshold is strongly correlated to the maximum Lyapunov exponent and we propose approximate relations between these quantities. The value of this threshold can be used to parametrize the space time complexity of CA.

Cellular automata (CA) are discrete dynamical systems that may exhibit complex space time patterns . It has been observed that CA may be synchronized by a stochastic coupling . We study all totalistic CA with four, five and six neighbors. For the CA considered, we find that a synchronization threshold is reached critically and that all CA with complex non periodic space time patterns have a positive threshold and conversely. We also find a strong relationship between the synchronization threshold and the maximum Lyapunov exponent (MLE) of CA .

Let us start considering the following asymmetric coupling for a continuous one-dimensional map $`f(x)`$ 

$$\begin{array}{cc}\hfill x^{}& =f(x),\hfill \\ \hfill y^{}& =(1p)f(y)+pf(x),\hfill \end{array}$$

with $`x=x(t)`$, $`x^{}=x(t+1)`$ (idem for $`y`$) and $`0p1`$. The function $`f(x)`$ depends in general by one (or more) parameter $`a`$; let assume that, for the chosen value of $`a`$, $`f(x)`$ is chaotic with Lyapunov exponent $`\lambda `$, and that $`x(0)y(0)`$. Then, $`x(t)`$ is always different from $`y(t)`$ for $`p=0`$, while for $`p=1`$ $`x`$ and $`y`$ synchronize in one time step. There exists a critical synchronization threshold $`p_c`$ for which both trajectories $`x(t)`$ and $`y(t)`$ become indistinguishable in the long time limit and

$$p_c=1\mathrm{exp}(\lambda ).$$
(1)   
In what follow we shall try to develop similar relations for CA. We begin with a brief review of the definition of maximum Lyapunov exponent for CA based on a linear expansion of the evolution rule. We then present a synchronization mechanism and show that the distance between two realizations goes to zero in a critical manner at $`p_c`$. The numerical experiments show a relation between $`p_c`$ and the maximum Lyapunov exponent which may be understood by considering several probabilistic CA. We restrict our study to one dimensional, totalistic Boolean cellular automata with four, five and six inputs, since their number is reasonably manageable and their evolution can be efficiently implemented .

A Boolean CA $`F`$ of range $`r`$ is defined as a map on the set of configurations $`\{𝐱\}`$ with $`𝐱=(x_0,\mathrm{},x_{N1})`$, $`x_i=0,1`$, and $`i=0,\mathrm{},N1`$ such that

$`𝐱^{}=F(𝐱)`$

where $`𝐱=𝐱(t)`$, $`𝐱^{}=𝐱(t+1)`$ and $`t=0,1\mathrm{}`$. The map $`F`$ is defined locally on every site $`i`$ by

$$x_i^{}=f(\{x_i\}_r)$$

where $`\{x_i\}_r=(x_i,\mathrm{},x_{i+r1})`$ is the neighborhood of range $`r`$ of site $`i`$ at time $`t`$, assuming periodic boundary conditions. For totalistic CA, the local function $`f`$ is symmetric and depends only on $`s`$ defined by

$$s(\{x_i\}_r)=\underset{j=0}{\overset{r1}{}}x_{i+j}.$$

That is $`x_i^{}=f(s(\{x_i\}_r))`$. It is useful to introduce the following operations between Boolean quantities: the sum modulo two (XOR), denoted by the symbol $``$, and the AND operation, which is analogous to the usual multiplication and shares the same symbol. These operations can be performed between two configurations component by component. We introduce the difference, or damage, $`𝐳=𝐱𝐲`$, whose evolution is given by $`𝐳^{}=F(𝐱)F(𝐲)`$ and we define the norm of $`𝐳`$ as $`|𝐳|=(1/N)_ix_iy_i`$.

A function $`f(x_i,\mathrm{},x_j,\mathrm{},x_{i+r})`$ is sensitive to its $`j`$-th argument for a given neighborhood $`(\{x_i\}_r)`$ if the Boolean derivative

$$\frac{f}{x_j}|_{\{x_i\}_r}=f(x_i,\mathrm{},x_j,\mathrm{})f(x_i,\mathrm{},x_j1,\mathrm{})$$

is 1. The Jacobian matrix $`J`$ of $`F`$ is an $`N\times N`$ matrix with components

$$J_{i,j}(𝐱)=\frac{f}{x_j}|_{\{x_i\}_r}.$$

The matrix $`J`$ is circular with zeroes everywhere except possibly on the main diagonal and the following $`r1`$ upper diagonals.

It is possible to “Taylor expand” a Boolean function around a given point using Boolean derivatives . To first order in $`|𝐳|`$ we have

$$F(𝐲)=F(𝐱)J(𝐱)𝐳$$
(2)   
where $``$ denotes the Boolean multiplication of a matrix by a vector. Compared to algebraic multiplication of a matrix by a vector, the sum and multiplication of scalars are replaced by the XOR and the AND operations respectively. Using Eq. (2) we may approximate the evolution of the damage configuration $`𝐳`$ by

$$𝐳^{}=J(𝐱)𝐳.$$

However, $`|𝐳|`$ grows at most linearly with $`t`$ since a damage cannot spread to more than $`r`$ neighbors in one time step: a fixed site $`i`$ at time $`t+1`$ can be damaged if at least one of its $`r`$ neighbors at time $`t`$ is damaged, but if more than one of the neighbors is damaged, the damage may cancel. Since

$$z_i^{}=\underset{j=i}{\overset{i+r1}{}}J_{i,j}(𝐱)z_j,$$

$`z_i^{}=1`$ if $`J_{i,j}(𝐱)z_j=1`$ on an odd number of sites. In order to account for all possible damage spreading we choose to consider each damage independently. If, at time $`t`$, $`m`$ damaged sites are present, we consider $`m`$ replicas each one with a different damaged site. On each replica, the damage evolves for one time step, without interference effects and so on.

This procedure is equivalent to choosing a vector $`𝝃(0)=𝐳(0)`$ that evolves in time according to

$$𝝃^{}=J(𝐱)𝝃$$
(3)   
where now the matrix multiplication is algebraic. The components $`\xi _i`$ are positive integers that count the number of ways in which the initial damage can spread to site $`i`$ at time $`t`$ on the ensemble of replicas.

We define the maximum Lyapunov exponent $`\lambda `$ of the cellular automaton $`F`$ by

$$\lambda (𝐱^0)=\underset{T\mathrm{}}{lim}\underset{N\mathrm{}}{lim}\frac{1}{T}\mathrm{log}\left(\frac{|𝝃^T|}{|𝝃^0|}\right),$$

where $`|𝝃|`$ may be taken as the Euclidean norm or as the sum of its components. The geometrical average $`\mu `$ of ones in the Jacobian matrix $`J`$ is defined by

$$\mu (𝐱^0)=\underset{T\mathrm{}}{lim}\underset{N\mathrm{}}{lim}\left(\underset{t=0}{\overset{T1}{}}\frac{1}{r}\underset{i,j}{}J_{i,j}(𝐱^t)\right)^{1/T}.$$

We show in Fig. 1 the points $`(r\mu ,\lambda )`$ of all the totalistic CA with $`r=4,5`$, and 6 that show complex non periodic space time patterns.

The process defined by Eq. (3) may be viewed as a deterministic directed bond percolation problem where a site $`i`$ at depth $`t`$ is wet if $`\xi _i(t)>0`$. The bonds exist where the components of $`J`$ are 1. A first approximation is obtained by replacing $`J`$ with a random matrix whose elements are zero except on the diagonal and the $`r1`$ following upper diagonals, where they are equal to one with probability $`\mu `$ . There is a critical value $`\mu _c(r)`$ below which the bond percolation process falls into the absorbing state so that the maximum eigenvalue of the product of random matrices is zero. We can further introduce a mean field approximation to the directed bond percolation process, which exhibits discrepancies only very near to $`\mu _c`$. In this case one can show that

$$\stackrel{~}{\lambda }=\mathrm{ln}(r\mu )$$
(4)   
is an upper bound to the MLE of the product of random matrices. The behavior of $`\stackrel{~}{\lambda }`$ is plotted in Fig. 1, and we note that it is a good approximation to the generic behavior of CA. We also report the value $`r\mu _c(r)`$ for which the maximum eigenvalue of the product of random matrices is zero, corresponding to the percolation threshold for the directed bond percolation problem. We found that $`r\mu _c(r)1.3`$ regardless of $`r`$, a fact that can be understood from the following argument. The percolation cluster has on the average $`\kappa `$ bonds per site with $`\kappa (r)=r\mu _c(r)`$; for this percolation model, one connection to a wet neighbor at the previous time step is sufficient, and further connections do not alter wetting.

We now discuss a synchronization mechanism for CA. Starting with two initial configurations chosen at random $`𝐱(0)`$ and $`𝐲(0)`$ we propose that

$$\begin{array}{cc}\hfill 𝐱^{}& =F(𝐱),\hfill \\ \hfill 𝐲^{}& =\overline{S^t(p)}F(𝐲)S^t(p)F(𝐱),\hfill \end{array}$$

where $`S^t(p)`$ is a Boolean random diagonal matrix with elements $`s_i^t(p)`$ that, at each time step, take the value one with probability $`p`$ and zero with probability $`1p`$; $`\overline{S(p)}=IS(p)`$ and $`I`$ is the identity matrix. On the average, $`y_i^{}`$ will be set to the value of $`x_i^{}=f(\{x_i\})`$ on a fraction $`p`$ of sites. In this way we introduced a stochastic synchronization mechanism over a deterministic process. This stochatic mechanism can be considered as a “random field” approximation of an intermittent coupling generated by a deterministic chaotic process. The evolution equation for the difference $`𝐳=𝐱𝐲`$ is

$$𝐳^{}=\overline{S(p)}\left[F(𝐱)F(𝐲)\right].$$
(5)   
The control and order parameters are $`p`$ and $`h(p)=lim_t\mathrm{}lim_N\mathrm{}|𝐳(t)|`$ respectively. We say that $`𝐱`$, the driver, and $`𝐲`$, the driven system, synchronize when $`h(p)=0`$. For $`p=0`$ both systems evolve independently, while for $`p=1`$ they synchronize in just one step; we expect then to find a synchronization threshold $`p_c`$. This behavior is shared by all the CA with complex non periodic space time patterns. All others synchronize for $`p0`$. This can be conversely expressed by saying that all CA that synchronize with a nontrivial $`p_c`$ exhibit complex non periodic space time patterns.

For totalistic linear rules, whose evolution is given by

$$f(\{x_i\}_r)=\underset{j=0}{\overset{r1}{}}x_{i+j},$$

the synchronization equation (5) is equivalent to the dilution (with probability $`1p`$) of the rule. For $`r=2`$ the dilution problem is equivalent to the line $`z=0`$ of Ref. , whose transition belongs to the universality class of directed percolation (DP) . For a similar synchronization mechanism, it has been recently claimed that the elementary CA rule 18 does not  and does  belong to the DP universality class. The presence of a single absorbing state and the absence of other conserved quantities (i.e. number of kinks) strongly suggests that the synchronization transition belongs to the DP universality class . The goal of this work was not that of computing critical exponents. However, since the critical point was located by means of the scaling law for the density of defects (in order to minimize finite size and time effects), we had the opportunity of computing the magnetic exponent $`\beta `$ for all the CA studied. Due to the large number of CA examined, this computation was performed semi-automatically, and the precision of the resulting exponent is quite low, nontheless we have not found any example of non-DP behavior. We performed bulk simulations for $`N=2000`$ and $`T=4000`$ and checked that the results do not change for doubling of halving these figures.

In Fig. 2 we plot the points $`(p_c,\lambda )`$ for all totalistic CA with $`r=4,5,6`$ and a non trivial value of $`p_c`$. All the numerical experiments were performed using a parallel algorithm that takes care of all the values of $`p`$ simultaneously .

Let us study the relation between $`p_c`$ and $`\lambda `$ by some random approximations. Near the synchronization threshold $`p_c`$, we may expand $`𝐲`$ around $`𝐱`$ with the help of Eq. (2) so that

$$𝐳^{}=\overline{S(p)}J(𝐱)𝐳.$$
(6)   
Eq. (6) may be written as

$$z_i^{}=\overline{s_i(p)}\left[J_{i,i}(𝐱)z_i\mathrm{}J_{i,i+r1}(𝐱)z_{i+r1}\right].$$
(7)   
During the time evolution of a particular CA a fixed value of $`\mu `$ is attained so that on average $`r\mu `$ derivatives inside the parenthesis on the r.h.s. of Eq. (7) are different from zero. A first approximation, model A, is obtained by replacing the derivatives with random variables $`m_i(\mu )`$ that are one with probability $`\mu `$ and zero with probability ($`1\mu `$). That is

$$z_i^{}=\overline{s_i(p)}\left[m_i(\mu )z_i\mathrm{}m_{i+r1}(\mu )z_{i+r1}\right].$$

On every site and every time step, $`r\mu `$ variables are chosen on average and if their sum is odd, then with probability $`p`$, $`z_i^{}=1`$. We then look for the synchronization threshold $`p_c(\mu )`$ and plot it as a function of $`\stackrel{~}{\lambda }=\mathrm{ln}(r\mu )`$, Eq. (4). This is the curve labeled $`A`$ in Fig. 2. The predicted values of $`p_c`$ are generally higher than those found for CA for the same value of $`\lambda `$ possibly due to correlations among the derivatives.

Since the typical “complex” CA pattern exhibit transient correlations (“triangles”), one can model them by choosing a fixed number $`kr`$ of derivatives equal to one. The simplest way is to take $`k=r\mu `$ with $`r\mu `$ an integer, model B, which . Then

$$z_i^{}=\overline{s_i(p)}\left[z_i\mathrm{}z_{i+k1}\right].$$

which is a dilution of the XOR with $`k`$ inputs. This process is expected to belong to the same universality class of directed percolation . The curve labeled B of Fig. 2 passes through all the calculated values for $`k=2,\mathrm{}6`$. We note that this second model is a better fit of CA behavior.

We can extend this last model allowing noninteger values of $`r\mu `$ by

$$z_i(t+1)=\overline{s_i^t(p)}m_i^t(\mu )\left[z_i\mathrm{}z_{i+k1}\right]$$

which we call model C<sub>k</sub>. In this way the average number of non zero derivatives is $`k\mu `$ with $`0\mu 1`$. Now $`k`$ is a free parameter, and this model can be useful to delimit the expected spread of $`(\lambda ,p_c)`$ points.

Since $`s_i^t(p)`$ and $`m_i^t(\mu )`$ are independent random variables we may write

$$z_i(t+1)=\overline{s_i^t(q)}\left[z_i\mathrm{}z_{i+k1}\right].$$

where $`(1q)=(1p)\mu `$. In this guise, this is model B with $`k`$ inputs. The synchronization threshold is given by $`p_c(k,\mu )=1(1q_c(k))/\mu `$ where $`q_c(k)`$ is the percolation threshold of the dilution of the XOR with $`k`$ inputs. Using the approximation $`\lambda =\mathrm{ln}(k\mu )`$ one has

$$p_c(k,\lambda )=1k\left(1q_c(k)\right)\mathrm{exp}(\lambda ).$$
(8)   
which bears some likeness to Eq. (1). The curves labeled $`C_2`$ and $`C_6`$ of Fig. 2 correspond to this last expression for $`k=2`$ and $`k=6`$ respectively. We note that the points $`(\lambda ,p_c)`$ of almost all the CA considered fall between these two curves.

For the totalistic one dimensional CA with $`r=4,5,6`$ we can safely say that all CA with a positive $`p_c`$ exhibit complex non periodic space time patterns and conversely. These CA also have a positive MLE . We also showed that the synchronization of CA is a critical phenomenon similar to directed percolation.

We proposed several approximations based on a combination of “linearization” of CA rules using Boolean Taylor expansions and stochasticity and showed the relation between the synchronization threshold and the MLE. In particular, model $`C`$ implies a relation similar to that found for continuous maps with the addition of a percolation threshold constant.

An analogous mechanism can be applied to coupled map lattices; in this case $`p`$ is the probability that a map $`y_i(t)`$ takes the value $`x_i(t)`$. One observes a synchronization transitions, but the critical value $`p_c`$ is not correlated to the usual MLE .

We acknowledge partial financial support from project DGAPA-UNAM IN103595.