Acoustic graphene network loaded with Helmholtz resonators: a first-principle modeling, Dirac cones, edge and interface waves

In this work, we study the propagation of sound waves in a honeycomb waveguide network loaded with Helmholtz resonators (HRs). By using a plane wave approximation in each waveguide we obtain a first-principle modeling of the network, which is an exact mapping to the graphene tight-binding Hamiltonian. We show that additional Dirac points appear in the band diagram when HRs are introduced at the network nodes. It allows to break the inversion (sub-lattice) symmetry by tuning the resonators, leading to the appearence of edge modes that reflect the configuration of the zigzag boundaries. Besides, the dimerization of the resonators also permits the formation of interface modes located in the band gap, and these modes are found to be robust against symmetry preserving defects. Our results and the proposed networks reveal the additional degree of freedom bestowed by the local resonance in tuning the properties of not only acoustical graphene-like structures but also of more complex systems.


Introduction
Over the last years, in the context of metamaterials, a plethora of sophisticated acoustic structures exhibiting unusual wave properties have been theoretically proposed and also experimentally studied. Examples include Dirac cones [1], unidirectional propagation [2][3][4], topological interface waves [5][6][7], robust localized modes [8], etc. Recently periodic acoustic/elastic media have been proven to be an excellent platform for studying various wave phenomena [9]. In particular, artificial structures exhibiting Dirac degeneracies have been widely studied due to the direct link between the conical dispersion and the topological properties of wave propagation [10][11][12].
Owing to the experimental advantages compared to other classical wave systems (such as easy implementations in macro scale [13,14] and direct measurements [15,16]), nowadays, acoustics is intensively used as a test-bed for topological wave physics.
The usual route to study these phenomena in continuous classical systems is to achieve particular designs which are analogs of discrete models with topological characteristics [17,18]. That way, it is expected that the proposed continuous wave systems will inherit the topological properties of these discrete models, especially regarding the robust transfer of protected states [19][20][21]. One of the first and most thoroughly studied models in the field of condensed matter is the honeycomb discrete tight-binding model, which was originally used to describe the electronic properties of graphene.
The honeycomb tight-binding model with a broken inversion (sub-lattice) symmetry of the unit cell is known to possess a quantum valley Hall (QVH) topological phase transition [22,23], and the corresponding valley-protected edge states have been shown to persist under particular types of imperfections [24][25][26][27]. To observe these phenomena in acoustics, different theoretical proposals and experimental demonstrations have Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
In most cases, a dispersion relation which mimics the one of graphene is sought by using either connected acoustic cavities [30] or sonic crystals [31] with a honeycomb geometry of the scatterers. However, to the best of our knowledge in all of these cases the corresponding discrete model is obtained by means of the k·p method (and thus it works well only close to a high symmetry point) [32], or by a tight-binding model whose coupling (hopping) coefficients are practically fitting parameters [33,34].
In this work, we study an acoustic honeycomb network composed by connected waveguides whose Bloch wave solutions satisfy an eigenvalue problem which is exactly the same as the Bloch Hamiltonian of graphene. Taking advantage of the simplicity of the structure and of the modeling, we study the effect of Helmholtz resonators (HRs) located at the network nodes. Our goal is to propose an alternative way that can break the inversion symmetry of the unit cell and create a gap around the Dirac point to exploit the corresponding valleyprotected edge and interface modes of the network. Our structure allows us to independently tune both the resonance frequency of HR and the Dirac point frequency. We use this flexibility of our design by setting these two frequencies at the same point and we study the combined effect of the resonance and of the symmetry breaking. Both edge and interface modes are identified and the influence of boundary conditions are investigated. We found that robust propagation against a corner of 120°can still be observed as in the case of valley-protected edges/interface modes.
The article is organized as follows: in section 2, the honeycomb acoustic network and its modeling are presented. The exact mapping of the acoustic network Hamiltonian to the one of graphene is obtained. Based on acoustic graphene network, the effect of HRs loaded to the acoustic network is studied. Section 3 shows the derivation and discussion of zigzag edge waves considering an inversion symmetry breaking network due to resonance stemming simply from HRs loaded to one of the graphene sublattice nodes. The existence of edge waves in two different edge configurations, that is, with/without loaded HR, is studied. In section 4, we further study the sound waves propagation on the interfaces by combining two symmetry breaking networks in section 3. The existence of interface modes, and their robustness against defect (a corner of 120°) is shown even in the appearance of the combined effect of the resonance and the symmetry breaking. Finally, conclusions are presented in section 5.

Model
A schematic presentation of the graphene network is shown in figure 1(a), where a unit cell containing two junctions (marked as A and B) is indicated by a red-dashed box. This network can be achieved by using connected acoustic waveguides. A three-dimensional (3D) view of the unit cell is shown in figure 1(b), where the distance between two consecutive junctions is L and the radius of the acoustic waveguides is R t . By placing the unit cell of figure 1(b) in a hexagonal lattice, the acoustic honeycomb network is constructed. The sound pressure at each junction is linked to the pressures of its neighboring junctions through the air channels. As we will show later, sound wave propagation in this network can be exactly mapped to the electron behavior in graphene under the tight-binding approximation.

Exact acoustic analog of graphene
To obtain the governing equations and the dispersion properties for the acoustic network, we adopt the methodology initially proposed in 2D lattices of tubes [35] and further developed in 2D acoustic channels of lattice structures [36,37]. We assume that sound wave propagation between junctions is well described as monomode propagation (plane wave approximation) as long as the radius of the air channels is much smaller than the distance between the two junctions A and B, i.e. R t = L and the frequencies of interest are much lower than the first cutoff frequency of the waveguides. Each unit cell is labeled using the normalized coordinates m and n (both integers) as shown in figure 1(a). Then, we employ the continuity of flux at each of the two junctions A, B and we obtain the following system of discrete equations for the sound pressure at site (m, n): In equations (1), a p m n , with α=A, B is the pressure at each junction, and ε α is the 'energy' term which is given by e e e = = .
, q y =3L k y /2 with k x , k y the wave vectors along the xand y-directions in the first Brillouin zone (BZ) as shown in figure 1(c). By substituting the wave solution into equations (1), the following governing equation can be derived, . By solving the eigenvalue problem described by equation (4) where ε g are the eigenvalues, we obtain the dispersion relation e = = + + F q q q q q , 3 4 cos cos 2 cos 2 , 5 which gives the mapping between Bloch wave number and frequency. It is very interesting to note that equation (4) is identical to the Hamiltonian of electrons in graphene under the first neighbor interactions, manifesting the fact that the honeycomb network is an exact acoustic analog of graphene. However, in contrast to many of the existent classical analogs of graphene where an effective tight-binding model is used to describe the properties of the system [30], the corresponding discrete model in equation (4) directly stems from the acoustic conservation laws.
To show more clearly that the honeycomb network is an exact acoustic analogue of graphene, we consider the network of L=0.4 m, R t =0.025 m, see figure 1(b). The dispersion relation given by equation (5) along the irreducible BZ is shown by solid blue lines in figure 1(c) as a function of ε g . In comparison, we also show the corresponding dispersion curves in dots obtained by finite element simulation and a good agreement with the analytic results is exhibited. The small discrepancy is due to three-dimensional nature of the junction with near field effects that increase with increasing R t /L. In the graphene dispersion curves, the Dirac points can be found at ε g =0, leading to the corresponding Dirac frequency f K =c/(4L)=215 Hz for the graphene network. Note also that the right-hand side of equation (5) does not depend on the geometrical properties of the particular system (e.g. length of the air channel). As such the discrete model is valid for any acoustic honeycomb network as long as monomode propagation in the air channels is justified.

Graphene network loaded with HRs
We are now interested in the effect of local resonances by side-loading HRs at the two junctions A and B of the graphene network as shown in figure 2(a). The sketch of the resonator is also shown in figure 2(a). At sufficiently low frequencies, those HRs are regarded as point scatterers, whose presence can be well described by modifying the flux conservation at each junction to include the flux entering the resonator. Then equations (1) remain unchanged, but the energy term ε α yields The function G α ( f, f α ) incorporates the geometrical characteristics of the HR and depends on both the sound wave frequency f and its individual resonance frequency f α (see appendix A). Considering the modifications on the energy term, the corresponding dispersion relation with the presence of resonators takes the following form e e = + + q q q 3 4 cos cos 2 cos 2 . 7 As an example, we consider the case when the two resonators are identical with a resonance frequency f A =f B =320 Hz. The dispersion curves are shown in figure 2(b) by blue solid lines. For comparison, the dispersion curves of the graphene network without resonators are also shown in figure 2(b) by red dashed lines.
As it can be seen, by adding HRs to the system, we have generated two more propagating branches and a band gap around the resonance frequency depicted by the green shaded area in figure 2. Note that, since the resonators are identical, the unit cell still preserves the lattice symmetry, hence the existence of the Dirac point at the corner of the BZ (K point). However, the appearance of the two extra propagating branches due to resonances leads to the presence of an additional Dirac point at the K point. The Dirac frequencies are determined by equation (7), which satisfies ε A ε B =0 at the K point. By combining with equation (6), it turns out that the Dirac frequencies depend also on the resonance frequency of the HRs (see appendix B). This can be seen in the dispersion curves in figure 2(c) where the resonance frequency is tuned to be f A =f B =f K . It shows that the resonance band gap has been moved around the resonance frequency f K while the two Dirac points are also shifted compared to figure 2(b). When the two resonators are different as shown in figure 2(d), the two Dirac points at the K point are lifted due to the breaking of the inversion symmetry of the unit cell. As a consequence, two full gaps marked by orange shaded areas appear. This is different from the gap depicted by green shaded area which stems from the resonances. It should be mentioned that there are different ways to break the inversion symmetry in order to open a gap at the K point of graphene/graphene-like structures. For example, in graphene/graphene-like systems, the inversion symmetry can be broken by introducing difference in A, B sublattices, such as, bianisotropic responses [23], refractive index [24], cylinder sizes [25], masses [26], bilayer designs [28], acoustic cavities [30], etc. However, in this paper we are interested in the case of breaking Dirac points by loading resonators to the junctions, which leads to interesting edge/interface waves in the corresponding gaps when boundaries are considered. In the next section, we will investigate these edge/interface waves which have not been widely reported in acoustic systems so far. = The red dashed lines in (b) correspond to the dispersion curves of the unloaded graphene network. The green areas represent the band gaps due to resonance, and the orange areas correspond to the gaps due to inversion symmetry breaking.

Zigzag edge waves
From hereon, we study the graphene network loaded with only one HR at junction A of the unit cell (single loaded graphene network (SLGN)), and the resonance frequency of the resonator is tuned to be » f f A K . The corresponding dispersion curves are shown in figure 3(a). The blue solid lines represent the analytical results from equation (7) for ε B =ε g (no resonator), and the black dot lines the finite element simulation. As it can be seen, the degeneracy of the Dirac point at the K point is lifted creating a band-gap, due to inversion symmetry breaking, similar to other inversion symmetry breaking systems [5,7,23,31]. On the other hand, a unique feature of our system is the appearance of an additional nearly flat band inside the generated band gap around the resonance frequency f A . It is interesting to note that this band becomes flatter as f A approaches f K . In particular, at the resonance frequency ( f A ) the pressure field at junction A (with resonator) should be equal to zero. Since we have chosen f A to be half of the Bragg frequency (Dirac point), the wavelength is 4L imposing a maximum pressure at junction B (no resonator). Thus, a localized solution (described in the appendix C), where pressure at only a single junction B of the whole lattice is non-zero and every other junction (A and B) is zero, is possible and induces the flat band.
To investigate the presence of edge waves when boundaries are considered, we consider a SLGN which is finite in the y-direction but infinite in the x-direction. A supercell of the SLGN containing N unit cells (inset of figure 3(a)) is depicted in one of the red-dashed boxes of figure 3(b). The two zigzag edges, e.g. top and bottom ends, are set to be zero pressure boundaries. For example, considering the supercell on site j (contains both 2m and 2m+1), the zigzag boundary conditions in the network are expressed as By combining with equation (1), the following equation connecting two consecutive supercells can be derived [ ] contains all the sound pressures at the junctions inside the supercell. D( f ), S( f ) are 2(N−1)×2(N−1) matrices depending only on the sound wave frequency, which can be obtained by summing up all the equations of motion of the N unit cells combining with the boundary conditions in equations (8). By applying an ansatz of Bloch wave solution, = v V e j q j 2i x , equation (9) leads to the following generalized eigenvalue problem By sweeping the frequency f, the corresponding Bloch wave vector q x can be obtained by equation (10), resulting in the dispersion curves for the supercell. It should be mentioned that the common approach to get dispersion curves by sweeping Bloch wave vector in the BZ is not suitable in our system when HRs are loaded to junction A of the network. The challenge is that the energy term ε A (equation (5)) contains the contribution part of the resonator G A , while ε B is not affected by resonance, leading to the breakdown of the eigenvalue problem of the form shown in equation (5). However, the eigenvalue equation (10) avoids the obstacle caused by resonance since its eigenvalue gives rise to the Bloch wave vector. In addition, from the eigenvalue problem in equation (10), not only the propagating modes, but also the evanescent modes can be obtained. This might lead to great advantages for the study of wave systems where dissipation is considered. The dispersion of a supercell with N=17 is shown in figure 3(c). We note here that the top edge is free of resonators, we call it free zigzag edge (FZE), while the bottom edge is loaded with resonators, we call it resonant zigzag edge (RZE). The dot lines correspond to the simulation results. As can be seen, there are two edge wave branches illustrated by blue dot lines around 185 and 270 Hz. These two branches support the edge waves on the RZE, the bottom edge of figure 3(b). This can be confirmed by the eigenmodes on the two branches. For example, we investigate the modes of the blue branches using a finite element method, the corresponding mode profiles are presented in figure 3(e). It shows that sound waves are localized on the RZE and are evanescent along the bulk.
On the other hand, edge waves exist also on FZE (on the top edge of figure 3(b)). These edges waves are located into the nearly flat band as it is shown by the close view of figure 3(d) where the edge branch of the FZE is indicated by a magenta dot curve. The eigenmode profile of the edge mode is also depicted in figure 3(e), where sound field is seen to be localized on the FZE and decaying into the bulk.

Zigzag interface waves
We now consider sound waves propagating along the interface constructed by two SLGNs. The study of the interface is interesting as it can lead to potential investigations of topological modes associated with the QVH effect in the current graphene network loaded with HRs. Regarding the geometry of the SLGN, there are two different kinds of interfaces: (1) Resonant zigzag interface with two consecutive junctions side-loaded by HR as shown in figure 4(c). (2) Free zigzag interface with two consecutive junctions with no resonators as shown in figure 5(c). The dispersion of interface modes is calculated in a similar way as discussed in section 3.
Considering the resonant zigzag interface, figure 4(a) displays the dispersion curves of the supercell, which contains N=17 unit cells with the resonant zigzag interface in the middle. As can be seen, the dispersion exhibits similar bulk modes features (gray dots), while six total branches of edge/interface modes are found, instead of three compared to figure 3(c). The four branches marked by red dots in figure 4(a) correspond to interface modes, and the other two branches close to the nearly flat region support edge modes on the two FZEs as shown in figure 4(b) by magenta dots. Compared to the case without interface ( figure 3(c)), the two RZE branches are absent as this RZE boundary is eliminated from the supercell of figure 4(c). In addition, one additional branch for edge waves on the FZE appears on the top of the near flat region (see figure 4(b)) due to the symmetry of FZEs on both ends of the supercell. Figure 4(c) shows three eigenmodes corresponding to three lower branches marked as red in figure 4(a) and magenta in figure 4(b). The modes in the red bands show the localization of sound pressure on the interface (middle of the supercell), exhibiting the feature of interface modes. The mode in the magenta band in figure 4(b) has a similar pressure field distribution as the FZE mode in figure 3(d) as they are both free zigzag modes on the open ends.
Regarding the free zigzag interface, figure 5(c), the band diagram is depicted in figure 5(a), where two edge wave branches on the RZEs (blue dots) and two interface branches (red dots) are observed in the band gaps. It should be noted that, unlike the results in figure 3(c), each of the RZE branches in figure 5(a) is a degeneracy of two branches due to the symmetry of the RZEs on both ends of the supercell. The two additional branches (red dots) correspond to the interface modes compared to figure 3(c). In figure 5(b) showing a close view around the flat band, we observe the absence of edge modes as expected since there is only one kind of boundary corresponding to RZE.
We also pick three eigenmodes corresponding to three different branches to show the pressure field distribution of each mode in figure 5(c). The modes of the interface branches (red bands in figure 5(a)) show a similar distribution as the red ones in figure 4(a) since pressures are localized on the interface. The mode in the RZE branch exhibits the localization of pressures on the RZE, which is consistent with the results of edge waves on the RZE as observed in figure 3.
Our results show that in the presence of HRs (resonance frequency equal to f K ) in the SLGN, there are 3 edge modes branches in the absence of interface and 6 edge/interface modes branches when interfaces are present. It has been reported in [11] that there are 3 branches of edges/interface modes due to the inversion symmetry breaking without resonance. However, our study of SLGN shows that the presence of HRs with resonance frequency close to f K injects one additional nearly flat band to the middle of the band gap, splitting the band gap into two. As a consequence, more branches for edge/interface modes appear in the gaps.
As it has been studied in the analogs of QVH system [5,7,23,31], the interface modes in the full gap appearing by breaking the Dirac cone at the K point are robust in a sense that they can travel through corners of 120 • with respect to its original direction of propagation. This phenomenon can also be observed in our SLGN. To verify the above analysis, we conduct simulations on a finite size SLGN of 34×27 unit cells. We create a FZE type of interface forming a 'Z' shape in the finite network. We preform numerical simulations of the finite structure in frequency domain assuming an incoming monochamatic wave at the beginning of the air channel  denoted by yellow arrow in figure 6. The other air channels at the edges are set to be zero pressure corresponding to the open boundary in figure 3(b). The results are shown in figure 6(a) for 208 Hz corresponding to the lower edge branch and 6(b) for 240 Hz to the upper edge branch. It shows that in both cases, sound waves travel along the interface path and do not propagate into the bulk or along the edges of the structure. More interestingly, the pressure field distributions of the interface modes indicate a robustness of sound waves propagation again the two 120 • corners along the interface path. Note that the presence of visco-thermal losses is found to change only quantitatively the propagation of the Z-shape interface mode. Disregarding the nearly flat band region in the band gap, the result suggests that the robustness of interface waves might stem from the topological property of the systems similar to the QVH systems, which can be further studied in a future work by taking into account the influence of the nearly flat band.

Conclusions
We have developed a discrete model to describe wave propagation in an acoustic graphene network including subwavelength resonators. Regarding the influence of resonances, we show that they are responsible for the appearance of additional branches and band gaps, and the tunable resonance frequency provides an extra degree of freedom to alter the symmetry of the system. For a strip geometry with zigzag boundaries, the breaking of Dirac cones at the K point due to resonance results in the existence of edge waves. We also show the appearance of interface modes that exhibit a robust propagation against the 120 • corners, similarly to QVH effect. This work is a preliminary step towards more complex topological systems with on-site resonances.

( )
Considering the total flux at junction B, the conservation of flux dictates that u 1 +u 2 +u 3 +u R =0. Replacing velocities with the corresponding pressures, and reproducing the process at junction A, we can obtain the governing equations (1). The corresponding energy terms ε A,B are given by equation (6). And the function G α ( f, f α ) is given by The entrance admittance Y α depends on the geometrical characteristics of the resonators (α=A, B) and is given by

Appendix B. Effects of the Helmholtz resonance on the graphene dispersion relation
According to equation(5) and the expression of F(q x , q y ), passbands are found when ò α or ò g are bounded within [−3 : 3]. In the absence of HRs, the 'energy' is defined by p =  fL c 3 cos 2 g ( )and thus there is no bandgap as it has to be for standard graphene. This function is plotted in figure A2 with the red dashed line. Note that when ò g is crossing the zero axis, it corresponds to a Dirac point. When the system is side loaded by identical HRs, due to the resonant shape of the entrance admittance of the HRs, the 'energy' ò α (see equqtion (6)) now crosses the zero axis at two different points (see figure A2 with the solid blue line) yielding two Dirac points in the considered frequency range.