{'fraction_non_alphanumeric': 0.05241888113531167, 'fraction_numerical': 0.017903714668358356, 'mean_word_length': 3.4133139776026544, 'pattern_counts': {'":': 0, '<': 44, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Limiting shapes for a non-abelian sandpile growth model and related cellular automata

We present limiting shape results for a non-abelian variant of the abelian sandpile growth model (ASGM), some of which have no parallel in the ASGM. One of our limiting shapes is an octagon. In our model, mass spreads from the origin by the toppling rule in Zhang's sandpile model. Previously, several limiting shape results have been obtained for the ASGM using abelianness and monotonicity as main tools. As both properties fail for our model, we use a new proof technique: in our main proof, we introduce several cellular automata to mimic our growth model.


Introduction
We consider the following setup. Start with a pile of mass n ≥ 1 at the origin of the rectangular grid Z d , and mass h < 1 at every other site, where n and h are real numbers. Mass may be moved by 'splitting piles', that is, one may take all the mass from one site, and divide it evenly among its neighbors. One can only split 'unstable' piles of mass at least 1, and one can not stop before all the piles have mass less than 1. We call T thepossibly random -set of sites where at least one split was performed.
Will the mass spread over all of Z d , or will T be a finite subset of Z d ? In the first case, how does it spread? In the last case, what is the size and shape of this set, depending on h and n? How do these answers depend on the order of splitting?
This splitting game is related to abelian sandpile growth models. In fact, the 'splitting pile' rule is the same as the toppling rule in Zhang's sandpile model [10], a sandpile model that has not been considered as a growth model before. A notorious difficulty of moving mass by the Zhang toppling rule is the non-abelianness, that is, the end result depends on the order of topplings [3,6].
In other studies of sandpile growth models [2,4,8], ample use was made of the freedom, by abelianness, to choose some convenient toppling order. In this way, information could be derived about limiting shapes, growth rates and about whether an explosion occurs or not. The term 'explosion' is introduced in [2]: if an explosion occurs, then the mass from the origin spreads over all of Z d , and every site topples infinitely often. For 'robust backgrounds', that is, values of h such that an explosion never occurs, one can examine the growth rate and existence of a limiting shape as n → ∞. Roughly speaking, the growth rate is the radius of T as a function of n, and if the set T , properly scaled, converges to a deterministic shape as n → ∞, then that is the limiting shape.
In [8], the main topic is the rotor router model. In the rotor router model, the mass consists of discrete grains, so that only discrete values of h are possible. All sites except the origin start empty, and at all sites there is a router which points to a neighbor. Instead of splitting, a pile that has at least one grain may give a grain to the neighbor indicated by the router. The router then rotates to point to the next neighbor, for a cyclic ordering of all the neighbors. It was proved that the limiting shape is a sphere, and the growth rate is proportional to n 1/d . The proof makes use of properties of Green's function for simple random walk. In this paper, we will demonstrate that the method of [8] can be adapted for our model, resulting in a growth rate proportional to n 1/d for all h < 0; see Theorem 3.4. In the abelian sandpile model (for background, see [1,9]), the mass also consists of discrete sand grains. Instead of splitting, a pile that consists of at least 2d grains may topple, that is, give one grain to each of its neighbors. Thus, 2d consecutive rotor router moves of one site equal one abelian sandpile toppling of that site. In [4], making use of this equivalence, it was shown that for the abelian sandpile growth model with h → −∞, the limiting shape is a sphere. Moreover, it was proved that for h = 2d − 2, the limiting shape is a cube. In [2], it was proved that for all h ≤ 2d−2, the growth rate is proportional to n 1/d .
For arbitrary splitting order, it may even depend solely on the order if there occurs an explosion or not; see the examples in [6], Section 4.1. In this paper, we focus primarily on the parallel splitting order, but several of our results are valid for arbitrary splitting order.
Another complicating property of our model is that unlike the abelian sandpile model, we have no monotonicity in h nor in n. In the abelian sandpile growth model, it is almost trivially true that for fixed h, T is nondecreasing in n, and for fixed n, T is nondecreasing in h. For our growth model however, this is false, even if we fix the parallel splitting order. One should take care when performing numerical simulations for this model. Since in each splitting, a real number is divided by 2d, one quickly encounters rounding errors due to limited machine accuracy. All simulations presented in this paper were done by performing exact calculations in binary format.
This article is organized as follows: after giving definitions in Section 2, we present our results and some short proofs in Section 3. Our main result is that for h explosive and with the parallel splitting order, the splitting model exhibits several different limiting shapes as t → ∞. We find a square, a diamond and an octagon. These results are stated in our main Theorem 3.1, which is proved in Section 4. Sections 5 and 6 contain the remaining proofs of our other results. Finally, in Section 7 we comment on some open problems for this model.

Definitions
In this section, we formally define the splitting model. While we focus primarily on the parallel order of splitting, some of our results are also valid for more general splitting order. Therefore, we give a general definition of the splitting model, with the splitting automaton (parallel splitting order) as a special case.
For n ∈ [0, ∞) and h ∈ (−∞, 1), η h n is the configuration given by For every t = 0, 1, 2, . . ., and for fixed n and h, η t is the configuration at time t, and the initial configuration is η 0 = η h n . We interpret η t (x) as the mass at site x at time t. We now describe how η t+1 is obtained from η t , for every t. Denote by U t = {x : η t (x) ≥ 1} the set of all unstable sites at time t. S t+1 is a (possibly random) subset of U t . We say that S t+1 is the set of sites that split at time t + 1. Then the configuration at time t + 1 is for all x defined by The splitting order of the model determines how we choose S t+1 , given U t , for every t. For example, if we have the parallel splitting order then we choose S t+1 = U t for every t. In this case, we call our model the splitting automaton. Some of our results are also valid for other splitting orders. In this paper we only consider splitting orders with the following properties: At every time t, S t+1 is non-empty unless U t is empty, and for every x that is unstable at time t, there exists a finite time t 0 such that x ∈ S t+t0 . For example, we allow the random splitting order where S t+1 contains a single element of U t , chosen uniformly at random from all elements of U t . With this splitting order, at each time step a single site splits, randomly chosen from all unstable sites. This splitting order is valid because U t increases slowly enough. Since only the neighbors of sites that split can become unstable, we have for any splitting order that U t+1 ⊆ U t + ∂U t , where with ∂X for a set X ⊂ Z d , we denote the set of sites that are not in X , but have at least one neighbor in X . But when every time step only a single site splits, at most 2d sites can become unstable, so that |U t+1 | ≤ |U t | + 2d.
We are interested in the properties of all the sites that split at least once until time t, as well as For a fixed splitting order, we say that η h n stabilizes if in the limit t → ∞, for every x, the total number of times that site x splits is finite. Note that if η h n does not stabilize, then every site splits infinitely often. We also remark that in order to show that η h n does not stabilize, it suffices to show that T = Z d . Namely, if T = Z d , then every site splits infinitely often. Otherwise, there is a site x and a time t such that x does not split at any time t > t, but each of its neighbors y i , i = 1, 2, . . . , 2d, splits at some time t i ≥ t. But then at time max i t i , x is not stable, because it received at least mass 1 2d from each of its neighbors. Therefore, x must split again.
We call η h n stabilizable if η h n stabilizes almost surely (The "almost surely" refers to randomness in the splitting order). As defined in [2], Definition 2.1. The background h is said to be robust if η h n is stabilizable for all finite n; it is said to be explosive if there is a N h < ∞ such that for all n ≥ N h , η h n is not stabilizable.
Remark We expect the splitting model for every h to be either robust or explosive, independent of the splitting order (see Section 7.1). However, even for a fixed splitting order we cannot a priori exclude intermediate cases where the background is neither robust nor explosive. For example, since the splitting model is not monotone in n, it might occur for some h that for every n, there exist n 1 > n 2 > n such that η h n1 is stabilizable, but η h n2 is not. Finally, we give our definition of a limiting shape. In [2,4,8], limiting shape results were obtained for T in the limit n → ∞, for robust background. In this paper however, we present limiting shape results for T t with n fixed, parallel splitting order and explosive background. We study the limiting behavior in t rather than in n. Accordingly, we have a different definition for the limiting shape.
Let C denote the cube of radius 1/2 centered at the origin {x ∈ R d : max i x i ≤ 1/2}. Then x + C is the same cube centered at x, and by V + C we denote the volume x∈V (x + C). Definition 2.2. Let V t , with t = 0, 1, . . . be a sequence of sets in Z d . Let S be a deterministic shape in R d , scaled such that max xi {x ∈ S} = 1. Let S and S denote respectively the inner and outer -neighborhoods of S. We say that S is the limiting shape of V t , if there is a scaling function f (t), and for all > 0 there is a t such that for all t > t , S ⊆ f (t) (V t + C) ⊆ S . If S is the limiting shape of T t , then we say that S is the limiting shape of the splitting automaton.

Main results
We have observed -see Figure 1 -that varying h has a striking effect on the dynamics of the splitting automaton. For large values of h, T t appears to grow in time with linear speed, resembling a polygon, but which polygon depends on the value of h. Figures 1 and 3 support the conjecture that for the parallel splitting order, as t → ∞, there are many possible different limiting shapes depending on h. Our main result, Theorem 3.1, is that there are at least three different polygonal limiting shapes, for three different intervals of h. For small values of h, the splitting model stabilizes; see Theorem 3.2. In between, there is a third regime that we were not able to characterize. It appears that for h in this regime, T t keeps increasing in time, but does not have a polygonal limiting shape. We comment on this in Section 7.  Let D be the diamond in R d with radius 1 centered at the origin; let Q be the square with radius 1 in R 2 centered at the origin, and let O be the octagon in R 2 with vertices (0,1) and ( 5 6 , 5 6 ), and the other six vertices follow from symmetry. See Figure 2 for these shapes. Theorem 3.1.
(1) The limiting shape of the splitting automaton on Z d , for 1 − 1 2d ≤ h < 1 and n ≥ 1, is the diamond D. The scaling function is f (t) = 1 t . (2) The limiting shape of splitting automaton on Z 2 , for 7/10 ≤ h < 40/57 and The limiting shape of the splitting automaton on Z 2 , for 5/7 ≤ h < 13/18 and n = 3, is the octagon O. The scaling function is f (t) = 5 3t . In the abelian sandpile growth model, h = 2d − 1 is the only possible explosive background value, and our proof of part 1 also works for that situation. However, the second two parts have no analog in the abelian sandpile growth model.
Our proof uses the method of mimicking the splitting automaton with a finite state space cellular automaton; we consider our three explicit results as introductory examples for this method. We expect that with this method, many more limiting shape results can be obtained.
Next, we characterize several regimes of h for the splitting model on Z d .
In the splitting automaton, for d ≥ 2, there exist constants C d < 1 − 3 4d+2 such that the background is explosive if h ≥ C d .
We give the proof of the first part here, because it is a very short argument. The proof of parts 2 and 3 will be given in Section 5, where we give the precise form of C d . We do not believe that this bound is sharp. From the simulations for d = 2 for example, a transition between an explosive and robust regime appears to take place at h = 2/3, while in our proof of part 3, C 2 = 13/19 ≈ 0.684.

Proof of Theorem 3.2, part 1
If n ≥ 1, then the origin splits, so then T is not empty. Now suppose that T is a finite set. Then there exist sites outside T that have a neighbor in T . Such a site received at least 1 2d , but did not split. For h ≥ 1 − 1 2d , this is a contradiction. Note that Theorem 3.2 does not exclude the possibility that there exists, for d fixed, a single critical value of h that separates explosive and robust backgrounds, independent of the splitting order. We only know this in the case d = 1, for which the first two bounds are equal.
We give another result that can be proved by a short argument: Proof. Because n ≥ 1, T is not empty. Suppose that T is not a rectangle. Then, as is not hard to see, there must exist a site that did not split, but has at least two neighbors that split. Therefore, its final mass is at least h + 1 d , but strictly less than 1. This can only be true for h < 1 − 1 d . In the case of the parallel splitting order, we additionally have symmetry. Thus, for h ≥ 1 − 1 d , T is a cube. We remark that the above proof also works for the abelian sandpile model, thus considerably simplifying the proof of Theorem 4.1 (first 2 parts) in [4].
Our final result gives information on the size and shape of T when h < 0. This theorem is similar to Theorem 4.1 of [8]. Let B r denote the Euclidean ball in R d with radius r, and let ω d be the volume of B 1 .

Limiting shapes in the explosive regime
In this section, we will prove Theorem 3.1. Each part of this theorem is stated for h in a certain interval, and the first two parts for n in a certain interval. That means the theorem is stated for uncountably many possible initial configurations. However, we will show that we do not need to know all the exact masses to determine T t for a certain t. For each part of the theorem, we will introduce a labeling of sites, using only a finite number of labels. It will suffice to know the labels of all sites at time t, to determine T t for all t ≥ t.
We will see in each case that the time evolution in terms of the labels is a lot more enlightening than in terms of the full information contained in η t . In each case, we can identify a certain recurrent pattern of the labels. Our limiting shape proofs will be by induction in t, making use of these recurrent patterns.
The label of a site at time t will depend on its own label at time t − 1 plus those of its neighbors at time t − 1. We will specify the labels at t = 0, and the transition rules for the labels. In other words, for each part of the theorem we will introduce a finite state space cellular automaton, that describes the splitting automaton for certain intervals of h and n in terms of the labels.
A cellular automaton is defined by giving its state space S, its initial configuration ξ 0 , and its transition rules. By ξ t (x), we denote the label of site x at time t. The state space will consist of a finite number of labels. The transition rules specify how the label of each cell changes as a function of its own current label and those of its neighbors. A cellular automaton evolves from the initial configuration in discrete time; each time step, all cells are updated in parallel according to the transition rules. We use the following notation for the transition rules. Let s and s, s 1 , s 2 , . . . denote labels in S (not necessarily all different). By we mean that if a cell has label s, and there is a permutation of the labels of its 2d neighbors equal to {s 1 , s 2 , . . . , s 2d }, then the label of this cell changes to s . By * , we will denote an arbitrary label. For example, if we have a transition rule s ⊕ s 1 , s 1 , * , . . . , * −→ s , then the label of a cell with label s will change into s if at least two of its neighbors have label s 1 , irrespective of the labels of the other neighbors.
We first give the proof for part 1, which is the simplest case. In fact, in this case the splitting model is equivalent to (1,d) bootstrap percolation: T t+1 is the union of T t with all sites that have at least one neighbor in T t . The proof we give below will seem somewhat elaborate for such a simple case. That is because we use this case to illustrate our method of labels and cellular automata.
We will need the following observation, which can be proved by induction in t: We call x an odd site if i x i is odd, otherwise we call x an even site. Then in the splitting automaton, even sites only split at even times, and odd sites only at odd times.
Proof of Theorem 3.1, part 1 We begin by defining the diamond cellular automaton. Definition 4.2.
The diamond cellular automaton has state space {e, , u}. We additionally use the symbol s to denote e or . In the initial configuration, every cell has label , only the origin has label u. The transition rules are: The above set of transition rules is sufficient to define the diamond cellular automaton, because, as we will demonstrate below, other combinations of cell and neighborhood labels do not occur.
Let G t , the growth cluster of the cellular automaton, be the set of all cells that do not have label at time t. We will first prove the limiting shape result for G t , and then demonstrate that if 1 − 1 2d ≤ h < 1, then G t is the same set as T t ∪ ∂T t for every t. Let D r be the diamond {x ∈ Z d : i |x i | ≤ r}. To prove the limiting shape result, we will show by induction in t that G t = D t , so that the limiting shape of G t is D, with scaling function f (t) = 1 t .
Our induction hypothesis is that at time t, ξ t is as follows (see figure 4): all sites x ∈ D t have label u if ( i x i − t) mod 2 = 0, and label e otherwise. All other sites have label . If this claim is true for all t, then we have G t = D t . u u e u u e u e u u e u u Figure 4. The induction hypothesis for the diamond cellular automaton at t = 2. Labels not shown are .
As a starting point, we take t = 0. At that time, the origin has label u, and all other cells have label . Therefore, the hypothesis is true at t = 0. Now suppose the hypothesis is true at time t. Then all sites with label e have 2d neighbors with label u, therefore by the second transition rule they will have label u at time t + 1. All sites with label u have 2d neighbors with label e or , therefore by the third rule they will have label e at time t + 1. All sites in D t+1 \ D t have label and at least one neighbor with label u. Therefore, by the fourth rule they will have label u at t + 1. Other labels do not change, by the first rule. This concludes the induction, and moreover shows that our set of transition rules is sufficient to define the diamond cellular automaton.
Finally, we show that if 1− 1 2d ≤ h < 1, then G t is the same set as T t ∪∂T t . To compare the configurations η t and ξ t , we give a mapping For a fixed diamond cellular automaton configuration ξ, define With this mapping, the initial configuration of the splitting automaton η h n is in M ξ0 d , for all n ≥ 1 and 1 − 1 2d ≤ h < 1. We will show by induction in t that η t is in M ξt d for all t. Suppose that at time t, η t is in M ξt d . We check one by one the transition rules: • (rule 1) If in the splitting automaton a site has mass h, and none of its neighbors split, then its mass does not change. This is true for all h < 1. • (rule 2) If in the splitting automaton a stable site has 2d unstable neighbors, then its receives at least 2d times 1 2d , therefore its mass will become at least 1. This is true for all h.
• (rule 3) If in the splitting automaton a site has mass at least 1, then it splits. If none of its neighbors splits, then its mass will become 0. This is true for all h. • (rule 4) If in the splitting automaton a site with mass h has at least one neighbor that splits, it receives at least 1 2d . Therefore, it will become unstable only if . This completes the induction.
Finally, by the following observations: • only the label u maps to mass 1 or larger, so a site is in T t if and only if it has had label u at least once before t, • the label of a site changes into another label if and only if at least one neighbor has label u, • if a site does not have label at time t, then it cannot get label at any time t ≥ t, we can conclude that G t of the diamond cellular automaton is the same set as T t ∪ ∂T t for the splitting automaton with 1 − 1 2d ≤ h < 1 and n ≥ 1. We now give the proofs of the remaining two parts; note that in these next two proofs, we are in dimension 2. We will need more elaborate cellular automata, in which there are several labels for unstable sites. For example, it is important to know whether the mass of a splitting site is below or above 4(1 − h): if its neighbor has mass h then in the first case it might not become unstable, but in the second case, it will.

Proof of Theorem 3.1, part 2
We begin by defining the square cellular automata. Definition 4.3.
The square cellular automaton has state space {e, , p, m, m , c, d}. We additionally use the symbol s to denote a label that is e, or p. In the initial configuration, every cell has label , only the origin has label d. The transition rules are: The above set of transition rules is sufficient to define the cellular automaton, because, as we will demonstrate below, other combinations of cell and neighborhood labels do not occur.
Recall that the growth cluster G t is the set of all cells that do not have label at time t. To prove the limiting shape result for the growth cluster of the square cellular automaton, we use induction. Let C r ∈ Z 2 be the square {(i, j) : |i| ≤ r, |j| ≤ r}. Let ζ r be the following configuration (see Figure 5): • All sites in C r−1 have label c.
• The labels in C r \ C r−1 are d, if (i − j) mod 2 = 0, and e otherwise.
• The labels outside C r are p if they have a neighbor with label e, and otherwise. Our induction hypothesis is that for every even t, ξ t = ζ t/2 . The initial configuration ξ 0 of the square cellular automaton is ζ 0 . Now suppose that at some even time t, ξ t = ζ t/2 . Then by using the transition rules, one can check that at time ξ t+2 will be ζ t/2+1 = ζ (t+2)/2 . This completes the induction. The shape result now follows: for every t, C t/2 ⊆ G t ⊆ C t/2+1 , so that the limiting shape of G t is the square Q, with scaling function f (t) = 2 t . Finally, to show that G t for the square cellular automaton is equal to T t ∪ ∂T t for the splitting automaton with 4 − 4h ≤ n ≤ 16 − 20h and 7/10 ≤ h < 40/57, we give a mapping For all h < 3/4, these intervals are nonempty, moreover, 4 − 4h > 1.
With this mapping, one may check that η h n is in M ξ0 s . By induction in t, we will show that η t is in M ξt s for all t. Suppose at time t, η t is in M ξt s . We check one by one the transition rules: • (rule 1) If in the splitting automaton a site has mass h, and none of its neighbors split, then its mass does not change. This is true for all h < 1. • (rules 2-5) If an unstable site splits, then by Lemma 4.1, its neighbors do not split.
Therefore, it will become empty. If a cell has split at least once, then from that time on it cannot receive sand from more than 4 neighbors before splitting itself. Therefore, no cell that split at least once can gain mass greater than 16 − 20h. This is true for all h < 1.
. This is true for all h < 1.
. This is true for all h < 1. For the final part, we first perform 8 time steps in the splitting automaton before we describe its further evolution as a finite state space cellular automaton. Otherwise, we would need many more labels and transition rules.

Proof of Theorem 3.1, part 3
We begin by defining the octagon cellular automaton. Definition 4.4.
The octagon cellular automaton has state space {e, , p, m, d, d , d!, c, c , q, q }. We additionally use the symbol s to denote a label that is e, or p, and the symbol u to denote a label that is any of the other.
The initial configuration is given in the table below. We show only the first quadrant (left bottom cell is the origin). The rest follows by symmetry. All labels not shown are . Our proof that the growth cluster of the octagon cellular automaton has O as limiting shape, is by induction. We will show that there is a pattern that repeats every 10 time steps. To describe this pattern, we introduce two subconfigurations that we call 'tile' and 'cornerstone'. They are given in Figure 6. We say the subconfiguration is at position (x, y), if the left bottom cell has coordinates (x, y).  To specify a configuration, we will only specify the cells x = (x, y) with y ≥ x ≥ 0; the rest follows by symmetry. • there is a cornerstone at position (5(i + 1), 5(i + 1)), • for every j = 0, . . . , i, there is a tile at position (5(i − j), 7 + 5i + j), • the leftmost tile is different, namely, χ i (0, 7 + 6i) = c , • for every cell (x, y) such that it is not in a tile or cornerstone, but every directed path from (x, y) to (0, 0) intersects a tile or cornerstone, χ i (x, y) = , • for every other cell (x, y), χ i (x, y) = e if (x + y) mod 2 = 0, and c otherwise.
Our induction hypothesis now is: In the octagon cellular automaton, at time 5 + 10i, ξ 5+10i = χ i . In words this says that every 10 time steps, an extra tile is formed.
The hypothesis can be verified for i ≤ 6, by performing 65 time steps starting from the initial configuration ξ 0 . We show the results of this computation in Figure 7, generated by a computer program of the octagon cellular automaton.
Suppose now that the hypothesis is true at time 5 + 10i, with i > 6. We will construct ξ 5+10(i+1) by performing 10 time steps starting from χ i . Observe that the label of a cell at time t + 10 depends only on its own label and that of all cells in x + D 10 ; we call this set of cells the '10-neighborhood' of x. By the definition of χ i , we have that for every i > 5 and for every x, there exists a cell y such that the labeling of the 10-neighborhood of x in χ i is identical to that of the 10-neighborhood of y in χ 5 . Therefore, the label of x in χ i+1 will be identical to that of y in χ 6 . Thus, we can construct ξ 5+10(i+1) from χ i , and we find that indeed, if ξ 5+10i = χ i then ξ 5+10(i+1) = χ i+1 .
Since G t is nondecreasing in t, we have that for every t there is an i such that G 5+10i ⊆ G t ⊆ G 5+10(i+1) . The radius of G 5+10i is 9 + 6i. This means that every 10 time steps, the radius increases by 6. We conclude that the limiting shape of G t is the octagon O, with scaling function f (t) = 5 3t . Finally, we prove that G t of the octagon cellular automaton is equal to T t+8 ∪ ∂T t+8 if n = 3 and h ∈ [5/7, 13/18).
We give a mapping For every x, if n = 3 and h ∈ [5/7, 13/18) then η 8 ∈ M ξ0 o . This can be verified by tedious, but straightforward inspection: In Figure 8 we give the configuration at t = 8 for the splitting automaton with n = 3 and h ∈ [5/7, 13/18). Next, we will prove by induction in t that η t+8 is in M ξt o for all t. We assume that for some t, η t+8 is in M ξt o . By examining every transition rule, we can then show that as a consequence, η t+9 is in M  We check one by one the transition rules: • (rule 1) If in the splitting automaton a site has mass h, and none of its neighbors split, then its mass does not change. This is true for all h < 1. • (rule 2) If an unstable site splits, then by Lemma 4.1, its neighbors do not split.
Therefore, it will become empty.
. This is true for all h < 1.
. This is true for all h < 1.
. This is true for all h < 1. . This completes the induction. Finally, by the following observations: • all sites in G 0 are in T 8 ∪ ∂T 8 , • only labels denoted as u map to values in [1, ∞), so for all t ≥ 8, a site is in T t if and only if it is in T 8 or it has had a label denoted as u at least once before t, • the label of a site with label changes into another label if and only if at least one neighbor has a label denoted as u, • if a site does not have label at time t, then it cannot get label at any time t ≥ t, we can conclude that if n = 3 and h ∈ [5/7, 13/18], then G t for the octagon cellular automaton is the same set as T t+8 ∪ ∂T t+8 .

Explosive and robust regimes
In this section, we prove parts 2 and 3 of Theorem 3.2.

Proof of Theorem 3.2, part 2
We will prove that for all n, all h < 1/2 and all t, |T t | ≤ n 1/2−h , where by |T | we denote the cardinality of a set T ⊂ Z d . It follows that so that for all h < 1/2, the background is robust. Let m 0 be the total mass in T t ∪ ∂T t at time 0, and let m t the total mass in T t ∪ ∂T t at time t. We have At time t, T t ∪ ∂T t contains a total mass of at least 1 2d times the number of internal edges in T t ∪ ∂T t . Namely, consider a pair of sites x and y connected by an internal edge. Each time that one of them splits, a mass of at least 1 2d travels to the other one. The number of internal edges in T t ∪ ∂T t is at least d|T t |. We demonstrate this by the following argument: Fix an ordering for the 2d edges connecting a site to its 2d neighbors, such that the first d edges of the origin are in the same closed orthant. For every site x in T t , all its edges are in T t ∪ ∂T t . If for every x ∈ T t we count only the first d edges, then we count each edge in T t ∪ ∂T t at most once, and we arrive at a total of d|T t |.
Therefore, at least a mass of d 1 2d |T t | = 1 2 |T t | remains in T t ∪ ∂T t . Moreover, since the sites in ∂T t did not split, the mass h at every site in ∂T t also remains in T t ∪ ∂T t . Therefore, we have Finally, we note that since up to time t no mass can have entered or left T t ∪ ∂T t , we have m 0 = m t . Putting everything together, we find from which the result follows.

Proof of Theorem 3.2, part 3
We first give the proof for d ≥ 3.
First, we need some notation. Denote by D r ⊂ Z d the diamond D r = {x : i x i ≤ r}, and by L r ⊂ Z d the layer {x : i x i = r}. Denote d k = (k, k, . . . k) ∈ Z d . Let Γ k,0 = d k , and for i = 1 . . . d, let Γ k,i be the set of sites in L dk+i that have i nearest neighbors in Γ k,i−1 . Observe that for every i, Γ k,i is not empty, and that Γ k,d = d k+1 .
Finally, we define and remark that C d ≤ 1 − 3 4d+2 , with equality only in the case d = 2. We will prove the following statement: Lemma 5.1. In the splitting automaton on Z d , if h ≥ C d , and n ≥ 2d(1 − h), then for every k ≥ 1, at time dk + 2, the sites in Γ k,1 split.
Theorem 3.2, part 3 follows from Lemma 5.1 combined with Theorem 3.3. Lemma 5.1 tells us that for every r, there is a site on the boundary of the cube {x : max i x i ≤ r} that splits at least once. By Theorem 3.3 and by symmetry, all sites in this cube split at least once. Therefore, lim t→∞ T t = Z d .
Proof. Note that at time t, no site outside D t can have split yet, so if a site in Γ k,i splits at time dk + i + 1, it does so for the first time. We will show that in fact the sites in Γ k,i do split at time dk + i + 1. Since we take the parallel splitting order, by symmetry, all sites in Γ k,i distribute the same mass when they split; denote by m(k, i) the mass distributed in the first split of a site in Γ k,i .
We will prove the lemma by induction. For k = 0, the lemma is true, because we chose n large enough. Now suppose it is true for some value k. Then at time dk + 2, the sites in Γ k,2 receive 2 2d m(k, 1) ≥ 1 d , because they each have two neighbors in Γ k,1 . Since they did not split before, their mass is now at least h + 1 d . Therefore, if h ≥ 1 − 1 d , they split at time dk + 3. This condition is fulfilled because C d ≥ 1 − 1 d . Continuing this reasoning, we find for all i = 2, . . . , d − 1 that at time dk + i + 1, the sites in Γ k,i split, because each site in Γ k,i has i neighbors in Γ k,i−1 , so it receives mass i 2d m(k, i − 1) ≥ i 2d ≥ 1 d . We calculate, using that for i = 2, ..., d, we have m(k, Proof. In the proof of Theorem 3.2, part 3, we have proved that at time 2k + 2, sites (k + 1, k) and (k, k + 1) split, making use of the fact that at time 2k + 1, site (k, k) splits. We did not take into account that more sites in L k might split at time 2k + 1.

Remark
We have extended this method further, obtaining even smaller bounds for h, but as we increase the number of sites we consider, the calculations quickly become very elaborate, and the bound we obtain decreases very slowly. The smallest bound we recorded was 0.683. 6. The growth rate for h < 0.
In this section, we prove Theorem 3.4. The proof will follow closely the method used for the proof of Theorem 4.1 in [8], based on the estimates presented in [8], Section 2. An important ingredient to obtain the bounds is that after stabilization, every site has mass at most 1, so that the mass n starting from the origin, must have spread over a minimum number n 1−h of sites. On the other hand, as shown in the proof of Theorem 3.1, part 2, we know: |T | ≤ n 1 2 −h . But crucial is the use of Green's functions and their asymptotic spherical symmetry, allowing to conclude more about the shape of the set of sites that split. Thus we can derive that T contains a ball of cardinality comparable to the coarse estimate n 1−h . Moreover, T is contained in a ball of cardinality close to n 1/2−h . The method for h < 0 does not depend on abelianness or monotonicity, therefore it can be adapted to the splitting model with arbitrary splitting order.
We start with introducing some notation. Denote by P 0 , E 0 as the probability and expectation operator corresponding to the Simple Random Walk X(t) starting from the origin. For d ≥ 3, define Defining the operator ∆ as From [7], we have ∆g(z) = −1 when x is the origin and ∆g(z) = 0 for all other x ∈ Z d . By u(x), we denote the total mass emitted from x during stabilization. Then ∆u(x) is the net increase of mass at site x during stabilization. For all x, let η ∞ (x) be the final mass at site x after stabilization. Since the final mass at each site is strictly less than 1, we have for all x ∈ Z d with δ 0,x = 1 if x is the origin and 0 for all other x. Moreover, since the final mass at each site x ∈ T is in [0, 1), we have for all x ∈ T Proof for the inner bound: For x ∈ Z d , |x| is the Euclidean distance from x to the origin. Let with e 1 = (1, 0, 0, ..., 0). From Lemma 2.2 of [8], we have: Furthermore, using (6.1), we have for all x ∈ Z d , Therefore, u − ξ d is superharmonic, which means that it reaches its minimum value on the boundary. Now, as in the proof of Theorem 4.1 of [8], the estimates in Lemmas 2.1 and 2.3 of [8] can be applied to conclude that there is a suitable constant c 2 such that u(x) is positive for all x ∈ B c1r−c2 . The proof for the outer bound is more involved that that in [8], because Lemma 4.2 from [8], which is valid for the abelian sandpile growth model with h ≤ 0, is not applicable for the splitting model. In essence, this lemma uses that if u(x) > u(y) for some x and y, then the difference must be at least 1, because mass travels in the form of integer grains. Clearly, we have no such lower bound in the splitting model.
We note that in [2], a different proof for the outer bound appeared which is valid for the abelian sandpile growth model with h < d. Unfortunately, we cannot adapt this proof for the splitting model either. We will comment on this in Section 7.
We therefore first present some lemma's which we need to prove the outer bound.
Lemma 6.1. Let h < 0, and take x 0 ∈ T adjacent to ∂T . There is a path Proof. We will first show that we can find a nearest neighbor path such that: Let x 1 be the nearest neighbor of x 0 that loses the maximal amount of mass among all the nearest neighbors of x 0 . If there is a tie, then we make an arbitrary choice. Because x 0 has at least one neighbor that does not split, ∆u(x 0 ) ≤ 2d−1 2d u(x 1 ) − u(x 0 ). Therefore: For k ≥ 1 take x k+1 = x k−1 to be the site that looses the maximal amount of mass among all the nearest neighbors of x k except x k−1 . If there is a tie, then we make an arbitrary choice. It is always possible to choose x k+1 . As long as x k is not the origin, then we get from (6.2) that: Thus we get a chain {x k , k = 0, 1, . . .} of nearest neighbors, that possibly ends at the origin.
We rewrite With this expression, and the fact that u 0 < u 1 by (6.3), it is readily derived by induction that u(x k−1 ) < u(x k ). Inserting this, we obtain u( Now it is left to show that the chain does end at the origin. We derive this by contradiction: suppose the chain does not visit the origin. Then the chain cannot end, because there is always a new nearest neighbor that loses the maximal amount of mass among all the new nearest neighbors. But the chain cannot revisit a site that is already in the chain, because u(x k ) is strictly increasing in k. But by (5.1), the chain cannot visit more than n 1/2−h sites. Therefore, the chain must visit the origin.
Define Q k (x) = {y ∈ Z d : max i |x i − y i | ≤ k} as the cube centered at x with radius k. Let be the average loss of mass of the sites in cube Q k (x), and Proof. From Proposition 5.3 of [6], we know for every x: Equation (6.1) tells that ∆u(y) = η ∞ (y) − h − (n − h)δ 0,y . Therefore Since Q k (x) has 2dk(2k + 1) d−1 internal bounds, we get: where a depends only on k, d and h.
Proof. For x / ∈ T (k) , there is at least one site y 0 ∈ Q k (x) that does not split. For l ≥ 1, take y l as the nearest neighbor of y l−1 that loses the maximal amount of mass among all the neighbors of y l−1 . Since y 0 does not split, we have Therefore u(y 1 ) < 2d(1 − h). For every l > 1, we have from (6.1) that 1 2d y∼y l u(y) < 1−h+u(y l ), therefore u(y l+1 ) < 2d(1−h)+2du(y l ). We know there are at most (2k +1) d sites in {y l } l=0 . Then: Proof of the outer bound: First, we wish to find an upper bound for u(x) for all x with c 1 r − 1 < |x| ≤ c 1 r, that does not depend on n. If x is not in T (k) , then we use Lemma 6.3.
This shows that u (k) − φ d is subharmonic on T (k) . So, it takes its maximal value on the boundary. We combine this information with some lemma's: • Lemma 2.4 of [8] gives that for all x, φ d (x) ≥ −a for some constant a depending only on d. • Lemma 6.3 gives that for every x ∈ ∂T (k) , u(x) < a . • Finally, from Lemma 2.2 of [8], there is ac 2 which only depends on , d and h, such that for x with c 1 r − 1 < |x| ≤ c 1 r, φ d (x) ≤c 2 .
The first two lemma's imply that for x ∈ ∂T (k) , u (k) (x) − φ d (x) ≤ a + a, an upper bound that does not depend on n. Therefore, since u (k) − φ d is subharmonic on T (k) , Combining this with the third lemma, we get: u(x) ≤c 2 + a + a, ∀x ∈ B c 1 r ∩ T k .
Therefore, there is an upper bound for u(x) that does not depend on n, for all x ∈ B c 1 r ∩ T k . From Lemma 6.3, we know also for x / ∈ T k , u(x) < a . Summarizing all, we obtain that for all x with c 1 r − 1 < |x| ≤ c 1 r, u(x) ≤C, withC a constant that does not depend on n.
To summarize, for all x with c 1 r − 1 < |x| ≤ c 1 r, u(x) ≤C, withC a constant that does not depend on n.
Now it remains to show that a site that splits, must lie at a bounded distance c 2 from B c 1 r . This follows from Lemma 6.1: From every site x 0 that splits, there is a path along which u(x) increases by an amount of at least − 2d 2d−1 h every step, and this path continues until the origin. Then along the way, this path must cross the boundary of B c 1 r , and there u(x) ≤C. Therefore, we can choose c 2 = − (2d−1)C 2dh .

Open problems
Based on numerical simulations, we present some tantalizing open problems.

A critical h?
In Theorem 3.2, we give two regimes for h for which we know that the splitting model is explosive resp. robust. In between, there is a large interval for h where we can prove neither. We conjecture however that the two behaviors are separated by a single critical value of h, and that this value does not even depend on the splitting order. In dimension 2, our simulations indicate that this critical h is 2/3. (1) For the splitting model on Z d , there exists a h c = h c (d) such that for all h < h c , the model is robust, and for all h ≥ h c , the model is explosive. (2) h c (2) = 2/3.

The robust regime.
We have proved Theorem 3.4 for all h < 0. We hoped to extend this result to all h < 1/2, by adapting the proof used in Section 3.1 of [2] for the abelian sandpile growth model (ASGM). However, the first step of this proof uses the fact that u n is nondecreasing in n, where u n is the total number of topplings that each site performs in stabilizing η h n . This follows from abelianness of the topplings. Since the splitting model is not abelian, we were not able to adapt this proof to work for our model. Nevertheless, we conjecture  We have only just started classifying the multitude of shapes of the splitting automaton, that one can observe by varying h. We are confident that our method is capable of generating many more limiting shape results. In some cases, we observe that varying n can make a difference, however, we expect the following to be true: Conjecture 7.3. For the splitting automaton on Z d , for every h ∈ [1 − 3 4d+2 , 1) there exists a n 0 such that for every n > n 0 , the limiting shape is a polygon, and depends only on h and d.
This conjecture is reminiscent of Theorem 1 on threshold growth in [5], but the splitting automaton is not equivalent to a two-state cellular automaton.
For smaller values of h, the behavior of the splitting automaton seems to be not nearly as orderly. In Figure 9, we show the behavior at h = 0.667, where we conjecture the model to be explosive. The shape of T seems to alternate between square and rounded. We are not sure whether a limiting shape exists for this value of h.