Main

Sound is probably the ‘simplest’ of all classical waves. It has no intrinsic spin and does not respond to magnetic fields; hence, the fundamental interactions underlying the quantum spin Hall effect (QSHE) and quantum Hall effect (QHE) do not apply to acoustic waves. However, by engineering the coupling in the z direction, we can create synthetic staggered flux and hence kz-preserved unidirectional edge modes in the xy plane in acoustic meta-crystals which have simple ‘static’ structures with no moving fluid and no dynamic modulation. The acoustic meta-crystals have Weyl points24,25,26 in the three-dimensional (3D) band structure, and the systems are acoustic analogues of the topological Haldane model2,23 for fixed values of kz.

Tight-binding model

To illustrate how the idea works, we start with a simple nearest-neighbour tight-binding model for an AA-stacked honeycomb lattice (Fig. 1a). Its acoustic implementation will be discussed later. The Hamiltonian H of this system consists of the intralayer part H0 and the interlayer part H1.

where a (b) and a (b) are the annihilation and creation operators on the sublattice sites, ɛ represents the on-site energy difference. Each lattice is specified by subscripts (i, k), where the first labels the position in each layer and the second labels the number of layers. The first term in equation (1) represents the sublattice on-site energy difference. The second term represents the intralayer hopping between nearest sublattices. The intralayer hopping tn is real and constant. The first Brillouin zone (BZ) of this lattice is shown in Fig. 1c. As this system is periodic along the z direction, kz is a good quantum number. For each fixed kz, and if we consider the dispersion and transport in the xy plane, the 3D system can be reduced to an effective two-dimensional (2D) system with a unit cell, as shown in Fig. 1b. The corresponding first BZ is illustrated in Fig. 1d, which represents a plane cut with the specified kz in the original first BZ in Fig. 1c.

Figure 1: Interlayer coupling induces different symmetry breakings in two-dimensional systems.
figure 1

a,b, AA stacking of a hexagonal lattice (a) and its unit cell (b). c, Reciprocal space of the system shown in a. d, Reduced first Brillouin zone for each fixed kz. Different interlayer couplings introduce different effects in the reduced two-dimensional lattice. e,f, Unequal interlayer coupling at different sublattice sites corresponds to broken inversion symmetry in 2D. g,h, Chiral interlayer coupling generates an effective gauge flux at each fixed kz. The red arrows in h represent the direction of positive phase hopping, and dotted and crossed black circles represent the direction of local flux. This is an analogue of the Haldane model.

We can now introduce different kinds of interlayer coupling, from which we choose two special examples (Fig. 1e, g, where hopping is non-zero only between connected sites). In Fig. 1e, the hopping amplitudes are different for the two different sublattices. In this case, the interlayer hopping part H1 is given by

where ta and tb (tatb and both are real) represent the interlayer hopping terms for sublattices a and b. The corresponding Bloch Hamiltonian H(k) is given by

where , a is the distance between the two sublattices, dh is the interlayer distance and k = (kx, ky, kz) is the Bloch wavevector. The eigenvalue of H(k) is given by . The first term under the square vanishes along the KH line. If [ɛ + (tatb)cos(kzdh)] is non-zero, the inversion symmetry of the reduced 2D hexagonal lattice is broken, as illustrated in Fig. 1f. When |tatb| > ɛ, there exist special values kz = ±arccos[ɛ/(tbta)]/dh, where H(x, y, kz) = H(−x, −y, kz). This implies that Dirac cones can form in the kxky plane for these values of kz. These special points can give rise to Weyl points in the 3D band structure, as we will discuss later.

Another interesting example is shown in Fig. 1g, in which we consider a chiral kind of interlayer coupling. We will see that if we consider the propagation in the xy plane for a fixed kz, this is effectively a realization of the topological Haldane model2. The interlayer coupling coefficients as indicated by the cyan bonds in Fig. 1g are denoted by tc, which are taken to be identical and real. For a fixed kz, the interlayer hopping becomes next-nearest-neighbour hopping in the reduced 2D system, with a complex hopping coefficient tce, where φ = ±kzdh depends on whether the hopping proceeds in the clockwise or anticlockwise direction. In Fig. 1h, we use red arrows to indicate the direction along which φ = −kzdh. After a complete loop of hopping in the direction indicated, the total phase accumulated is −3kzdh. This means the gauge flux enclosed by this loop is −3kzdh. In Fig. 1h, we use dotted and crossed black circles to represent the sign of local flux. Although the total flux inside each unit cell is zero, we have non-zero local flux. The chiral interlayer coupling introduces Peierls phase27 (Supplementary Information II) for the hopping parameters in 2D for any non-zero kz, thereby achieving a staggered synthetic gauge flux in a ‘static’ system. This is different from the standard paradigm of using dynamical perturbation to induce synthetic gauge flux17. For a system with both the coupling shown in Fig. 1e and that in Fig. 1g, one can then tune the system across a topological transition point by changing the value of (tatb)/tc or kz. The phase diagram is shown in Supplementary Information I.

Breaking of inversion symmetry

Let us now consider real acoustic systems. We start with a periodic array of acoustic cavities linked together by tubes, as shown in Fig. 2a (top view) and b (side view). The resonance cavities can be viewed as ‘meta-atoms’ and the hopping strength between the meta-atoms can be tuned by changing the radius of the connecting tubes. The light blue colour in Fig. 2b denotes the area where the hard boundary condition is applied and the system is filled with air. The intralayer couplings are set to be equal by giving all horizontal connecting tubes the same radius (w0). The interlayer couplings (along the z direction) are set to different values by choosing different radii (w1w2) for different sublattice sites. Here we consider the mode (Fig. 2c) whose pressure is described by a sinusoidal function along the z direction and does not vary in the horizontal plane. We note here that different mode profiles give basically the same results, except that the working frequency will be different. In Fig. 2d, we show the band dispersions of this mode in the kxky plane with different values of kz. At kz = 0 (red lines in Fig. 2d), the dispersion in the kxky plane has a gap at the point owing to inversion symmetry breaking. As this reduced 2D system still has mirror symmetry with respect to the xz plane, the Berry curvature is an odd function in the reciprocal space and the two red bands are topologically trivial, with zero Chern numbers. At kz = 0.623, there is a Dirac point at , as shown by the black lines in Fig. 2d, consistent with the tight-binding model’s prediction that degeneracy for the reduced 2D system can be recovered at some finite kz value. The value of kz, where the system has a Dirac cone in the kxky plane, is slightly different from the tight-binding prediction. This is because the introduction of the connecting tubes changes the resonance frequency of the cavities, which is equivalent to ɛ ≠ 0 in equation (2). In Fig. 2e, we show the band dispersion (black curves) along the KH direction (z direction). The difference in frequency between the two bands along kz reflects the strength of inversion symmetry breaking as a function of kz. As shown in Fig. 2e, the band dispersion near the Dirac point in Fig. 2d is also linear along the z direction, indicating that the degeneracy point is a Weyl point24,25,26,28 in the 3D band structure.

Figure 2: Different interlayer coupling coefficients induce effective breaking of in-plane inversion symmetry.
figure 2

a,b, Top and side views, respectively, of a unit cell of the acoustic system. c, Real part of the eigenpressure field of a cavity mode for a single resonator at 2,144 Hz, where red and blue represent positive and negative pressure, respectively. d, Band dispersion of the cavity mode in c for different values of kz, where the black curves have a Dirac point at for a special value of kz and the degeneracy is lifted (red curves) otherwise. The Chern number of the two red bands are both C = 0. e, Band dispersion along the KH direction. The two bands cross at a special point (kz = 0.623 with normalized unit π/(hc + l)). This special point corresponds to a Weyl point of this acoustic system.

Synthetic gauge flux

We now consider the realization of effective acoustic gauge flux. In fact, creating effective acoustic gauge flux for a fixed kz with chiral interlayer coupling can go beyond the tight-binding description. Figure 3a, b shows the top and side views of a unit cell of an acoustic system that exhibits synthetic gauge flux. The chiral coupling is characterized by the relative rotation angle θ between the holes on the upper and lower boundaries of the planar (xy plane) waveguides. When θ = 0, the coupling is non-chiral and the synthetic gauge flux vanishes. The strength of the gauge flux depends on the rotational angle as well as the radius of the connecting tube r1 (Supplementary Information IV).

Figure 3: Chiral interlayer coupling induces a synthetic gauge flux.
figure 3

a,b, Top view and side views, respectively, of a unit cell of the acoustic system. The dashed circles and solid circles in a represent the holes opened at the upper and lower sides of the in-plane (xy) sound waveguide, respectively. θ represents the rotation angle of the connecting tubes. c, Band dispersion of the lowest mode for two values of kz, where the black/red curves represent bands without/with the effective gauge flux. The Chern numbers of the two red bands are C = +1 (lower band) and C = −1 (upper band). The sign of Chern numbers can be changed by reversing the sign of kz. d, Band dispersion along the KH direction in the reciprocal space. K and H points are Weyl points of the 3D band structure. The strength of the synthetic gauge flux reaches its maximum near kz = 0.5π/(hc + l). e, Projection band (grey) along the x direction with kz = 0.5π/(hc + l). Red and blue curves represent the surface states localized at the lower and upper boundaries of a ribbon, respectively (see Supplementary Fig. 5, in which the eigenpressure fields of the states marked by the black stars are shown). f, The surface wave for kz = −0.5π/(hc + l) propagates in the clockwise direction around the corners and the defect of a finite system without being backscattered. The purple star marks the position of our sound source, the black arrows illustrate the direction of propagation of the sound wave inside the system and the purple arrows at the bottom show the coupling of the sound wave from the boundary of the system to the outside.

We consider the lowest-order acoustic mode and, as before, the waveguide is filled with air. In Fig. 3c, red/black curves indicate the band dispersion in the reduced 2D BZ at kz ≠ 0/kz = 0, representing the system with/without the effective gauge flux. The degeneracy at is lifted by the gauge flux. The effect of this synthetic gauge flux can also be seen from the Chern number of each isolated band, which is found to be +1/ −1 for the lower/upper band when kz is positive. The strength of the gauge flux can be seen from the width of the gap at the point. In Fig. 3d, we show the band dispersion along the KH direction. The two bands are required to be degenerate at the K and H points by a combination of time-reversal and C6 rotational symmetry. Along the KH line, the two states repel and the width of the gap reaches its maximum near the middle of the KH line.

A non-zero Chern number for a non-zero kz implies the existence of a topologically protected chiral edge mode in the boundary between this system and a topologically trivial system inside the common gap region. We construct a hexagonal ribbon with finite width along the y direction and periodic along the x direction (see Supplementary Information V). We use hard boundaries to confine the sound wave in the y direction to within the acoustic system. The hard boundary condition can be regarded as a trivial band gap with zero decay length. Figure 3e shows the projection band (grey) along the x direction and the dispersions of two surface states with kz = 0.5π/(hc + l). The red and blue colours denote surface states localized at the lower and upper boundaries in the y direction, respectively (see Supplementary Information V). The two degeneracy points in Fig. 3d correspond to Weyl points24,25,26,28 in the 3D band structure. If we consider a surface BZ spanned by kx and kz, the allowed boundary modes for a given excitation frequency will trace out trajectories connecting two Weyl points that are analogous to ‘Fermi arcs’ in electronic Weyl semi-metals25,28,29 (see Supplementary Information VI). If the width of the ribbon is large enough (larger than the decay length of the edge state), the edge state localized at one boundary cannot be scattered backwards as long as kz is still preserved. In Fig. 3f, we show the property of the edge state. The edge state is excited on the left boundary (marked by the purple star) inside the gap region, and propagates clockwise around the corners and the defect without being backscattered. Black and purple arrows are drawn to show the direction of propagation of the sound wave. The direction, either clockwise or anticlockwise, depends on the sign of kz. We supplemented our full wave simulation for an infinite system (with periodic boundary conditions) with simulations using a 3D finite-sized tight-binding model (Supplementary Information III). Different from the one-way edge states of 2D acoustic systems published recently20, the time-reversal symmetry in our 3D system is preserved. We note that the time-reversal partner of a clockwise state for a particular kz is an anticlockwise state at −kz.

Weyl points

This synthetic gauge flux is related to the Weyl points. In 3D, Weyl point dispersion is governed by the Weyl Hamiltonian H(k) = ∑ kivijσj, i, j {x, y, z} (refs 24, 25, 26, 28, 30, 31, 32), where vij are the group velocities and σj are the Pauli matrices. Weyl points have associated topological charges (or chirality, c = sgn[det(vij)] = ±1), which can be regarded as monopoles of Berry flux26. In Fig. 4a, b, we show the Weyl points of the two acoustic systems considered in Figs 2 and 3, respectively. The systems considered in Fig. 2 have mirror symmetry with respect to the xz plane, which ensures that the two mirror-symmetric Weyl points in the reciprocal space have charges with opposite signs, and hence the net charge in the horizontal light blue plane in Fig. 4a is zero. For any 2D band with an arbitrary (but fixed) value of kz, the net Berry flux vanishes and its Chern number is zero. In contrast, the system with chiral coupling (for example, Fig. 3) does not possess mirror symmetry and the remaining C6 symmetry ensures all the Weyl points on the same kz plane have charges with the same sign. Meanwhile, the net charge of Weyl points inside the first BZ must vanish29, which means there must be at least two planes with different kz possessing Weyl points of different charges. For example, the kzdh/π = 0 and kzdh/π = ±1 planes in Fig. 4b carry net topological charges of +2 and −2 respectively. Thus, for an arbitrary fixed kz lying between these two planes, the net Berry flux through the reduced 2D BZ is 2π and the Chern number is ±1, with the sign determined by the sign of kz (refs 24, 25, 29). The chiral coupling guarantees the non-zero Chern number, which corroborates the existence of a synthetic gauge flux.

Figure 4: Weyl points in the reciprocal space.
figure 4

a,b Weyl points for the acoustic systems studied in Figs 2 and 3, respectively. The red and blue spheres represent Weyl points with positive and negative topological charges, respectively.

Our idea of manipulating acoustic waves can also be extended to electromagnetic wave systems. The vector property of electromagnetic waves as well as the possibility of breaking time-reversal symmetry using magnetic fields in electromagnetic systems can offer more flexibility in introducing interesting phenomena such as analogues of chiral anomaly33. The results shown in Figs 2 and 3 demonstrate separately and respectively the consequences of on-site coupling difference and chiral coupling. If both symmetry-breaking mechanisms are simultaneously incorporated in the design of the acoustic meta-crystal, we can in principle enter all the regimes in the phase diagram of the Haldane model as we vary kz. The effective gauge flux induced by chiral coupling may stimulate new ideas for manipulating sound wave propagation, and may have implications in fields such as sound signal processing, sound energy harvesting, and noise protection.

Methods

All simulations were performed using the commercial solver package COMSOL Multiphysics. The 3D geometry was implemented by imposing periodic boundary conditions on some specified boundaries. The systems were filled with air (density ρ = 1.3 kg m−3 and speed of sound ν = 343 m s−1). Eigenmode calculations were carried out to find the band structures as well as eigenmodes in Figs 2c–e and 3c–e. Frequency domain calculations were carried out to find the kz conserved backscattering immune transport behaviour of the edge mode shown in Fig. 3f.

Figure 2a and b, respectively, show the top view and side view of the unit cell studied in Fig. 2. The parameters used were rc = 3 cm, a = 9 cm, w0 = w2 = 0.6 cm, w1 = 1 cm, hc = 8 cm and l = 3 cm. The light blue colour in Fig. 2b denotes the area where the hard boundary condition was applied, and the white colour indicates the area where the periodic boundary condition was applied with a given Bloch wavevector. Figure 2c shows the real part of the pressure field of a cavity mode found at 2,144 Hz, where the parameters of the cavity were the same as before—that is, rc = 3 cm and hc = 8 cm. The hard boundary condition was applied over all the cavity boundaries and the eigenmode calculation was then performed to find the eigenpressure distribution, where the red/blue colour indicates positive/negative local pressure.

Figure 3a and b, respectively, show the top and side views of the unit cell studied in Fig. 3. The parameters of the unit cell were r0 = 3 cm, rt = 5.5 cm, a = 8 cm, r1 = 1.2 cm, θ = 90°, l = 5 cm, hc = 8 cm. The dashed circles and solid circles in Fig. 3a represent the holes opened at the upper and lower sides of the sound waveguide in the xy plane, respectively. θ represents the rotation angle of the connecting tubes. In Fig. 3b, the light blue colour marks the area where the hard boundary condition was applied, and solid white circles represent areas where Floquet periodic boundary conditions were applied with a given Bloch wavevector. Floquet periodic boundary conditions were also applied to the side walls (they are transparent to expose the interior structure of the unit cell). To calculate the projection band shown in Fig. 3e, we constructed a ribbon with finite length (16 unit cells) along the y direction. Floquet periodic boundary conditions were then applied to the side walls in the x and z directions. The remaining boundaries were all set as hard boundaries. The Bloch wavevector along the z direction was fixed at kz = 0.5π/(l + hc) in this calculation. In Fig. 3f, the purple star marks the position of our source. To couple waves inside our system, we adopted the plane wave radiation boundary condition for one of the side boundaries of a unit cell with a plane wave excited at this port and a working frequency of 1,360 Hz. All the side boundaries in the negative y direction were also set as plane wave radiation boundaries to couple waves outside of our system, and we used purple arrows to indicate the direction of wave propagation through these ports. Here the ‘plane wave radiation boundary’ plays two roles in our simulation. The one at the side boundary serves as the source and those in the negative y direction at the bottom are used to couple the wave from our system to the outside so that the surface wave will not go around all the side boundaries and back to the source port boundary. The plane wave radiation boundary, which serves as source, can also be replaced by other kinds of source, such as a line source. All the boundaries on the upper and lower sides in the z direction are set as Floquet periodic boundary conditions with a Bloch wavevector along the z direction given by kz = −0.5π/(l + hc), where the ‘+’ or ‘−’ sign of kz determines the direction (clockwise or anticlockwise) of surface wave propagation. All the remaining boundaries are set as hard boundaries.