Quasiperiodic granular chains and Hofstadter butterflies
rrsta.royalsocietypublishing.org
Research
Article submitted to journal
Subject Areas:
Nonlinear waves, nonlinear dynamics,dynamical systems,condensed-matter physics
Keywords:
Granular chains, nonlinear latticesystems, condensed-matter physics,quasiperiodicity, Hofstadter butterfly,localization
Author for correspondence:
Mason A. Portere-mail: [email protected]
Quasiperiodic granular chainsand Hofstadter butterflies
Alejandro J. Martínez , Mason A. Porter ,and P. G. Kevrekidis Oxford Centre for Industrial and Applied Mathematics,Mathematical Institute, University of Oxford, OxfordOX2 6GG, UK Department of Mathematics, University of California,Los Angeles, CA 90095, USA Department of Mathematics and Statistics, Universityof Massachusetts, Amherst, MA, 01003, USA
We study quasiperiodicity-induced localization ofwaves in strongly precompressed granular chains.We propose three different setups, inspired by theAubry–André (AA) model, of quasiperiodic chains;and we use these models to compare the effectsof on-site and off-site quasiperiodicity in nonlinearlattices. When there is purely on-site quasiperiodicity,which we implement in two different ways, we showfor a chain of spherical particles that there is alocalization transition (as in the original AA model).However, we observe no localization transition in achain of cylindrical particles in which we incorporatequasiperiodicity in the distribution of contact anglesbetween adjacent cylinders by making the angleperiodicity incommensurate with that of the chain.For each of our three models, we compute theHofstadter spectrum and the associated Minkowski–Bouligand fractal dimension, and we demonstratethat the fractal dimension decreases as one approachesthe localization transition (when it exists). Finally, ina suite of numerical computations, we demonstrateboth localization and that there exist regimes ofballistic, superdiffusive, diffusive, and subdiffusivetransport. Our models provide a flexible set ofsystems to study quasiperiodicity-induced analogs ofAnderson phenomena in granular chains that one cantune controllably from weakly to strongly nonlinearregimes. c (cid:13) The Authors. Published by the Royal Society under the terms of theCreative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author andsource are credited. a r X i v : . [ n li n . PS ] J a n r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A ..................................................................
1. Introduction
Quasicrystals are solids whose structure is ordered but not periodic [1]. For many years, it wasthought that quasicrystals could only be produced artificially. However, almost a decade ago,the first natural quasicrystal was found in Russia [2]. A common type of quasicrystal ariseswhen atoms are arranged so that they possess symmetries, such as -fold symmetry, that areforbidden to periodic crystals . A famous two-dimensional (2D) example is a Penrose tiling [4].More generally, quasiperiodic structures can exist in any number of dimensions as structureswith a broken translational symmetry. In 1D, the most common models used in the study ofquasiperiodic systems are a Fibonacci quasicrystal [5] and the Aubry–André (AA) model [6].These models are topologically equivalent to each other [7], in the sense that it is possible tocontinuously deform one into another without closing any bulk gaps.A key feature of quasiperiodic potentials, arising prominently in the study of Schrödingerequations [6], is the transition from a “metallic” phase, in which all eigenstates are delocalized,to an “insulating” phase, in which they are localized. See, e.g., the analysis in [8,9] and theexperiments in [11]. It is of considerable interest to extend these studies in various ways,including to nonlinear systems (see, e.g., the work of [10] and references therein), many-bodysystems, discrete systems, and settings with controllable interactions. For example, there havebeen relevant investigations in both bosonic and fermionic settings [12,13].In the present paper, we study the effects of quasiperiodicity in strongly precompressedgranular chains [14,15], in which neighboring particles interact with each other according toa Hertzian potential. Our aim is to illustrate that localization of eigenmodes can occur inquasiperiodic granular chains and to explore the conditions — with respect to both modelsand experimental setups — in which it occurs. We illuminate these conditions by comparing avariety of different models with one or both of on-site and off-site quasiperiodic structures. Asdemonstrated in a wealth of research over more than three decades [16], granular chains areextremely versatile, as one can adjust interaction potentials; readily tune them between almostlinear, weakly nonlinear, and strongly nonlinear regimes by applying different precompressionstrengths; construct them using particles of different sizes, shapes, and material properties; andso on. This has yielded a wealth of insights about a diverse set of physical phenomena, includingacoustic realizations of many concepts from condensed-matter physics [17]. Most centrally forour discussion, this includes the dynamics of wave transport and localization in disorderednonlinear systems in both theoretical [18,19] and experimental [20] studies. Other phenomena(both theoretical and applied) from condensed-matter and quantum physics that have beenrealized in granular crystals include an analog of a Ramsauer–Townsend resonance in the formof a square well [21], switching and rectification [22], and others. More broadly, granular chainsprovide a wonderful playground that enables systematic exploration of the role of lattice structure(e.g., material heterogeneity [23,24]) and fundamental dynamic (e.g., rogue waves [25]) andthermodynamic (e.g., equipartition [26]) phenomena.The remainder of our paper is organized as follows. In Section 2, we briefly review the AAmodel. In Section 3, we present three models of 1D quasiperiodic lattices: two with on-sitequasiperiodicity and one with off-site quasiperiodicity. In Section 4, we linearize the governingequations of our three models. In Section 5, we demonstrate that such models may possess aHofstadter butterfly structure. In Section 6, we examine energy transport and localization in ourmodels. We conclude and suggest some interesting directions for future research in Section 7.
2. A brief review of the Aubry–André model
The prototypical form of the Aubry–André (AA) model at the level of a tight-binding model is EΨ n = Ψ n +1 + Ψ n − + λV ( nq + r ) Ψ n , (2.1) By the crystallographic restriction , crystals can only have certain rotational symmetries: -fold, -fold, -fold, or -foldsymmetry [3]. r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. where Ψ is a wavefunction, n indexes the lattice site, E is energy, and V is a potential. We supposethat the on-site energy is modulated by a lattice distortion with wave vector q ∈ R , which isincommensurate with π and has phase r . We also suppose that V ( x ) = V ( x + 2 πn ) , where x ∈ R and n ∈ Z .In [6], Aubry and André proved several fundamental properties of the eigenmodes of (2.1).They showed that a ground state exists and that it undergoes a transition from analyticity for λ < λ c to nonanalyticity for λ > λ c for certain λ c ( q ) and when q is not a Liouville number. Thatis, q ∈ R \ Q and there exist γ and r such that (cid:12)(cid:12)(cid:12)(cid:12) q π − p p (cid:12)(cid:12)(cid:12)(cid:12) > γ p r (2.2)is satisfied for any rational number p /p .Aubry and André also showed, using a perturbative approach, that the analyticity-breakingcauses the eigenmodes of (2.1) to have very rich spatial properties. Specifically, one can write theeigenmodes of the modulated system as Ψ n ( k ) = e ikn + λ ∞ (cid:88) m = −∞ v m e im ( qn + h ) mq + k ) − cos( k )] , (2.3)where v m are the coefficients of the Fourier expansion of V . The eigenmodes and eigenvaluesfor λ = 0 are given by e ikn and k ) , respectively. For the series in (2.3) to be convergent, oneneeds to satisfy a Diophantine condition. That is, there exist two positive constants K and β suchthat (cid:12)(cid:12)(cid:12)(cid:12) kπ + m q π − n (cid:12)(cid:12)(cid:12)(cid:12) > Km β (2.4)for any integers m and n .For the sake of simplicity, consider the special case in which the phase r = 0 and V ( ξn ) =cos(2 πξx ) , where ξ equals the golden ratio (1 + √ / . So, Eq. (2.1) takes the form: EΨ n = Ψ n +1 + Ψ n − + λ cos(2 πξn ) Ψ n . In this case, Aubry and André showed that for λ > λ c = 2 , allof the eigenmodes of (2.1) are exponentially localized, as in the Anderson model (in which thepotential V is disordered rather than quasiperiodic) [27], with the same characteristic localizationlength ζ = 1ln (cid:16) λ (cid:17) . (2.5)That is, Ψ n decays asymptotically as e − n/ζ as n → ∞ . However, when λ < , most eigenmodes aregiven by extended, modulated plane waves. Interestingly, this implies that the loss of analyticityis also associated with a transition at λ c to spatial localization of the eigenmodes, where λ c = 2 in this specific example. This transition is called a localization transition or Aubry–André transition .This result is a generic phenomenon in Schrödinger lattices, and it is thus relevant for a diversevariety of physical systems [11,28–30], including photonic lattices, Bose–Einstein condensates,and many others. Additionally, the spectra of the corresponding Schrödinger operators (with twoor even more frequencies) is a topic of intense mathematical interest; see, e.g., [31] and referencestherein.
3. Implementing the AA model in granular chains
Several recent studies have generalized conventional granular chains in various ways. Theyhave yielded several interesting insights, and they promise to result in a host of others in thecoming years [16,17]. One type of a generalized granular chain is a cradle system [32], in whichparticles are attached to linear oscillators, enabling the use of on-site potentials in a way thatis independent of particle–particle interactions. Several potential realizations of such a settingwere given in [33] (although they have yet to be implemented experimentally, to the best ofour knowledge). Another fascinating variant arises from examining a chain of particles with r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. internal resonators, such as by coupling a secondary particle inside a principal one. This leadsto a locally-resonant granular chain, which is sometimes called a mass-in-mass (or mass-with-mass , if the secondary particle is external) chain [34–36]. Additionally, the use of particles withnon-spherical geometries can drastically modify particle–particle interactions. For instance, withcylindrical particles, although one has the same functional relation between the force and thedisplacements as with spherical particles, one can tune the magnitude of such interactions bychanging the contact angle between adjacent cylinders [37,38].We will consider granular chains with all of the above types of variations. The equations ofmotion in our general setting are ¨ u n = α n [ δ n + u n − − u n ] / − α n +1 [ δ n +1 + u n − u n +1 ] / (cid:124) (cid:123)(cid:122) (cid:125) Hertzian interaction − β n u n (cid:124) (cid:123)(cid:122) (cid:125) elastic restitution − γ n ( u n − v n ) (cid:124) (cid:123)(cid:122) (cid:125) mass-in-mass interaction , (3.1) ¨ v n = γ n ( u n − v n ) , (3.2)where u n is the displacement of the n th particle (where n ∈ { , , . . ., N } ) measured from itsequilibrium position in the initially compressed chain, v n is the displacement of the n th interiormass (when one particle is located inside another), and δ n = (cid:18) F A n (cid:19) / (3.3)is a static displacement for each particle that arises from the static load F = const. There is aHertzian interaction between a pair of particles only when they are in contact, so each particle isaffected directly only by its nearest neighbors and experiences a force from a neighbor only whenit overlaps with it. This yields [ x ] + = (cid:40) x , if x > , if x ≤ . (3.4)In our subsequent discussions (see Section 3), we will consider various special cases of (3.1) and(3.2), depending on the type of particle that we use to construct chains. Specifically, we work withthree different models: two of them have spherical particles, and one has cylindrical particles.For our analysis and computations, we assume that α n , β n , and γ n vary sinusoidally in spaceaccording to the following formulas: α n = ¯ α + ¯ α cos(2 πnξ ) , (3.5) β n = ¯ β [1 + cos(2 πnξ )] , (3.6) γ n = ¯ γ [1 + cos(2 πnξ )] , (3.7)where ξ is the golden mean √ (unless we explicitly state otherwise), ¯ α > ¯ α ≥ , and ¯ β , ¯ γ ≥ .For simplicity, we separately examine the effects of (3.5)–(3.7). This leads to three different modelsin which it may be possible to observe an AA transition. In Eqs. (3.5)–(3.7), ¯ α , ¯ β , and ¯ γ are thequasiperiodic parameters that determine the strengths of the modulations for the different termsin Eqs. (3.1) and (3.2). (a) On-site quasiperiodicity: Two different variants of the AA model usingspherical particles In this section, we discuss the effect of an on-site quasiperiodicity on the dynamics of a granularchain by considering chains of spherical particles with local potentials. We set ¯ α = 0 in Eq. (3.5), r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. on - s it e po t e n c i a l (a)(b) Figure 1.
Schematics of (a) model Ia and (b) model Ib. so the coupling parameter α n is given by α n = A n /m n , where A n = 4 E n − E n (cid:16) R n − R n R n − + R n (cid:17) / (cid:104) E n (1 − ν n − ) + E n − (1 − ν n ) (cid:105) , (3.8)and the n th particle has elastic modulus E n , Poisson ratio ν n , radius R n , and mass m n . We assumethat the particles are identical, so E n = E , ν n = ν , R n = R , and m n = m . This, in turn, implies that α n = ¯ α , and we let ¯ α = 1 without loss of generality. (i) Model Ia: ¯ β (cid:54) = 0 and ¯ α = ¯ γ = 0 Suppose that the particles in the chain are attached to a mechanical restitution mechanism, suchthat there is a linear force in the equations of motion (3.1,3.2). In the limit of small angles, thiscan describe the well-known Newton’s cradle, a system in which particles are aligned in onedimension and are suspended from a ceiling by strings so that the particles collide with each otheralong one dimension and oscillate. In the top panel of Fig. 1, we show a schematic of this system.Studies of this system have focused primarily on waves that arise by releasing one of the particlesat one end with some velocity [33]. This produces a transfer of energy across the chain in the formof traveling waves [32]. Mulansky and Pikovsky studied disorder in closely related (nonlinearlycoupled, and locally linear or nonlinear) oscillator systems [39], showing numerically and using afractional-nonlinear-diffusion approach that energy transport is subdiffusive. This helps furthermotivate the study of modulated systems, such as the AA model, in granular chains. In our case,the equation of motions are ¨ u n = [ δ + u n − − u n ] / − [ δ + u n − u n +1 ] / − ¯ β [1 + cos(2 πnξ )] u n , (3.9)where ¯ β > . In other words, the Hookean spring constants are positive for all n . r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. (ii) Model Ib: ¯ γ (cid:54) = 0 and ¯ α = ¯ β = 0 We now consider chains composed of particles that include an internal degree of freedom(i.e., mass-in-mass particles) [34]. Previous studies have focused on the generation [35,36], ofsuch lattices, their traveling-wave solutions [40,41], and their (bright and dark) breather-likeexcitations [42,43]. These works illustrate that coherent structures and their dynamics are enrichedsignificantly by the presence of the internal degree of freedom (DOF). To give just one example,incorporating an internal DOF in the particles can lead to nonlocal solitary waves with non-vanishing tails (so-called “nanoptera”), which have been observed experimentally in woodpilegranular chains [44]. In our setting, we envision embedding a particle in the interior of each hostparticle, such that a particle and its interior mass are coupled via a linear restitution mechanism(such as a Hookean spring).The equations of motion, upon quasiperiodic modulation of the mass-in-mass (MiM)resonator, are ¨ u n = [ δ + u n − − u n ] / − [ δ + u n − u n +1 ] / − ¯ γ [1 + cos(2 πnξ )] ( u n − v n ) , (3.10) ¨ v n = ¯ γ [1 + cos(2 πnξ )] ( u n − v n ) , (3.11)where ¯ γ > . (b) Off-site quasiperiodicity: The AA model with cylindrical particles Another approach for incorporating quasiperiodicity in a granular chain is by tuning particle–particle interactions. In existing experimental setups, to implement an AA model, one canuse chains of cylindrical particles (rather than spherical ones). We are motivated by recentexperiments [37], where it was demonstrated that cylindrical particles offer more flexibilitythan spherical particles for tuning particle–particle interactions, as one can change the contactangle between contiguous cylinders. Moreover, spatial and even temporal (periodic) variation ofcontacts between cylindrical particles has been proposed as an efficient way for implementingvarious functionalities, including that of an acoustic transistor in [38]. (i) Model II: ¯ α (cid:54) = 0 and ¯ β = ¯ γ = 0 To give equations of motion for a chain of cylindrical particles, we first need to know the form ofthe Hertzian coefficient in this case. For identical cylinders, the interaction coefficients are [46] α n ( φ n ) = Ym g ( φ n ) (cid:104) g ( φ n ) (cid:112) g ( φ n ) (cid:105) / , (3.12)where Y depends on the physical parameters of the particles, m is the mass of a particle, and φ n is the contact angle between cylinders n − and n and it is defined mod π/ . Explicit forms of Y and g i are Y = 2 E √ R − ν ) ,g ( φ n ) = (cid:115) φ n ) K (cid:104) e ( φ n ) (cid:105) π − / ,g ( φ n ) = 4 πe ( φ ) ,g ( φ n ) = (cid:18) a b E (cid:104) e ( φ n ) (cid:105) − K (cid:104) e ( φ n ) (cid:105)(cid:19) (cid:16) K (cid:104) e ( φ n ) (cid:105) − E (cid:104) e ( φ n ) (cid:105)(cid:17) , where E is the elastic modulus, ν is the Poisson ratio, and R is the radius of the circular crosssection of the cylinders. The functions K and E , respectively, are the complete elliptical integrals r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. π /4 π /202468 0 50 100 150 2000 π /4 π /2 0 50 100 150 2000 π /4 π /2 (a) (b)(c) Figure 2. (a) Interaction coefficient α n for model II as a function of the contact angle φ n between adjacent cylinders.In panels (b) and (c), we show the contact-angle distributions for two cases: (i) ¯ α = 3 and ¯ α = 1 and (ii) ¯ α = 3 and ¯ α = 3 . We use arrows to represent the mapping process that we describe in the text. of the first and second kinds [45]. They are K ( k ) = (cid:90) π/ dθ (cid:112) − k sin θ ,E ( k ) = (cid:90) π/ (cid:112) − k sin θdθ , where e = (cid:112) − ( b/a ) is the eccentricity of the contact area, a is the semi-major axis lengthof the ellipse, and b is its semi-minor axis length. One can approximate the quotient b/a by [(1 − cos φ n ) / (1 + cos φ n )] / .This yields the following equations of motion: ¨ u n = α n ( φ n )[ δ n + u n − − u n ] / − α n +1 ( φ n )[ δ n +1 + u n − u n +1 ] / . (3.13)One can then control the interactions between particles by changing φ n [38]. This raises thefollowing question of what the distribution of contact angles { φ n } has to be to obtain α n ( φ n ) =¯ α + ¯ α cos(2 πξn ) . We address this issue by numerically inverting Eq. (3.12), so that quasiperiodicvariation of α n yields a quasiperiodic variation of angles in the interval ( φ min , φ max ) . In Fig. 2, weshow the contact-angle distributions generated in two cases: (i) when ¯ α = 3 and ¯ α = 1 and (ii) ¯ α = 3 and ¯ α = 3 . We also note that α ( φ n ) → ∞ as φ n → and that α n ( φ n ) has a lower bound at φ n = π/ .
4. Linear approximation
Depending on the strength of the precompression that we apply to a granular chain and themagnitude of the strains that arise in (or are exerted on) the chain, one can expand the Hertzianforce in a power series about the equilibrium state. This process reduces the equations of motionfor the granular chain to ones that resemble those for a Fermi–Pasta–Ulam–Tsingou (FPUT) chain[17]. In particular, if the precompression is strong enough — specifically, if δ n (cid:29) | u n − − u n | for r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. Model Ia Model Ib
Figure 3. (Top) Linear spectrum as a function of the quasiperiodicity parameter and (bottom) inverse participation ratio.We show our results for model Ia in the left column and for model Ib in the right column. all n — the dominant terms are the linear ones, so we linearize Eq. (3.13) to get ¨ u n = B n u n − + B n +1 u n +1 − ( B n + B n +1 + β n + γ n ) u n + γ n v n , (4.1) ¨ v n = γ n ( u n − v n ) , where B n = 32 α n δ / n . (4.2)By considering (without loss of generality at the level of this linear approximation) a complexrepresentation of the wavefunctions, u n = φ n e iωt and v n = ψ n e iωt , and defining E = − ω , weobtain Eφ n = B n φ n − + B n +1 φ n +1 − ( B n + B n +1 + β n + γ n ) φ n + γ n ψ n , (4.3) Eψ n = γ n ( φ n − ψ n ) , an eigenvalue problem that we can solve numerically by diagonalization. Using this lineardescription, we can now address the issue of a localization (“metal–insulator”) transition forsuitable incommensurate periodic coefficient variations of different types. (a) Linear spectrum and localization transition In addition to the linear spectrum, which we obtain by solving (4.3), we also compute the inverseparticipation ratio (IPR) P − = (cid:80) n (cid:104) h ( u n , ˙ u n ) + h ( v n , ˙ v n ) (cid:105)(cid:2)(cid:80) n h ( u n , ˙ u n ) + h ( v n , ˙ v n ) (cid:3) ∈ [0 , (4.4)as a measure of the amount of localization of the eigenmodes. For modal analysis, we use h ( u n , ˙ u n ) = u n and h ( v n , ˙ v n ) = v n . A value of P − = 1 accounts for modes when only oneparticle is vibrating. In contrast, a mode is fully extended if P − → as N → ∞ . This provides a r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. qualitative understanding of the nature of the linear modes, and a transition in the IPR also givesa way to quantitatively describe the AA transition.In Fig. 3, we show the linear spectrum and IPR as a function of the quasiperiodicity parameter(which is ¯ β in model Ia and ¯ γ in model Ib) for a chain of N = 200 particles for models Iaand Ib. We observe that these two models have a complex structure of bands and gaps, withsome frequencies that appear isolated in the gaps and others that form bands that appear tocluster. Isolated frequencies are associated with modes that are similar to impurity modes. Similarstructures of bands and gaps have been observed in other physical systems, such as in optics(see, e.g., [11,28]). Interestingly, we observe from the IPR that AA transitions occur in granularchains, most prominently in model Ia, where the transition is effectively identical to that inthe original AA model. This is a consequence of modulating only the on-site potential with anexternal mechanism, so linearizations of the two systems yield the same equations. If B n = 1 and γ n = 0 for all n , the transition occurs at ¯ β = 2 . This differs starkly from the localization propertiesof the linear modes in the Anderson model, where low-frequency linear modes remain extendedindependently of the amount of disorder [18]. In Fig. 4(a), we show the transition to localization inthe fundamental mode of model Ia. For model Ib, we double the number of modes in the system,because we double the number of DOF by incorporating the internal particles. This system hasa very rich spectrum, where the upper part (half of the modes) has the same structure as inmodel Ia, but there is also a bottom part (the other half of the modes) associated with modesthat do not undergo the localization transition and consequently are extended independently ofthe modulation. This is straightforward to explain by writing the system (4.1) in terms of in-phase( x n = u n + v n ) and out-of-phase ( y n = u n − v n ) variables. This yields ˙ x n = f ( x n , x n ± , y n , y n ± ) and ˙ y n = g ( x n , x n ± , y n , y n ± ) − γ n y n . Only the equations for ˙ y n are affected explicitly by thequasiperiodicity. In this case, γ n enters as a prefactor of y n , instead of γ n as in model Ia. Thisexplains why the localization transition occurs at ¯ γ = 1 instead of at ¯ γ = 2 . Note that modes inthe upper part of the spectrum also correspond to out-of-phase modes (between u n and v n ),whereas the bottom part of the spectrum is associated with in-phase modes. The latter do not seethe quasiperiodicity in practice (because they effectively satisfy the original Hertzian dynamicswithout the MiM contribution), so they are generically extended.In contrast to models Ia and Ib, model II does not have a localization transition; instead,we observe that all modes are extended, except for the ones that are associated with isolatedfrequencies in the gaps. In Fig. 4(b), we show the IPR as a function of ¯ α for model II with ¯ α = 3 .This suggests that, without an on-site potential, one cannot observe this sort of transition ingranular chains of cylindrical particles. In the future, it will be particularly worthwhile to explorethe generality of this conclusion. Specifically, a relevant question is whether it is generically thecase that it is impossible for inter-site interactions, modulated by one or more frequencies, toinduce a localization transition in a granular chain.
5. Hofstadter butterfly
Another property of the AA model’s spectrum is its fractal nature. To explore this, we compute thespectrum as a function of ξ (see Fig. 5), and we observe a structure that is known as a Hofstadterbutterfly . The butterfly is a footprint of the spectrum’s fractality, and one can see its statisticalself-similarity in the figure.The Hofstadter butterfly was first predicted in 1976 [47] in a completely different system: Blochelectrons on 2D lattices and in the presence of an orthogonal magnetic fields. In typical crystals,one needs to use magnetic fields that are at least of the order of several thousands of tesla toobserve a Hofstadter butterfly. As a result, it took until 1997 — in a microwave system [48] — thata Hofstadter butterfly was observed experimentally. A Hofstadter butterfly was then observedin graphene in 2013 [49] and using interacting photons in superconducting qubits in 2017 [50].The possibility to also observe Hofstadter butterflies in granular chains is very exciting, given thesimplicity and controllability of the latter system. r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. Mode number
Figure 4. (Left) Localization transition as a function of ¯ β for the fundamental mode in model Ia. (Right) Inverseparticipation ratio P − as a function of ¯ α ∈ [0 , and ¯ α = 3 for model II with N = 100 cylindrical particles. Figure 5.
Linear spectrum, in the form of a Hofstadter butterfly, as a function of ξ for (left) model Ia, (center) model Ib,and (right) model II. To characterize the self-similarity of the spectrum in the different cases, we computethe Minkowski–Bouligand fractal dimension ( D m ), which (by Moran’s theorem) is the sameas the Hausdorff dimension ( D h ) for strictly self-similar fractals [51,52]. The procedure tonumerically compute D M is known as box-counting , and we proceed as follows. First, we mapthe Hofstadter spectrum into a square of × pixels, we then partition the square into boxesof characteristic size (side length) l B , and finally we count the number N B of boxes that includeat least one point of the spectrum. We do this procedure for different values of l B ; if ln( N B ) scales linearly with ln( l B ) , then the fractal dimension D M satisfies the relation N B ∝ l D M B . Inpractice, one computes D M as the best fit to N B ∝ l D M B . For this, we compute a linear regression r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. ● ● ● ●●● ●●●●●●●●●●●● ●●● ● ● ● l B N B ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ □□□□□□□□□□□□□□□□□□□□□□□□□□□□□ α n , β n , γ n D m ● Model I ▲ Model II □ Model III
Model IaModel IbModel II (a) (b)
Figure 6. (a) Example of the box-counting process that we use to compute the fractal dimension for model Ia and ¯ β = 2 .(b) Minkowski–Bouligand fractal dimension D M as a function of the quasiperiodic parameter ¯ β , ¯ γ , and ¯ α for models Ia,Ib, and II, respectively. We use the parameter value ¯ α = 3 for model II. of the logarithm of the data using gradient descent. In Fig. 6(a), we show an example of the boxcounting for the Hofstadter spectrum of model Ia at the localization transition (i.e., when ¯ β = 2 ).To show how the band gaps are filled with boxes, we have superimposed the N B boxes over thespectrum for different values of l B .We now compute the fractal dimension D M as a function of the quasiperiodic parameters ¯ α , ¯ β , and ¯ γ for our three models. We expect D M to be between and , because D M = 1 for a lineand D M = 2 for a plane. In Fig. 6(b), we show the results of our computations. We observe thatthe minimum fractal dimension occurs at the same point as the localization transition for modelsIa (at ¯ β = 2 ) and Ib (at ¯ γ = 1 ). We calculate that D M (cid:39) . for model Ia and D M (cid:39) . for modelIb. It is interesting to note the similar non-monotonic dependence of the fractal dimension on themodel parameters in models Ia and Ib. Presumably, this arises from the aforementioned similarityof the former model and the out-of-phase excitations of the latter model. In contrast, in model II,for ¯ α = 3 and ¯ α ∈ (0 , , we observe that the fractal dimension decreases monotonically as ¯ α increases.
6. Energy transport and localization
Another important issue, which one can examine in several ways, is how energy transport isaffected by the quasiperiodicity [53,54]. For instance, one can compute a second moment m ofthe energy distribution as a function of time to quantitatively characterize the temporal evolutionof the energy distribution’s width [18–20,55–60]. Following recent studies in disordered granularchains [20], we study the evolution of the energy distribution { H n } Nn =1 immediately after theimpact of a striker against the first particle ( n = 1 ). We write the particles’ energy as H n = ˙ u n v n (cid:26) α n δ n + u n − − u n ] / + 2 α n +1 δ n +1 + u n − u n +1 ] / (cid:27) (6.1) + β n u n + γ n u n − v n ) . The total energy is H = (cid:80) n H n , which is a conserved quantity. Note that, in the contextof experiments, one should also consider dissipation which is neglected here; see a relevantsummary of models thereof in [17]. Importantly, we expect that our results will be robust enough r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. to be observable experimentally even in the presence of a small amount of dissipation. Thisclaim is supported by recent experimental results on the observation of other linear phenomena(e.g., Anderson-like localization [20] and an analog of a Ramsauer–Townsend resonance [21]) ingranular chains.We then compute the second moment m ( t ) = (cid:80) n ( n − H n (cid:80) n H n , (6.2)and we estimate a scaling relationship between m and t . When the scaling is approximatedreasonably as a power law — for which, in exact form, m ( t ) ∼ t η as t → ∞ for some exponent η — one can categorize as ballistic the case when η = 2 , superdiffusive the one when η ∈ (1 , , diffusive when η = 1 , and subdiffusive when η ∈ (0 , . If η = 0 , we say that there is no diffusion;in other words, all energy remains localized. In a perfectly homogeneous granular chain, it isknown that energy transport is ballistic. However, when disorder is added to the chain, thedynamics can change drastically, and one can observe different energy transport regimes [18–20].These previous studies have focused on the interplay between disorder and nonlinearity. Here, incontrast, we show that even in a strongly precompressed (i.e., almost linear) chain, one can obtainany desired exponent diffusion η ∈ [0 , for the energy transport. However, we find that on-sitepotentials are essential to have localization.One advantage of working with quasiperiodic chains instead of disordered ones is that wedo not need to compute averages over a large number of realizations to obtain robust insights.Our quasiperiodic chains are produced in a deterministic way, so given a parameter ξ , one getsone specific chain. This enables us to cover the whole parameter space with considerably fewercomputations than when studying disordered chains. We are also interested in characterizingenergy transport in realistic frameworks, so we set N = 21 , which gives a long enough chain toqualitatively capture the nature of transport, at least in several recent experimental and theoreticalexplorations [18,20].To integrate (3.1), we use a fifth-order explicit Runge–Kutta (RK5) method with a step sizeof d t = 0 . . We set δ n = 1 for all n . We also set u ( t ) = u N +1 ( t ) = 0 for all t , so we have fixedboundary conditions on both sides. We use the stopping criterion for our simulations that either T = 20 or that energy reaches the boundary opposite to the one impacted by the striker. We usethe former condition to stop the code in cases in which all of the energy is trapped in the form oflocalized states. In other words, there is no diffusion. This is expected, for instance, in model Ia for ¯ β > (i.e., after the localization transition occurs). However, it is an uncommon scenario in modelII, which does not have a localization transition. In Fig. 7, we show our results for energy transportin our three models. In models Ia and Ib, we explore the parameters ranges ¯ β ∈ [0 , and ¯ γ ∈ [0 , ,respectively, so we can compare energy transport on both sides of the localization transition. Inmodel II, we consider ¯ α = 3 and ¯ α ∈ [0 , . For all three models, and for ξ ∈ [0 , , we can tunethe energy transport over from subdiffusive to ballistic behaviors. This allows a great deal ofcontrol of the energy-transport properties, and it is remarkable that we are able to do this using adeterministic model. In model Ia, we also observe localization (for which η = 0 ) in contrast withobservations in disordered granular crystals [18–20]. This suggests that the inclusion of on-sitepotentials is crucial for this localization phenomenon.
7. Conclusions
In this article, we introduced different types of quasiperiodic granular chains that were inspiredby the work by Aubry and André in condensed-matter physics and by recent developments(cradle and mass-in-mass systems) in the context of granular lattice systems. We studied thelocalization and spectral properties of such chains. To achieve each type of quasiperiodic chain,we incorporated spatial modulation (which is incommensurate with the chain’s period) into oneof the physical parameters. We proposed three models: models Ia and Ib use spherical particles,in which quasiperiodicity enters via an on-site potential — either in the form of a local oscillator r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A ..................................................................
210 21 21 d i ff u s i on e xpon e n t a) b) c) Figure 7.
Diffusion exponent γ for (a) model Ia, (b) model Ib, and (c) model II as a function of the quasiperiodic parametersand ξ . (as in a cradle) or in the form of a local resonator (for a mass in mass system) — and model II usescylindrical particles. In models Ia and Ib, we demonstrated the existence of an analog of the well-known AA transition. However, in model II, we showed that, without an on-site potential andwith quasiperiodicity affecting only inter-site interactions, a localization transition cannot occur.We also computed the Hofstadter spectrum for each of the models and studied their fractalproperties by computing the Minkowski–Bouligand fractal dimension. In models Ia and Ib, weshowed that the minimum fractal dimension D M of the spectrum coincides with the point atwhich the localization transition occurs. For model II, we observed that the spectrum’s fractaldimension decreases monotonically as a function of the quasiperiodic parameter.Finally, we numerically studied energy transport by exciting the granular chains with a striker.We demonstrated the existence of different regimes — ranging from ballistic to subdiffusivetransport — as well as localization. In contrast to prior work, which achieved such control using acombination of disorder and adjusting the nonlinearity strength [18–20], we were able to controlthe energy transport using a deterministic model in a strongly precompressed (and almost linear)granular chain.Naturally, it will be particularly valuable to implement some of these ideas in laboratoryexperiments, as one can then further explore the role that granular systems can play in the study ofquasiperiodic operators [48–50]. Achieving an experimental realization of a Hofstadter butterflyin a granular chain would also be very exciting in its own right. Among the settings that wehave proposed in this paper, model Ib (i.e., the chain of mass-in-mass particles) is the clearestcandidate for observing a localization transition (in the out-of-phase variables), given that thecradle system has yet to be experimentally realized. Arguably, the woodpile setup of [44] mayalso constitute an excellent playground for such studies. However, a key consideration, givenexperimental limitations, is that one seeks to build a chain with as few particles as possible suchthat one can (still) observe the relevant phenomenology.There are also numerous open issues to explore computationally and theoretically. Examplesinclude the effect of nonlinearity (e.g., through larger excitation amplitudes) on these modesand their localization, how these phenomena differ for granular crystals in different numbersof dimensions, and others. Relevant extensions will be considered in future studies.Funding. AJM acknowledges support from CONICYT (BCH72130485/2013). PGK gratefully acknowledgessupport from the US-AFOSR under FA9550-17-1-0114.
Acknowledgements.
We thank R. Chaunsali, Alain Goriely, Robert MacKay, and Francisco J. Muñoz forhelpful comments. We also thank the editors of this special issue for their invitation to contribute an article toit.
References
1. Janot C. 1994
Quasicrystals: A Primer . 2nd edition, Clarendon Press, Oxford, UK. r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A ..................................................................
2. Bindi L, Steinhardt PJ, Yao N, Lu PJ. 2009. Natural quasicrystals.
Science , 1306.3. Lifshitz R. 2003. Quasicrystals: A matter of definition.
Foundations of Physics , 1703.4. Penrose R. 1978. Pentaplexity. Eureka , 16.5. Levine D, Steinhardt PJ. 1984. Quasicrystals: A new class of ordered structures. Phys. Rev. Lett. , 2477.6. Aubry S, André G. 1980. Analyticity breaking and Anderson localization in incommensuratelattices. Ann. Isr. Phys. Soc. , 133.7. Kraus YE, Zilberberg O. 2012. Topological equivalence between the Fibonacci quasicrystal andthe Harper model. Phys. Rev. Lett. , 116404.8. Grempel DR, Fishman S, Prange, RE. 1982. Localization in an incommensurate potential: Anexactly solvable model.
Phys. Rev. Lett. , 833.9. Aullbach, C, Wobst A, Ingold, G-L, Hänggi P, Varga I. 2004. Phase-space visualization of ametal-insulator transition, New J. Phys. , 70.10. Flach S, Ivanchenko M, Khomeriki R. 2012. Correlated metallic two-particle bound states inquasiperiodic chains. EPL , , 66002.11. Lahini Y, Pugatch R, Pozzi F, Sorel M, Morandotti R, Davidson N, Silberberg Y. 2009.Observation of a localization transition in quasiperiodic photonic lattices. Phys. Rev. Lett. ,013901.12. Iyer S, Oganesyan V, Refael G, Huse DA. 2013. Many-body localization in a quasiperiodicsystem,
Phys. Rev. B , , 134202.13. Mastropietro V. 2015. Localization of interacting fermions in the Aubry–André model. Phys.Rev. Lett. , , 180401.14. Nesterenko VF. 2001 Dynamics of Heterogeneous Materials . Springer-Verlag, Heidelberg,Germany.15. Sen S, Bang J, Avalos E, Doney R. 2008. Solitary waves in the granular chain.
Phys. Rev. ,21.16. Porter MA, Kevrekidis PG, Daraio C. 2015. Granular crystals: Nonlinear dynamics meetsmaterials engineering.
Phys. Today (11), 44.17. Chong C, Porter MA, Kevrekidis PG, Daraio C. 2017. Nonlinear coherent structures ingranular crystals. J. Phys.: Condens. Matter , 413003.18. Martínez AJ, Kevrekidis PG, Porter MA. 2016. Superdiffusive transport and energylocalization in disordered granular crystals. Phys. Rev. E , 022902.19. Achilleos V, Theocharis G, Skokos Ch. 2016. Energy transport in one-dimensional disorderedgranular solids. Phys. Rev. E , 022903.20. Kim E, Martínez AJ, Phenisee SE, Kevrekidis PG, Porter MA, Yang J. 2018. Direct measurementof Superdiffusive energy transport in disordered granular chains. Nat. Commun. (in press). Seehttps://arxiv.org/abs/1705.08043.21. Martínez AJ, Yasuda H, Kim E, Kevrekidis PG, Porter MA, Yang J. 2016. Scattering of wavesby impurities in precompressed granular chains.
Phys. Rev. E , 052224.22. Boechler N, Theocharis G, Daraio C. 2011. Bifurcation based acoustic switching andrectification. Nat. Mater. , 665.23. Starosvetsky Y, Jayaprakash KR, Arif Hasan M, Vakakis AF. 2017. Topics on the NonlinearDynamics and Acoustics of Ordered Granular Media, World Scientific , Singapore.24. Lindenberg K, Harbola U, Romero AH, Lindenberg K. 2011. Pulse propagation in granularchains.
AIP Conf. Proc. , , 9725. Przedborski M, Sen S, Harroun TA. 2017. Fluctuations in Hertz chains at equilibrium. Phys.Rev. E , , 032903.26. Han D, Westley M, Sen S. 2014. Mechanical energy fluctuations in granular chains: Thepossibility of rogue fluctuations or waves. Phys. Rev. E , , 032904.27. Anderson PW. 1958. Absence of diffusion in certain random lattices. Phys. Rev. , 1492.28. Martínez AJ, Molina, MI. 2012. Surface solitons in quasiperiodic nonlinear photonic lattices.
Phys. Rev. A , 013807.29. Roati G, D’Errico C, Fallani L, Fattori M, Fort C, Zaccanti M, Modugno G, Modugno M,Inguscio M. 2008. Anderson localization of a non-interacting Bose–Einstein condensate. Nature , 895.30. Yuce C. 2014. PT symmetric Aubry–André model.
Phys. Lett. A , 2024.31. Goldstein M, Schlag W, Voda M. 2017. On the spectrum of multi-frequency quasiperiodicSchrödinger operators with large coupling, arXiv:1708.09711. r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A ..................................................................
32. James G. 2011. Nonlinear waves in Newton’s cradle and the discrete p-Schrödinger equation.
Math. Models methods Appl. Sci. , 2335.33. James G, Kevrekidis PG, Cuevas J. 2013. Breathers in oscillator chains with Hertzianinteractions, Physica D , , 39.34. Huang H, Sun C, Huang G. 2009. On the negative effective mass density in acousticmetamaterials. Ing. J. Eng. Sci. , 4.35. Kevrekidis PG, Vainchtein A, Serra Garcia M, Daraio C. 2013. Interaction of traveling waveswith mass-with-mass defects within a Hertzian chain. Phys. Rev. E , , 042911.36. Bonanomi L, Theocharis G, Daraio C., 2015. Wave propagation in granular chains with localresonances. Phys. Rev. E , , 033208.37. Khatri D, Ngo D, Daraio C. 2012. Highly nonlinear solitary waves in chains of cylindricalparticles. Granular Matter , 63.38. Li F, Chong C, Yang J, Kevrekidis PG, Daraio C. 2014. Wave transmission in time- and space-variant helicoidal phononic crystals. Phys. Rev. E , , 053201.39. Mulansky M, Pikovsky A. 2013. Energy spreading in strongly nonlinear disordered lattices. New J. Phys. , 053015.40. Xu H, Kevrekidis PG, Stefanov A. 2015. Traveling waves and their tails in locally resonantgranular systems, J. Phys. A , , 195204.41. Vorotnikov K, Starosvetsky Y, Theocharis G, Kevrekidis PG. 2018. Wave propagation in astrongly nonlinear locally resonant granular crystal. Physica D , , 27.42. Liu L, James G, Kevrekidis P, Vainchtein A. 2016. Strongly nonlinear waves in locally resonantgranular chains. Nonlinearity , , 3496.43. Liu L, James G, Kevrekidis P, Vainchtein A. 2016. Breathers in a locally resonant granular chainwith precompression. Physica D , , 27.44. Kim E, Li F, Chong C, Theocharis G, Yang J, Kevrekidis PG. 2015. Highly nonlinear wavepropagation in elastic woodpile periodic structures, Phys. Rev. Lett. , , 118002.45. Digital Library of Mathematical Functions (release 1.0.10). 2015. National Institute ofStandards and Technology. Available at http://dlmf.nist.gov/ .46. Johnson KL. 1987. Contact Mechanics . Cambridge University Press, New York, USA.47. Hofstadter DR. 1976. Energy levels and wavefunctions of Bloch electrons in rational andirrational magnetic fields.
Phys. Rev. B , 2239.48. Kuhl U, Stöckmann H-J. 1998. Microwave realization of the Hofstadter butterfly. Phys. Rev.Lett. , 3232.49. Dean CR, Wang L, Maher P, Forsythe C, Ghahari F, Gao Y, Katoch J, Ishigami M, Moon P,Koshino M, Taniguchi T, Watanabe K, Shepard KL, Hone J, Kim P. 2013. Hofstadter’s butterflyand the fractal quantum Hall effect in moiré superlattices. Nature , 598.50. Roushan P, Neill C, Tangpanitanon J, Bastidas VM, Megrant A, Barends R, Chen Y, ChenZ, Chiaro B, Dunsworth A, Fowler A, Foxen B, Giustina M, Jeffrey E, Kelly J, Lucero E,Mutus J, Neeley M, Quintana C, Sank D, Vainsencher A, Wenner J, White T, Neven H,Angelakis DG, Martinis J. 2017. Spectroscopic signature of localization with interacting photonsin superconducting qubits.
Science , 1175.51. Falconer K. 1990
Fractal geometry: mathematical foundations and applications . John Wiley & Sons,Chichester, UK.52. Schroeder M. 1991
Fractals, Chaos, Power Laws: Minutes from an Infinite Paradise . W. H. Freemanand Company, New York, USA.53. Kramer B, MacKinnon A. 1993. Localization: Theory and experiment.
Rep. Prog. Phys. , 1469.54. Laptyeva TV, Ivanchenko MV, Flach S. 2014. Nonlinear lattice waves in heterogeneous media. J. Phys. A , 493001.55. Datta PK, Kundu K. 1995. Energy transport in one-dimensional harmonic chains. Phys. Rev. B , 6287.56. Lepri S, Schilling R, Aubry S. 2010. Asymptotic energy profile of a wave packet in disorderedchains. Phys. Rev. E , 056602.57. Naether U, Rojas-Rojas S, Martínez AJ, Stützer S, Tünnermann A, Nolte S, Molina MI, VicencioRA, Szameit A. 2013. Enhanced distribution of a wave-packet in lattices with disorder andnonlinearity. Opt. Express , 927.58. García-Mata I, Shepelyansky DL. 2009. Delocalization induced by nonlinearity in systemswith disorder. Phys. Rev. E Europhys. Lett. , 30001. r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A ..................................................................
60. Rojas-Rojas S, Morales-Inostroza L, Naether U, Xavier GB, Nolte S, Szameit RA, Vicencio R,Lima G, Delgado A. 2014. Analytical model for polarization-dependent light propagation inwaveguide arrays and applications.
Phys. Rev. A90