Theory of holey twistsonic media

Rotating two overlapping lattices relative to each other produces the well known moiré interference patterns and has surprisingly led to strongly correlated superconductivity in twisted bilayer graphene. This seminal effect that is associated with electrons occupying ﬂ at dispersion bands has stimulated a surge of activities in classical wave physics such as acoustics to explore equivalent scenarios. Here, we mimic twisted bilayer physics by employing a rigorous sound wave expansion technique to conduct band engineering in holey bilayer plates, i.e

T he moiré effect is possibly best known from the shimmering effect seen when two geometrically regular arrangements are superimposed as in the case of overlaid textiles, fabrics, parallel lines, etc. These rippled appearances seem particularly impressive at acute angles, at which the overlaid and seemingly competing fabrics stand out at large intervals, i.e., waves with noticeable moiré wavelength. In an apparent unrelated area comprising atomically thin materials such as graphene, silicene, borophene, etc., or their assembled van der Waals heterostructures, a twist among them has shown to give rise to unprecedented electronic properties. Specifically, rotating two overlapping graphene monolayers at a so-called magic angle, gives rise to the collapse of their Dirac cones into flat dispersion bands, which is associated to an array of intriguing properties, such as 2D magnetism, Mott-insulating phases and unconventional superconductivity [1][2][3][4][5][6] .
Generally speaking, lately, artificial sonic and phononic lattices have been used widely to explore with acoustic and elastic waves exotic topological phases and properties, which are hard and sometimes even out of reach to demonstrate in electronic settings. An attractive motivation compared to their electronic counterparts is their easy fabrication and tunability, which often reveal novel and unexpected effects to lead to entirely differing analogous connection to the original physical context. Chern insulators, valley-Hall phases and higher-order topological insulators are a few of many arenas that have been conquered with classical wave acoustics [7][8][9][10][11][12][13][14][15][16] . Among the latest efforts, several groups have already demonstrated that structured and twisted bilayer plates indeed host striking sonic, vibrating and photonic similarities compared to twisted bilayer graphene physics [17][18][19][20][21][22] .
In this paper, we employ a theoretical approach to study how the twist angle, the plate separation and thickness are the responsible actors to enable moiré acoustic dispersion engineering in twistsonic media. The so-called mode-matching technique (MMT) enables us to conduct an entirely semi-numerical study of acoustic waves interacting with this complex twisted structure, on a much faster and storage-efficient basis compared to commercially available finite-element-method solvers like COMSOL. Our findings showcase how the band flatness, the spectral location of the bandgaps and the moiré interference patterns can be studied efficiently, thus providing useful insight into this contemporary branch of physics.

Results and discussion
Theory. Here we present a numerical technique that allows one to obtain the band structure of the twistsonic medium that consists of two rigid plates of thickness h 1 and h 2 that are separated by a thin gap of thickness h g , as shown in Fig. 1. Both plates are perforated by round holes that are periodically distributed in a unit cell of size a 0 . The unit cell contains N holes of radius R 0α located at the positions R 1 α and R 2 α where α = 1, 2, . . . , N. We separate the twistsonic medium into five different regions: region I refers to the first plate onto which sound waves impinge, region II refers to the field inside the through-holes of the first plate, region III refers to the gap separating the plates, region IV refers to the field inside the holes of the second plate and in region V sound is able to emerge into free-space as shown in Fig. 1b.
This separation of regions is the backbone of our theoretical formalism that is based on the MMT, which comprises a modal expansion of the involved pressure and velocity fields in these zones 23,24 . We assume that the surrounding medium is air with mass density and speed of sound ρ air = 1.225 kg/m 3 and c air = 343 m/s, respectively, and that the plates are made of steel or brass, in which the perfect rigid body approximation is very accurate. Thus, if an incident pressure wave P in of unitary amplitude impinges the holey twistsonic medium, the wave will partially be reflected through diffracted sound and partially be funneled through the holey twistsonic medium. The pressure in regions I, III, and V is expanded in terms of plane waves whereas inside the α th subwavelength hole (regions II and IV), only the fundamental cavity eigenmode is excited. The pressure fields in those regions read where G is the reciprocal lattice vector, K is the parallel wave vector, k 0 = ω/c air is the free-space wavenumber, φ a = h 1 , φ b = h 1 + h g and φ c = h 1 + h g + h 2 are out-of-plane phase contributions and q G ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi k 2 0 À jK þ Gj 2 q . The radiation and hole impedances are written as Z G = k 0 /q G and Z = k 0 /q 0 , respectively, whereas r G (t G ) is the reflection (transmission) coefficient and A α , B α are the modal cavity expansion coefficients. Likewise, continuity across the boundaries applies to the z-component, i.e., the normal component of the acoustic velocities that are As detailed in the Methods section, having all regional fields defined, we are now able to impose continuity of the respective fields across the respective interfaces. Moreover, the MMT is applied by projecting diffracted Bloch modes over the normal velocity fields and the fundamental cavity mode over the pressure field. This process introduces an in-plane phase contribution that defines the intra-layer hole-to-hole interaction within the moiré supercell. Additionally, at each interface site, we define a set of Fig. 1 The twistsonic medium. a The rigid bilayer is perforated with N round holes of radius R 0 , which are arranged in a honeycomb lattice comprising a primitive unit cells of size a 0 and the moiré super cell of size a. b The side view of the structure enables one to distinguish the individual regions used for the MMT. modal velocities: The terms introduced in the system, all have a physical interpretation. I α is the irradiation term over the fundamental cavity mode of the α th hole of the first plate, H 1ð2Þ α is the overlap function, ϵ 1(2) is related to the acoustic bouncing back and forth inside the holes, G V 1ð2Þ refers to the coupling of the two sides of each hole, G 1ð2Þ αβ measures the coupling between the holes of each plate, and finally, ψ 1ð2Þ αβ and Φ 1ð2Þ αβ account for coupling across the plates. Upon solving the unknown modal coefficients of Eq. (4), both the cavity expansion and the scattering coefficients can be determined. The latter can be specifically expressed as Moreover, in the absence of incoming sound, i.e., I α = 0, computing the zero determinant of Eq. (4) as a function of both frequency and the parallel momentum, allows us to determine the dispersion relation in dependence of the geometrical parameters. More importantly though, also in dependence of the rotation angle of the holey twistsonic medium.
Numerical results. In numerically treating the holey twistsonic medium, the upper layer is considered perforated by N holes that are placed at sites R 1 α ¼ n 1 a 1 þ n 2 a 2 , forming a honeycomb lattice of period a 0 and lattice vectors a 1;2 ¼ a 0 ½ ± cosðπ=3Þ; sinðπ=3Þ. Beyond this, the N holes of the bottom layer are placed at R 2 α ¼ m 1 b 1 þ m 2 b 2 where the vectors b 1,2 are the rotated counterparts to a 1,2 to account for the twisting. All holes in the study have circular cross section of radius R 0α = 0.25a 0 . Note that the moiré pattern in twisted bilayer is not always commensurate with the original lattice period. Among all possible twisting angles we can work with, the commensurate angles θ ¼ arcos½ð3m 2 þ 3m þ 1=2Þ=ð3m 2 þ 3m þ 1Þ, for some integers m, are the ones defining the size of the resulting superlattice and its underlying sonic moiré interference pattern. Hence, the twisting induced moiré lattice constant becomes In what follows, we show how the size of the supercell alters the dispersive nature of the holey twisted structure, however, additionally, we discuss how the gap separation and the plate thickness change the band flatness and its spectral location, respectively. In this study, we concentrate on twistsonic media comprising three twist angles, which means we have to construct three different superlattices. With respect to the above given formula, we fix those three configurations, as summarized in Table 1. Figure 2 displays multiple computed dispersion relations of twistsonic configurations with constant hole radius and plate thickness as captured in the figure. The three rows indicate computations at constant unit cell size, whereas the rows represent a successive shrinking of the gap separation, i.e., h g = 0.8a 0 , to h g = 0.5a 0 . The solutions of the system correspond to vanishing determinants, whereas the dispersionless soundlines, i.e., the lattice singularities are accompanied with large determinant values that are not solutions to the problem. Reducing the twist angle, i.e., increasing the superlattice size leads to an obvious change of the moiré interference pattern that is indicated by strong modifications of the dispersionless sound lines seen throughout all examples in Fig. 2. Moreover, at the K-points in the band diagrams, we predict distinctive Dirac cones whose respective touching points reside at identical frequencies (red circles). A more distinctive feature of the computations is the compression of these Dirac cones when the two involved plates are pushed together. E.g., for m = 2 at around ka 0 =2π ¼ 0:215 we see that the flatness of the dispersion bands of the cones appear more pronounced when h g is steadily decreased. This finding shows that not only twisting towards the magic angle can lead to band flatness, but a more accessible approach, namely to tune the plate separation, can lead to an equivalent effect. In the following subsection, we further compare the MMT with the finiteelement method and find that the MMT displays superior advantages, both in terms of the computational speed and the storage size.
In the same spirit with the foregoing study in Fig. 2, we now fix the gap separation as indicated in the caption of Fig. 3, but next to varying the twist angle (m = 1-3) we also increase the plate thicknesses, h = h 1 = h 2 . What immediately stands out in all computed band diagrams for the chosen twist angles, is that an increase of the thickness of both plates leads to a redshift (towards lower frequencies) of the dispersion branches. We In order to visualize the acoustic moiré interference fringes, we compute the pressure fields at spectral regions marked by red circles in Fig. 2. We chose to inspect the insonified side of the holey twistsonic media, i.e., region I in Fig. 1, onto which sound waves impinge, thus the first pressure field in Eq. (1) is to be considered in connection with its complex reflection coefficient in Eq. (6). In doing so, we first solve the equations in Eq. (4) to obtain the modal fields. Afterwards, we calculate the reflection coefficient r G in Eq. (6), upon which we introduce it in the expression of P I in Eq. (1) for the chosen plane. At the nearest vicinity of the upper holey plate (z = −0.001a 0 ) of the three structures under study, the backreflected sound unequivocally displays the fingerprints of the periodic moiré fringes as the computations in Fig. 4 showcase. We have taken a square-shaped domain (10a 0 × 10a 0 ) in each case, to display the scattered moiré features for comparison. A remarkable attribute in using holey plates as artificial twisted structures can be discerned upon close inspection of the field plots, which is the said complex interplay between twisted induced moiré pattern and acoustic energy concentrated within the holes.
Comparing the MMT with the finite-element method. In this section, we discuss the advantages of the MMT after comparison with the widely used commercial finite-element-method solver COMSOL.
(i) The computational speed of the MMT is faster. We demonstrate that one of the advantages of the MMT is the computational speed. As shown in Fig. 5, the band diagrams of identical structures (see caption) are calculated with the MMT (implemented using the Julia Programming Language) and by using the commercial finite-element-method (FEM) simulation software COMSOL Multiphysics. The obtained bands show good agreement except for a negligible frequency shift. To compare the computational speed, along the wave vector axis both methods have identical intervals, i.e., 50 steps. Since the number of frequencies in COMSOL depends on the available number of eigensolution, in the MMT we chose 500 points. We ran the MMT and the FEM COMSOL computations one by one using the same computer, which is a Dell Precision 7920 workstation with a Intel(R) Xeon(R) Silver 4210 CPU and 640 GB of RAM. The computational time using COMSOL is 8077s, whereas that of the MMT is only 407s, which is almost 20 times faster.
(ii) The storage size of the MMT code is smaller. The storage size of the FEM COMSOL file is about 1.27 × 10 7 kB, while the storage size of the MMT code is only 416 kB, which is about 3.28 × 10 5 times smaller.

Conclusions
Holey twistsonic media, i.e., the acoustic counterpart of twisted bilayer graphene have been studied semi-numerically thanks to a modal expansion approach. In contrast to widely used commercial finite element solvers, our technique produces relatively fast dispersion relations results without sacrificing storage space. Holey plates are highly flexible and tunable structures for metamaterials, topological and twistsonic applications. Hence, we believe that our tool should serve as solid basis to conduct experimental studies along this line. We showed how these acoustic twisted media not only host these flat moiré dispersion   bands, but moreover, that they can be straightforwardly altered through the involved geometry.

Methods
The semi-analytical method used to compute the dispersion relations as well as the pressure fields is the Mode Matching Technique, whose details are derived in this section.
Let us assume that we have two acoustically rigid plates of thickness h 1 and h 2 placed in the xy plane at z = 0. Both plates are periodically perforated by N round holes of radius R 0α and located at the positions R 1 α and R 2 α where α = 1, 2, . . . , N in a unit cell of size a 0 . These two plates are separated by a thin layer of thickness h g . We assume that the plates are made of steel or brass, in which the perfect rigid body approximation is very accurate. If the system is irradiated by an incident wave P in of unitary amplitude and propagating along the z axis with wavenumber k ¼ K þ q 0ẑ , the same scattering properties are obtained in different frequency regimes by scaling all the geometrical parameteres with the same factor.
In order to apply the Mode Matching Technique, five different regions are taking into account: region I refers to the first plate onto which sound waves impinge, region II refers to the field inside the through-holes of the first plate, region III refers to the gap separating the plates, region IV refers to the field inside the holes of the second plate and in region V sound is able to emerge into free-space. We assume that the surrounding medium is air, with the mass density ρ air and speed of sound c air .
In regions I, III and V the pressure is expanded in terms of plane waves whereas inside the holes the field is written as a linear combination of the corresponding waveguide eigenmodes: where G is the reciprocal lattice vector, K is the Bloch wave vector, k 0 = ω/c air is the incident wave vector, φ a = h 1 , φ b = h 1 + h g and φ c = h 1 + h g + h 2 are the phases we introduce in order to center the axis in the bottom or upper part of the layers, . The radiated and hole impedance are written as Z G = k0/ q G and Z = k 0 /q 0 , respectively, whereas r G (t G ) is the reflection (transmission) coefficient and A α , B α are the modal cavity expansion coefficients. Based on the subwavelength scales of the holes, it is assumed that only their fundamental mode is exited, so q 0 = c 1 /c 2 k 0 .
In the following, we show the normal velocity within the different regions, which is obtained deriving the pressure of each region with respect to the z-axis.
The Mode Matching Technique is applied by projecting the Bloch modes with the normal velocity field and the cavity modes for the pressure field: • Pressure: First, we impose the continuity of the pressure at the upper (z = 0) and bottom (z = h) part of the α holes: By subtituting the pressure field in the previous expressions: After doing the integrals we get: • Normal velocity: we impose the continuity of the normal velocity through along the entire unit cell.
By subtituting the normal velocity field in the previous expressions: After doing the integrals we get: After applying the mode matching technique through the pressure and the normal velocity field, we obtained a system of 8 equations and 8 unknowns wich are: r G ; τ G ; Γ G ; t G ; A II α ; B II α ; A IV α and B IV α . In order to solve the system, we are going to first eliminate four of these unknowns by reducing the system from 8 to 4 equations. For doing that, we will combine the result of applying the mode matching through the normal velocity of one region, with the result of applying the method to the pressure field as we can see in the following equations: Z A II