This project explores the sandpile model in 2D and 3D. Using simulations in the Wolfram Language, it examines how initial configurations, ranging from a single unstable cell at the center to layered, dome-like shapes or difference in the behavior of 2D and 3D sandpiles starting on disk or sphere towers, evolve until they reach stability and the. These types of simulations can be used as a simple analogy for processes such as star formation in gas clouds, where local instability can lead to large-scale changes, even though they are not strictly physically accurate.

Introduction

The Sandpile Model

The sandpile model is a classic cellular automaton known for modeling self-organized criticality, a tendency of some dissipative systems to gravitate toward critical states through cascades of events called avalanches. It helps explain how unpredictable, large-scale behaviors naturally emerge in real-world phenomena.
Consider a grid where grains of sand are dropped, one by one, into one of the cells; if a cell carries as many or more grains of sand as it has neighbors, it is considered unstable, and topples, spreading one grain to each of them. This localized redistribution can subsequently push adjacent cells beyond their own capacities, triggering cascading chain reactions where the size and frequency of the resulting avalanches follow a predictable power-law distribution.

The Relation to Star Formation

In this work, we explore the sandpile model as a simple conceptual analogy for the process of star formation within a gas cloud. Although the model does not capture the full physics of forces such as gravity or fluid dynamics, it is valuable because it illustrates how complex, large-scale patterns can emerge from seemingly simple local interactions.
The direct relationship between the model and star formation is established through the accumulation of matter and local instability, the propagation of energy waves, and simulation using a realistic density gradient.
The paper “Rotational Sandpile Models for Spiral Pattern Formation”, by Zou (2026), establishes a precedent for using the Abelian sandpiles to in cosmological modeling applied to spiral galaxy formation using rotations. In this project, we apply sandpile models to exploring related questions about star formations and supernovae. We explore how the sandpile model propagates particles of matter (grains of sand), given an initial setup analogous to a gas cloud or star in that it is a radially decreasing density gradient with a circular or spherical boundary. The sandpile model distributes matter along this gradient by toppling until it reaches a fixed point.
We represent these density gradients as the sums of binary disk matrices or increasing radii in 2D and 3D. When simulated on these initial conditions, the sandpile model will attempt to distribute matter so as to stabilize the unstable initial states. In some aspects, such as the fact that the model distributes matter outwards to minimize instability, the model is conceptually analogous to the dynamics at play in the formation of a star from a gas cloud, or the dynamics within a star primed to become a supernova.

Circularity and Sphericity

We choose to measure the circularity of sandpile states throughout sandpile model simulations as a proxy for the homogeneity/heterogeneity of the spreading of matter through space given the simulation initial conditions, as this allows for the assessment of the impact of initial conditions on the final symmetry.
To extend our analysis to three dimensions of space, we also compute the sphericity of 3D sandpile states for understanding the propagation dynamics of avalanches. Since the real gas cloud, where stars are born, is dense at the center and thin at the edges, calculating these metrics helps to understand, in quantitative terms, how physical density differences control the propagation of matter and energy through the system.

Sandpile Models on 2D and 3D Arrays

To understand how the sandpile model applies to star formation in gas clouds, first, it is important to define the operation of a sandpile on an array and visualize its progression toward stability.

Sandpile 2D

In the sandpile model on a 2D array, when a cell reaches a total of four grains of sand, it topples, triggering an avalanche effect, and distributes one grain of sand to each of its immediate neighbors. We can capture this behavior in a function using list convolution:
Define a 2D toppling step:
In[]:=
sandStep[s_]:=s+ListConvolve[{{0,1,0},{1,-4,1},{0,1,0}},UnitStep[s-4],2,0]
Applying this function to an array of integers ≥ 0 computes one toppling step of a sandpile model simulation. Lets’s define an initial sandpile to topple.
Initialize a 2D pile:
In[]:=
init=ArrayPad[{{1000}},20];​​frames2D=FixedPointList[sandStep,init];
This makes it possible to visualize the entire progression of the model and how the avalanche unfolds, starting with just a single unstable cell.
Animate 2D toppling:
Out[]=
Ultimately, when the model reaches stability, it is possible to observe the distribution of grains across the grid, as well as the number of sand grains in each of the cells.

Sandpile 3D

The sandpile model can also be represented using a 3D array. In this case, instead of the 4 grains a cell needs to collapse, 6 grains are required, because they are also distributed to the cells directly above and below the one receiving the sand grains.
Define 3D toppling step:
In[]:=
sandStep3D[s_]:=s+ListConvolve[​​{​​{{0,0,0},{0,1,0},{0,0,0}},​​{{0,1,0},{1,-6,1},{0,1,0}},​​{{0,0,0},{0,1,0},{0,0,0}}​​},​​UnitStep[s-6],​​2,​​0​​];
Using this function allows to compute the step-by-step toppling process in a 3D array, and it can be initialized as follows:
Initialize 3D pile:
In[]:=
init3D=ArrayPad[{{{1000}}},5];​​frames3D=FixedPointList[sandStep3D,init3D,2000];
Animate 3D toppling:
Out[]=
As a final result, the model achieves stability, and, just as in the 2D visualization, it is possible to observe the distribution of grains in each cell. Activity unfolds at the center of the sandpile model, and the grains expand toward the outermost layers, increasing the area of the resulting pile.
In the example above, the interior of the sandpile is not visible; however, it is important to be aware of the internal configurations, as the sandpile is not completely filled, and its interior undergoes various changes until all cells become stable. We can view the interior of the 3D sandpile states by shaving off half of the state arrays.
Animate 3D slice view:
Out[]=
Another way to visualize the interior of the sandpile model in 3D is by taking Pac-Man style cutout of the result.
Apply a Pac-Man style cutout to a 3D array:
Visualize Pac-Man style cutout:
In this way, it is possible to visualize the symmetry both inside and outside the model.
After defining the sandpile models and computing the toppling process to observe the results, the impact of the initial configurations can be calculated based on the final visualization of the simulation.

Computing Sandpile State Circularity and Sphericity

Avalanches are a key component in understanding the overall dynamics of matter propagation. Therefore, it is calculated the circularity and sphericity metrics to gain a better understanding of how initial configurations affect the final visualization.

Circularity

Define circularity:
Measure final 2D circularity:
Ultimately, the closer the ratio is to 1, the more circular the 2D sandpile is.
However, this result will not always be so circular; in the initial stages, when starting from just a single unstable cell, the visualization tends to be much more square-shaped due to the shape of the grid cells.
Extract the circularity values from each frame:
Compare the evolution of the circularity:
Thus, one can observe the model’s tendency toward circularity, driven by the effect of avalanches expanding the model’s edges and altering its shape to something more closely resembling a circle.
However, it is not possible to perform this calculation as we have done above for a 3D sandpile model, as the calculation for sphericity differs.

Sphericity

Define sphericity:
Since the final result of the 3D sandpile is composed of multiple blocks, a simpler calculation of sphericity would yield a very low result due to an apparent staircase effect formed on the outer layer.
Measure final 3D sphericity:
As with the 2D model, the closer the result is to 1, the more spherical the 3D sandpile configurations are.
Another similarity to the 2D version is that, just as with circularity, the initial sphericity is very low and tends to become closer to that of a sphere as the sandpile approaches stability.
Extract the sphericity values from each frame:
Compare the evolution of the sphericity:
Thus, the variation in the model’s sphericity is evident, as is its approach to a “perfect” sphere as more grains are added and stability is reached.
In the universe, there are certain events that can be demonstrated using the sandpile model (such as the birth of a star within a gas cloud and the occurrence of a supernova); however, they do not begin with a single particle (or grain of sand), but rather with a matter gradient in which the inner layers are denser and the outer layers are less dense.

Concentric Sandpile Configurations

One way to circumvent the problem of starting with just a single particle is to create an initial configuration for the sandpile model where the center is the most unstable (containing the most grains of sand), and where the layers become less dense (containing fewer grains of sand) as they move outward from the center.

Sum of Concentric Disks Towers with Regularly Spaced Radii

On a 2D grid, this sandpile gradient is represented by towers of concentric disks with regularly spaced radii.
Define disk tower:
It is possible to apply this function to a 2D array and visualize the sandpile states until it becomes stable.
Initialize disk tower:
Inspect disk tower:
This visualization is a 2D analogy for events related to stellar dynamics, such as star formation within a gas cloud and the occurrence of a supernova, since, in both cases, there is a higher concentration of matter in the system’s interior.
However, these configurations are not stable; there is a displacement of matter from the center outward that alters the sandpile’s shape, demonstrating that, in the previously mentioned space events, the systems’ appearance fluctuates over time.
As mentioned earlier, calculating the model’s circularity provides insight into how the sandpile’s initial state shapes its final form.
Circularity of the final sandpile disk tower state:
This calculation allows us to measure the achieved circularity of the sandpile simulation over time, accounting for the physical limitations of the grid. In a regular square simulation grid, the classical toppling of a single sand pile from a single tower does not spread into a perfect circle, but rather into structures with odd, flattened sides and square corners due to the shape of the grid blocks.
Extract the circularity value from each frame:
Compare the evolution of the circularity:
Unlike the evolution of circularity when starting with just a single unstable cell, the model formed by towers of concentric discs with regularly spaced radii begins much more circular but exhibits fluctuations in this rate before reaching stability.
Animate disk tower toppling:
Next, we apply this approach to a sandpile model setup in three dimensions.

Sum of Concentric Sphere Towers with Regularly Spaced Radii

In a 3D grid, this gradient of material density differences is represented by spherical towers with regularly spaced radii.
Define sphere tower:
Using this function on a 3D array, each step the model takes until it ceases to be unstable is computed.
Initialize sphere tower:
Therefore, it is possible to generate a visualization of the sandpile’s initial state, but it is also interesting to visualize the interior to understand the distribution of the layers.
Inspect sphere tower:
In this way, the density differences between the layers more closely approximate the reality of stellar dynamics, because we can think of a sandpile model setup in which the initial condition is a radially decreasing unstable density gradient as an abstract representation of a gas cloud, or the interior of a star.
Following the trend of other ways of representing sandpile models, the final representation is directly affected by the avalanche effects that occur until stability is reached.
Sphericity of the final sandpile sphere tower state:
Thus, the closer the value was to 1, the more spherical the model became, just like the other representations highlighted in this study.
Extract the sphericity value from each frame:
Compare the evolution of the sphericity:
An interesting behavior observed is the lack of changes in the model’s sphericity rate over time, which is always very close to the value of 1.
This occurs because all the activity takes place in the innermost layers of the sphere and is not reflected in the outermost layers, given that the avalanches failed to cause a change in the system’s size before it stabilized.
Animate sphere tower toppling:
Performing sandpile simulations on this setup, and comparing them in 2D and 3D may yield some insight into the way the underlying space of the model affects its dynamics, and the propagation of activity within or outside of the density gradient.
This is represented by the system’s need for a large number of sand grains to generate waves that drive outward activity, analogous to the stellar events studied as a sandpile-like phenomenon.

Concluding Remarks

In this project, we explored the sandpile model as a simple conceptual analogy for stellar dynamics, such as the formation of a star within a gas cloud and the process of a supernova. The application of the model to initial conditions consisting of a radially decreasing density gradient revealed that the system behaves substantially differently in 2D and 3D.
The simulations demonstrated that 2D unstable density gradients, represented by towers of disk matrices, readily result in sand grains spreading outside of the initially specified unstable region. In contrast, when the model is applied in 3D as a stack of balls, the grains stay contained within the initially specified unstable region. The measurement of circularity and sphericity metrics proved fundamental to illustrate these differences and to help understand, in quantitative terms, how physical density differences control the propagation of matter and energy through the system.
This raises an interesting scientific question: To what extent does the dimensionality of space affect the dynamics of forming stars and of supernovae? In future work, the scope will be expanded to attempt to recreate this phenomenon on a large scale. The intention is to see whether the interaction of several such simulations, running simultaneously and rotating, can, for example, simulate a galaxy, thereby helping to further understand spiral pattern formation and the dynamics of the universe.

Acknowledgments

I would like to thank Stephen Wolfram for the research idea and for creating this program. I would also like to thank everyone involved in the WSRI, especially my mentor Phileas Dazeley-Gaist, as well as Júlia Campolim de Oste, Lucas Squillante and all the amazing people I met during this time period.
To my professors at the Federal Rural University of Pernambuco, Jones Albuquerque, Silvana Bocanegra, and Elizabeth Tscha, for all the help that brought me to WSRI.
And last but not least, I would like to thank all my family and friends who helped make it possible for me to take part in this experience.

References

1. Wolfram (2002), A New Kind of Science, Wolfram Media.
2. Bak, Tang, and Wiesenfeld (1988), “Self-organized criticality”, Physical review A.
3. Zou (2026), “Rotational Sandpile Models for Spiral Pattern Formation”, https://www.researchgate.net/publication/406980794_Rotational_Sandpile_Models_for_Spiral_Pattern_Formation

AI Disclosure

The following generative AI tool was used in this project: Gemini 3.1 Pro. It was used for debugging code and editing prose. All code and written content was reviewed, understood and approved by the author.

CITE THIS NOTEBOOK

Sandpile Model Dynamics on Radial Density Gradients​
by Maria Luiza Cordeiro​
Wolfram Community, STAFF PICKS, July 16, 2026
​https://community.wolfram.com/groups/-/m/t/3762924