Analytical description of the diffusion in a cellular automaton with the Margolus neighbourhood in terms of the two-dimensional Markov chain

Cellular automata (CA) are the powerful tool for modelling the matter transport [1]. The attractive feature of CA is that contemporary algorithms of parallel computing can be naturally embedded in them. Various cellular algorithms were successfully implemented on field-programmable gate arrays (FPGA) [2, 3] and on graphic cards [4, 5].

In the study of diffusion phenomena in physics, chemistry, and biology, CA are positioned as one of the effective tools for modelling systems of many particles. The simplest CA with boolean alphabet modelling the diffusion are the asynchronous CA with naive diffusion and the synchronous CA with Margolus neighbourhood [6] or briefly Margolus CA (MCA). The last one is the object of our study. Exactly this type of the CA attracted our attention since synchronous CA are more advantageous in terms of the computing speed (see, e.g., comparison in [7]). Note that the MCA is the representative of partitioning CA.

The most common macroscopic characteristics of the diffusion transfer is the diffusion coefficient ${\mathcal{D}}_{c}$ that is defined as the proportionality factor in Fick‘s law [8, 9]. The diffusion coefficient is usually invariant when the movements of particles are caused just by random fluctuations rather than by any fields. Nevertheless, the processes with the time dependent diffusion coefficient are also of interest, e.g., in the ambipolar diffusion phenomena [10, 11, 12]. The most popular model of the diffusion is the kinetic model, which is based on a differential equations that can be solved, for example, using difference schemes [13, 14, 15]. In such approach, the diffusion coefficient is an explicit parameter of the model. For cellular automata, the diffusion coefficient is not presented in the model explicitly. Calculating the diffusion coefficient for the CA with naive diffusion is a trivial task since it directly models the random walk of particles that is well-studied. However, it is not trivial for partitioning CA including MCA. We show here that the diffusion in the two-dimensional MCA can be described in terms of one two-dimensional Markov chain as opposed to the one-dimensional random walk. In our paper, we study the exact solution of this two-dimensional Markov chain whose characteristics are crucial for the mathematical description of the diffusion in MCA.

While multidimensional Markov chains arise in a number of applications in service theory [16], genetic networks [17], Gibbs sampling [18], and other areas [19], there are very few exactly solvable multidimensional Markov chains that do not degenerate to the one-dimensional one. Therefore, the exact solutions to the two-dimensional Markov chain under consideration can be of interest for specialists in applied mathematics and useful for applications other than cellular automaton theory.

The paper is organized as follows. In Section II, we describe the diffusion model within the framework of the two-dimensional cellular automata with Margolus neighbourhood. Also, we explain the way of obtaining the diffusion coefficient for this model based on the considering one-dimensional movement of a single particle. In Section III, the stochastic process describing the one-dimensional movement of a single particle is given. The idea of introducing the auxiliary random variable associated with the direction of movement is proposed. This approach allows us to consider the process under consideration as the two-dimensional Markov chain. In Section IV, the formal mathematical definition of this chain and basic notations are given. In Section V, we solve the master equation for the probability distribution of the particle coordinate using the probability generating function method. The moments of this distribution are derived. In Section VI, the obtained results are discussed within the framework of the diffusion model. The exact analytical formula for the diffusion coefficient is derived, and properties of the cellular automaton with the Margolus neighbourhood are discussed. In Section VII, we conclude with some remarks.

In the two-dimensional cellular automata with Margolus neighbourhood, the diffusion is described based on the following algorithm [6]:
1) The plane is divided into identical square cells forming a cell grid. One cell contains a boolean variable that takes the value 0 or 1 corresponding to the absence or presence of a particle in the cell respectively.
2) The cell grid is divided into blocks of $2\times 2$ cells. Such partition can be even or odd (see Fig. 1). The even partition corresponds to the case when the left lower cell of each block has even vertical and horizontal numbers. The odd partition corresponds to the case when the left lower cell of each block has odd vertical and horizontal numbers.
3) Beginning from an initial state, on each discrete time step the following rule is applied to each block independently. We rotate the block $90^{\circ}$ clockwise or counterclockwise with the probability of $0<p\leq 0.5$ for each outcome. The rule is applied to blocks of the even partition on even time steps and it is applied to blocks of the odd partition on odd time steps. If each time step corresponds to the same time interval $\Delta t$ and $p=\mathop{\rm const}\nolimits\,$, then such MCA simulates the diffusion with the time independent diffusion coefficient ${\mathcal{D}}_{c}$.

The classic MCA uses $p=0.5$, i.e. blocks always rotate. In such approach, the diffusion rate is controlled with the length of the time interval corresponding to one discrete time step and with the size of a cell. However, this model can be applied to the multicomponent diffusion problem. In this case, the cell contains few bits that correspond to particles of different substances. It can be treated as few layers of the cell grid. Then, we can not directly apply the described rule with $p=0.5$ for every layer of particles if the diffusion rates for these substances are not equal. The difference in diffusion rates can be realized using different $p$ for different substances. The only drawback of such approach is that the diffusion coefficient ${\mathcal{D}}_{c}={\mathcal{D}}_{c}(p)$ nonlinearly depends on $p$. The other way is using the same $p=0.5$ for each layer but skip the two consecutive time steps for one or few layers with a probability equal to $p_{s}$ (the rotation rule is not applied to any blocks of the respective layer on the skipped time step). In this approach, the diffusion coefficient linearly depends on $p_{s}$. The differences between these two approaches will be discussed in Section VI.

In order to model the diffusion using the MCA, the diffusion coefficient corresponding to the particular MCA must be determined. It can be shown (see, e.g., [20]) that

$${\mathcal{D}}_{c}(p)=k(p)\displaystyle\frac{\Delta x^{2}}{\Delta t},$$

(1)

where ${\mathcal{D}}_{c}(p)$ is the diffusion coefficient, $\Delta x$ is the side length of the cell, $\Delta t$ is the time interval corresponding to the discrete time step, and $k(p)$ is the proportionality factor. However, relation (1) does not help to obtain the exact value of the multiplier $k$. Moreover, $k$ nonlinearly depends on $p$ as it was mentioned earlier. One way is to determine $k$ is to compute the specific diffusion problem using the numerical realization of the MCA and to compare the result with the analytical or numerical solutions of the respective diffusion equation. Such approach allows one to find $k$ approximately. However, it is hard to estimate error of obtained value. Moreover, this error depends on the choice of the specific diffusion problem. The other way is to obtain the distribution of the particle position in the MCA. Assuming that the correlation of these distributions for two particles tends to zero with the growth of time steps number, the distribution of the position of the single particle is equal to the Green function of the diffusion equation that is given by

$$G(\vec{r},\tau)\sim\exp\left(-\displaystyle\frac{\vec{r}\,^{2}}{4{\mathcal{D}}_{c}\tau}\right),$$

(2)

where $\vec{r}$ is the radius-vector of the particle position at time $t=\tau$ relative to the initial position at time $t_{0}=0$, ${\mathcal{D}}_{c}$ is the diffusion coefficient. For the isotropic and homogeneous problem, one can consider the projection of the particle position on one axis (let it be $x$). As $t\to\infty$, the distribution of the particle position tends to the normal distribution with the probability density function given by

$$P_{X_{n}}(x)\sim\exp\left(-\displaystyle\frac{(x-M_{X_{n}})^{2}}{2D_{X_{n}}}\right),\qquad t=n\Delta t,\qquad n\in{\mathbb{N}},$$

(3)

where $X_{n}$ is the random position of the particle at the $n$-th time step, $M_{X_{n}}$ is the expectation of $X_{n}$, and $D_{X_{n}}$ is the dispersion of $X_{n}$. Then, for $\Delta x=1$, $\Delta t=1$, we have

$$\begin{gathered}M_{X_{n}}\to 0,\qquad n\to\infty,\\ \displaystyle\frac{D_{X_{n}}}{2n}\to{\mathcal{D}}_{c},\qquad n\to\infty.\end{gathered}$$

(4)

Hereinafter, ${\mathcal{D}}_{c}$ will be given for $\Delta x=1$, $\Delta t=1$ unless otherwise stated. The trick is to find the distribution $P_{X_{n}}(x)$ or its moments analytically. In [21], the authors tried to obtain the diffusion coefficient using this approach for $p=0.5$. However, they assumed that $X_{n}$ is the one-dimensional Markov chain, i.e. the simplest random walk. Later on, using the software realization of MCA with higher number of cells and time steps, it was shown that this assumption led to that the diffusion coefficient obtained in [21] for 2D MCA is approximately three times higher than the real one (see, e.g., [20, 22, 23]). However, the exact distribution of $X_{n}$ and the diffusion coefficient of MCA were not obtained since then. In this work, we present the exact distribution of $X_{n}$ in MCA and obtain the exact analytical formula of ${\mathcal{D}}_{c}(p)$ for an arbitrary $p$ including the special case $p=0.5$.

In this Section, we will describe the rules of the stochastic movement of a single particle in MCA. As it was mentioned in the previous Section, we consider the movement along the $x$-axis that is directed along the side of a cell. Notice that if a particle is in an odd column at an odd time step, then its $x$-coordinate can either increase by 1 or do not change according to the rotation rule of the MCA (see Fig. 2). The same rule works for a particle in a even column at an even time step. Let the variable $d$ show the direction of possible movement of the particle along the $x$-axis. If $d=1$, then $x$ can either increase by 1 or do not change. If $d=-1$, then $x$ can either decrease by 1 or do not change. It is easy to see that $d=-1$ corresponds to a particle in an even column at an odd time step and to a particle in an odd column at an even time step.

Notice that the direction $d$ does not change when the $x$-coordinate changes and the direction $d$ reverts its sign when the coordinate $x$ does not change. The first outcome is realized with the probability of $p$ at every time step while the second outcome is realized with the probability $(1-p)$. Note that the sequence of the $x$-coordinates of the single particle has a memory, i.e. the current $x$-coordinate of the particle depends the history of its movements. It the main issue in the analytical description of partitioning cellular automata.

Next, we give the formal definition of the two-dimensional Markov chain described in Section III.

Let $X_{t}:\Omega\to{\mathbb{Z}}$ and $\Delta_{t}:\Omega\to\{-1,+1\}$ be sequences of random variables where $t\in{\mathbb{Z}}_{+}$ and $\Omega$ is a set of outcomes. We associate the values of $X_{t}$ with the discrete $x$-coordinate of the particle on a discrete time step $t$, and $\Delta_{t}$ is associated with its direction. The probability of transitions for the chain under consideration, which were described in Section III, are given by

$$\begin{gathered}{\rm Prob}(X_{t+1}=X_{t}+\Delta_{t};\Delta_{t+1}=\Delta_{t})=p,\\ {\rm Prob}(X_{t+1}=X_{t};\Delta_{t+1}=-\Delta_{t})=1-p.\end{gathered}$$

(5)

The semicolon in (5) implies the logical conjunction. Let us denote

$$P_{t}(x,d)={\rm Prob}(X_{t}=x;\Delta_{t}=d).$$

(6)

The transition rules (5) in notations (6) yield the following master equation:

$$P_{t+1}(x,d)=pP_{t}(x-d,d)+(1-p)P_{t}(x,-d).$$

(7)

Hereinafter, $x\in{\mathbb{Z}}$, $d\in\{-1,+1\}$, and $t\in{\mathbb{Z}}_{+}$. We consider the problem with the deterministic initial position of the particle, i.e. ${\rm Prob}(X_{0}=0)=1$. Then, the initial condition for (7) reads

$$P_{0}(x,+1)=\varepsilon\cdot\delta_{x0},\qquad P_{0}(0,-1)=(1-\varepsilon)\cdot\delta_{x0},$$

(8)

where $\delta_{ij}$ is the Kronecker delta.

Within the framework of the diffusion problem, it is natural to consider the symmetric case $P_{0}(0,+1)=P_{0}(0,-1)$. Thus, the following initial condition will be considered:

$$P_{0}(x,+1)=\displaystyle\frac{\delta_{x0}}{2},\qquad P_{0}(x,-1)=\displaystyle\frac{\delta_{x0}}{2}.$$

(9)

Since the distribution over $x$ is of primary interest to us, we will study the marginal distribution given by

$$P_{t}(x)=P_{t}(x,+1)+P_{t}(x,-1).$$

(10)

In order to construct the analytical description of the two-dimensional Markov chain under consideration, we will solve the master equation (7) with the initial condition (9) in the next Section.

The exact solution of (7), (9) can be derived using the method of the probability generating function [24] for the marginal distribution (10). Such probability generating function reads

$$G_{t}(z)=\sum_{x=-\infty}^{\infty}P_{t}(x)z^{x}.$$

(11)

Let us also introduce the following probability generating functions:

$$G_{t}^{+}(z)=\sum_{x=-\infty}^{\infty}P_{t}(x,+1)z^{x},\qquad G_{t}^{-}(z)=\sum_{x=-\infty}^{\infty}P_{t}(x,-1)z^{x}.$$

(12)

In view of (10), one readily gets

$$G_{t}(z)=G_{t}^{+}(z)+G_{t}^{-}(z).$$

(13)

The master equation (7) can be written as

$$\left\{\begin{array}[]{l}P_{t+1}(x,+1)=pP_{t}(x-1,+1)+(1-p)P_{t}(x,-1),\cr P_{t+1}(x,-1)=pP_{t}(x+1,-1)+(1-p)P_{t}(x,+1).\end{array}\right.$$

(14)

Multiplying the system (14) by $z^{x}$ and summing over $x\in{\mathbb{Z}}$, we derive the following equation for the probability generating functions (12):

$$\begin{gathered}g_{t+1}(z)=b(z)g_{t}(z),\qquad g_{0}(z)=\displaystyle\frac{1}{2}\begin{pmatrix}1\\ 1\end{pmatrix},\\ g_{t}(z)=\begin{pmatrix}G^{+}_{t}(z)\\ G^{-}_{t}(z)\end{pmatrix},\qquad b(z)=\begin{pmatrix}pz&1-p\\ 1-p&pz^{-1}\end{pmatrix},\end{gathered}$$

(15)

where $g_{0}(z)$ is obtained from (9), (12).

The solution of the system (15) has the following form:

$$g_{t}(z)=u(z)\begin{pmatrix}\lambda_{1}^{t}(z)&0\\ 0&\lambda_{2}^{t}(z)\end{pmatrix}u^{-1}(z)g_{0}(z),\qquad u(z)u^{\top}(z)=u^{\top}(z)u(z)=1,$$

(16)

where the solution of the associated matrix eigenvalue problem yields

$$\begin{gathered}\lambda_{1}(z)=\displaystyle\frac{p}{2}\left(z+\displaystyle\frac{1}{z}\right)+\displaystyle\frac{m(z)}{2},\qquad\lambda_{2}(z)=\displaystyle\frac{p}{2}\left(z+\displaystyle\frac{1}{z}\right)-\displaystyle\frac{m(z)}{2},\\ u(z)=\begin{pmatrix}c_{1}(z)(1-p)&c_{2}(z)(1-p)\\ \displaystyle\frac{c_{1}(z)}{2}\left(m(z)-r(z)\right)&-\displaystyle\frac{c_{2}(z)}{2}\left(m(z)+r(z)\right)\end{pmatrix},\qquad m(z)=\sqrt{p^{2}\left(z+\displaystyle\frac{1}{z}\right)^{2}-8p+4},\\ r(z)=p\left(z-\displaystyle\frac{1}{z}\right),\qquad c_{1}(z)=\sqrt{\displaystyle\frac{m(z)+r(z)}{2(1-p)^{2}m(z)}},\qquad c_{2}(z)=\sqrt{\displaystyle\frac{m(z)-r(z)}{2(1-p)^{2}m(z)}}.\end{gathered}$$

(17)

The solution (16), (17) can be rewritten as

$$\begin{gathered}G^{+}_{t}(z)=\displaystyle\frac{1}{2m(z)}\left[-\lambda_{1}^{t}(z)\lambda_{2}(z)+\lambda_{1}(z)\lambda_{2}^{t}(z)+(1-p+pz)\left(\lambda_{1}^{t}(z)-\lambda_{2}^{t}(z)\right)\right],\\ G^{-}_{t}(z)=\displaystyle\frac{1}{2m(z)}\left[\lambda_{1}^{t+1}(z)-\lambda_{2}^{t+1}(z)+(1-p-pz)\left(\lambda_{1}^{t}(z)-\lambda_{2}^{t}(z)\right)\right].\end{gathered}$$

(18)

Raw moments of $X_{t}$ for the given $\Delta_{t}$ read

$$\mu^{(n)\pm}_{t}=\sum_{x\in{\mathbb{Z}}}x^{n}P_{t}(x,\pm 1).$$

(19)

The moments (19) can be obtained from

$$\mu^{(n)\pm}_{t}=\displaystyle\frac{d^{n}G_{t}^{\pm}(z)}{dz^{n}}\Big{|}_{z=1}.$$

(20)

Using (20), (18), one readily gets

$$\begin{gathered}\mu^{(1)+}_{t}=-\mu^{(1)-}_{t}=\displaystyle\frac{p}{4-4p}\left(1-(2p-1)^{t}\right),\\ \mu^{(2)+}_{t}=\mu^{(2)-}_{t}=\displaystyle\frac{p^{2}}{4(1-p)^{2}}\left(-1+\displaystyle\frac{2(1-p)}{p}t+(2p-1)^{t}\right).\end{gathered}$$

(21)

Then, the expectation $M_{X_{t}}$ and the dispersion $D_{X_{t}}$ of the marginal distribution (10), which are given by

$$M_{X_{t}}=\sum_{x\in{\mathbb{Z}}}xP_{t}(x),\qquad D_{X_{t}}=\sum_{x\in{\mathbb{Z}}}(x-M_{X_{t}})^{2}P_{t}(x),$$

(22)

read

$$M_{X_{t}}=0,\qquad D_{X_{t}}=\displaystyle\frac{p^{2}}{2(1-p)^{2}}\left(-1+\displaystyle\frac{2(1-p)}{p}t+(2p-1)^{t}\right).$$

(23)

While the formula (18) is convenient for the derivation of moments, it does not allow one to directly obtain the distribution (10). The explicit representation of the probability generating function as the power series in $z$ is a bit tricky. The analytical form for the marginal distribution (10) is given by the following theorem

The proof of this theorem is given in Appendix A.

Now, from (4) and (23), we can obtain the diffusion coefficient ${\mathcal{D}}_{c}(p)$:

$${\mathcal{D}}_{c}(p)=\lim_{t\to\infty}\displaystyle\frac{D_{X_{t}}}{2t}=\displaystyle\frac{1}{2}\displaystyle\frac{p}{1-p}.$$

(30)

In Fig. 3, 4, 5, the probability distribution from (24), (25), (26), (27) (or (29) for $p=\frac{1}{2}$), and (28) is plotted along with the following normal probability distribution:

$$f_{t}(x)=\sqrt{\displaystyle\frac{1-p}{2t\pi p}}\exp\left(-\displaystyle\frac{(1-p)x^{2}}{2tp}\right).$$

(31)

Note that in view of (3), (4) the distribution $P_{t}(x)$ tends to $f_{t}(x)$ as $t\to\infty$ at least for $p\in[0;\frac{1}{2}]$. The moments of the normal distribution (31) are given by (30), (23).

From Fig. 3, 4, it is clear that the distribution is very close the normal one even for small $t$.

Note that the probability $p>\displaystyle\frac{1}{2}$ does not make sense for MCA. However, the definition (5) of the two-dimensional Markov chain under consideration admits $p\in(0;1)$ ($p=0$ and $p=1$ are trivial cases), and Theorem 29 holds for $p\in(0;1)$. For completeness, we have shown the probability distribution for $p>\displaystyle\frac{1}{2}$ in Fig. 5. Apart from the fact that $p>\displaystyle\frac{1}{2}$ has no physical meaning regarding the diffusion, we can note that the probability distribution $P_{t}(x)$ is a nonmonotonic function with respect to $x\geq 0$ for $p>\displaystyle\frac{1}{2}$.

Next, we will discuss the differences between the two types of MCA mentioned in Section II. Let the MCA parameterized by $p$ be termed the first type MCA. The second type MCA corresponds to the $p=\displaystyle\frac{1}{2}$ but it is parameterized by $p_{s}$. The parameter $p_{s}$ determines the probability that the rotation rule does not applies to any block of the odd and even partition at the next two time steps. In Section II, we described it as the probability of skipping two consecutive time steps. We will limit ourselves to the simple case when the probability $p_{s}$ is checked only at odd time steps. Note that $p_{s}$ corresponds to the global rule while the probability $p$ corresponds to the local rule since it is checked for each individual block independently.

We have shown that the diffusion coefficient ${\mathcal{D}}_{c}(p)$ (30) in the first type MCA nonlinearly depends on $p$. On the other hand, the diffusion coefficient $\tilde{{\mathcal{D}}_{c}}(p_{s})$ in the second type MCA linearly depends on $p_{s}$ and is given as follows:

$$\tilde{{\mathcal{D}}_{c}}(p_{s})=(1-p_{s}){\mathcal{D}}_{c}(0.5).$$

(32)

Both of these MCA can be adjusted to the same diffusion coefficient by the appropriate choice of $p$ and $p_{s}$, i.e. they have the same macroscopic behaviour at large times. However, they behave differently at small times $t$. Since the first and the third moments of the distribution given by (24), (25), (26), (27), (28) equal to zero, the core information about the process behaviour is carried by the dispersion. The behaviour of the first type MCA corresponds to the dispersion $D_{X_{t}}(p)$ given by (23) while the behaviour of the second type MCA corresponds to the dispersion $\tilde{D}_{X_{t}}(p_{s})$ given by

$$\tilde{D}_{X_{t}}(p_{s})=(1-p_{s})\cdot D_{X_{t}}(0.5).$$

(33)

In Fig. 6, the time dependence of $D_{X_{t}}(p)$ and $\tilde{D}_{X_{t}}(p_{s})$ are shown.

We see that the dispersion $D_{X_{t}}(p)$ tends to its large time asymptotic value faster than $\tilde{D}_{X_{t}}(p_{s})$. Moreover, for small $t$, it is also much closer to this value that determines the diffusion coefficient in MCA especially for $p<\displaystyle\frac{1}{2}$. It means that the first type MCA shows its macroscopic diffusion behaviour at smaller $t$. That can be treated as a better time resolution of such MCA. Note that the comparison of the dispersions yields just an lower estimate for the time resolution since we do not take into account correlation between the movement of different particles in MCA. However, there is a reason to believe that this correlation is even higher for the second type MCA with $p_{s}\neq 1$ since the skipping of the two consecutive time steps is performed for the whole cellular grid as opposite to the first type MCA, where it is checked for the every block whether it rotates or not on the current time step. Hence, the difference mentioned above may be even greater than one obtained from our estimates.

Let us compare our analytical results with some empirical results of other authors.

In [20], the author has compared the results of modelling with MCA and solutions of PDE numerically. The diffusion coefficient ${\mathcal{D}}_{c}(0.5)=0.512$ was obtained. The error of 5% was declared in [20]. Their numerical results agrees with the exact value of $0.5$ obtained in this work.

In [23], the authors have empirically obtained the dependence of ${\mathcal{D}}_{c}(p)$ (in relative units) on $p$. They constructed the following regression model for this dependence:

$$p=-0.35\left(\displaystyle\frac{{\mathcal{D}}_{c}(p)}{{\mathcal{D}}_{c}(0.5)}\right)^{2}+0.86\displaystyle\frac{{\mathcal{D}}_{c}(p)}{{\mathcal{D}}_{c}(0.5)}.$$

(34)

Their experimental values for (34) fit well to the curve (30) taking into account the inaccuracy of their method of matching the solutions of MCA and PDE. Hence, they could obtain the regression model that is close to the exact formula if they would use the rational function approximation, which is widely used in many applications [25], instead of the polynomial model.

Note that along with the MCA with the Boolean alphabet the one with the integer alphabet is also used [26, 7]. Such approach allows one to lower the concentration noise [27] when modelling reaction-diffusion systems. In order to control the diffusion coefficient independently of reaction rates in such MCA, the rotation rule applies not to the whole integer value but to its percentage. Applying the rotation rule to the percentage of particles is equivalent to the changing the block rotation probability $p$ by the same factor in terms of the diffusion coefficient. Thus, the dependence of the diffusion coefficient on the percentage $\xi$ of particles that are involved in the rotation has the same form as the dependence (30) up to changing $\xi=2p$.

We should note that, in some works related to the multicomponent MCA, reseachers falsely suppose that the dependence ${\mathcal{D}}_{c}(p)$ should be linear for a wide range of $p$. This assumption originates from the work [28], where the authors claimed that ${\mathcal{D}}_{c}(p)$ is linear on $p$ even when the rotation probability is independent for each block. Actually, they considered the modified version of the MCA that is neither the MCA with boolean alphabet, nor the multicomponent MCA. Our exact analytical formula (30) shows that the linear approximation is inaccurate when $p$ is substantially lower than $\frac{1}{2}$. Note that it also agrees with the cited numerical results [23].

We have considered the one-parameter stochastic process describing the one-dimensional movement (projection on the $x$-axis) of a single particle in two-dimensional MCA. In order to find the distribution of the $x$-coordinate of the particle at every discrete time step (the probability distribution of the random variable $X_{t}$), we have introduced the additional random variable $\Delta_{t}$ associated with the direction of the movement of the particle. Using our approach, we have reduced this process to the two-dimensional Markov chain of $X_{t}$ and $\Delta_{t}$ and have obtained the desired distribution by the probability generating function method. In particular, we have derived the analytic formula for the probability distribution (Theorem 29) and for two first moments (23) of $X_{t}$. The probability distribution is expressed in compact form in terms of Jacobi polynomials.

The primary results of this work are the Theorem 29 and formula (30). The formula (30) gives the diffusion coefficient of the one-parameter two-dimensional MCA. It establishes the correspondence between the cellular automaton model of the diffusion and the kinetic theory. Note that in our formalist diffusion coefficient is given by the asymptotic behaviour of the dispersion of $X_{t}$ for the distribution under consideration. It is shown that the formula (30) agrees with earlier empirical results of other authors in the specified cases. Thus, the formula (30) is a fundamental mathematical result for modelling the diffusion using MCA. Theorem 29 gives the probability distribution of a single particle in MCA that is given by the exact solutions to the two-dimensional Markov chain (7), (8). While (30) is of main applied interest, Theorem 29 is a more general fundamental result. It allows us to study specific features of various realizations of MCA discussed in Section 6. Moreover, the exactly solvable two-dimensional Markov chain are also valuable for the fundamental mathematics. In view of Theorem 29, the empirically proven fact that MCA describes the diffusion can be presented as a specific property of Jacobi polynomials (see Fig. 3, 4).

Due to the potential interest for other applications, it is worth to study other properties of the process under consideration. Note that the formulae derived in this work are valid for the range $p\in(0;1)$ while only $p\in(0;\frac{1}{2}]$ have a physical meaning in MCA. Therefore, the more detailed study of this process for the parameter $p\in(\frac{1}{2};1)$ is planned in a separate future work. Also, it is of interest to generalize our approach for the three-dimensional MCA, especially the formula (30).

We are thankful to M.L. Gromov and N.A. Shalyapina for stimulating discussions.