Structural aging of a cohesive and amorphous granular solid under cyclic loading

Abstract

We investigate how cyclic loading evolves the structure and deformation behavior of a granular raft composed of particles floating at an air-oil interface. The raft has a disordered particle packing structure and is cohesive due to capillary interactions between particles. Under uniaxial cyclic loading with a small strain amplitude, the raft's packing structure experiences an aging process characterized by logarithmically increasing packing fraction and decreasing structural heterogeneity. The observed structural change is due to particle dynamics that are organized around morphologically evolving voids in the raft. The raft is then subjected to quasi-static tension or compression tests until failure. In comparison with non-aged rafts, the rafts that experienced cyclic loading show a higher strength, higher stiffness, and lower ductility, along with qualitatively different features, such as a stress overshoot in the loading curve.

---

1 Introduction

The deformation of amorphous solids is an important process throughout numerous industries, geophysical processes, and everyday life.^1–3 Processes like squeezing toothpaste, spreading mayonnaise, and scrubbing a foaming soap exemplify ductile deformation, where a material's constituent particles continuously rearrange without system-spanning fracture. Conversely, shattering glass, breaking chalk, and large-scale landslides exemplify brittle failure, characterized by the formation of shear bands and a steep drop in the loading curve.^4–13

Tuning the ductility of amorphous solids is important for many engineering applications. A common approach is thermal annealing, which uses sequential temperature changes to manipulate a material's structure and its associated energy landscape. Past studies have uncovered the roles of the annealing temperature, duration, and quenching rate on the ductility of various materials.^6,8,14–20 As the interplay between structure, elasticity, and plasticity is complex in amorphous solids,^2,21,22 a change in ductility is often accompanied by concurrent changes in other deformation parameters, such as the elastic modulus and yield strength.^{4–8,12,13,23} However, this complex interplay also enables the manipulation of a disordered structure through mechanical deformation.

An amorphous solid can be mechanically annealed toward lower energy states via cyclic strain.^14,24 The lower energy serves to enhance a material's stability and generally decrease its ductility.^4–6,8,14 Numerical simulations have identified several important factors that influence the behaviors of amorphous solids during cyclic shear.^{9,11,14,25–27} For example, large shearing amplitudes may induce shear bands,^9,23 while small amplitudes can lead the particle system to a localized energy minimum, characterized by reversible particle trajectories.^25,28–31 However, studies on this topic often rely on simulations of idealized particle systems, which discount important complexities that stem from realistic particle interactions such as 3-body effects. Additionally, many of these works focus on pure or simple shear, but do not address alternative or mixed loading modes that induce volumetric strains.

While experiments measuring the change of bulk behaviors of amorphous solids such as metallic glasses have been informative,^32,33 obtaining the necessary particle-scale information for a mechanistic understanding is challenging. This renders experiments with macroscopic particles useful, such as granular particles and colloids, as the particle packing structure and dynamics during cyclic loading can be directly measured.^30,34–40 Existing experimental work with macroscopic particles has already revealed novel features in cyclic deformation, such as a perpetual creeping motion with free or stress-controlled boundaries, potentially resulting from complexities in the energy landscape.^38,40–42 However, most existing experiments on this topic used particles with steric repulsive interactions, which fail to consider more complicated effects such as long-range attractions, and possible many-body effects that are generally associated with molecular and atomic solids.²⁷ As such, we emphasize the need to understand the structure and deformation behaviors stemming from realistic particle interactions and mixed loading modes using macroscopic particles that can be directly observed.

We use a model experimental system that consists of a disordered and two-dimensional particle raft floating at an air–liquid interface.^22,43–47 The particle interaction includes a contact-based short range repulsion, friction, and a long-range capillary attraction resulting from buoyancy and surface tension, which has a functional form close to that of atomic interactions.⁴⁸ The attraction may also have many-body effects if the liquid surface distortion due to the presence of many nearby particles cannot be superposed.⁴⁹ Carrying over from our previous study, we configured a uniaxial tension test with the cohesive raft having two free boundaries.^22,43 This offers a loading mode that contains both volumetric and deviatoric strains, and permits changes to the particle packing fraction.

In this work, we measure the dynamics and structural evolution in a cohesive granular raft subject to cyclic uniaxial loading protocols. Different from previous cyclic shear protocols that are volume-preserving, we observe a logarithmic increase in the particle packing fraction and a decrease in the structural heterogeneity in the raft, which resemble long-term aging in thermal systems. We connect these overall structural changes to microscopic particle dynamics that are organized around void collapsing behavior during cyclic loading. Finally, we discuss how the structural changes in the cyclic loading impact the mechanical behaviors of the raft when it is subject to subsequent tension or compression deformation until failure.

2 Experimental methods

Using a uniaxial deformation apparatus,^43,50–52 we tested a granular raft composed of a monolayer of polydisperse spherical particles floating at an air–liquid interface, illustrated in Fig. 1a. The raft is cohesive due to many-body capillary attractions between particles. The particles are made of closed-cell Styrofoam, with a density of approximately 15 kg m⁻³ and polydisperse diameters within the range d = 1.0 ± 0.1 mm.⁴³ The liquid is a mineral oil, which has a density of ρ = 870 ± 10 kg m⁻³ and a kinematic viscosity of ν = 13.5 cSt. Its surface tension is γ = 27.4 ± 0.7 dyn cm⁻¹, corresponding to a capillary length of mm. This sets a small bond number of Bo = d²/4l_c² ≈ 0.08, at which l_c approximates the interaction range of the capillary attraction between nearby particles.^48,49 The elasto-capillary number, γ/Yd, is on the order of 10⁻⁵ since the Young's modulus (Y) of Styrofoam is on the order of 1 MPa, which means that particles do not deform significantly due to capillary or applied forces.

Fig. 1 (a) Schematic of the experimental setup with an image of the particle packing. (b) A schematic demonstrating the procedure for a single loading cycle. (c) Images demonstrating a tensile loading test. (d) Images demonstrating a compression loading test.

The granular raft has a rectangular shape with height L₀ = 120d and width W₀ = 60d, corresponding to roughly 7000 particles. The two shorter sides of the rectangle are bonded to floating boundaries made of hollow carbon fiber tubes, visualized in Fig. 1a. The two longer sides are free, however, so as not to induce resistance to motion. To generate the initial packing, we first manually laid a thin column of particles connecting the two floating boundaries. Then we sprinkled more particles near this column, prompting the particles to join the columnar assembly via capillary attraction. The generated rafts are dense and strictly monolayer.

For the loading experiments, we fixed the position of the bottom boundary and used two soft and thin beams to drive the top boundary with controlled displacement δ. We measured the force, F, exerted on the top boundary by the particles by monitoring the small deflection of the two beams, as described in our previous work.⁴³ The speed of the moving boundary is δ = 2 µm s⁻¹, setting a capillary number, Ca = νρδ/γ ≈ 1 × 10⁻⁶, meaning that the viscous drag on the particles is negligible in comparison to the surface tension-induced particle attraction. The strain rate under this speed is δ/L₀ = 1.7 × 10⁻⁵ s⁻¹. We can therefore assume that the induced deformation is quasi-static.

Two types of experiments were conducted. The first is a uniaxial small amplitude cyclic deformation experiment driven by the displacement-controlled boundary. Each deformation cycle begins at δ = 0, after which the raft is extended to δ = δ_A. It is then compressed until δ = −δ_A, before finally returning to δ = 0, as shown in Fig. 1b. The number of loading cycles, N, was set to 1000 for observing long-term behavior as required for understanding the aging of amorphous solids. Because performing thousands of loading cycles can be time consuming, we increased the driving speed to . A cyclic strain amplitude of ε_A = δ_A/L₀ = 0.3% is selected, the choice of which will be discussed in Section 3.

In the second experiment, the raft is subjected to uniaxial tension (δ_t > 0) or compression (δ_c < 0) until it fails, as illustrated in Fig. 1c and d, respectively. The strain is again induced by the displacement-controlled boundary, under assumed quasi-static motion. To specifically understand the influence of cyclic loading on the mechanical behavior of the raft, tension and compression tests were performed both on the rafts directly following construction and on the “aged” rafts that had experienced cyclic loading with N = 1000. In addition to monitoring the force F, we imaged and tracked the trajectories of all particles in the system, while only using particles that are at least five particle diameters away from the boundaries for analysis. Particle centers are tracked using an image analysis algorithm with sub-pixel accuracy.^43,50

3 Effect of strain amplitude and number of cycles

Our preliminary results suggest that for ε_A > 1%, cyclic loading rapidly leads to the formation of a system-spanning shear band, which develops into a fracture at N ≈ 100, see an experimental snapshot in Fig. 2. For ε_A ≤ 1%, however, no such shear band forms, allowing the structure of the raft to evolve without global failure. The global strain of 1% roughly coincides with the yield strain of the rafts during tensile tests as reported in our previous study⁴³ and in Section 7 of this study, placing our choice of ε_A = 0.3% in the quasi-elastic regime. This choice is further strengthened by another preliminary experiment in which a cyclic loading amplitude of ε_A = 0.3% was performed up to N = 5000, where no indications of failure were observed in the structure of the raft. The comparison between N = 100 and N = 5000 for ε_A = 0.3% is also shown in Fig. 2.

Fig. 2 Snapshots of granular rafts subjected to a variety of loading amplitudes at several different loading cycle numbers.

4 Global dynamics during cyclic loading

To quantify the degree of structural change between consecutive cycles, we use the non-affine displacement, D²_min, calculated for each particle i by subtracting a best affine fit of its surrounding displacement field from the actual displacements of its neighbors j.⁵³ The formula iswith r_ji = r_j − r_i as the relative position vector between particles i and j, n_j is the number of neighbors within 1.25d of i, and E is the best fit affine transformation matrix.^53,54 We consider the time interval Δt to be one cycle and time t to indicate the beginning of each cycle. The system and ensemble average, 〈D²_min〉, decreases nearly logarithmically as a function of the loading cycle N, see Fig. 3a, indicating a slow-down in dynamics.

Fig. 3 Dynamic and structural quantities during cyclic loading with strain amplitude ε_A = 0.3% over the course of N = 1000 cycles. (a) The system-averaged D²_min vs. the cycle number. (b) Normalized particle mean-squared displacement vs. the cycle number. The orange curve is the original MSD and the blue is the cage-relative MSD. Slope of one provided for reference. (c) The global packing fraction vs. the cycle number. (d) The system-averaged bond order parameter Ψ₆ vs. the cycle number. (e) Distribution of Q_k for N = 0 (yellow) and N = 1000 (blue). The solid curve is a fit to the Gaussian distribution. Inset demonstrates the Voronoi cell anisotropy vectors for a single Delaunay triangle. (f) The standard deviation of Q_k vs. the cycle number. In all panels, the data points and shades respectively represent the mean and the test-to-test fluctuations from 19 repetitions.

To quantify longer-term particle dynamics, we calculate the mean-squared displacement (MSD), 〈Δr_i²(t)〉 = 〈|r_i(t) − r_i(0)|²〉, with respect to N = 0. As shown in Fig. 3b, the normalized MSD has two regimes: a super-diffusive regime with a slope above 1 for N < 10 and a sub-diffusive regime with a slope below 1 for N ≥ 10.

To better understand the observed super-diffusivity, we calculate a cage-relative MSD by subtracting the mean displacement of a particle's neighborhood from the displacement of that particle, written as

This eliminates the effect of rigid body motion and leaves only the local affine and non-affine displacements. The result in Fig. 3b shows a consistent sub-diffusive behavior for all shear cycles. For N < 10, the super-diffusivity in the original MSD is due to the rigid body motion of groups of particles. These groups may reside near locations of intense particle rearrangements, possibly originating from residual stresses due to the initial raft construction.

As previously noted, the structural evolution of our rafts distinguishes itself from the results of isochoric cyclic shearing simulations, as the presence of the free boundaries permits density change. The density is tracked through the packing fraction ϕ calculated from a Delaunay triangulation of the particle centers. For each triangle k, we calculate its area A_k and the sum of the circular sector areas, A^k_sc, due to the overlap area between k and its three constituent particles (highlighted in green in Fig. 4b).

Fig. 4 Void identification and void area calculation. (a) A binarized image showcasing an identified void. (b) The original experimental image corresponding to (a), overlaid with its Delaunay triangles highlighted in red. An example triangle k is highlighted with black borders, and the associated particles are highlighted in orange. The circular sectors of the particles are highlighted in green. The insets illustrate the areas A^k_bp, A^k_sc, and A_k that are calculated for each triangle k. (c) Distributions of void areas based on the number of constituent particles making up the boundary of a void across all experiments taken at N = 1000. (d) System-averaged void size across all experiments against the cycle number. The shaded region indicates standard error across repeated experiments.

The packing fraction is then calculated as . Over the course of 1000 loading cycles, Fig. 3c shows that ϕ consistently increases beyond the first few loading cycles, with a clear long-term logarithmic trend emerging for N > 100.

To detect higher-order structural changes accompanying the densification, we characterize local ordering through the widely used bond orientation parameter for each particle i,

with θ_j being the angle between the vector r_ji and a chosen reference direction, and n_j denoting the first shell of neighbors from the triangulation.

Under this definition, a hexagonally ordered arrangement of neighbors around particle i results in Ψ_6,i ≈ 1, which decreases with increasing disorder. Similar to the packing fraction, 〈Ψ₆〉 gradually increases, with a logarithmic trend for N > 100, as shown in Fig. 3d, indicating increased ordering in the raft.

We then characterize the packing structure's spatial fluctuations based on the Delaunay triangulation and the corresponding Voronoi tessellation. Specifically, we calculate the weighted divergence of the Voronoi cell anisotropy vector,^43,55 Q_k, which can be calculated for each Delaunay triangle k as

Q_k ≡ (∇·c)(A_k/〈A〉), (4)

where c is the Voronoi cell anisotropy vector pointing from a particle's center to the centroid of its Voronoi cell, as illustrated in the inset of Fig. 3e. The divergence in triangle k is evaluated based on the c vectors of the three constituent particles,⁵⁵ which is further weighted by the ratio between the triangle area A_k and the average triangle area 〈A〉 to ensure a zero mean. Under this definition, a triangle covering a tightly packed region should have its c vectors pointing outward, resulting in Q_k > 0, and a triangle covering a loosely packed region should have c pointing inward and Q_k < 0.

The distributions of Q_k for N = 0 and N = 1000 are shown in Fig. 3e, and both are Gaussian-like around Q_k = 0. For larger |Q_k|, the distributions have tails that trend above a Gaussian fit, signaling the existence of regions with highly heterogeneous local packing fractions in the raft. When comparing the distributions of Q_k from N = 0 and N = 1000, a decrease in the standard deviation, std(Q_k), is observed for larger values of N, which, similar to ϕ and Ψ₆, behaves logarithmically for N > 100 as in Fig. 3f.

The long-term logarithmic structural evolution is reminiscent of the aging of a thermal glass below the glass transition temperature, which evolves toward lower energy states characterized by densification, more homogeneous structure, and increased ordering.

Such observations lead to our description of the small amplitude oscillatory deformation as mechanical aging. In terms of dynamics, the decrease in non-affine motion (Fig. 3a) is similar to previous molecular dynamics simulation results,⁹ which suggest that particles should evolve toward configurations with lower potential energies. The slow and long-term increase in the cage-relative MSD (Fig. 3b) is also reminiscent of the MSD observed in aging thermal glasses.^56,57 We note that similar correspondences between mechanical driving and thermal activation have been recognized in other disordered solids, particularly metallic glasses. Atomistic studies have demonstrated that external stress can markedly accelerate relaxation dynamics – effectively mimicking an increase in temperature – through facilitated exploration of a complex, fractal-like energy landscape.^58,59 Moreover, stress-temperature scaling relations have been proposed in steady-state flow,⁶⁰ indicating the existence of quantitative mappings between mechanical and thermal variables. These precedents suggest that the resemblance uncovered here is not incidental but may reflect a more general principle of glassy dynamics. In this light, the present experimental system could, in the future, provide a controlled platform to probe and possibly quantify the correspondence between mechanical forcing and effective thermal aging in granular glasses. A systematic exploration of such mappings, however, is beyond the scope of this manuscript and warrants further study.

5 Void-driven microscopic dynamics

The trends displayed in the previous section point toward a microscopic mechanism that may increase the stability of the system, leading to lower values of 〈D²_min〉 and subdiffusive dynamics, as indicated by Fig. 3a and b. The same conclusion can be drawn from the 〈Ψ₆〉 parameter from Fig. 3d, as particles that arrange in a hexagonal packing can be considered particularly stable in two dimensions. We can further speculate about the governing mechanism through indications of increasing packing fraction and decreasing spatial heterogeneity (Fig. 3c–e), which points to the possibility that loosely packed regions are especially prone to rearrangements that result in local densification and spatial homogenization.

To quantify this possible interplay between structure and deformation, we identify loosely packed regions as voids from experimental images by distinguishing the particle phase from the background oil phase through binarization, with an example shown in Fig. 4a. To evaluate the area of these voids, we again rely on the Delaunay triangulation, which partitions each void into triangles, where the vertices of the Delaunay triangles fall on the boundary particles surrounding a void. For each triangle k, we calculate three areas, including the triangle area A_k, the sum of the particle circular sector areas, A^k_sc, and the enclosed area, A^k_bp, as the area between three constituent particles if they are closely packed. All area types are depicted in Fig. 4b. The void area A_v is then

where N_vt denotes the number of constituent triangles for each void. Under this definition, a void surrounded by three closely packed particles (as in Fig. 4b) will have zero area, as there is no room for its area to further reduce.

In Fig. 4c, we show the distributions of normalized void areas, Ã_v = 4A_v/πd², which are categorized by the number of their constituent particles N_bp. For voids with larger N_bp, both the average and the variance of the void areas are larger.

Note that although voids with N_bp = 3 should have Ã_v = 0, some non-zero values exist due to measurement noise, and they should be neglected as they are only small fractions of a particle's projected area.

We then track the system-averaged void areas as N increases, see Fig. 4d. Initially, 〈Ã_v〉 increases slightly, before decreasing nearly logarithmically after N = 100 to a value significantly lower than the initial average void size. This behavior is consistent with the packing fraction curve in Fig. 3c, which has an initial decrease followed by a logarithmic increase after N = 100.

To understand how voids influence the dynamics of the raft during structural aging, we examine the particle displacement field of the raft, with an example shown in Fig. 5a. Specifically, we compute the average relative displacement field around each void by spatially binning the particle displacement as a function of the particle position relative to the centroid of each void and subtracting the average background displacement. We then discount all bins with less than 20 data points.

Fig. 5 Dynamics of particles near voids. (a) An experimental image overlaid with particle trajectories (blue) and large identified voids (red). Trajectories are from N = 0 to N = 1000. (b) Relative displacement field around a void from N = 0 to N = 1000 with N_bp = 4, averaged across all experiments. (c) Average relative displacement field around a void from N = 0 to N = 1000 with N_bp > 4 across all experiments. (d) Particles colored by D²_min,L with large voids highlighted in red. (e) The averaged D²_min,L of particles as a function of their distance to the closest void, plotted for different void sizes.

The resulting displacement field, when calculated for N_bp = 4, only exhibits small radially outward displacements roughly one particle diameter away from the center, as shown in Fig. 5b. When N_bp > 4, however, a noticeable tendency for particles to move toward the void's center is observed, as shown in Fig. 5c. A slightly higher displacement magnitude is observed in the direction of the free boundaries, highlighting the effect of their presence.

To further quantify the dynamics of the raft near voids, we examine the non-affine particle motion close to the identified voids, with an example visualized in Fig. 5d. As a void typically evolves over many loading cycles, we evaluate eqn (1) using particle positions from the initial and the final cycle of an experiment, yielding D²_min,L, with the subscript L introduced to distinguish from the frame-by-frame D²_min. For each void, we establish a local reference at the void's centroid and average the D²_min,L of the nearby particles at different distances, d_v. The resulting normalized average, D̅²_min,L/d², shown in Fig. 5e, indicates that the non-affinity generally rises close to voids and is higher for larger voids. All void sizes seem to have the same radius of influence of d_v/d ≈ 3, at which D̄²_min,L converges to the background value. Voids with Ã_v < 0.29, however, tend to have slightly lower D̄²_min,L when measured close to the void than the background value measured far away. Based on Fig. 4c, the area threshold of Ã_v ≤ 0.29 roughly corresponds to N_bp ≤ 4, and higher bins are mostly populated by voids with N_bp > 4. It follows that voids with areas above the Ã_v = 0.29 threshold contribute significantly more to the change in structure of the raft, which may be due to lower rearrangement energy barriers in their close proximity. If so, we may expect that a particle located in a region with a high density of voids may be particularly susceptible to rearrangements that result in densification, some instances of which can be seen in Fig. 5a and d.

The results presented in this section help explain the densification and lowered 〈D²_min〉 observed in Section 3. A void with N_bp > 4 tends to see nearby particles move toward its center, decreasing its area and contributing to the overall densification of the raft. As smaller voids induce less non-affine motion, the shrinkage of the voids over time also results in smaller D²_min in the system on average.

6 Evolution of void morphology

Although the dynamics of particles near voids have been explored, we have yet to investigate the way that a void's morphology changes as a result of those dynamics. Therefore, we further quantify the evolution of voids by tracking each void across consecutive loading cycles, driven by the non-affine displacements of nearby particles. For a time series of binarized 2D images (Fig. 4a), we stack them into a 3D spatial-temporal image, in which each void occupies a continuous volume that can be easily identified, eventually yielding a time history of its area Ã_v. We also track events involving several voids merging into one, and a single void splitting into several voids. In this case, we calculate the summed area of all contributing voids for an event, also denoted by Ã_v, so that we can quantify the change in void area, ΔÃ_v, due to a given event.

We classify the morphological changes of the voids into four types of events. We use Common-type to denote a change in a void's area as its constituent particles rearrange, but without change to its connectivity in the 3D spatial-temporal image. For this analysis, we do not include voids with N_bp = 3.

Moving on, the Split-type denotes an event in which a single void splits into multiple voids, while the Merge-type corresponds to multiple voids merging into a single void. Finally, the Other-type represents a complex morphology change, where multiple voids reorganize into multiple new voids, in a way that cannot be decomposed into simple merging or splitting. For all four types, we calculate ΔÃ_v using a small loading cycle interval of ΔN = 2.

Under this definition, we sampled events throughout a cyclic loading test and from 14 tests with different initial packing structures. Fig. 6a shows a frequency count of the change in void area ΔÃ_v for the four types of events, revealing that Common-type events are the most common by several orders of magnitude. Similar numbers of Merge-type and Split-type events are observed, followed by the least frequent Other-type events.

Fig. 6 The change of void morphology during cyclic loading. (a) Frequency count of the area change for different void behavior types. (b) The total void area change (left axis, turquoise) and average individual void area change (right axis, dark blue) due to specific event type. (c) The change in void area vs. current void area for Common-type voids. The shaded regions indicate standard error across experiments. (d) A histogram of the orientations of the principal axes for voids with N_bp > 4 (right half) and N_bp = 4 (left half). Data points taken at N = 0 are shown in dark blue, and data points taken at N = 1000 are shown in turquoise.

For a given type of void event, we then calculate the average area change per void, , and the total area change, ΔÃ^tot_v, see Fig. 6b. On average, Split-type and Merge-type events result in negative and positive area change, respectively. Common-type events have a slightly negative ΔÃ_v on average and the associated magnitude is significantly smaller than that of Merge-type and Split-type events, which are similar to each other. Despite this, Common-type events are so much more frequent that their contribution to the total area change, ΔÃ^tot_v, is greater than any other event type. Additionally, as Merge-type and Split-type events occur at similar frequencies and with opposite signs in , their contributions to ΔÃ^tot_v nearly cancel out. Finally, other-type events have a negative , but their low frequencies prevent any significant ΔÃ^tot_v.

As Common-type events have the greatest contribution to the densification of the raft, we examine the dependence of their void area change, ΔÃ_v, on their current size Ã_v. We binned all applicable voids by Ã_v, and calculated the bin-averaged ΔÃ_v for voids sampled at various N, which is shown in Fig. 6c. The results show that for the smallest area bin, which is almost entirely composed of voids with N_bp = 4, a small, positive, area change is observed. This indicates that these voids grow in small increments more than they shrink, consistent with the small radially outward displacement fields observed in Fig. 5b. It is worth pointing out that the area change due to a void switching from N_bp = 4 to N_bp = 3 is not considered Common-type and is therefore not represented in Fig. 6. For larger voids, however, the change in area is negative and nearly identical across different values of Ã_v. The radially inward displacement field from Fig. 5c corroborates this observation. This helps explain the initial observation of decreasing structural heterogeneity, as initially loose regions of the raft will slowly collapse toward the center of these regions.

The sign switch in between the two smallest area bins provides further evidence of differing void behavior between N_bp > 4 (Ã_v > 0.29) voids and N_bp ≤ 4 (Ã_v ≤ 0.29) voids. This aligns with the tendency of N_bp > 4 voids to contribute significantly more to the raft's dynamics than N_bp ≤ 4 voids, which was first observed in Fig. 5e. There is also a dependence of on N, showing that |ΔÃ_v| is smaller for larger N, similar to the dynamic slow-down observed in Fig. 3a. This suggests that most of the structural change occurs in the earlier stages of cyclic loading.

Finally, void orientation can be investigated as a result of aging. To determine the orientation of a single void, the pixels composing a void from an experimental image are used to identify the major axis of an ellipse with the same normalized second central moments. The orientation, θ, is then defined as the angle that this axis makes with another axis running parallel to the fixed boundaries of the raft. We separate the orientation of voids into two categories, N_bp = 4 and N_bp > 4, to distinguish whether the behaviors are different, motivated by the differences between N_bp = 4 and N_bp > 4 voids observed in Fig. 5 and Fig. 6c. These are further divided into groups of N = 0 or N = 1000 to observe the effects of mechanical aging. The resulting histograms of void orientations are shown in Fig. 6d.

We see that all void sizes initially show a preference toward θ = 0, especially for the case of voids with N_bp = 4. At N = 1000, voids with N_bp > 4 tend to be oriented perpendicularly more often than at N = 0, which is likely a result of the void collapsing with the free lateral raft boundaries shrinking inward. Voids with N_bp = 4 at N = 1000, however, are found most frequently at θ = π/4 and θ = 3π/4, coinciding with the direction of maximum shear stress in a uniaxially deformed material, and also similar to the orientation of the shear bands observed in the next section.

The re-orientation of these small, potentially stable voids toward two distinct directions (π/4 and 3π/4) is reminiscent of the bistable elements, or hysterons, thought to underpin the aging behavior of cyclically driven disordered systems.⁶¹ Recent theoretical work has shown that frustrated interactions between such hysteretic elements can lead to the persistent logarithmic aging observed in cyclically driven systems, which continuously explores new configurations with reduced dissipation and less participating elements.²⁴ While voids with N_bp > 4 shrink and voids N_bp = 3 remain rigid, whether voids with N_bp = 4 serve as physical realization for such elements remains an intriguing direction for future study.

7 Implications of aging on raft failure

To investigate how the structural changes induced by cyclic loading affect the yielding behavior of the raft, we subject aged rafts (loaded to N = 1000), as well as non-aged rafts (without any cyclic loading history), to uniaxial tension and uniaxial compression. We measure F_i/γW₀ against the global strain ε = δ_i/L₀, where i = t for tension and i = c for compression. We also visualize the deformation behavior of the rafts by overlaying experimental images with particles colored by their respective D²_min values, which are calculated over a strain interval Δε = 0.001. The results are displayed in Fig. 7.

Fig. 7 Uniaxial tension and compression tests on the non-aged and aged rafts. (a) Stress–strain curves for the non-aged (light orange) and aged (dark orange) rafts in tensile tests, averaged over 16 and 14 repetitions for each experiment, respectively. (b) Stress–strain curves for the non-aged (light blue) and aged (dark blue) rafts in compression tests, averaged over five and eight repetitions for each experiment, respectively. (c) A visualization of the particle D²_min in the non-aged and aged rafts undergoing compression and tension, shown for various global strains. (d) The average void area vs. global strain for the cases in (a) and (b) using the same color scheme. The shaded regions in (a), (b), and (d) indicate test-to-test fluctuations.

In tensile tests, both aged and non-aged rafts display a similar tensile stiffness, measured as the slope of the initial increment of the loading curve, until the onset of yielding at ε ≈ 0.01. However, the strength of the aged rafts, defined as the maximum value of F_t/γW₀, is slightly larger than that of the non-aged rafts, which is shown in Fig. 7a. Also associated with the aged rafts is a faster decay in the tensile stress following ε > 0.01. These results indicate that aged rafts are slightly stronger and more brittle during tensile deformation. The visualization of two example rafts (Fig. 7c), representing aged and non-aged rafts, reveals shear banding that is narrower in the aged case through more concentrated regions of high D²_min.

While the loading curves were similar up to yielding in the tensile tests of the aged and non-aged rafts, there are clear differences between the two types in compression tests as shown in Fig. 7b. The relation between F_c/γW₀ and ε indicates a significantly higher stiffness for the aged rafts, as well as a significantly higher yield strength at ε ≈ 0.03, which coincides with a unified yield strain value for a variety of disordered solids.⁶² Beyond the yielding point, the two loading curves are rather distinct. While the stress for the non-aged rafts continuously increases until it plateaus, a significant overshoot exists in the curve for the aged raft, followed by a stress drop toward the same plateau as in the non-aged case, indicating a possibly history-independent behavior. The stress drop in the aged rafts before the plateau is indicative of a brittle behavior, while the continuous stress increase in the non-aged rafts indicates ductile behavior. Correspondingly, the shear bands tend to be narrower in the aged rafts, with the particles inside having higher D²_min, as seen in Fig. 7c.

Another difference between tension and compression tests manifests in the orientation of the shear band during failure. The orientation of the shear band during compression is above π/4, and below π/4 during tensile tests, suggesting that volumetric effects and friction may be important for the yielding process.^63,64

As a final analysis, we investigate the behavior of voids during the uniaxial loading tests by examining the average void area 〈Ã_v〉 vs. ε, which can be calculated using eqn (5) and normalized in the same way as in Section 5. In all tests, 〈Ã_v〉 increases with ε, showing that the raft generally dilates regardless of the deformation mode. The plot of 〈Ã_v〉 against ε is shown in Fig. 7d. For tension, this is a straightforward outcome due to the imposed extensional volumetric strain. The dilation observed in the compression tests, however, indicates that the deviatoric part of the imposed uniaxial deformation dominates the volumetric part, as deviatoric deformation often induces dilation in dense granular materials. This suggests that the densification of the raft played a role in its lowered ductility, as greater dilation was required to reach steady-state behavior for the aged rafts, which is a signature of brittle-like failure.

8 Conclusions

We have demonstrated that a cohesive granular raft subject to small amplitude cyclic loading experiences structural aging with a long-term logarithmic increase in the raft's density. This density change is possible due to the presence of free boundaries and uniaxial loading conditions, which promotes changes to volume. Other parameters tied closely to the structure of the raft, such as the bond orientation parameter and structural heterogeneity, also change logarithmically with cycle number, consistent with the dynamic slow-down in the cage-relative MSD and 〈D²_min〉.

To address the source of structural evolution during aging, we monitor voids, defined as discrete areas associated with gaps in the particle packing structure. We see that the area of voids also decreases logarithmically in the long term, but that large and small voids exhibit very different behaviors. Large voids tend to displace nearby particles toward their center during aging, shrinking the void's area and increasing the density of the raft. Additionally, using D̄²_min,L as a metric for structural evolution, we see that larger void areas tend to induce more structural change in their close vicinity. Further analysis suggests that although large scale rearrangements individually contribute the most to the area change of a void, it is the smaller rearrangements that collectively dominate the density change of the raft. Thus, the structural change of the raft is associated closely with the shrinking of large voids, and the density change is most influenced by small but commonplace rearrangements to those voids.

Smaller voids, which can be classified as having four constituent boundary particles, behave very differently. On average, these voids tend to grow slightly, but on a magnitude appreciably smaller than the area change associated with larger voids. Despite this small contribution to total raft density, we see that these voids are not inactive, as a distribution of their orientations changes to align primarily with an angle of π/4 radians with the horizontal after 1000 loading cycles, despite initially aligning with the horizontal axis. Larger voids, meanwhile, see very little change to the distribution of their orientations.

The structural changes to the raft culminate in altered deformation properties and failure mechanisms. Aged rafts undergoing tensile tests have a higher strength than non-aged rafts, but also exhibit more brittle behavior. A similar trend is observed under compression, with the added observation of increased stiffness and a stress overshoot, indicating brittle failure with dilation. This additional dilation is due to the increased density of the raft, and thus the shrinking of large voids. Regardless of preparation history, however, all rafts under compression tend to a history-independent behavior at large strains. Investigating the area of voids during these uniaxial tests confirms particularly strong dilation during the compression of aged rafts. Also observed is a shear band close to π/4 radians for all failure tests, but with a higher or lower angle depending on whether the test was under compression or tension, respectively.

It should be noted that our results may depend on particle diameter. Previously, we showed that rafts having particles d = 3.3 mm (larger than the capillary length l_c) deform with a more brittle behavior, while rafts with d = 0.7 mm (smaller than l_c) deform similarly to our current rafts with d = 1.0 mm particles.⁴³ The role of particle diameter and therefore the corresponding capillary interactions in cyclic loading should be examined in future studies.

Our experiment has potential applications in the study of yielding disordered solids, as we can intentionally trigger a stress overshoot during yielding via well-controlled cyclic loading. Our methods for monitoring the raft's microscopic structure and dynamics could be used to better understand how the preparation history of an out-of-equilibrium disordered solid influences the nature of its yielding transition.^3,4,65 That said, the structural aging of rafts under volume-conserving deformation should be studied and compared with our results. The capillary-induced particle attraction and the particle friction also render our experimental system suitable to study whether complex particle interactions give rise to rugged energy landscapes of an amorphous solid, exhibiting possibly perpetual irreversibility.^38,40,54

Conflicts of interest

There are no conflicts to declare.

Data availability

Data for this article are available on the Deep Blue Data Repositories of the University of Michigan. https://doi.org/10.7302/tf47-ch94.

Acknowledgements

The authors would like to acknowledge the funding from the National Science Foundation grant MRSEC/DMR-1720530 and grant CMMI-2519512. We also thank Andrea J. Liu for helpful discussions.

References

1. D. Bonn, M. M. Denn, L. Berthier, T. Divoux and S. Manneville, Rev. Mod. Phys., 2017, 89, 035005.
2. A. Nicolas, E. E. Ferrero, K. Martens and J.-L. Barrat, Rev. Mod. Phys., 2018, 90, 045006.
3. L. Berthier, G. Biroli, L. Manning and F. Zamponi, Nat. Rev. Phys., 2025, 7, 1–18.
4. M. Ozawa, L. Berthier, G. Biroli, A. Rosso and G. Tarjus, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 6656–6661.
5. M. Popovi, T. W. de Geus and M. Wyart, Phys. Rev. E, 2018, 98, 040901.
6. H. J. Barlow, J. O. Cochran and S. M. Fielding, Phys. Rev. Lett., 2020, 125, 168003.
7. M. Singh, M. Ozawa and L. Berthier, Phys. Rev. Mater., 2020, 4, 025603.
8. W.-T. Yeh, M. Ozawa, K. Miyazaki, T. Kawasaki and L. Berthier, Phys. Rev. Lett., 2020, 124, 225502.
9. N. V. Priezjev, J. Non-Cryst. Solids, 2022, 590, 121697.
10. A. D. Parmar, S. Kumar and S. Sastry, Phys. Rev. X, 2019, 9, 021018.
11. P. Leishangthem, A. D. Parmar and S. Sastry, Nat. Commun., 2017, 8, 14653.
12. J. Pollard and S. M. Fielding, Phys. Rev. Res., 2022, 4, 043037.
13. R. Sharma and S. Karmakar, Nat. Phys., 2025, 21, 253–261.
14. P. Das, A. D. Parmar and S. Sastry, J. Chem. Phys., 2022, 157, 044501.
15. N. V. Priezjev, J. Non-Cryst. Solids, 2018, 479, 42–48.
16. M. Utz, P. G. Debenedetti and F. H. Stillinger, Phys. Rev. Lett., 2000, 84, 1471.
17. M. Fan, M. Wang, K. Zhang, Y. Liu, J. Schroers, M. D. Shattuck and C. S. O'Hern, Phys. Rev. E, 2017, 95, 022611.
18. J. Shen, Y. Huang and J. Sun, J. Mater. Res., 2007, 22, 3067–3074.
19. E. Bouchbinder and I. Procaccia, Phys. Rev. E, 2013, 87, 042310.
20. J. Ketkaew, W. Chen, H. Wang, A. Datye, M. Fan, G. Pereira, U. D. Schwarz, Z. Liu, R. Yamada and W. Dmowski, et al., Nat. Commun., 2018, 9, 3271.
21. G. Zhang, H. Xiao, E. Yang, R. J. Ivancic, S. A. Ridout, R. A. Riggleman, D. J. Durian and A. J. Liu, Phys. Rev. Res., 2022, 4, 043026.
22. H. Xiao, G. Zhang, E. Yang, R. Ivancic, S. Ridout, R. Riggleman, D. J. Durian and A. J. Liu, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, e2307552120.
23. N. V. Priezjev, J. Non-Cryst. Solids, 2021, 566, 120874.
24. D. Shohat, P. Baconnier, I. Procaccia, M. van Hecke and Y. Lahini, Proc. Natl. Acad. Sci. U. S. A., 2026, 123, e2515075123.
25. N. C. Keim and P. E. Arratia, Phys. Rev. Lett., 2014, 112, 028302.
26. P. K. Jana and N. V. Priezjev, J. Non-Cryst. Solids, 2023, 600, 121996.
27. H. Li, H. Xiao, T. Egami and Y. Fan, Mater. Today Phys., 2024, 49, 101582.
28. Y. Jin, P. Urbani, F. Zamponi and H. Yoshino, Sci. Adv., 2018, 4, eaat6387.
29. N. C. Keim and D. Medina, Sci. Adv., 2022, 8, eabo1614.
30. S. Slotterback, M. Mailman, K. Ronaszegi, M. Van Hecke, M. Girvan and W. Losert, Phys. Rev. E, 2012, 85, 021309.
31. J. R. Royer and P. M. Chaikin, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 49–53.
32. M. Zhang, L. H. Dai and L. Liu, J. Appl. Phys., 2014, 116, 053522.
33. Y. Wang, M. Zhang and L. Liu, Scr. Mater., 2015, 102, 67–70.
34. N. W. Mueggenburg, Phys. Rev. E, 2005, 71, 031301.
35. P. Richard, M. Nicodemi, R. Delannay, P. Ribiere and D. Bideau, Nat. Mater., 2005, 4, 121–128.
36. J. Zhang, T. Majmudar, A. Tordesillas and R. Behringer, Granul. Matter, 2010, 12, 159–172.
37. S. Farhadi, A. Z. Zhu and R. P. Behringer, Phys. Rev. Lett., 2015, 115, 188001.
38. H. Xiao, A. J. Liu and D. J. Durian, Phys. Rev. Lett., 2022, 128, 248001.
39. Y. Zhao, Y. Zhao, D. Wang, H. Zheng, B. Chakraborty and J. E. Socolar, Phys. Rev. X, 2022, 12, 031021.
40. Y. Yuan, Z. Zeng, Y. Xing, H. Yuan, S. Zhang, W. Kob and Y. Wang, Nat. Commun., 2024, 15, 3866.
41. J. B. Knight, C. G. Fandrich, C. N. Lau, H. M. Jaeger and S. R. Nagel, Phys. Rev. E, 1995, 51, 3957.
42. Y. Yuan, W. Kob and H. Tanaka, Proc. Nat. Acad. Sci. U. S. A., 2026, 123, e2528600123.
43. H. Xiao, R. Ivancic and D. Durian, Soft Matter, 2020, 16, 8226–8236.
44. K.-Î. Tô and S. R. Nagel, Soft Matter, 2023, 19, 905–912.
45. B. C. Druecke, R. Mukherjee, X. Cheng and S. Lee, Phys. Rev. Fluids, 2023, 8, 024003.
46. S. Protière, Annu. Rev. Fluid Mech., 2023, 55, 459–480.
47. X. Yan, T. C. Watson and H. Xiao, arXiv, preprint, arXiv:2602.17539, 2026.
48. M. Nicolson, Math. Proc. Camb. Philos. Soc., 1949, 288–295.
49. M.-J. Dalbe, D. Cosic, M. Berhanu and A. Kudrolli, Phys. Rev. E, 2011, 83, 051403.
50. J. M. Rieser, Deformation of two-dimensional amorphous granular packings, University of Pennsylvania, 2015.
51. M. Harrington and D. J. Durian, Phys. Rev. E, 2018, 97, 012904.
52. M. Harrington, H. Xiao and D. J. Durian, Granul. Matter, 2020, 22, 1–9.
53. M. L. Falk and J. S. Langer, Phys. Rev. E, 1998, 57, 7192.
54. W. Li, J. M. Rieser, A. J. Liu, D. J. Durian and J. Li, Phys. Rev. E, 2015, 91, 062212.
55. J. M. Rieser, C. P. Goodrich, A. J. Liu and D. J. Durian, Phys. Rev. Lett., 2016, 116, 088001.
56. B. Seoane and F. Zamponi, Soft Matter, 2018, 14, 5222–5234.
57. A. P. Hammond and E. I. Corwin, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 5714–5718.
58. P. Cao, M. P. Short and S. Yip, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 13631–13636.
59. C. Liu and Y. Fan, Phys. Rev. Lett., 2021, 127, 215502.
60. P. Guan, M. Chen and T. Egami, Phys. Rev. Lett., 2010, 104, 205701.
61. N. C. Keim, J. D. Paulsen, Z. Zeravcic, S. Sastry and S. R. Nagel, Rev. Mod. Phys., 2019, 91, 035002.
62. E. D. Cubuk, R. Ivancic, S. S. Schoenholz, D. Strickland, A. Basu, Z. Davidson, J. Fontaine, J. L. Hor, Y.-R. Huang and Y. Jiang, et al., Science, 2017, 358, 1033–1037.
63. J. Bardet, Comput. Geotechnol., 1990, 10, 163–188.
64. K. Karimi and J.-L. Barrat, Sci. Rep., 2018, 8, 1–10.
65. R. L. Moorcroft, M. E. Cates and S. M. Fielding, Phys. Rev. Lett., 2011, 106, 055502.

---