Elastic Wannier-Stark ladders and Bloch oscillations in 1D granular crystals

The study of energy localization in materials has increasingly drawn attention of researchers in recent years.^1,2 On these lines, granular crystals, which are tightly packed arrangements of granular particles interacting as per nonlinear contact law, have served as an attractive playground for energy localization due to their tunability.^3–5 Such systems are fertile test beds of studying elastic wave propagation as a result of intermixing effects of dispersion, nonlinearity, defect, and disorder. The localization of wave—in the form of the spatially localized linear or nonlinear resonant modes—has thus been shown in the studies related to breathers,^6–8 defects,^9,10 topological interfaces,¹¹ and disorder.^12–15 

Here, we explore a new way of achieving wave localization in granular systems by the formation of the elastic Wannier-Stark ladders (WSLs) and Bloch oscillations (BOs). BOs were originated from quantum mechanics referring to the spatially localized oscillations of electrons in a periodic potential with an externally applied electrical field. After those were predicted by Bloch and Zener,¹⁶ this phenomenon remained controversial for more than 60 years.¹⁷ Chynoweth et al. observed WSLs in crystals with an applied electric field and verified the Bloch-Zener model.¹⁸ These WSLs referred to the equally-spaced electron energy bands and represented the frequency-domain counterparts of BOs. After the initial experimental verifications of electric BOs in semiconductor superlattice,¹⁹ analogous systems have been studied and demonstrated for plasmonic waves,²⁰ light waves,^17,21,22 and acoustic waves.^23–27 WSLs and BOs were also observed in elastic systems for torsional^28,29 and surface³⁰ waves. However, experimental studies on elastic BOs are still limited to specific experimental designs that are not easily tunable, and this fact hinders further studies and potential applications to some extent.

In this research, we numerically and experimentally demonstrate the existence of WSLs and the corresponding time-resolved BOs by using a tunable one-dimensional (1D) granular system consisting of a chain of uniform cylinders.³¹ According to the Hertz contact law,³² the contact stiffness between two adjacent cylinders depends on the contact angle between them, and this fact brings extreme tunability to our system. Previous research has shown that such cylindrical granular crystals provide an enhanced freedom in tailoring stress wave propagation in solids.^31,33–39 For the current study, we can therefore vary the alignment of contact angles in a tunable fashion and introduce a gradient in contact stiffness along the chain—mimicking the effect of an externally applied electric field for BOs of electrons. This results in the tilting of frequency band structure—an essential feature for the formation of the WSLs. We use a discrete element model to analyze the dynamics of the system and show that it has a good agreement with the experimental observations. Furthermore, we numerically demonstrate that such systems can provide pathways for further tailoring of elastic BOs.

The experimental setup is shown in Fig. 1. The vertical 1D chain is composed of N = 21 cylindrical granular particles made of fused quartz with density ρ = 2187 kg/m³, elastic modulus E = 72 GPa, and Poisson's ratio ν = 0.17. All the cylinders are identical with the same radius R = 9 mm and height H = 18 mm. We use 3D-printed enclosures [white color casing in Fig. 1(b) and enlarged view in Fig. 1(c)] to hold the cylindrical particles and stack them vertically. All the enclosures, along with the cylinders inside, can be rotated freely about their central axis. Taking advantage of the calibration markers on the outer surface of the enclosures, we are able to achieve a fine control over the rotation angle, and thus maintain contact angles with a resolution as small as $2.5°$⁠.

In the present study, we manipulate the contact angles between the particles to create a stiffness gradient along the length of the chain. To this end, we first take a dimer chain, i.e., a chain with two alternating contact angles, ${90°, 40°}$⁠, and then modify the same judiciously to achieve a gradient. We set the first six contact angles to be ${90°, 5°, 90°, 7.5°, 90°,15°}$⁠, followed by alternating ${90°, 40°}$ until the end of the chain. Note that we count the particles from the actuation side [i.e., from the bottom of the chain in Fig. 1(a)]. The primary reason behind choosing these set of angles is to achieve a large stiffness gradient at the beginning of the chain, so that WSLs and BOs are located near the actuator for an effective excitation.

We compress the system with a static force F₀ = 100 N. We mount a static force sensor at the upper bracket to measure this force. The piezoelectric actuator at the bottom employs two different types of excitation signals: (i) a linear chirp signal (2 to 30 kHz) for band structure estimation, and (ii) Gaussian-modulated sinusoidal pulse (GMSP) at a target central frequency for time-history response. The full-width-half-magnitude width of the GMSP signal is 0.82 ms. We also place a piezoelectric force sensor at the top end of the chain—in direct contact with the uppermost particle—to measure the transmitted dynamic force [see Fig. 1(a)]. The maximum transmitted dynamic force is about 0.55 N, which is approximately 0.55% of the initial static compression. We use a Laser Doppler Vibrometer (LDV) to measure the velocity of each particle one by one and obtain system's response. Here, the laser passage holes designed on each enclosure [Fig. 1(c)] allow us to point the laser at $45°$ from the vertical direction and measure a component of particle velocity. We synchronize the velocity profiles measured from all particles with respect to the actuation signals to form a spatio-temporal map of BOs.

We model the dynamics of our system by the discrete element method,⁴⁰ in which we treat the chain as a system of lumped masses interconnected by nonlinear springs. This nonlinear interaction arises from the nature of the contact between cylinders, and it is governed by the Hertz's contact law.³² If α be the contact angle between two adjacent cylinders, the contact area would be elliptical for $0°<α<90°$ and circular for $α=90°$⁠. The nonlinear force-displacement relationship is thus given as³³ 

where F is the contact force, and $z=δ+Δu$ is the total relative displacement of two adjacent cylindrical particles, and it comprises of the initial static compression δ and the relative dynamic displacement Δu [see Fig. 1(d)]. The contact stiffness coefficient β(α) is of the following form:

where the eccentricity ε of the elliptical contact area is $1−(r2r1)2$ with r₁ and r₂ as the semi-major and semi-minor axes, respectively. $K(ε)$ and $E(ε)$ are the complete elliptic integrals of the first and second kind, respectively, and $r2r1≈(1+cos (α)1−cos (α))−23$⁠.

Therefore, we write the equation of motion for the jth cylinder in the chain as

Here, u_j is the dynamic displacement of the jth particle, and δ_j is the initial static compression between the jth and (j + 1)th particles due to the precompressive force. β_j is the simplified notation for β(α_j) used hereafter [see its dependence on the contact angle in Fig. 1(e)]. The mass of each cylinder is m = ρπR²H, and [z]₊ denotes max(0, z). The final term represents the viscous force that is dependent on velocity $uj̇$ and damping coefficient τ. Under a large initial static compression, we neglect the effect of gravity in this short chain. This is because the additional static compression by gravity at the bottom of the chain is about 2% of the applied precompression, which leads to only 1.1% variation in band-edge frequencies investigated later.

To obtain the dispersion characteristics of the system, we employ an analytical formulation by linearizing Eq. (3). This is because the initial static force (displacement) is large compared to instantaneous dynamic force (displacement) in the system; therefore, we have $|uj+1−uj|≪δj$⁠. Moreover, neglecting the effect of viscous damping on the eigenfrequencies of our system, we suppress the damping term. We thus deduce

where k_j is the linearized contact stiffness between the jth and (j + 1)th particles and is given as $kj=32βjδj1/2=32βj1/2F01/3$ with F₀ being the initial static compression force.⁴¹ 

Now, if we have two alternating contact angles along the chain, the system is a “dimer” system with two alternating contact stiffness values, say k_s (soft contact) and k_h (hard contact). Bloch's theorem⁴² can be employed on this periodic system to obtain a representative frequency band structure shown in Fig. 2(a). This includes two bands, namely, acoustic and optical bands,⁴³ which are flat along the spatial extent of the chain, and the band-edge frequencies are given by⁴³ 

Our system, however, as we have mentioned earlier, is a variant of this dimer chain, and we keep one contact angle uniform along the chain (i.e., same k_s) but vary the other angle to achieve a gradient in k_h as a function of space—shown in the upper inset of Fig. 2(b). Although this system is no longer a periodic system, an approximate dispersion behavior can still be estimated by calculating the local band-edge frequencies. To this end, we take a dimer unit cell with two stiffness values at different space locations and calculate the band-edge frequencies from Eq. (5). This approach is similar to the one taken in the seminal paper by Zener.¹⁶ As a result, the optical band is tilted and becomes space dependent [Fig. 2(b)].

We carry out numerical simulations by directly solving the set of equations given by Eq. (3). To account for the boundary conditions in our finite system, we write the equation of motion for the first and last particles as

where u₀ is the displacement of the actuator, k_a is the linear contact stiffness between the actuator and the first particle, and k_e is the linear contact stiffness between the last particle and the dynamic force sensor, which is assumed to be fixed. The linear contact assumption is valid for high initial static compression. To match the mode shapes of the boundary modes (front and end) observed in experiments, we set the stiffness coefficients k_a = 4.2 × 10⁷ N/m and k_e = 1.5 × 10⁸ N/m. Moreover, we choose the damping factor τ = 0.0125 ms to match the decay of short input Gaussian pulse in the experimental results to be discussed later. δ_a (δ_e) is the initial static compression between the actuator (dynamic force sensor) and the first (last) cylinder.

We rewrite the equations in the form of first-order state space and solve them using the ODE45 solver in MATLAB. We obtain the velocity response with a sampling time step of 10^–6s. The excitation conditions are linear chirp signal and GMSP for the band structure and time-history analyses, respectively, given at the front of the chain in the form of the actuator displacement u₀ in Eq. (6a).

WSLs are the frequency-domain counterparts of BOs. To directly observe those, we calculate the frequency band structure of the system under a chirp excitation by performing the Fast Fourier Transformation (FFT) on time-history response at each cylinder location. The power spectral density (PSD) of the velocity response obtained by direct numerical simulation is shown in Fig. 3(a). The colored areas denote the presence of energy and the white area indicates the presence of local band gaps. For the chosen contact angles in the current setup (detailed in Sec. II), we observe a titled optical band in the beginning of the chain, which matches well with the analytically obtained local band-edge frequencies [cf. Eq. (5)] shown as dashed black lines. Within this inclined optical band, we observe WSLs in the form of three localized resonance peaks, marked as WSL1 (20.42 kHz), WSL2 (18.26 kHz), and WSL3 (15.79 kHz). In Fig. 3(b), we show the experimentally obtained spectrum. We remarkably notice the presence of all three WSLs along with the titled optical band. These are WSL1 (22.80 kHz), WSL2 (19.60 kHz), and WSL3 (17.20 kHz) with relative errors of 11.66%, 7.34%, and 8.93%, respectively. Though we see the experimental results match numerics to a reasonable accuracy, an up-shift of WSL frequencies is observed in the experiments. This may be due to several factors, including (1) inaccuracy in the standard values of material properties, (2) deviation in the contact law exponent, and (3) stiffening caused by viscoelastic effects.^6,34 In Figs. 4(c) and 4(d), we plot the mode shapes at the frequencies corresponding to WSL1 and WSL2, respectively. We observe that both numerical and experimental results agree well in terms of predicting the energy localization (high amplitude) at particles {2, 3} and particles {4, 5} for WSL1 and WSL2, respectively. This localization will be utilized in Sec. IV B for exciting BOs at different spatial locations depending on the input frequency.

To observe BOs in the time domain, we measure time-history of the velocity at each particle under the GMSP excitation centered at a WSL frequency. Figures 4(a) and 4(b) show numerical and experimental time-resolved BOs at WSL1 frequency. We observe that as this short pulse is injected into the system, it does not travel deep and is primarily localized at the front end of the chain. Specifically, we observe more lasting vibrations at the second and third particles (enclosed within dashed black lines), and this complies well with our earlier prediction in Fig. 3(c) that the modal energy is localized at the same spatial location in the chain. However, given the fact that we are using a short pulse with limited energy, this localized energy decays due to dissipation in the experimental setup. We then change the central frequency of input GMSP at WSL2 and plot the results in Figs. 4(c) and 4(d). We clearly notice that the localization of energy is deeper into the chain—to the fourth and fifth particles. This also complies well with the localization observed in Fig. 3(d) for WSL2. This therefore marks a direct observation of the time-resolved elastic BOs, in which we have shown a remarkable change in spatial localization of elastic energy as a function of input frequency excitation in the same system.

In this section, we show how this contact-based granular chain can be further tuned to change the dispersion characteristics of the system and achieve desired spatial localization at a single operating frequency f_c. We numerically show that the contact angles in the chain can be fine-tuned in such a way that a linear [Fig. 5(a)] and quadratic [Fig. 5(b)] tilt in the optical band is obtained, which again matches well with analytically obtained curves (dashed black lines). Here, the terms linear and quadratic refer to the approximate shape of the lower boundary of the optical band. The contact angle distributions along the chain to achieve such gradients of the optical band are shown in the upper insets. These are back calculated from the desired stiffness distribution by fitting a polynomial to the contact law given by Eq. (2). We note that this includes a minute variation of angles as small as $0.3°$⁠. This fine control of contact angles was not possible with the current experimental setup. Thus, we limit our investigation only to numerical study now, while we postulate that experimental verification may be possible in the future with an improved design. We send a GMSP from the front end at a center frequency f_c = 18 kHz. Time-history plots of velocity in Figs. 5(c) and 5(d) show distinct spatial and temporal characteristics of BOs. More specifically, for the linear optical band case, we witness the wave localization between the 10th and the 21th particles, while for the quadratic case, the localization spot has shifted to the space between the 7th and the 15th particles. Similarly, the time period of the BOs has changed from 2.43 ms to 1.56 ms [by comparing Figs. 5(c) and 5(d)]. Lesser amplitude of velocity observed in Fig. 5(c) is due to the fact that the localization region is farther away from the actuator, and only a limited amount of energy is pumped to excite the localized mode using evanescent waves. Therefore, we confirm that by changing the contact angles in this type of granular system, we can shift BO in space and change its period in time.

In this work, we have studied wave localization phenomena in a 1D granular chain consisting of vertically stacked cylindrical particles. By tuning the contact angles between the cylinders, we introduced a gradient in contact stiffness along the chain, which resulted in a spatially tilted optical band and formed the Wannier-Stark ladders (WSLs). We used the frequencies of these WSLs to excite spatially localized resonant modes and observed their corresponding time-resolved elastic Bloch oscillations (BOs) through numerics and experiments. Using numerical experiments, we also showed the possibility of further tailoring these BOs in terms of altering their spatial and time responses if a more accurate experimental setup can be designed. We anticipate that this study triggers further studies on WSLs and BOs, especially in the context of nonlinear effects, which can easily be invoked in the current system. Moreover, such localized modes might be useful for engineering applications, for example, in energy harvesting using piezoelectric materials, where these resonances could potentially enhance mechanical-to-electrical energy conversion.

We thank Hiromi Yasuda and Hyunryung Kim at the University of Washington, Professor Eunho Kim in Chonbuk National University of Korea, Professor Chris Chong at Bowdoin College, and Professor Panayotis Kevrekidis at the University of Massachusetts Amherst for useful discussions and technical help. We would like to acknowledge the financial support from the NSF (CAREER-1553202) and the AFOSR (FA9550-17-1-0114). Last, we thank the editors for the opportunity to contribute to this special topic section.