Effect of Disorder on Bulk Sound Wave Speed : A Multiscale Spectral Analysis

Disorder of size (polydispersity) and mass of discrete elements/particles in randomly structured media (e.g. granular matter like soil) has numerous effects on the materials’ sound propagation characteristics. The influence of disorder on energy and momentum transport, the sound wave speed and its low pass frequency filtering characteristics is the subject of this study. Goal is understanding the connection between the particle-micro-scale disorder and dynamics and the system-macro-scale wave propagation, which can be applied to non-destructive testing, seismic exploration of buried objects (oil, mineral, etc.) 5 or to study the internal structure of the Earth. To isolate the longitudinal P-wave mode from shear and rotational modes, a one-dimensional system of equal size elements/particles is used to study the effect of mass disorder alone via (direct and/or ensemble averaged) real time signals, signals in Fourier space, energy and dispersion curves. Increase in mass disorder (where disorder has been defined such that it is independent of the shape of the probability distribution of masses) decreases the sound wave speed along a granular chain. Energies associated with the eigenmodes are conserved, independent of time, and can 10 be used to obtain better quality dispersion relations for disordered chains; these dispersion relations confirm the decrease in pass-frequency and wave speed with increasing disorder acting opposite to the wave acceleration close to the source.


Introduction
Sound wave propagation through matter has been an extensive area of research (for a textbook example, see Aki and Richards, 2002) being applied to study earthquakes or the internal structure of the Earth, as well as oil, gas or mineral ex-ploration (seismology).Waves can be used for dissecting the human body without using blades, revealing material properties through nondestructive testing (ultrasonics), studying the structure of lattices or designing metamaterials.There are numerous applications and uncountable problems which still need to be solved, where the challenge has always been resolving the finest structures of matter using wave propagation.Hence, steps are being taken in the direction of micromechanics of seismic waves; see, e.g., O'Donovan et al. (2016).
Disordered, heterogeneous and random media cause multiple scattering of seismic waves that are dispersed, attenuated and localized in space (Sato, 2011;Scales and Van Vleck, 1997).The phenomenon of multiple scattering causes the formation of the so-called "coda", which is the tail of a propagating wave pulse.While coda was earlier treated as noise (Weaver, 2005), now it has given way to coda wave interferometry with multiple applications (Snieder et al., 2002).The coda has been studied in detail in laboratory experiments with uniaxial or triaxial devices, e.g., pulse propagation across glass beads (Jia et al., 1999) and sintered glass beads (Güven, 2016), indicating extreme sensitivity towards system preparation and configuration and getting washed out on ensemble averaging with only the coherent part of the signal remaining.In Van Der Baan (2001), it was shown that macroscopic or seismic waves governed by the classical wave equation did not exhibit localization at lower frequencies, but this idea got repudiated by Larose et al. (2004), where weak localization (a mesoscopic phenomenon, precursor to wave localization; Sheng, 2006) was experimentally observed at frequencies as low as 20 Hz, indicating the inadequacy of the classical wave equation.
In recent years, wave propagation through granular materials has attracted a lot of attention.Granular material is heterogeneous with many discretized units and can be used for other modeling geometrically heterogeneous media (Matsuyama and Katsuragi, 2014).The studies done using ordered and disordered lattices for wave propagation (Gilles and Coste, 2003;Coste and Gilles, 2008, etc.) have helped the understanding of wave propagation in granular materials through dispersion relations, frequency filtering, etc. Scaling laws allow the relation of various physical parameters such as density, pressure, coordination number, etc., with the moduli, forming an effective medium theory (EMT) for granular matter (Makse et al., 2004).Nesterenko (1983) showed the existence of localized wave packets propagating in a nonlinear granular chain (onedimensional granular material) under the condition of "sonic vacuum" (in the limit of zero acoustic wave speed and vanishing confining pressure), thus forming supersonic solitary waves; such concepts have been exploited immensely to develop various kinds of metamaterials such as those used for shock and energy trapping (Daraio et al., 2006), an acoustic diode (Boechler et al., 2011) or for understanding and studying jamming transitions in granular matter (van den Wildenberg et al., 2013;Upadhyaya et al., 2014).Some of the open questions and developments related to wave propagation in unconsolidated granular matter, such as higher harmonics generation, nonlinear multiple scattering, soft modes, rotational modes, etc., have been addressed by Tournat and Gusev (2010).However, in the following, the focus of attention will not be on solitons and unconsolidated granular matter; hence, there will be no sonic vacuum during analyses (no opening and closing of contacts of particles).
A striking characteristic of consolidated granular matter is that grain-grain forces are arranged and correlated in a linear manner known as force chains (Somfai et al., 2005).Similar to the force chains, Pasternak et al. (2015) showed the existence of moment chains in granular media, i.e., correlations of grain-grain mutual rotations.These chains are mesoscopic structures and are just one of the many microrotational effects of granules.Cosserat continuum theory can be used to model these micropolar and/or microrotational effects, as discussed in detail by Pasternak and Mühlhaus (2005).
The force chains and granular chains which carry the large forces of the system supposedly support faster sound transmission across granular matter (Ostojic et al., 2006).In Owens and Daniels (2011), it was seen from experiments with two-dimensional photo-elastic disks that vibration propagates along the granular chains, visualized by the brightness due to compression between the particles; however, the exact mechanisms of propagation of the vibrations are still a matter of ongoing research.Our system under investigation will be a single one of such granular chains (Fig. 1); it will assist in isolating the P wave or the longitudinal excitation from all other kinds of excitations (S wave, rotational wave, etc).In Merkel et al. (2010) it was seen that inclusion or removal of rotation does not significantly affect the longitudinal mode in an ordered granular crystal.However, the situation is different when rotations become prominent and other wave modes cannot be ignored (see Yang and Sutton, 2015;Merkel and Luding, 2017, and the references therein).
Although it is very simplistic, a polydisperse granular chain can have various kinds of disorder: size, mass or stiffness disorder (Lawney and Luding, 2014).The size or mass disorder has a much stronger contribution towards disorder than stiffness because mass ∝ radius 3 , whereas stiffness ∝ radius 1/3 (Achilleos et al., 2016).Hence, only mass disorder for the disordered granular chain has been chosen.However, there are processes when stiffness disorder cannot be ignored, for instance, the processes when the repulsive interaction force between the fragments or elements of the material being modeled has different stiffness during compression and tension (bilinear oscillator; Dyskin et al., 2014), infinite stiffness during compression (impact oscillator; Dyskin et al., 2012 andGuzek et al., 2016) or negative stiffness (Pasternak et al., 2014, andEsin et al., 2016).
In Sect. 2 an impulse propagating across a granular chain is modeled.A similar model was used in Marketos and O'Sullivan (2013), Lawney and Luding (2014) and Otsubo et al. (2017).Section 2.8 concerns the dispersion relation for wave propagation across a granular medium, Sect.2.9 concerns the group velocity and Sect.2.10 concerns a novel way of computing the dispersion relation in terms of moments of eigenmodal energy.In Sect.3, the quantities mentioned in Sect. 2 are computed and the observations are discussed.Section 4 summarizes and concludes the observations made in Sect. 3 with Sect. 2 as the foundation and an outlook of the ongoing as well as possible future research work on wave propagation in granular matter is given.
2 Modeling a general one-dimensional chain A one-dimensional chain of N + 2 particles is considered (Fig. 2).Each particle i has mass m (i) and contact stiffness κ (i,j ) with respect to a neighboring particle j .The tilde symbols are used for dimensional quantities.The interaction force experienced by neighboring particles i and j is At rest : During wave propagation : with the contact stiffness κ (i,j ) and the particle overlap |, with the radius r and coordinates x of the centers of the particles.The Hertzian and linear models are given by β = 1/2 and β = 0, respectively (Lawney and Luding, 2014).This force acting from , corresponding in one dimension to be positive and negative for j < i and j > i, respectively.This force resembles the framework of the discrete element method where the overlap of particles substitutes their deformations at the contacts, which would be much more difficult and time consuming to resolve with a finite element model of deformable bodies.Assuming that the chain is precompressed by an external applied force F o , the characteristic overlap of the particles in static equilibrium ( o ), when all the contact stiffness ( κ (i,j ) ) of particles are chosen as κ o (mean characteristic contact stiffness), is thus defined as where the unit of κ depends on β.

Nondimensionalization
A length scale can be chosen such that the scaled particle overlap δ (i,j ) = δ (i,j ) / yields There are several length scales that can be chosen, e.g., the particle size, the driving amplitude or the initial overlap of the particles in static equilibrium.The latter is chosen for computations here so that (i,j ) = o ≡ 1 if all κ (i,j ) = κ o .Other dimensionless quantities are the mass b = m/ M 1 , where M 1 is the first moment of the mass distribution of the particles, as shown in Appendix C (the unscaled average mass of the particles); the dimensionless displacement u = u/ l; and the dimensionless spring constant κ = κ/ κ o .
The characteristic timescale becomes which gives us the dimensionless time t = t/ t c .The displacement of particle i from its equilibrium position o , so that the overlap becomes, δ (i,j ) = (i,j ) − (u (j )  − u (i) ).Finally, the interaction forces scale as 2.2 Equation of motion: nonlinear (Hertzian) The equation of motion for any particle i (except the boundary particles at either end of the chain) by using Eqs.( 3), (4) and nondimensionalization (Sect.2.1) can be written as which can also be written as .

Equation of motion: linear
The repulsive interaction force can be expressed as a power series and can be expanded about the initial overlap (i,j ) due to precompression.
For small displacements from the equilibrium condition (during wave propagation), using the definition of δ (i,j ) and after ignoring higher-order nonlinear terms, we arrive at − u (i) .(10) 438 R. K. Shrivastava and S. Luding: Effect of disorder on bulk sound wave speed Inserting the force relation (Eq.10) in Eq. ( 7), we get the general, linearized equation of motion: − u (i) .(11) For Hertzian nonlinear repulsive interaction force between the particles, the scaled stiffness κ (i,j ) and initial overlap (i,j ) are given as follows (see Appendix B for details): and The Hertzian nonlinear repulsive interaction force is appropriate for spherical particles (Landau and Lifshitz, 1970).Equation ( 11) can be written in a linearized form as Since we are interested only in mass disorder, we can choose all coupling stiffness (κ (i,j ) ) as 1.Then, Eq. ( 14) for individual particles can be written as The factor 1 1+β becomes 1 for the linear contact model (β = 0) and it becomes 2/3 for the Hertzian contact model (β = 1/2).It can be observed that the factor 1 1+β has only multiplicative influence on the physical parameters.Since in our system of equations (Eq.15) only mass disorder is present, the masses of the particles get multiplied by this factor ( 1 1+β ).For further analysis, β = 0 has been chosen so that This results in N equations which eventually can be expressed in matrix form: where M is a diagonal mass matrix with entries b (1) , b (2) , b (3) , . .., b (N) and zero otherwise; K is a matrix with diagonal entries −(κ (i+1,i) + κ (i−1,i) ) = −2, superdiagonal (κ (i+1,i) ) and subdiagonal (κ (i−1,i) ) entries +1 and zero otherwise for κ = 1.The term f is the external force which depends on the specified driving.Introducing A = −M −1 K, Eq. ( 17) can then be written as 2.4 Analysis in real space or spatial Fourier space Using an ansatz for real space and another ansatz for spatial Fourier space in Eq. ( 18) (the calligraphic fonts from now onwards will depict the spatial Fourier transform counterparts of the real space parameters), one has with wavenumber k and ku) dt du as double Fourier transform (spatial as well as temporal) ansatz.Equation ( 20) is a familiar eigenvalue problem.The eigenvalues ω 2 j and eigenvectors s (j ) of the matrix A give the eigendomain of the granular chain that are independent of the external driving.The square roots of the eigenvalues, ω j , are the natural frequencies of the chain.The set of eigenvectors can be orthonormalized to obey the orthonormality condition: with δ ij being the Kronecker delta symbol.The S matrix or the eigenbasis matrix can be constructed with s (j ) as the columns of the matrix, which can be used to transform back and forth from the real domain to the eigendomain.The columns (s (j ) ) of the matrix S are sorted such that the corresponding eigenvalues ω j are in increasing order.The vector of eigenmode amplitudes is A matrix G consisting of eigenvalues ω j along the diagonal (in increasing order) is formulated such that G = S −1 AS which allows the transformation of Eq. ( 17) into the eigendomain as which defines h and H implicitly.The differential equations (Eq.23) are decoupled and can be solved to give where C (1) is a diagonal matrix with C (1) j,j = sin(ω j t), C (2) is a diagonal matrix with C (2) j,j = cos(ω j t), and z P (t) or Z P (t) are the particular solutions of the differential equations, which depend on h or H and, hence, depend on the external driving force f or F. The vectors a or A and b or B are determined by the initial conditions from the initial displacement and The terms a and b or A and B are column vectors with column elements a j and b j or A j and B j , associated with a particular eigenfrequency (ω j ).The solution in real space can be obtained by the transformation mentioned in Eq. ( 22), which can be applied on Eq. ( 24) to give

Initial conditions: impulse driving
The initial conditions required to solve various special cases are the initial displacements (u o ) and initial velocities (v o ) in real space and V o and U o in spatial Fourier space.Besides the sinus driving used in Lawney and Luding (2014), we apply an impulse driving initial condition.For an impulse driving mode, the boundary conditions are as follows: An impulse driven chain has an impulse imparted to the first particle, i = 1 with initial velocity v o .Since the focus of our study is not on the occurrence of sonic vacuum (Nesterenko, 1983), the initial impulse (v o ) should be chosen small enough to avoid opening of contacts.Using Eqs. ( 25), ( 26) and ( 27) the initial conditions for the impulse driven chain, i.e., f = 0 (no driving present), and which implies that displacements and velocities of all particles p are given analytically by S pj S 1j sin(ω (j ) t) ω (j ) and In wavenumber space (spatial Fourier transform), the initial condition is specified by V o (k), which can be a sine or cosine function in terms of wavenumber (k).Using Eq. ( 27) and V o (k), we get and thus 2.6 Mass distribution, disorder parameter (ξ ), ensemble averaging & binning The mass distribution of the monodisperse chain has been selected randomly from normal (f and binary (f (bi) (b)) distributions whose standard deviations give the measure of the disorder of mass in the chain (ξ ).For instance, the normal distribution is given by High disorder means that the difference between the lightest particle and the heaviest particle is very large.It was observed in Lawney and Luding (2013) that the three distributions showed quantitatively similar behavior if the first two moments of the distributions were the same.Here, the first two moments of the aforementioned three distributions have been matched.The mathematical details of the distributions are given in Appendix C. Ensembles of chains with different realizations for a particular disorder and distribution have been taken into consideration.Angular brackets will be used to denote ensemble-averaged physical quantities such as u , E tot , etc.The first five moments of the three distributions for different disorder (standard deviation) ξ = 0, ξ = 0.1, ξ = 0.2, ξ = 0.35, ξ = 0.5 and ξ = 0.8 are given in Table 1 (500 ensembles scaled), Table 2 (500 ensembles unscaled) and Table 3 (10 000 ensembles).

Participation ratio & localization length
The participation ratio (P j ) (introduced in Bell and Dean, 1970, and used previously in Allen andKelner, 1998, Zeravcic et al., 2009) is a crucial tool in determining the localization length ( L j ) associated with a particular eigenmode.This localization length can be seen as the length beyond which elastic waves with a particular frequency become evanescent, i.e., they decay exponentially in a disordered system (Mouraille, 2009).It is instrumental in determining the length within which the elastic waves become confined in space and is dependent on the frequency and thus the eigenmode (Anderson, 1958).The participation ratio of eigenmode j is defined as Table 1.Scaled moments of ensemble-averaged distributions (500 ensembles) used for the one-dimensional chain (256 elements long).
Table 2. Unscaled moments of ensemble-averaged distributions (500 ensembles) used for the one-dimensional chain (256 elements long).
with the normalization condition on the eigenmodes For one dimension, the localization length is defined as L = P j d, where d is the particle center distance in equilibrium, i.e., under precompression.The localization length can now be nondimensionalized by the internal particle scale of separation ∼ d to give L j ∼ = P j .As discussed and pointed out in Allen and Kelner (1998), the localization length of the lowest eigenmode is often attributed to the length of the chain (which would be regarded as a force chain in our analysis), and hence it becomes important to find the localization length of an ordered chain, ξ = 0 as reference.

Distribution
Disorder The eigenvalues of this matrix are ω 2 j = 4sin 2 j π 2N and its eigenvectors are s (j ) = {sin j π N , sin 2j π N , sin 3j π N . . .sin (N−1)j π N }.After respecting the normalization condition and the definition of the participation factor, the localization length of the lowest eigenmode (P 1 ) can be analytically calculated from the eigenvectors as , and hence for N = 256, P 1 = 170.667≈ 171.

Dispersion
The analytical expression for the dispersion relation in an ordered chain of particles or elements with linear contact forces are given by (Brillouin, 1946;Tournat et al., 2004;Lawney and Luding, 2014) where the wavenumber can be nondimensionalized by the microscopic particle scale of separation ( d) and frequency by κ o M 1 giving the nondimensional dispersion relation with π = 2 for ordered chains with ξ = 0. Equation (39) holds for propagative as well as evanescent waves.The positive roots of this relation correspond to propagative waves and the imaginary roots to evanescent waves (Tournat et al., 2004).This expression also holds for longitudinal wave propagation in three-dimensional granular packings (Mouraille and Luding, 2008) and in one-dimensional chains (Lawney and Luding, 2014).From the dispersion relation, it can be noted that disorder creates a maximum permissible frequency ( π ) for propagating waves, frequencies below π are propagative and the frequencies above π are evanescent.The dispersion relation (Eq.39) for ordered chains which is the dispersion relation for propagative waves.
www.nonlin-processes-geophys.net/24/435/2017/ Nonlin.Processes Geophys., 24, 435-454, 2017 442 R. K. Shrivastava and S. Luding: Effect of disorder on bulk sound wave speed 2.9 Total energy dispersion From Eq. (A5) it can be observed that the total energy of the eigenmodes is constant with respect to time as given by By taking the first moment of this eigenmodal total energy representation about frequency, a dominant frequency related to a particular wavenumber can be obtained.Moments of the eigenmodal total energy representation are defined as The dominant frequency is given by the first moment, The dominant frequency can be measured by averaging over all eigenmodes for a single realization with A j (k) as a multiplicative factor which depends on the Fourier initial condition V o (k) (Eq.32).The dispersion relation for the propagating waves can be obtained by taking ensemble averages of this dominant frequency ( (k) ), which will be plotted in Fig. 10b below for different disorder strengths (500 ensembles).

Group velocity
The group velocity is given by for both propagative waves and evanescent waves.It can be obtained by differentiating Eq. ( 40) that where π = π (ξ ) depends on disorder, as we will see below.

Results and discussions
The analytical expressions derived in the previous sections are computed for chains that are N = 256 particles long.The impulse imparted to the first particle is v o = 0.05.The time step utilized for the output is t = 0.0312 and the maximum time up to which the computations have been carried out is t max = 256 such that the pulse has just about reached the 256th particle.As can be seen from Tables 1 and 3, the scaled average mass of the particles has been kept M 1 = 1 and ξ = 0.0, 0.1, 0.2, 0.35, 0.5 and 0.8 disorder parameters (standard deviation; see Appendix C) have been used for analysis.Using the analytical solution of the linearized system (Eq.31), ensembles of 500 and 10 000 chains along with representative single realizations will be shown in this section.

Nonlinear (Hertz) and linear space-time responses
Equation ( 8) with Eqs. ( 12) and ( 13) has been solved numerically with Verlet integration to get space-time responses of particles with nonlinear (Hertzian) repulsive interaction forces.The time step used for the numerical integration is t = 0.00038147.Figure 3 shows the space-time responses calculated numerically for the nonlinear equation of motion (Eq.8) and the space-time responses calculated for the linearized equation of motion (Eq.14) using a small initial velocity v o = 0.05.The space-time responses are obtained for a single realization of a granular chain without ensemble averaging.The nonlinear space-time responses coincide with the linear space-time responses for small enough v o , confirming that the solution given by Eq. ( 31) is also appropriate for particles with nonlinear repulsive interaction forces for small displacements.
In order to quantify the limitations of the linear spacetime responses obtained from Eq. ( 31), Fig. 4 is plotted.The difference between the maximum value (u peak ) of the spacetime responses for Hertzian and linear repulsive interaction forces (u peak(linear) ) is chosen as a parameter to judge the appropriateness of the linear space-time response for the nonlinear equation of motion (Eq.8).The difference increases nonlinearly irrespective of particle position and disorder parameter of the granular chain, with very good agreement for v o < = 0.05.

Displacement response of three mass distributions
Only the mass disorder of the particles in the chain with length 256 is taken into consideration and κ is chosen as 1 (Sect.2.3). Figure 5 shows the displacement as a function of time of the 150th particle (Fig. 1a and c) and of the 220th particle (Fig. 1b and d), which are placed before and after the reference localization length (the maximum possible, L max = 171, Sect.3.7) for two disorder parameters ξ = 0.1 and ξ = 0.5 with three mass distributions (normal, uniform and binary).For weak disorder (ξ = 0.1), it is observed that the displacement wave packets are perfectly superposed, affirming the conclusion in Lawney and Luding (2013) and Lawney and Luding (2014) that the shape of the distribution has no effect on the propagating pulse if the first two moments of the distribution are the same (Table 1).For stronger disorder (ξ = 0.5), the wave packets are not collapsing per- fectly (Fig. 1c and d).As can be seen in Table 1, there is a numerical mismatch between the unscaled moments of the distributions, leading to a dissimilarity between the second scaled moments ( M 2 ).This also causes the real standard deviation (disorder; ξ ) which has been numerically calculated ( ; Table 1) to deviate a little bit from its intended value.It can also be observed from Fig. 5c and d that the pulse shapes of binary distribution and uniform distribution are closer to each other in comparison to normal and binary or normal and uniform, since the scaled second moments ( M 2 ) of binary and uniform distributions for ξ = 0.5 are closest to each other (Table 1).Similar conclusions about similarity, dissimilarity and closeness can be drawn about pulse shapes of different distributions for different disorder parameters (ξ (intended), (numerically obtained)) on the basis of moments of the mass distribution.For larger ξ , higher moments (as listed in Tables 1, 2 and 3) have to be considered (Ogarko and Luding, 2013), but discussing these higher moments and their consequences go beyond the scope of this study.

Displacement response for different disorder parameters (ξ )
Mechanical waves propagating through disordered media or granular media such as soil (on the receiver end) can be divided into two parts, the coherent part and the incoherent part (Jia et al., 1999;Jia, 2004).The coherent part is the leading wave packet and is self-averaging in nature (it maintains its shape after ensemble averaging) and it is used for determining the bulk sound wave velocity.In contrast, the incoherent part is the scattering, non-self-averaging part, which is  strongly system-configuration-dependent, also known as the coda or tail of the mechanical wave.Figure 6 contains the displacement of two particles (150th and 220th particles) used in Sect.3.2, for consistency.Here, attention has been given to the effect of the mass distribution on the time of arrival or flight, and hence the wave velocity of the initial wave packet.Figure 6a and b contain the displacements of the 150th and 220th particle and after L , Sect.3.7) for single realizations, 500 and 10 000 ensembles.The leading wave packet is the same for 500 and 10 000 ensembles in both figures, i.e., the coherent part of the wave which maintains its shape after averaging.The coda is more or less pronounced at 150 or 220, respectively, and vanishes due to ensemble averaging.Figure 6c and d show the displacement response of the 150th and the 220th particle with respect to time for chains with different mass disorder.The speed of the coherent wave packet (from source to receiving particle, peak of signal) is increasing with disorder.Higher disorder leads to higher coherent wave (peak) speed, irrespective of the localization length (L max ).However, this increase in wave speed can also be attributed to sound wave acceleration near the source, as pointed out by Mouraille et al. (2006), and may not be generalized as effect of mass disorder in the chain, as investigated in the next section.

Coherent wave speed and disorder
Tables 4, 5 and Fig. 7 contain the velocity of the peak of the coherent wave, the velocity of the rising part of the coherent wave packet when the displacement of the particle has attained 5, 10, 70 and 90 % of the peak value, and the first time when the displacement of the particle becomes 0 after it has attained the peak value of the coherent wave (zero crossing), all constituting the coherent wave packet.The velocities were determined through velocity picking (particle position divided by the time of arrival).The particles used for computing the velocities were 130, 150, 200 and 220 (Table 4; 2 before localization length and 2 after localization length, L max = 171).It can be observed that irrespective of the rising part of the coherent wave packet and the peak (Fig. 7a,  b, c, d), the wave velocity increases with disorder.However, for zero crossing (Fig. 7e), the velocity decreases with an increase in disorder, and the same can be said for the part Nonlin.Processes Geophys., 24, 435-454, 2017 www.nonlin-processes-geophys.net/24/435/2017/ of the coherent wave packet which lies after the peak value; this can be attributed to the increased spreading of the wave packet with an increase in disorder.Notably, the speed measured at particle 130 is larger (smaller) if the earlier (later) parts of the signal are considered.Figure 7f shows the velocity of the peak value of the coherent wave packet of all the particles of the granular chain for different disorder parameters and it also exhibits a similar kind of acceleration of signal or mechanical wave near the source to what was observed in Mouraille et al. (2006).This acceleration is caused by self-demodulation of the initial impulse imparted to the granular chain and the noteworthy point is that the velocity increases while the acceleration decreases with an increase in disorder.Due to this observation we cannot generalize the effect of disorder on wave speed.The sudden rise in velocity of the peak value in Fig. 7f after the 250th particle is due to the boundary effect as well as due to the presence of the coherent wave front of the traveling wave around that position as the maximum time window used is t max = 256.For practical purposes, we remark the wave speed measured varies by a few percent up or down, dependent on which part of the signal is used for measurement.
To understand the effect of disorder on wave speed without taking into account this "source effect", the velocity based on the time taken by the pulse for propagating a common distance of seven particles has been computed in Table 6 and the results have been plotted in Fig. 8.The reason for selecting such a low common distance of separation was to keep the effect of the source as minimal as possible; the sets of points (130 to 137, 150 to 157, 220 to 227 and 240 to 247 particles) were used with different reference points of the coherent wave front (5, 10, 70, 90 %, peak value and zero crossing).From Fig. 8 and Table 6 it can be observed that the same trend is followed, except the velocity computed using the zero crossing reference point, which is more or less constant with little fluctuations (Table 6).Figure 8a shows a consistent increase of velocity as it is the closest to the source (dominated by the source effect); however, as the set of particles is selected farther away from the source, the velocity decreases slightly and then increases with increasing ξ (Fig. 8b and c). Figure 8d exhibits a consistent decrease of velocity with an increase in ξ because the set of particles (247-240) are far from the source (source effect is weakest).From Fig. 8d, it can be concluded that higher disorder results in a decrease in wave velocity.

Frequency response and dispersion
In Fig. 9a and c a fast Fourier transform (FFT) with respect to time is carried out on the displacement response of a 256element-long chain for disorder, with ξ = 0.01 and ξ = 0.35, respectively (when an impulse of v o = 0.05 has been applied to the first particle) to observe the frequency content with distance, (the sampling frequency is ω sample = 2π t ) and responses up to half of the sampling frequency were taken into account to avoid aliasing (Nyquist criterion; Shannon, 1949).The first five particles have been excluded from the Fourier transform to avoid an overwhelming driving signal effect.Figure 9a exhibits the existence of a cut-off frequency (ω = 2) above which the waves become evanescent.The bending of the intensity with distance (particle number), especially at large distances, is attributed to dispersion and the finite time window.Using the group velocity (p is the particle number) for an ordered chain (ξ = 0), π = 2 and using Eq. ( 45), the frequency envelope is which is the red curve plotted in Fig. 9a.A spatial as well as temporal two-dimensional FFT is carried out for a single realization of a 256-element-long chain with disorder ξ = 0.01 and ξ = 0.35 to observe the dispersion relation (Fig. 9b and d; ω vs. k).Two-dimensional FFT has been used previously for one-dimensional and threedimensional polydisperse granular packings for obtaining dispersion relations (Luding and Mouraille, 2008;Lawney and Luding, 2014;O'Donovan et al., 2015), but strong frequency filtering due to the disordered system resulted in ambiguous dispersion relations (flat bands and absence of certain frequencies below the cut-off frequency, which indicates the nonpropagative bands due to the presence of defect modes).This can also be observed from Fig. 9b and d.Equation (40) (the dispersion relation for an ordered chain) has been plotted in Fig. 9b, which gives a perfect fit for the denser regime in the figure.However, for the disordered chain, ξ = 0.35, as proposed earlier in Sect.2.9, the dispersion relation is better given by (k) by ensemble averaging the dominant frequencies with respect to different wavenumbers.The term (k) for 500 ensembles with disorder ξ = 0.35 has been plotted in Fig. 9d (the green curve).For low frequencies the green curve perfectly superposes the dense regime in the displacement's temporal and spatial Fourier transform; for higher frequency (ω > 1.5) due to the appearance of a flat band (defect mode) the intensity is not visible near the green curve, which holds true for low and intermediate frequencies and wavenumbers.

Total energy dispersion in disordered chains
(k) from Eq. ( 43), which was plotted for ξ = 0.35 in Fig. 9b, is plotted for ξ = 0.1, 0.2, 0.35, 0.5 and 0.8 in Fig. 10a.It is observed that the maximum permissible frequency ( π ) above which the waves become evanescent decreases with increasing disorder.The slope of ω vs. k curves indicates the wave speed which clearly can be observed to be decreasing with increasing disorder, confirming what was observed in Sect.3.4.

Participation ratio & localization length
Figure 11 shows the participation ratio ( P ), i.e., the localization length ( L , from Sect.2.7) for binned 500 ensembleaveraged realizations of chains (with 0.0781 as frequency bin size) and with different disorder parameters ξ .The lowest frequencies have the same localization length L max = 171, independent of the disorder of the chain; see Sect.2.7.Towards higher frequency the localization length decays to zero more rapidly with increasing disorder, characterized by a particular frequency (crossover or pass frequency) where p = L max /2.Unlike infinitely long chains, where L ∝ ω −2 (as suggested in Azbel, 1983), the finite disordered chains for higher frequencies have L ∝ ω −q where q 2, decreasing with increasing disorder.For understanding the effect of different types of disorder on the pass frequencies (ω pass ) associated with the localization length (L), ω pass(1/2) (the frequency associated with L max /2), has been plotted in Fig. 12a, and L, associated with ω = 1 scaled with L max , has been plotted in Fig. 12b, selected from the dashed-line crossings in Fig. 11a.Both quantities exhibit a decreasing trend with an increase in disorder, as characterized by the empirical fits: where ω o = 0.3502 and ξ ω = 0.2536, and where L o = 2.7668 and ξ L = 0.3165, with some errors of the order of ±5 %.Both curves saturate for large ξ values.More data and a closer analysis are necessary for improving this analysis and putting a better theoretical basis to the fits.

Total energy of eigenmodes
The density of states or density of vibrational modes is an important quantifying factor in studying the vibrational properties of materials such as jammed granular media (Schreck et al., 2014), etc.However, it tells us only about the number of vibrational modes but does not paint the complete picture of spectral properties of momentum and energy transport.Equation (A5) gives us the energies of individual eigenmodes and shows that the energy is constant with respect to time. Figure 13a plots the ensemble-averaged density of states for 500 mass-disordered granular chains with frequency bins of size 0.0781.The peak of the density is decreasing with increasing disorder and shifting to smaller ω. Figure 13b gives the ensemble-averaged energy spectrum for the same frequency bins used in Fig. 13a (500 realizations), giving an energy distribution over frequency.The shape of the energy distribution is wider for larger ξ ; the energy distribution becomes more sharp, shifting to smaller ω with increasing disorder.In both plots, Fig. 13a and b, the tails are broader for larger ξ , where the shapes in panel (a) are independent of driving, while the shapes in panel (b) depend on the initial condition.

Conclusions
An impulse driven wave propagating through a precompressed mass-disordered granular chain has been studied.Motivation comes from the existence of force chains which form the backbone network for mechanical wave propagation in granular materials such as soil.The scaled standard deviation of the mass probability distribution of the elements or particles of the granular chain has been identified as the relevant disorder parameter (ξ ; see Sect.2.6), as suggested already in Lawney and Luding (2014).Chains with normal, binary and uniform mass distributions have quantitatively identical signal transmission characteristics as long as the first two moments of the mass distribution are the same and ξ is not too large.Interestingly, on first sight, the dependence of wave speed on magnitude of disorder looks nonmonotonous.This surprising increase of wave speed for weak disorder, and decrease for stronger disorder, is due to two different effects overlapping: the increase of wave speed takes place close to the source (see Fig. 7), i.e., our one-dimensional granular chain has the ability to model the physics of accelerating waves, as observed in complex higher dimensional granular structures (Mouraille et al., 2006).The competing mechanism of decreasing wave speed with disorder is only clearly observed when the velocities are measured via the travel time while maintaining constant separation far away from the source (Fig. 8 and Table 6).The group velocity given by Eq. ( 45) also shows a decrease in wave speed with an increase in disorder.When the travel time is measured from the source, the two mechanisms overlap and interfere, causing the nonmonotonous behavior, but possibly allowing for the tuning of particular propagation characteristics in short chains.
As another main result, Eq. ( 43) gives an effective, weighted dispersion relation as the normalized first moment of eigenmodal (total) energies with frequency.This gives a much better signal to noise ratio for ω vs. k in comparison to two-dimensional FFT of displacement or velocity signals, reported previously (Mouraille et al., 2006;O'Donovan et al., 2015).The upper (maximum permissible) limit frequency due to the discreteness of the system slightly decreases with increasing disorder, ξ , and waves consistently propagate slightly slower with increasing disorder if scaled by mean mass, i.e., an effect of ξ .From the energy content one also observes (in disordered systems) that waves above a low-frequency pass band (ω pass(1/2) ) become evanescent after they have traversed a localization length, L = L(ξ, ω), associated with a particular pass frequency (ω = 1) for which (yet) no analytical prediction is known to the authors (Otsubo et al., 2017).
The energy analysis presented in this article can be used for understanding pulse propagation in disordered, weakly or strongly nonlinear granular chains and its attenuation, widening and acceleration (experimentally and numerically investigated in Langlois and Jia, 2015).It would also be interesting to understand the effect of damping on the eigenmodes, velocity of the propagating wave, changes in frequency filtering and the energy of the eigenmodes.Also, a different kind of averaging (micro-macro transition) should be developed using frequency bands to develop a master equation for propagation (or localization) of total energy in terms of wavenumber and frequency at different regimes of disorder, nonlinearity and material properties.Such macromodels, taking into account multiple scattering, dispersion, attenuation, etc., will allow for modeling of realistic wave propagation in granular materials such as soil on large scales.
Data availability.Data have been generated using the aforementioned theoretical model.The readers can reproduce it by using the equations mentioned in their respective sections.
www.nonlin-processes-geophys.net/24/435/2017/ Nonlin.Processes Geophys., 24, 435-454, 2017 summation).Hence, E tot (ω j ) = 1 2 a 2 j ω 2 j .Now, by replacing u, v, a with their spatial Fourier transform counter parts U, V and A (calligraphic) by using the ansatz in spatial Fourier space as in Eq. ( 19) for Eq. ( 18), we obtain the harmonic total energy (in terms of wavenumber): Appendix B: Hertz contact model If a Hertzian repulsive force is taken into consideration between particles (Landau and Lifshitz, 1970;Lawney and Luding, 2014), then κ (i,j ) = Y (i,j r i r j r i + r j 1/2 , (B1) where E i and ν i are the elastic modulus and Poisson's ratio, respectively, of particle i.If the particles are made up of the same material, Y (i,j ) and ν become same for all the contacts, The characteristic stiffness of the contact is The characteristic initial overlap becomes The characteristic time is Nonlin.Processes Geophys., 24, 435-454, 2017 www.nonlin-processes-geophys.net/24/435/2017/

Figure 1 .
Figure 1.A granular or force chain from a network (schematic).

Figure 2 .
Figure 2. Prestressed chain of granular elements during dynamic wave propagation.

Figure 3 .
Figure 3.The displacement as a function of time is shown for the 100th particle in a chain of particles with disorder parameter ξ = 0.1 in (a) and the 130th particle in a chain of particles with disorder parameter ξ = 0.35 in (b) for impulse v o = 0.05.

Figure 4 .Figure 5 .
Figure 4.The nonlinear increase in u (p) diff (a parameter which shows dissimilarity between linear and nonlinear space-time responses) with initial impulse velocity (v o ).The value v o = 0.05 is in the zone where linear and nonlinear space-time responses are almost identical.

Figure 6 .
Figure 6.Displacements of 150th (a, c) and 220th (b, d) particle with respect to time for different ensembles (a, b) of disorder ξ = 0.5 and of multiple disorder parameters (c, d); for 500 ensembles.

Figure 7 .
Figure 7. Coherent wave velocities determined through velocity picking.The peak of the coherent wave packet's velocity (d) and the rising part of the packet (a, b, c) as well as the falling part (e) are taken into consideration.Panel (f) displays the peak velocity for all the particles in a granular chain with different ξ .Note that (a, b, e) have a different axis range as (c, d).

Figure 8 .
Figure 8. Wave speed for common distance of separation (seven particles or elements) with different disorder ξ .Note that the vertical range in (a, c) is different from (b, d).Red line is 70 % of the peak value and blue line gives the velocities when the peak is considered.

Figure 9 .
Figure 9. Panel (a) is the temporal Fourier transform of displacement of particles for normal distribution and disorder parameter ξ = 0.01 (single realization) with group velocity (v g ) depicting the propagation of the wave front, and (b) is the temporal as well as spatial Fourier transform (two-dimensional FFT, single realization) calculated for obtaining the dispersion relation of a chain, while (k) gives the true ensemble-averaged (500) dispersion relation from Eq. (43).Panels (c, d) are the higher-disorder ξ = 0.35 counterparts of (a, b), respectively.

Figure 13
Figure 13.(a) Density of states and (b) energies of the binned frequencies for different disorder ξ for all bin sizes in ω.

Figure A1 .
Figure A1.Multiplicative Factor a j for normal distribution obtained after ensemble averaging (500).

Table 4 .
Scaled coherent wave velocity picking for different particles before and after localization length for a disordered chain with normal distribution (256 elements long, 500 ensembles).

Table 5 .
Unscaled coherent wave velocity picking ( M 1 ) for different particles before and after localization length for a disordered chain with normal distribution (256 elements long, 500 ensembles).

Table 6 .
Coherent wave velocity calculated from the time taken by the pulse to travel a common distance of separation (seven particles or elements) with time calculated in reference to 5, 10, 70 and 90 % of the peak value and the peak value of the coherent wave packet.