# Evaluating the Quality of Local Structure Approximation Using Elementary Rule 14

Henryk Fukś, Francis Kwaku Combert

► To cite this version:

Henryk Fukś, Francis Kwaku Combert. Evaluating the Quality of Local Structure Approximation Using Elementary Rule 14. 24th International Workshop on Cellular Automata and Discrete Complex Systems (AUTOMATA), Jun 2018, Ghent, Belgium. pp.43-56, 10.1007/978-3-319-92675-9_4. hal-01824874

HAL Id: hal-01824874

https://inria.hal.science/hal-01824874

Submitted on 27 Jun 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Evaluating the Quality of Local Structure Approximation Using Elementary Rule 14

Henryk Fukś and Francis Kwaku Combert

Department of Mathematics and Statistics, Brock University,
St. Catharines, ON, Canada
hfuks@brocku.ca, fc15uy@brocku.ca

**Abstract.** Cellular automata (CA) can be viewed as maps in the space of probability measures. Such maps are normally infinitely-dimensional, and in order to facilitate investigations of their properties, especially in the context of applications, finite-dimensional approximations have been proposed. The most commonly used one is known as the local structure theory, developed by H. Gutowitz et al. in 1987. In spite of the popularity of this approximation in CA research, examples of rigorous evaluations of its accuracy are lacking. In an attempt to fill this gap, we construct a local structure approximation for rule 14, and study its dynamics in a rigorous fashion, without relying on numerical experiments. We then compare the outcome with known exact results.

**Keywords:** rule 14, local structure approximation, invariant manifolds

## 1 Introduction

One-dimensional elementary cellular automata (CA) can be viewed as maps in the space of probability measures over bi-infinite binary sequences (to be called *configurations*). This can be understood as follows. Suppose that we start with a large set of initial configurations drawn from a certain distribution (for example, from the Bernoulli distribution). Let us now suppose that we apply a given cellular automaton rule to all these configurations. The resulting set of configurations is usually no longer described by Bernoulli distribution, but by some other distribution. We can thus say that the CA rule transforms the initial probability measure into some other measure, and when we apply the local rule again and again, we obtain a sequence of measures, to be called the *orbit* of the initial measure.

This approach, however, is not without difficulties. In order to fully describe a probability measure over bi-infinite binary sequences, one needs to specify infinitely many *block probabilities*, that is, probabilities of the occurrence of 0, 1, 00, 01, 10, 11, 000, etc – in short, the probabilities of occurrence of all possible binary words. This means that the CA rule treated as a map in the space of probability measures is an *infinitely-dimensional map*.

Infinite-dimensional maps are difficult to investigate, even numerically, thus from the early days of CA research, efforts were made to find a way to approximate them by finite-dimensional maps. In a seminal paper [8], published over 30
years ago, H. Gutowitz et al. proposed such an approximation, which they called the *local structure theory*. It was an application of a well known idea of Bayesian extension, widely used in statistical physics as a basis of so-called mean-field theories, finite-cluster approximations, and related methods.

Since 1987 the local structure theory has been widely used in CA research, as witnessed by a large number of citations of [8]. This could be somewhat surprising, given that relatively few rigorous results are known about the local structure theory. Usually, the authors using this method simply construct a finite-dimensional map or recurrence equations following the recipe given in [8], and declare that these possess orbits approximating the dynamics of the actual CA or related system which they investigate. Judgments on the quality of the approximation are usually made based on numerical iterations of local structure maps and numerical simulations of the CA in question. Numerical results are thus compared with other numerical results.

In recent years, however, partial orbits of Bernoulli measures have been computed for some selected elementary CA [7], making a somewhat more rigorous approach possible. The goal of this paper is to provide an example of a CA rule for which some block probabilities are known exactly, and for which local structure equations can be analyzed rigorously, without relying exclusively on numerical iterations. This way, the quality of the approximation could be evaluated in a solid and rigorous fashion, without worrying about numerical errors, finite size effects, etc.

We selected elementary CA rule 14 as the most promising example for such study. It has several interesting features: exact probabilities of blocks of length up to three are known for the orbit of the symmetric Bernoulli measure under this rule, and some of these block probabilities exhibit non-trivial behaviour - for example, convergence toward the steady state as a power law with fractional exponent. At the same time, rule 14 conserves the number of pairs 10 [4], and the existence of this additive invariant provides a constrain simplifying local structure equations, making them easier to analyze. Since block probabilities of length 3 are known for this rule, we will construct local approximation of level 3 and investigate its dynamics not only by simple numerical iterations, but by finding invariant manifolds at the fixed point and determining the nature of the flow on these manifolds.

One should stress here that in what follows we will use only very minimal formalism. More formal details about the construction of probability measures over infinite bisequences and the construction of local structure maps for arbitrary rules (both deterministic and probabilistic) can be found in [5], where the reader will also find more references on these subjects.

## Preliminary remarks about rule 14

Consider the fully discrete dynamical system (called *cellular automaton*) where $s_i(n) \in \{0,1\}$ is the state of site $i \in \mathbb{Z}$ at time $n \in \mathbb{N}$, with dynamics defined by $s_i(n+1) = f(s_{i-1}(n), s_i(n), s_{i+1}(n))$. The function $f: \{0,1\}^3 \to \{0,1\}$
is called the *local rule*. In this paper, we will consider $f$ which is defined by
$f(x_0, x_1, x_2) = x_1 + x_2 - x_1x_2 - x_0x_1 - x_0x_2 + x_0x_1x_2$, and we call the above
*rule 14*, following the numbering scheme of Wolfram [11].

Usually, the initial state at $n = 0$ is drawn from the Bernoulli distribution, where each site $s_i(0)$ is either in state 1 with probability $\rho$, or in state 0 with probability $1 - \rho$, independently of each other, where $\rho \in [0, 1]$. When $\rho = 1/2$, we call this symmetric Bernoulli distribution.

A classical problem in cellular automata theory is to compute the probability of the occurrence of a given binary string **a** in a configuration obtained after *n* iterations of the rule, assuming that the initial configuration is drawn from the Bernoulli distribution. Such probability will be denoted by $P_n(\mathbf{a})$ and called *block probability*. It is easy to show that if the initial distribution is Bernoulli, then the probability of occurrence of **a** is independent of its position in the configuration. We will call such block probabilities *shift invariant*.

The set of shift-invariant block probabilities $P_n(\mathbf{a})$ for all binary strings $\mathbf{a}$ defines a shift-invariant probability measure on the set of infinite binary sequences, but we will not be concerned with the formal construction of such measures here. Interested reader can find all relevant details and references in [5].

Consider now a configuration in which $s_i(n + 1) = 1$. By using the definition of rule $f$, one can easily figure out that $s_i(n + 1)$ is determined entirely by the triple $(s_{i-1}(n), s_i(n), s_{i+1}(n))$, and that the only possible values of $(s_{i-1}(n), s_i(n), s_{i+1}(n))$ producing $s_i(n+1) = 1$ are $(0, 0, 1)$, $(0, 1, 0)$ or $(0, 1, 1)$. This means that probability of obtaining 1 at time $n+1$ is equal to the sum of probabilities of occurrence of blocks 001, 010, and 011 at time $n$, $P_{n+1}(1) = P_n(001) + P_n(010) + P_n(011)$. One can carry out a similar reasoning for longer blocks. For example, a pair of 1s, that is, $s_i(n+1) = 1$ and $s_{i+1}(n+1) = 1$, can appear only and only if at the previous time step $n$ the lattice positions $i-1, i, i+1, i+2$ assumed values 0,0,1,0 or 0,0,1,1, i.e., $(s_{i-1}(n), s_i(n), s_{i+1}(n), s_{i+2}(n)) = (0, 0, 1, 0)$ or $(s_{i-1}(n), s_i(n), s_{i+1}(n), s_{i+2}(n)) = (0, 0, 1, 1)$. This yields $P_{n+1}(11) = P_n(0010) + P_n(0011)$.

Obviously, one can write analogous equations for probabilities of any binary
block, obtaining an infinite system of difference equations. The complete set of
such equations for blocks of length up to 3 for rule 14 is shown below.

$$
\begin{align*}
P_{n+1}(0) &= P_n(000) + P_n(100) + P_n(101) + P_n(110) + P_n(111), \\
P_{n+1}(1) &= P_n(001) + P_n(010) + P_n(011), \\
P_{n+1}(11) &= P_n(0010) + P_n(0011), \\
P_{n+1}(00) &= P_n(0000) + P_n(1000) + P_n(1100) + P_n(1101) + P_n(1110) + P_n(1111), \\
P_{n+1}(01) &= P_n(0001) + P_n(1001) + P_n(1010) + P_n(1011), \\
P_{n+1}(10) &= P_n(0100) + P_n(0101) + P_n(0110) + P_n(0111), \\
P_{n+1}(000) &= P_n(00000) + P_n(10000) + P_n(11000) + P_n(11100) + P_n(11110) \\
&\qquad + P_n(11110) + P_n(11111), \\
P_{n+1}(001) &= P_n(00001) + P_n(10001) + P_n(11001) + P_n(11010) + P_n(11011), \\
P_{n+1}(010) &= P_n(10100) + P_n(10101) + P_n(10110) + P_n(10111),
\end{align*}
$$
$$
\begin{align*}
P_{n+1}(011) &= P_n(00010) + P_n(00011) + P_n(10010) + P_n(10011), \\
P_{n+1}(100) &= P_n(01000) + P_n(01100) + P_n(01101) + P_n(01110) + P_n(01111), \\
P_{n+1}(101) &= P_n(01001) + P_n(01010) + P_n(01011), \\
P_{n+1}(110) &= P_n(00100) + P_n(00101) + P_n(00110) + P_n(00111), \\
P_{n+1}(111) &= 0. \tag{1}
\end{align*}
$$

One thing which is immediately obvious is that not all of these equations are in-
dependent because the block probabilities themselves are not independent. Block
probabilities must satisfy so-called *Kolmogorov consistency conditions*, which are
in fact just additivity conditions satisfied by a measure induced by block proba-
bilities. For example, we must have $P_n(1)+P_n(0) = 1$, $P_n(01)+P_n(00) = P_n(0)$,
etc. Consistency conditions can be used to express some block probabilities by
others. One can show that for binary strings, among probabilities of blocks of
length $k$, only $2^{k-1}$ are independent [5], in the sense that one can choose $2^{k-1}$
block probabilities which are not linked to each other via consistency conditions.
For blocks of length up to 3, there are 14 block probabilities, $P_n(0)$, $P_n(1)$,
$P_n(00)$, $P_n(01)$, $P_n(10)$, $P_n(11)$, $P_n(000)$, $P_n(001)$, $P_n(010)$, $P_n(011)$, $P_n(100)$,
$P_n(101)$, $P_n(110)$, and $P_n(111)$. Among them only $2^{3-1} = 4$ are independent.
While there is some freedom in choosing which ones are to be treated as inde-
pendent, we will choose the following four, $P_n(0)$, $P_n(00)$, $P_n(000)$, and $P_n(010)$.
This is called the *short block representation*, and a detailed algorithm for choos-
ing block this way is described in [5]. Here it is sufficient to say that short block
representation ensures that the blocks selected as independent are the shortest
possible ones.

Using consistency conditions, one can now express the remaining blocks of
length up to 3 in terms of $P_n(0)$, $P_n(00)$, $P_n(000)$, and $P_n(010)$, as follows:

$$
\begin{align*}
P_n(1) &= 1 - P_n(0), \\
P_n(01) &= P_n(0) - P_n(00), \\
P_n(10) &= P_n(0) - P_n(00), \\
P_n(11) &= 1 - 2P_n(0) + P_n(00), \\
P_n(001) &= P_n(00) - P_n(000), \\
P_n(011) &= P_n(0) - P_n(00) - P_n(010), \\
P_n(100) &= P_n(00) - P_n(000), \\
P_n(101) &= P_n(0) - 2P_n(00) + P_n(000), \\
P_n(110) &= P_n(0) - P_n(00) - P_n(010), \\
P_n(111) &= 1 - 3P_n(0) + 2P_n(00) + P_n(010). \tag{2}
\end{align*}
$$

Using the above substitutions one can reduce eqs. (1) to the following set of four
equations,

$$
\begin{align}
P_{n+1}(0) &= 1 - P_n(0) + P_n(000), \tag{3} \\
P_{n+1}(00) &= 1 - 2P_n(0) + P_n(00) + P_n(000), \nonumber
\end{align}
$$
$$ 
\begin{aligned}
P_{n+1}(000) &= 1 - 3P_n(0) + 2P_n(00) + P_n(000) + P_n(010) - P_n(01000), \\
P_{n+1}(010) &= P_n(0) - 2P_n(00) + P_n(000).
\end{aligned}
$$

Note that the above cannot be iterated, because on the right hand side, in addition to the four aforementioned independent probabilities, we have probability $P_n(01000)$, the probability of the block of length 5.

Fortunately, in spite of the above problem, if the initial Bernoulli measure is symmetric, exact expressions for probabilities $P_n(0)$, $P_n(00)$, $P_n(000)$ and $P_n(010)$ for rule 14 (that is, the solution of eqs. (3)) can be obtained by combinatorial methods. We will quote the relevant results below, omitting the proof, which can be found in [6].

**Proposition 1 (Fukś et al. 2009).** For elementary rule 14, if the initial configuration is drawn from symmetric Bernoulli distribution, the probabilities of block of length up to 3 are given by

$$ P_n(0) = \frac{1}{2} \left( 1 + \frac{2n-1}{4^n} C_{n-1} \right), \qquad (4) $$

$$ P_n(00) = 2^{-2-2n}(n+1)C_n + \frac{1}{4}, \qquad (5) $$

$$ P_n(000) = 2^{-2n-3}(4n+3)C_n, \qquad (6) $$

$$ P_n(010) = 2^{-2-2n}(n+1)C_n, \qquad (7) $$

where $C_n$ is the n-th Catalan number, $C_n = \frac{1}{n+1} \binom{2n}{n}$.

Note that although the above proposition provides probabilities of $P_n(0)$, $P_n(00)$, $P_n(000)$ and $P_n(010)$ only, the remaining probabilities of blocks of length up to 3 can be easily computed using eqs. (2).

Although we know exact solution of eqs. (3), we can also attempt to obtain an approximate solution by approximating the “problematic” block probability $P_n(01000)$. There exists a method for approximating longer block probabilities by probabilities of shorter blocks. This method is called the Bayesian extension, and it is known to produce block probabilities satisfying consistency conditions [5]. Applying the Bayesian extension to $P_n(01000)$, one obtains

$$ P_n(01000) \approx \frac{P_n(010)P_n(100)P_n(000)}{P_n(10)P_n(00)}. \qquad (8) $$

In the above, by definition, the fraction on the right hand side is considered to be zero whenever its denominator is equal to zero. Using eqs. (2) we can now express $P_n(01000)$ in terms of our four independent block probabilities,

$$ P_n(01000) \approx \frac{P_n(010)(P_n(00) - P_n(000))P_n(000)}{(P_n(0) - P_n(00))P_n(00)}. \qquad (9) $$

If we replace $P_n(01000)$ in eqs. (3) by the above approximation, we will obtain the system of four coupled difference equations,

$$ x_{n+1} = -x_n + z_n + 1, \qquad (10) $$
$$y_{n+1} = -2x_n + y_n + z_n + 1, \qquad (11)$$

$$z_{n+1} = 1 + z_n + v_n - 3x_n + 2y_n - \frac{v_n(y_n - z_n)z_n}{y_n(x_n - y_n)}, \qquad (12)$$

$$v_{n+1} = x_n - 2y_n + z_n, \qquad (13)$$

where for brevity we introduced variables $x_n = P_n(0)$, $y_n = P_n(00)$, $z_n = P_n(000)$, and $v_n = P_n(010)$. Equations (10)–(13) will be referred to as local structure equations of level 3, following nomenclature of [8, 5]. The designation “level 3” pertains to the fact that we used block probabilities of length up to 3.

## 2 Exact solutions vs. local structure approximation

How does the orbit of local structure equations (10)–(13) compare with known exact solutions given by eq. (4)–(7)? In order to find this out, we will assume that the initial probability measure is symmetric Bernoulli, meaning that $x_0 = P_0(0) = 1/2$, $y_0 = P_0(00) = 1/4$, $z_0 = P_0(000) = 1/8$, and $v_0 = P_0(010) = 1/8$. Figure 1 shows differences between exact probabilities obtained from eq. (4)–

<figure><img src="image_2.png" /><figcaption>Fig. 1. Differences between exact and approximate values of block probabilities as a function of $n$. Two differences are shown, $P_n(00) - y_n$ (lower curve) and $P_n(000) - z_n$ (upper curve).</figcaption></figure>

(7) and values obtained by iterating local structure equations (10)–(13). Two
differences are shown, $P_n(00) - y_n$ (lower curve) and $P_n(000) - z_n$ (upper curve). In both cases we can see that the difference tends to zero as $n \to \infty$. Values of $P_n(0) - x_n$ and $P_n(010) - v_n$ (not shown) exhibit similar behaviour.

This indicates that even though the local structure approximation of level 3 does not produce exact values of block probabilities at finite $n$, it seems to become exact in the limit $n \to \infty$. To verify this, let us first note that from eq. (4)-(7) we obtain

$$ \lim_{n \to \infty} P_n(0) = \frac{1}{2}, \quad \lim_{n \to \infty} P_n(00) = \frac{1}{4}, \quad \lim_{n \to \infty} P_n(000) = 0, \quad \lim_{n \to \infty} P_n(010) = 0. $$

We will denote these values by $(x^*, y^*, z^*, v^*) = (\frac{1}{2}, \frac{1}{4}, 0, 0)$. One can easily verify that $(x^*, y^*, z^*, v^*)$ is a fixed point of eqs. (10)-(13). In what follows, we will investigate stability of this fixed point. We will prove that the following property holds.

**Proposition 2.** If the dynamical system given by eqs. (10)-(13) is iterated starting from initial conditions $x_0 = 1/2$, $y_0 = 1/4$, $z_0 = 1/8$, and $v_0 = 1/8$, then $\lim_{n \to \infty} (x_n, y_n, z_n, v_n) = (x^*, y^*, z^*, v^*) = (\frac{1}{2}, \frac{1}{4}, 0, 0)$.

This means that the local structure map approximates the exact probabilities remarkably well, converging to the same fixed point as the exact values. We will prove Proposition 2 by reducing local structure equations to two dimensions and by computing local manifolds at the fixed point.

## Reduction to two dimensions

Close examination of equations (10)-(13) reveals some obvious symmetries. First of all, it is easy to check that $x_{n+1} - y_{n+1} = x_n - y_n$. Since $x_0 - y_0 = \frac{1}{4}$, we have $x_n - y_n = \frac{1}{4}$ for all $n$, thus

$$ x_n = y_n + \frac{1}{4}. \qquad (14) $$

Further simplification is possible. Note that $v_{n+1} - y_{n+1} = 3(x_n - y_n) - 1 = 3 \cdot \frac{1}{4} - 1 = -\frac{1}{4}$. This implies that for any $n > 0$, $v_{n+1} = y_{n+1} - \frac{1}{4}$, or, equivalently, that for any $n > 1$,

$$ v_n = y_n - \frac{1}{4}. \qquad (15) $$

Note that this does not hold for $n = 0$, because in this case $v_0 = y_0 - 1/8$. Now, using eqs. (14) and (15), we can reduce our dynamical system to two dimension, as eqs. (11) and (12) become

$$ \begin{aligned} y_{n+1} &= -2(y_n + \frac{1}{4}) + y_n + z_n + 1, \\ z_{n+1} &= 1 + z_n + (y_n - \frac{1}{4}) - 3(y_n + \frac{1}{4}) + 2y_n - \frac{(y_n - \frac{1}{4})(y_n - z_n)z_n}{y_n((y_n + \frac{1}{4}) - y_n)}. \end{aligned} $$
After simplification we obtain, for $n \ge 1$,

$$y_{n+1} = \frac{1}{2} - y_n + z_n, \qquad (16)$$

$$z_{n+1} = z_n - \frac{(4y_n - 1)(y_n - z_n)z_n}{y_n}, \qquad (17)$$

where we start the recursion at $n = 1$, taking $y_1 = 3/8$, $z_1 = 7/32$. The last two values were obtained by direct computation of $y_1$ and $z_1$ from eqs. (11) and (12) for $n = 0$, by substituting $x_0 = 1/2$, $y_0 = \frac{1}{4}$, $z_0 = v_0 = 1/8$ on the right hand side.

We will prove the following result.

**Proposition 3.** If the dynamical system described by eqs. (16) and (17) is iterated starting at $y_1 = 3/8$, $z_1 = 7/32$, then

$$\lim_{n \to \infty} (y_n, z_n) = \left(\frac{1}{4}, 0\right). \qquad (18)$$

In order to prove the above proposition let us first denote $\mathbf{x} = \begin{bmatrix} y \\ z \end{bmatrix}$. In this notation, eqs. (16) and (17) define two-dimensional map

$$\mathbf{F}(\mathbf{x}) = \begin{bmatrix} \frac{1}{2} - y + z \\ z - \frac{(4y-1)(y-z)z}{y} \end{bmatrix}. \qquad (19)$$

It is easy to check that the map $\mathbf{F}$ has the fixed point $\mathbf{x}^* = \begin{bmatrix} \frac{1}{4} \\ 0 \end{bmatrix}$. In order to prove Proposition 3, all we need is to show that $\mathbf{x}^*$ is asymptotically stable (or at least semi-stable in the relevant domain).

The Jacobian matrix of $\mathbf{F}$ evaluated at the fixed point $\mathbf{x}^*$ is given by

$$A = \begin{bmatrix} -1 & 1 \\ 0 & 1 \end{bmatrix},$$

and its eigenvalues are $-1$ and $1$. Since these eigenvalues have an absolute value equal to $1$, the fixed point $\mathbf{x}^*$ is a non-hyperbolic fixed point and one cannot determine its stability by eigenvalues alone. We will investigate its stability by resorting to the center manifold theory.

Let $P$ be the matrix of column eigenvectors of $A$, and let $P^{-1}$ be its inverse,

$$P = \begin{bmatrix} 1 & \frac{1}{2} \\ 0 & 1 \end{bmatrix}, \quad P^{-1} = \begin{bmatrix} 1 & -\frac{1}{2} \\ 0 & 1 \end{bmatrix}.$$

We will first move the fixed point to the origin and simultaneously diagonalize the linear part of $\mathbf{F}$. The following change of variables accomplishes this task,

$$\mathbf{X} = P^{-1}(\mathbf{x} - \mathbf{x}^*), \qquad (20)$$
where the components of the new variable **X** will be denoted by **Y** and **Z**. Eq. (20) thus yields

$$Y = y - \frac{1}{4} - \frac{z}{2}, \qquad (21)$$

$$Z = z. \qquad (22)$$

Change of variables from **x** to **X** transforms the dynamical system $\mathbf{x}_{n+1} = \mathbf{F}(\mathbf{x}_n)$ into the system

$$\mathbf{X}_{n+1} = P^{-1}\mathbf{F}(P\mathbf{X}_n + \mathbf{x}^*) - P^{-1}\mathbf{x}^*. \qquad (23)$$

This yields, after simplification,

$$Y_{n+1} = -Y_n + \frac{1}{2} \frac{(4Y_n + 2Z_n)(Y_n - \frac{1}{2}Z_n + \frac{1}{4})Z_n}{Y_n + \frac{1}{2}Z_n + \frac{1}{4}}, \qquad (24)$$

$$Z_{n+1} = Z_n - \frac{(4Y_n + 2Z_n)(Y_n - \frac{1}{2}Z_n + \frac{1}{4})Z_n}{Y_n + \frac{1}{2}Z_n + \frac{1}{4}}. \qquad (25)$$

One can immediately see that the above system has (0,0) as a fixed point, and that its linear part is given by $Y_{n+1} = -Y_n$, $Z_{n+1} = Z_n$. As mentioned earlier, there is nothing we can say about the stability of (0,0) by examining the linear part alone, except that in the vicinity of (0,0) the **Y** variable is changing its sign at each iteration. We will use the method outlined in [10] to find the invariant manifold corresponding to -1 eigenvalue. We will call this manifold the *flip manifold* and denote it by $W^f$.

Let us assume that $W^f$ has the equation $Z = h(Y)$, where $h$ in the vicinity of 0 is given by the series $h(Y) = a_2Y^2 + a_3Y^3 + a_4Y^4 + a_5Y^5 + ...$. Note that the series starts from the quadratic term, and this is because the manifold $Z = h(Y)$ must be tangent to the **Y** axis (we already diagonalized our dynamical system).

The condition for invariance of $W^f$ requires that the relationship $Z_n = h(Y_n)$ remains valid in the next time step, meaning that $Z_{n+1} = h(Y_{n+1})$. Let us rewrite eqs. (24) and (25) as

$$Y_{n+1} = G_1(Y_n, Z_n), \qquad (26)$$

$$Z_{n+1} = G_2(Y_n, Z_n), \qquad (27)$$

where

$$G_1(Y, Z) = -Y + \frac{1}{2} \frac{(4Y + 2Z)(Y - \frac{1}{2}Z + \frac{1}{4})Z}{Y + \frac{1}{2}Z + \frac{1}{4}}, \qquad (28)$$

$$G_2(Y, Z) = Z - \frac{(4Y + 2Z)(Y - \frac{1}{2}Z + \frac{1}{4})Z}{Y + \frac{1}{2}Z + \frac{1}{4}}. \qquad (29)$$

Condition $Z_{n+1} = h(Y_{n+1})$ now becomes $G_2(Y, Z) = h(G_1(Y, Z))$, and, by taking $Z = h(Y)$, it yields

$$G_2(Y, h(Y)) = h(G_1(Y, h(Y))). \qquad (30)$$
This means that if we expand $G_2(Y, h(Y)) - h(G_1(Y, h(Y)))$ into the Taylor series with respect to $Y$, all coefficient of the expansion should be zero. Such expansion, done by the Maple symbolic algebra system, yields

$$ 
\begin{aligned}
G_2(Y, h(Y)) - h(G_1(Y, h(Y))) = & (2a_3 - 4a_2)Y^3 + \\ & (-4a_3 - 4\left(4 + \frac{1}{2}a_2\right)a_2 + 16a_2 + 4a_2^2)Y^4 + \\ & \left(2a_5 - 4a_4 - 4\left(4 + \frac{1}{2}a_2\right)a_3 - 8a_2a_3 + 16a_2^2 + 16a_3\right) \\ & - \left(-4a_3 - 4\left(4 + \frac{1}{2}a_2\right)a_2 + 16a_2\right)a_2)Y^5 + O(Y^6).
\end{aligned}
$$

Coefficients in front of $Y^3, Y^4, Y^5, \dots$ must be zero, yielding the system of equations for $a_2, a_3, a_4, \dots$,

$$ 
\begin{aligned}
0 &= 2a_3 - 4a_2, \\
0 &= -4a_3 - 4\left(4 + \frac{1}{2}a_2\right)a_2 + 16a_2 + 4a_2^2, \\
0 &= 2a_5 - 4a_4 - 4\left(4 + \frac{1}{2}a_2\right)a_3 - 8a_2a_3 + 16a_2^2 + 16a_3 \\
  &\quad - \left(-4a_3 - 4\left(4 + \frac{1}{2}a_2\right)a_2 + 16a_2\right)a_2, \\
  &\dots
\end{aligned}
\quad (31) $$

Solving the above system one obtains $a_2 = 4, a_3 = 8, a_4 = 32, a_5 = -32$, etc. The flip manifold $W^f$ is, therefore, given by

$$ Z = h(Y) = 4Y^2 + 8Y^3 + 32Y^4 - 32Y^5 + O(Y^6). \quad (32) $$

By substituting $Z_n$ by $h(Y_n)$ on the right hand side of eq. (24) and Taylor expanding again one obtains the equation describing the dynamics on the flip manifold $W^f$,

$$ Y_{n+1} = -Y_n + 8Y_n^3 + 32Y_n^4 + O(Y_n^5). \quad (33) $$

The above equation has 0 as a fixed point, and we need to determine its stability. Recall that a fixed point $\bar{x}$ of $x_{n+1} = f(x)$ is said to be asymptotically stable if there exist $\delta > 0$ such that for any $x_0$ satisfying $|x_0 - \bar{x}| < \delta$ we have $\lim_{n \to \infty} x_n = \bar{x}$. We will use the following general test for asymptotic stability [9].

**Theorem 1 (Murakami 2005).** Let $\bar{x}$ be a fixed point of $x_{n+1} = f(x)$. Suppose that $f \in C^{2k-1}(\mathbb{R})$, $f'(\bar{x}) = -1$, $f''(\bar{x}) = 0$ for $j \in \{2, 3, \dots, k-1\}$, and that $f^{(k)}(\bar{x}) \neq 0$. If $k$ is odd and $f^{(k)}(\bar{x}) > 0$, then $\bar{x}$ is asymptotically stable.

In our case, for eq. (33), $f(x) = -x + 8x^3 + 32x^4 + O(x^5)$, $\bar{x} = 0$, $f'(\bar{x}) = -1$, $f''(\bar{x}) = 0$, and $f^{(3)}(\bar{x}) = 48$, thus the theorem applies, meaning that zero is asymptotically stable fixed point of eq. (33).
<figure><img src="image_3.png" /><figcaption>Fig. 2. The flip manifold $W^f$ and the center manifold $W^c$ in transformed $(Y, Z)$ coordinates (top) and original $(y, z)$ coordinates (bottom). Points represent numerically computed orbits of a sample point on $W^c$ (○) and $W^f$ (●).</figcaption></figure>

We need to perform a similar analysis for the eigenvalue 1 and the corresponding center manifold $W^c$. Let us assume that $W^c$ has equation $Y = g(Z)$, where $g$ in the vicinity of 0 is given by the series $g(Z) = b_2Z^2 + b_3Z^3 + b_4Z^4 + b_5Z^5 + \dots$. The condition for invariance of $W^c$ requires that $Y_n = g(Z_n)$ remains valid at the next time step, $Y_{n+1} = g(Z_{n+1})$. Using our previous notation this means that $G_1(Y, Z) = g(G_2(Y, Z))$, which, by substituting $Y = g(Z)$, yields

$$G_1(g(Z), Z) = g(G_2(g(Z), Z)). \qquad (34)$$

As before, by expanding $G_1(g(Z), Z) = g(G_2(g(Z), Z))$ into the Taylor series and setting all coefficient of the expansion to be zero we obtain, using Maple, $b_2 = \frac{1}{2}$, $b_3 = -\frac{1}{2}$, $b_4 = -4$, $b_5 = -\frac{3}{2}$, etc. The equation of the center manifold is,
therefore,

$$Y = \frac{1}{2}Z^2 - \frac{1}{2}Z^3 - 4Z^4 - \frac{3}{2}Z^5 + O(Z^6). \qquad (35)$$

By substituting $Y_n$ by $g(Z_n)$ on the right hand side of eq. (25) and Taylor expanding again one obtains the equation describing the dynamics on the center manifold $W^c$,

$$Z_{n+1} = Z_n - 2Z_n^2 + 6Z_n^3 - 6Z_n^4 + O(Z_n^5). \qquad (36)$$

In order to determine the stability of 0 in the above difference equation, let us first define semistability. A fixed point $\bar{x}$ of $x_{n+1} = f(x)$ is said to be asymptotically semistable from the right if there exist $\delta > 0$ such that for any $x_0$ satisfying $x_0 - \bar{x} < \delta$ we have $\lim_{n \to \infty} x_n = \bar{x}$. One can show [3] that if $f'(\bar{x}) = 1$ and $f''(\bar{x}) < 0$ then $\bar{x}$ is asymptotically stable from the right. In our case, for eq. (36), we have $f(x) = x - 2x^2 + 6x^3 - 6x^4 + O(x^5)$, $\bar{x} = 0$, $f'(\bar{x}) = 1$ and $f''(\bar{x}) = -4 < 0$, thus for eq. (36), zero is asymptotically semistable from the right.

Figure 2 shows manifolds $W^f$ and $W^c$ together with sample orbits generated numerically by iterating eqs. (24) and (25). Direction of the flow is indicated by arrows. Note that $W^c$ is asymptotically semistable only on the right (for $Z > 0$), and unstable on the left (for $Z < 0$). The left-sided instability is irrelevant for us, since $Z$ represents the probability of 000 block, thus it must always be positive.

Since 0 is asymptotically stable on $W^f$, and asymptotically semistable on $W^c$, we conclude that for $Z_0 > 0$, $\lim_{n \to \infty}(X_n, Z_n) = (0, 0)$, or, equivalently, $\lim_{n \to \infty}(x_n, z_n) = (1/4, 0)$, as claimed in Proposition 3. Proposition 2 follows automatically. □

### 3 Quality of local structure approximation

We have demonstrated so far that for rule 14, the local structure approximation of level 3 reproduces correctly the limiting values of probabilities of blocks of length up to 3. What about the rate of convergence to these limiting values? In order to find this out, let us consider rates of convergence to zero of $P_n(000)$ and its approximation $z_n$. We know that $P_n(000) = 2^{-2n-3}(4n+3)C_n$, where $C_n = \frac{1}{n+1}\binom{2n}{n} = \frac{(2n)!}{n!(n+1)!}$. Using Stirling's formula for large $n$, $n! \sim \sqrt{2\pi n}\left(\frac{n}{e}\right)^n$, the Catalan number $C_n$ can be approximated as

$$C_n = \frac{1}{n+1} \frac{(2n)!}{(n!)^2} \sim \frac{1}{n+1} \frac{\sqrt{4\pi n} \left(\frac{2n}{e}\right)^{2n}}{\left(\sqrt{2\pi n} \left(\frac{n}{e}\right)^n\right)^2} = \frac{1}{n+1} \frac{2^{2n}}{\sqrt{\pi n}} = \frac{1}{n+1} \frac{4^n}{\sqrt{\pi n}},$$

meaning that $P_n(000)$ converges toward zero as a power law $P_n(000) \propto n^{-1/2}$, where $x \propto y$ means the ratio $x/y$ tends to a positive number as $n \to \infty$.

Let us now examine convergence of $z_n$ to 0. We do not have a formula for $z_n$, but we can generate $z_n$ numerically, by iterating the local structure equations. Figure 3 shows the graph of $z_n$ vs. $n$ in log-log coordinates together with the
graph of $P_n(000)$ vs. $n$. We can see that both graphs appear to be almost straight lines, confirming that both $z_n$ and $P_n(000)$ behave as $n^\alpha$ for large $n$. The difference is in the value of the exponent $\alpha$. For $P_n(000)$ the exponent (computed as a slope of the upper line in Figure 3) is $\alpha \approx -1/2$, whereas for $z_n$ the exponent (computed as a slope of the lower line) is $\alpha \approx -1$.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Plot of $P_n(000)$ (upper line) and its local structure approximation $z_n$ (lower line) as a function of $n$ in log-log coordinates.</figcaption></figure>

The value of the exponent $\alpha \approx -1$ can be explained as follows. The starting point of the local structure approximation orbit, $y_1 = 3/8$, $z_1 = 7/32$, lies almost on the center manifold $W_c$. The convergence toward the fixed point is, therefore, dominated by eq. (36), which, if we keep only leading terms, becomes $Z_n = Z_n - 2Z_n^2$. Although this equation is not solvable in a closed form, we can obtain its asymptotic solution using the standard technique used in the theory of iterations of complex analytic functions. We can namely conjugate the map $Z \to Z - 2Z^2$ with appropriate Möbius transformation, which moves the fixed point to $\infty[1, 2]$. In our case, the Möbius map will simply be the inverse, meaning that we change variables in the equation $Z_n = Z_n - 2Z_n^2$ to $u_n = 1/Z_n$, obtaining

$$u_{n+1} = u_n + 2 + \frac{4}{u_n - 2} \qquad (37)$$

Since $u_n \to \infty$, the above can be approximated for large $n$ by $u_{n+1} = u_n + 2$, which has the solution $u_n = 2t + u_0$, or, going back to the original variable, $Z_n = \frac{1}{2t+1/z_0}$. The result $z_n = Z_n \propto t^{-1}$ immediately follows.

In conclusion, one could thus say that the local structure approximation correctly reproduces not only the coordinates of the fixed point but also the
type of convergence toward the fixed point (as a power law). It fails, however, to reproduce the correct value of the exponent in the power law. This in agreement with the commonly reported results of investigations of critical phenomena: mean-field type theories cannot reproduce values of fractional exponents in power laws.

It would be interesting and beneficial to extend results of this paper to non-symmetric initial Bernoulli measures. Numerical evidence suggests that local structure approximation remains exact in the limit of $n \to \infty$ in such cases, but to be sure one would need to generalize eqs. (4)–(7) to non-symmetric initial measure. This, in principle, should be possible, and will be attempted in the future.

**Acknowledgement:** H.F. acknowledges financial support from the Natural Sciences and Engineering Research Council of Canada (NSERC) in the form of Discovery Grant.

## References

1. Beardon, A.F.: Iteration of Rational Functions. Springer, New York (1991)

2. Devaney, R.L.: An Introduction to Chaotic Dynamical Systems. Addison-Wesley, Reading, Massachusetts (1993)

3. Elyadi, S.N.: An Introduction to difference equations. Springer, New York (1999)

4. Fukś, H.: Remarks on the critical behavior of second order additive invariants in elementary cellular automata. Fundamenta Informaticae 78, 329–341 (2007)

5. Fukś, H.: Construction of local structure maps for cellular automata. J. of Cellular Automata 7, 455–488 (2013), (Updated version available at arXiv:1304.8035)

6. Fukś, H., Haroutunian, J.: Catalan numbers and power laws in cellular automaton rule 14. Journal of cellular automata 4, 99–110 (2009)

7. Fukś, H.: Orbits of Bernoulli measures in cellular automata. In: Meyers, R.A. (ed.) Encyclopedia of Complexity and Systems Science, pp. 1–19. Springer, Berlin, Heidelberg (2017), https://doi.org/10.1007/978-3-642-27737-5_676-1

8. Gutowitz, H.A., Victor, J.D., Knight, B.W.: Local structure theory for cellular automata. Physica D 28, 18–48 (1987)

9. Murakami, K.: Stability for non-hyperbolic fixed points of scalar difference equations. J. of Math. Analysis and Applications 310(2), 492–505 (2005)

10. Perko, L.: Differential Equations and Dynamical Systems. Springer, New York (2008)

11. Wolfram, S.: Cellular Automata and Complexity: Collected Papers. Addison-Wesley, Reading, Mass. (1994)