Observation of Hofstadter butterfly and topological edge states in reconfigurable quasi-periodic acoustic crystals

The emergence of a fractal energy spectrum is the quintessence of the interplay between two periodic parameters with incommensurate length scales. crystals can emulate such interplay and also exhibit a topological bulk-boundary correspondence, enabled by their nontrivial topology in virtual dimensions. Here we propose, fabricate and experimentally test a reconfigurable one-dimensional (1D) acoustic array, in which the resonant frequencies of each element can be independently fine-tuned by a piston. We map experimentally the full Hofstadter butterfly spectrum by measuring the acoustic density of states distributed over frequency while varying the long-range order of the array. Furthermore, by adiabatically changing the phason of the array, we map topologically protected fractal boundary states, which are shown to be pumped from one edge to the other. This reconfigurable crystal serves as a model for future extensions to electronics, photonics and mechanics, as well as to quasi-crystalline systems in higher dimensions. Hofstadter’s butterfly is a fractal pattern which pictorially represents the behavior of electrons under an applied magnetic field in a 2D lattice as a pair of butterfly wings. Here, the authors recreate this pattern by measuring the acoustic density of states in a fine-tuned one-dimensional acoustic array.

T he almost Mathieu operator H ϕ;θ ¼ T þ T y þ λ cosð2πϕX þ θÞ was introduced by Peierls in 1933 1 , in his study of electron dynamics under magnetic fields. It acts on a 1D lattice, where T and X are the translation and position operators, respectively. When λ = 2, the operator connects to Harper's equation 2 , which expresses the 2-dimensional (2D) magnetic Schrödinger operator in the Landau gauge, with ϕ representing the magnetic flux through the primitive cell. Ever since the work of Hofstadter from 1976 3 , material scientists are well aware of the spectral properties of this operator, which were further investigated with semianalytical methods by Aubry and Andre 4 , and with rigorous methods by many mathematical physicists 5,6 . It is known, for example, that for all irrational values of ϕ, the spectrum is a Cantor set whose size measures to Spec H ð Þ j j¼j 4 À 2λj. Furthermore, for λ > 2, the eigenmodes are localized in space with exponentially decaying profiles, while the eigenmodes are extended for λ > 2. At the transition point λ = 2, the spectrum is fractal and its measure vanishes, hence it consists mostly of gaps. As we shall see, all these gaps are topological in the sense that they all support nontrivial bulk-boundary correspondences, a statement that actually remains valid even when the potential and the couplings are replaced by generic quasiperiodic functions 7 .
Fractal spectra encode interesting patterns of self-similarity, and large efforts have been devoted to their experimental observation. Several groups have proposed indirect methods to chase the butterfly spectrum, for example, by measuring quantization Hall conductance in magnetotransport 8,9 or by using ultra cold atoms as a tool of mapping Haper Hamiltonian [10][11][12] to obtain evidence of fractal spectrum. The first realization of a butterfly spectrum was carried out by measuring microwave as well as acoustic transmittance through an array of scatterers 13,14 . Years later, several groups independently reported evidence of butterfly spectra by investigating the electronic behavior of graphene super-lattices under the modulation of long-range periodic potential and magnetic field perpendicular to graphene superlattices [15][16][17] , and fractal magnetic Bloch states were revealed 18 . A more recent work used nine superconducting qubits as a tool to demonstrate the spectroscopic signature of this energy spectrum 19 .
In Harper's equation, the parameter θ represents the electron quasimomentum in a chosen direction. As such, by varying θ over the full interval [0, 2π], one essentially recovers the physics of electrons on a 2D lattice subjected to a perpendicular magnetic field 20 . These quantum systems display the Integer Quantum Hall Effect (IQHE), hence it is not difficult to understand why the almost Mathieu operator possess a nontrivial bulk-boundary correspondence. Nevertheless, the possibility of implementing the operator with quasiperiodic 1D platforms opens up new ways to observe and study the associated topological phenomena. A number of theoretical and experimental works have recently induced quasiperiodicity in a broad variety of systems 7,21-27 to explore topological phase transitions and edge states. Even more, topological phases in higher dimensions 23,[28][29][30][31] have been predicted and the those characterized by second class Chern number in 2D quasiperiodic crystals observed experimentally 21,22 .
Despite all these advances, reconfigurable experimental platforms were missing to construct the exact quasiperiodic lattice and to emulate the topological effects, yet they are much needed for a thorough exploration of their physics, such as the observation of the Hofstadter butterfly spectrum or pumping the edge states from one boundary to the other. In this paper, we present our design and 3D-printed realization of a reconfigurable 1D acoustic metamaterial, made of resonators with tunable onsite frequency. The reconfigurable design allows us to construct arbitrary phason configurations of the quasiperiodic array, thus, it enables experimentally mapping Hofstadter butterfly spectrum and topological edge states for sound in a straightforward way by simply tuning the piston of resonators. Compared to the previous works which mapped the fractal bulk spectrum in the classical system 13,14 through simulating "onsite potential" by a step function, our structure can mimic the variation of the "onsite potential" in Harper/Aubry-André equation in a "sinusoidal" manner. In addition, both adiabatic pumping (mapping edge states spectrum) and nonadiabatic pumping (mapping Hofstadter butterfly spectrum) can be easily emulated in our reconfigurable acoustic platform, which are very hard to achieve in the solid states of matter.

Results
Reconfigurable acoustic metamaterial design. The array consists of tunable acoustic resonators connected via circular channels, called bridges hereafter. The height of each resonator can be adjusted by rotating a piston inserted inside a thread, which is printed both on the inner side of the resonator and on the surface of the piston itself. The reconfigurable resonators are shown in Fig. 1a, and the geometry of the structure was optimized to observe the acoustic spectrum with high resolution (see Methods). For example, in order to realize high and stable quality factors for the resonant modes (~50), the unperturbed height of the resonators is chosen as h 0 = 40 mm and the modulation depth as δh 0 ¼ 0:12h 0 . The diameter of the connecting bridges affects the size of the fractal gaps in the butterfly spectrum, therefore an effort was put into optimizing the diameter value to maximize the size of these gaps. The details of the structure optimization can be found in Supplementary Note 1.
The first-order acoustic pressure mode of the single resonator oscillating in the axial direction of the cylinder (~4.5 kHz) (shown in Supplementary Fig. 1c), with one node at the center of the cavity, is one of the interests for our experiment. bridges are placed at the bottom of the resonators to minimize the effect on the coupling strength κ of the varying height through the pistons at the top. As long as the modulation depth of height δh is not very large, this assumption holds well, as experimentally confirmed by testing the frequency response of two-coupled resonators (for details see S.M. Section. II). Therefore, if δh h 0 ( 1, and the variation of the resonators' height follows the modulation where n is the site index of the resonators, the recursive equations describing an infinite chain of coupled resonators assume the form where ω 0 ffi πc h 0 is the resonator unperturbed onsite angular frequency, c is the speed of sound, γ ffi δh h 0 ω 0 , and φ n is the amplitude of the collective resonant mode at site n. The parameter θ is a phase factor often referred to as phason. The value of the coupling constant κ is determined by the position and the diameter of the connecting bridges and it can be adjusted by moving the bridges up or down as well as by varying their diameter. The numerical value of κ was experimentally determined from the hybridization of the modes in a dimer configuration (Supplementary Note 2), and the position and geometry of the bridges was adjusted until κ = γ/2, hence tuning the system at the transition λ = 2, where the spectrum is mostly composed of gaps.
Bulk-boundary correspondence. We can encode the data of the collective mode in the vector jφi ¼ P n φ n jni, in which case Eq. (1) reduces to the eigensystem problem for the almost Mathieu operator H ϕ,θ . It is known that, if ϕ is irrational, the spectrum of H ϕ,θ is actually independent of θ, but this is not the case if ϕ is rational. Hence it is more convenient and appropriate to work with the reunion of the spectra Spec ϕ ð Þ ¼ ∪ θ SpecðH ϕ;θ Þ, since, as we shall see, the topological gaps occur in Spec(ϕ). For example, the individual spectra Spec(H ϕ,θ ) can display additional gaps for which it is not possible to define topological invariants, as shown by the midgap in the Fig. 2a. As subsets of the real axis, Spec (ϕ) are known to depend continuously of ϕ and, as such, by sampling ϕ over rational values we can compute an accurate representation of the spectral butterfly, as shown in Fig. 2b. As one can see, in all cases there is a fractal network of spectral gaps. If ϕ is irrational, then between any pair of such gaps there is an infinite number of additional gaps, hence counting the gaps in an orderly fashion is futile. Nevertheless, each gap can be uniquely labeled 32 by the value of the integrated density of states (IDS) inside the gaps, defined as the number of eigenvalues below that gap divided by the length of the system, as the length is taken to infinity. For our context, IDS is known 30 to be quantized as Furthermore and very importantly, there is the constraint 0 IDS 1. The integer m is connected to the first Chern number C, a fact that can be seen from Streda's formula ∂IDS ∂ϕ ¼ C (The details of applying this formula in the aperiodic case are presented in ref. 33 ). Given the constraints on IDS, one can immediately see that C ≠ 0 for all gaps except for two cases, namely for the gaps above and below the spectral butterfly. In other words, every true gap of the butterfly is topological. The Chern numbers are continuous (hence constant) functions of ϕ as long as one navigates inside a spectral gap 34 , and hence does not cross any spectrum of the butterfly. As such, the Chern number can be evaluated for each gap using rational approximants p/q of ϕ (the analytic form of the Chern number for both irrational and rational ϕ is presented in Supplementary Note 3). e.g., by turning Eq. (2) into the Diophantine equation 35,36 where j is the index of the minigap counted from the lowest frequency gap. Alternatively, one can integrate the Berry curvature of the gap projection P G over the magnetic Brillouin zone (detailed in Supplementary Note 3). Both methods lead to the same answers and some of the computed Chern values are shown in Fig. 2b for (p, q) = (1, 6). As was shown by Hatsugai 37 , the edge of a halved IQHE system carries a topological invariant whose numerical value coincides with the bulk Chern number. Furthermore, this edge topological invariant counts the difference between the positively and negatively dispersing edge bands 38 . In our context, this means that, at the edge of a half acoustic quasiperiodic chain, we should observe a number of chiral bands equal to the value of bulk Chern number as θ is varied. This is indeed confirmed by our explicit calculations reported in Fig. 2b. Importantly, the stability of this bulk-boundary correspondence against disorder was established in ref. 30 .
Furthermore, we demonstrate in Fig. 2c the spectacular change in the spatial profile of the resonant modes as the system transitions between the localized and de-localized regimes by tuning κ, for example, changing the diameter of the bridges d c . When κ f 0 < 0:06 (corresponding to the case λ > 2), the eigenstates are localized and exponentially decay, in the example of d c = 8 mm, the acoustic pressure waves in the array are indeed localized (Fig. 2c). However, when κ f 0 > 0:06 (corresponding to the case λ < 2), as shown in Fig. 2e, for d c = 12 mm, the acoustic pressure waves extend at every site along the array. Contrary to what one will expect, based on the fact that topological invariants are usually carried by the extended states, the topological bulkboundary correspondence is not affected by these qualitative differences. The reason is that the topological bulk Chern number, defined above, involves also the virtual dimension of the phason 6 .
Observation of the Hofstadter butterfly spectrum. In order to observe the fractal structure of the butterfly spectrum, the number of degrees of freedom of the system has to be large, namely the number of resonators should be sufficient to resolve the fractal features. However, the resonance peaks of the system have finite bandwidth due to unavoidable radiation and material loss (e.g., introduced by holes drilled to probe the system and the resin absorption). In addition, acoustic viscosity also contributes to the loss since the size scale of the resonator is in the centimeter. Therefore, it is impossible to resolve the higher-order minigap as the corresponding peaks overlap with each other if the bands are too close (with the separation less than~70 Hz). Here we choose a number of cells q = 37, as a compromise between the requirements of high resolution and time needed to take measurement for different configurations of the array, while being aware of the finite Q factor of each resonance that would fundamentally limit the observation of higher-order gaps.
For the first time, we mapped the Hofstadter butterfly through the acoustic density of states (DOS) instead of measuring the acoustic wave transmittance/reflectance for each value of p. We extracted DOS by measuring the frequency responses of each resonator in the quasiperiodic array, in such way the boundary effects on the bulk spectrum can be reduced by excluding DOS of the boundary sites from the spectrum. In the measurement, the speaker is placed at the bottom of each resonator, and the directional microphone is placed at the top of the piston screwed into the same resonator. Both speaker and microphone are placed next to the probe holes, which are intentionally introduced into the geometry as seen in Fig. 1a. The frequency response range from 3600 to 5800 Hz was measured for every resonator, except for the boundary ones. The normalized DOS is obtained through these data (see Methods Section for details). As shown in Fig. 3a, the phason is set as π 3 , the integer p varies from 1 to q = 37, and its variation changes the spatial periodicity of the sinusoidal function in the array. For example, if p = 8, the sinusoidal wave in the array undergoes eight cycles, and the corresponding measured DOS, shown in Fig. 3b, clearly reveals the main gap m 0 in the butterfly spectrum, where the indexing method of Hofstadter's work is adopted 3 . The higher-order fractal gaps c −1 in the c −1 region, and one in the r −1 sub-region of the c −1 region are also clearly observed. These results allow us to conclude that at least three orders of fractal gaps in the spectrum can be resolved in our setup.
Normalized DOS for all p values have been measured, and their logarithmic scale amplitude is presented in Fig. 3c, with the brighter part of the spectrum indicating higher DOS in that region. As expected, the spectrum is overall symmetric in the horizontal direction within the range of tolerance due to experimental errors. However, a slight asymmetry present in the vertical direction can be observed, and it is attributed to the fact that the height variation has a minor effect on the hopping strength κ, which has also been verified by our first-principle calculations (shown by red colored dots in Fig. 3d). To better resolve the minigap in the spectrum, the positions of the resonance peaks are extracted from each site for all p values and depicted in blue colored circles in Fig. 3d, and match the simulation results in a precise way. The legality of taking the open boundary in the bulk spectrum measurement comes from several aspects. First of all, frequency responses at boundary sites were not taken into account to generate the butterfly spectrum, therefore, we do not observe many edge states existed inside the bandgaps in the butterfly spectrum (there are few in the lower right main gap of the spectrum). Second, the bulk states, due to the existence of loss in the acoustic structure, cannot propagate over a large distance when they are excited by the source, thus, the effects of boundary conditions on the exited bulk states are lessen because of the reduced propagation length of bulk states.
Pumping topological edge states. To explore the fractal gaps and edge states associated with them in our acoustic quasicrystal, we fix the magnetic flux to ϕ = 1/6. In order to resolve the fractal edge state spectrum, a large number of cells are required. We choose p = 4 and q = 24 in our experimental setup. The field profiles distributed in the arrays shown in Fig. 4a reveals the topological pumping of boundary mode from one side of the array to the other side, as the phason θ varies in the range indicated by red colored lines in Fig. 4b. The phason θ acts similar to a momentum vector, although we present the variation of the arrays in real space. Although the field profiles of the arrays are depicted schematically, they closely reflect the results from firstprinciple calculations. Only half of the energy spectrum calculated from COMSOL simulations is presented in Fig. 4b due to the spectrum symmetry. Because of the bulk-edge correspondence, one and two gapless edge states are expected in the m 0 and c −1 gaps at each boundary.
Compared to the work 27 , there is no limitation on the number of phason configuration, though we found that taking measurement for 48 phason values is sufficient to clearly resolve the edge spectrum with better resolution than the one in ref. 27 . Followed the same procedure as one in measuring butterfly spectrum, the acoustic DOS for each value was obtained. The reconstructed energy spectrum for all phasons is plotted in Fig. 4c, and it is found consistent with our first-principle calculations. One and two edge states are observed in m 0 and c −1 gaps, respectively, as indicated by dashed lines in Fig. 4c. When compared to the (1) at λ = 2, with ϕ sampled over the rational values p/q, q = 611, p = 0,…, q. c-e Spatial profiles of the resonant modes (picked at 4531 Hz) simulated with COMSOL Multiphysics for the 1D quasiperiodic array. The arrays in c-e have the bridges with different diameters d c , which correspond to the cases of λ > 2, λ = 2 and λ < 2, respectively. For all three cases, the parameters were fixed at ϕ ¼ 8 37 and θ ¼ π 3 . The colorbar shows the normalized magnitude of acoustic pressure theoretical results, the crossing of these edge spectra shifts to higher frequency because of the effect of different boundary conditions, specifically, perfectly matched layer boundary conditions are applied in the simulations, while in reality the boundary cells are open to the ambient air environment, which emulates the open boundaries set in the numerical calculation 39 . However, the edge states emerge irrespective of the boundary conditions, and always retain their gapless and asymmetric character with respect to the phason variation, which vividly demonstrates the robustness of the topological edge states.
The logarithmic scale field profiles of the edge states shown in the Fig. 4d, e confirm their localization at the boundaries both in simulations and in the experimental results, respectively. The edge state profile lines in Fig. 4e are color coded to agree with colors in Fig. 4c. In addition, as can be seen from Fig. 4e, the edge states present in the m 0 gap (yellow colored) are found to be more localized than those present in the c −1 gap (blue colored), although the edge states in both gaps have the same distance to the nearest (upper) bulk bands, which excludes the effect of decay length due to bulk-edge frequency difference. Interestingly, a power-law behavior is discovered in the decay profile of the edge states, which reveals the modal property of the quasiperiodic crystals, as predicted theoretically in ref. 40 .

Discussion
In this work, we experimentally mapped the Hofstadter butterfly and fractal edge spectra using a 3D-printed reconfigurable acoustic quasicrystal, and first time used DOS instead of the transmittance/reflectance of the acoustic waves to construct the butterfly spectrum. Our approach can be easily extended to other systems, such as LC circuits, or arrays of silicon ring resonators 41 , which may possess higher-quality factors, thus offering the ability of mapping richer features in the spectrum. For example, if the parameter ϕ is an irrational number, or the integer p (coprime with q) is very large, the very high order fractal gaps and their associated edge states can be observed in the spectrum for the case of high-quality factor resonators [see Supplementary Fig. 4]. Thus, the larger the quality factor of the modes, the finer fractal features may be observed. In addition, the idea of reconfigurable design can be generalized to systems emulating higher dimensions, for example by mapping the effects of second or third Chern number in 4-and 6-dimensional spaces, implemented in 2-and 3-dimensional quasicrystals, respectively. Such systems may enable the exploration of even more subtle phenomena, which may be unfeasible in periodic topological systems. Our work opens up the possibility of realizing topological devices with robust performance stemming from the fundamental physics of quasiperiodic topological systems.

Methods
Experimental details. The reconfigurable quasiperiodic array in Fig. 1b  B9Creator v1.2 3D printer was used to print the reconfigurable model. All resonators and pistons were made with acrylic-based light-activated resin, a type of plastic that hardens when exposed to UV light. The shell of the resonator was printed with a sufficient thickness (2-3 mm) to ensure a hard wall boundary condition. Two connected cells were printed at a time and the cells were designed in a way to interlock tightly with each other.
The frequency generator and Fast Fourier Transform (FFT) spectrum analyzer scripted in LabVIEW were used in the data processing. For the details of measurement, the speaker was placed at the bottom port and the microphone at the top port of the same site. Both the speaker and microphone were closely touched with the ports to achieve the maximum coupling between source and the system, as well as maximum amplitude of signal. The frequency generator was used to run a sweep from 3600 to 5800 Hz in 20 Hz intervals and with the dwell time of 0.5 s which is enough for the FFT spectrum analyzer to obtain the stable amplitude responses φ(f) at each frequency. For the calculation of normalized DOS, field distributions φ(i, f) are obtained by repeating this process for each site i 2 ½2; q À 1. The data for each site was normalized on the total volume of signal summed over frequencies and sites, as well as on the free space amplitude response between the microphone and the speaker, ð Þ=φ air f ð Þ. The signal spectrum for an array of q−2 resonators S n f ð Þ ¼ To observe the field distribution excited by the speaker, the speaker was placed at the port of the site of interest, and the microphone over each site of the array to measure the magnitude at the desired frequency. In all acoustic measurements, the noise floor in the signal is less than −120 dB, which is much less than the signal level in the desired frequency interval.
Numerical details. For the first-principle calculations, the finite element solver Multiphysics COMSOL 5.2a and the Acoustic module were used to perform fullwave simulation. In the acoustic propagation wave equation, the speed of sound was set as c = 343.2 m/s, and density of air as ρ ¼ 1:225 kg=m 3 . The dimensions of the structure are the same as the fabricated ones. For eigenvalue calculations, the continuity conditions were imposed along the ends of the quasiperiodic array. For the calculations of edge states spectrum, perfect matched layer boundary conditions were applied on the boundaries of the array. Coupled model equations were used to fit the physical parameters of the cells, calculate the Hofstadter butterfly spectrum and topological edge spectrum. For the two coupled resonators, the fitted onsite frequency is f 0 ∼ 4502 Hz, coupling strength is κ ∼ 270 Hz. Mapping the edge spectrum using adiabatic pumping of phason θ for the case (p, q) = (4, 24). a Schematics of edge state transfer during phason θ pumping; the range of this pumping is indicated by the orange line in Fig. 4b. b Band structure obtained from first-principle calculations. c Band spectrum mapped from the measured normalized density of states (DOS) for each θ, the colorbar shows the magnitude of DOS in arbitrary unit. d Acoustic pressure distributions of edge states found from first-principle simulations; these field profiles, from the bottom panel to the top panel, correspond to the edge states marked by dots in Fig. 4b and are arranged from the lower frequency to the higher frequency, the colorbar represents the magnitude of acoustic pressure. e Measured acoustic pressure amplitude distributions S n in logarithmic scale at different sites of the array from 1 to 24, these field profiles correspond to the edge states marked in dots in Fig. 4c