New approach to evaluate the asymptotic distribution of particle systems expressed by probabilistic cellular automata

We propose some conjectures on the asymptotic distribution of the probabilistic Burgers cellular automaton (PBCA), which is defined by a simple rule of particle motion with a probabilistic parameter. Asymptotic distribution of configurations converges to a unique steady state for PBCA. We propose a new and widely-applicable approach to analyze probabilistic particle systems and apply it concretely to PBCA and its extensions. We introduce a conjecture on the distribution and derive the asymptotic probability expressed by the GKZ hypergeometric function. If the space size goes into infinity, we can evaluate the relationship between the density and flux of particles for infinite space. Moreover, we propose two extended systems of PBCA and analyze their asymptotic behavior.


Introduction
Cellular automata (CA) are dynamical systems with discrete time, discrete space, and a finite set of state values. Their dynamics are generally determined by a simple rule that depends on the values of the neighboring space sites. This simple abstraction has interesting mathematical properties and has been studied in theoretical and applied research. For example, systems called "elementary cellular automata" (ECA) were classified according to the behavior of their solutions and have been extensively studied by many researchers [1][2][3][4].
There exist 256 independent rules for ECA. ECA with one of the rules, Rule 184, is known as a non-trivial particle system. It is also called Burgers cellular automaton (BCA) since it can be derived from Burgers equation by using the ultradiscretization method, which was discovered in the field of integrable systems [4]. There exists a threshold of the density of particles, and the asymptotic behavior of solutions drastically changes between the regions of lower and higher density. Thus, a phase transition occurs for BCA at the threshold. If we consider that the dynamics of BCA represents a transportation system with cars, this phase transition can be interpreted as a primitive model of occurrence of traffic jams.
We can introduce probabilistic parameters into the deterministic CA; this topic has also been extensively researched in theory and practice. For example, the asymmetric simple exclusive process (ASEP) is a well-known standard statistical model for a simple stochastic particle system. It is a random walk model of multiple particles, where each particle moves to the neighboring sites with an excluded volume effect. Various exact evaluations of statistical results have been completed for ASEP [5][6][7][8][9]. Sasamoto et al. revealed exact relations between ASEP and the orthogonal polynomials [10].
Moreover, several probabilistic CA represent realistic dynamical systems. For example, Nagel and Schreckenberg proposed a quite efficient dynamical model to investigate the physics of traffic jams, known as Nagel-Schreckenberg (NS) model [11]. They applied their model to real freeway traffic and obtained good accordance between the observed data and their theoretical estimation [12].
In previous publications, the author and his co-workers analyzed the asymptotic behavior of probabilistic CA with 4 neighbors or with a higher order of conserved quantities [13,14]. There exist some bilinear equations for probabilities of local patterns in the asymptotic solutions of the systems. Using these equations, they derived a theoretical expression of the relationship between density and flux for the asymptotic behavior, which is called a fundamental diagram (FD).
In this paper, we focus on a probabilistic extension of BCA (PBCA). It is partially equivalent to the "totally" asymmetric simple exclusion process (TASEP) obtained by restricting the motion of particles of ASEP to the only one direction [6,7,15]. However, PBCA and TASEP update particle positions differently. Although one of the particles to be updated is chosen at every time step in TASEP, the "parallelupdate" is used for PBCA: motions of all particles from the current time step to the next time step are determined simultaneously.
We report a new analysis to understand the asymptotic behavior of parallelupdated PBCA with the periodic boundary condition. We can consider PBCA as a one-dimensional random process and derive a transition matrix for the process. Assuming the random process is ergodic, we propose a conjecture on the asymptotic distribution of the system. Using the conjecture, we can derive the asymptotic probability of each configuration of particles in space sites and can derive FD from their expected values. Moreover, we give the expression of FD of PBCA by a hypergeometric function proposed by Gelfand, Kapranov, and Zelevinsky (GKZ). It is called the GKZ hypergeometric function and is obtained by extending the hypergeometric function of a single variable to multiple variables [16]. They obey the specific forms of differential equations, and their contiguity relation can be obtained in the form of New approach for particle systems of probabilistic cellular automata matrix [17][18][19]. The FD of PBCA with an infinite number of space sites is calculated by using the limit of the contiguity relations. This limiting case coincides with that of previous research based on a specific ansatz [11]. Furthermore, we propose two types of extensions of PBCA and derive the asymptotic probability of configurations and FD using a similar conjecture for PBCA.
The contents of this paper are as follows. In Sect. 2, we introduce definition and properties of PBCA, and we present a new type of analysis of PBCA. In Sect. 3, we propose two types of extensions of PBCA and present the results based on the similar conjecture of PBCA. One of the extensions is a system of which motion rule of particle depends on four neighbors. The other is an extension of PBCA to a coupled system. In Sect. 4, we give the conclusion. In Appendix A, we show calculation of PBCA using the GKZ hypergeometric function when the space size goes into infinity.

Definition of particle system
A probabilistic Burgers cellular automaton (PBCA) is defined by the following maxplus equation: where Subscript j is an integer site number, and superscript n is integer time. Real constant satisfies 0 ≤ ≤ 1 . The value of probabilistic parameter a n j is determined independently for every (j, n). We assume a periodic boundary condition for the space sites with a period L, that is, u n j+L = u n j . From the above evolution equation, it is easily shown that Here, u n j is the number of particles at site j and time n. Thus, the sum of all state values over L is conserved for n and is determined by the initial data. The rule of particle motion for PBCA can be expressed as follows: q n j = min(a n j , u n j , 1 − u n j+1 ), a n j = 1 (prob. ) 0 (1 − ) . (2) The particle at jth site moves to the right with probability only if no particle exists at ( j + 1)th site. Figure 1 shows an example of time evolution in a PBCA.

Asymptotic distribution of PBCA
Supposing that a set of values over all L sites corresponds to a configuration of a random process, PBCA is a one-dimensional random process on L C m configurations if m = ∑ L j=1 u n j , that is, if m particles move through the sites. Figure 2a shows a histogram of all configurations from n = 0 to n = 1000 obtained by a numerical computation. Figure 2b-d show those from n = 0 to n = 10,000, 100,000, and 1,000,000, respectively. These figures suggest that the distribution New approach for particle systems of probabilistic cellular automata of configurations converges after enough time steps. Moreover, the heights of bins can be divided into classes after enough time steps. We made various numerical calculations on the histogram and an approximately unique steady state is always obtained. Therefore, we can assume that PBCA is ergodic. To illustrate our assumption and to explain the relation among heights of classes, we show the exact results for some small L and m.
For L = 4 and m = 2 , a set of all configurations (which is denoted by ) is and transition probabilities are obtained as where a ij is a transition probability from the ith configuration to the jth configuration. The characteristic equation for the matrix has a simple root 1, and the absolute value of other roots is always less than 1. Subsequently, time evolution of the random process determined by this matrix has a unique steady state, and an eigenvector of the matrix for eigenvalue 1 becomes the asymptotic distribution of configurations. The corresponding eigenvector is and components of this eigenvector give the ratios of heights of the histogram.  and its eigenvector for eigenvalue 1 as a simple root of the characteristic equation is We propose the following conjecture from the exact results obtained for small L and the numerical results for large L and m. To give the conjecture, let us introduce some notations. Define as a set of all configurations for L and m, as the number of patterns s 1 s 2 ⋯ s k included in x considering the periodic boundary condition, and p(x) as the probability of x in the steady state.
Conjecture: For any x ∈ , we have where C is a normalization constant satisfying ∑ x∈ p(x) = 1. Considering the case of L = 4 and m = 2 as an example, we have Considering another case of L = 11 , and m = 6 , we have 000111 001110 011100 ⋮ 110100 101001 010101 101010 By the above-mentioned conjecture, we can categorize the probability by #10 . Then, the probability of a configuration x with #10(x) obtained in the limit n → ∞ for the space size L and the number of particles m is where N L,m (k) is the number of configurations with #10 = k defined by

Fundamental diagram of PBCA
A fundamental diagram (FD) of a general particle system is a diagram that shows the relation between the density and flux of the particles averaged over all sites in the limit n → ∞ . Since PBCA is also a particle system, we can derive its FD. Density is defined by m/L, where L and m are the number of sites and particles, respectively. Let Q L, (m) be the expected values of flux of the steady state, where is the hopping probability of a particle. Then, Q L, (m) of PBCA is an expected value of #10∕L , since only particles next to an empty site can move with probability . Thus, we have Figure 3 shows FD obtained by (3) and that by the numerical calculation. The former is shown by small black circles ( • ) and the latter by white circles ( ○ ). Their good coincidence can be observed from this figure.
Utilizing GKZ hypergeometric function, we can evaluate Q L, (m) in the limit of L → ∞ preserving m = L as which coincides with the result of previous research [11]. We show the proof of this derivation in the Appendix A [17][18][19].

Extended systems and their properties
In this section, we introduce two systems that are extensions of PBCA. We evaluate their asymptotic probability of configurations and derive their FD, as we have shown in the previous section. Note that we also assume the ergodicity for both systems.

Extension to 4 neighbors
Let us consider the probabilistic CA expressed by a max-plus equation, where The probabilistic parameters a n j and b n j are defined by Since this equation is also in the conservation form, the number of particles ( = ∑ j u n j ) is preserved during the time evolution. From the max-plus expression of flow q n j , we can easily show the motion rule of particles denoted by 1 as follows.
, a n j )). Note that any particle in the configuration other than 100 and 101 can not move. We call this system EPBCA1. In the case of = , any particle in the pattern "10" moves with probability . Thus, PBCA is included in EPBCA1 as a special case. Figure 4 shows an example of time evolution of the system.
We can obtain an exact form of the transition matrix and of eigenvector for eigenvalue 1 for small values of the space size L and the number of particles m. Then, we give a conjecture for asymptotic distribution of EPBCA1 from those results. Let us introduce two examples that support our conjecture.
Suppose a set of configurations for L = 8 and m = 4 . Size of is 70. To make the expression of the transition matrix compact, we define classes of configurations up to cyclic rotation and divide into the equivalence classes. The number of the equivalence classes is 10, and their representative elements are as follows.
where Since the mean flux Q L, , (m) is an expected value of ( #100 + #101)∕L , it becomes Figure 5 shows FD obtained by (7)and that by the numerical calculation. The former is shown by small black circles ( • ) and the latter by white circles ( ○ ). Their good coincidence can be observed from this figure. Although we have not yet obtained the limit distribution as L → ∞ , the profile of Q does not change drastically, for larger L.

Extension to 2 types of particles
Let us consider a probabilistic CA defined by the following system of max-plus equations. Assume u n j , v n j ∈ {0, 1} and probabilistic parameters a n j and b n j are defined by The above evolution rule can be interpreted as a probabilistic system of two types of particles, considering that u n j and v n j are the number of particles of type A and B at site j and time n, respectively. Moreover, let us assume these two types of particles do not exist at the same site and at the same time: (u n j , v n j ) ≠ (1, 1) . If we set this condition for the initial data, it is always satisfied during the time evolution according to the above evolution rule. Subsequently, we express the configuration of particles (u, v) = (0, 0) , (1, 0), and (0, 1) by symbols , , and , respectively. Then, the rule of particle motion is given as the follows: We call the above system EPBCA2. If v n j ≡ 0 , the motion rule for u n j of EPBCA2 reduces to that of PBCA, and that for v n j trivial. Thus, PBCA is included in this system as a special case. Figure 6 shows an example of time evolution in the EPBCA2. Any particle does not pass over the other particles according to the evolution rule. Let us define a "sequence" as an array of particles from a given configuration preserving their relative positions. For example, the sequence in configuration is , and this sequence may change into the sequences , , and during the time evolution. However, sequence can not change into since particles do not u n+1 j = u n j + q n j−1 − q n j , q n j = min(a n j , u n

Fig. 6
Example of time evolution in the EPBCA2 for = 0.4 and = 0.8 . White, grey, and black squares express (u, v) = (0, 0) (no particle), (1, 0) (only particle A exists), and (0, 1) (only particle B exists), respectively New approach for particle systems of probabilistic cellular automata pass over one another. Similarly, sequence may change into but not , , , or . If we consider a transition matrix of configurations for given initial data, the configurations are restricted to those obtained from the initial data. Therefore, if a set of configurations includes a configuration x = x 1 x 2 … x L , we restrict to be constructed from any combination of 's, 's, and 's preserving their numbers and the order of sequence included in x up to cyclic rotation. Subsequently, we give two examples of transition matrices and their eigenvectors.
The first example is that includes configuration . The size of is 84. To make the expression of the transition matrix compact, we define classes of configurations up to cyclic rotation and divide into the equivalence classes. The number of equivalence classes is 12, and their representative elements are as follows: The transition matrix derived by transition probabilities from each representative element to the others is and its eigenvector for eigenvalue 1 is Another example is that includes configuration . The size of is 140, the number of equivalence classes is 20, and their representative elements are The transition matrix is and the eigenvector for eigenvalue 1 is , , , . and … , and let n be that of … and … . For example, if configuration is given, then k = 5 , k = 4 , n = 6 , and n = 3 . Note that m + m + n + n = L . Using these notations, we introduce the following conjecture.

Conjecture: The probability of configuration x in the steady state is
where C is a normalization constant satisfying ∑ x∈ p(x) = 1. From the above conjecture, the probability of any configuration x for k , k , n , and n in the steady state is given by where N(k , k , n , n ) denotes the number of all configurations included in with k , k , n , and n and is calculated as Considering the motion rule of particles A and B, the flux in the steady state is New approach for particle systems of probabilistic cellular automata Figure 7a shows FD calculated by (9) and Fig. 7b shows the numerical result. In these figures, A and B are densities of particles A and B, respectively. Figure 8 shows a comparison of FDs of Fig. 7 for = 0.5 . Small black circles ( • ) are obtained by (9) and white circles ( ○ ) by the numerical calculation. Their good coincidence can be observed from this figure.

Concluding remarks
We introduced some conjectures on the asymptotic distribution of PBCA and its extensions, assuming that these systems are ergodic. In the conjectures, the following points are the most important.
• Probability of any configuration in the steady state depends only on the number of some specific local patterns included in the configuration. For example, it depends only on #10 for PBCA. • Since the probabilities do not depend on the location of patterns in the configuration, they are equal to each other if the numbers of patterns included in the configurations are also equal. For example, the probabilities of configurations are the same if their #10 are the same in the case of PBCA.
Based on the conjecture, we derived the probability of configuration in the steady state and FD of the system for arbitrary size of L.
On the other hand, FD of PBCA for L → ∞ is reported in the previous research as using an ansatz on relations of probabilities of local patterns [11]. Using the GKZ hypergeometric function, we evaluate the limit of FD of PBCA (3) and confirm that the diagram coincides with the above result.
Furthermore, since FD of the extended systems is also expressed by some kind of the GKZ hypergeometric function, we expect that FD in the limit of infinite size can be derived similarly as PBCA. In particular, for EPBCA2, for the initial condition of sequence … … and a probabilistic parameter = 1 , the rule of particle motion in the steady state becomes as follows. This motion rule is the same as that of the L-R system obtained by the ultradiscrete Cole-Hopf transformation of SECA84 with a quadratic conserved quantity [14]. The FD of L-R system in the limit of infinite space size is evaluated and derived in a simple form depending on the density of particles. Therefore, FD of EPBCA2 in the limit of L → ∞ includes FD of SPCA84 as a special case.
New approach for particle systems of probabilistic cellular automata Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/ licenses/by/4.0/.

Definition of the GKZ hypergeometric function and its properties
In this section, we define the GKZ hypergeometric function and its properties. First, let us define B k (x) by It satisfies the following contiguity relations.
For an m × n matrix A = (a ij ) and an n-dimensional vector , general solution to is expressed by where is a special solution to the equation. Using A, , and , the GKZ hypergeometric function ( ; ) is defined by where is an n-dimensional vector and 0,j , k j and x j are jth element of , , and , respectively. Although also depends on the matrix A, A is fixed below and we omit the dependency. Using contiguity relations, the following two properties for GKZ hypergeometric function are derived.
First, for a given vector J + ( ) and J − ( ) are defined by subsets of {0, 1, … , n} as Then, we obtain Therefore, satisfies a differential equation Second, for an operator j = x j x j and defined by satisfies Using ∈ ℂ × and the ith row vector = (a i1 , … , a in ) of A, define D( ) by Since Thus, is obtained.

FD of PBCA expressed by the GKZ hypergeometric function
In this subsection, we express FD of PBCA by the GKZ hypergeometric function. The FD of PBCA (3) for the infinite space size is where density is constant and To evaluate the above limit, we choose the following matrix A and vector in the definition of the GKZ hypergeometric function in A.1.
On the other hand, derivative of (17) with respect to x 3 is Therefore, from (21)

3
New approach for particle systems of probabilistic cellular automata