Critical fluctuations of elastic moduli in jammed solids

Abstract

We investigate sample-to-sample fluctuations of the shear modulus in ensembles of particle packings near the jamming transition. Unlike the average modulus, which exhibits distinct scaling behaviours depending on the interparticle potential, the fluctuations obey a critical exponent that is independent of the potential. Furthermore, this scaling behaviour has been confirmed in two-dimensional packings, indicating that it holds regardless of spatial dimension. Using this scaling law, we discuss the relationship predicted by heterogeneous-elasticity theory between elastic-modulus fluctuations and the Rayleigh scattering of sound waves across different pressures. Our numerical results provide a useful foundation for developing a unified theoretical description of the jamming critical phenomenon.

---

1. Introduction

When particles are compressed, the system jams and acquires rigidity at a certain packing fraction. The phenomenon, called the jamming transition, is ubiquitous in diverse materials such as foams, colloidal suspensions, pastes, and granular materials, and forms a basis for understanding properties of disordered systems, spanning from structural glasses to biological systems.^1,2 Over the past 30 years, various studies, including theoretical,^3,4 numerical^5,6 and experimental^7,8 approaches, have characterised basic properties of this phenomenon. The most central observation is made in critical phenomena; various physical quantities exhibit power-law scaling near the transition.^5,6 This critical behaviour manifests in elastic moduli,^6,9 vibrations,^10,11 transport properties,^12–14 and rheology.¹⁵ Furthermore, at the transition point, the number of contacts exactly matches the number of degrees of freedom in the system, a state known as isostaticity.⁶ It is now well established that the jamming scaling law is controlled by the contact number excess to the isostatic value.^3,4

One of the central challenges in statistical physics is the understanding of phase transitions and critical phenomena.^16–19 In this context, theoretical tools such as mean-field theory and, more importantly, the renormalization group have been developed, and it has become clear that fluctuations play a crucial role in critical behaviour. In particular, as a system approaches a critical point, both the fluctuations and the associated correlation length scales diverge, leading to various anomalies at the critical point. In the field of jamming transition, some works have followed this line of research, including scaling collapse and finite-size scaling,^9,20,21 diverging length scales,^3,10,22,23 and Widom scaling hypothesis.²⁴ Thus, understanding how fluctuations in physical quantities behave near the jamming critical point remains an important issue. Despite its relevance, the study of fluctuations in jamming systems is still limited, and further development is anticipated.^25–28

Beyond their relevance in statistical physics, fluctuations in elastic moduli play an essential role in the theoretical description of the anomalous properties of amorphous solids. Among the various theoretical frameworks,^29–34 the heterogeneous-elasticity theory (HET)^35–39 is one of the most well-established approaches for describing the anomalous behaviour of amorphous materials such as jammed solids. HET postulates that the local shear modulus exhibits spatial fluctuations throughout the system, and extensive numerical studies have indeed confirmed the existence of such elastic heterogeneity in amorphous systems.^40–42 Within this theoretical framework, the so-called disorder parameter γ, which quantifies the variance of local shear moduli, acts as a key control parameter.³⁶ In the standard formulation of HET, the distribution of local shear moduli is assumed to be Gaussian, and its mean value and standard deviation serve as inputs for computing macroscopic observables. This formulation successfully accounts for a variety of universal features of amorphous solids. For example, the low-frequency vibrational modes in excess over the Debye prediction (the boson peak) can be captured within HET, with both their characteristic frequency and intensity expressed as functions of γ. Likewise, the theory predicts a Rayleigh-type frequency dependence of the acoustic attenuation, Γ∼Ω^d+1 (Γ is the attenuation coefficient and Ω is the propagation frequency), which is also observed in experiments on glasses.^43,44 The prefactor of this relation, corresponding to the strength of Rayleigh scattering, is again determined by γ, linking vibrational and transport anomalies to the underlying elastic heterogeneity. These theoretical predictions, spanning both vibrational spectra^14,45,46 and acoustic attenuation,^14,47,48 have been extensively tested through numerical simulations.

Another major theoretical framework developed to describe and predict the properties of amorphous solids is the effective medium theory (EMT), which is based on spring-network models.⁴⁹ While HET treats the spatial fluctuations of local shear moduli as the sole mesoscale input, EMT instead relies on microscopic parameters, namely the mean coordination number z and the contact strain e, often referred to as the pre-stress. Despite this difference in the scale of description, EMT successfully reproduces the same key features of vibrational and scattering phenomena as HET, including the emergence of the boson peak and the Rayleigh-type acoustic attenuation. Given that both theories capture similar physical phenomena from different levels of description, it is natural to expect that they are connected via an appropriate coarse-graining procedure. However, the precise nature of this connection remains elusive. One major reason for this persistent gap is the limited numerical characterisation of fluctuations of elastic modulus. From this perspective, a detailed characterisation of local shear modulus fluctuations is of particular theoretical significance, particularly for bridging the gap between these two frameworks.

In this work, we characterise the sample-to-sample fluctuations of elastic moduli in ensembles of jammed packings by means of numerical simulation. These fluctuations are believed to encode information equivalent to the disorder parameter in HET.^50,51 We focus on two prototypical models of jammed amorphous solids: harmonic spheres and Hertzian spheres. Our results reveal that, while the average elastic modulus exhibits a model-dependent scaling behaviour, its fluctuations diverge in a manner independent of the interaction potential. This critical scaling is also observed in two-dimensional packings. Finally, we discuss the theoretical implications of these findings, with particular emphasis on their relevance to the HET prediction on acoustic attenuation in amorphous solids.

2. Methods

2.1. System description

We consider randomly jammed packings in a three-dimensional box with periodic boundary conditions. Particles interact with a finite-range purely repulsive potential⁶where r_ij is the distance between particles i and j, σ_ij = (σ_i + σ_j)/2 is the sum of radii of two particles, and H(x) is the Heaviside step function. Two models of jammed particles are considered, harmonic spheres (α = 2) and Hertzian spheres (α = 5/2). For both models, we use 50 : 50 binary mixtures of particles whose diameters are σ and 1.4σ with identical mass m. The number of particles in the system is denoted as N. Masses, lengths, and energies are measured in units of m, σ, and ε, respectively. The number density is expressed as ρ = N/L^d, where d = 3 is the spatial dimension.

2.2. Numerical generation of jammed packings

We start the preparation of packings by minimising the energy of random configurations at a sufficiently high packing fraction φ_init = 1.0 with the FIRE algorithm.⁵² The algorithm is terminated when the condition is reached. Then, we perform the global compression/decompression to generate packings at a given target pressure p.⁵³ For harmonic spheres, pressures range over p = 10⁻², 10⁻³, …, 10⁻⁶ for all system sizes. For Hertzian spheres, the highest pressure is p = 10⁻² for all system sizes, and the lowest pressure is p = 10⁻⁸ for N = 2000. For the case of N = 32 000, the lowest pressure is p = 10⁻⁷. A packing ensemble prepared with the compression/decompression-only protocols may contain samples that have finite pressure but are unstable to shear.⁵⁴ We discard such samples with a negative shear modulus when constructing the ensembles. The sizes of ensembles used in this study are given in Appendix A. After packings are generated, we recursively remove rattler particles whose number of contacts N_c satisfies N_c ≤ d.⁵⁵

2.3. Calculation of stressed and unstressed shear moduli

Here, we recap the linear response formulation of shear moduli used in this study.^56–58 Because the system is isotropic, we measure the elastic response exclusively for a shear strain in the x direction with a strain gradient in the y direction.

The elastic modulus is decomposed into two components, the affine and non-affine modulus: G = G_A − G_N.^59,60 The affine modulus (Born–Huang term) is the second derivative of total potential energy E, defined as

with respect to the affine strain, evaluated as^{56–58,61–63}where V is the volume of the simulation box, f_ij = −v′(r_ij) and k_ij = v″(r_ij) are the first and second derivatives of the potential v(r_ij) with respect to r_ij. The distance r_ij is given by the norm of the displacement vector r_ij = (x_ij, y_ij, z_ij). The sum is taken for all bonds 〈ij〉 in the system. We note that the second term in the summation of eqn (3) is a correction term for the finite pressure.^64,65

The second term G_N is due to non-affine relaxation. This term is evaluated with the Hessian matrix^56–58

where the i-th block of the affine force Ξ (3N-dimensional vector) is provided asand the Hessian matrix (size 3N × 3N) is defined as
where r_i = (x_i, y_i, z_i) is the coordinate of particle i. In amorphous solids, affine deformation causes an additional non-affine displacement field δu = [script letter H]⁻¹·Ξ.^56–59,66 This displacement field δu (3N-dimensional vector) is usually obtained though the decomposition of Ξ with the normal modes of [script letter H].^56–58 However, this method is computationally hard when the size of matrix becomes large. Therefore, we obtain the non-affine displacements δu via the FIRE minimisation of a cost functioninstead of performing the direct diagonalization of the Hessian matrix.⁶⁷ In this minimisation, the termination condition is set to . This iterative method yields displacements that are identical to the direct solution of the linear equation [script letter H]·δu = Ξ within numerical tolerance.

In the harmonic limit, the energy change ΔE caused by relative displacements {u_i} from the mechanically equilibrated position {r_i} takes the form^3,4,9,61,68

where u_ij,‖ and u_ij,⊥ denote the components of the displacement parallel and perpendicular to r_ij, respectively. This may be written compactly using the Hessian matrix: ΔE = u·[script letter H]·u. Because the contact force f_ij = −v′(r_ij) is always repulsive in the jammed packings, the original state is called the stressed system. On the other hand, we have also studied the shear modulus of the unstressed system,⁴ where we remove f_ij in the Hessian matrix while maintaining k_ij. This case corresponds to the system where we replace the interactions between particles with relaxed (unstretched) springs of each stiffness k_ij. In practice, we omit the first derivatives in the Hessian matrix in the cost function (7) when calculating the shear modulus of the unstressed system. We denote the shear modulus of the stressed system by G₁ and that of the unstressed system by G₀.

3. Behaviour of average shear modulus

First, we present the critical behaviour of the ensemble-averaged shear modulus. As reported in numerous studies,^6,9,69,70 the shear modulus exhibits the following critical scaling with pressure p:

G ∼ p^{(α−3/2)/(α−1)}. (9)

In addition, the excess contact number δz = z − z_iso, which quantifies the distance from the isostatic point at the onset of jamming,^4,6 is known to scale as⁶:

δz ∼ p^1/(2α−2). (10)

Here, z_iso = 2d(1 − 1/N) is the isostatic contact number that satisfies the Maxwell criterion, and the 1/N correction arises from the global zero-frequency translational modes of the finite-size systems.²⁰ Combining the two scaling relations yields G ∼ δz for harmonic spheres (α = 2), and G ∼ δz² for Hertzian spheres (α = 5/2). As clearly shown in Fig. 1, our numerical data are consistent with these predictions: harmonic spheres exhibit the scaling 〈G₁〉 ∼ δz, while Hertzian spheres obey 〈G₁〉 ∼ δz², where 〈·〉 denotes the ensemble average.

Fig. 1 Critical scaling of the average shear moduli 〈G₁〉 and 〈G₀〉/2 for harmonic and Hertzian packings with N= 32000. The solid and dashed lines represent 〈G〉 ∼ δz and 〈G〉 ∼ δz², respectively.

We now turn to the unstressed modulus G₀. Similar to the stressed modulus G₁, the unstressed modulus G₀ follows the same scaling laws for harmonic and Hertzian spheres, as shown in Fig. 1. In addition, according to the EMT for spring networks, the stressed shear modulus G₁ is expected to be twice the unstressed modulus G₀ near the jamming point.⁴⁹ To test this prediction, we plot G₀/2 for both types of packings in Fig. 1. Near the jamming transition, 〈G₀〉/2 collapses onto 〈G₁〉, indicating good quantitative agreement with EMT. However, further from the transition, 〈G₁〉 exceeds 〈G₀〉/2, suggesting that the pre-stress e at higher pressure deviates from its critical value e_c.⁴⁹ Overall, Fig. 1 clearly demonstrates that both 〈G₀〉 and 〈G₁〉 obey the same scaling law within each model, while the scaling exponent itself depends on the interaction potential.

4. Divergence of critical fluctuations

In this section, we examine the sample-to-sample probability distributions within the ensemble used to compute the average quantities discussed above. The discussion is divided into two subsections, focusing respectively on the contact number and the shear modulus. In each case, we analyse the probability distribution of the observable of interest and quantify its fluctuations.

4.1. Excess contact number

We begin by presenting the sample-to-sample probability distribution of the excess contact number δz. Fig. 2(a) shows the distribution P([small delta, Greek, circumflex]z) computed from the ensemble of Hertzian packings with N = 32 000. To facilitate comparison across different pressures, we rescale the variable as [small delta, Greek, circumflex]z = δz/〈δz〉 − 1. As illustrated in Fig. 2(a), the distribution is symmetric and approximately Gaussian in shape. As the pressure decreases and the system approaches the jamming transition, the distribution broadens.

Fig. 2 Probability distributions of (a) the excess contact number δz and (b) the shear modulus G for the ensemble of Hertzian packings (N = 32 000) at various pressures. To enable comparison between different pressures, each quantity X is rescaled as X̂ = X/〈X〉 − 1. In panel (b), the distributions of the stressed modulus G₁ are shown as solid curves, whereas those of the unstressed modulus G₀ appear as dashed curves.

We next quantify the pressure dependence of these fluctuations across the ensemble. To do so, we adopt a disorder quantifier for a generic observable X, defined as²⁸

where σ_X is the standard deviation. We use this measure to characterise the fluctuations in the excess contact number (X = δz).

Fig. 3 displays χ_δz for both interaction types across various system sizes. Harmonic spheres exhibit scaling behaviour χ_δz ∼ δz^−2/5, consistent with the findings of Giannini and coworkers,²⁸ albeit subject to finite-size effects. On the other hand, contact-number fluctuations χ_δz of the Hertzian spheres follow a similar scaling behaviour with δz. Giannini and coworkers identified distinct regimes within the jammed phase based on the behaviour of χ_δz. They reported that χ_δz becomes independent of system size at high pressures, where the scaling becomes steeper, with a more negative exponent. According to their analysis, this crossover occurs at δz ≈ 7 × 10⁻¹, and they referred to the corresponding high-δz region as the glassy regime (δz ≳ 7 × 10⁻¹). Our results for both harmonic and Hertzian spheres are consistent with these observations, although we do not explore this aspect in further detail here.

Fig. 3 Critical fluctuations of the excess contact number δz. The dashed line represents the scaling χ_δz ∼ δz^−2/5. The shaded region marks the glassy regime, δz ≳ 7 × 10⁻¹, where this scaling behaviour is no longer expected to hold.²⁸

4.2. Shear modulus

We now turn to the sample-to-sample fluctuations of the shear modulus. Fig. 2(b) displays the probability distributions of the shear moduli computed from the ensemble of Hertzian packings. Both the stressed modulus G₁ and the unstressed modulus G₀ are shown, represented by solid and dashed curves, respectively.

Compared with the excess contact number δz (Fig. 2(a)), both moduli exhibit broader distributions at the same system size. As the pressure is lowered and the system approaches the jamming transition, these distributions become broader, mirroring the trend observed in P(δz). The distinction between G₁ and G₀ is also evident in Fig. 2(b). At a given pressure, G₀ displays a narrower distribution than G₁. Although both distributions broaden as pressure decreases, G₀ consistently remains narrower than G₁, yet both are wider than the distribution of δz.

The shapes of the distributions also differ notably. At high pressure (p = 10⁻²), both G₀ and G₁ exhibit symmetric distributions. However, at lower pressures (e.g., p = 10⁻⁴), the peak of P(Ĝ₁) shifts towards positive values (Ĝ₁ > 0), and the distribution becomes skewed. This asymmetry becomes more pronounced as the system approaches the jamming point (p = 10⁻⁶). By contrast, P(Ĝ₀) retains a symmetric, Gaussian-like shape throughout this pressure range.

Such non-Gaussian behaviour in the distribution of G₁ has been previously reported,⁵⁰ and the skewness is attributed to spatially localised vibrational modes at low frequencies, which are suppressed in the unstressed system.¹¹ Moreover, outliers are found only in the distribution of the stressed modulus. To quantify this, we examine the minimum and maximum values of Ĝ₀ and Ĝ₁ within the ensemble. While these rare samples are not visible in Fig. 2(b) due to limited statistics, we observe that |minĜ₁|/maxĜ₁ ≈ 2.5 to 4, indicating the presence of samples with anomalously small G₁. Consequently, a Shapiro–Wilk test applied to the raw data returns p-values close to zero, confirming strong deviation from normality. In contrast, for the unstressed modulus G₀, we find |minĜ₀|/maxĜ₀ ≈ 1, indicating an absence of outliers. The Shapiro–Wilk test on G₀ returns p-values of approximately 0.8 or larger, even without removing any outlier data, suggesting consistency with a Gaussian distribution.

We now quantify the sample-to-sample fluctuations in the shear moduli. As noted above, the distribution of G₁ is markedly different from Gaussian and includes significant outliers, making the standard definition of the fluctuation χ_G problematic. To address this, a more robust, median-based measure has been proposed:⁷¹

which avoids the need for explicit outlier-exclusion protocols while remaining statistically meaningful. We adopt this definition for both stressed and unstressed moduli (X = G₀, G₁).

Fig. 4 shows the fluctuation of the shear modulus as a function of δz for both moduli. For the unstressed modulus G₀, the fluctuations exhibit a clear scaling relation: χ_G ∼ δz^−1/2. Finite-size effects appear negligible in this case, consistent with the narrower distribution of G₀ compared to G₁ (Fig. 2(b)). For the stressed modulus G₁, the fluctuations approach the same scaling χ_G ∼ δz^−1/2 as the system size increases. At higher pressures (i.e., larger δz), the scaling breaks down, consistent with the identification of the glassy regime at δz ≳ 7 × 10⁻¹ by Giannini and coworkers.²⁸ An important conclusion from Fig. 4 is that, in contrast with the average shear modulus (Fig. 1), the fluctuation χ_G appears independent of the interaction potential exponent α. For harmonic spheres, the scaling exponent is in agreement with previous findings.^21,28

Fig. 4 Critical fluctuations of the shear moduli. The unstressed modulus G₀ and stressed modulus G₁ are shown in panels (a) and (b), respectively. The dashed lines indicate the scaling χ_G ∼ δz^−1/2. The legend is identical to that in Fig. 3. The shaded region denotes the glassy regime.²⁸

5. Results on two-dimensional packings

Finally, we perform a similar analysis on two-dimensional packings. We use ensembles of binary disks interacting via harmonic and Hertzian potentials described in Section 2.⁶ As in three dimensions, the ensemble sizes used for this calculation are listed in Appendix A. Fig. 5 shows the jamming scaling of fluctuations of the excess contact number (panel (a)), the unstressed shear modulus (panel (b)), and the stressed shear modulus (panel (c)). Analysis of the excess contact number δz reveals that the contact-number fluctuation χ_δz in two dimensions follows a different critical exponent than in three dimensions: χ_δz ∼ δz^−3/5 (Fig. 5(a)), based on our numerical data. The dimensional difference in contact-number fluctuations was also reported in the literature²⁶ using a slightly different metric. By contrast, our analysis of shear moduli demonstrates the shear-moduli fluctuations in two dimensions diverge in the same manner as in three-dimensional systems, χ_G ∼ δz^−1/2 (Fig. 5(b) and (c)). This observation is consistent with previous results for harmonic packings,²¹ which suggested σ_G1/〈G₁〉 ∼ (pN²)^−1/4 ∼ δz^−1/2 in both two and three dimensions. Our results therefore indicate that shear-modulus fluctuations are dimensionally independent, whereas contact-number fluctuations are not. Notably, the robustness with respect to the interaction potential holds for fluctuations of both quantities.

Fig. 5 Critical fluctuations in two-dimensional disk packings of harmonic and Hertzian interactions. Panels (a), (b), and (c) display data of the excess contact number, the unstressed modulus, and the stressed modulus, respectively. In panel (a), the dashed line in black represents χ_δz ∼ δz^−3/5. Dashed lines in panels (b) and (c) represent χ_G ∼ δz^−1/2. Circle markers represent harmonic disks and square markers represent Hertzian disks. System sizes for each interaction potential are given in Appendix A.

6. Conclusion and discussions

In this paper, we have investigated the sample-to-sample fluctuations of the shear modulus and contact number in two models of jammed packings. We have found that, while the average values exhibit critical exponents that depend on the interaction potential, the critical exponents of fluctuations display potential-independent values with respect to δz. Specifically, the fluctuation indicators scale as χ_δz ∼ δz^−2/5 for the excess contact number and χ_G ∼ δz^−1/2 for the shear modulus in three dimensions. The exponent of the shear modulus is identical for both stressed and unstressed systems. While the exponent of excess contact number is dimension-dependent, the exponent of shear moduli is identical for two- and three-dimensional packings.

As stated in the Introduction, the divergence of fluctuations implies the divergence of associated correlation lengths. Our results indicate that the exponent governing the length scale associated with contact-number fluctuations depends on spatial dimension, whereas the length scale associated with elastic-modulus fluctuations exhibits the same scaling behaviour in both two and three dimensions. The exponent controlling the divergence of elastic-modulus fluctuations coincides with that of the correlation length l_c ∼ δz^−1/2, which characterises low-frequency vibrational modes and diverges at the jamming transition.¹⁰ Crucially, the shear modulus of jammed packing is significantly influenced by non-affine responses to the deformation, which are known to be spatially organised over l_c under global shear.⁷² Just as the four-point susceptibility reflects the growth of domain size of dynamic heterogeneity in glass-forming liquids, the divergence of χ_G manifests the increasing correlation length of the non-affine displacement fields. As the system approaches the jamming transition, these displacements become increasingly correlated over larger distances, leading to greater variability in the shear modulus between different samples. This length scale l_c related to elastic-modulus fluctuations also manifests itself in vibrational modes^10,23,73 and responses to local perturbation,²² and these behaviours can also be captured within EMT.^49,74 The fact that χ_G ∼ δz^−1/2 ∼ l_c holds across different interaction potentials and spatial dimensions strongly suggests that the elastic criticality is robustly controlled by this mesoscopic length scale, regardless of local microscopic details. By contrast, the physical implication of the length scale associated with contact-number fluctuations remains unclear. However, our results highlight the fact that those fluctuations obey different scaling exponents in two and three dimensions, while the average contact number follows a single exponent. This suggests that the spatial correlation of the contact network or geometric constraints may involve different characteristic features in dimensions.

Finally, we situate our results on elastic fluctuations in the context of the HET, which provides a theoretical framework for describing sound attenuation in amorphous solids. According to HET, the attenuation coefficient Γ and the propagation frequency Ω are related through the Rayleigh scaling³⁵

where the disorder parameter γ is determined by spatial fluctuations of the local elastic modulus. Here, ω₀ denotes the characteristic frequency associated with the elastic correlation length.^75,76 Several specific forms of ω₀ have been proposed, some taking it as the boson peak frequency ω_BP,⁵¹ or evaluating it as ω_G = c_T/a₀,⁵⁰ where c_T is the sound speed of transverse wave and a₀ = ρ^−1/d is the characteristic interparticle spacing. It should be noted that HET does not provide a fully quantitative prediction of the absolute magnitude of Γ, and a system-dependent rescaling factor is generally required to achieve quantitative agreement with numerical data.^50,77

Detailed numerical simulations have elucidated the jamming scaling laws governing sound attenuation.¹⁴ In the low-frequency regime of the Rayleigh scattering, it follows

for both two and three dimensions (d = 2,3), where B represents a constant of the attenuation strength. Here, ω* denotes the onset frequency of the plateau in the vibrational density of states, and it follows the same jamming scaling law as ω_BP: ω* ∼ ω_BP ∼ δz.^10,11 Consequently, if one adopts the HET prediction with ω₀ = ω_BP, the predicted prefactor γ becomes a constant γ ∼ δz⁰ in two and three dimensions. By contrast, if one instead chooses ω₀ = ω_G, whose jamming scaling is ω_G ∼ δz^1/2, the corresponding prefactor γ becomes dimension-dependent: γ ∼ δz⁻¹ in two dimensions and γ ∼ δz^−3/2 in three dimensions.

Recent simulations of soft-core and Lennard-Jones systems have demonstrated that this disorder parameter could be estimated from sample-to-sample fluctuations^50,51,78 by assuming the spatial correlations of local elastic moduli are short-ranged and saturate upon coarse-graining.^40,41,61,79 This approach also has the merit of sidestepping the direct evaluation of local elastic moduli,^41,47 and it is claimed that spatial fluctuations of shear modulus γ can be replaced by the sample-to-sample fluctuations χ_G².^50,51 Our present result of jammed packings shows that χ_G² ∼ δz⁻¹ regardless of the spatial dimension, and this is indeed aligned with numerical simulations of acoustic attenuation in two-dimensional packings¹⁴ if we select ω_G as ω₀ in eqn (13) (note that the jamming scaling of the disorder parameter γ varies on the selection of the normalisation frequency ω₀).

Overall, our findings highlight that jamming scaling laws play a central role in governing elastic fluctuations and their impact on the acoustic attenuation. Nonetheless, the numerical character of the disorder parameter employed in HET remains poorly characterised; clarifying which fluctuations it captures and their scaling near jamming is necessary before HET can be applied to the attenuation prefactor. To resolve these issues, systematic studies at lower pressures and larger system sizes, direct measurements of spatial correlations of local elastic moduli, and quantitative comparisons with acoustic scattering data will be crucial.^77,80–84 Pursuing these directions should clarify how universal the observed scalings are and under which physical conditions the HET-based description becomes predictive.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data relevant to this work can be accessed at the University of Osaka Institutional Knowledge Archive.⁸⁵

Appendix

A Enumeration of the number of packings

In this appendix, we present the sizes of packing ensembles that are used in the present study. Table 1 shows the numbers of three-dimensional packings, whereas Table 2 shows the numbers of two-dimensional packings. All cases of the number of particles N, pressures p, and interaction potentials are presented in the two tables.

Table 1 The number of samples of three-dimensional packings. The numbers are presented for each number of particles N and each pressure p

N/p	10^−2	10^−3	10^−4	10^−5	10^−6	10^−7	10^−8
Harmonic
1000	11 958	11 825	11 472	10 693	10 691	—	—
4000	10 192	10 157	10 054	9832	9820	—	—
8000	5800	5796	5751	5700	5700	—	—
32 000	1530	1530	1527	1513	1510	—	—
64 000	690	690	688	681	681	—	—
128 000	407	407	407	407	405	—	—
Hertzian
2000	8047	8044	8014	7943	7987	7960	7935
32 000	1800	1800	1799	1794	1767	1742	—

---

Table 2 The number of samples of two-dimensional packings. The numbers are presented for each number of particles N and each pressure p

N/p	2 × 10^−2	10^−2	10^−3	10^−4	10^−5	10^−6	10^−7	10^−8
Harmonic
1000	10 945	10 556	10 312	9628	9145	9075	9126	—
4000	6013	6017	5981	5863	5510	4718	4551	—
8000	3709	3708	3704	3664	3724	3538	3421	—
32 000	1330	1330	1330	1329	1317	1298	1163	—
64 000	696	696	696	695	694	648	426	—
Hertzian
1000	7992	7986	7946	7820	7589	7135	6964	5699
32 000	—	576	576	576	574	575	556	289

---

Acknowledgements

KS acknowledges Yusuke Hara for useful discussions in the early stages of this work. This work is supported by JSPS KAKENHI Grant Number 25H01519 and JST ERATO Grant Number JPMJER2401.

Notes and references

1. A. J. Liu and S. R. Nagel, Annu. Rev. Condens. Matter Phys., 2010, 1, 347–369,  DOI:10.1146/annurev-conmatphys-070909-104045.
2. M. van Hecke, J. Phys.: Condens. Matter., 2010, 22, 033101,  DOI:10.1088/0953-8984/22/3/033101.
3. M. Wyart, S. R. Nagel and T. A. Witten, Eur. Lett., 2005, 72, 486–492,  DOI:10.1209/epl/i2005-10245-5.
4. M. Wyart, L. E. Silbert, S. R. Nagel and T. A. Witten, Phys. Rev. E, 2005, 72, 051306,  DOI:10.1103/PhysRevE.72.051306.
5. C. S. O’Hern, S. A. Langer, A. J. Liu and S. R. Nagel, Phys. Rev. Lett., 2002, 88, 075507,  DOI:10.1103/physrevlett.88.075507.
6. C. S. O’Hern, L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. E, 2003, 68, 011306,  DOI:10.1103/physreve.68.011306.
7. T. G. Mason, J. Bibette and D. A. Weitz, Phys. Rev. Lett., 1995, 75, 2051–2054,  DOI:10.1103/PhysRevLett.75.2051.
8. T. S. Majmudar, M. Sperl, S. Luding and R. P. Behringer, Phys. Rev. Lett., 2007, 98, 058001,  DOI:10.1103/PhysRevLett.98.058001.
9. W. G. Ellenbroek, E. Somfai, M. van Hecke and W. van Saarloos, Phys. Rev. Lett., 2006, 97, 258001,  DOI:10.1103/physrevlett.97.258001.
10. L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. Lett., 2005, 95, 098301,  DOI:10.1103/physrevlett.95.098301.
11. H. Mizuno, H. Shiba and A. Ikeda, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, E9767–E9774,  DOI:10.1073/pnas.1709015114.
12. N. Xu, V. Vitelli, M. Wyart, A. J. Liu and S. R. Nagel, Phys. Rev. Lett., 2009, 102, 038001,  DOI:10.1103/physrevlett.102.038001.
13. V. Vitelli, N. Xu, M. Wyart, A. J. Liu and S. R. Nagel, Phys. Rev. E, 2010, 81, 021301,  DOI:10.1103/physreve.81.021301.
14. H. Mizuno and A. Ikeda, Phys. Rev. E, 2018, 98, 062612,  DOI:10.1103/physreve.98.062612.
15. P. Olsson and S. Teitel, Phys. Rev. Lett., 2007, 99, 178001,  DOI:10.1103/physrevlett.99.178001.
16. H. E. Stanley, Introduction to Phase Transitions and Critical Phenomena, Oxford University Press, 1971.
17. J. J. Binney, N. J. Dowrick, A. J. Fisher and M. E. J. Newman, The Theory of Critical Phenomena: An Introduction to the Renormalization Group, Oxford University Press, 1992.
18. N. Goldenfeld, Lectures on Phase Transitions and the Renormalization Group, CRC Press, 1992.
19. H. Nishimori and G. Ortiz, Elements of Phase Transitions and Critical Phenomena, Oxford University Press, 2010.
20. C. P. Goodrich, A. J. Liu and S. R. Nagel, Phys. Rev. Lett., 2012, 109, 095704,  DOI:10.1103/physrevlett.109.095704.
21. C. P. Goodrich, S. Dagois-Bohy, B. P. Tighe, M. van Hecke, A. J. Liu and S. R. Nagel, Phys. Rev. E, 2014, 90, 022138,  DOI:10.1103/physreve.90.022138.
22. E. Lerner, E. DeGiuli, G. Düring and M. Wyart, Soft Matter, 2014, 10, 5085–5092,  10.1039/c4sm00311j.
23. L. Yan, E. DeGiuli and M. Wyart, Europhys. Lett., 2016, 114, 26003,  DOI:10.1209/0295-5075/114/26003.
24. C. P. Goodrich, A. J. Liu and J. P. Sethna, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 9745–9750,  DOI:10.1073/pnas.1601858113.
25. D. Hexner, A. J. Liu and S. R. Nagel, Phys. Rev. Lett., 2018, 121, 115501,  DOI:10.1103/physrevlett.121.115501.
26. D. Hexner, P. Urbani and F. Zamponi, Phys. Rev. Lett., 2019, 123, 068003,  DOI:10.1103/physrevlett.123.068003.
27. H. Ikeda, J. Chem. Phys., 2023, 158, 056101,  DOI:10.1063/5.0127064.
28. J. A. Giannini, E. Lerner, F. Zamponi and M. L. Manning, J. Chem. Phys., 2024, 160, 034502,  DOI:10.1063/5.0176713.
29. V. L. Gurevich, D. A. Parshin and H. R. Schober, Phys. Rev. B, 2003, 67, 094203,  DOI:10.1103/physrevb.67.094203.
30. S. Franz, G. Parisi, P. Urbani and F. Zamponi, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 14539–14544,  DOI:10.1073/pnas.1511134112.
31. T. Damart, A. Tanguy and D. Rodney, Phys. Rev. B, 2017, 95, 054203,  DOI:10.1103/physrevb.95.054203.
32. M. Baggioli and A. Zaccone, J. Phys.: Condens. Matter, 2022, 34, 215401,  DOI:10.1088/1361-648X/ac5d8b.
33. F. Vogel and M. Fuchs, Phys. Rev. Lett., 2023, 130, 236101,  DOI:10.1103/PhysRevLett.130.236101.
34. F. Vogel, P. Baumgärtel and M. Fuchs, Phys. Rev. X, 2025, 15, 011030,  DOI:10.1103/PhysRevX.15.011030.
35. W. Schirmacher and G. Ruocco, Heterogeneous Elasticity: The Tale of the Boson Peak, in Low-TemperatureThermal and Vibrational Properties of Disordered Solids, ed. M. A. Ramos World Scientific (Europe), 2022, ch. 9, pp. 331–373 DOI:10.1142/9781800612587_0009.
36. W. Schirmacher, Europhys. Lett., 2006, 73, 892–898,  DOI:10.1209/epl/i2005-10471-9.
37. W. Schirmacher, G. Ruocco and T. Scopigno, Phys. Rev. Lett., 2007, 98, 025501,  DOI:10.1103/physrevlett.98.025501.
38. A. Marruzzo, W. Schirmacher, A. Fratalocchi and G. Ruocco, Sci. Rep., 2013, 3, 1407,  DOI:10.1038/srep01407.
39. W. Schirmacher, T. Scopigno and G. Ruocco, J. Non-Cryst. Solids, 2015, 407, 133–140,  DOI:10.1016/j.jnoncrysol.2014.09.054.
40. M. Tsamados, A. Tanguy, C. Goldenberg and J.-L. Barrat, Phys. Rev. E, 2009, 80, 026112,  DOI:10.1103/PhysRevE.80.026112.
41. H. Mizuno, S. Mossa and J.-L. Barrat, Phys. Rev. E, 2013, 87, 042306,  DOI:10.1103/physreve.87.042306.
42. H. Mizuno, L. E. Silbert and M. Sperl, Phys. Rev. Lett., 2016, 116, 068302,  DOI:10.1103/physrevlett.116.068302.
43. G. Baldi, V. M. Giordano, G. Monaco and B. Ruta, Phys. Rev. Lett., 2010, 104, 195501,  DOI:10.1103/physrevlett.104.195501.
44. G. Baldi, V. M. Giordano, B. Ruta, R. Dal Maschio, A. Fontana and G. Monaco, Phys. Rev. Lett., 2014, 112, 125502,  DOI:10.1103/physrevlett.112.125502.
45. H. Mizuno, S. Mossa and J.-L. Barrat, Phys. Rev. B, 2016, 94, 144303,  DOI:10.1103/physrevb.94.144303.
46. K. Shiraishi, H. Mizuno and A. Ikeda, J. Chem. Phys., 2023, 158, 174502,  DOI:10.1063/5.0142648.
47. H. Mizuno, S. Mossa and J.-L. Barrat, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 11949–11954,  DOI:10.1073/pnas.1409490111.
48. L. Wang, L. Berthier, E. Flenner, P. Guan and G. Szamel, Soft Matter, 2019, 15, 7018–7025,  10.1039/c9sm01092k.
49. E. DeGiuli, A. Laversanne-Finot, G. Düring, E. Lerner and M. Wyart, Soft Matter, 2014, 10, 5628–5644,  10.1039/c4sm00561a.
50. G. Kapteijns, D. Richard, E. Bouchbinder and E. Lerner, J. Chem. Phys., 2021, 154, 081101,  DOI:10.1063/5.0038710.
51. S. Mahajan and M. Pica Ciamarra, Phys. Rev. Lett., 2021, 127, 215504,  DOI:10.1103/PhysRevLett.127.215504.
52. J. Guénolé, W. G. Nöhring, A. Vaid, F. Houllé, Z. Xie, A. Prakash and E. Bitzek, Comput. Mater. Sci., 2020, 175, 109584,  DOI:10.1016/j.commatsci.2020.109584.
53. C. P. Goodrich, PhD thesis, University of Pennsylvania, 2015, https://repository.upenn.edu/handle/20.500.14332/28587.
54. S. Dagois-Bohy, B. P. Tighe, J. Simon, S. Henkes and M. van Hecke, Phys. Rev. Lett., 2012, 109, 095703,  DOI:10.1103/physrevlett.109.095703.
55. P. Charbonneau, E. I. Corwin, G. Parisi and F. Zamponi, Phys. Rev. Lett., 2015, 114, 125504,  DOI:10.1103/physrevlett.114.125504.
56. C. Maloney and A. Lemaître, Phys. Rev. Lett., 2004, 93, 195501,  DOI:10.1103/physrevlett.93.195501.
57. C. E. Maloney and A. Lemaître, Phys. Rev. E, 2006, 74, 016118,  DOI:10.1103/physreve.74.016118.
58. A. Lemaître and C. Maloney, J. Stat. Phys., 2006, 123, 415–453,  DOI:10.1007/s10955-005-9015-5.
59. S. Alexander, Phys. Rep., 1998, 296, 65–236,  DOI:10.1016/S0370-1573(97)00069-0.
60. N. W. Ashcroft and N. D. Mermin, Solid State Physics, Saunders Collage Publishing, 1976.
61. H. Mizuno, K. Saitoh and L. E. Silbert, Phys. Rev. E, 2016, 93, 062905,  DOI:10.1103/physreve.93.062905.
62. J. F. Lutsko, J. Appl. Phys., 1989, 65, 2991,  DOI:10.1063/1.342716.
63. T. H. K. Barron and M. L. Klein, Proc. Phys. Soc., 1965, 85, 523,  DOI:10.1088/0370-1328/85/3/313.
64. J. P. Wittmer, H. Xu, P. Polińska, C. Gillig, J. Helfferich, F. Weysser and J. Baschnagel, Eur. Phys. J. E, 2013, 36, 131,  DOI:10.1140/epje/i2013-13131-y.
65. J. P. Wittmer, H. Xu, P. Polińska, F. Weysser and J. Baschnagel, J. Chem. Phys., 2013, 138, 12A533,  DOI:10.1063/1.4790137.
66. A. Zaccone and E. Scossa-Romano, Phys. Rev. B, 2011, 83, 184205,  DOI:10.1103/physrevb.83.184205.
67. Y. Hara, H. Mizuno and A. Ikeda, Soft Matter, 2023, 19, 6046–6056,  10.1039/d3sm00566f.
68. W. G. Ellenbroek, M. van Hecke and W. van Saarloos, Phys. Rev. E, 2009, 80, 061307,  DOI:10.1103/physreve.80.061307.
69. H. A. Makse, N. Gland, D. L. Johnson and L. M. Schwartz, Phys. Rev. Lett., 1999, 83, 5070–5073,  DOI:10.1103/PhysRevLett.83.5070.
70. P. Wang, S. Zhang, P. Tuckman, N. T. Ouellette, M. D. Shattuck and C. S. O’Hern, Phys. Rev. E, 2021, 103, 022902,  DOI:10.1103/physreve.103.022902.
71. G. Kapteijns, E. Bouchbinder and E. Lerner, Phys. Rev. E, 2021, 104, 035001,  DOI:10.1103/physreve.104.035001.
72. G. Düring, E. Lerner and M. Wyart, Soft Matter, 2013, 9, 146–154,  10.1039/c2sm25878a.
73. M. Shimada, H. Mizuno, M. Wyart and A. Ikeda, Phys. Rev. E, 2018, 98, 060901,  DOI:10.1103/physreve.98.060901.
74. M. Wyart, Europhys. Lett., 2010, 89, 64001,  DOI:10.1209/0295-5075/89/64001.
75. W. Schirmacher, C. Tomaras, B. Schmid, G. Baldi, G. Viliani, G. Ruocco and T. Scopigno, Condens. Matter Phys., 2010, 13, 23606,  DOI:10.5488/cmp.13.23606.
76. W. Schirmacher, J. Non-Cryst. Solids, 2011, 357, 518–523,  DOI:10.1016/j.jnoncrysol.2010.07.052.
77. C. Caroli and A. Lemaître, Phys. Rev. Lett., 2019, 123, 055501,  DOI:10.1103/PhysRevLett.123.055501.
78. G. Szamel and E. Flenner, J. Chem. Phys., 2025, 163, 050903,  DOI:10.1063/5.0280663.
79. A. Shakerpoor, E. Flenner and G. Szamel, Soft Matter, 2020, 16, 914–920,  10.1039/c9sm02022e.
80. C. Caroli and A. Lemaître, J. Chem. Phys., 2020, 153, 144502,  DOI:10.1063/5.0019964.
81. G. Szamel and E. Flenner, J. Chem. Phys., 2022, 156, 144502,  DOI:10.1063/5.0085199.
82. S. Mahajan and M. Pica Ciamarra, Phys. Rev. E, 2024, 109, 024605,  DOI:10.1103/PhysRevE.109.024605.
83. E. Flenner and G. Szamel, Sci. Adv., 2025, 11, eadu6097,  DOI:10.1126/sciadv.adu6097.
84. E. Flenner and G. Szamel, J. Phys. Chem. B, 2025, 129, 1855–1863,  DOI:10.1021/acs.jpcb.4c07545.
85. K. Shiraishi and H. Mizuno, Dataset: Critical fluctuations of elastic moduli in jammed solids, 2025 DOI:10.60574/103296.

---