Chaos in granular matter consisting of non-convex particles

1 Introduction

The behavior of granular materials is fundamentally influenced by particle shape, particularly for non-convex geometries that may interlock (Tran et al., 2024; Aponte et al., 2025). While spherical particles have been extensively studied, the dynamics of non-convex particles remain poorly understood. To phenomenologically capture the behavior of non-convex particles prone to interlocking, dipods and tripods have been introduced recently (Maier et al., 2025). These pods are formed by rigidly connected spheres, abstracting the true shape of real-world particles. Thereby, the need for the (computationally expensive) shape representation of physical particles is avoided. Di- and tripod geometries are characterized by a dimensionless interlocking parameter ξ / D, which governs the transition between free-flowing and interlocking behavior. Understanding how particle shape affects stacking dynamics is critical for optimizing industrial powder handling, predicting granular flow in complex systems, and advancing Discrete Element Method (DEM) simulation strategies for non-convex particulate matter. Our present study focuses on a systematic and fundamental understanding of particle shape effects in quasi-two-dimensional pod systems, with the vision to lay the foundation for more applied studies in the near future.

The influence of particle shape on granular behavior has been extensively documented across multiple scales and systems. Early studies by Olson et al. (2002) demonstrated that grain anisotropy influences ordering, stability, and collective dynamics in granular systems. More recently, Murphy et al. (2019) showed that particle shape affects mesoscale plasticity and stress redistribution during avalanching in granular assemblies. The role of geometric parameters extends to force transmission networks, as Xiao (2023) investigated how particle aspect ratio affects force transmission and pressure dip formation in Three-dimensional (3D) sandpiles. Using DEM simulations, Cúñez et al. (2024) demonstrated that particle shape (i.e., the behavior of spheres versus that of cubes) strongly affects segregation in granular flows, highlighting the role of geometry in collective behavior. Furthermore, Lee et al. (2025) highlighted how particle-scale anisotropy controls bulk stress, packing, and flow properties in sheared granular media, emphasizing the fundamental connection between microscopic geometry and macroscopic material response.

Non-convex particle geometries introduce an additional layer of complexity through “geometric cohesion” and interlocking mechanisms (Tran et al., 2024). Aponte et al. (2024) and Wang et al. (2025b) examined this “geometric cohesion” in Two-dimensional (2D) and 3D systems of star-shaped particles and its effect on collective stability, establishing a framework for understanding shape-induced interactions. This framework was recently expanded by Kim et al. (2025), who identified a percolation transition in entangled networks where rigidity arises purely from geometric interlocking. Experimental validation of these concepts was provided by Aponte et al. (2025), whose study of Platonic polypods showed that non-convex arms induce “geometric cohesion” and solid-like behavior in particle columns. The computational implications of similar - but abstracted - particle geometries were explored by Maier et al. (2025), who demonstrated that interlocking necessitates multiple simulation runs. Crucially, they explicitly addressed the inherent simulation uncertainty, an aspect often overlooked by previous studies.

The complex dynamics of granular systems often exhibit nonlinear and chaotic characteristics that challenge predictive modeling. Hill et al. (1999) showed that segregation-driven patterns in granular flows emerge from the competition between chaos and order. This theme was further developed by Banigan et al. (2013), who identified chaotic dynamics in granular jamming using nonlinear dynamical systems analysis, a connection recently refined by Ness and Cates (2020) who linked the onset of chaos to absorbing-state transitions near jamming. Advanced computational techniques have proven valuable in characterizing these dynamics: Small et al. (2013) applied attractor reconstruction to DEM simulations, revealing chaotic dynamics under large strain compression. More recently, Li et al. (2024) investigated the routes to chaos in coupled granular-structural systems, identifying specific vibration thresholds for chaotic motion. These findings underscore the sensitivity of granular systems to initial conditions and the importance of statistical approaches when analyzing their behavior.

In industrial applications, particle shape significantly impacts flow behavior and system uncertainty. Reddy et al. (2021) demonstrated that dumbbell and disc mixtures alter discharge rates and stagnant zone formation in silos due to geometric interlocking, while Wang et al. (2025a) showed how highly non-convex particles dictate flow regimes and repose angles in rotating drums. Complementing this, Song et al. (2025) investigated scale-up rules for elongated particles in rotating drums, demonstrating that mixing rates and contact forces are governed by the interplay between Froude number and drum size.

Addressing the challenge of calibrating DEM models for complex powder systems, Kronlachner et al. (2023) introduced a block-movement-based calibration approach for fine cohesive powders to capture flow dynamics. Notably, their work addressed variability in continuously operated experiments (and in a quasi steady-state situation) as a calibration feature.

To fully understand the chaotic nature of interlocking and stacking, we examine the variability across multiple realizations of a simple dropping event. Therefore, we will show that recent advances in Machine Learning (Fransen et al., 2025) and white-box discovery of governing equations (Han et al., 2025) offer powerful new tools for data interpretation. Most important, we use the word “chaos” in the sense of the so-called “Sensitive Dependence on Initial Conditions (SDIC)”, which is commonly also known as the “butterfly effect”: infinitesimally small differences in initial particle positions and velocities are amplified by nonlinear contact interactions, resulting in an inability to predict in the final state (Lorenz, 1963; Strogatz, 2015). For the systems studied here, this behavior is expected to be present during settling of our pods, and the resulting distribution of final states is treated as its measurable statistical signature (Lai and Tél, 2011).

Despite progress in understanding granular systems, critical gaps remain: (i) the relationship between the interlocking parameter ξ / D and the transition from free-flowing to interlocking behavior is poorly characterized, (ii) the sensitivity of stacking configurations to initial conditions and the emergence of chaotic behavior requires a more systematic investigation, and (iii) rational guidelines for determining the required number of simulation repetitions as a function of particle shape are absent. These gaps hinder the optimization of DEM simulations for non-convex particles and limit predictive capabilities for industrial powder handling systems.

Leveraging data-driven approaches to analyze the variability across multiple realizations with slightly different initial conditions is the focus of the present study. Our model captures the phenomenological behavior of a certain group of particles; however, it does not aim at representing a realistic particle geometry: the connecting “cylinder” is purely visual. In the simulation, it (i) is empty space that allows other particles to pass through freely, and at the same time (ii) connects spheres rigidly. The specific goals are to: (i) systematically investigate shape effects on stacking behavior of dipods and tripods using DEM simulations, (ii) identify critical shape parameters (ξ / D) that trigger transitions to chaotic behavior, (iii) employ feature engineering and symbolic regression to reveal key dependencies between initial particle configurations and final stacking outcomes, (iv) establish a framework for quantifying stacking behavior using the vertical position of the mean Center of Mass (COM) as a simple metric, and (v) provide rational guidelines for selecting shape parameters and estimating the required number of simulation repetitions in DEM-based studies of non-convex particles.

2 Materials and methods

In the following subsections, we outline the pod approach, the associated simulation setup used during our investigations, and the postprocessing applied to gain our results.

2.1 Pod approach

To investigate interlocking, we focused on parcels, more specifically on so-called “pods”, which are made of a small number of spheres being rigidly connected and equally spaced to the COM of the whole pod. Importantly, the connection to the COM is purely mathematical, thus, other particles can pass through freely. Such pods were introduced by Radl et al. (2024) and previously investigated for their interlocking capabilities in case of tetrapods, i.e., pods in which four spheres in a tetrahedral arrangement are rigidly connected, by Radl et al. (2024), Maier et al. (2025). Pods have another beneficial feature: they can degenerate to a sphere (in case we set ξ / D = 0) while conserving mass by inflating the associated diameter D *. This fact gives rise to the opportunity to transition from a convex particle, i.e., a sphere, to a non-convex pod in a continuous fashion.

2.1.1 Dipods

The simplest pod, a dipod, consists of two spheres rigidly connected by a straight rod. The sphere's distance to the COM of the pod can be described by the pod radius ξ *, while the individual spheres have a diameter D *. Thus, the geometric constant governing their behavior is the interlocking parameter ξ / D. An illustration of dipods is shown in Figure 1A. The associated, purely geometric, critical parameters for interlocking behavior of a dipod are shown in Table 1.

TABLE 1. Critical shape parameters for a dipod, defining transitions between different overlap and interlocking regimes.

2.1.2 Tripods

A tripod extends the pod concept to three spheres rigidly connected to a central point. The three spheres are positioned at equal distances from the COM, described by the pod radius ξ *, and arranged in a planar configuration with 120° angular spacing. Each individual sphere has a diameter D *. Similar to dipods, the interlocking parameter ξ / D governs the geometric behavior and interlocking characteristics of tripods. An illustration of tripods is shown in Figure 1B. The associated, purely geometric, critical shape parameters for interlocking behavior of a tripod are shown in Table 2.

TABLE 2. Critical shape parameters for a tripod, defining transitions between different overlap and interlocking regimes.

2.2 Simulation setup

We performed DEM simulations to investigate the packing behavior of dipods and tripods in a quasi-2D geometry using LIGGGHTS® 3.8.0. The simulation protocol consisted of three main stages: (i) randomly orienting the pods, (ii) pod insertion, and (iii) gravitational settling. All pods were assumed to be almost frictionless to isolate shape effects. The exact procedure is outlined in the following sections. An illustration of the setup is shown in Figure 1. All particle parameters relevant for the simulation are shown in Table 3. The timestep t * was chosen based on 10% of the Rayleigh time step (Thornton and Randall, 1988). The sliding friction coefficient μ r s [pw] was set to 1 × 1 0 − 4, and not exactly 0, to avoid excessively long simulation times, while still allowing for almost frictionless settling. The classical Hertz-Mindlin contact model was used to model particle-particle and particle-wall interactions (Hertz, 1881). Animations for both dipods and tripods using this setup are provided in the Supplementary Material.

TABLE 3. DEM parameters used in the simulations.

2.2.1 Particle insertion and initialization

A total of N pods were inserted at an initial packing fraction of φ 0 = 0.1. Each particle was randomly positioned and randomly rotated in the x and y directions. This low initial packing fraction ensured minimal particle overlap during the insertion phase and allowed for unconstrained initial particle configurations. The limited particle counts (N ≤ 10) were intentional, enabling microscopic-level analysis of particle arrangements. The application to a larger ensemble of particles is explicitly identified as future work, on which our present foundational study can build.

2.2.2 Simulation box geometry

The simulation box was designed to enforce a quasi-2D constraint while allowing realistic particle dynamics. This choice reflects the intrinsic symmetry of our particle geometries (i.e., the pods): since dipods are collinear, and tripods are planar by construction, confining the simulation domain to a quasi-2D region is physically appropriate as no 3D effects originate from the pod shape itself. A fully 3D domain would primarily add out-of-plane degrees of freedom (which was not the focus of our present contribution), and which can be easily studied in follow-up work. Our simulation box geometry was divided into two regions: an insertion box, where particles are introduced, and the remaining settling domain. The insertion box dimensions were defined as follows:

• Depth (z-direction): Equal to the sphere diameter D *, with confining walls positioned at z = ± D * / 2.
• Lateral dimensions (x y-plane): Square cross-section with side length calculated to achieve the target initial packing fraction φ 0 = 0.1. Consequently, the side length varied depending on the number of pods N and their respective volumes. Configurations with large ξ / D and low N resulted in the algorithm's inability to place the pods without overlap. These cases were excluded from further analysis and are left blank in the corresponding figures.
• Height (y-direction): The box center was positioned as such, that it was a side length above the ground, i.e., the bottom wall, to provide sufficient space for gravitational settling.

The bottom wall was an infinite plane located at y = 0 (i.e., “the ground”) to serve as the deposition surface. The narrow depth of the box (equal to one sphere diameter) restricted particle motion to the x y-plane, creating a quasi-2D system while preserving the rotational degree of freedom.

We also investigated a setup with lateral constraining walls, where the dimensions were chosen such that twice the number of inserted pods could fit on the ground. Since this setup did not yield significantly different results compared to the infinite plane configuration, we focus our discussion on the latter. The corresponding figures and tables for the wall-bounded setup are provided in the Supplementary Material.

2.2.3 Settling dynamics

Following insertion and orientation, particles underwent frictionless gravitational settling. The particles were subjected to gravitational acceleration in the negative y-direction, causing them to fall and interact with the bottom wall and neighboring particles. The settling process continued either until E k i n , i * ≤ E pot , 0 * ⋅ 1 0 − 6, where E k i n , i * represents the instantaneous kinetic energy and E pot , 0 * is the initial potential energy, or until 2 × 10⁶ steps had been completed.

2.3 Postprocessing

To extract meaningful physical insights from the simulation data, we implemented a structured postprocessing workflow. This process began with the systematic analysis of outputs generated by the LIGGGHTS® simulations. Following data validation, we defined the center of mass and identified attractors, established sample size requirements, and performed feature engineering. Finally, we employed symbolic regression to derive the governing equations that describe the system’s behavior.

We interpret chaos as SDIC that manifests during the transient settling phase of our pods: repeated simulations with near-identical macroscopic setup but slightly perturbed initial microstates diverge, i.e., their trajectories differ significantly, and may settle into different final configurations. In our present work, the transient settling phase spans from particle insertion through settling until one of the two stopping criteria was reached (i.e., 1: E k i n , i * ≤ E pot , 0 * ⋅ 1 0 − 6, or 2: 2 × 10⁶ steps). Final states are the stable configurations thereafter.

2.3.1 Center of mass (COM) calculation & attractor definition

The initial step in processing the LIGGGHTS® simulation data involved enforcing three key validation criteria: (i) pods actually were in contact with the bottom wall (in some simulations, the simulation terminated before this was actually the case), (ii) sphere-sphere overlap was constrained to be less than 1% in the initial configurations, and 5% in the final configurations, and (iii) the kinetic energy calculated from the vertical velocity component (v y *) was less than 1 0 − 4 of the initial potential energy (this was useful to validate the small number of simulations that reached the second stopping criterion).

For each validated configuration, we computed the vertical coordinate of the system’s center of mass, y COM. The mean COM values of the final configurations were then interpreted as attractors, with position tolerances set to 1.00 × 10⁻⁵ m for dipods and 1.00 × 10⁻⁴ m for tripods. To facilitate comparative analysis, these coordinates were normalized relative to the attractor corresponding to the minimum height configuration (i.e., “attractor 0”), wherein all pods rest individually in their most stable ground positions.

Given that the objective of this work is to characterize the macroscopic outcome of a practical packing process, the vertical COM position (y COM) is used as a physically meaningful scalar parameter to quantify the packing height. For the quasi-2D geometry considered here, the vertical COM position is a sensitive macroscopic indicator of interlocking and stacking, and distinct packing modes are reflected by qualitatively different y COM distributions. Although y COM does not fully resolve the underlying contact network, it provides a simple descriptor for transitions between qualitatively different final packing outcomes.

2.3.2 Sample size considerations

To ensure statistical reliability of the mean COM measurements, we determined the required number of simulation repetitions using confidence interval analysis following the approach of Maier et al. (2025). The required sample size n s t a t for a confidence interval at a two-sided confidence level of 1 − α with width x ̄ ± ( x ̄ ⋅ E ) can be calculated using Equation 1.

n s t a t = t 1 − α 2 , n − 1 ⋅ s E ⋅ x ̄ 2 (1)

Within this definition, t 1 − α 2 , n − 1 is the critical value from the student-t-distribution, s is the sample standard deviation, x ̄ is the sample mean, and E is the relative error. The value for t reaches a constant value for large sample counts ( n ≥ 30 ). Although the underlying distribution of individual COM measurements may exhibit non-normal characteristics such as skewness or multi-modality, the Central Limit Theorem (CLT) ensures that the sampling distribution of the mean approximates a normal distribution given our large sample size. This justifies our use of this confidence interval estimation method (Montgomery and Runger, 2018). In our study, each simulation was typically run 1000 times (i.e., n = 1000) with different seed values for the random number generator that defines the initial state.

For all considerations in the present work, we targeted a relative error E of 5% with a confidence level of 95% (i.e., α = 0.05, E = 0.05), yielding t 0.975 , 999 = 1.96. This ensured that with 95% confidence, the real mean μ of the population with regard to the target variable (in this case y C O M ̄) is within x ̄ ± 5 %.

In addition to the sample size, the probability of classes of attractors was investigated. For this purpose, the attractors were normalized using the “attractor 0” and classified as follows: (i) baseline y COM < 1.01; (ii) low-rise 1.01 ≤ y COM < 1.05; (iii) mid-rise 1.05 ≤ y COM < 1.15; (iv) high-rise 1.15 ≤ y COM. Following this, the 95% confidence interval for each class was computed. This allows for a qualitative assessment of the behavior in case of a multimodal distribution, thus, offering additional information besides just a number of samples needed to match the mean. The data were analyzed up to ξ / D = 2.0 as the statistical behavior stabilized, and no additional information was expected to be gained by extending the analysis to higher ξ / D values.

2.3.3 Feature engineering

To prepare the dataset for subsequent analysis, we first filtered the attractors, retaining only those with at least 10 occurrences and sources associated with more than one distinct label. The dataset was then partitioned stratified by global attractor label, allocating 0.8 of the data to the training set using deterministic sampling via a fixed seed.

It is important to note that the comprehensive feature set was computed exclusively based on the initial configuration to enable predictive modeling. Each feature was documented by type (continuous or discrete/counts) and fallback value. These features were organized into six distinct categories: (i) statistical and aggregate measures, including basic statistics and per-molecule aggregates projected onto the XY-plane (XY); (ii) spatial and distance-based metrics, such as pairwise distances, nearest-neighbor statistics, and the Radial distribution function (RDF); (iii) structural descriptors, encompassing local structure indices, Bond Orientational Order (BOO) parameters, cluster analysis metrics, and shape descriptors; (iv) topological and graph-based features, including graph connectivity measures, Topological Data Analysis (TDA) proxies, and interaction fields; (v) geometric properties, such as moments of inertia in the XY, symmetry descriptors, and angle-based features; as well as (vi) multi-scale features, incorporating multi-resolution analysis and fixed-size aggregate context.

To prevent data leakage, normalization parameters were derived exclusively from the training set. Three distinct normalization strategies were applied: (i) Z-score normalization for statistical measures and dimensionless quantities; (ii) prior dimensional scaling followed by Z-score normalization for physical quantities requiring dimensional analysis; and (iii) no normalization for raw counts, identifiers, categorical features, or features yielding null values after dimensional scaling.

Feature selection was performed solely on the training dataset through a multi-stage procedure. Initially, all constant features and those containing Not a Number (NaN) values were removed. Next, features were ranked according to their relevance to the target variable using sigmoid-modulated conditional mutual information maximization, employing a greedy Conditional Mutual Information Maximization (CMIM) selection algorithm for regression. This approach assigns a discounted redundancy penalty (minimum weight ≈ 0.3) to redundant but highly informative features, rather than a fixed penalty (Fleuret, 2004; Vergara and Estévez, 2010; Souza et al., 2022; Gong et al., 2022). The feature set was further refined using a Quantum-behaved Parallel Swarm Ant Lion Optimization (QPSALO)-style wrapper method (Zhou et al., 2022). The final feature subset identified during training was then consistently applied to the test dataset to maintain methodological rigor.

2.3.4 Symbolic regression

Building on the curated feature set, we employed symbolic regression to distill interpretable, physically meaningful governing equations from the data. Unlike black-box models such as neural networks, symbolic regression offers transparency and facilitates the discovery of natural laws (Schmidt and Lipson, 2009; Udrescu and Tegmark, 2020; Han et al., 2025).

The symbolic regression process was implemented in Julia, using 40 iterations, a population size of 1000, and a tournament selection size of 5 (Cranmer, 2023; Bezanson et al., 2017). The top-performing equation was evaluated on the test dataset in Python, with performance metrics reported as Coefficient of Determination (R ²) and Root Mean Square Error (RMSE).

To ensure the model’s physical interpretability and guard against overfitting, we analyzed the derived equation using partial dependence plots (PDPs) and Shapley Additive Explanations (SHAP) analysis (Lundberg and Lee, 2017; Friedman, 2001). The Partial Dependence Plots (PDPs) revealed the influence of individual features on the predicted value, indicating whether relationships were direct, indirect, or non-monotonic. Meanwhile, the SHAP analysis highlighted the formation of feature clusters, offering insights into the non-monotonic behavior of certain variables.

3 Results

This section presents the main findings of our study, structured in two parts: first, we examine the behavior of dipod systems, followed by an analysis of tripod configurations. Each analysis is structured by the number of the involved pods N.

3.1 Dipods

3.1.1 Distribution of dipods

The distribution of attractors for dipod systems, as illustrated in Figure 2, exhibits a clear dependence on the geometric parameter ξ / D. As ξ / D increases, the system transitions through distinct regimes, revealing both universal behaviors and configurations specific to the number of dipods present. An overview of all distribution data gained during this study for dipods is shown in the Supplementary Figure S1.

3.1.1.1 Two-Dipod systems

For systems comprising two dipods, the attractor distribution is characterized by a right-skewed, mono-modal profile, with only two dominant attractor types observed (excluding “attractor 0”). This skewness arises because configurations within the same attractor type progressively flatten with increasing ξ / D. Notably, stacking behavior emerges abruptly at ξ / D = 0.8, which also holds true for larger systems up to eight dipods. This threshold aligns with the geometric requirement for interlocking, which theoretically begins at ξ / D = 1 / 2 (see Table 1). The discrepancy between these two values (the observed and theoretical thresholds) is a sampling resolution artifact, as we chose the sampling increment equal to 0.1 in the ξ / D space. Stacking ceases at ξ / D = 1, where the gap between spheres becomes sufficiently large to allow one sphere to pass through. For small systems (i.e., N = 2, and N = 3 as discussed in the next paragraph), stacking ceases at ξ / D = 1, where the inter-sphere gap becomes large enough for one sphere to pass through freely, precluding stable interlocked configurations. For larger systems, however, this upper bound does not mark a simple cessation of stacking. Instead, the additional geometric degrees of freedom introduce further complexity, and new attractor types emerge beyond this point.

3.1.1.2 Three-Dipod systems

The introduction of a third dipod transforms the distribution from mono-modal to bimodal, with two additional distinct attractor types emerging at ξ / D values of 0.8 and 1.2. This shift underscores the limited diversity of attractor types relative to the total number of attractors (Table 4), suggesting that geometric constraints predominantly govern the accessible configuration space.

TABLE 4. Total number of identified attractors as a function of the number of pods ( N ) for dipod and tripod systems.

3.1.1.3 Systems with four or more dipods

In systems with more than three dipods, the attractor distribution becomes multimodal. A key observation is the gradual increase in the mean center of mass (COM) height with the addition of more dipods (Figure 3A). The additional transition points at ξ / D ≈ 1.0, 1.4, and 1.9 (emerging with increasing N and clearly visible in Table 5) are consistent with new collectively interlocking configurations that only become geometrically accessible once sufficient particles are present to form them. With N = 2 only the simplest interlocking mode is possible; larger N expands the configuration space, introducing additional attractor states at characteristic ξ / D values. However, a precise geometric derivation of each threshold is beyond the scope of this present study and remains an open question.

TABLE 5. Probability of “attractor 0” for all investigated scenarios (i.e., ξ / D, N combinations) for dipods. Cell shading intensity corresponds to the probability value.

3.1.1.4 Friction effect for fixed ξ / D values

Having established how system size shapes the frictionless attractor distribution, we now examine the role of friction at fixed ξ / D. Specifically, at a fixed value of ξ / D = 0.8, our data for simulations using friction (i.e., the “with-friction panels” (iii) and (iv) in Figure 4) show a larger number of discrete y COM outcomes than the “frictionless panels” (i) and (ii). This is true for N = 4 and N = 8. However, the overall spread of the outcome (i.e., y COM values) remains comparable. This is consistent with a qualitatively stronger SDIC when considering friction effects: we emphasize that Figure 4 illustrates only a qualitative effect of friction at this particular value for ξ / D, and hence should not be interpreted as a general claim that friction always increases the number of discrete y COM outcomes. Interestingly, the maximum of the y COM values for both setups (i.e., almost frictionless and with friction) are similar. Since we have performed only a limited number of higher-friction simulations, we treat ξ / D ≈ 0.8 as a specific crossover value for (almost) frictionless dipods, rather than a universal boundary for pods with any value of friction.

3.1.1.5 Probability of “attractor 0” for dipod systems

The geometric boundaries of dipod systems also significantly influence the probability of “attractor 0” (Table 5). Each time an additional sphere fits between dipods, the probability of “attractor 0” spikes, whereas the introduction of new dipods reduces its likelihood at the same ξ / D. Conversely, the probability of non-stacking is lowest just before a new sphere can fit, highlighting the system’s sensitivity to geometric constraints.

3.1.2 Sample size considerations

Following the concept outlined in Section 2.3.2, we calculated the required sample size for all investigated pairs of (ξ / D, N), presenting the results in Table 6. Given the non-normal underlying data distribution, we validated the use of t-statistics by bootstrapping the mean 1000 times to approximate its sampling distribution and confirming normality via the Shapiro-Wilk test (Efron and Tibshirani, 1993; Shapiro and Wilk, 1965). Pairs of (ξ / D, N) not meeting this requirement are marked with *. Such non-conforming pairs do not appear systematically across ξ / D and N, but are scattered across the whole data range (see Table 6). Thus, such cases were treated as stochastic sampling artifacts; since they do not cluster in any particular region of the parameter space, the general validity of the CLT still holds for our dataset. The calculated sample sizes reveal a clear bifurcation into a deterministic regime (for ξ / D < 0.8) and a chaotic regime where the variance becomes significant. The increased spread and multimodality of the final COM values in this region is consistent with SDIC-driven divergence, i.e., the famous “butterfly effect” that is used to describe chaotic behavior of a system. Also, this behavior is very different from chaotic fluid flow, e.g., turbulence, since we do reach a non-trivial steady state position of our pods for long times. While variability generally scales with ξ / D, the required sample size exhibits a complex, non-monotonic dependence on particle count within the chaotic regime. This analysis reveals that choosing a sample size of 100 is adequate for preliminary investigations, providing a balance between computational efficiency and statistical reliability. However, for high-precision studies, especially in the chaotic regime, larger sample sizes are advisable to capture the full extent of variability, including rare events.

TABLE 6. Required sample size for all investigated scenarios involving dipods. Cell shading intensity is proportional to the sample size. Values marked with * indicate that the CLT conditions were not met. Entries marked with - indicate that the standard deviation was zero, preventing the calculation of the sample size.

Looking at the classified dipod data shown in Table 7, derived using the method outlined in Section 2.3.2, a few key observations can be made: (i) at ξ / D = 0.7, the structure collapses to baseline for all N (P B ≈ 1); (ii) as ξ / D and N increase, the probability shifts away from the baseline toward high-rise, with the transition already visible at ξ / D = 0.8 and becoming pronounced at ξ / D = 0.9; (iii) for larger systems (approximately N ≥ 7), high-rise dominates most cases for ξ / D ≥ 1.2, while mid-rise remains a persistent intermediate regime. Some low-N cases show localized non-monotonic re-entrance toward baseline; and (iv) low-rise outcomes are rare throughout and never dominant.

TABLE 7. Full classified outcome probabilities for dipod cases up to ξ / D = 2.0.

Statistical validation confirms that the generated dataset of 1000 samples comfortably exceeds the maximum calculated requirement of 103 samples across all configurations, while ensuring that the confidence intervals of ±0.02–0.03 are sufficiently narrow to reliably distinguish between outcome classes throughout the entire parameter space.

3.1.3 Symbolic regression for attractor basin identification of two dipod systems

To identify the basins of attractors for two-dipod systems, we employed the feature engineering, selection, and symbolic regression framework described in Sections 2.3.3 and 2.3.4. This process yielded a closed-form equation (Equation 2) for two-dipod systems, utilizing the variables defined in Table 8. The derived equation achieved an R ² of 0.839 and an RMSE of 0.125.

y C O M = ξ / D E + 0.997 E = 5.676 − cos s t d _ x _ p e r _ m o l _ m a x F F = 1.869 − r i p l e y _ L r _ m i n u s _ r _ m e a n ⋅ 1.100 − v o r o n o i _ p e r i m e t e r _ m e a n d i s t a n c e s _ m e a n − 4.029 (2)

TABLE 8. Selected features used for symbolic regression analysis of two-dipod systems, along with their physical and geometric interpretations.

3.1.3.1 Partial dependence analysis

The partial dependence plots (PDPs) presented in Figure 5 reveal the relative effect of each feature on the predicted center of mass (y COM). Both the geometric parameter ξ / D and the spatial metric ripley_Lr_minus_r_mean exhibit strongly increasing, convex relationships with the predicted COM (train–test correlation = 1.000). These results confirm that higher spatial regularity and elongated particle geometries consistently promote stacking behavior. Conversely, voronoi_perimeter_mean displays a decreasing trend (correlation = 0.978), indicating that smaller Voronoi perimeters—reflective of dense clustering in the 2D projection—are associated with elevated COM values. Of note, std_x_per_mol_max shows a strongly increasing trend in the training set but a decreasing trend in the test set, suggesting potential instability and overfitting for this particular feature.

3.1.3.2 SHAP feature importance

The SHAP analysis (Figure 6A) quantifies the global importance of the features identified. voronoi_perimeter_mean emerges as the dominant predictor, followed by ripley_Lr_minus_r_mean and ξ / D. This hierarchy suggests that local packing structure (captured by Voronoi and Ripley metrics) is the primary determinant of the attractor basin, surpassing the influence of pure geometric constraints. The SHAP summary plot (Figure 7A) further reveals distinct clusters of SHAP values, reflecting the multimodal nature of the system. These clusters correspond to discrete regimes of packing density, each associated with a specific attractor basin.

3.1.3.3 Model performance and implications

Collectively, this analysis underscores that transitions between attractors are physically driven by local packing density. The model effectively leverages the Voronoi perimeter as a proxy for stacking: stacked dipods overlap in the 2D plane, generating small, densely packed Voronoi cells. The derived equation successfully assigns all samples to the attractor type with a higher mean COM, as seen in Figure 8A and its corresponding Table 9A, demonstrating the transferability of identified trends. However, the error is primarily driven by two rare intermediate attractors (1.388 and 1.477), representing only ≈ 11 % of samples but exhibiting significantly higher RMSE values. This indicates that while the model captures the main stacking regimes, it struggles to resolve these transitional configurations. As will be introduced in Section 3.2.3, this behavior stems from the complexity inherent to our data set, and is not directly related to model choice. Especially for low N, the model error is driven by rare attractors, while the general trend is well captured.

TABLE 9. Root Mean Square Error (RMSE) of the symbolic regression model for (A) two-dipod and (B) two-tripod systems, broken down by attractor basin.

3.2 Tripods

3.2.1 Distribution of tripods

When transitioning from dipods to tripods, the attractor distribution retains its multimodal character, but shifts it to higher values of ξ / D, as illustrated in Figure 9. This shift arises from the distinct geometric boundaries inherent to tripods compared to their dipod counterparts. While tripod systems do not reach the same stacking height as dipod systems, they exhibit a greater probability for stacking. The number of individual attractors differs markedly (Table 4). The overall shallower distribution with a tenfold increase in geometric tolerance, as stated in Section 2.3.1, explains the gradual growth in attractor count for systems comprising more than five tripods. An overview of all distribution data gained during this study for tripods is shown in the Supplementary Figure S2.

3.2.1.1 Two-Tripod system

In tripod systems, stacking initiates at a ξ / D value of 1.0, at which point one sphere fits precisely within the tripod cavity (Table 2). Additionally, the cavity constrains the spheres, resulting in greater similarity between initial and final states, as trapped spheres cannot escape. For two-tripod systems, four distinct attractor types are observed. A notable departure from two-dipod systems is the emergence of sphere-like attractors, as shown in Figure 9 (subplot for two-tripod systems). The distribution remains right-skewed and although data for ξ / D values above 1.2 are not available due to the limitations of our approach outlined in Section 2.2.2, it is hypothesized that sphere-like attractors persist at higher ξ / D values, given their theoretical feasibility.

3.2.1.2 Three-Tripod system

The addition of a third tripod transforms the distribution into a bimodal profile, analogous to that observed in three-dipod systems. A novel configuration appearing at this system is the tube-like attractor as depicted in Figure 9 (subplot for three tripods).

3.2.1.3 Mean COM of tripod systems

An analysis of the mean COM for each combination of tripod number and ξ / D value (Figure 3B) reveals that both multimodality and stacking height are far less pronounced than in dipod systems. The maximum mean COM for tripod systems (1.20) is less than half the maximum mean COM observed for dipod systems (1.55) under equivalent conditions.

3.2.1.4 Probability of “attractor 0” for tripod systems

The geometric boundaries of tripods result in a generally lower probability of “attractor 0” across the system, and a shifted increase in “attractor 0” to greater ξ / D values compared to their dipod-system counterparts, as shown in Table 10. At a ξ / D value of about 1.15, the gap between spheres widens sufficiently to accommodate one sphere; at ξ / D = 1.8, two spheres fit between two spheres of the tripod. These geometric thresholds induce a slight increase in the probability of “attractor 0”, though it remains markedly lower than in dipod systems.

TABLE 10. Probability of “attractor 0” for all investigated scenarios (i.e., ξ / D, N combinations) for tripods. Cell shading intensity corresponds to the probability value.

3.2.2 Sample size considerations

Similar to the dipod analysis, we determined the required sample size for tripods across all (ξ / D, N) configurations, as detailed in Table 11. The statistical validation followed the same bootstrapping procedure to ensure the applicability of t-statistics. The calculated sample sizes indicate a transition to a chaotic regime at a higher geometric parameter value compared to dipods, specifically around ξ / D = 1.0. Below this threshold, the system behaves largely deterministically. Notably, the maximum required sample size for tripods is 38, which is significantly lower than the 103 samples required for dipods. This reduction in variability suggests that the additional sphere in the tripod structure imposes stronger geometric constraints, limiting the diversity of stable configurations and thus the variance in the center of mass. Consequently, a sample size of 100 remains more than sufficient for high-precision studies of tripods, providing an even greater safety margin than for dipods.

TABLE 11. Required sample size for all investigated scenarios involving tripods. Cell shading intensity is proportional to the sample size. Values marked with * indicate that the CLT conditions were not met. Entries marked with - indicate that the standard deviation was zero, preventing the calculation of the sample size.

Looking at the classified tripod data using the method outlined in Section 2.3.2, a few key observations can be made (see Table 12): (i) Below the critical ξ / D of approximately 1.0, the structure collapses to baseline for all N; (ii) At and above ξ / D ≈ 1.0, increasing ξ / D and N shifts outcomes from baseline toward mixed mid-rise/high-rise states, with high-rise becoming dominant at larger N and higher ξ / D (reaching about 0.62 at N = 10, ξ / D = 2.0); (iii) A pronounced re-entrant baseline pocket around ξ / D ≈ 1.2 for lower N (e.g., N = 3–5) interrupts otherwise increasing high-rise tendency and suggests a geometric resonance effect; and (iv) Mid-rise persists as a robust intermediate regime across most post-critical conditions, while low-rise appears at moderate but secondary levels and is never dominant.

TABLE 12. Full classified outcome probabilities for tripod cases up to ξ / D = 2.0.

Statistical validation confirms that the generated dataset of 1000 samples comfortably exceeds the maximum calculated requirement of 38 samples across all configurations, ensuring that the confidence intervals of ±0.02–0.03 are sufficiently narrow to reliably distinguish between outcome classes throughout the entire parameter space.

3.2.3 Symbolic regression for attractor basin identification of two tripod systems

To identify the attractor basins for two-tripod systems, we applied the same symbolic regression framework as before. This analysis yielded the closed-form equation shown in Equation 3. This equation utilizes the variables defined in Table 13 and achieved an R ² of 0.818 with an RMSE of 0.073. While the R ² is approximately 2.5% lower than in the two-dipod case, the RMSE is roughly 40% lower. This substantial improvement in error metrics aligns with the significantly narrower distribution of the COM observed for tripods, as discussed in Section 3.2.2. Notably, the model for tripods is more complex, incorporating eight features compared to the five used for the dipod system.

y C O M = log G H G = a v g _ c l u s t e r i n g ξ / D I + 0.964 I = K 0.435 + 0.588 p o i n t _ g r o u p _ s y m m e t r y H = 0.493 J J = K + L ⋅ i n e r t i a _ t e n s o r _ x x K = i n e r t i a _ t e n s o r _ x y ⋅ s k e w n e s s _ x L = l a r g e s t _ c l u s t e r _ c i r c u l a r i t y ⋅ n n _ m e d i a n (3)

TABLE 13. Selected features used for symbolic regression analysis of two-tripod systems, along with their physical and geometric interpretations.

3.2.3.1 Partial dependence analysis

The partial dependence plots (Figure 10) reveal three distinct behavioral categories among the features in the two-tripod system. First, avg_clustering, and point_group_symmetry both exhibit positive correlation with the predicted center of mass. avg_clustering shows a strong, slightly convex increasing trend, indicating that higher clustering coefficients are associated with greater stacking heights. point_group_symmetry, by contrast, shows a weaker, positive linear trend.

Second, inertia_tensor_xx, nn_median, and ξ / D display negative correlations with the predicted COM. inertia_tensor_xx follows a convex decreasing curve, suggesting that its negative impact diminishes at higher values. Both ξ / D and nn_median exhibit approximately linear decreasing trends, confirming that larger geometric spacing and nearest-neighbor distances generally reduce stacking height.

Third, largest_cluster_circularity, inertia_tensor_xy, and skewness_x exhibit non-monotonic behaviors. Notably, largest_cluster_circularity shows a distinct peak, indicating an optimal cluster shape that maximizes stacking height, which is consistent with the sphere-like attractors observable for tripod systems. The complex shapes of inertia_tensor_xy and skewness_x likely capture specific orientational configurations. The close alignment of training and test curves across most features suggests the model robustness, although confidence intervals widen at the extremes of avg_clustering and largest_cluster_circularity.

3.2.3.2 SHAP feature importance

The SHAP analysis (Figure 6B) identifies avg_clustering as the most influential predictor, followed by largest_cluster_circularity and inertia_tensor_xx. This hierarchy indicates that, for tripods, the degree of clustering and the shape of the clusters are more critical determinants of the attractor basin than the pure geometric parameter ξ / D, which ranks fourth in importance.

The SHAP summary plot (Figure 7B) reveals a striking separation in the SHAP values for avg_clustering, where samples split into two distinct clusters with a significant gap. This indicates that avg_clustering acts almost as a binary switch, strongly differentiating between two primary regimes, likely corresponding to stacked and unstacked configurations. largest_cluster_circularity also shows clustering but with more internal variance, while features like inertia_tensor_xx and ξ / D exhibit more continuous distributions of SHAP values. This suggests that these latter features provide fine-tuning within regimes rather than defining the primary state.

3.2.3.3 Model performance and implications

The analysis implies that the transition between attractors in tripod systems is fundamentally distinct from that in dipod systems. While dipods rely on continuous local packing density (Voronoi), tripods are governed by a discrete interlocking mechanism, captured by the binary nature of avg_clustering. The derived equation effectively functions as a robust classifier between non-stacked states (i.e., “attractor 0”) and stacked configurations, driven by this dominant feature, which a marked difference from the two-dipod equation, where it struggled to distinguish exactly in this region, whereas the two-tripod equation conversely struggles with the highest COM configurations. Specifically, Figure 8B and Table 9B highlight the attractors at 1.463 and 1.497 exhibiting elevated RMSE values of 0.213 and 0.141, respectively. Although these complex interlocking states comprise only ≈ 11.5 % of the dataset, they challenge the model’s precision. This limitation likely arises because once tripods interlock, their clustering values converge; thus, secondary geometric features (like inertia_tensor_xx) provide incremental refinement but are insufficient to fully resolve the subtle orientational differences between these distinct high-stacking modes. To investigate the limits of predictability, we compared symbolic regression against cross-validated baseline models. Specifically, we considered linear, tree-based, and random forest-based regression models. Predictability decreases sharply for systems with more than two tripods, regardless of the model, and the same holds for dipod systems. Additionally, models trained on pooled data across all tripod counts performed comparably to those trained on individual counts. Results and hyper-parameters are provided in Supplementary Figure S13 and Supplementary Table S9. Across all baselines, performance is adequate for simple configurations (i.e., low N) but degrades systematically as the attractor distribution widens, suggesting complexity grows sharply with particle count.

4 Discussion

Our study established a unified framework for probing the effect of non-convex geometries, linking their geometric properties (Tables 1,2) to a certain mechanical behavior. By quantifying the results in terms of their COM distribution, we have demonstrated that the interlocking parameter ξ / D serves as a critical predictor for both dipod and tripod systems. This parameter governs the transition from an effectively deterministic settling to a chaotic behavior.

4.1 Geometric control and topological predictors

Our analysis identifies ξ / D not merely as a shape parameter, but as a control switch for the transition in the flow behavior. The correspondence between the interlocking threshold in our setup ( ξ / D ≈ 0.8 ) and the onset of chaos confirms that geometric exclusion drives the shift from deterministic to multi-stable packing, directly dictating the computational cost required to capture system variability. Specifically, we observed that for dipods, significant variability emerges beyond ξ / D = 0.8, necessitating up to ≈ 100 simulations in the most chaotic regimes. In contrast, tripods exhibit a delayed transition to chaotic behavior at ξ / D = 1.0 and require significantly fewer samples (i.e., in our case 38), likely due to increased geometric constraints. As no formal chaos diagnostics (e.g., Lyapunov exponents) were computed in our present contribution, we propose that SDIC (i.e., the “butterfly effect”) is the most probable mechanism that causes the chaotic behavior of our system, a claim that is supported by our ensemble statistics (Ott, 2002). Consequently, we recommend a baseline of 100 simulation repetitions for such non-convex particle systems to robustly capture the underlying variability.

Furthermore, our symbolic regression analysis highlights the superiority of advanced topological descriptors over simple geometric features for predicting final configurations. The dominance of the feature “voronoi_perimeter_mean” in dipod systems reveals that the three-dimensional stacking behavior is mechanically rooted in a local two-dimensional density parameter: the system “stacks” by minimizing the projected void space between centroids. This implies that the flow of complex granular matter is better captured by its statistical topology (i.e., Voronoi- or Ripley-based metrics) than by macroscopic parameters such as the packing fraction. We expect that this understanding is important to efficiently predict packing states from initial conditions by, e.g., data-driven approaches based on machine learning.

In contrast to the continuous packing dynamics of dipods, the flow of the tripod systems investigated by us is affected by a “key-and-lock” interlocking mechanism. We identify the emergence of an average clustering coefficient as the dominant feature that governs the flow behavior. Specifically, we find that this feature acts as a binary “on/off” parameter in our symbolic regression model, indicating that tripod stability is not governed by local density minimization, but by the formation of distinct contacts. This distinction suggests a fundamental shift in granular flow behavior as the complexity of particle shape increases: from a “packing-dominated” regime (observed for dipods) where stability arises from maximizing density (i.e., “Voronoi minimization”), to an “interlocking-dominated” regime (observed for tripods), where stability arises from discrete geometric constraints. This transition mirrors findings in entangled granular media, where increasing non-convexity shifts the system from fluid-like caging (Torquato and Stillinger, 2010) to solid-like interlocking (Murphy et al., 2016).

4.2 Outlook

Future work should extend this framework to higher-order X-pods, such as tetrapods, hexapods, or octapods, to systematically explore how the degree of chaos scales with particle complexity. A key question is whether higher-order symmetries lead to increasingly deterministic packing behaviors, or introduce new modes of disorder. Additionally, incorporating particle flexibility, e.g., allowing pods to compress and expand, and moving from inherently transient to (pseudo) steady-state simulation setups would provide a more comprehensive understanding of real-world granular materials. Direct experimental comparisons are expected to be most feasible in low-ξ / D regimes where neighboring spheres can physically contact, enabling the physical fabrication of pods. For larger ξ / D values, the lack of intra-pod contacts makes fabrication increasingly impractical, rendering one-to-one experimental benchmarking across the full parameter space difficult or even impossible. Beyond these practical constraints, it should be noted that the pods studied here are abstract, non-convex model particles rather than physical objects. Direct experimental comparison is therefore not precluded, but would deviate from our intended concept of using abstracted particle geometries to model the packing behavior of granular matter. Accordingly, our present work must be understood as a systematic exploration of the geometric parameter space of abstract non-convex particles. Such particles are not meant to represent the shape of real-world particles, but to phenomenologically capture their flow and stacking behavior.

On the theoretical front, future efforts should also focus on implementing symbolic regression for systems with higher particle counts or more diverse, non-regular particles. Our initial efforts to do so failed, suggesting that accurately predicting the interlocked states of large particle ensembles remains challenging with current methods. This points towards the need for new descriptors that can better resolve subtle orientational differences in systems with N > 2.

By expanding our dataset to include larger particle counts, future work could aim at deriving generalized predictive laws for final configurations beyond regular particles (i.e., particles that can be described by a simple geometrical shape). Ultimately, leveraging established links between geometric parameters and macroscopic behavior would offer a pathway for the inverse design of granular materials. This would enable the engineering of particles with tailored packing densities and other (mechanical) properties, making optimal particles for industrial or architectural applications a reality.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript. Portions of the post-processing code and selected textual sections were drafted or refined with the assistance of GitHub Copilot (Claude Sonnet 4.5, Gemini 3 Pro). All outputs were reviewed and verified by the authors.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fceng.2026.1766550/full#supplementary-material

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Nomenclature

Acronyms

2D Two-dimensional 3D Three-dimensional BOO Bond Orientational Order CLT Central Limit Theorem COM Center of Mass CMIM Conditional Mutual Information Maximization DBSCAN Density-Based Spatial Clustering of Applications with Noise DEM Discrete Element Method HPC High-performance computing N/A Not available NaN Not a Number PDP Partial Dependence Plot p Particle pp Particle-particle pw Particle-wall QPSALO Quantum-behaved Parallel Swarm Ant Lion Optimization R ^**2** Coefficient of Determination RDF Radial distribution function RMSE Root Mean Square Error SDIC Sensitive Dependence on Initial Conditions SHAP Shapley Additive Explanations TDA Topological Data Analysis w Wall XY XY-plane

Latin characters

D * Sphere diameter [m] E Relative error [–] e Coefficient of restitution [–] E k i n , i * Kinetic energy at time i [J] E pot , 0 * Initial potential energy [J] m i * Mass of particle i [kg] N Number of pods [–] s Sample standard deviation [–] n Sample size [–] n s t a t Statistically required sample size [–] t * Time (or time step) [s] t Critical value from student t-distribution [–] v i * Velocity of particle i [m/s] v y * Velocity component in y-direction [m/s] x i * Position of particle i [m] x Cartesian coordinate in x-direction [–] x ̄ Sample mean [–] Y * Young’s modulus [J/m³] y COM Normalized height of COM [–] y COM * Absolute height of COM [m] y COM ̄ Normalized mean height of COM [–] y Cartesian coordinate in y-direction [–] z Cartesian coordinate in z-direction [–]

Greek characters

α Significance level (confidence interval) [–] μ Population mean [–] μ r r Rolling friction coefficient [–] μ r s Sliding friction coefficient [–] ν Poisson ratio [–] ρ p * Particle density [kg/m³] φ Packing fraction [–] φ 0 Initial packing fraction [–] ξ * Pod radius [m] ξ / D Interlocking parameter (i.e., pod radius to sphere diameter ratio) [–]

Footnotes

Edited by: Chuan-Yu Wu, University of Surrey, United Kingdom

Reviewed by: Sina Zinatlou Ajabshir, University of Salerno, Italy

Chao Zheng, China University of Petroleum, China