Quasiperiodic granular chains and Hofstadter butterflies

Quasicrystals are solids whose structure is ordered but not periodic [1]. For many years, it was thought 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 arises when atoms are arranged so that they possess symmetries, such as $5$-fold symmetry, that are forbidden to periodic crystals¹¹1By the crystallographic restriction, crystals can only have certain rotational symmetries: $2$-fold, $3$-fold, $4$-fold, or $6$-fold symmetry [3].. A famous two-dimensional (2D) example is a Penrose tiling [4]. More generally, quasiperiodic structures can exist in any number of dimensions as structures with a broken translational symmetry. In 1D, the most common models used in the study of quasiperiodic 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 to continuously deform one into another without closing any bulk gaps.

A key feature of quasiperiodic potentials, arising prominently in the study of Schrödinger equations [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 the experiments 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-body systems, discrete systems, and settings with controllable interactions. For example, there have been relevant investigations in both bosonic and fermionic settings [12, 13].

In the present paper, we study the effects of quasiperiodicity in strongly precompressed granular chains [14, 15], in which neighboring particles interact with each other according to a Hertzian potential. Our aim is to illustrate that localization of eigenmodes can occur in quasiperiodic granular chains and to explore the conditions — with respect to both models and experimental setups — in which it occurs. We illuminate these conditions by comparing a variety of different models with one or both of on-site and off-site quasiperiodic structures. As demonstrated in a wealth of research over more than three decades [16], granular chains are extremely versatile, as one can adjust interaction potentials; readily tune them between almost linear, weakly nonlinear, and strongly nonlinear regimes by applying different precompression strengths; construct them using particles of different sizes, shapes, and material properties; and so on. This has yielded a wealth of insights about a diverse set of physical phenomena, including acoustic realizations of many concepts from condensed-matter physics [17]. Most centrally for our discussion, this includes the dynamics of wave transport and localization in disordered nonlinear systems in both theoretical  [18, 19] and experimental [20] studies. Other phenomena (both theoretical and applied) from condensed-matter and quantum physics that have been realized in granular crystals include an analog of a Ramsauer–Townsend resonance in the form of a square well [21], switching and rectification [22], and others. More broadly, granular chains provide 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]) and thermodynamic (e.g., equipartition [26]) phenomena.

The remainder of our paper is organized as follows. In Section 2, we briefly review the AA model. In Section 3, we present three models of 1D quasiperiodic lattices: two with on-site quasiperiodicity and one with off-site quasiperiodicity. In Section 4, we linearize the governing equations of our three models. In Section 5, we demonstrate that such models may possess a Hofstadter butterfly structure. In Section 6, we examine energy transport and localization in our models. We conclude and suggest some interesting directions for future research in Section 7.

The prototypical form of the Aubry–André (AA) model at the level of a tight-binding model is

where $\Psi$ is a wavefunction, $n$ indexes the lattice site, $E$ is energy, and $V$ is a potential. We suppose that the on-site energy is modulated by a lattice distortion with wave vector $q\in\mathds{R}$, which is incommensurate with $2\pi$ and has phase $r$. We also suppose that $V(x)=V(x+2\pi n)$, where $x\in\mathds{R}$ and $n\in\mathds{Z}$.

In [6], Aubry and André proved several fundamental properties of the eigenmodes of (1). They showed that a ground state exists and that it undergoes a transition from analyticity for $\lambda<\lambda_{c}$ to nonanalyticity for $\lambda>\lambda_{c}$ for certain $\lambda_{c}(q)$ and when $q$ is not a Liouville number. That is, $q\in\mathds{R}\setminus\mathds{Q}$ and there exist $\gamma$ and $r$ such that

is satisfied for any rational number $p_{1}/p_{2}$.

Aubry and André also showed, using a perturbative approach, that the analyticity-breaking causes the eigenmodes of (1) to have very rich spatial properties. Specifically, one can write the eigenmodes of the modulated system as

where $v_{m}$ are the coefficients of the Fourier expansion of $V$. The eigenmodes and eigenvalues for $\lambda=0$ are given by $e^{ikn}$ and $2\cos(k)$, respectively. For the series in (3) to be convergent, one needs to satisfy a Diophantine condition. That is, there exist two positive constants $K$ and $\beta$ such that

for any integers $m$ and $n$.

For the sake of simplicity, consider the special case in which the phase $r=0$ and $V(\xi n)=\cos(2\pi\xi x)$, where $\xi$ equals the golden ratio $(1+\sqrt{5})/2$. So, Eq. (1) takes the form: $E\Psi_{n}=\Psi_{n+1}+\Psi_{n-1}+\lambda\cos(2\pi\xi n)\Psi_{n}$. In this case, Aubry and André showed that for $\lambda>\lambda_{c}=2$, all of the eigenmodes of (1) are exponentially localized, as in the Anderson model (in which the potential $V$ is disordered rather than quasiperiodic)[27], with the same characteristic localization length

That is, $\Psi_{n}$ decays asymptotically as $e^{-n/\zeta}$ as $n\rightarrow\infty$. However, when $\lambda<2$, most eigenmodes are given by extended, modulated plane waves. Interestingly, this implies that the loss of analyticity is also associated with a transition at $\lambda_{c}$ to spatial localization of the eigenmodes, where $\lambda_{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 diverse variety of physical systems [11, 28, 29, 30], including photonic lattices, Bose–Einstein condensates, and many others. Additionally, the spectra of the corresponding Schrödinger operators (with two or even more frequencies) is a topic of intense mathematical interest; see, e.g., [31] and references therein.

Several recent studies have generalized conventional granular chains in various ways. They have yielded several interesting insights, and they promise to result in a host of others in the coming years [16, 17]. One type of a generalized granular chain is a cradle system [32], in which particles are attached to linear oscillators, enabling the use of on-site potentials in a way that is independent of particle–particle interactions. Several potential realizations of such a setting were given in [33] (although they have yet to be implemented experimentally, to the best of our knowledge). Another fascinating variant arises from examining a chain of particles with internal resonators, such as by coupling a secondary particle inside a principal one. This leads to 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, 35, 36]. Additionally, the use of particles with non-spherical geometries can drastically modify particle–particle interactions. For instance, with cylindrical particles, although one has the same functional relation between the force and the displacements as with spherical particles, one can tune the magnitude of such interactions by changing the contact angle between adjacent cylinders [37, 38].

We will consider granular chains with all of the above types of variations. The equations of motion in our general setting are

where $u_{n}$ is the displacement of the $n$th particle (where $n\in\{1,2,{\ldots},N\}$) measured from its equilibrium position in the initially compressed chain, $v_{n}$ is the displacement of the $n$th interior mass (when one particle is located inside another), and

is a static displacement for each particle that arises from the static load $F_{0}=\text{const}$. There is a Hertzian interaction between a pair of particles 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

In our subsequent discussions (see Section 3), we will consider various special cases of (6) and (7), depending on the type of particle that we use to construct chains. Specifically, we work with three different models: two of them have spherical particles, and one has cylindrical particles. For our analysis and computations, we assume that $\alpha_{n}$, $\beta_{n}$, and $\gamma_{n}$ vary sinusoidally in space according to the following formulas:

where $\xi$ is the golden mean $\frac{\sqrt{5}+1}{2}$ (unless we explicitly state otherwise), $\bar{\alpha}_{1}>\bar{\alpha}_{2}\geq 0$, and $\bar{\beta}$, $\bar{\gamma}\geq 0$. For simplicity, we separately examine the effects of (10)–(12). This leads to three different models in which it may be possible to observe an AA transition. In Eqs. (10)–(12), $\bar{\alpha}_{2}$, $\bar{\beta}$, and $\bar{\gamma}$ are the quasiperiodic parameters that determine the strengths of the modulations for the different terms in Eqs. (6) and (7).

Depending on the strength of the precompression that we apply to a granular chain and the magnitude of the strains that arise in (or are exerted on) the chain, one can expand the Hertzian force in a power series about the equilibrium state. This process reduces the equations of motion for 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 $\delta_{n}\gg|u_{n-1}-u_{n}|$ for all $n$ — the dominant terms are the linear ones, so we linearize Eq. (18) to get

where

By considering (without loss of generality at the level of this linear approximation) a complex representation of the wavefunctions, $u_{n}=\phi_{n}e^{i\omega t}$ and $v_{n}=\psi_{n}e^{i\omega t}$, and defining $E=-\omega^{2}$, we obtain

an eigenvalue problem that we can solve numerically by diagonalization. Using this linear description, we can now address the issue of a localization (“metal–insulator”) transition for suitable incommensurate periodic coefficient variations of different types.

Another property of the AA model’s spectrum is its fractal nature. To explore this, we compute the spectrum as a function of $\xi$ (see Fig. 5), and we observe a structure that is known as a Hofstadter butterfly. The butterfly is a footprint of the spectrum’s fractality, and one can see its statistical self-similarity in the figure.

The Hofstadter butterfly was first predicted in 1976 [47] in a completely different system: Bloch electrons 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 to observe a Hofstadter butterfly. As a result, it took until 1997 — in a microwave system [48] — that a Hofstadter butterfly was observed experimentally. A Hofstadter butterfly was then observed in 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 the simplicity and controllability of the latter system.

To characterize the self-similarity of the spectrum in the different cases, we compute the Minkowski–Bouligand fractal dimension ($D_{m}$), which (by Moran’s theorem) is the same as the Hausdorff dimension ($D_{h}$) for strictly self-similar fractals [51, 52]. The procedure to numerically compute $D_{M}$ is known as box-counting, and we proceed as follows. First, we map the Hofstadter spectrum into a square of $480\times 480$ pixels, we then partition the square into boxes of characteristic size (side length) $l_{B}$, and finally we count the number $N_{B}$ of boxes that include at 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}\propto l_{B}^{D_{M}}$. In practice, one computes $D_{M}$ as the best fit to $N_{B}\propto l_{B}^{D_{M}}$. For this, we compute a linear regression of the logarithm of the data using gradient descent. In Fig. 6(a), we show an example of the box counting for the Hofstadter spectrum of model Ia at the localization transition (i.e., when $\bar{\beta}=2$). To show how the band gaps are filled with boxes, we have superimposed the $N_{B}$ boxes over the spectrum for different values of $l_{B}$.

We now compute the fractal dimension $D_{M}$ as a function of the quasiperiodic parameters $\bar{\alpha}_{2}$, $\bar{\beta}$, and $\bar{\gamma}$ for our three models. We expect $D_{M}$ to be between $1$ and $2$, because $D_{M}=1$ for a line and $D_{M}=2$ for a plane. In Fig. 6(b), we show the results of our computations. We observe that the minimum fractal dimension occurs at the same point as the localization transition for models Ia (at $\bar{\beta}=2$) and Ib (at $\bar{\gamma}=1$). We calculate that $D_{M}\gtrapprox 1.69$ for model Ia and $D_{M}\gtrapprox 1.82$ for model Ib. It is interesting to note the similar non-monotonic dependence of the fractal dimension on the model parameters in models Ia and Ib. Presumably, this arises from the aforementioned similarity of the former model and the out-of-phase excitations of the latter model. In contrast, in model II, for $\bar{\alpha}_{1}=3$ and $\bar{\alpha}_{2}\in(0,3)$, we observe that the fractal dimension decreases monotonically as $\bar{\alpha}_{2}$ increases.

Another important issue, which one can examine in several ways, is how energy transport is affected by the quasiperiodicity [53, 54]. For instance, one can compute a second moment $m_{2}$ of the energy distribution as a function of time to quantitatively characterize the temporal evolution of the energy distribution’s width [18, 19, 20, 55, 56, 57, 58, 59, 60]. Following recent studies in disordered granular chains [20], we study the evolution of the energy distribution $\{H_{n}\}_{n=1}^{N}$ immediately after the impact of a striker against the first particle ($n=1$). We write the particles’ energy as

The total energy is $H=\sum_{n}H_{n}$, which is a conserved quantity. Note that, in the context of experiments, one should also consider dissipation which is neglected here; see a relevant summary of models thereof in [17]. Importantly, we expect that our results will be robust enough to be observable experimentally even in the presence of a small amount of dissipation. This claim 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]) in granular chains.

We then compute the second moment

and we estimate a scaling relationship between $m_{2}$ and $t$. When the scaling is approximated reasonably as a power law — for which, in exact form, $m_{2}(t)\sim t^{\eta}$ as $t\rightarrow\infty$ for some exponent $\eta$ — one can categorize as ballistic the case when $\eta=2$, superdiffusive the one when $\eta\in(1,2)$, diffusive when $\eta=1$, and subdiffusive when $\eta\in(0,1)$. If $\eta=0$, we say that there is no diffusion; in other words, all energy remains localized. In a perfectly homogeneous granular chain, it is known that energy transport is ballistic. However, when disorder is added to the chain, the dynamics can change drastically, and one can observe different energy transport regimes [18, 19, 20]. These previous studies have focused on the interplay between disorder and nonlinearity. Here, in contrast, we show that even in a strongly precompressed (i.e., almost linear) chain, one can obtain any desired exponent diffusion $\eta\in\left[0,2\right]$ for the energy transport. However, we find that on-site potentials are essential to have localization.

One advantage of working with quasiperiodic chains instead of disordered ones is that we do 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 $\xi$, one gets one specific chain. This enables us to cover the whole parameter space with considerably fewer computations than when studying disordered chains. We are also interested in characterizing energy transport in realistic frameworks, so we set $N=21$, which gives a long enough chain to qualitatively capture the nature of transport, at least in several recent experimental and theoretical explorations [18, 20].

To integrate (6), we use a fifth-order explicit Runge–Kutta (RK5) method with a step size of $\mathrm{d}t=0.01$. We set $\delta_{n}=1$ for all $n$. We also set $u_{0}(t)=u_{N+1}(t)=0$ for all $t$, so we have fixed boundary 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 use the former condition to stop the code in cases in which all of the energy is trapped in the form of localized states. In other words, there is no diffusion. This is expected, for instance, in model Ia for $\bar{\beta}>2$ (i.e., after the localization transition occurs). However, it is an uncommon scenario in model II, which does not have a localization transition. In Fig. 7, we show our results for energy transport in our three models. In models Ia and Ib, we explore the parameters ranges $\bar{\beta}\in[0,4]$ and $\bar{\gamma}\in[0,3]$, respectively, so we can compare energy transport on both sides of the localization transition. In model II, we consider $\bar{\alpha}_{1}=3$ and $\bar{\alpha}_{2}\in[0,3]$. For all three models, and for $\xi\in[0,1]$, we can tune the energy transport over from subdiffusive to ballistic behaviors. This allows a great deal of control of the energy-transport properties, and it is remarkable that we are able to do this using a deterministic model. In model Ia, we also observe localization (for which $\eta=0$) in contrast with observations in disordered granular crystals [18, 19, 20]. This suggests that the inclusion of on-site potentials is crucial for this localization phenomenon.

In this article, we introduced different types of quasiperiodic granular chains that were inspired by 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 the localization 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 one of 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 (as in a cradle) or in the form of a local resonator (for a mass in mass system) — and model II uses cylindrical 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 and with 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 fractal properties by computing the Minkowski–Bouligand fractal dimension. In models Ia and Ib, we showed that the minimum fractal dimension $D_{M}$ of the spectrum coincides with the point at which the localization transition occurs. For model II, we observed that the spectrum’s fractal dimension 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 subdiffusive transport — as well as localization. In contrast to prior work, which achieved such control using a combination of disorder and adjusting the nonlinearity strength [18, 19, 20], we were able to control the 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 laboratory experiments, as one can then further explore the role that granular systems can play in the study of quasiperiodic operators [48, 49, 50]. Achieving an experimental realization of a Hofstadter butterfly in a granular chain would also be very exciting in its own right. Among the settings that we have proposed in this paper, model Ib (i.e., the chain of mass-in-mass particles) is the clearest candidate for observing a localization transition (in the out-of-phase variables), given that the cradle system has yet to be experimentally realized. Arguably, the woodpile setup of [44] may also constitute an excellent playground for such studies. However, a key consideration, given experimental limitations, is that one seeks to build a chain with as few particles as possible such that one can (still) observe the relevant phenomenology.

There are also numerous open issues to explore computationally and theoretically. Examples include the effect of nonlinearity (e.g., through larger excitation amplitudes) on these modes and their localization, how these phenomena differ for granular crystals in different numbers of dimensions, and others. Relevant extensions will be considered in future studies.

AJM acknowledges support from CONICYT (BCH72130485/2013). PGK gratefully acknowledges support from the US-AFOSR under FA9550-17-1-0114.

We thank R. Chaunsali, Alain Goriely, Robert MacKay, and Francisco J. Muñoz for helpful comments. We also thank the editors of this special issue for their invitation to contribute an article to it.