Formations of force network and softening of amorphous elastic materials from a coarsen-grained particle model

Amorphous materials, such as granular substances, glasses, emulsions, foams, and cells, play significant roles in various aspects of daily life, serving as building materials, plastics, food products, and agricultural items. Understanding the mechanical response of these materials to external forces is crucial for comprehending their deformation, toughness, and stiffness. Despite the recognition of the formation of force networks within amorphous materials, the mechanisms behind their formation and their impact on macroscopic physical properties remain elusive. In this study, we employ a coarse-grained particle model to investigate the mechanical response, wherein local physical properties are integrated into the softness of the particles. Our findings reveal the emergence of a chain-like force distribution, which correlates with the planar distribution of softness and heterogeneous density variations. Additionally, we observe that the amorphous material undergoes softening due to the heterogeneous distribution of softness, a phenomenon explicable through a simple theoretical framework. Moreover, we demonstrate that the ambiguity regarding the size ratio of the blob to the force network can be adjusted by the amplitude of planar fluctuations in softness, underscoring the robustness of the coarse-grained particle model.

Amorphous materials, including granular substances, glasses, emulsions, foams, and cells, find widespread applications as building materials, plastics, food products, and thermal insulators in our daily lives.The mechanical properties of these materials are commonly characterized by parameters such as Young's modulus and bulk modulus, which are pivotal for material design.Despite their adherence to classical mechanics at the level of individual compositions, these amorphous materials often display intricate macroscopic behaviors.For instance, certain components within these materials may exhibit fluid-like behavior, such as hourglass phenomena and avalanches, while also demonstrating solid-like characteristics like clogging hoppers and pipes 1-5 .Therefore, elucidating the complex mechanical responses of amorphous materials is essential for leveraging their diverse functionalities effectively.
In amorphous elastic materials, external forces propagate through localized pathways known as force networks or stress networks [6][7][8][9][10][11][12][13][14] .These force networks have been extensively investigated as crucial components for understanding complex dynamic responses.Over the years, force visualization experiments have been conducted to explore the mechanical behavior of granular materials [6][7][8][9][10][11] .For instance, it has been observed that compressing an unjammed state results in a shear jamming state when the force network percolates in all directions 8 .Furthermore, studies have indicated that the force network becomes isotropic through the vibration of granular materials 9 .
Despite investigations into the features of force networks, the mechanism underlying their formation remains unclear.Furthermore, the relationship between the force network and macroscopic properties, such as the bulk modulus of amorphous materials, remains ambiguous.Understanding these aspects requires studying amorphous materials on a scale larger than their individual components.For instance, the falling dynamics of grains have been likened to Rayleigh-Taylor instability phenomena in fluids [15][16][17][18][19] , with typical length scales of these macroscopic flows several orders of magnitude larger than the grain size.To elucidate such macroscopic www.nature.com/scientificreports/behaviors, a coarse-grained model is often considered effective as a meso-scale statistical mechanics approach.This modeling technique has traditionally been employed to explain phenomenological behaviors in various systems, including polymers, liquid crystals 20 , and phase separations 21 .Recently, coarse-grained models have successfully demonstrated the fracture of amorphous materials 22,23 .In these models, it is assumed that density is heterogeneous, and the shear modulus strongly depends on density.Application of shear leads to enhanced density fluctuations, ultimately resulting in fracture.
Recent molecular dynamics simulations have revealed that the local elastic moduli exhibit spatial distribution in the random packing state [24][25][26] .Moreover, it has been observed that fluctuations in local elastic moduli increase as the jamming threshold is approached 25 .Additionally, features of sound wave propagation in amorphous solids have been found to correlate with the heterogeneous distribution of local elastic moduli 24 .Based on these findings, we hypothesize that the formation mechanism of force networks is also linked to the distribution of local elastic moduli.We characterize the local elastic moduli as the softness of a disk, incorporating factors such as local packing fraction, local structure, and interactions.Subsequently, we investigate the relationship between the force network and the planar distribution of softness using a coarse-grained particle model.

Models and methods
We assume that microscopic physical properties should be averaged at mesoscales, such as the typical length of force networks.Accordingly, we partition the elastic material into blobs with a certain diameter.Within each blob, we characterize the local elastic moduli as softness, considering factors such as local packing fraction, local structure, and interactions refer to Fig. 1.By employing this concept, we investigated a two-dimensional soft disk model wherein the softness of each disk varies.As the blobs encompass the heterogeneous planar distribution of local packing fraction as softness, they are arranged in a close-packed hexagonal lattice configuration.
Each disk denoted as i possesses an elastic parameter G i indicative of its softness.A total of 65536 soft disks with a diameter D are organized within a triangular grid measuring 256 × 256.Periodic boundary conditions are imposed in both the x and y directions.The value of D is fixed at 1, while the positions of the disks along the x and y axes are indexed by the number of layers n x and n y .
We generated a heterogeneous distribution of disk elasticities using the Cahn-Hilliard-Cook equation 21,27 expressed as follows: where t denotes time, Ŵ(r) represents a field determining the softness of particles at position r and ξ stands for the correlation length of fluctuations of Ŵ(r) .Assuming local equilibrium and employing the fluctuation-dissipation relation, the correlation function of the dimensionless noise g(r, t) can be described as where i, j = x, y , δ ij is Kronecker delta and δ(r) denotes a delta function.We conducted numerical simulations of Eq. (1) on a 256 × 256 square lattice.As the order parameter is conserved, Ŵ(r) remains 0. Each lattice point (n x , n y ) corresponds to a blob located at that position.The relationship between Ŵ(r) and G i is defined as G i = �G i � + �Ŵ(r)/(Ŵ max − Ŵ min ) , where G i is the mean value of G i , represents the strength of fluctua- tions of the softness, and Ŵ max and Ŵ min are maximum and minimum of Ŵ in the entire system.We set the mean softness parameter �G i � = 1 .Since Ŵ max ≈ −Ŵ min , a range of values for G i becomes 1 − �/2 ≤ G i ≤ 1 + �/2 .Consequently, ξ represents the correlation length of the fluctuations of G i and represents the amplitude of these fluctuations.Additionally, ξ and are treated as variable independent parameters in this simulation, although the critical point argument for hyperuniformity suggests a potential interdependence between them 28 . (1) Figure 1.Schematic of a coarse-grained particle model for a granular material, serving as an example of amorphous elastic materials.It depicts the hierarchical structure of granular systems, which comprises individual grains, several grains forming a blob, the arrangement of these blobs, and the resulting force network.The two images on the left showcase grains, with scale bars indicating 1 mm.The image on the right visualizes the force using a photoelastic disk, with a scale bar of 5 mm.Additionally, for reference, the hierarchical structure of polymers (including blobs, spring-beads model, and polymer network in gels) is depicted, suggesting an analogy to granular systems.
Taking into account the micro-deformation within the system, the interaction between the disks can be approximated using a harmonic potential.Disks i and j interact via a pairwise harmonic potential: where G ij = (G i + G j )/2 and r ij is a center-to-center distance between disks i and j.Then we obtain where j signifies the summation over disk j, which is in contact with disk i. F i represents the force vector act- ing on disk i, while rij denotes the unit vector pointing from disk i to disk j.We specify that the interaction is purely repulsive, thus f(r) = �(r) where is the Heaviside function.Additionally, we assume a loose packing in granular systems or soft components like emulsions and bubbles.Consequently, compression of the blob and overlap between pairs of disks are permissible.To streamline the model, friction between the disks is disregarded.
We compressed the system by 0.1% in both the x and y directions.Each particle follows the normalized overdamped equation.
where dr i represents the displacement of disk i during each simulation time step dt.We set dt = 0.001.We consider the system to have reached a steady state when the maximum velocity of all particles becomes less than 1.0× 10 −6 .We conducted 5 runs with the same ξ and , but with the different planar distribution of G.We confirmed that the results remained essentially unchanged regardless of the initial conditions.

Result Formation of the force network
Initially, we examined the force distribution.Figure 2a and b illustrate the planar distribution of the top 10% of G ij and the top 10% of the force F ij between disks i and j, respectively.The structure is generated using Eq. ( 1) with ξ = 10 and = 1 at t = 1000.For large sections of the pattern, there appears to be a significant correlation between F ij and G ij .However, the distribution of G ij exhibits a cluster-like behavior, while that of F ij forms a thin network, which is a typical characteristic of amorphous materials [6][7][8][9][10][11][12][13][14] .Figure 3 presents an image plot of the probability of F ij and G ij at ξ = 10 and = 1.It indicates that F ij is qualitatively correlated with G ij , although F ij is widely dispersed around a straight line.Fitting by a linear function yields a slope of 5.8 × 10 −4 .If the displace- ments of the disks were uniform, F ij and G ij would be perfectly proportional, and the slope would be 1.0 × 10 −3 .Hence, the smaller slope suggests spatially inhomogeneous displacement.Additionally, we compute the Pearson correlation coefficient P between F ij and G ij .The obtained value of P = 0.61 indicates a correlation between F ij and G ij , albeit not a strong one.This finding is consistent with the difference in planar distribution between F ij and G ij .Further insight into the deviation of the correlation between F ij and G ij will be provided in the next subsection.
We also examined the force distribution when the planar distribution of G ij remains the same, but only the amplitude of the fluctuation is varied.It was observed that the appearance of the planar distribution of F ij changes minimally.Figure 3b shows the correlation of F ij with the same combination of i and j at = 0.1 and = 1.The data scatter is significantly narrower, and it was found that the slope is 10.0, identical to the ratio of .Additionally, the Pearson correlation coefficient P = 0.994 for F ij between = 0.1 and = 1 is close to 1.This suggests that the planar distribution of F ij is primarily determined by the planar distribution of G ij , rather than .Thus, the characteristic length of the force network is solely determined by ξ , which represents the characteristic length of G fluctuations.
(3) In broad sections of the pattern, there appears to be a significant correlation between F ij and G ij .Meanwhile, the distribution of G ij exhibits a cluster-like pattern, whereas the distribution of F ij forms a network.www.nature.com/scientificreports/

Sandwiched system
To explore the mechanism underlying the reduction in correlation between F ij and G ij , numerical simulations were conducted using a model system.As depicted in Fig. 2, regions sandwiched by high G ij exhibit the larger F ij .To simplify the analysis, a binary system was considered, comprising regions with high ( G i = 1.5 ) and low ( G i = 0.5 ) G i values.Figure 4a illustrates the planar distribution of G i , with high G i regions positioned at the upper and lower centers, referred to as a "sandwiched system".In this configuration, n x and n y are set to 64, with the width and height of each high G i region being 8 and 16, respectively.The system is subjected to a compres- sion of 0.1%.Figure 4b illustrates the spatial distribution of the top 10% of F ij within the sandwiching system.The rec- tangular areas enclosed by solid lines denote the high G i regions.The force is distributed prominently not only within the high G i regions but also within the low G i region sandwiched by the high G i region.Previous studies have highlighted the crucial role of density change δρ in force distribution, given the strong dependence of bulk modulus on density in a jammed state 29 .Thus, we calculated the local density around disk i, denoted as ρ i , utilizing the equation: where ri represents the mean value of r ij averaged over the nearest neighbor disk j.Subsequently, the density change is defined as δρ i = ρ i − ρ 0 , with ρ 0 (= 1) denoting the initial density.Figure 5a   The rectangular area enclosed by the solid line represents the high G i region.The force is prominently distributed not only within the high G i region but also within the low G i region that is sandwiched by the high G i region.
profile as a function of x at the center along the y-axis ( n y = 32).The region sandwiched by the high G i regions is enclosed by dotted lines, as shown in Fig. 5.It is observed that the density within the sandwiched region surpasses that of other low G i regions.Consequently, elasticity within the sandwiched region effectively increases after compression, leading to significant force distribution.
Here we extend the findings derived from the model systems to the fluctuated system, as depicted in Fig. 2. Figure 5b shows δρ plotted against G i , where δρ represents the mean value within the range G i − 0.005 < G i ≤ G i + 0.005 , with error bars indicating the standard deviation.It is observed that the density increases in the lower G i region, consistent with the density change observed in the sandwiched model.This observation suggests that the force network is constructed by connecting the higher G i regions.

Total potential energy
Although the underlying microscopic mechanism behind force chains has been unveiled, the impact of the force network on macroscopic physical properties remains ambiguous.In this study, we quantify the total potential energy within the systems, where the second derivative of the potential energy with respect to the volume corresponds to the bulk modulus.Figure 6a shows the variation of the total potential energy U with respect to .Different symbols such as circles, squares, triangles, and inverted triangles represent U values corresponding to ξ = 1, 3, 5, and 10, respectively.It is observed that U predominantly relies on , with a minor dependency on ξ .This implies that the material undergoes softening when is substantial.The influence of ξ on U will be addressed subsequently.
To explore the relationship between the variation in U and the structural changes, we computed the potential change δU = U 0 − U , where U 0 represents the potential energy in a homogeneous system ( = 0).Figure 6b displays δU plotted against .Our analysis reveals that δU exhibits a proportional relationship with 2 for each ξ value.Given our previous demonstration of how the planar distribution of G i induces density heterogeneity, we further investigated the correlation between δU and the standard deviation of density σ , as depicted in Fig. 6c.Notably, we observed that δU across all systems with varying and ξ values conforms to a scaling law δU ∼ σ 2 .This finding suggests that density change plays a crucial role in determining the total potential energy.
It is essential to recognize that the size of the blob is arbitrary and subject to ambiguity.The influence of this ambiguity on the coarse-grained particle model needs to be considered.Our investigation revealed that the size of the force network is determined by ξ , representing the characteristic length of planar softness inhomogeneity relative to the blob size.As such, changes in the definition of the blob size should correspondingly affect ξ .For instance, in a system with ξ = 10, increasing the blob size by a factor of 10 should theoretically yield equivalence to a system with ξ = 1.However, despite the scale change, slight discrepancies in δU between ξ = 10 and ξ = 1 are observed, as shown in Fig. 6b.To address this discrepancy empirically, we introduce a rescaling factor � r = (0.015ξ + 0.985)� , where r represents the rescaled .Specifically, we select ξ and values such that r matches for different combinations, such as (ξ , �) = (1, 1) and (ξ , �) = (10, 0.881) .We find that δU collapses across different ξ values when r remains consistent, despite variations in the ξ and combinations, as shown in Fig. 6d.Thus, although the blob size remains ambiguous, achieving system equivalence through r rescaling demonstrates the robustness of this coarse-grained particle model.

Theoretical approach
In this subsection, we present a theoretical elucidation of the mechanism behind the decrease in total potential energy with .Initially, we examine a simplified scenario involving a one-dimensional arrangement of three particles, labeled as disks i, j, and k, arranged sequentially.The interactions between disk i and j, and between disk j and k, are denoted by U ij and U jk , respectively.The size of the disks is normalized to 1.We assume that the disk i The density change δρ profile as a function of x at the center along the y-axis ( n y = 32) within the sandwiched system.Despite uniform G i values, notably higher density is observed in the region sandwiched by the high G i regions (within the dotted lines).(b) δρ as a function of G i within the fluctuated system, as shown in Fig. 2. Error bars denote standard deviations.It is notable that δρ is larger in regions with the lower G i , while conversely, it is smaller in regions with the higher G i . is harder than the disk k.Under compression by α , the distance between disk i and j contracts to α/2 − δ , while the distance between disk j and k contracts to α/2 + δ .Subsequently, the total potential energy can be expressed as U = U ij (α/2 − δ) + U jk (α/2 + δ) .Expanding U to the second-order terms with respect to δ , we obtain where U ′ and U ′′ represent the first and second derivatives of U with respect to the distance, respectively.Meanwhile, by considering the balance of forces, we derived Expanding Eq. ( 8) to the first term with respect to δ , we obtain Subsequently, Eq. ( 9) is substituted into Eq.( 7), resulting in: Given that the second derivative of U is positive from the stability of the system, the total potential energy decreases δ = 0 .Finally, we substitute our potential, where www.nature.com/scientificreports/Thus, we obtained δU ∝ � 2 in the simplified model.Next, we explore the scaling law δU ∼ � 2 r ∼ σ 2 in a two-dimensional system.The total potential energy can be expressed as the summation of Eq. ( 3) as follows In a homogeneous system with a compression rate of α = 0.001 , the distance r ij between particles i and j uni- formly becomes D − αD , yielding a total potential energy of U 0 = �G ij �(αD) 2 N(Z/2)/2 = 0.098304 .Here, Z(= 6) denotes the number of particles in contact with each particle.Subsequently, we consider the heteroge- neous case where G ij and r ij are described as G ij = �G ij � + δG ij and r ij = D − αD + δr ij , where δG ij and δr ij are deviations from the average.Consequently, i j δG ij = 0 and i j δr ij = 0 .Given that δr ij ∼ 10 −4 , much smaller than αD and δG ij , the term δr 2 ij becomes negligible.Thus, δU can be approximated as follows : In this context, when the correlation length ξ exceeds 1, the value of G j for disk j, which is in contact with disk i, closely resembles G i .Consequently, we posit that δr ij is approximately equal to δr i and δG ij is approximately equal to δG i , where δr i and δG i denote the mean distance and the mean G ij averaged across neighboring disk j.
Given the observation that higher values of G ij correspond to lower density or larger δr ij , and conversely, smaller δr ij corresponds to lower G ij , we establish the following relationship : Here, r represents the difference of G i , thus Finally, we derive Expressed as σ = 1 N i (ρ 2 i − �ρ�) 2 , the deviation in density can be formulated.Given the smallness of δr ij and α , the expression can be approximated as follows Combining Eqs. ( 15), (17), and (19), it is obtained that δU ∝ σ 2 .

Discussion
In our study, we compare our findings with prior vibration experiments, wherein the force network achieves homogeneity through vibration 9 .It is established in literature that the packing fraction of granular systems increases under vibration, approaching random close packing 1 .Notably, the state of random close packing is identified as hyperuniform, characterized by nearly uniform compressibility across all locations 28,[30][31][32] .Our simulations corroborate that when ξ is small, the propagation force attains uniformity, aligning well with experimental observations.Furthermore, previous studies have reported that the system exhibits a significant response to shear deformation near the jamming point, despite the lack of direct coupling between shear deformation and density changes, unlike compression.For instance, shear-induced decreases in the effective moduli of aggregates have been reported 33 .Additionally, fractures often occur following enhancements in density fluctuations induced by shear stress 22,23 .Investigating the behavior of systems characterized by inhomogeneous G under shear deformation would be of interest.
Additionally, we highlight the significance of the formation of the force network in controlling stress-induced phenomena.Research has shown that stress can suppress liquid-liquid phase separation in cells 34 , and it can also influence the absorption and desorption of guest molecules in metal-organic frameworks (MOFs) 35,36 .Therefore, the characteristics of the force networks are expected to have a substantial impact on the functionality of Vol:.( 1234567890

Conclusion
Many amorphous elastic materials, including glasses, emulsions, foams, and cells, exhibit complex macroscopic behaviors and are extensively utilized in everyday applications, necessitating both mechanical and chemical stability.Despite the recognized importance of the force network, the mechanisms governing its formation and its relationship with macroscopic properties remain unclear.Consequently, a qualitative understanding of the macroscopic properties of amorphous materials is still lacking.In such large systems, a phenomenological approach, such as a simple coarse-grained model, is generally effective.In this study, we investigate the formation of the force network and the variation in total potential energy using a coarse-grained particle model incorporating local packing fraction, particle shape, and microstructure into the softness G of the disk.We observe the formation of a force network resembling that observed in amorphous elastic materials, driven by the correlation length of fluctuations of G. Additionally, we find that amorphous elastic materials undergo softening due to the amplitude of the fluctuations of G, with softening scaling by the square of the density deviation, explicable through a simple theory.Furthermore, we observe that the ambiguity in blob size suggests the robustness of the coarse-grained particle model.In this study, simulations using a coarse-grained particle model suggest that alterations in the effective local elastic modulus, prompted by changes in density, play a significant role.To validate local density changes and changes in local modulus of elasticity, it is essential to incorporate feedback from molecular dynamics simulations of compression.Additionally, we operated under the assumption that deformation is small in this simulation, thus applying the harmonic potential.However, it's crucial to acknowledge that the potential may not remain harmonic under significant deformation.Although the decrease in total potential energy remains consistent, the relationship between total potential energy and density change could potentially vary.Therefore, exploring the potential dependence of our findings is a critical avenue for future research.

Figure 2 .
Figure 2. (a) planar distribution of the top 10% of G ij with ξ = 10 and = 1 .(b) planar distribution of the top 10% of F ij .In broad sections of the pattern, there appears to be a significant correlation between F ij and G ij .Meanwhile, the distribution of G ij exhibits a cluster-like pattern, whereas the distribution of F ij forms a network. https://doi.org/10.1038/s41598-024-59498-2

Figure 3 .
Figure 3. (a) Image plot of the probability of F ij and G ij at ξ = 10 and = 1 .Upon fitting with a linear function, the obtained slope is 5.8 × 10 −4 , which is smaller than the slope when the displacements are uniform (1.0 × 10 −3 ).It is observed that the data exhibit wide dispersion.(b) Image plot of the probability of F ij in a system with = 0.1 and F ij with = 1 for the same combination of i and j.The planar distribution of G ij remains unchanged.It is noted that the scattering of data is significantly narrower.The colors correspond to the probability.

Figure 4 .
Figure 4. (a) The planar distribution of G i , with a portion of the low G i (green) region being sandwiched by the high G i (purple) region.(b) The planar distribution of the top 10% of F ij .The rectangular area enclosed by the solid line represents the high G i region.The force is prominently distributed not only within the high G i region but also within the low G i region that is sandwiched by the high G i region.

Figure 5 .
Figure 5. (a)The density change δρ profile as a function of x at the center along the y-axis ( n y = 32) within the sandwiched system.Despite uniform G i values, notably higher density is observed in the region sandwiched by the high G i regions (within the dotted lines).(b) δρ as a function of G i within the fluctuated system, as shown in Fig.2.Error bars denote standard deviations.It is notable that δρ is larger in regions with the lower G i , while conversely, it is smaller in regions with the higher G i .

Figure 6 .
Figure 6.(a) The total potential energy U as a function of , showcasing various ξ values.Circle, square, triangle, and inverted triangle symbols represent U with ξ = 1, 3, 5, and 10, respectively.(b) The dependence of potential change δU on .Notably, δU exhibits a proportional relationship with 2 for each ξ .(c) The potential change δU as a function of σ , where σ denotes the standard deviation of ρ .It is noteworthy that δU across all systems with varying and ξ values adheres to a scaling law δU ∼ σ 2 .(d) δU is plotted against r , illustrating collapsed behavior across different ξ values when r remains consistent, despite variations in the ξ and combinations.
these stress-induced phenomena.Consequently, material design based on this coarse-grained particle model is considered necessary.