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

The one-parameter two-dimensional cellular automaton with the Margolus neighbourhood is analyzed based on the considering the projection of the stochastic movements of a single particle. Introducing the auxiliary random variable associated with the direction of the movement, we reduce the problem under consideration to the study of a two-dimensional Markov chain. The master equation for the probability distribution is derived and solved exactly using the probability generating function method. The probability distribution is expressed analytically in terms of Jacobi polynomials. The moments of the obtained solution allowed us to derive the exact analytical formula for the parametric dependence of the diffusion coefficient in the two-dimensional cellular automaton with the Margolus neighbourhood. Our analytic results agree with earlier empirical results of other authors and refine them. The results are of interest for the modelling two-dimensional diffusion using cellular automata especially for the multicomponent problem.

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 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 onedimensional 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.

II. DIFFUSION IN THE MARGOLUS CELLULAR AUTOMATON
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 × 2 cells. Such partition can be even or odd (see 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 • clockwise or counterclockwise with the probability of 0 < p ≤ 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 ∆t and p = const , then such MCA simulates the diffusion with the time independent diffusion coefficient 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 D c = 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 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 where r is the radius-vector of the particle position at time t = τ relative to the initial position at time t 0 = 0, 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 → ∞, the distribution of the particle position tends to the normal distribution with the probability density function given by where X n is the random position of the particle at the n-th time step, M Xn is the expectation of X n , and D Xn is the dispersion of X n . Then, for ∆x = 1, ∆t = 1, we have Hereinafter, D c will be given for ∆x = 1, ∆t = 1 unless otherwise stated. The trick is to find the distribution P Xn (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 D c (p) for an arbitrary p including the special case p = 0.5.

III. ONE-DIMENSIONAL MOVEMENT OF A SINGLE PARTICLE
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.

IV. MASTER EQUATION
Next, we give the formal definition of the two-dimensional Markov chain described in Section III.
and Ω 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 ∆ t is associated with its direction. The probability of transitions for the chain under consideration, which were described in Section III, are given by The semicolon in (4.1) implies the logical conjunction. Let us denote The transition rules (4.1) in notations (4.2) yield the following master equation: Hereinafter, x ∈ Z, d ∈ {−1, +1}, and t ∈ Z + . We consider the problem with the deterministic initial position of the particle, i.e. Prob(X 0 = 0) = 1. Then, the initial condition for where δ 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: Since the distribution over x is of primary interest to us, we will study the marginal distribution given by In order to construct the analytical description of the two-dimensional Markov chain under consideration, we will solve the master equation (4.3) with the initial condition (4.5) in the next Section.

V. PROBABILITY GENERATING FUNCTION
The exact solution of (4.3), (4.5) can be derived using the method of the probability generating function [24] for the marginal distribution (4.6). Such probability generating function reads Let us also introduce the following probability generating functions: In view of (4.6), one readily gets The master equation (4.3) can be written as (5.4) Multiplying the system (5.4) by z x and summing over x ∈ Z, we derive the following equation for the probability generating functions (5.2): where g 0 (z) is obtained from (4.5), (5.2).
The solution of the system (5.5) has the following form: where the solution of the associated matrix eigenvalue problem yields The solution (5.6), (5.7) can be rewritten as Raw moments of X t for the given ∆ t read The moments (5.9) can be obtained from (5.10) Using (5.10), (5.8), one readily gets Then, the expectation M Xt and the dispersion D Xt of the marginal distribution (4.6), which are given by While the formula (5.8) is convenient for the derivation of moments, it does not allow one to directly obtain the distribution (4.6). 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 (4.6) is given by the following theorem where P α,β n [x] are the Jacobi polynomials, and C k n = n! k!(n − k)! are the binomial coefficients.

(5.19)
The proof of this theorem is given in Appendix A.
Note that the probability p > 1 2 does not make sense for MCA. However, the definition (4.1) of the two-dimensional Markov chain under consideration admits p ∈ (0; 1) (p = 0 and p = 1 are trivial cases), and Theorem V.1 holds for p ∈ (0; 1). For completeness, we have shown the probability distribution for p > 1 2 in Fig. 5. Apart from the fact that p > 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 ≥ 0 for p > 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 = 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 We see that the dispersion D Xt (p) tends to its large time asymptotic value faster thañ D Xt (p s ). Moreover, for small t, it is also much closer to this value that determines the diffusion coefficient in MCA especially for p < 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 = 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 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 D c (p) (in relative units) on p. They constructed the following regression model for this dependence: Their experimental values for (6.5) fit well to the curve (6.1) 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 [7,26]. 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 ξ of particles that are involved in the rotation has the same form as the dependence (6.1) up to changing ξ = 2p.
We should note that, in some works related to the multicomponent MCA, reseachers falsely suppose that the dependence D c (p) should be linear for a wide range of p. This assumption originates from the work [28], where the authors claimed that 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 (6.1) shows that the linear approximation is inaccurate when p is substantially lower than 1 2 . Note that it also agrees with the cited numerical results [23].

VII. CONCLUSION
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 ∆ 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 ∆ 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 V.1) and for two first moments (5.13) 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 V.1 and formula (6.1). The formula (6.1) 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 (6.1) agrees with earlier empirical results of other authors in the specified cases. Thus, the formula (6.1) is a fundamental mathematical result for modelling the diffusion using MCA.
Theorem V.1 gives the probability distribution of a single particle in MCA that is given by the exact solutions to the two-dimensional Markov chain (4.3), (4.4). While (6.1) is of main applied interest, Theorem V.1 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 V.1, 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 ∈ (0; 1) while only p ∈ (0; 1 2 ] have a physical meaning in MCA. Therefore, the more detailed study of this process for the parameter p ∈ ( 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 (6.1). ACKNOWLEDGEMENT We are thankful to M.L. Gromov and N.A. Shalyapina for stimulating discussions. Let us write the probability generating function (5.8) as where q(z) = p z + 1 z , λ = 4(2p − 1), and m(z) is defined in (5.7). Next, we use the following properties: and arrive at the following formulae: Using the relation m 2 (z) = q 2 (z) − λ, the binomial theorem, and changing the order of summation, we transform (A3) to the power series in q(z): where the numbers Q, R are given by Using the binomial theorem for the power of q(z), changing the order of summation, in view of (5.3), we obtain the following identities: G 2n+1 (z) = 1 2 2n+2 n j=0 n l=j p 2l+1 a n,l C l−j 2l+1 z 2j+1 + where a n,l = (−1) n−l λ n−l R(n − 1, n − 1 − l) + R(n, n − l) , b n,l = (−1) n−l λ n−l Q(n, n − l), r n,l = (−1) n−l λ n−l Q(n − 1, n − 1 − l) + Q(n, n − l) , d n,l = (−1) n−l λ n−l R(n, n − l).

(A8)
In order to simplify (A8), let us prove the lemma below.
Lemma A.1. The following relation holds true: Assuming n to be a half-integer number, one readily gets (A16) from (A17) and (A18).