Compression-driven jamming in porous cohesive aggregates

Abstract

I investigate the compression-driven jamming behavior of two-dimensional porous aggregates composed of cohesive, frictionless disks. Three types of initial aggregates are prepared using different aggregation procedures, namely, reaction-limited aggregation (RLA), ballistic particle–cluster aggregation (BPCA), and diffusion-limited aggregation (DLA), to elucidate the influence of aggregate morphology. Using distinct-element-method simulations with a shrinking circular boundary, I numerically obtain the pressure as a function of the packing fraction ϕ. For the densest RLA and the intermediate BPCA aggregates, a clear jamming transition is observed at a critical packing fraction ϕ_J, below which the pressure vanishes and above which a finite pressure emerges; the transition is less distinct for the most porous DLA aggregates. The jamming threshold depends on the initial structure and, when extrapolated to infinite system size, approaches ϕ_J = 0.765 ± 0.004 for RLA, 0.727 ± 0.004 for BPCA, and 0.602 ± 0.023 for DLA, where the errors denote the standard error. Above ϕ_J, the pressure follows P ≈ A(ϕ − ϕ_J)², which implies that the bulk modulus K of jammed aggregates is proportional to ϕ − ϕ_J. Rigid-cluster analysis of jammed aggregates shows that the average coordination number within the largest rigid cluster increases linearly with ϕ − ϕ_J. Taken together, these relations suggest that the elastic response of compressed porous aggregates is analogous to that of random spring networks.

---

1 Introduction

Soft-matter systems consisting of large numbers of discrete particles, such as granular materials, emulsions, foams, and cellular tissues, jam into mechanically rigid solids once their packing fraction, ϕ, exceeds a critical value.¹ This jamming transition has been explored extensively over the past few decades.^2–8 In athermal assemblies of frictionless particles with purely repulsive contacts, macroscopic observables, including pressure, bulk modulus, and shear modulus, display critical behavior in the vicinity of the jamming point ϕ_J.^4,6,9,10

Although much of the foundational work has dealt with purely repulsive, frictionless systems, interparticle cohesion and friction are ubiquitous in real materials. Granular packings that are frictional but otherwise repulsive jam at significantly lower ϕ than their frictionless counterparts.¹¹ The impact of attractive forces on the mechanical response has likewise been examined in a series of studies.^12–17 For example, Koeze et al.¹⁵ showed that the bulk modulus (and hence the pressure) of two-dimensional cohesive, frictionless packings follows a qualitatively different ϕ-dependence from that of purely repulsive systems, a trend later confirmed in three-dimensional systems near jamming by Yoshii and Otsuki.¹⁶

Cohesive forces do more than modify contact mechanics; they also drive aggregation before compression. Under dilute conditions, attractive particles spontaneously form porous clusters. Such aggregates are widespread in nature, ranging from colloidal gels to cosmic dust,^18,19 and their morphology is dictated by the coagulation pathway.²⁰ Most jamming studies of cohesive, frictionless particles, however, start from initially dispersed states and therefore overlook how the pre-compression structure influences the subsequent mechanical response.

Here, I numerically investigate compression-driven jamming in two-dimensional porous aggregates composed of cohesive, frictionless disks. Three types of initial aggregates are prepared using different aggregation procedures, and I perform numerical compression tests using the distinct element method (DEM).²¹ By analyzing the pressure as a function of the packing fraction ϕ, I determine the jamming threshold ϕ_J and its dependence on the initial aggregate morphology. I further identify a critical power-law relation between the pressure and ϕ − ϕ_J above ϕ_J for each aggregate type. These results indicate that the elastic response of compressed porous aggregates is analogous to that of random spring networks.

2 Models

2.1 Particle interaction

I consider systems of monodisperse disk-shaped particles initially confined by a circular wall of radius R_wall and simulate quasi-static compression by gradually reducing R_wall. Each particle has unit mass, m_disk = 1, and unit radius, R_disk = 1. The number of particles is set to N = 1024, 4096, or 16 384 for each simulation. The equation of motion is given bywhere r_i is the position of particle i, and F_i is the total force acting upon it. Rotational degrees of freedom are neglected in this study. The timestep for numerical integration is fixed at Δt = 0.02. For two particles i and j in contact, the interparticle force, F_ij, is given by the sum of the elastic term F_e and dissipative term F_d:

F_ij = F_e + F_d. (2)

I consider a simple dissipation model that is widely used to describe energy dissipation in foam bubbles:²²

where v_i is the velocity of particle i, and the dissipation coefficient for damping (C_d) is set to C_d = 1. In this study, F_d was introduced to facilitate the energy and structural relaxation of the system. Rheological properties, such as viscosity, are beyond the scope of this study, and I expect that the specific choice of F_d does not significantly affect the main conclusions.

The elastic term F_e is a function of the compression length (i.e., penetration depth) between two particles, δ_ij = 2R_disk − |r_i − r_j|. I assume that a contact between particles i and j forms when they collide (δ_ij = 0). A repulsive force acts when the compression length is positive, whereas an attractive force acts when it is negative. In this study, I consider two types of elastic interaction, namely, the unbreakable and the breakable contact models (Fig. 1). For the unbreakable contact model, I do not allow a contact to break once it has formed, and F_e is given by

where δ_ij = 2R_disk − |r_i − r_j| denotes the compression length (i.e., penetration depth) between two particles, and n_ij = (r_i − r_j)/|r_i − r_j| is the normal unit vector. The spring constant (k) is set to k = 10.

Fig. 1 Schematics of elastic interaction models. (a) Unbreakable contact model. F_e = 0 when the two particles are separated (black line). A contact is formed when the particles collide (red filled circle). Once a contact has formed, F_e is given by eqn (4) (red line). A repulsive force acts when the compression length is positive, whereas an attractive force acts when it is negative. In this model, a contact never breaks once it has formed. (b) Breakable contact model. F_e = 0 when the two particles are separated (black line). A contact is formed when the particles collide (red filled circle). Once a contact has formed, F_e is given by eqn (5) (red line). In this model, the contact breaks when δ reaches −2δ_c (black open circle).

In contrast, for the breakable contact model, I consider the weakening and eventual breakage of contacts when δ_ij reaches prescribed threshold values. In this study, I assume a simple model:

and a contact breaks when δ_ij = −2δ_c. The critical length δ_c is set to δ_c = 0.01.

The force exerted on particle i interacting with the wall, F_wall, is

F_wall = F_e,wall + F_d,wall, (6)

F_d,wall = −C_d(v_i − v_wall), (7)

where v_wall = v_walln_wall is the velocity of the wall at the contact point. The speed of the wall, v_wall, is defined as v_wall = −dR_wall/dt. For the unbreakable contact model, F_e,wall is given by

F_e,wall = kδ_walln_wall, (8)

where δ_wall = R_disk − (R_wall − |r_i|) is the compression length with wall, and n_wall = −r_i/|r_i| is the normal unit vector. For the breakable contact model, F_e,wall is given byand a contact breaks when δ_wall = −δ_c. The packing fraction of the system is ϕ = N/R_wall², and the pressure, P, is defined as follows:where the summation is taken over all particle–wall contacts.

In the present study, I perform numerical simulations using the unbreakable contact model to obtain the main results (Section 3). The impact of contact breakage is discussed in Section 4.1.

2.2 Initial aggregate structure

I employ three types of aggregates prepared via distinct procedures as initial configurations: reaction-limited aggregation (RLA), ballistic particle–cluster aggregation (BPCA), and diffusion-limited aggregation (DLA). These aggregates are generated in continuous (off-lattice) space. The corresponding structures are illustrated in Fig. 2(a), (b), and (c), respectively. Below, I briefly describe these aggregation protocols (see also ref. 23).

Fig. 2 Three types of aggregates prepared using distinct aggregation procedures. Panels (a)–(c) show the initial structures of DLA, BPCA, and RLA aggregates with N = 16 384, respectively.

2.2.1 Reaction-limited aggregation (RLA). RLA, also known as the Eden growth model, is a sequential growth process in which particles attach to randomly selected sites on a pre-existing cluster. RLA-type cluster formation typically occurs when the probability of particle attachment upon collision is low, and the overall growth is limited by the attachment kinetics. Such structures are widely observed in nature, ranging from bacterial colonies²⁴ to topological-defect turbulence in nematic liquid crystals.²⁵ The structural properties of RLA clusters have been extensively investigated.^26,27 An initial RLA aggregate with N = 16 384 is shown in Fig. 2(a), and the initial wall radius is set to R_wall,ini = 200.

2.2.2 Ballistic particle–cluster aggregation (BPCA). BPCA is a sequential growth process that occurs via collisions between a particle and a cluster. In this process, a particle launched from a distant point moves ballistically and adheres to the surface of a growing cluster upon collision. BPCA-type cluster formation typically arises in dilute media where the mean free path of particles exceeds the aggregate size. BPCA has been widely used as a model for cosmic dust aggregates^28–30 and atmospheric aerosols.³¹ An initial BPCA aggregate with N = 16 384 is shown in Fig. 2(b), with the initial wall radius set to R_wall,ini = 270.

2.2.3 Diffusion-limited aggregation (DLA). DLA is a sequential growth process driven by collisions between a particle and a cluster. In this process, a randomly walking particle starting far from the origin collides with a growing cluster and adheres to its surface. DLA-type cluster formation is typically observed when growth occurs within a diffusion field.³² Such structures are widely observed in nature, ranging from bacterial colonies³³ to dendritic crystals formed via electrochemical deposition.³⁴ The structural properties of DLA clusters have been extensively investigated.^35,36 An initial DLA aggregate with N = 16 384 is shown in Fig. 2(c), with the initial wall radius set to R_wall,ini = 600.

2.3 Compression and relaxation

In this study, I obtain statically compressed aggregates through a two-step numerical procedure described below (Fig. 3). In the first stage, the wall radius R_wall is continuously decreased, producing a series of configurations C with packing fractions ϕ (Fig. 3(b) and (c)). In the second stage, the generated configurations C(ϕ) are relaxed at constant R_wall to new configurations C_relax(ϕ) (Fig. 3(d) and (e)). The pressure P = P(ϕ) is then measured in the final static configuration.

Fig. 3 Snapshots of a compressed aggregate obtained using a two-step numerical procedure. (a) Initial configuration of BPCA-16384-1. (b) and (c) configurations at ϕ = 0.694 and 0.717 during the first compression stage. (d) and (e) Configurations at ϕ = 0.694 and 0.717 after relaxation. The red particles in panels (d) and (e) form the largest rigid cluster (see Section 3.2).

During the first compression stage, the wall radius R_wall decreases exponentially with time according to a characteristic timescale τ_comp:

R_wall(t) = R_wall,ini exp(−t/τ_comp), (11)

where R_wall,ini is the initial wall radius, and I set τ_comp = 10⁶. In the second relaxation stage, the wall radius is held constant (R_wall = const.), and the aggregate is allowed to relax over a waiting time of τ_relax = 2 × 10⁵.

3 Results

In this study, I prepare three types of initial aggregates (RLA, BPCA, and DLA) for three aggregate sizes (N = 16 384, 4096, and 1024). For each combination of type and size, I generate five aggregates with different initial configurations. Each simulation run is labeled, for example, as RLA-16384-1 or DLA-1024-5.

3.1 Pressure and jamming point

I analyze the pressure as a function of the packing fraction. As P is proportional to k in static configurations, I consider the normalized pressure P/k in this study. Fig. 4 shows the relationship between (P/k)^1/2 and ϕ near the jamming transition. I find that there exists a critical packing fraction ϕ_J such that the pressure vanishes for ϕ < ϕ_J and increases with ϕ for ϕ > ϕ_J. The value of ϕ_J depends on the initial aggregate configuration.

Fig. 4 Normalized pressure P/k after relaxation. The figure shows the relationship between (P/k)^1/2 and ϕ near the jamming transition. (a)–(c) Results for N = 16 384. (d)–(f) Results for N = 4096. (g)–(i) Results for N = 1024. Data with positive P are plotted using filled symbols, whereas data with negative P are plotted as −P using open symbols.

Fig. 4(a)–(c) show the results for RLA, BPCA, and DLA aggregates with N = 16 384, respectively. For RLA and BPCA, I find a linear relationship between (P/k)^1/2 and ϕ for ϕ > ϕ_J. These results indicate that P/k is approximately proportional to (ϕ − ϕ_J)² in these cases. Results for N = 4096 and 1024 are also shown in Fig. 4(d)–(f) and (g)–(i), respectively. The linear relationship between (P/k)^1/2 and ϕ is less clear for some DLA aggregates (e.g., DLA-16384-4), although other DLA aggregates (e.g., DLA-16384-2 and DLA-4096-1) do exhibit an approximately linear trend.

I estimate the value of ϕ_J for each aggregate by least-squares fitting, assuming the following critical scaling relation with exponent α = 2:

The least-squares fits are performed using data in the range 0.005 < (P/k)^1/2 < 0.04 for RLA and BPCA aggregates, and 0.005 < (P/k)^1/2 < 0.02 for DLA aggregates. The fitting parameters for all runs are summarized in Table 1. I acknowledge that, in the immediate vicinity of ϕ_J, i.e., in the region where (P/k)^1/2 < 0.005, the values of P/k exhibit considerable scatter, and thus there is a limit to how close to ϕ_J I can take the fitting range.

Table 1 Fitting parameters (ϕ_J and A) for the packing-fraction dependence of the normalized pressure P/k (eqn (12)). The least-squares fits are performed using data in the range 0.005 < (P/k)^1/2 < 0.04 for RLA and BPCA aggregates, and 0.005 < (P/k)^1/2 < 0.02 for DLA aggregates

Aggregate ID	ϕ J	A
RLA-16384-1	0.7516	1.173
RLA-16384-2	0.7526	1.121
RLA-16384-3	0.7546	1.258
RLA-16384-4	0.7550	1.095
RLA-16384-5	0.7561	1.222
RLA-4096-1	0.7435	1.068
RLA-4096-2	0.7261	0.951
RLA-4096-3	0.7339	0.915
RLA-4096-4	0.7305	1.022
RLA-4096-5	0.7368	0.821
RLA-1024-1	0.7266	0.871
RLA-1024-2	0.7064	0.738
RLA-1024-3	0.7060	0.444
RLA-1024-4	0.7041	0.736
RLA-1024-5	0.7128	0.588
BPCA-16384-1	0.7074	0.659
BPCA-16384-2	0.7218	0.819
BPCA-16384-3	0.7226	0.820
BPCA-16384-4	0.7185	0.708
BPCA-16384-5	0.7163	0.789
BPCA-4096-1	0.7138	0.756
BPCA-4096-2	0.6953	0.481
BPCA-4096-3	0.7058	0.737
BPCA-4096-4	0.7091	0.886
BPCA-4096-5	0.7044	0.660
BPCA-1024-1	0.6739	0.398
BPCA-1024-2	0.6885	0.548
BPCA-1024-3	0.6987	0.698
BPCA-1024-4	0.6774	0.542
BPCA-1024-5	0.6937	0.544
DLA-16384-1	0.6119	0.148
DLA-16384-2	0.5824	0.055
DLA-16384-3	0.5630	0.057
DLA-16384-4	0.6054	0.126
DLA-16384-5	0.5773	0.055
DLA-4096-1	0.6207	0.261
DLA-4096-2	0.5814	0.101
DLA-4096-3	0.5938	0.086
DLA-4096-4	0.5183	0.063
DLA-4096-5	0.6186	0.119
DLA-1024-1	0.6201	0.204
DLA-1024-2	0.5846	0.084
DLA-1024-3	0.5803	0.074
DLA-1024-4	0.4607	0.050
DLA-1024-5	0.5465	0.148

---

To elucidate the effect of finite system size N on the jamming point ϕ_J, I investigate how ϕ_J depends on N for the three types of aggregates. Fig. 5 shows the system-size dependence of ϕ_J. I find that ϕ_J is approximately given by the following empirical relation:

ϕ_J = ϕ_J,∞ − cN^−1/2, (13)

where ϕ_J,∞ is the jamming point for an infinitely large aggregate and c is a fitting constant. From the fitting shown in Fig. 5, I obtain ϕ_J,∞ = 0.765 ± 0.004 for RLA, 0.727 ± 0.004 for BPCA, and 0.602 ± 0.023 for DLA, where the quoted errors denote the standard errors. These results highlight that the initial aggregate morphology strongly influences the jamming point and underscore the importance of the pre-compression aggregation history in cohesive granular materials.

Fig. 5 System-size dependence of the jamming point ϕ_J. Symbol shapes (circle, triangle, diamond, square, and pentagon) correspond to those in Fig. 4. The solid lines represent the best fits obtained by least-squares fitting to eqn (13).

I briefly discuss the bulk modulus of the compressed aggregates. The bulk modulus K is defined as

Given that the pressure follows the critical scaling described by eqn (12), the bulk modulus is expressed as

Eqn (15) implies that K varies continuously across ϕ_J.

The dependence of pressure on the filling factor in cohesive granular materials differs qualitatively from that in purely repulsive granular systems, where F_e is defined as

The elastic properties of granular systems composed of frictionless particles with linear repulsive interactions have been extensively studied in previous works.^3,4 It is known that the pressure scales as P/k = A(ϕ − ϕ_J) for ϕ > ϕ_J,⁴ resulting in a discontinuous jump of the bulk modulus K at ϕ_J. For granular systems composed of frictionless Hertzian repulsive particles, one instead finds P/k = A(ϕ − ϕ_J)^3/2 for ϕ > ϕ_J.⁴ In Section 3.2, I perform a rigid-cluster analysis to understand the origin of this difference.

In my simulations, a linear relationship between (P/k)^1/2 and ϕ − ϕ_J holds with high accuracy for RLA and BPCA with N = 16 384 (Fig. 4(a) and (b)). Therefore, in the analyses presented in Sections 3.2 and 3.3, I focus on these results.

3.2 Rigid cluster analysis

I perform a rigid-cluster analysis to identify the largest rigid cluster in each final static configuration. Rigid clusters are defined as those whose zero-frequency vibrational modes correspond only to rigid-body motions, and I use the pebble-game algorithm³⁷ to identify them. The red particles in Fig. 3(d) and (e) form the largest rigid cluster. An aggregate confined by a boundary wall enters a jammed state when a rigid cluster has more than two contact points with the boundary wall. This rigidity percolation is evident in Fig. 3(e); the largest rigid cluster has many contacts with the wall for ϕ > ϕ_J.

I define N_cluster as the number of particles that constitute the largest rigid cluster, and the fraction of the largest rigid cluster as f_cluster = N_cluster/N. Fig. 6(a) and (b) show 1 − f_cluster as a function of ϕ − ϕ_J. Results for the RLA and BPCA aggregates with N = 16 384 are shown in Fig. 6(a) and (b), respectively. I find that 1 − f_cluster decreases sharply at ϕ_J in both cases, and that it can be approximated by

for ϕ > ϕ_J. The dashed lines in Fig. 6(a) and (b) represent the fits obtained using data in the range 0.01 < ϕ − ϕ_J < 0.06. The fitting parameters (f₀ and s) are summarized in Table 2.

Fig. 6 (a) and (b) Fraction of the largest rigid cluster, f_cluster = N_cluster/N, as a function of ϕ − ϕ_J. Panels (a) and (b) show the results for the RLA and BPCA aggregates with N = 16 384, respectively. Symbol shapes (circle, triangle, diamond, square, and pentagon) correspond to those in Fig. 4. Filled symbols indicate data for ϕ ≥ ϕ_J, whereas open symbols indicate data for ϕ < ϕ_J. The dashed lines represent the fits obtained using data in the range 0.01 < ϕ − ϕ_J < 0.06 (eqn (17)). (c) and (d) Coordination number within the largest rigid cluster, Z_eff, as a function of ϕ − ϕ_J. Panels (c) and (d) show the results for the RLA and BPCA aggregates with N = 16 384, respectively. (e) and (f) Coordination number of the whole aggregate, Z, as a function of ϕ − ϕ_J. Panels (e) and (f) show the results for the RLA and BPCA aggregates with N = 16 384, respectively.

Table 2 Fitting parameters (f₀ and s) for the packing-fraction dependence of the fraction of the largest rigid cluster (eqn (17)). The least-squares fits are performed using data in the range 0.01 < ϕ − ϕ_J < 0.06

Aggregate ID	f 0	s
RLA-16384-1	0.12	0.031
RLA-16384-2	0.11	0.034
RLA-16384-3	0.12	0.029
RLA-16384-4	0.10	0.034
RLA-16384-5	0.08	0.039
BPCA-16384-1	0.27	0.037
BPCA-16384-2	0.17	0.040
BPCA-16384-3	0.16	0.040
BPCA-16384-4	0.21	0.042
BPCA-16384-5	0.19	0.038

---

I then calculate the coordination number within the largest rigid cluster, Z_eff (Fig. 6(c) and (d)). The largest rigid cluster is jammed at ϕ_J, and my numerical results show that Z_eff − Z_iso ≈ 0 at the jamming point, where Z_iso = 4 is the isostatic coordination number for two-dimensional aggregates. I also find that Z_eff − Z_iso increases linearly with ϕ − ϕ_J for both RLA and BPCA aggregates:

Z_eff − Z_iso ≈ m(ϕ − ϕ_J), (18)

where m is the slope in Fig. 6(c) and (d). For linear repulsive systems, Z_eff − Z_iso is known to scale as Z_eff − Z_iso = m(ϕ − ϕ_J)^1/2.⁴ Thus, the linear relationship shown in Fig. 6(c) and (d) is a distinctive feature of cohesive granular systems.

For aggregates of linear cohesive particles, K is approximately proportional to ϕ − ϕ_J (eqn (15)). In addition, since ϕ − ϕ_J increases linearly with Z_eff − Z_iso, K is proportional to Z_eff − Z_iso. This linear relationship is consistent with that found for random spring networks.⁶ Ellenbroek et al.⁶ investigated the elastic response of two types of spring networks: one derived from real packings of linear repulsive particles, and another obtained by randomly cutting bonds in a highly connected network derived from a jammed packing. They showed that K is proportional to Z_eff − Z_iso and remains continuous at Z_iso for random spring networks, whereas a discontinuous jump of K at Z_iso appears for networks derived from repulsive particle packings.

Although both Ellenbroek et al.⁶ and I investigate the elastic response of random spring networks, the procedures used to construct these networks differ. As described above, they reduced the coordination number Z by randomly cutting springs in the contact network obtained from a jammed packing. In contrast, in the present study I increase Z by compressing a system of cohesive particles and thereby evolving the contact network. It is intriguing that similar critical scaling is observed in both cases, despite the fact that the evolution of Z proceeds in opposite directions.

3.3 Interparticle compression length

After relaxation, the force acting on a particle is determined by the elastic term (see eqn (4)), which is proportional to the interparticle compression length. Here, I define the mean and root-mean-square of the interparticle compression length δ as δ_mean and δ_rms, respectively. Fig. 7(a) and (b) show δ_mean and δ_rms as functions of ϕ − ϕ_J for RLA and BPCA aggregates with N = 16 384, respectively. The dashed and solid lines represent the fits obtained via least-squares fitting, and I find that the following relations describe their dependence on the packing fraction:

δ_mean = B(ϕ − ϕ_J)^β, (19)

δ_rms = C(ϕ − ϕ_J)^γ. (20)

Fig. 7 δ _mean and δ_rms as functions of ϕ − ϕ_J. Panels (a) and (b) show the results for the RLA and BPCA aggregates with N = 16 384, respectively. Filled and open symbols correspond to δ_mean and δ_rms, respectively. The dashed and solid lines represent the fits obtained via least-squares fitting (eqn (19) and (20)).

I obtain β = 2.20, B = 9.03, γ = 1.61, and C = 2.90 for the RLA case, and β = 2.00, B = 3.06, γ = 1.39, and C = 1.13 for the BPCA case. Notably, β is close to 2 in both cases, indicating that δ_mean is approximately proportional to P/k. In other words, the mean interparticle force F_mean = (k/2)δ_mean scales linearly with the macroscopic pressure P. This proportionality is consistent with the Rumpf equation:³⁸

where Z is the coordination number of the aggregates. The dependence of Z on the filling factor is shown in Fig. 6(e) and (f). For ϕ ≈ ϕ_J, Z ≈ 4 and therefore F_mean ≈ (π/ϕ_J)P.

I note that eqn (19) and (20) are valid only in the regime where δ_mean/δ_rms ≪ 1. Fig. 7 indicates that the ratio δ_mean/δ_rms increases with ϕ, because δ_mean and δ_rms are governed by different exponents, β and γ. However, by definition, the mean must be less than or equal to the root mean square; thus, δ_mean/δ_rms should always be smaller than 1.

Since δ can take both positive and negative values, the root-mean-square value, δ_rms, characterizes the typical amplitude of δ in the system. The scaling of δ_rms with ϕ − ϕ_J (eqn (20)) is therefore expected to be related to the geometric properties of the contact network, although the connection remains unclear at present.

4 Discussion

4.1 Breakable contact model

In this study, I mainly investigate the compressive behavior of aggregates whose constituent particles interact via unbreakable contacts (Fig. 1(a)). Although this idealized and simple setting is useful for capturing the essence of cohesive systems, the unbreakable contact model may be unrealistic in some situations. I therefore discuss the impact of contact breakage on the compressive behavior in this section.

Fig. 8 shows the normalized pressure as a function of ϕ. Filled circles indicate the results for the unbreakable contact model, whereas open squares indicate those for the breakable contact model with δ_c = 0.01. The pre-compression aggregate is DLA-4096-1 for both contact models. It is clear that the pressure is essentially independent of the choice of contact model around the jamming point, whereas the two curves start to deviate from each other when ϕ exceeds a certain value (in this case, ϕ > 0.633). The critical scaling around the jamming point, P/k = A(ϕ − ϕ_J)², is therefore a common feature of both the unbreakable and breakable contact models.

Fig. 8 Normalized pressure P/k as a function of ϕ. Filled circles indicate the results for the unbreakable contact model, whereas open squares indicate those for the breakable contact model with δ_c = 0.01. The pre-compression aggregate is DLA-4096-1 for both contact models.

When δ ≥ −δ_c holds for all interparticle contacts within the aggregate, F_e in the breakable contact model exactly coincides with that in the unbreakable contact model. In contrast, if δ < −δ_c for some contacts, the values of F_e for those contacts differ between the two models. Fig. 9 shows the frequency distribution of δ, p(δ). Here, I define p(δ) such that it satisfies the normalization condition . Fig. 9(a)–(c) correspond to the unbreakable contact model, and Fig. 9(d)–(f) correspond to the breakable contact model. For the unbreakable contact model, the distribution is approximately symmetric with respect to δ = 0. At ϕ = 0.628, the distribution for the breakable contact model (Fig. 9(d)) is identical to that for the unbreakable contact model (Fig. 9(a)). This is consistent with the fact that the pressure at ϕ = 0.628 is independent of the contact model.

Fig. 9 Frequency distribution of δ, p(δ). Panels (a)–(c) show the results for the unbreakable contact model, and panels (d)–(f) show the results for the breakable contact model. Vertical dashed lines at δ = −δ_c indicate the critical compression length in the breakable contact model (eqn (5)).

The distributions p(δ) for the two models deviate from each other at ϕ = 0.636 and 0.646. In these cases, a fraction of the contacts in the unbreakable contact model (Fig. 9(b) and (c)) have δ < −δ_c. In contrast, in the breakable contact model (Fig. 9(e) and (f)), no contacts with δ < −δ_c are present. Although interparticle contacts are broken at δ = −2δ_c in the breakable contact model (Fig. 1(b)), contacts with −2δ_c < δ < −δ_c are mechanically unstable and are not expected to exist in static configurations.¹⁶

4.2 Compression of randomly dispersed cohesive particles

In this study, I investigated the compression behavior of aggregates using the DEM and determined the critical scaling of the pressure (eqn (12)). Although ϕ_J varies depending on the initial structure, the exponent α = 2 appears to be common. Furthermore, my analysis of Z_eff and a comparison with random spring networks⁶ suggest that this scaling reflects the properties of the interparticle interaction model (Section 3.2).

This naturally raises the question: how does the compression behavior change if the particles are initially distributed randomly rather than forming aggregate structures? To address this question, I perform additional simulations. Using the unbreakable contact model, I carry out compression simulations of a system of randomly dispersed cohesive particles with N = 4096 and an initial wall radius R_wall,ini = 128 (corresponding to an initial packing fraction of 0.25). The initial configuration is shown in Fig. 10(a), and the configuration at ϕ = 0.666 is shown in Fig. 10(b).

Fig. 10 (a) Initial configuration of randomly dispersed cohesive particles. (b) Snapshot of a compressed system at ϕ = 0.666. (c) Normalized pressure P/k as a function of ϕ. Data with positive P are plotted using filled symbols, whereas data with negative P are plotted as −P using open symbols.

Fig. 10(c) shows the normalized pressure P/k after relaxation as a function of ϕ. I find that a jamming transition occurs at ϕ_J ≈ 0.63. As in the case of aggregates, the data near the jamming transition are approximately described by P/k = A(ϕ − ϕ_J)². This result further supports the idea that the critical scaling with exponent α = 2 originates from random networks of linear springs.

It should be noted, however, that this result for the randomly dispersed system differs from the conclusions of earlier studies that carried out simulations under similar conditions.^14,15 Koeze et al.¹⁵ performed compaction simulations of randomly dispersed cohesive particles and investigated the relationship between the bulk modulus K and ϕ. The interparticle interaction they used is such that two particles i and j interact when δ_ij > −2δ_c, and the corresponding force is given by eqn (5). In contrast to the breakable contact model adopted in the present study, they did not take into account any hysteresis between connected and separated states. In addition, they employed periodic boundary conditions and compressed the system by repeatedly applying an affine transformation followed by relaxation. They found a critical scaling with exponent α ≈ 3.5 and reported that ϕ_J depends on δ_c. The differences between their results and those of the present study indicate that both the exponent α and the jamming point ϕ_J depend on the details of the interparticle attraction model and/or on the compression protocol.

5 Summary

In this study, I investigated the compression-driven jamming behavior of two-dimensional porous aggregates composed of cohesive, frictionless disks. Three types of initial aggregates were prepared using different aggregation procedures (RLA, BPCA, and DLA) in order to clarify how the initial morphology affects the jamming behavior. I showed that the system exhibits a clear jamming transition at a critical packing fraction ϕ_J that depends strongly on the initial structure. By extrapolating to infinite system size, I obtained ϕ_J = 0.765 ± 0.004 for RLA aggregates, 0.727 ± 0.004 for BPCA aggregates, and 0.602 ± 0.023 for DLA aggregates. Above the transition, the pressure follows a power-law scaling, P = A(ϕ − ϕ_J)², which implies that the bulk modulus K is approximately proportional to ϕ − ϕ_J.

To understand the microscopic origin of this behavior, I performed a rigid-cluster analysis. The jamming transition is accompanied by the percolation of the largest rigid cluster that makes multiple contacts with the confining wall. At the jamming point, the effective coordination number within the largest rigid cluster satisfies Z_eff ≈ Z_iso = 4, and for ϕ > ϕ_J I found that Z_eff − Z_iso increases linearly with ϕ − ϕ_J (Fig. 6(c) and (d)). Consequently, the bulk modulus is proportional to Z_eff − Z_iso, in contrast to frictionless systems with purely repulsive, linear interactions, for which Z_eff − Z_iso is known to scale as (ϕ − ϕ_J)^1/2.⁴ This linear trend implies that the elastic response of compressed porous aggregates is analogous to that of random spring networks.⁶

I also examined the effect of contact breakage by comparing an unbreakable contact model with a breakable contact model, in which contacts are weakened and eventually broken under finite tensile loading (Fig. 1(b)). The results show that the critical scaling of the pressure near the jamming point, P = A(ϕ − ϕ_J)², is essentially identical for both models, and that differences appear only at higher packing fractions where tensile contacts become unstable in the breakable model (Fig. 8).

In summary, this work demonstrated that the jamming landscape of cohesive soft materials is strongly modulated by the aggregate formation pathway, thereby bridging concepts from colloidal gelation, porous-media mechanics, and granular physics. Future extensions to three-dimensional systems, frictional interactions, and shear- or tensile-driven dynamics will be invaluable for establishing a comprehensive, morphology-aware jamming framework relevant to a broad range of soft-matter physical contexts.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data for this article are available at Zenodo (https://doi.org/10.5281/zenodo.15722246).

Acknowledgements

This research is supported by JSPS KAKENHI Grant (JP24K17118, JP25K00025, and JP25K01063). Numerical computations were in part carried out on the general-purpose PC cluster at the Center for Computational Astrophysics, National Astronomical Observatory of Japan.

References

1. M. van Hecke, J. Phys.: Condens. Matter, 2010, 22, 033101.
2. H. A. Makse, D. L. Johnson and L. M. Schwartz, Phys. Rev. Lett., 2000, 84, 4160–4163.
3. C. S. O'Hern, S. A. Langer, A. J. Liu and S. R. Nagel, Phys. Rev. Lett., 2002, 88, 075507.
4. C. O'Hern, L. Silbert, A. Liu and S. Nagel, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2003, 68, 011306.
5. T. S. Majmudar, M. Sperl, S. Luding and R. P. Behringer, Phys. Rev. Lett., 2007, 98, 058001.
6. W. G. Ellenbroek, Z. Zeravcic, W. van Saarloos and M. van Hecke, EPL, 2009, 87, 34004.
7. A. J. Liu and S. R. Nagel, Annu. Rev. Condens. Matter Phys., 2010, 1, 347–369.
8. G. Parisi and F. Zamponi, Rev. Mod. Phys., 2010, 82, 789–845.
9. A. Zaccone and E. Scossa-Romano, Phys. Rev. B:Condens. Matter Mater. Phys., 2011, 83, 184205.
10. A. Zaccone, J. Appl. Phys., 2025, 137, 050901.
11. A. P. Santos, D. S. Bolintineanu, G. S. Grest, J. B. Lechman, S. J. Plimpton, I. Srivastava and L. E. Silbert, Phys. Rev. E, 2020, 102, 032903.
12. G. Lois, J. Blawzdziewicz and C. S. O'Hern, Phys. Rev. Lett., 2008, 100, 028001.
13. L. M. Lopatina, C. J. Olson Reichhardt and C. Reichhardt, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2011, 84, 011303.
14. D. J. Koeze and B. P. Tighe, Phys. Rev. Lett., 2018, 121, 188002.
15. D. J. Koeze, L. Hong, A. Kumar and B. P. Tighe, Phys. Rev. Res., 2020, 2, 032047.
16. K. Yoshii and M. Otsuki, J. Phys. Soc. Jpn., 2025, 94, 084801.
17. A. T. Grigas, A. Fisher, M. D. Shattuck and C. S. O'Hern, Phys. Rev. E, 2024, 109, 034406.
18. C. P. Johnson, X. Li and B. E. Logan, Environ. Sci. Technol., 1996, 30, 1911–1918.
19. J. Blum and G. Wurm, Annu. Rev. Astron. Astrophys., 2008, 46, 21–56.
20. P. Meakin, Rev. Geophys., 1991, 29, 317–354.
21. P. A. Cundall and O. D. L. Strack, Geotechnique, 1979, 29, 47–65.
22. D. J. Durian, Phys. Rev. E:Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1997, 55, 1739–1751.
23. P. Meakin, J. Sol-Gel Sci. Technol., 1999, 15, 97–117.
24. S. Kaczmarczyk, F. Koza, D. Śnieżek and M. Matyka, Phys. Rev. E, 2025, 111, 014406.
25. K. A. Takeuchi and M. Sano, J. Stat. Phys., 2012, 147, 853–890.
26. K. A. Takeuchi, J. Stat. Mech.:Theory Exp., 2012, 2012, 05007.
27. N. Kobayashi and H. Yamazaki, J. Phys. Soc. Jpn., 2018, 87, 014005.
28. K. Wada, H. Tanaka, T. Suyama, H. Kimura and T. Yamamoto, Astrophys. J., 2009, 702, 1490–1501.
29. S. Arakawa, H. Tanaka, E. Kokubo, D. Nishiura and M. Furuichi, Astron. Astrophys., 2023, 670, L21.
30. S. Krijt, S. Arakawa, M. Oosterloo and H. Tanaka, Mon. Not. R. Astron. Soc., 2024, 534, 2125–2133.
31. D. W. Schaefer and A. J. Hurd, Aerosol Sci. Technol., 1990, 12, 876–890.
32. T. A. Witten, Jr. and L. M. Sander, Phys. Rev. Lett., 1981, 47, 1400–1403.
33. M. Matsushita and H. Fujikawa, Phys. A, 1990, 168, 498–506.
34. D. Grier, E. Ben-Jacob, R. Clarke and L. M. Sander, Phys. Rev. Lett., 1986, 56, 1264–1267.
35. S. Ohta, J. Phys. Soc. Jpn., 2009, 78, 013605.
36. S. Arakawa, J. Phys. Soc. Jpn., 2024, 93, 024401.
37. D. J. Jacobs and M. F. Thorpe, Phys. Rev. Lett., 1995, 75, 4051–4054.
38. H. C. H. Rumpf, Chem. Ing. Tech., 1970, 42, 538–540.

---