Long-range nonspreading propagation of sound beam through periodic layered structure

Linear spreading of a wave packet or a Gaussian beam is a fundamental effect known in evolution of quantum state and propagation of optical/acoustic beams. The rate of spreading is determined by the diffraction coefficient D which is proportional to the curvature of the isofrequency surface. Here, we analyzed dispersion of sound in a solid-fluid layered structure and found a flex point on the isofrequency curve where D vanishes for given direction of propagation and frequency. Nonspreading propagation is experimentally observed in a water steel lattice of 75 periods (~1 meter long) and occurs in the regime of anomalous dispersion and strong acoustic anisotropy when the effective mass along periodicity is close to zero. Under these conditions the incoming beam experiences negative refraction of phase velocity leading to backward wave propagation. The observed effect is explained using a complete set of dynamical equations and our effective medium theory. Recently dispersion-free propagation of spatiotemporally-modulated optical pulses over long distances has been demonstrated, but acoustic analogues are less well-explored. Here, anomalous dispersion and strong anisotropy is demonstrated for the long-distance non-spreading propagation of acoustic waves through an underwater metamaterial.

T he dependence of the index of refraction on wavelength, n = n(λ) known as dispersion, is responsible for a wide variety of effects in optics and acoustics. In particular, dispersion is the main reason for broadening of a wave packet propagating in a dispersive medium. Only in a hypothetical dispersionless medium does a signal propagate without changing its shape. However, even in this case a collimated beam of initial width σ 0 conically spreads. Spreading is due to the presence of the Fourier components with wave vectors making angle Δθ~λ/(σ 0 n) with direction of propagation.
It was demonstrated recently that initially spatiotemporaly modulated optical pulses may propagate in a linear medium without dispersion and diffraction 1 . In these 3D light bullets, the effects of dispersion and diffraction are not manifested and evolution of the packets is reduced to translational displacement without form change. Initial localization of the wave packet in 3D is achieved due to transversal modulation by a Bessel function and due to longitudinal and temporal modulations by an Airy function. Each function represents an exact solution of the corresponding wave equation which exhibits shape-invariant propagation [2][3][4][5] .
In acoustics, the area of dispersion/diffraction-free propagation of signals is less developed. Generation and diffraction-free propagation of a bottle-shaped localized acoustic pulses in a homogeneous medium was demonstrated by Zhang et al. 6 . A helical structure metamaterial possessing a large effective index of refraction for audible sound was proposed by Zhu et al. 7 . Despite high values of n eff~1 0, the index of refraction was practically independent of frequency, i.e. a slow sound wave did not suffer from dispersion spreading. Dispersionless propagation of mechanical gravity waves with nonlinear spectrum was observed on the surface of a water tank 8 . In all these experiments, as well as in the experiments with light pulses, the nonspreading propagation becomes possible due to a specially prepared initial wave packet.
In practical applications of acoustics in industry, biology, or medicine the most common signal generated by a transducer is a Gaussian beam. Divergence of a beam reduces sound energy concentration, resolution of acoustic imagining, and depth of penetration of the sound wave into a scanning medium. These and other drawbacks of acoustic devices may be improved by a robust and inexpensive method of collimation of acoustic beams.
In this paper, we report long-distance and nonspreading propagation of an acoustic Gaussian beam through a periodic stack of solid-fluid layers with hyperbolic dispersion. Strong suppression of the dispersion broadening of a Gaussian beam is predicted and observed for a set of angles of incidence and frequencies of sound when the dispersion coefficient vanishes exactly. This effect is due to some peculiarities in the dispersion law of the mixed acoustic mode propagating in strongly anisotropic layered environment.

Results
Evolution of a Gaussian beam. Propagation of a Gaussian beam in a dispersive medium is a standard problem of wave dynamics. Let a monochromatic source generate a pressure beam at x = 0 of width σ 0 along axis y, pðx ¼ 0; y; tÞ ¼ p 0 exp ik 0y y À y 2 =4σ 2 0 À iω 0 t . The beam propagates in an anisotropic elastic medium with nonlinear dispersion ω = ω (k x , k y ). The direction of the z-axis is chosen in such a way that the wave vector has two nonzero components, k 0 = (k 0x , k 0y , 0). To emphasize the effect of nonspreading propagation we assume that the beam is homogeneous along z, i.e. its width along z is much greater than along y. Under this condition, z-dependence of pressure in the wave can be ignored.
The distribution of pressure in the region x > 0 is obtained by integration over k y of the spatial Fourier components of p(x = 0, y, t) with phase-accumulation factor exp½iðk x x þ k y yÞ. Asymptotical behavior of this integral is obtained by quadratic expansion of the diffraction relation k x = k x (k y , ω 0 ) near k 0y (see, e.g. ref. 9 ) pðx; y; tÞ ¼ p 0 e iðk 0 ÁrÀω 0 tÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 À i Dx Here σðxÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi is the coordinate-dependent width of the beam, Vðk 0y Þ ¼ ∂k x =∂k y j k y ¼k 0y is the inclination parameter, and Dðk 0y Þ ¼ ∂ 2 k x =∂k 2 y j k y ¼k 0y is the diffraction coefficient of the medium [9][10][11] . At short distances, x < σ 2 0 =jDðk 0y Þj, the diffraction can be neglected. At longer distances, x ) σ 2 0 =jDðk 0y Þj the width of the beam grows linearly with x and the beam approaches a conical shape with the opening angle Δθ ¼ arctanðjDjÞ. Thus, the diffraction coefficient D defines the rate of beam spreading in spatial domain. The diffraction is normal if D < 0, otherwise it is anomalous.
The inclination parameter V(k 0y ) gives the direction of the transverse shift off the center of the Gaussian beam. The center of the beam follows the straight line y = −Vx, which is parallel to the vector of group velocity V g ¼ ð ∂ω ∂k x ; ∂ω ∂k y Þ. In an anisotropic medium the group V g and phase V ph ¼ ðω=kÞk velocities are not parallel.
Linear diffraction spreading does not occur if D(k 0y ) = 0. This condition can be realized if for certain frequency and direction of propagation the corresponding isofrequency surface has a region with zero curvature. In particular, in a medium with hyperbolic/ elliptic dispersion where ωðkÞ ¼ This can be realized if the dynamic mass ρ x = 0, i.e. if the elastic medium exhibits very strong anisotropy, ρ y /ρ x → ∞. In what follows we demonstrate that the dynamic mass density ρ x may be very small in a periodic layered water-steel lattice that gives rise to long-range nonspreading propagation of an acoustic beam.
Nonspreading propagation of an optical beam due to zero curvature of isofrequency surface was realized in a 3D photonic crystal 12 . Later, anomalous diffraction and diffractionless propagation of optical beams in discrete optical waveguides were studied in refs. 9,10,13 . Anomalies in refraction of sound related to curvature have been observed in ref. 14 . A general consideration of anomalous diffraction in metamaterials is proposed in ref. 11 .
Diffractionless propagation of acoustic beams was predicted 15,16 and observed in the experiment on acoustic transmission through perforated layered elastic structures 17 . It was shown that for these structures, the isofrequency curves have hyperbolic topology. While the elastic structures studied in refs. [15][16][17] were layered, they possess 3D periodicity due to perforation of solid plates. Our goal is to demonstrate that the diffraction coefficient D(k 0y ) vanishes in a layered solid-fluid structure with one-dimensional periodicity.
Effective dynamical parameters of a periodic layered system. Let us consider propagation of sound waves in a periodic binary system with layers of materials a and b and lattice spacing d = a + b. In the experiment, material a is water and b is steel. The geometrical and physical parameters of the structure are given in Fig. 1. In the long-wavelength limit kd ≪ 1, a layered structure behaves like a homogeneous anisotropic fluid with effective mass densities ρ x and ρ y = ρ z and effective elastic modulus B. Frequency-dependent effective parameters (EPs) of a layered structure can be calculated using the method of homogenization 18 for a binary structure with high acoustic contrast between the constituents. Generalization of the results 18 to a structure with arbitrary contrast is given in Supplementary Note 1: Here k aðbÞ ¼ (2). However, this approximation turns out to be sufficient to semi-quantitatively analyze the effect of nonspreading propagation. Note that the EPs depend on ω and k accounting for temporal and spatial dispersions. In the quasi-static limit, ω → 0 formulas (2) are reduced to Rytov's results 19 , Analysis of isofrequency curves. The dispersion relation with the EPs (2) reproduces well the long-wavelength region of the band structure obtained for the water-steel layers. In Fig. 1 the first three transmission bands are plotted for sound Bloch wave with Bloch vector k x varying within the first half of the Brillouin zone and conserved k y = 89 m −1 propagating in water-steel periodic lattice. Exact dispersion curves account for transverse and longitudinal modes are shown in black and those obtained in the effective medium approximation in red. The second band has anomalous dispersion where ∂ω/∂k x < 0. It starts at frequency f = 106.8 kHz. When the angle γ = 0 the band edge is red-shifted to the frequency f 0 = 104.8 kHz (see Supplementary  Fig. 3).
Near the frequency f 0 the dispersion relation, f = f(k x , k y ) is represented by a saddle surface shown in the inset to Fig. 2. Horizontal planes cutting the saddle surface above and below the level of f 0 give two sets of isofrequency curves topologically similar to hyperbolas. These two sets are separated by the isofrequency contour f(k x , k y ) = f 0 as shown in Supplementary  Fig. 4. The separatrix is not a straightline since the curves are not canonical hyperbolas. Two orthogonal parabolas passing through the saddle point are highlighted in red in the inset to Fig. 2. Projections of these parabolas on the horizontal plane form a local set of coordinates for the hyperbolic-like curves. Each coordinate line is the axis of symmetry of the hyperbolic-like curves. In Fig. 2 black and blue curves represent the two sets of hyperbolic-like isofrequency curves. The green line represents an isofrequency curve for f = 24 kHz with elliptic-like (normal) dispersion corresponding to the first transmitting band in Fig. 1.
For propagation along the lattice axis, k y = 0, the band gap occupies the interval f 1 < f < f 2 , where the gap edges f 1 = 46.7 kHz and f 2 = 77.5 kHz are reached at the edge of the Brillouin zone, k x = π/d = 314 m −1 . The isofrequency curves in the first passing band (f < f 1 ) are closed since equation f = f(k x , k y ) has real solutions for any direction of propagation. For the frequencies slightly above f 1 a narrow interval of the directions of propagation almost parallel to the axis, k y << k x becomes inaccessible. This interval increases with f. Thus, the topological transition from closed to opened isofrequency curves occurs at f = f 1 . Evolution of  The two blue curves are the plots for f = 105.1 kHz, which belongs to the second transmission band with anomalous dispersion. The solid blue curve is obtained from the exact dispersion relation derived for a solid-liquid layered structure taking into account longitudinal and transverse modes in metal plates (see Supplementary Note 2). The dash-dotted blue line is the result that follows from the effective medium theory, where the transverse mode is neglected. Because of this approximation and also since the parameter k x d ≈ 0.5 is not sufficiently small to provide the quantitative validity of the effective medium theory, the solid and dash-dotted blue lines deviate from each other.
Topologically, both blue curves are similar to hyperbolas. However, they are not given by a canonical second-order hyperbolic equation. The dash-dotted blue curve is obtained from Eqs. (2) and (3). Since the EPs in Eq. (2) depend on k y , Eq. Thus, the incoming wave suffers negative phase velocity refraction at the water-lattice interface and backward propagation inside the lattice. Similar effects were predicted for perforated layered structures 16 .
Refraction of a sound beam of f = 105.1 kHz at water-lattice interface is demonstrated by Supplementary Movie 1. It is seen that the wave fronts in water propagate towards the lattice, i.e. along positive direction of axis x. However, entering the lattice under the angle of incidence θ = 10°the x-projection of their velocity is reversed. It is a signature of negative phase velocity refraction. It is interesting to note that negative group velocity refraction cannot be observed at normal incidence. Unlike this, negative phase velocity refraction and backward wave propagation exist even for normal incidence. This Supplementary Movie 2 demonstrates backward wave propagation in steel-water structure for normal incidence of acoustic beam with frequency f = 104.5 kHz lying in the second transmission band with anomalous dispersion. However, the dispersion is normal in the third band (f = 155.4 kHz) that results in positive phase velocity refraction (see Supplementary Movie 3).
The flex point on the dash-dotted blue curve in Fig. 2 lies at k y ≈ 71 m −1 (see Fig. 2). At this point the dynamic density ρ x ≈ 32 kg/m 3 is much less than the densities of the lattice constituents. It does not vanish exactly (unlike the case of elliptic/hyperbolic dispersion) since the EPs in Eq. (3) are k y -dependent. We may conclude that the nonspreading propagation of acoustic beam occurs under the condition of very strong acoustic anisotropy, ρ y / ρ x ≈ 40 (see Supplementary Note 1 and Supplementary Fig. 1). One more visible signature is that at the flex point, group and phase velocities are almost orthogonal. Also the dispersion curve of the second transmission band in Fig. 1 is quite flat, that leads to the inequality |V g | << |V ph |, which is related to strong scattering at the interfaces 20 .
The isofrequency curves in Fig. 2 obtained from the exact solution of equations of elasticity exhibit discontinuous cuts. Sound waves with frequency and direction corresponding to the cut do not propagate though a given layered structure. The effect of zero acoustic transparency through a solid plate surrounded by a fluid was recently predicted by Quotane et al. 21 . This effect is due to a narrow Fano resonance caused by interference between a discrete trapped mode of a solid plate and a continuum of acoustic modes of a fluid. In the experiment, the direction of propagation of this Fano resonance is realized at must be ruled out. Of course, the blue dash-dotted line does not show Fano resonance since the EPs are obtained neglecting the transverse mode.
Using Snell's law sin θ= sin γ ¼ c a =V ph;ref the angle of incidence is calculated to be θ ≈ 10°. This angle is counted from the positive direction of axis x since the incident wave propagates in water where V g;inc Á V ph;inc ¼ c 2 a > 0. Note, that the angle of incidence calculated using the effective medium theory is 9°that is close to the exact result (see Supplementary Note 2).
Numerical modeling and experimental results. Numerically modeled transmission of a sound beam of frequency f = 105.1 kHz and angle of incidence θ = 10°through the steel-water structure is shown in Fig. 3. The width of the beam in the x−y plane practically does not change along its pass through a 100period long lattice since for the selected parameters the diffraction coefficient D ¼ ∂ 2 k x =∂k 2 y vanishes. However, in the 3D map of sound intensity shown in Fig. 4, the width of the beam in the y−z plane grows linearly because the corresponding diffraction coefficient ∂ 2 k x =∂k 2 z ≠ 0 (see Supplementary Note 3). The spreading of the beam eventually results in interaction between the beam and the top-bottom boundaries that manifested in unwanted multiple reflections from the acrilic holders. However, the net effect of these reflections is minimal since the impedance mismatch between the acrylic and water produces a low reflection coefficient of~8%. Also only a small part of the acoustic beam expands over the edges of 17.8 × 17.8 cm plates, i.e. very small part of the signal that reaches the receiving transducer suffers from multiple reflections.
The directions of phase and group velocities in the incident and refracted beams are indicated by arrows. The direction of the group velocity in the refracted beam can be calculated from Snell's law using the definition of the group index of refraction 16 Since in Fig. 2 the derivative ∂k x /∂k y is negative at the flex point (k y = 0.89 cm −1 with V ðxÞ g > 0), the beam exhibits positive refraction for group velocity, n g > 0. The angle of refraction calculated from Eq. (4) is γ g = 35 ∘ .
This angle defines the geometry of the fabricated lattice. Due to strong refraction of the acoustic beam at water-lattice interface, the steel plates in the structure are vertically shifted as shown in Fig. 5a, b. The input Gaussian beam of width σ 0 ≈ 2 cm propagates a distance 75= cos 35 % 88 cm through 75-period long lattice. A comparison of the beam spreading in free water versus the lattice is given in the plot inset of Fig. 5b. At 71 cm from the emitter, the beam spreads from ≈2.7 to ≈10 cm through free water. However, the same input spreads to only ≈3.5 cm when propagating through the lattice. Thus, the effect of diffraction broadening is strongly suppressed in the lattice supporting the contention that the diffraction coefficient D is close to zero.
The condition of diffractionless propagation, D = 0, establishes a dependence among three variables, k x , k y , and f. This dependence can be reduced to the relation between the angle of incidence θ and frequency f (see Supplementary Note 3). The dependence f = f(θ) is plotted in Supplementary Fig. 7. The frequency f = 105.1 kHz is selected as a point on the curve f = f(θ) corresponding to relatively small angle of incidence (θ = 10°), fitting the second transmission band close to the Γ point, and lying away from Fano resonance.
The intensity of the signal (at the center of the beam) transmitted through the lattice is reduced by 930 times as compared to the signal intensity at the emitting transducer. The ratio of the intensities averaged over the beam profile is about 500. The reduction of the transmitted signal is due to strong reflection at the first and last metal plates, to the viscous losses and scattering inside the lattice, and to losses outside the lattice. The latter factor is related to the strongly inhomogeneous acoustic field generated by the transducer displaced 8 cm away from the first plate. Sound intensity measured near the first plate is ≈5.5 times lower than at the transducer. The averaged signal  intensity also slightly decays (≈2 times) on the way (4 cm) from the last plate to the receiving hydrophone. Thus, we estimate that, passing through the lattice, the signal intensity is reduced by 500/ 5.5 × 2 = 45 times. Numerical modeling shown in Fig. 3 gives the decay of the signal intensity by 36 times. This includes the losses due to impedances mismatch and to viscosity. We attribute 45/ 36 ≈ 1.3 times higher losses in the experiment mostly to scattering at period fluctuations and absorption of sound by small amount of air bubbles which are always present in water. The transmission can be optimized by replacing stainless steel plates by plates of softer material (see Supplementary Note 4 and Supplementary Fig. 8).
The regime of backward propagation of sound inside the lattice was demonstrated by measuring the phases of the wave at six successive unit cells. The results are presented in Fig. 6 and the details are given in Supplementary Note 5. Experimental and numerical results are shown by squares and dots connected by solid and dashed line, respectively. Red and blue symbols show the phases of two passing waves with frequencies lying in the first (red) and second (blue) passing bands. In the experiment the phases of the Bloch wave that exhibit nonspreading propagation at f = 112.7 kHz (second passing band) gradually decrease that serves as a signature of backward propagation. Unlike this, the phases of the wave with normal dispersion with f = 52.7 kHz (first passing band) gradually increase. Numerically calculated phase behaves similarly. Since in the experiment the nonspreading propagation was observed at the frequency ≈7 kHz higher than the theory predicts, only qualitative agreement is expected. Note that for backward propagation the experimental and numerical phases change much slower than for forward propagation. This occurs because in both cases the frequencies are close to the upper edge of the bandgap. Exactly at the band edge spatial variation of the phase is compensated by its temporal variation since a propagating Bloch wave becomes a standing wave.
Discussion. The divergenceless propagation results in a collimated beam as shown in Fig. 5b. The emerging beam from the lattice does not exhibit strong spreading over the 6 cm from the end of the lattice to the end of the scanned area. From Fig. 5b inset, it can also be seen that propagation through the lattice maintains a strong intensity gradient at the beam edges as compared with the weaker intensity gradient observed at comparable distance without the lattice. Collimation of acoustic beam and energy funneling were previously observed due to diffraction of incoming beam at corrugated surfaces [22][23][24] . This method partially reduces transverse diffraction of acoustic beam. Resonant interaction of the incoming sound with the surface acoustic wave excited at a corrugated surface also reduces angular spreading of the beam, while numerical simulations 22 show presence of lateral diffraction maxima. Lateral diffraction disappears if the corresponding isofrequency surface has a flat region 11,14,25 . Usually, the curvature never vanishes completely, leading to some spreading. Experimental realization of a flat band in acoustics was demonstrated in a 2D phononic crystal 25 of steel rods in water, and more recently in perforated Plexiglas plates in air 17 . Simultaneous collimated propagation of acoustic and electromagnetic waves in phoxonic crystal of silicon rods in air was numerically predicted 26 . Dissipative losses were not taken into account and high level of collimation is achieved due to existence of a single flex point on isofrequency line.
In these periodic structures where collimation of a sound beam was predicted or observed, sound usually propagates through narrow apertures that strongly enhance viscous losses. Therefore, the length of the collimated beam did not exceed 10 cm. The scheme proposed here is free from such drawback. Much less energy losses allow much longer collimated propagation. Also the angular width of the beam is much less since it is required that curvature vanishes exactly only at a single point. From practical point of view the 1D periodic lattice proposed here is much easier to fabricate.
The exit position of the outgoing beam exhibits strong dependence on frequency. Experimentally, nonspreading propagation occurs at 112.17 kHz. Even small variation of the frequency leads to noticeable vertical displacement of the output beam, as shown by the colored beam profiles in Fig. 5c. This displacement is due to dispersion of the refraction coefficient Eq. (4). The graph in Fig. 5c shows calculated vertical displacement caused by small variation of the refraction coefficient (4), and the angle of refraction, with frequency. Calculated and measured displacements are in reasonable agreement.
While the experimental results confirm most of our theoretical predictions for nonspreading propagation, there is a disagreement over the frequency. The effect is observed at frequency 7 kHz higher than the theory predicts. We attribute this difference to partial vibration coupling and subsequent minor mode shifting of the lattice modes due to the two acrylic plates which hold the metal plates on the top and bottom (see Fig. 5a). These holders are not sufficiently massive to rule out sound transmission, therefore the stiffness of the system as a whole turns out to be higher than that of the steel-water lattice.
We have demonstrated long-range nonspreading propagation of sound through a layered solid-liquid structure with onedimensional periodicity. Diffraction, usually leading to linear divergence of the beam, is avoided due to vanishing of the diffraction coefficient. Since this coefficient is proportional to the curvature of the isofrequency curve, the possibility of diffractionless propagation depends on the presence or absence of a flex point on the curve and the direction of propagation is defined by the position of the flex point. The isofrequency curves are obtained from the exact solution of the wave equation in the layered structure and from the developed effective medium theory, which accounts for temporal and spatial dispersions of the EPs. Both methods give close results. We have shown that the effect of nonspreading transmission is accompanied by two other anomalies: negative phase velocity refraction and strong acoustic anisotropy. The reported results may find applications in many areas where robust transmission of narrow acoustic beams is required.

Methods
Experiment. Our sample is assembled from 75 pieces of 17.78 cm square shape grade 409 stainless steel plates with thickness b = 0.891 mm. The plates are 10 mm apart, center-to-center with 0.1 mm tolerance. The plates are held in place by specially fabricated acrylic top and bottom plates.
The acoustic experiment was conducted in a large deionized water tank (150 × 90 × 60 cm), in a bistatic setup. Signal generation was achieved with a Lecroy WaveStation2012 Arbitrary Waveform Generator set to 20 V peak-to-peak, sweeping from 80 to 120 kHz in 80 s. A 25 mm diameter Ultran plane wave NCG100-D25 immersion transducer served as the emitter. Detection was done using a Müller-Platte 1 mm needle hydrophone connected to a Tektronix MDO 3024b in spectrum analyzer mode. The emitting transducer was rotated in the x−y plane by a Thorlabs CR1-Z6 angular translation stage for variable incident angles. The needle hydrophone was translated in the x−y plane by a two-dimensional 1500 mm by 1200 mm Newmark eBelt linear translation stage controlled by a Newmark NSG-G2-X2 translation stage controller. Calibration of the angle and subsequent alignment was performed for each experiment with angular tolerance of 0.5. An automated raster scan was completed by a pre-prepared MATLAB code with data acquired at each point for 90s. The incident angle with respect to the lattice plane was 10°. For the incident scanned area, the beam profile was measured 60 mm from the incident lattice plate in a 9 cm × 4 cm area in 3 and 10 mm steps, respectively (Fig. 5). The output signal was measured starting 20 mm from the last plate with a scanned area of 9 cm × 4 cm in 3 and 10 mm intervals.
Theory. The EPs Eq. (2) are obtained from two dynamical conditions that define the response of elastic medium on external perturbation.
First, the average acceleration of the unit cell in a Bloch wave equals to that calculated for the equal volume of the effective medium under action of propagating plane wave of same amplitude a i ¼ À Here p uc (x, y) is the distribution of pressure within the unit cell, p h (x, y) is the distribution of pressure in a homogeneous effective medium, and ρ(x) is the coordinate-dependent density within the unit cell. Second, the average deformations are also equal Àd=2 p uc ðx; yÞ λðxÞ dx Here λ(x) is the coordinate-dependent bulk modulus within the unit cell. Equations (5) and (6) define the EPs not only for a binary structure but for any distribution of density ρ(x) and elastic modulus λ(x) within the unit cell. Derivation of the EPs ρ x , ρ y , and B is given in Supplementary Note 1. The dispersion relation for fluid-solid lattice is obtained from the condition of matching the solution of the wave equations for longitudinal mode in layer a (fluid) with the solution for longitudinal and transverse modes in layer b (solid). Together with the Bloch theorem for periodic medium, this leads to a set of 6 × 6 linear homogeneous equations. Numerical simulations are conducted by using an acoustic model in commercial software (COMSOL MULTIPHYSICS).