Our numerical model of granular materials is the aggregate of two-dimensional disks. We introduce the force between the two disks, i and j, in contact as the sum of elastic and viscous forces. The elastic force consists of normal and tangential parts as k _n ξ _ij and k t ξ i j ⊥ , respectively, where k _n (k _t ) is the normal (tangential) stiffness. Here, ξ _ij ≡ R _i + R _j − r _ij > 0 with the disk radii, R _i and R _j , and center-to-center distance, r _ij , represents an overlap between the particles, while ξ i j ⊥ is the relative tangential displacement at contact. On the other hand, the viscous force consists of normal and tangential parts as η n ξ ̇ i j and η t ξ ̇ i j ⊥ , respectively, where η _n (η _t ) is the viscosity in the normal (tangential) direction.

In our numerical model, the strength of contact forces is determined by the microscopic stiffness and viscosity, i.e., k _n , k _t , η _n , and η _t . To control the strength of contact forces, we introduce the following dimensionless parameters [13], ρ K ≡ k t k n , ρ D ≡ η t η n , ϵ n ≡ η n m k n . (1) Here, the stiffness ratio ρ _K quantifies the strength of tangential elastic forces, while the damping ratio ρ _D represents the relative magnitude of tangential viscous forces. In addition, the damping factor ϵ _n , together with ρ _D , controls the strength of energy dissipation.

The dimensionless parameters, Eq. 1, were suggested by Kruyt [13] to quantify the influence of microscopic properties on the sound in granular crystals. In Ref. [36], we have studied the role of ρ _K in the sound characteristics of two-dimensional “disordered” granular disks, where the contact damping was absent, i.e., ρ _D = ϵ _n = 0. In this paper, we will focus on the effects of ρ _D and ϵ _n on the sound properties of disordered granular disks. Note that we do not introduce the tangential elastic force and viscous forces, i.e., ρ _K = ρ _D = ϵ _n = 0, when we prepare disordered configurations of the disks (Section 2.2). However, we introduce these forces when we simulate small oscillations of the disks around their equilibrium positions, where ρ _K is fixed to unity (Section 2.3).

We prepare disordered configurations of granular disks by the same method as in Refs. [36, 40]. Our system is a 50:50 binary mixture of N = 32,768 disks, where every disk has the same mass m and different diameters, d _S and d _L = 1.4d _S . A repulsive force between the disks, i and j, in contact is given by the elastic force in normal direction, f _ij = k _n ξ _ij . We randomly distribute the N disks in a L × L square periodic box and fully minimize elastic energy of the system with the FIRE algorithm [41]. We stop the energy minimization if the maximum acceleration of the disks becomes less than 1 0 − 9 d 0 / t 0 2 with the mean disk diameter, d ₀ ≡ (d _S + d _L )/2, and time unit, t 0 ≡ m / k n . In the following, we scale every length and time by d ₀ and t ₀, respectively. Since our system is bi-dispersed, disk positions after the energy minimization are disordered, where we denote their disordered configurations as { r _i (0)} (i = 1, …, N). Note that the packing fraction of the disks is fixed to 0.9 and thus our system is far above the jamming transition [42–44]. Nevertheless, there are typically 0.2% of rattlers that do not contribute to the mechanical contact network.

To simulate sound in the disordered granular media, we introduce linear equations of motion of the granular disks [45, 46]. Now, we introduce the tangential elastic force and viscous forces between the disks in contact. Because the disordered configurations, { r _i (0)} (Section 2.2), are mechanically stable, we describe small oscillations of the disks around { r _i (0)} by the following equation [36, 40, 47], M q ̈ t = − D q t − B q ̇ t . (2) Here, t denotes time and q t ≡ u i t , θ i t T (3) is a 3N-dimensional displacement vector which consists of translational displacement in the xy-plane, u _i (t) ≡ r _i (t) − r _i (0), and angular displacement of each disk, θ _i (t). In the linear equations of motion (Eq. 2), M , D , and B are 3N × 3N mass matrix, dynamical matrix, and damping matrix, respectively [36, 40]. The mass matrix is diagonal and consists of mass and moment of inertia of each disk. The dynamical matrix is defined by second derivatives of elastic energy [47–51], whereas the damping matrix is given by second derivatives of dissipation function [40, 52]. In SM, we derive Eq. 2 and show explicit forms of M , D , and B .

In Eq. 2, the elastic forces between the disks in contact are given by − D q ( t ) . To calculate each element of D , we define the elastic energy as the sum of harmonic potentials stored in normal and tangential directions (see SM). Note that the initial configurations, { r _i (0)}, are mechanically stable even if we introduce the tangential component of elastic energy [36]. The normal (tangential) component of elastic energy is characterized by the normal (tangential) stiffness, k _n (k _t ). On the other hand, the viscous forces between the disks in contact are given by − B q ̇ ( t ) , where q ̇ t ≡ u ̇ i t , θ ̇ i t T (4) is the time derivative of the 3N-dimensional displacement vector, q (t). Each element of B is defined by the dissipation function [40], which is also decomposed into normal and tangential components (see SM). As in the case of the elastic energy, the normal (tangential) component of the dissipation function is characterized by the normal (tangential) viscosity, η _n (η _t ).

The linear equations of motion (Eq. 2) are equivalent to the so-called “spring-dashpot model”, i.e., a canonical model of granular materials for DEM simulations [36]. However, the matrices, D and B , are given by the initial equilibrium positions, { r _i (0)}, such that the interactions, − D q ( t ) and − B q ̇ ( t ) , are calculated based on the initial contact network. Since the displacements, q (t), are so small, such harmonic approximations of the elastic energy and dissipation function are valid (SM) [45, 46]. Note that the static friction between the disks in contact is modeled by the tangential elastic force, k t ξ i j ⊥ . The static friction coefficient is infinite, i.e., the dynamical (Coulomb) friction is not implemented, because our dynamical matrix is given by the elastic energies and cannot describe the dynamical friction.

To examine sound properties of the granular disks, we employ a similar method as in Refs. [36, 40, 47]. We numerically integrate the linear equations of motion, Eq. 2, under periodic boundary conditions. Initial velocities of the disks are given by a sinusoidal standing wave, q ̇ 0 = A ⁡ sin k ⋅ r i 0 T , (5) where A and k are amplitude and wave vectors, respectively [37, 40, 47]. As shown in SM, we use different amplitude and wave vectors to demonstrate three different types of elastic wave, i.e., longitudinal (L), transverse (T), and rotational (R) modes. The latter (R mode) represents micropolar rotations of the disks [36], which are not relevant in frictionless systems [40, 47].

By using numerical solutions of the linear equations of motion (Eq. 2), we visualize time evolution of the sinusoidal standing wave. Figure 1 displays snapshots of our numerical simulation at t/t ₀ = (A) 0, (B) 2, (C) 4, and (D) 6. Each disk (circle) is colored according to its angular velocity, θ ̇ i ( t ) (i = 1, …, N), where θ ̇ i ( t ) increases from −A _θ (blue) to A _θ (red). In this figure, the wave vector is k = (k, 0) (as indicated by the arrow) with the wave number, k ≃ 0.29 d 0 − 1 , where the dimensionless parameters are given by ρ _K = 1, ρ _D = 0.2, and ϵ _n = 0.1. As can be seen, the initial standing wave (Figure 1A) is attenuated with time (Figures 1B, C) and eventually vanishes in a long time limit (Figure 1D). Such the wave attenuation is caused not only by scattering (due to the disordered configuration of the disks) [36, 47] but also by energy dissipation (due to the viscous forces between the disks in contact) [40]. We can also observe similar wave attenuation when we visualize the time evolution of translational velocities of the disks, u ̇ i ( t ) (i = 1, …, N) (data are not shown).

From numerical solutions of Eq. 2, we obtain the data of disk velocities, q ̇ ( t ) (Eq. 4). We apply Fourier transforms to the velocities as u ̇ k t , θ ̇ k t = ∑ i = 1 N u ̇ i t , θ ̇ i t e − I k ⋅ r i t (6) with the imaginary unit I, where the disk position r _i (t) is also obtained from the numerical solutions of Eq. 2. Then, we introduce the L and T modes as u ̇ k ‖ t ≡ u ̇ k t ⋅ k ̂ k ̂ , (7) u ̇ k ⊥ t ≡ u ̇ k t − u ̇ k ‖ t , (8) respectively, where k ̂ ≡ k / k is a unit vector parallel to the wave vector.

The normalized VAFs of L, T, and R modes are defined as C L k , t = 〈 u ̇ k ‖ t ⋅ u ̇ − k ‖ 0 〉 〈 | u ̇ k ‖ 0 | 2 〉 , (9) C T k , t = 〈 u ̇ k ⊥ t ⋅ u ̇ − k ⊥ 0 〉 〈 | u ̇ k ⊥ 0 | 2 〉 , (10) C R k , t = 〈 θ ̇ k t θ ̇ − k 0 〉 〈 | θ ̇ k 0 | 2 〉 , (11) respectively. Figure 2 shows the time evolution of (A) C _L (k, t), (B) C _T (k, t), and (C) C _R (k, t), where the dimensionless parameters, ρ _K , ρ _D , and ϵ _n , are as in Figure 1. To calculate each VAF, we use two different wave vectors (as listed in Table 1 in the SM) and average each VAF over the two samples (wave vectors). As can be seen, the oscillation of L mode is faster than that of T mode (Figures 2A, B). Furthermore, the oscillation of R mode is much faster than those of L and T modes (note the different horizontal scale in Figure 2C), meaning that the sound speed of micropolar rotations of the disks is much larger than acoustic sound speeds. In addition, the R mode decays much faster than the acoustic L and T modes, implying stronger scattering and energy dissipation of rotational motions. Notice that the VAFs are entirely damped (“overdamped”) without oscillations if the parameters, ρ _D and ϵ _n , are sufficiently large.

To extract sound characteristics of the granular disks, we fit a damped oscillation, C α k , t = e − γ α k t ⁡ cos ω α k t , (12) to the data of normalized VAFs (α = L, T, R) [36, 40, 47]. In Eq. 12, the frequency ω _α (k) represents the dispersion relation, while the coefficient γ _α (k) quantifies the sound attenuation of each mode. For each wave number k, we adjust the fitting parameters, ω _α (k) and γ _α (k), to see perfect agreement between the data of normalized VAFs and damped oscillations. The solid lines in Figure 2 are the damped oscillations; Eq. 12, fitted to the data of normalized VAFs (circles). We also confirm perfect agreement between the data and Eq. 12 for all the dimensionless parameters, ρ _K , ρ _D , and ϵ _n , used in simulations (data are not shown).

When the VAFs are overdamped, the data of normalized VAFs cannot be described by the damped oscillation, Eq 12. Therefore, there is a need to introduce a criterion for the fitting parameters, ω _α (k) and γ _α (k). Here, we employ the Ioffe-Regel (IR) limit for the criterion, where Eq. 12 is meaningful only if the condition π γ α k ω α k < 1 (13) is satisfied [18]. The ratio πγ _α (k)/ω _α (k) is a monotonically increasing function of the wave number k. Hence Eq. 13 is equivalent to k < k α * , where the limit wave number k α * is defined as π γ α ( k α * ) / ω α ( k α * ) = 1 . In SM, we show our results of the limit wave number, k α * . The limit wave number is a monotonically decreasing function of the strength of contact damping, ρ _D and ϵ _n . In the following, we only show the results of ω _α (k) and γ _α (k) in k < k α * .

We analyze the dispersion relation of each mode, ω _α (k) (α = L, T, R), extracted from the data of VAFs and clarify its dependence on the contact damping. Figure 3A displays ω _α (k) as functions of the scaled wave number kd ₀, where we control the damping ratio ρ _D as listed in the legend (see SM for their dependence on the damping factor ϵ _n ). In this figure, the dispersion relations of L and T modes exhibit ordinary acoustic branches, where ω _L (k) and ω _T (k) increase from zero with the wave number. The dispersion relation of L mode is larger than that of T mode, i.e., ω _L (k) > ω _T (k), over the whole range of k. This means that the oscillation of the VAF of L mode is always faster than that of T mode (as shown in Figures 2A, B). In addition, these dispersion relations are quite insensitive to ρ _D (Figure 3A) and ϵ _n (SM). This trend is consistent with the theoretical prediction of granular crystals [13] though ω _T (k) is cut-off at the limit wave number k T * which decreases with the increase of ρ _D (see SM).

In contrast to the acoustic dispersion relations, the dispersion relation of the R mode exhibits a characteristic optical-like branch (Figure 3A) [7]. The influence of contact damping is significant at large wave number, where ω _R (k) dramatically decreases if we increase either ρ _D (inset to Figure 3A) or ϵ _n (SM). It is interesting that the viscous forces in normal and tangential directions, characterized by ϵ _n and ρ _D , respectively, have a similar effect on ω _R (k) because micropolar rotations of the disks are driven only by tangential forces. In addition, ω _R (k) in Figure 3A is not cut-off though ω _T (k) is at k T * . Therefore, compared with the T mode, the oscillation of R mode is long-lived, implying that micropolar rotations are not strongly coupled with transverse motions of the disks.

We examine sound attenuation of each mode by the attenuation coefficients, γ _α (k) (α = L, T, R), extracted from the data of VAFs. Figure 3B displays γ _α (k) as functions of the scaled wave number, where we vary ρ _D as listed in the legend (see SM for the effect of ϵ _n on γ _α (k)). In the continuum limit, k → 0, the attenuation coefficients of the L and T modes are quadratic in the wave number, i.e., γ _L (k), γ _T (k) ∼ k ², as indicated by the dashed line. The quadratic growth of the attenuation coefficients is typical of viscoelastic media [40] and is also predicted by the lattice theory of granular crystals [13]. In addition, regardless of the wave number, both γ _L (k) and γ _T (k) increase with the increase of strength of contact damping (see also SM).

In contrast, the attenuation coefficient of the R mode remains constant, γ _R (k) ∼const., in the continuum limit and is much larger than those of acoustic modes, i.e., γ _R (k) ≫ γ _L (k), γ _T (k), for small wave numbers. Therefore, the decay of the normalized VAF of R mode is much faster than those of L and T modes (as shown in Figure 2). However, if the damping ratio is relatively small (ρ _D < 0.4), the attenuation coefficient of the R mode becomes smaller than those of acoustic modes, i.e., γ _R (k) < γ _L (k), γ _T (k), at large wave number. We also observe this phenomenon if the damping factor ϵ _n is small enough (SM). Because the lattice theory of granular crystals implies that the attenuation of R mode is always stronger than those of acoustic modes [13], our results suggest that the weak attenuation of R mode at large wave number is specific to disordered granular disks. Similar to the acoustic modes, γ _R (k) significantly increases with the increase of strength of contact damping regardless of the wave number (see also SM).

We quantify sound speed of each mode by phase speed defined as c _α (k) ≡ ω _α (k)/k (α = L, T, R). Figure 3C shows c _L (k) and c _T (k) as functions of the scaled wave number, where we vary the parameter ρ _D as listed in the legend (see SM for different values of ϵ _n ). As can be seen, the phase speeds of acoustic L and T modes converge to constants (horizontal lines) in the continuum limit, k → 0. The continuum limit, c _α (0) (α = L, T), is insensitive to the strength of contact damping, ρ _D and ϵ _n . Therefore, the viscous forces between the disks in contact do not affect macroscopic speeds of sound [40].

The phase speed of L mode, c _L (k), is a monotonically decreasing function of the wave number, meaning that the dispersion relation, ω _L (k), becomes sub-linear at large wave number (see Figure 3A). The phase speed of T mode, c _T (k), also decreases from the continuum limit, c _T (0), when the wave number k increases from zero. However, further increasing k, we observe that c _T (k) starts increasing and generates a small “dip” at intermediate wave number (as indicated by the arrow in Figure 3C). Such a small dip in the phase speed is characteristic of (energy conserving) disordered media [47, 53] and has been considered to be a sign of the boson peak in vibrational density of states [54–58]. It is believed that the boson peak is a consequence of elastic heterogeneities in disordered materials [59–61]. Therefore, the small dip in c _T (k) is unique to our study on disordered granular disks, i.e., is not expected to exist in c _T (k) of granular crystals [7, 10–13]. Note that c _L (k) exhibits no dip, which is in sharp contrast with the results of disordered frictionless disks [40]. Moreover, c _L (k) is relatively insensitive to the strength of contact damping, ρ _D and ϵ _n . Similarly, c _T (k) is not affected by the damping ratio ρ _D , though it is cut-off at k T * which decreases with the increase of ρ _D (see SM). In addition, the magnitude of small dip in c _T (k) decreases with the increase of damping factor ϵ _n (see SM) as previously found in the model of disordered frictionless disks [40].

As shown in Figure 3D, the phase speed of the R mode, c _R (k), exhibits asymptotic behavior, i.e., c _R (k) ∼ k ⁻¹ (dashed lines), in the continuum limit, k → 0. Thus, the optical-like branch in the dispersion relations becomes flat, ω _R (k) ∼const., when the wave number is small enough, Figure 3A. The influence of contact damping is significant at large wave number, where c _R (k) is lowered with the increase of contact damping (see the inset).

We compare our numerical results with theoretical predictions of granular crystals. In particular, we focus on the dependence of sound characteristics on the strength of contact damping, ρ _D and ϵ _n .

The lattice theory of granular crystals [13] predicts that, in the continuum limit, k → 0, the dispersion relations of the L and T modes depend on the stiffness ratio as ω L ( k ) t 0 ≃ 1 + ρ K k d 0 and ω T ( k ) t 0 ≃ ( 1 + ρ K ) / 2 k d 0 , respectively. Our results of ω _L (k) and ω _T (k) are consistent with this prediction, except for 25%–30% smaller prefactors. In SM, we show our results of ω _L (k) and ω _T (k) at the smallest wave number, Δk ≡ 2π/L, where both are insensitive to ρ _D and ϵ _n , i.e., as in Figure 3A, and are described as ω L ( Δ k ) t 0 = ω L d 1 + ρ K Δ k d 0 and ω T ( Δ k ) t 0 = ω T d ( 1 + ρ K ) / 2 Δ k d 0 with the prefactors, ω L d ≃ 0.741 and ω T d ≃ 0.712 .

The lattice theory also predicts that, in the continuum limit, the dispersion relation of the R mode is controlled by the dimensionless parameters according to ω R k ∝ ρ K − a 0 ρ D ϵ n 2 , (14) where a ₀ is a dimensionless constant [13]. To examine the theoretical prediction, we plot our numerical results of the dispersion relation at the smallest wave number, ω _R (Δk), as a function of the damping ratio ρ _D in Figure 4A. In this figure, the dashed line is the theoretical prediction; Eq. 14, approximated to the data of ω _R (Δk) for ϵ _n = 0.02 (where a ₀ ≃ 1.3 × 10⁻² with a prefactor, 4.59). It is apparent that our results are qualitatively different from the theoretical prediction; Eq. 14 is convex upward, while the data are convex downward. Our findings suggest that the influence of disordered disk configurations on the R mode cannot be neglected even in the continuum limit.

In the continuum limit, the attenuation coefficients of the L and T modes are predicted to be quadratic in the wave number as γ L ( k ) t 0 ≃ ϵ n ( 1 + 4 ρ D ) ( k d 0 ) 2 and γ T ( k ) t 0 ≃ ϵ n { ( 1 + ρ D ) / 2 } ( k d 0 ) 2 , respectively [13]. Furthermore, the attenuation coefficient of the R mode is predicted to be proportional to the strength of contact damping as γ _R (k) ∝ϵ _n ρ _D in the continuum limit [13]. Therefore, all the attenuation coefficients predicted by the lattice theory are proportional to ϵ _n and are linear in ρ _D , meaning that the sound is not attenuated if the contact damping is absent, i.e., γ _α (k) = 0 if ϵ _n = ρ _D = 0 (α = L, T, R). However, in disordered media, the sound is also attenuated by scattering even if the contact damping does not exist, i.e., γ _α (k) > 0 even if ϵ _n = ρ _D = 0 [40]. The scattering of the acoustic L and T modes is represented by the scaling, γ _α (k) ∝ k ³, i.e., like Rayleigh scattering in two dimensions [1]. Assuming that the attenuation coefficient of the R mode is also finite, γ _R (k) > 0, in the limit, ϵ _n = ρ _D = 0, we modify the theoretical predictions to take the influence of disorder into account as γ L k t 0 = γ L d ϵ n 1 + r L ρ D k d 0 2 + q L k d 0 3 , (15) γ T k t 0 = γ T d ϵ n 1 + r T ρ D k d 0 2 + q T k d 0 3 , (16) γ R k t 0 = γ R d ϵ n 1 + r R ρ D + q R , (17) where γ α d , r _α , and q _α (α = L, T, R) are introduced as fitting parameters.

Adjusting the parameters, we see good agreements between Eqs. 15–17 and numerical results. Figures 4B–D show dimensionless functions, (B) f L ( k ) ≡ γ L ( k ) t 0 / { ϵ n ( k d 0 ) 2 } − q L k d 0 / ϵ n − γ L d , (C) f T ( k ) ≡ γ T ( k ) t 0 / { ϵ n ( k d 0 ) 2 } − q T k d 0 / ϵ n − γ T d , and (D) f R ( k ) ≡ γ R ( k ) t 0 / ϵ n − q R / ϵ n − γ R d , at the smallest wave number, Δk. In these figures, the dashed lines represent (B) γ L d r L ρ D , (C) γ T d r T ρ D , and (D) γ R d r R ρ D , indicating Eqs 15–17. As can be seen, all the numerical results are nicely collapsed on the dashed lines though the data of f _R (Δk) in (D) considerably deviate up to 8.1% from γ R d r R ρ D (dashed line) as the strength of contact damping, ρ _D and ϵ _n , increases. The dimensionless parameters in Eqs. 15–17 are given by γ L d ≃ 0.440 , r _L ≃ 0.259, q _L ≃ 0.104, γ T d ≃ 0.155 , r _T ≃ 0.700, q _T ≃ 0.061, γ R d ≃ 2.16 × 1 0 − 2 , r _R ≃ 5.27 × 10², and q _R ≃ 3.72 × 10⁻². This means that, even if we modify the constants in the theoretical predictions, our numerical results cannot be explained. Therefore, structural disorder significantly alters the sound attenuation in granular media and the improved damping relations, Eqs. 15–17, have to be explained by advanced theory in future.