Emerging contact force heterogeneity in ordered soft granular media

Under external perturbations, inter-particle forces in disordered granular media are well known to form a heterogeneous distribution with filamentary patterns. Better understanding these forces and the distribution is important for predicting the collective behavior of granular media, the media second only to water as the most manipulated material in global industry. However, studies in this regard so far have been largely confined to granular media exhibiting only geometric heterogeneity, leaving the dimension of mechanical heterogeneity a rather uncharted area. Here, through a FEM contact mechanics model, we show that a heterogeneous inter-particle force distribution can also emerge from the dimension of mechanical heterogeneity alone. Specifically, we numerically study inter-particle forces in packing of mechanically heterogeneous disks arranged over either a square or a hexagonal lattice and under quasi-static isotropic compression. Our results show that, a hexagonal packing exhibit a more heterogeneous inter-particle force distribution than a square packing does. For both packing lattices, preliminary analysis shows the consistent coexistence of outliers (i.e., softer disks sustaining larger forces while stiffer disks sustaining smaller forces) in comparison to their homogeneous counterparts, which implies the existence of nonlocal effect. Further analysis on the portion of outliers and on spatial contact force correlations suggest that the hexagonal packing shows more pronounced nonlocal effect over the square packing under small mechanical heterogeneity. However, such trend is reversed when assemblies becomes more mechanically heterogeneous. Lastly, we confirm that, in the absence of particle reorganization events, contact friction merely plays the role of packing stabilization while its variation has little effect on inter-particle forces and their distribution.


Introduction
Upon external perturbations, inter-particle contact forces in disordered granular media are well known to form, both experimentally (Majmudar and Behringer, 2005) and numerically (Radjai et al., 1996), a spatially heterogeneous distribution with filamentary patterns (i.e., force chains). These forces and together with the distribution are well known to play a pivotal role in determining how granular media collectively behave (e.g., shear banding (Kawamoto et al., 2018) and solid-liquid phase transitioning (Andreotti et al., 2013;Li and Andrade, 2020)) and interact with external stimuli (e.g., intruder impact (Clark et al., 2015) and wave propagation (Zhai et al., 2020)). Understanding them is therefore relevant to many applications in engineering (e.g., designing adaptive devices (Daraio et al., 2006;Wang et al., 2021)) and (geo-) physics (e.g., mitigating geophysical hazards (Johnson and Jia, 2005)).
Numerous studies have shown that features of inter-particle forces and the distribution depend non-trivially and sensitively on the specific L. Li et al. to quantify contact forces in rigid particle packings. Very recently, these two numerical methods are also being extended to study the collective compaction behavior of highly deformable (Vu et al., 2019;Cantor et al., 2020) or compressible (Vu et al., 2021) particle packings, or that of bi-mixtures of rigid and deformable particle packings . These experimental and numerical studies found that, for disordered disk or sphere packings, both the normal and tangential (frictional) contact forces show exponential distributions for strong forces (i.e., those above the mean) and show power-law distributions for weak forces (i.e., those below the mean) (Majmudar and Behringer, 2005;Radjai et al., 1996). Further, the observation of exponentially distributed strong normal forces seems to be insensitive to particle shape variation based on studies investigating rigid polyhedron (Azéma et al., 2009) and deformable ellipse packings (Wang et al., 2021b), although the associated scaling exponent depends on particle shape. In addition, strong normal forces were found to gradually switch from showing an exponential distribution to showing a Gaussian distribution as a packing's size polydispersity decreases (Voivret et al., 2009), with the Gaussian distribution being discovered in ordered packings of frictionless and rigid disks (van Eerd et al., 2007). Lastly, it was found that, as far as normal contact forces at the boundaries of sphere packings are concerned, they show an exponential distribution whose shape is insensitive to contact friction (Blair et al., 2001). However, nearly all studies to date have been focusing on granular media composed of mechanically homogeneous (being either rigid or deformable) particles, leaving the aspect of mechanical heterogeneity a rather uncharted area. The first attempt to aim at exploring the aspect of mechanical heterogeneity, to the best of our knowledge, dates back to 1986 where a set of experiments were performed to investigate force transmissions in bi-mixtures of plexiglass and rubber particles arranged over a hexagonal lattice (Travers et al., 1986). This experimental study suggested that -though only qualitatively -heterogeneous contact forces can also be induced by only mechanical heterogeneity. Unfortunately, further investigations along this aspect have since remained largely undeveloped. As a result, it still remains unclear how contact forces are distributed in mechanically heterogeneous granular media.
In this paper, we attempt to study quantitatively via simulations inter-particle forces and the distribution within granular media composed of mechanically heterogeneous particles. We consider it an interesting and important problem, not only because of its relevance to many geophysical applications (where geo-materials can be highly heterogeneous mechanically), but due to a more fundamental aspect of opening a potential avenue of engineering novel granular media through a bottom-up perspective. Such bottom-up engineering may be achieved, in the future, through a tactical combination of mechanical heterogeneity and geometric heterogeneity, which in turn allows us to actively control the contact force distribution of granular media (e.g., achieving a homogeneous inter-particle force distribution). As a point of departure, in this paper we numerically approach this problem in its minimal dimension possible: we isolate the dimension of mechanical heterogeneity by considering packing of mono-sized disks (plane-strain cylinders) arranged over two canonical lattices: a square lattice and a hexagonal lattice. We also assume small deformation within every disk such that point-like contacts and contact forces can still be rationally defined. Lastly, we assume that no particle reorganization occurs upon quasi-static loading. Here the phrase ''particle reorganization" refers to particle movements that involve dynamical events and potential finite deformations that can also lead to the loss of inter-particle contacts. Examples are abrupt dynamical frictional slips between two contacting disks that may also cause large disk rotations and the loss of contacts between the two disks. Essentially, we restrain our scope to studying disk assemblies where every disk is able to achieve static equilibrium given the surrounding contact tractions (normal and frictional tractions) in the limit of small deformation. With the scope being defined, we wish to explore the following three questions: • At the system scale, what does the inter-particle force distribution look like and how does it differ between the two packing lattices? • At the particle scale, how does the variation of the mechanical property of a disk alters the amount of force the disk sustains and how does it differ between the two packing lattices? • At both the system and the particle scale, does contact friction play a role for either lattice?
The rest of the paper is organized as follows. In section. 2, we briefly introduce an implementation of a 2D FEM multi-body contact mechanics algorithm whose details together with benchmark tests are presented in Supplementary Material. In section. 3, we apply the implemented algorithm to model ordered packing under quasi-static isotropic compression and analyze inter-particle forces on both the system scale and the particle scale. We also discuss the effect of contact friction. In section. 4, we conclude with a brief summary of our findings and the inspired outlook for future work.

Modeling methodology
Within our scope of small deformation and no particle reorganization in the quasi-static limit, the implementation is greatly simplified. The central idea is to find iteratively the displacement field such that the resulting contact traction together with all other boundary conditions, equilibrate each solid body in a granular system under consideration. In turn, when equilibration is not possible, we take it as a sign of particles undergoing reorganization (e.g., induced by frictional instabilities) being inevitable. 1 We adopt the classical penalty formulation to model contacts between solid bodies. Using the penalty formulation allows us to readily detect contact between two adjacent disks. We pay extra attention to pick the appropriate values for penalty parameters (the normal contact stiffness and the tangential contact stiffness ) that can capture reasonably well the contact physics but at the same time prevent numerical instabilities from happening. During the early stage of the implementation, we consulted the book (Laursen, 2013) and the paper (Simo and Laursen, 1992) for general theoretical perspectives. The source code is publicly accessible through https://github.com/liuchili/2D-FEM-multibody-contactmechanics.git. Details of the implementation are discussed from a topdown perspective in Supplementary Material. Interested readers can consult the Supplementary Material for a quicker understanding of the implementation.

Modeling ordered packings of disks under quasi-static isotropic compression
In this section, we use the developed implementation to model ordered packing of mono-sized disks. We first introduce the model setup, the calibration of contact parameters, and the preparation of initial configurations. After that, we discuss the simulation results concerning inter-particle forces and the distribution on both the system scale and the particle scale.

Virtual experiment setup
We consider mono-sized disks (with radius = 5 mm) arranged spatially over two canonical packing lattices: a square one with 625 disks (see Fig. 1(a)) and a hexagonal one with 711 disks (see Fig. 1(b)). These two assemblies are confined in two similar-sized rectangular domains respectively and subjected to quasi-static isotropic compression under plane-strain condition. We choose such a system size (500 ∼ 1000 particles) to be consistent with the commonly adopted system size in L. Li et al. The left and right boundaries are designed to have a zig-zag shape for their surfaces touching the disks, which aims at providing a more homogeneous contact force distribution along the boundary. The right boundary is under the same constraints as the right boundary in the square packing, while the left boundary is constraint entirely from moving, to mimic the fixed rigid wall used as the left boundary in the square packing. real experiments (e.g., photoelasticity) that investigate the particlescale and meso-scale physics of granular materials (Bassett et al., 2015;Papadopoulos et al., 2016). For the square packing, the bottom and left boundaries are treated as stationary rigid walls, while the right and top boundaries are treated as rectangular-shaped solids and are both subjected to a constant compressive force and are constrained along the direction perpendicular to the direction of , as shown in Fig. 1(a).
For the hexagonal packing, all boundaries are treated essentially the same as for the square packing case, except that the inner surfaces (those touching the disks) of the left and the right solids are changed to have zig-zag shapes, as shown in Fig. 1(b). We make such a change to increase the homogeneity of contact forces for the hexagonal packing in the reference configuration where every disk has the same mechanical properties.

Input of mechanical properties
We consider disks whose mechanical properties are close to those of rubber-like materials. Since rubber is mostly incompressible ( rubber = 0.5), we consider a ''plane-strain-equivalent" material whose and satisfy = 8∕9 rubber and = 1∕3, where rubber is the Young's modulus of the rubber. This relation can be deduced (see (Vu et al., 2019) for detail) by matching the strain energy density between a rubber-like material and its ''plane-strain-equivalent" counterpart and demanding the in-plane principle stretches to be the same. In our virtual experiments, we remain = 1∕3 unchanged, and we sample for each disk from a truncated Gaussian distribution with a mean mean = 2.75 MPa, a lower bound min = 0.5 MPa (which is close to the material used in Vu et al. (2019)), and an upper bound max = 5 MPa (which is close to the material used in Hurley et al. (2014)). We use the truncated Gaussian distribution for its convenient approximation of both a uniform distribution (by picking a relatively large standard L. Li et al. deviation) and a Dirac-delta distribution (by picking a relatively small standard deviation). Specifically, we vary the standard deviation = 0.125, 0.25, 0.5, 1, 2, 4, and 32, and for each we sample 10 different configurations to get meaningful statistics. Fig. 2 shows the normalized probability density function for each standard deviation . When = 32 we are essentially sampling from a uniform distribution from = 0.5 MPa to = 5 MPa, whilst when = 0.125 we are essentially sampling from a peaked distribution that is very close to the Dirac-delta distribution ( mean ).

Determination of penalty parameters
It is important to pick the appropriate penalization parameters ( and ) for a contact problem simulation, especially in our cases where contacts happen among disks with different mechanical properties. Ideally, we will need to pick values of and that are as large as possible to approximate as close as possible the physical contact laws, which require no normal inter-disk penetration and no tangential inter-disk slip when frictional tractions are below the thresholds set by normal tractions and the contact friction. Also, values of and needs to be larger when the considered contacting solids are stiffer (e.g., having a larger Young's modulus). In our cases, if the values of and are large enough to physically capture the contact mechanics between stiffest disks ( = 5 MPa), and at the same time if such values are not overly large so that the contact interaction softest disks ( = 0.5 MPa) is free from numerical instabilities, we will be able to accurately model contact mechanics between any disks that are between the softest and the stiffest. One more factor to consider is to pick the appropriate number of elements/nodes per disk that is computationally feasible for us. In light of these considerations, we perform displacement-controlled (with = 0.002 mm) isotropic compression tests on a single disk with 10 loading steps ( Fig. 3(a)), considering both = min and = max , and considering both a dense mesh with 2321 nodes (top figure in Fig. 3(b)) and a coarse mesh with 167 nodes (bottom figure in Fig. 3(b)). The goal is to find appropriate values of and that can quantitatively capture contact forces using the coarse mesh by comparing to results obtained from the dense mesh. We find that when is larger than 200 MPa/mm (taking = ), for = max the resulting contact force no longer changes appreciably, at least for the range of considered loading steps. We then apply the same and to a case using the coarse mesh and find good agreement (Fig. 3(c)). However, this value of 200 MPa/mm is too large for cases with = min to converge, and after calibration we find a value of 150 MPa/mm is a suitable choice, as (1) it can converge simulations with = min giving accurate results, and (2) it reduces negligibly the accuracy for simulations with = max , at least for contact forces smaller than 30 N.
Lastly, we note that due to the symmetry of the isotropic compression configuration, we find the above results insensitive to the specific value of (we tried with = 0 and = 0.5). Based on the above discussions, for our virtual isotropic compression experiments, we use = = 150 MPa/mm, and we apply an incremental load of = 250 N with three loading steps. At the last loading step with = 750 N the resulting contact force will be around 30 N which is in the accuracy range of using = = 150 MPa/mm. Lastly, we pick a contact friction = 0.5 (a good choice for rubber-like material similar to (Li et al., 2019)) between disks, and we set the contact friction between disks and the four boundaries to be zero.

Preparation of initial configurations
For each packing lattice, similar to the procedure adopted in Liu and Sun (2020), we prepare an initial configuration by stacking disks with slightly larger radii 5.01 mm and relaxing the configuration with the four boundaries being held fixed. We apply this protocol to generate a very small overlap between two disks. Such overlap serves as a good initial guess for our implementation to find the solution iteratively, and it aids the convergence of our computations. We assign mean to every disk and relax the configuration multiple times with = = 150 MPa/mm, to get sufficiently small inter-disk overlap. Multi-step relaxation can be realized by setting , −1 tot = after we finished the first − 1 relaxation steps, updated the configuration along the way and before we start the th relaxation step. As expected, the inter-disk overlap becomes smaller and smaller as indicated by the distribution of maximum shear stress within each disk (Fig. 4). For both packing lattices, we encounter convergence issues at the fourth relaxation step, which implies that the inter-disk overlap has become sufficiently small. Accordingly, for each packing lattice, we take the relaxed and updated configuration after the third relaxation as the initial configuration of our virtual experiments.

Contact force heterogeneity: intensity
We first investigate the system-level heterogeneity variations of contact forces in both packing lattices as is varied. We first quantify such system-level heterogeneity by computing the standard deviation of all contact force magnitudes, terms as [ ] , in a simulated configuration. To obtain [ ], we compute ↔ (the vectorial form of ↔ ) by ↔ = ( ↔ − ↔ )∕2 between every two contacting disks, where ↔ and ↔ are the summation of forces on active nodes of solid ( ) with respect to solid ( ) and that on active nodes of solid ( ) with respect to solid ( ). We take a minus sign between these two quantities as  they satisfy, in theory, ↔ + ↔ = because of Newton's third law. We take an average between these two quantities to reduce possible bias on evaluating ↔ due to FEM discretization. Fig. 5(a) shows the variation of [ ] with the variation of for both packing lattices and for the three applied loads = 250 N, = 500 N and = 750 N, in a semi-log plot. Each data point shows the averaged value of [ ] considering ten configurations independently sampled from a truncated normal distribution with a given , together with the error bar indicating the variation. For both packing lattices, [ ] first increases with the increase of and plateaus when goes beyond 2 (the dashed vertical line). However, [ ] from a square lattice increases much faster when is smaller than two compared to that of a hexagonal lattice. In addition, the difference of [ ] between the two packing lattices increases with the increase of external load , with [ ] of a square lattice being consistently greater than that of a hexagonal lattice for all . Alternatively, if we quantify the systemlevel heterogeneity using [ ∕ 0 ] instead of [ ] , as shown in Fig. 5(b), we find that contact forces in a hexagonal lattice is actually more heterogeneous than those in a square lattice. Here [ ∕ 0 ] indicates the standard deviation of [ ∕ 0 ] where [ 0 ] represents the collection of contact force magnitudes from the reference configuration in which all disks share the same ( → 0). This observation can be made more conclusively by looking at the probability distribution of normal contact force (which is discussed in more detail in the next paragraph), as shown in Figs. 6(a)-(h), where the hexagonal packing shows a more widespread distribution than the square packing does. The apparent discrepancy between [ ] and [ ∕ 0 ] may be explained by the fact that a hexagonal lattice allows for larger coordination number (⟨ ⟩ = 6) over a square lattice (⟨ ⟩ = 4). As a result, in an absolute term ( [ ] ), a contact in a square lattice sustains on average a larger contact force compared to a contact does in a hexagonal lattice, thus giving greater contact force heterogeneity under the same applied load. Additionally, a larger coordinate number allows for more options for contact forces to distribute in space, thus in a relative term ( [ ∕ 0 ] ), enabling greater contact force heterogeneity. Moreover, we find that, once normalized by [ 0 ], the relative contact force heterogeneity [ ∕ 0 ] become rather insensitive to the considered range of external loading, as shown by the collapse of curves in Fig. 5(b). Lastly, to aid visualization, we show the spatial distribution of maximum in-plane shear strain ( max ) of one configuration sampled from = 32 and subjected to = 750 N from both the square lattice (Fig. 5(c)) and the hexagonal lattice (Fig. 5(d)). We can clearly observe that some disks are under much larger shear deformation than others, in a similar way to how cylinders shine with different intensities in a photo-elastic experiment. However, we note that, unlike in a photo-elastic experiment where brighter cylinders indicate locations of larger contact forces, in our virtual experiments ''brighter" (in terms of shear strains) disks instead suggest locations of smaller contact forces. We present a zoomed-in plot of a small area of Fig. 5(c) and Fig. 5(d), respectively, as shown in Figs. 5(e)(f). It can be observed that larger contact forces (represented by longer and thicker solid black lines) take place generally between disks with smaller shear strains. In addition, the direction of each solid black line is aligned with the direction of the corresponding contact force. Its deviation from the direction of the branch vector (a vector connecting the center of mass of two contacting disks) suggests the existence of frictional forces that arise from non-zero contact friction.
We next move to analyze the probability distributions of contact forces from our virtual experiments and compare the results with the classical ones obtained from disordered packing of rigid disks (termed as DPRD hereafter for simplicity). We first compute the normal contact force magnitudes [ ] and tangential (frictional) contact forces magnitudes [ ] from [ ]. Since we implement the contact law on a stress level instead of on a force level as in ordinary DEM (Cundall and Strack, 1979), for a contact between solid ( ) and solid ( ), we compute ↔ and ↔ in the following way: Suppose that we have ↔ (the vectorial form of ↔ ) already computed following the procedure described in the preceding paragraph, and suppose that due to small deformation we can approximate the contact normal ↔ by the branch vector ↔ = ( ) − ( ) using ↔ = ↔ ∕‖ ↔ ‖ 2 , where ( ) and ( ) are the centroid positions of solid ( ) and solid ( ) in the undeformed configuration. With these quantities at hand, we can compute ↔ = . Then, following the convention adopted in Radjai et al. (1996), we plot the normalized probability distributions, ∕⟨ ⟩ and ∕⟨ ⟩, for both packing lattices subjected to all three loading steps and with = 32 (Fig. 6(a) and Fig. 7(a)), = 1 (Fig. 6(b) and Fig. 7(b)), = 0.5 (Fig. 6(c) and Fig. 7(c)) and = 0.25 (Fig. 6(d) and Fig. 7(d)). Here ⟨ ⟩ and ⟨ ⟩ are the mean normal contact force magnitude and the mean tangential (frictional) force magnitude, respectively. In addition, similar to data presented in Figs. 5(a)(b), data shown in every sub-figure of Figs. 6 and 7 are results from the average of ten configurations independently sampled with a given . First, let us focus on the normalized probability distribution of . Generally speaking, the hexagonal packing shows a more widespread probability distribution of ∕⟨ ⟩ over the square packing, indicating the distribution of contact force magnitudes for the hexagonal packing is more scattered. This general observation is consistent with the observation from Fig. 5(b) that [ ∕ 0 ] from the hexagonal packing is greater than that from the square packing. For the square packing, the probability distribution curve narrows toward the mean (where the probability peaks) as decreases, indicating the convergence to a homogeneous contact force distribution where = ⟨ ⟩ for every contact. For the hexagonal packing, although the narrowing trend also appears, the probability distribution seems to converge to a different type of distribution as decreases. This distribution has multiple peaks ( Fig. 6(a)) whose normal forces correspond to 0 (which is defined per contact) obtained from the reference configuration. Note that, in the reference configuration, unlike the square packing where 0 is the same for all contacts, 0 is not the same for all contacts for the hexagonal packing. This is due to the fact that, unlike the square lattice case, the loading direction (a square type) is not aligned with the lattice direction (a hexagonal type), thus leading to different contact forces among contacts. Interestingly, when we instead plot the probability distribution of ∕ 0 , data from the hexagonal packing no longer show any peak when us small, as indicated for example by Fig. 6(e) in comparison to Fig. 6(a) and by Fig. 6(f) in comparison to Fig. 6(b). Further, data from the hexagonal packing rapidly converge to data from the square packing as decreases, and they both show exponential-like tails for normalized normal force above the mean ( > 0 ) and also for normalized normal force below the mean ( < 0 ), as shown in Figs. 6(e), (f), (g) and (h). We further fit exponential functions of the form exp(− ∕ 0 ) to those tails with characteristic exponents (Fig. 6(f)) weak (for normalized normal forces below the mean) and strong (for normalized normal forces above the mean), and we plot the variation of weak and strong along for both packing lattices, as shown in Fig. 6(i). It can be observed that for both packing lattices, the magnitudes of both weak and strong rapidly decay as increase and almost saturate when goes beyond one. Interestingly, these exponential tails of the normal force distributions observed in our hexagonal but mechanically heterogeneous packings are different from the Gaussian-like tails observed in hexagonally-arranged rigid frictionless packings (van Eerd et al., 2007). The apparent discrepancy observed in the normal force distribution for these two packings may be explained by different protocols from which forces are generated: In our work normal forces are generated from only one loading type (i.e., compressing uniformly along the x and y direction), while in van Eerd et al. (2007) normal forces are generated statistically from all force ensembles as long as the sampled forces satisfy overall stresses and balance each particle. It may be interesting to see whether the normal distribution will have a Gaussian-like tail if we further included normal forces generated from multiple other loading types (e.g., a uniaxial compression).
Next, we focus on discussing the normalized probability distribution of . In contrast to , the hexagonal packing shows a less widespread probability distribution of than the square packing (Figs. 7(a)(b)(c)). However, as increases, results from the hexagonal packing gradually converge to those from the square packing, exhibiting a latticeindependent trend. In particular, when is relatively large (i.e., = 32), results from both packing lattices show nearly identical trend for above the mean (Figs. 7(d)); this trend has a tail decaying faster than an exponential one which is commonly observed in DPRD. For below the mean the probability keeps increasing as approaches zero, which is in qualitative agreement with the power law scaling observed in DPRD for below the mean. In all, by comparing the probability distributions of and for both packing lattices, we conclude that, in the absence of particle reorganization, plays a more dominant role over in determining the degree of contact force heterogeneity of an ordered packing. We close this section by presenting two figures visualizing the contact force distributions for both the square lattice (Fig. 8) and the hexagonal lattice ( Fig. 9) obtained from configurations sampled from different and subjected to different . Contact forces in both figures are scaled with the same constant. We can clearly observe the emergence of preferred locations of contact forces as increases, in a way reminiscent of force chains observed in DPRD.

Contact force heterogeneity: orientation
In natural disordered granular media, it is well known that heterogeneity emerges not only as heterogeneous force intensities, but also as heterogeneous force orientations, and that both heterogeneities give rise to the ability of a natural disordered granular packing to resist an external loading. Accordingly, in this section, we shift our attention to analyze the contact force heterogeneity in terms of the orientation, instead of the intensity as discussed in section. 3.5.1.
We propose a parameter to quantify the force orientation heterogeneity in a granular packing. We define as the deviation from the direction of a contact force, , to that of the contact force in the reference configuration, 0 . In the reference configuration where every disk has the same mechanical property, 0 is either horizontal or vertical in a square packing, while in a hexagonal packing 0 can be horizontal, or 60 • from being horizontal, or 120 • from being horizontal. We compute both the average, [ ] , and the standard deviation, [ ] , of for both packing lattices considering different values of , as shown in Fig. 10(a) and Fig. 10(b). Both results suggest that for both packing lattices the level of force orientation heterogeneity increases with the increase of mechanical heterogeneity (i.e., the increase of ). In addition, when is small, the force orientation heterogeneity of the square packing is weaker than that of the hexagonal packing, but it quicks converges to that of the hexagonal packing as increases. These two observations can also be made by computing the probability distribution of at different values of , as shown from Fig. 10(c) to Fig. 10(f) where values of are 0.25 MPa, 0.5 MPa, 1 MPa and 32 MPa, respectively. The observed generally stronger force orientation heterogeneity in the hexagonal packing may be explained by the larger coordination number the hexagonal packing possesses in comparison to the square packing, which gives the hexagonal packing more possibilities for contact force orientations.

Particle-scale gain of contact forces
Observing the interesting patterns presented in Figs. 8 and 9, it is then natural to wonder about a possible correlation between the mechanical property of a disk and the amount of contact force that disk sustains. Or in other words, does the mechanical property of a disk (in this work just the Young's modulus) play a role in determining the amount of contact force that disk receives? For disks being laterally confined into a one-dimensional chain, the answer is trivially that mechanical properties play no role, since geometrical constraints allows for a single path for contact forces to locate in space to balance the externally applied load regardless of how soft or how compressible a constituent disk is. However, when disks are arranged in 2D (and, of course, 3D) arrays the scenarios become much less straightforward, since there are multiple potential paths allowing for contact forces to form networks to balance the externally applied load. In particular, different mechanical properties of contacting disks lead to non-affine deformation from one disk to another (see for example Fig. 5(c)), which can cause stress redistribution (and subsequently contact force redistribution) that otherwise vanishes in disk packing with homogeneous mechanical properties. In this regard, we define the following two quantities to investigate the possible correlation between Young's modulus and the sustained contact force of a disk: L. Li et al. • = ( ∕ mean − 1) × 100% defined as the relative ''stiffness" variation for a disk with respect to its reference state mean . • = ( ∕ 0 − 1) × 100% defined as the corresponding relative contact force gain/loss for the disk with respect to its referenced state 0 , where is the total contact force sustained by that disk with a given Young's modulus , and 0 is defined similarly but is obtained from the reference configuration in which = mean for every disk.
Solely from an energy point of view without any consideration of the potential spatial correlation among constituent disks, under the same external load, should be positively correlated to , i.e., stiffer disks gain larger contact forces, which corresponds to less work done by the L. Li et al. Fig. 8. Visualizations of the computed inter-particle forces in the square packing whose configurations are sampled from = 32(a)(b)(c), = 1(d)(e)(f) and = 0.5(g)(h)(l). Contact forces are represented by solid black lines whose length and thickness are scaled according to the magnitudes of contact forces. applied external load and consequently less strain energy stored by the disks. 2 First, we focus on discussing the correlation between and for the square packing. We pick five scenarios ( = 0.25, 0.5, 1, 2 and 32) under = 250 N ( Fig. 11(a)), and for each scenario, we plot all and data obtained from the ten independently sampled configurations (upper panel of Fig. 11(b)). We color data in blue if and appear in the first and third quadrants of each sub-figure, i.e., they are positively correlated as > 0, and we color data in red if otherwise ( < 0 residing in the second and fourth quadrants). Two qualitative observations can be made. First, there indeed exist a positive correlation (data in blue) between and when is relatively small ( < 2). However, there are also ''outliers" (data in red) regardless of the particular value of , i.e., softer disks ( < 0) ending up gaining larger forces ( > 0) while stiffer disks ( > 0) ending up gaining smaller forces ( < 0). Second, as increases, the positive correlation becomes weaker and the portion of ''outliers" become greater. The above two observations suggest the existence of nonlocal effect, i.e., how much force a disk sustains depends not only on how stiff that disk is, but also how stiff the surrounding disks are. We accordingly perform a first-order nonlocality check, labeling the stiffness of a disk by a new quantity that considers immediate 2 We assume that energy potentially dissipated through friction is not significant compared to the stored strain energy.
contacting disks: where is the number of contacting disks ( = 4 for a square packing), and is the corresponding Young's modulus of a contacting disk. For disks near the boundaries, we simply set to be infinity (1∕ → 0). Each quantity associated with a contact inside the summation on the right-hand side of Eq. (3.1) can be viewed as an effective stiffness from two springs linked in a serial, each of whose stiffness equals the Young's modulus of a corresponding contacting disk. thus on average quantifies how ''stiff" a disk is by incorporating the stiffness of nearby disks. Accordingly, we use a new quantity = ( ∕ ,mean − 1) × 100% and plot its correlation with (shown in the second panel of Fig. 11(b)). Here ,mean is simply the value of in the reference configuration where = mean for every disk. As shown in the second panel of Fig. 11(b)), compared to simply using , appears to show better correlation with , especially when is relatively small ( < 2). More specifically, correlated data (colored in blue) are less scattered and the portion of ''outliers" (colored in red), particularly of those lying in the fourth quadrant ( > ,mean but < 0 ), reduces considerably. However, when further increases, the relevance of to decreases, though still performs slightly better than . For the hexagonal packing, similarly, at the same five values (Fig. 11(c)), we plot the correlation between and (top panel of Fig. 11(d)), and that between and (bottom panel of Fig. 11(d)). The general trends are very similar to those of the square packing case L. Li et al. but the correlations are much stronger, especially when is relatively large (e.g., = 32). The observed much stronger correlation in a hexagonal packing may be explained by its larger coordination number (⟨ ⟩ = 6) over that of a square packing (⟨ ⟩ = 4), which statistically promotes energetically favored load-bearing paths formed by stiffer disks to balance the externally applied load, thereby implying a less pronounced nonlocal effect. For both packing lattices, however, we do observe the consistent existence of ''soft outliers" (those with < 0 but > 0) regardless of the use of or . Such observations suggest that there are non-negligible spatial correlation among constituent disks, and that such correlations appear to be long-range extending beyond a first-order nonlocality.
To shed further light on nonlocality, we proceed to quantify the spatial force correlations within these ordered packings. Following (Majmudar and Behringer, 2005), we compute the two-point correlation function ( ) = ⟨ ( ) ( + )⟩ , where is the sum of the magnitudes of the contact forces on a particle, 3 and is the position of its centroid. In order to investigate the correlation along different directions, we do not average over the angle. We pick the same four scenarios as earlier ( = 0.25, 0.5, 1 and 32) under = 250 N, and for each scenario, we plot the correlation (averaged over the ten independently sampled configurations) against the normalized radial distance along various angles (Fig. 12). In particular, the left column (Figs. 12(a), (c), (e), and (g)) corresponds to the square packing, while the right column (Figs. 12 (b), (d), (f), and (h)) corresponds to the hexagonal packing. In each plot, the insets reveal the polar contours of spatial correlation. The immediate observation is that, in both arrangements, the spatial force correlation indeed extends far beyond the first neighbors. This is consistent with the observation of force chain-like structures in Section 3.5.1. Despite the isotropic loading conditions, and in contrast to observations in DPRD, the correlation is highly anisotropic, and reflects the orientation of contacts in each arrangement. Interestingly, in the case of the square packing, the anisotropy appears to be independent of the degree of mechanical heterogeneity. On the contrary, the hexagonal packing with low mechanical heterogeneity exhibits a spatial distribution of forces whose correlation is more pronounced in the horizontal direction, while the same packing with large mechanical heterogeneity shows a spatial correlation that is equally pronounced along the six contact directions inherent to the arrangement. We postulate that this is due to the directional bias of the contact forces, inherent to the homogeneous hexagonal packing, which becomes less pronounced as heterogeneity increases. Finally, in Fig. 13(a), we plot the evolution of the normalized correlation lengtĥ= ∕ (where denotes the particle diameter) as a function of the material heterogeneity ( ). Note that the correlation length is obtained by fitting an exponential correlation kernel ( ) = exp(− ∕ ) to the angle-averaged correlation data, and is indicative of the characteristic size of particle chains and clusters that are responsible for force transmissions. In an isotropic (angle-averaged) sense, the hexagonal packing exhibits a slightly larger correlation length than the square packing in the low heterogeneity regime (e.g., < 2 MPa), while the opposite is true in the large heterogeneity regime (e.g., ≥ 2 MPa). This phenomenon is more L. Li et al. pronounced in the case where the correlation length is computed along the horizontal direction, as shown in Fig. 13(b). These observations are in line with our expectations of the portion of ''outliers" (discussed in Section 3.5.3) with respect to the total number of disks computed and shown in Fig. 13(c), where the hexagonal packing shows a larger portion of ''outliers" over the square packing when is relatively small (e.g., < 2 MPa), while the opposite is true when goes beyond 2. Finally, as expected, both packings show an overall decaying correlation length upon increasing material heterogeneity.

Effect of contact friction
We have also tried varying the contact friction to be either smaller or larger than 0.5, and we find has little effect on our findings so long as we are able to get converged solutions. In our case non-convergence happens when we reduce , and it happens more frequently with the increase of . On the contrary, when > 0.5 we are always able to find converged solutions. These non-converging scenarios imply reorganization events (e.g., large rotations of particles induced by frictional instabilities) which can also lead to loss of contacts. 4 For instance, we find that for < 0.35 we were unable to find converged solutions for hexagonal configurations sampled from = 32. It is possible in these scenarios that a relatively soft particle loses all contacts, while the surrounding particles are relatively stiffer and 4 We rule out the possibility of numerical instabilities (e.g. those induced by poor mesh qualities) since our simulations converge for certain values of .
are connected in a way capable of sustaining the external load. When reorganizations happen under smaller they allow greater flexibility among particles to explore and form more energetically favored loadbearing paths, which may then lead to a more heterogeneous contact force distribution and a stronger correlation between and . Extending our findings to these scenarios would either require an extension our current implementation to be dynamic, or experimentations with techniques (Hurley et al., 2014;Marteau and Andrade, 2017) that can measure inter-particle forces among different types of materials.

Summary and outlook
In this work, we numerically explore the effect of mechanical heterogeneity on inter-particle forces in soft granular media using a 2D FE contact mechanics algorithm. Specifically, we study two canonical packing lattices -a square lattice and a hexagonal packing latticesunder quasi-static isotropic compression. For both packing lattices, we show that a heterogeneous inter-particle force distribution emerges as the Young's moduli of constituent disks gradually deviate from being homogeneous, despite disks being arranged orderly in the absence of geometric heterogeneity.
At the system level, under the same level of mechanical heterogeneity, we observe that the hexagonal packing lattice shows a more heterogeneous inter-particle force distribution than the square packing lattice does. This observation may be explained by the larger coordination number of the hexagonal packing lattice that promotes more load-bearing paths. Correspondingly, we find that, on the one hand, for normal force well above the mean, the probability distribution from the for the square packing subjected to = 250 N. We pick five values and plot the correlation between and (top panel of (b)) and that between and . Data are shown by overlaying results from the ten independently sampled configurations for each . Data colored in blue locating in the first and third quadrants, represent disks whose stiffness correlate positively to the contact force they sustain, i.e., > 0 or > 0. Data colored in red locating in the second and fourth quadrants, represent ''outliers" including softer disks ( < 0 or < 0) gaining larger forces ( > 0) and stiffer disks ( > 0 or hexagonal packing shows a longer tail than that from the square packing. For normal forces well below the mean, the probability distribution from both lattices show tails dipping toward zero. However, when the probability distribution is plotted by normalizing normal forces against their reference values (i.e., 0 ), it shows a much weaker latticedependent trend with both forces below and above the mean exhibiting exponential tails. On the other hand, tangential (frictional) forces well above the mean in both lattices show tails decaying faster than exponential ones that typically appear in DPRD. Finally, both studied systems exhibit long-range spatial force correlation, which is consistent with our observations of emerging force chain-like structures.
At the particle scale, both packing lattices show beyond-first-order spatial correlation (i.e., nonlocal) effect in the sense that the amount of contact force a disk sustains can be determined neither by considering the stiffness of that disk alone nor by further considering the stiffness of immediately contacting disks. Specifically, for both packing lattices we also observe the coexistence of ''outliers": softer disks gaining larger forces and stiffer disks gaining smaller forces. Additional analysis on spatial force correlation confirms that the spatial correlation effect is indeed beyond first order. Further, the analysis suggests that the hexagonal packing lattice shows strong nonlocal effect over the square packing when is relatively small (i.e., < 2) and the opposite holds when goes beyond 2, which is in line with the observation of the portion of ''outliers" in a hexagonal packing being larger than that of the square packing when < 2 and being smaller otherwise.
Concerning the effect of contact friction, we find that so long as no particle reorganization occurs, our findings are insensitive to the particular value of . However, our simulations show that more mechanically heterogeneous packing requires a higher to prevent reorganizations from happening, suggesting a role of packing stabilization played by contact friction.
Looking forward, we present several potential future research directions. The first direction concerns a deeper understanding aiming at correlating the amount of contact force a disk sustains to the mechanical characterization of that disk. Our simulations show that even in simple ordered packing, it is challenging to account for the range of nonlocal effect. We believe that it is promising to tackle this challenge through a combination of network theory and machine learning. Specifically, network theory allows us to extract communities (Bassett et al., 2015;Karapiperis and Andrade, 2021), thus implicitly taking into account the nonlocal effect. These communities serve as excellent sources from which one can extract multiple descriptors associated to a single disk in a way similar to (Cubuk et al., 2015). From these descriptors we could train machine learning (ML) algorithms to identify relevant descriptors, first as a labeling problem (i.e., correctly identifying ''soft" and ''stiff" disks) and later as a regression problem (i.e., quantitatively predicting the amount of contact forces). These trained ML models may be used to further explore possible finite size (boundary) effects. The second direction concerns extending our findings to regimes where Fig. 12. Spatial two-point correlation of the total magnitude of the force acting on a particle, evaluated along specific directions, and plotted as a function of the normalized distance. Plots (a), (c), (e), and (g) correspond to = 32, 1, 0.5, 0.25 MPa respectively for the square packing. Similarly, plots (b), (d), (f), and (h) correspond to = 32, 1, 0.5, 0.25 MPa respectively for the hexagonal packing. Insets show the polar contours of correlation for each case. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) particle reorganizations occur. Reorganization events allow granular media to explore more packing configurations, thus potentially leading to more energy-favored loading paths. It would be interesting to investigate how reorganization events change the inter-particle forces and their distribution. Such investigations can be achieved by extending our implementation to be fully dynamic (and possibly to account for finite kinematics), or through experiments using techniques capable of measuring inter-particle forces between particles made of different materials (Hurley et al., 2014;Marteau and Andrade, 2017). A final promising direction relates to the development of analytical models of force transmission in these mechanically heterogeneous packings. In this regard, it would be worthwhile studying how theories such as L. Li et al. for the square (blue) and the hexagonal (red) packing. Error bars indicate one standard deviation from the ten independently sampled configurations for each . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) the q-model (Liu et al., 1995) could be extended and adapted for these systems.
In all, by exploring the ''disordered" space in terms of particles' mechanical properties, our work offers a fresh perspective to the classical understanding of inter-particle forces in granular media. It is in the hope of the authors that this work will promote further investigations along this direction.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.