Superdiffusive Transport and Energy Localization in Disordered Granular Crystals
Abstract
We study the spreading of initially localized excitations in onedimensional disordered granular crystals. We thereby investigate localization phenomena in strongly nonlinear systems, which we demonstrate to be fundamentally different from localization in linear and weakly nonlinear systems. We conduct a thorough comparison of wave dynamics in chains with three different types of disorder: an uncorrelated (Andersonlike) disorder and two types of correlated disorders (which are produced by random dimer arrangements), and for two families of initial conditions: displacement perturbations and velocity perturbations. We find for strongly precompressed (i.e., weakly nonlinear) chains that the dynamics strongly depends on the initial condition. In particular, for displacement perturbations, the longtime asymptotic behavior of the second moment has oscillations that depend on the type of disorder, with a complex trend that is markedly different from a power law and which is particularly evident for an Andersonlike disorder. By contrast, for velocity perturbations, we find that a standard scaling (for some constant ) applies for all three types of disorder. For weakly precompressed (i.e., strongly nonlinear) chains, and the inverse participation ratio satisfy scaling relations and , and the dynamics is superdiffusive for all of the cases that we consider. Additionally, when precompression is strong, the inverse participation ratio decreases slowly (with ) for all three types of disorder, and the dynamics leads to a partial localization around the core and the leading edge of the wave. For an Andersonlike disorder, displacement perturbations lead to localization of energy primarily in the core, and velocity perturbations cause the energy to be divided between the core and the leading edge. This localization phenomenon does not occur in the sonicvacuum regime, which yields the surprising result that the energy is no longer contained in strongly nonlinear waves but instead is spread across many sites. In this regime, the exponents are very similar (roughly and ) for all three types of disorder and for both types of initial conditions.
pacs:
63.50.x, 45.70.n, 46.40.Cd, 87.15.hjI Introduction
The study of wave propagation in disordered lattice and continuum systems has been an important and popular research theme during the past several decades. Some of the most prominent recent studies on these topics Schwartz:Nature2007 (); Billy:Nature2008 (); Roati:Nature2008 (); Gurarle:PRL2008 (); Conti:NatPhys2008 (); Efremidis:PRL2008 (); Hu2008 (); Gremaud:PRL2010 () have generalized to weakly nonlinear settings the ideas of P. W. Anderson, who showed theoretically that the diffusion of waves is curtailed in linear random media (where the randomness arises from defects or impurities) Anderson:PR1958 (); Kramer:RPP1993 (). This interplay between disorder and nonlinearity — which often arises in the presence of lattice discreteness — is of considerable interest to a vast array of ongoing studies, as is evidenced by the recent reviews Flach:ArxivRep20141 (); Flach:ArxivRep20142 () (see also the numerous references therein). The set of different physical scenarios in which Anderson localization has been investigated is staggering: it ranges all the way from electromagnetism Schwartz:Nature2007 () and acoustics Hu2008 () to subjects such as quantum chromodynamics andersonqcd ().
As in the above studies, we are interested in waves in disordered media, but we depart from the earlier work in a very important way: we seek to explore order–disorder transitions with a particular emphasis on strongly nonlinear media. This contrasts sharply with the linear and weakly nonlinear media in which Andersonlike models have traditionally been studied Flach:ArxivRep20141 (); Flach:ArxivRep20142 (). Our approach is motivated predominantly by the strong (and increasing) interest in granular crystals Nesterenko:book (); Sen:PR2008 (); pgk:2011 (), which (as we discuss below) are very important both for the study of fundamental nonlinear phenomena and for numerous engineering applications. The examination of disordered systems in general — and of Andersonlike phenomena in particular — is a key challenge in the study of nonlinear chains Kuzovkov:PhysScr2011 (); Skokos:PRE2009 (); Mulansky:PRE2009 ().
Onedimensional (1D) granular crystals, which consist of closely packed chains of elastically colliding particles, are a paradigmatic system for the investigation of chains of strongly nonlinear oscillators. Their strongly nonlinear dynamic response has inspired numerous studies of the interplay between nonlinearity and discreteness Nesterenko:book (); Sen:PR2008 (); pgk:2011 (). One can construct granular crystals using materials of numerous types and sizes, and their properties are thus extremely tunable, tractable, and flexible Nesterenko:book (); Sen:PR2008 (); Daraio:PRE2006 (); Coste:PRE1997 (). This also makes them very wellsuited for investigating the effects of structural and material heterogeneities on nonlinear wave dynamics. Studies have examined the role of defects Sen:PRE1998 (); PhysRevE.57.2386 (); Hascoet:EPJB2000 (); Hong:APL2002 (); Theocharis:PRE2009 () (including in experimental settings Job:PRE2009 (); Many ()), interfaces between two different types of particles Nesterenko:PRL2005 (); Daraio:PRL2006 (), decorated and/or tapered chains Doney:PRL2006 (); Harbola:PRE2009 (), chains of diatomic and triatomic units Porter:PRE2008 (); Porter:PhysicaD2009 (); Herbold:Mechanica2009 (); Molinari:PRE2009 (); staros1 (); staros2 (); staros3 (); boechler4 (), and quasiperiodic and random configurations Sokolow:AnnPhys2007 (); Chen:PhysicaB2007 (); Fraternali:MAMS2010 (). The tunability of granular crystals is valuable not only for fundamental studies of their underlying physics but also in potential engineering applications — including shock and energy absorbing layers Daraio:PRL2006 (); Hong:PRL2005 (); Fraternali:MAMS2010 (), sound focusing devices and delay lines Spadoni:PNAS2010 (), actuators Khatri:SPIE2008 (), vibration absorption layers Herbold:Mechanica2009 (), sound scramblers Daraio:PRE2005 (); Nesterenko:PRL2005 (), and acoustic switches and logic gates daraio_natcom (). Because one can model granular chains as a type of FermiPastaUlam (FPU) lattice, they have also been employed in studies of phenomena such as equipartition (see, e.g., szel2013 (); zhang2015 ()).
As was illustrated recently, localization in strongly nonlinear systems can have a fundamentally different character from localization in linear and weakly nonlinear systems Ponson:PRE2010 (). Importantly, one can use the setting of granular crystals to explore a regime (the socalled “sonic vacuum”) in which no linear counterpart whatsoever exists Nesterenko:book (). It is our goal in the present paper to investigate this regime and nearby regimes in detail and to conduct what we believe is the first systematic study of the differences between localization is linear, weakly nonlinear, and strongly nonlinear systems. There are numerous types of disorder in a granular chain, and — as we demonstrate in this paper — it matters whether the disorder is uncorrelated (as in the original Anderson model) or correlated. Moreover, there are multiple types of possible correlations in disordered arrangements, and we illustrate using randomized arrangements of dimers (see Ref. Ponson:PRE2010 () for an example arrangement that was studied in the context of granular crystals) that seemingly small differences in disorder can have a large impact on the dynamics of wave propagation in strongly nonlinear systems. Moreover, because granular chains are a type of FPU system FPU () — so nonlinearities arise from the potentials that connect adjacent nodes of the lattice — they are fundamentally different from the nonlinear Schrödinger (NLS) and Klein–Gordon (KG) lattices in which disordered configurations have been extensively studied recently Flach:ArxivRep20141 (); Flach:ArxivRep20142 (). In fact, as we will demonstrate in the present paper, this difference leads to much more rapid transport in disordered granular chains than what occurs in either NLS or KG lattices. This fundamental difference is one of the main findings of our work: the dynamics of strongly nonlinear, disordered granular crystals includes regimes with superdiffusive transport.
Much of the significant volume of work involving nonlinear disordered lattices has focused on the effect of weak nonlinearity on the wellestablished Anderson model Flach:ArxivRep20141 (); Flach:ArxivRep20142 (); Pikovsky:PRL2008 (); Tietsche:EPL2008 (); Johansson:EPL2010 (); Laptyeva:EPL2010 (); GarciaMata:PRE2009 (). One of the most remarkable findings in this body of work is the fact that a small amount of nonlinearity in a disordered lattice can induce interaction between Anderson modes, which eventually can lead to a subdiffusive delocalization process. Interestingly, this phenomenon emerges in a highly nontrivial way: even when both disorder and nonlinearity separately tend to localize energy, they also “interfere” with each other’s transportgeneration processes and consequently destroy the Andersonlocalization mechanism. In the context of strongly nonlinear disordered lattices, a noteworthy recent effort is that of Mulansky:NJP2013 (). The authors of that paper examined disordered lattices — where disorder is introduced via either a linear or a nonlinear onsite term — in which the coupling leads to a strongly nonlinear setting. They found that initially localized wave packets tend to spread in a subdiffusive way. However, the spreading in nonlinearly coupled linear oscillators is slow in comparison to purely subdiffusive behavior. We believe that the subdiffusive behavior in their setting is a consequence of the local potential, because (as we demonstrate in our paper) the dynamics tends to be superdiffusive when only strongly nonlinear interactions are present.^{1}^{1}1We note in passing that superdiffusive behavior was reported very recently in a socalled “pseudotwodimensional” random dimer rodrigo (). Both the setting and qualitative behavior of the system in Ref. rodrigo () are different from ours in fundamental ways. They consider purely linear dynamics, whereas we consider both linear and (especially) nonlinear dynamics. Additionally, they found superdiffusive transport for lattices with shortrange correlations but subdiffusive transport for uncorrelated disorders (such as the one in the Anderson model), whereas we find superdiffusive transport in lattices with either correlated or uncorrelated disorder.
The remainder of our paper is organized as follows. In Sec. II, we describe the fundamental equations that characterize a disordered granular chain with Hertzian interactions, and we examine different approximations that depend on the amount of precompression. In Sec. III, we present three different types of disorder and study their correlation properties, which will prove to be of crucial importance for the qualitatively different transport dynamics that can occur. In Sec. IV, we briefly discuss the influence of impurities in homogeneous chains on the modes that emerge and on the dynamics more generally. In Sec. V, we present our computational results. We describe the fundamental differences between the different types of disorder, and we discuss the properties of the linear spectrum and the different types of impuritylike modes that appear for each type of disorder. We also study the transport and localization properties for both linear and nonlinear waves for each type of disorder, and we demonstrate with numerical simulations that localization (as either breathers or traveling waves) is no longer possible in a strongly nonlinear regime. Additionally, for strongly precompressed chains and for initially localized displacement excitations, we demonstrate that the second moment exhibits a complicated trend that differs markedly from a power law. By contrast, we observe superdiffusive transport for all of the other configurations and initial conditions. We summarize our conclusions and discuss future challenges in Sec. VI.
Ii Disordered Granular Crystals
ii.1 Equations of Motion
One can describe a 1D crystal of spherical particles as a chain of nonlinear coupled oscillators with Hertzian interactions between each pair of particles Nesterenko:book (); Sen:PR2008 (); pgk:2011 (). Hertzian forces are applicable to a wide variety of materials Johnson:Book () (including steel, aluminum, brass, bronze, and many more). The equations of motion in this setting are
(1) 
where is the displacement of the th particle (where ) measured from its equilibrium position in the initially compressed chain, is the mass of the th particle, and
(2) 
is a static displacement for each particle that arises from the static load . The parameter is given by
(3) 
where the elastic modulus of the th particle is , the Poisson ratio of the th particle is , and the radius of the th particle is . A Hertzian interaction between a pair of particles occurs only when they are in contact, so each particle is affected directly only by its nearest neighbors and experiences a force from a neighbor only when it overlaps with it. This yields the bracket
(4) 
in Eq. (1). The exponent and prefactor in Eq. (1) are consequences of the elastic nature of the particle interactions and of the particle geometry Nesterenko:book (); Sen:PR2008 (). Other particle shapes, such as ellipsoids Ngo:PRE2011 () and cylinders Khatri:Granular2012 (), can also exhibit Hertzian interactions.
The boundary conditions of Eq. (1) are given by considering and . If one of the radii in Eq. (3) is infinite, then one obtains the interaction coefficient between an elastic plate and an elastic sphere:
(5) 
where is the elastic modulus, is the Poisson ratio of the elastic plates at the boundaries, and the suffixes and , respectively, indicate the left and right boundaries of the chain. Consequently, the equations of motion for spheres and are
(6)  
(7) 
Equation (1) does not include effects from restitution or dissipation, so we assume that we can neglect energy that dissipates into internal degrees of freedom. Most investigations of granular crystals make these assumptions, and a conservative (and Hamiltonian) description of granular crystals has been extremely useful for numerous comparisons of theoretical and computational results to laboratory experiments Nesterenko:book (); Sen:PR2008 (), including in the presence of disorder Ponson:PRE2010 (); Fraternali:MAMS2010 (). The proper physical form of dissipation is not known and is still a subject of ongoing debate. See Refs. Sen:PR2008 (); nester2007 (); Carretero:PRL2009 (); vergara () for recent discussions of dissipative forces in granular crystals. Note additionally that we will not worry about incorporating proper restitution forces, as we conduct our simulations in domains of sizes that ensure that the waves that we examine do not reach the domain boundary during the reported time horizon.
Even for homogeneous chains, Eq. (1) includes several interesting features that are not present in other lattice models (such as the wellknown nonlinear Schrödinger (NLS) and KleinGordon (KG) lattices Flach:PR2008 ()). From a structural perspective, the present model is a type of FPU lattice FPU (). It exhibits important differences from NLS and KG lattices, which typically include both a linear coupling and an onsite nonlinearity. However, there are also respects in which granular chains differ fundamentally from traditional FPU models Sen:PR2008 (). In particular, when there is no precompression (which is known as the “sonicvacuum regime” Nesterenko:book ()), the sound speed goes to and the system becomes purely nonlinear (i.e., linearizing it simply yields ). This allows solutions like compactons to occur in PDE limits of Eq. (1). Additionally, because compactons are not exact solutions of the original granular chain (which has a fundamentally discrete nature), traveling waves in strongly nonlinear regimes exhibit a superexponential decay at the edge of the distribution instead of having compact support chaterjee (); Pikovsky:PhysicaD2006 (). In the presence of precompression, which yields a linear term in Eq. (1), Refs. cho1 (); cho2 () illustrated both analytically and numerically that energylocalizing states can arise in the form of dark breathers.
ii.2 Precompression Regimes
By changing the magnitude of the static load relative to the displacements between particles, one can tune the strength of the nonlinearity in Eq. (1). To do this, we approximate the force using a powerseries expansion, which is known to be suitable for a strongly compressed or weakly nonlinear chain Nesterenko:book (). We thereby distinguish three different regimes, which we now discuss.
ii.2.1 “Linear” Regime ()
In this regime, we linearize Eq. (1) about the equilibrium state in the presence of precompression to obtain
(8) 
where
(9) 
Note that we have neglected the higherorder terms (even the quadratic ones) in the expansion for very weak strains (i.e., small relative displacements). This linear limit corresponds to a chain of coupled harmonic oscillators. We represent the solutions to Eq. (8) as complex wavefunctions to obtain a complete set of eigenfunctions of the form , where is the eigenfrequency (so the eigenvalue is ). In the homogeneous case — i.e., for for all — we obtain plane waves , and the dispersion relation,
(10) 
where is the mass of a particle, gives a single acoustic branch. The frequency satisfies the bounds and , so the group velocity in this homogeneous case is
(11) 
The maximum of the group velocity is , and it occurs when . From , we are able to write expressions for several quantities. For instance, given an initially localized excitation , wave spreading takes place within a cone . Therefore, in our simulations, we consider systems that have at least spheres, where is the integration time and we recall that the ceiling function is . This consideration allows us to avoid boundary effects when we study dynamics.
For an arbitrary arrangement of spheres in a granular chain, the eigenvalue problem associated with Eq. (8) takes the generic form
(12) 
where is the eigenvalue and , where are the elements of the diagonal matrix of masses and
(13) 
are the elements of a tridiagonal symmetric matrix. Note that we have used fixed boundaries at both ends of the chain (see Eqs. (6) and (7)). For a disordered chain (see Sec. III for the different types of disorder that we study), both matrices have random entries. Consequently, is an asymmetric tridiagonal matrix with random entries.
ii.2.2 “Weakly Nonlinear” Regime ()
This amounts to a particular case of the FPU model FPU () that includes the socalled “” and “” terms from two of the forms of nonlinearity in the original FPU model. One interesting feature of this regime is that smallamplitude intrinsic localized modes (ILMs, which are also often called “discrete breathers”) Campbell:Today2004 (); Flach:PR2008 () of the bright type (i.e., on top of a nonvanishing background) do not exist ^{2}^{2}2However, on top of a nonvanishing background, dark breathers can arise cho1 (); cho2 () from the linear limit. in the absence of disorder because of the specific relations between the parameters , , and in Eq. (15). This phenomenon was discussed in Ref. Theocharis:PRE2009 () based on the consideration of modulational instabilities (MIs) of linear waves due to nonlinearity. An MI is a generic mechanism to generate such localized waves from linear waves at band edges of a linear spectrum. However, to have an MI, it is necessary that , which is not satisfied in the present case. Nevertheless, introducing impurities in a granular chain leads to the emergence of breatherlike “defect” solutions that bifurcate from linear impurity modes Theocharis:PRE2009 ().
ii.2.3 “Strongly Nonlinear” Regime ()
When precompression is sufficiently weak in comparison to the strains (and for vanishing precompression), one can no longer approximate Eq. (1) by truncating a Taylor expansion. In general, for materials in which the sound speed goes to or remains very small, it is not reasonable to use a standard linear approximation as a starting point for a perturbative analysis Nesterenko:book (). This is particularly interesting from the point of view of transport and localization theory in nonlinear disordered systems, because almost all of the research in the field has focused on the influence of nonlinearity for disordered systems in which the linear spectrum is initially either full of or partially full of localized states Flach:ArxivRep20141 (); Flach:ArxivRep20142 (). Consequently, understanding the interplay between disorder and nonlinearity in the sonicvacuum limit brings new theoretical challenges, and — as we shall see — it also produces a fundamentally distinct form of dynamics.
As a starting point towards developing a theory for transport and localization in granular crystals, we nevertheless start by focusing our efforts in a standard way by extending the linear theory to the nonlinear regime.
(N)  

ii.3 Physical Parameters
We take advantage of the numerous experimental investigations of granular crystals Sen:PR2008 () to incorporate physically meaningful values for the parameters in Eq. (1). We suppose that all the spheres are made of steel, and we use the parameters given in Theocharis:PRE2009 () unless we specify otherwise. In particular, the elastic modulus is GPa, the Poisson ratio is , and the density is kg/m. We also suppose that the elastic plates at the boundaries have the same mechanical properties as the spheres, so and . In this paper, we examine disordered bidisperse granular chains, and we choose the radii of the spheres to be mm and , where . Note that reduces the system to the case of a homogeneous chain. To explore different nonlinear regimes, we use in the range between N and N. As we discussed in Subsection II.2, the amount of nonlinearity in the dynamics of each bead depends on the ratio . In Table 1, we show the initial value of the ratio for a homogeneous chain with an initially localized displacement excitation with m.
Iii Types of Disordered Configurations
In this article, we study three qualitatively different configurations of disordered granular chains: an Andersonlike configuration in which adjacent sites are uncorrelated and two types of random dimer models (RDMs) that include correlations across sites.
We consider bidisperse granular chains, so each chain consists of some configuration that includes two possible types of spheres: type 1 has radius , and type 2 has radius , where . We also suppose that all of the spheres are made from the same material, so their elastic properties are the same. That is, , , and the density is the same. We take their masses and to be different. Note that reduces the system to the case of a homogeneous chain. We consider three different ways of distributing the particles to produce disorder: (1) an Andersonlike distribution that amounts to uncorrelated disorder, (2) a random dimer distribution that follows a choice in Ref. Dunlap:PRL1990 () (RDM1), and a random dimer distribution that follows the choice in Ref. Ponson:PRE2010 () (RDM2). We show all three types of disorder in Fig. 1, and we note that both RDM1 and RDM2 are correlated types of disorder. See Izrailev2012125 () for a review of localization in systems with correlated disorder.
Specifically, we construct our three families of disordered chains as follows:

Anderson (A): For each of the particles in the chain, choose radius with a probability of and radius with a probability of . (Importantly, note that the choice for each particle is independent of all other particles.)

Random dimer model 1 (RDM1): For each of the dimers in the particle chain, we choose the configuration (i.e., both particles have a radius of ) with a probability of and the configuration with a probability of . (In the literature, this family specifically is what is usually meant by the term “random dimer model” Dunlap:PRL1990 ().)

Random dimer model 2 (RDM2): For each of the dimers in the particle chain, we choose the configuration with a probability of and the configuration with a probability of . (That is, we are choosing the orientation of the dimer, which we imagine to be a spin with two possible states Ponson:PRE2010 ().)
Because these granular chains include two types of spheres (and are oriented horizontally, so we can ignore gravity), there are three types of sphere–sphere interactions:

(between two spheres of radius ),

(between two spheres of radius ),

(between spheres of different radii).
One can characterize the disorder using two parameters. The parameter defines the extent of disorder. Thus, and are fully ordered cases, and is the most disordered case ^{3}^{3}3Note, however, that there are specific chain configurations — such as a periodic sequence of dimers with alternating spin orientations in RDM2 — for which gives a maximal order with respect to higherorder correlations. Upon averaging over many configurations with the same value of , such situations contribute little to the expected dynamics due to their low probability of occurrence.. The other parameter is , which is deterministic and defines the strength of the disorder by affecting the inertia (via the mass) of the particles and the magnitude of the interactions and .
It is worth remarking that in the original Anderson configuration, the radius of the th particle is , where corresponds to some uncorrelated sequence, such that and is the disorder strength. The “Anderson” model that we study (which is more precisely designated as “Andersonlike”) is an example of a “random binary alloy” Evangelou:PRB1993 () that has the same correlation properties as the original Anderson model. In the most general case, a random binary alloy can also include correlations due to dimer terms like the ones in RDM1 and RDM2 Evangelou:PRB1993 (). We study the correlation properties for each type of disorder in the next subsection.
iii.1 Correlations
Let be a random vector generated by the rules that we described above. Without loss of generality, we label each entry of as or . Thus, the Anderson chain has vector components of
(16) 
where . For the dimer models, we have
(17) 
where , for RDM1, and for RDM2.
Let be the cyclic permutation of that satisfies
(18) 
To characterize the amount of correlation in each case, we calculate the Pearson correlation coefficient
(19) 
where is the covariance between and , the vector has all elements equal to , the standard deviation of the vector is , and the mean of is given by the expectation . For our computations, it is convenient to write the covariance as
(20) 
Note that some statistical properties, such as the mean and standard deviation, are independent of permutations (i.e., and ), so we can write Eq. (19) as
(21) 
where one can calculate the term in parentheses as the sum of conditional probabilities that depend on the type of disorder.
In the next three subsubsections, we calculate the correlation coefficients for each type of disorder in the thermodynamic (i.e., ) limit.
iii.1.1 Anderson
In a granular chain with an Andersonlike disorder, the mean value of is , and the standard deviation is . Both quantities depend on the probability , but the correlation
(22) 
is independent of . When or , the correlation becomes because the granular chain is homogeneous. Equation (22) implies that the Anderson disorder is an uncorrelated type of disorder. This is true exactly in the thermodynamic limit (i.e., as the number of particles ). However, it is also true in an average sense for finite systems, which implies that the mean of the Pearson correlations for a large number of finite systems approaches the value of the correlation for a single system as . In Fig. 2(a), we show the Pearson correlation coefficient as a function of the relative distance between spheres. The black curve shows the mean value, which tends to (as we just discussed).
iii.1.2 Random Dimer Model 1 (RDM1)
As in the Anderson case, the mean value of for the RDM1 granular chain is , and the standard deviation is . However, because an RDM1 granular chain consists of a sequence of dimers, there is now a shortrange correlation. The Pearson correlation coefficient is
(23) 
which we note is again independent of . Consequently, the chain has the above shortrange correlation between secondnearestneighbors neighbors for any .
In Fig. 2(b), we show the Pearson correlation coefficient as function of the relative distance between spheres.
iii.1.3 Random Dimer Model 2 (RDM2)
The RDM2 granular chain has rather different statistical properties from the other two types of disordered chains.
The mean value of is , and the standard deviation is . Both the mean and the standard deviation are independent of the probability , because affects only the orientation of the the dimer; the numbers of values and values are unchanged. This type of disorder includes a longrange correlation that one can tune with the parameter . The Pearson correlation coefficient is
(24) 
An interesting special case occurs when , as the correlation reduces to a shortrange anticorrelation:
(25) 
Other interesting limits are the ordered diatomic chains that arise for and . Because the orientation of the dimer units is constant in these limits, there is a perfect correlation between particles that are an even distance apart and a perfect anticorrelation between particles that are an odd distance apart:
(26) 
In Figs. 2(c,d), we show the Pearson correlation coefficient as a function of the relative distance between spheres. We use different values of for the two panels.
Iv Impurities in a Homogeneous Granular Chain
Reference Theocharis:PRE2009 () confirmed for granular chains the general notion that either localized or resonant modes arise when an otherwise homogeneous system (a socalled “host” chain) includes impurities. The nature of such modes depends on the relation between the parameters of the impurities and those of the other spheres in the host chain. If an impurity mass is smaller (respectively, larger) than the rest of the particles, then the associated mode is localized (respectively, resonant). Moreover, it is possible to extend localized linear modes into the weakly nonlinear regime using a continuation procedure.
Impurities also break the translational symmetry of a chain, which implies that scattering processes around the impurities play a significant role in the dynamics. This becomes increasingly important as the number of impurities in a chain increases. To emphasize the role of impurities in the transport and localization properties of a system, we highlight the socalled “random dimer model” Dunlap:PRL1990 (), which we call RDM1 in the present article. In Ref. Dunlap:PRL1990 (), it was shown for the Schrödinger lattice with onsite energy distributed in an RDM1 way that — even when almost all of the linear modes are spatially localized — there is always one mode that is extended for a certain value (which depends on the strength of the disorder) of the wavenumber . Furthermore, for finite 1D systems, there is a set of modes for wavenumbers near (in particular, for wavenumbers , with as ) that have a localization length that is larger than the length of the system Datta:PRB1995 ().
A similar effect from double impurities has been observed in acoustic chains with harmonic interactions Datta:PRB1995 (). However, due to the acoustic characteristics of the linear spectrum, the frequency linear mode is extended in either a homogeneous or a disordered system (independently of the type of disorder) and for any system length. Therefore, even for an Andersonlike disorder configuration, modes with and (as ) have a localization length that is larger than the size of the system Datta:PRB1995 (); Lepri:PRE2010 (). We expect reflectionless modes to emerge in strongly compressed granular chains — i.e., in the linear regime (see Sec. II.2.1), in which a harmonic approximation of the interactions is suitable.
V Numerical Results
In general, it is difficult to precisely determine localization properties in disordered systems — primarily because most tests are based on the asymptotic behavior of particular observables (e.g., energy). From a practical perspective, one needs to consider long chains (and large volumes in larger dimensions) and very long integration times, and (from a theoretical perspective) one should let both time and system size go to infinity Kramer:RPP1993 (). Such scenarios are difficult to achieve experimentally, and even numerical simulations pose considerable difficulties Skokos:PRE2009 (). In particular, one is often interested in the asymptotic behavior of the energy distribution. Hence, to conduct longtime simulations without significant (and unphysical) variation in a system’s total energy, it is necessary to employ carefullychosen numericalintegration schemes. Additionally, because we are examining disordered systems and we thus need to average over a large number of realizations of a particular type of disorder to obtain appropriate statistical power, it is also necessary to employ sufficiently fast numericalintegration schemes that are also particularly accurate in their energy conservation. We thus use a symplectic integrator from Refs. Laskar:Celest2001 (); Skokos:PRE2009 (); skokoserratum ().
We also rely on indirect methods to develop intuition about the asymptotic behavior of disordered granular chains. One such method is to study the structure of the linear spectrum and the extent of localization of the linear modes. For instance, in the classical Anderson model in a 1D electronic system Anderson:PR1958 (), all of the linear modes are localized exponentially for any amount of disorder. This leads to an absence of diffusion that manifests as a saturation of the second moment of the probability distribution as a function of time. In other words, excitations remain spatially localized. By contrast, as we mentioned in Sec. IV, the RDM1 Dunlap:PRL1990 () behaves differently from the Anderson model in this respect, as the former includes extended modes that cause the second moment to grow as a function of time.
In our ensuing discussions, we investigate the influence of the three different types of disorder on the structure of the linear spectrum and the presence of localized states in both the bulk and the surface of a granular chain. We subsequently investigate transport and dynamical localization in the bulk for disordered Hertzian chains (1).
v.1 Direct Diagonalization of Eq. (12)
There are various ways of measuring localization in linear modes. In finite systems, it is useful to calculate the inverse participation ratio (IPR) Kramer:RPP1993 ()
(27) 
where some certain distribution. For modal analysis, we use , which allows one to measure the fraction of particles whose displacement of position from equilibrium differs markedly from . We can thereby measure the extent of localization. For instance, a plane wave with all sites equally excited satisfies as the number of particles . By contrast, a strongly localized wave satisfies , and exactly when only one sphere is vibrating (i.e., when for all and ).
Calculating the IPR makes it possible to directly obtain a qualitative understanding of the nature of the linear modes. In Fig. 3, we show the spectrum and the extent of localization (i.e., its IPR) associated with the linear modes for one realization of each of the three types of disorder. In Fig. 4, we show the mean value of the IPR over 100 realizations of each type of disordered chain as a function of the probability parameter and the size parameter . In both figures, we have sorted the modes from smallest frequency to largest frequency. Diagonalizing Eq. (12) directly yields the displacement distribution of the particles in the chain that are associated with the different modes. In this section, we use these displacement distributions to compute the IPRs that we show in Figs. 3 and 4. We also evaluate Eq. (27) using the energydensity distribution (given by Eq. (29), as we will discuss in Sec. V.3), and we obtain qualitatively similar results. For each type of disorder, we will use the energydensity distribution (see Sec. V.4) to characterize the dynamical localization.
We first consider Andersonlike disorder. For frequencies [see Eq. (10) and Fig. 3], we observe a complicated gap structure that includes isolated frequencies between the two bandedge frequencies. In the frequency range , there is also a small region in which has multiple peaks with values that are close to . These peaks are associated with singlenode impuritylike modes, in which the energy oscillates primarily around one particle. As was discussed in Ref. Theocharis:PRE2009 (), linear localized modes are bound to small particles for a single impurity, and the frequency of these modes is larger than the lower edge frequency of the homogeneous chain. Additionally, for a given precompression force , the frequency depends only on the strength of the impurity, and it thus depends only on the size parameter . There are also modes with that are related to double impurities. More precisely, is slightly smaller than because the mode does not consist exactly of two particles that vibrate, as there is also a tail that decays as a function of space. Modes with a lower IPR are associated with different local configurations. For example, a mode with two small masses that vibrate with a large amplitude and are separated by a large mass that oscillates with a small amplitude has . Additionally, modes that have spheres that effectively participate in the system dynamics, for , have , and one can make analogous statements for other values of .
One can interpret the probability parameter as a measure of the density of small impurities (i.e., particles with radius ) in a host chain of particles with radius . As , the granular chain is composed almost exclusively of spheres with radius , and its few small impurities generate impurity modes whose frequencies are larger than . The rest of the spectrum consists mostly of an acoustic branch that is bounded above by . This explains why the Anderson chain with in Fig. 4 has an IPR whose maximum occurs near the maximum mode number (i.e., it is close to the frequency edge ). When decreases, the fraction of particles with radius increases, and the population of modes with frequencies between and increases as well. In particular, the maximum value of in Fig. 4 in the Anderson case (which occurs for ) is about , which implies that most localized linear modes are double impuritylike modes instead of single impuritylike modes. However, the frequency of these modes does not change for a fixed value of , and it is close to the frequency edge at .
Another interesting feature of the Anderson model in granular chains is that the frequency mode is extended for all values of and . In other words, it is independent of the amount of disorder and of the relative sizes of the two types of particles Datta:PRB1995 (); Lepri1 (). Near , there is a nontrivial region in the  parameter space in which one observes extended modes in a finitesize chain. One expects the area of this region to vanish as the system size Datta:PRB1995 (). However, the presence of this extended mode opens a channel for the transportation of energy even in a disordered chain.
For an RDM1 chain, the frequency structure is similar to that of an Anderson chain. However, there are several highfrequency modes, which each have frequency between and , that form an almost flat structure in plots of frequency versus mode number (see Fig. 3). These frequencies are related to quasidegenerate modes, which have almost the same frequency as each other, and such modes arise more often in an RDM1 chain than in an Andersonlike chain. As in the Andersonlike chain, an RDM1 chain also includes some highly localized linear modes that are related to double impurities. Nevertheless, the main difference arises in the distribution, which for an RDM1 chain includes an extra minimum near a frequency of that depends on the parameters and . For example, when , we obtain roughly kHz for and the physical parameters described in Sec. II.3. This is related to extended modes that are centered at a nonzero frequency. Furthermore, as one can see from Fig. 4, the IPR tends to be smaller for most values of and in an RDM1 chain in comparison with an Anderson case. This occurs because the impurities in RDM1 chains are twice as large as those in Anderson chains, which implies in turn that RDM1 chains have large impurity modes.
An RDM2 chain exhibits completely different — and rather remarkable — features in its spectrum and IPR distribution from the other two types of disordered chains. To explain these differences, it is important to interpret the RDM2 system as a perturbation of a perfectly ordered diatomic chain instead of as a perturbation of a monoatomic one. In fact, most of the eigenvalues for an RDM2 chain occur between the frequency edges of the ordered diatomic chain (i.e., within its pass bands). The rest of the eigenvalues are organized predominantly into almost flat distributions within the band gaps (see Fig. 3). An RDM2 chain tends to have more degenerate modes than an RDM1 chain. RDM2 chains also have very interesting localization properties. In Fig. 4, for example, we observe that the distributions are (on average) almost independent of the degree of disorder (i.e., on the parameter ). We also see from Fig. 3 that most of the degenerate modes are also equally localized. In other words, they have almost the same value of . To explain the features of the IPR, observe that there exist a few single impuritylike modes with , but most of the localized modes consist of two (associated with ), three () or four () vibrating particles. Additionally, the RDM2 disorder is symmetric with respect to by construction (so, e.g., and are equivalent situations).
v.2 Spreading and Partial Localization Due to Disorder and Nonlinearity
Force distributions are particular useful in granular crystals, because it is easier and more reliable to measure forces than to measure energy. Moreover, examining forces as a function of time allows one to indirectly measure spreading and localization. Thus, in this section, we examine how the force evolves at specific spots in the chain and also how the force distributions are affected by changes in the precompression in homogeneous and Andersonlike disordered chains.
In Fig. 5, we show example force distributions from applying an excitation that consists of an initially localized displacement at the center of a homogeneous chain. For a strongly compressed chain (e.g., for N), the initial excitation spreads along the chain, and the dynamics arise from the decomposition of the Kronecker into linear modes. However, the spreading is slightly asymmetric, because the nonlinearity cannot be neglected entirely. Increasing the nonlinearity in the system by decreasing the precompression leads to a lessening of the distribution width due a decrease in the system’s sound speed. One directly observes this effect in the force distribution, and one can also see it indirectly by examining the force at different places in the chain as a function of time. For instance, the time that takes to detect fluctuations in the force at particles 30 and 60 sites away from the position of the initial excitation becomes longer as one decreases the precompression. Additionally, in the sonicvacuum regime, solitary waves emerge clearly, and the energy is divided mainly into two pulses that move in opposite directions.
In the presence of disorder (see Fig. 6), we observe that the spatial force distribution changes abruptly (i.e., even for a small amount of disorder) from the distribution in an associated homogeneous chain. When linear effects are dominant (e.g., at N), the force distribution has a maximum near the position of the initial excitation, and it decays exponentially away from this point. Near the central position of the distribution, the temporal force dynamics includes largeamplitude, persistent oscillations that exist for long times. The forces in particles that are a few sites away from the center (e.g., see the particles that are and sites away from the center in Fig. 6) exhibit oscillations whose amplitudes are ordersofmagnitude lower.
When we increase the effective nonlinearity in a granular chain — in particular, in the weakly nonlinear situation, such as the one in Eq. (14) — resonances of linear modes are induced by nonlinear shifts of the frequencies Flach:ArxivRep20141 (). This leads to a nonlinear mechanism of energy exchange between the localized and extended modes in the spectrum (see Section V.1), which in turn implies that energy that was previously stored in localized modes can now be carried through the system by being transferred either to other localized modes that are spatially close to the original one or to extended modes. In short, there is more transport. Consequently, the force is distributed among a larger number of particles in the chain. This effect is analogous to phenomena that have been observed in disordered NLS and KG lattices Flach:ArxivRep20141 (), and analogous dynamics has also been observed experimentally in the context of waveguide arrays Heinrich:NJP2012 (); Naether:OL2013 (). Remarkably, the localization goes away completely when the precompression goes to , and instead a pure spreading process occurs. In other words, the localization phenomenon, in which nearly all of the energy at vanishing precompression would be partitioned into localized traveling waves (which each have a support on only a few site of the chain) hinch (), is modified drastically because the presence of disorder.
v.3 Energy Distribution and Second Moment
As we stated previously, characterizing whether or not dynamics is localized — and which particular transport properties can characterize localization in a quantitative way — is a difficult task Kramer:RPP1993 (), and it has been examined from many different perspectives by several authors. Such methods include (1) computing a localization length LL1 (); LepriLya (), which gives information on how fast the distributions decay; (2) computing finitetime Lyapunov exponents Lyapunov1 (); Johansson:EPL2010 () to study KAM tori and chaotic dynamics; (3) directly estimating scaling properties of the energy distribution Lepri:PRE2010 (); Mulansky:NJP2013 (); and (4) calculating moments of distributions that are associated with the dynamics Dunlap:PRL1990 (); Datta:PRB1995 (); Kopidakis:PRL2008 (); Lepri:PRE2010 (); Naether:OL2013 (). The calculation of moments has been especially popular, and it is particularly common to investigate the growth of the second moment as a function of time, as this gives information about the width of a distribution. However, the exclusive use of the second moment as a singleparameter description is problematic and can lead to a misunderstanding of a system’s actual dynamics Flach:ArxivRep20141 (); Flach:ArxivRep20142 (), particularly in strongly nonlinear situations. Consequently, following Flach:ArxivRep20141 (); Flach:ArxivRep20142 (), in the present work, we examine dynamics by computing not only the second moment but also the IPR (see Section V.4).
Proceeding with our analysis, we note that the total energy of the system is conserved by the dynamics. We are thus interested in the energy distribution’s second moment
(28) 
where is the energy density of the th particle and is the position of the center of the distribution. The energy density of the th particle is given by
(29) 
where
(30) 
is the particle’s kinetic energy and the potential energy depends on the model. For example, in the linear limit, the potential energy is
(31) 
In the weakly nonlinear regime, the potential energy is
(32) 
In the strongly nonlinear regime of a Hertzian potential, the potential energy is
(33) 
The two last terms in the righthand side of the Hertzian potential energy of Eq. (V.3) have minus signs, so in the weakly nonlinear limit and in the linear limit. The first term in Eq. (V.3) gives only a trivial contribution to the total energy, because it corresponds to the (constant) background energy associated with the precompression. The last term in Eq. (V.3) is a telescopic series when one considers all , and the boundaries do not play any significant role because we are interested in the bulk dynamics. The displacement and the momentum at the edges of the chain are both exactly for all times.
In the linear regime and in the absence of disorder, the only possible situation after a very long time is for the system to thermalize Lepri1 (); Lepri2 (), so one obtains equipartition of energy between the different degrees of freedom. As a result (and as is wellknown), the asymptotic spreading dynamics in a homogeneous chain is ballistic (i.e., as ) regardless of whether the initial condition is a local displacement perturbation (i.e., ) or a local velocity perturbation (i.e., ) Landau:StatPhys (); Datta:PRB1995 (). However, introducing either disorder or nonlinearity can drastically change transport properties Flach:ArxivRep20141 (). For example, attempting to estimate a scaling relationship for the second moment now typically produces a different exponent: as , where .
However, one can expect even more complicated phenomena, so in particular it is not always meaningful to fit the spreading of the second moment to a power law with a single exponent Lepri:PRE2010 (). When there is reasonable powerlaw scaling, the behavior is called “superdiffusive” when , “diffusive” when , and “subdiffusive” when . There is no diffusion when . Following the work by Lepri et. al. Lepri:PRE2010 (), we attempt to identify the situations in which it is reasonable to construe the second moment as having a powerlaw scaling by using as a diagnostic the logarithmic derivative,
(34) 
where we calculate as a mean over some number of different realizations of the disorder. We expect that when as , but that can exhibit oscillations when the dynamics are more complicated. In our numerical computations, we estimate the logarithmic derivative using the finitedifference approximation , where we discretize time as in our numerical integration. The criterion that we use to state when has a powerlaw scaling is
(35) 
with a small parameter and an arbitrary time within our observation horizon. We thereby separate the cases in which oscillations of the numerical data for the second moment are admissible as statistical fluctuations from the ones in which oscillations are larger than statistical fluctuations.
It is also useful to compute the spectral density associated with the dynamics, as that allows one to identify which frequencies are involved in the dynamics Naether:PRE2013 (). We use the spatiotemporal displacement distribution to calculate the normalized spectral density
(36) 
where
and we use the time points to partition the interval into uniform subintervals.
In Sec. V.1, we discussed the effects of disorder in strongly precompressed chains of spheres, and we showed that disorder splits the spectrum into a lowfrequency region (in which the modes are extended) and a highfrequency region (in which modes tend to be localized). We now seek to explore the interplay between disorder and nonlinearity in both the stronglyprecompressed (i.e., weakly nonlinear) regime and the strongly nonlinear regime (whose limiting case is a sonic vacuum). We integrate Eq. (1) numerically using a “” algorithm Laskar:Celest2001 (); Skokos:PRE2009 (); skokoserratum (), which is a symplectic integrator that allows one to conserve energy for long temporal evolution. Using , the relative error in energy is between and (depending on the simulation parameters) using a reasonably small time step of s.
v.3.1 DisplacementPerturbation Initial Conditions
In Fig. 7, we show the spatiotemporal energy distribution and the spectral density for both homogeneous and Andersonlike disordered chains for different levels of precompression and for a displacementperturbation initial condition , with m. For N, we see that the main contribution to the dynamics comes from the linear modes (as we discussed previously). For the homogeneous case, maxima at nonzero frequencies give the bandedge frequencies, where the linear spectrum is denser than it is near frequency. When one decreases the precompression, the band width decreases, and the spreading of waves from the linear modes becomes slower because the sound speed also decreases. One observes clear nonlinear pulses in the dynamics, and the speed of these pulses is larger than the sound speed for sufficiently small precompression. (See, for instance, the panels in Fig. 5 with N.) However, for , the localized initial condition splits into traveling pulses that propagate in opposite directions. This occurs because all of the frequencies of the linear spectrum tend to for .
The chain with Andersonlike disorder exhibits more complicated dynamics than the homogeneous chain. In Fig. 7, we observe Andersonlike localization for strong levels of precompression. Spikes in the spectral density indicate the modes that contribute the most to the dynamics. The highest spike is located at a high frequency, so the main contribution comes from a localized mode (see Fig. 3) that is presumably close to the position of the initial excitation. As we can see from the lowfrequency spikes in the spectral density, the localization process occurs on top of a diffusive background pattern that arises primarily because of extended modes. As we consider weaker precompression, we observe a narrower frequency band near frequency in the spectral density, analogous to our observations for homogeneous chains. Although there have been many efforts to study the interplay between disorder and nonlinearity — and their effect on spreading dynamics — most prior research has concentrated on weakly nonlinear settings. In fact, the majority of prior work has concentrated on NLS and KG lattices (see, e.g. Ref. Flach:ArxivRep20141 () and references therein). It has been observed in these settings that transport is typically subdiffusive. A notable example in which neighboring lattice sites are not coupled linearly was investigated recently in Ref. Mulansky:NJP2013 (), who considered strongly nonlinear lattices in which both the onsite and the intersite interactions are nonlinear. However, those systems also exhibits subdiffusive spreading. We believe that the contribution of the onsite nonlinearity is crucial for obtaining subdiffusive spread in lattices with Andersonlike disorder, as the energyspreading exponents in Mulansky:NJP2013 () differ considerably from the ones that we identify in the present work. Although effects from nonlinearity and disorder can separately localize energy — and, indeed, that is their general predilection, as we can see in Figs. 7(c,q) — exactly the opposite can occur in some situations that include both of these factors [see Figs. 7(o,s)]. In particular, we find when both disorder and nonlinearity are present that it is possible for spreading to be enhanced rather than for the two features to conspire to create additional localization.
For granular chains in the strongly nonlinear regime, neither localization in the form of intrinsic localized modes nor exact localization as traveling nearlycompact waves is possible, as each of these structures is destroyed by disorder. It is also impossible to localize in an Andersonlike way, as such localization is suppressed by nonlinearity and the absence of a linear limit. Instead, the energy spreads among the particles in a peculiar but characteristic way: strongly localized (and nearly compact) waves are still present at the edges of the energy distribution during the spreading process at N [see Fig. 7(o)]; however, for N, the disorder induces multiple scattering events, which causes the wave amplitudes to decrease [see Fig. 7(s)].
v.3.2 VelocityPerturbation Initial Conditions
To analyze the dynamics for an initial velocity perturbation, we consider , and we set these perturbations to have the same energy as with the initial displacement perturbation . To get as a function of (or viceversa) one needs to solve . For example, to express the velocity perturbation in terms of the displacement perturbation, we write
(37) 
Thus, in our numerical simulations, we set m/s, which is the value that we obtain for a homogeneous chain with N and m.
In Fig. 8, we show spatiotemporal energy distributions and spectral density for both homogenous and Andersonlike chains using an initially localized velocity perturbation. The main — and fundamental — difference compared to what we observed using displacementperturbation initial conditions (see Fig. 7) comes from the spectral density. When there is strongly precompression, we observe that the distribution of modes that are excited by the velocityperturbation initial condition is denser near frequency than it is elsewhere. In the disordered case, this implies that the mean contribution to the dynamics comes from extended modes rather than localized modes. This contrasts starkly with our observations using displacement perturbations, and it leads to dynamics in which the energy spreads much faster than for displacement excitations. Moreover, for velocity perturbations, the energy that diffuses in the background is comparable to the amount of energy that remains localized. For N, we observe in both homogeneous and Andersonlike chains that a solitary wave propagates faster than the spreading pattern [see Figs. 8(i,k)]. For weaker precompression, the solitary wave still propagates in the homogeneous chain, but its amplitude decays in an Andersonlike chain. In particular, when , the solitary waves are delocalized due to scattering with defects in the disordered chain, and the energy pattern that emerges is qualitatively similar to what was observed in Ponson:PRE2010 () for transport of solitary waves in the RDM2 case in a highdisorder regime.
To visualize what happens to the energy from the dynamics in Andersonlike chains, we average the energy distribution at s over 100 realizations. In Fig. 9, we show this energy distribution using a logarithmic scale for both displacement excitations and velocity excitations. We examine how the distribution changes depending on the strength of nonlinearity. Specifically, we observe that the energy distribution grows exponentially near the edges of a chain. For and N, this occurs in a narrow region (fewer than 30 sites) of the chain, and energy is localized at the edge of the distribution because traveling waves survive the disorder. This phenomenon is considerably less prominent for displacementperturbation initial conditions than for velocityperturbation initial conditions. For N (sonicvacuum regime), we also observe exponential growth of the energy distribution near the chain edges. In this case, however, it occurs over a wider region (about 100 sites for displacement excitations and about 150 sites for velocity excitations), and the exponential growth has a considerably lower exponent than in the chains with nonzero precompression.
v.4 Transport Arising from Nonlinearity
To quantitatively characterize transport and localization processes, we conduct longtime simulations — up to s — in chains with spheres. We use long chains to avoid boundary effects during the entire numerical integration; no waves reach the boundary of the system within the simulation time. We compute the second moment and the IPR as functions of time for the three types of disorder, and we average our results over 500 realizations of a chain configuration in each case. To confirm our numerical results, we conduct several tests. For example, we compare our results from with those using a Runge–Kutta scheme with a very small time step (between s and s), and we obtain quantitatively the same results for the same realization of disorder ^{4}^{4}4Importantly, using allows much longer simulation times and a significant improvement in energy conservation in comparison to using a Runge–Kutta scheme.. We also test the scheme using smaller time steps ( s and s) and larger system sizes ( and particles), and we again obtain the same results. In the current section, we compute the IPR [see Eq. (27)] using the energy distribution instead of the displacement distribution. In other words, , and and are also based on the energy distribution.