# Long-range nonspreading propagation of sound beam through periodic layered structure ### 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 left(i{k}_{0y}y-{y}^{2}/4{sigma }_{0}^{2}-i{omega }_{0}tright)). The beam propagates in an anisotropic elastic medium with nonlinear dispersion ω = ω(kxky). The direction of the z-axis is chosen in such a way that the wave vector has two nonzero components, k0 = (k0xk0y, 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 ky of the spatial Fourier components of p(x = 0, yt) 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 kx = kx(kyω0) near k0y (see, e.g. ref. 9)

$$p(x,y,t)=frac{{p}_{0}{mathrm{{{e}}}}^{i({{bf{k}}}_{{bf{0}}}cdot {bf{r}}-{omega }_{0}t)}}{sqrt{1-ifrac{Dx}{2{sigma }_{0}^{2}}}}exp left[-frac{{left(y+xVright)}^{2}}{4{sigma }^{2}(x)}left(1+ifrac{Dx}{2{sigma }_{0}^{2}}right)right].$$

(1)

Here (sigma (x)=sqrt{{sigma }_{0}^{2}+{left(Dx/4{sigma }_{0}right)}^{2}}) is the coordinate-dependent width of the beam, (V({k}_{0y})=partial {k}_{x}/partial {k}_{y}{| }_{{k}_{y} = {k}_{0y}}) is the inclination parameter, and (D({k}_{0y})={partial }^{2}{k}_{x}/partial {k}_{y}^{2}{| }_{{k}_{y} = {k}_{0y}}) is the diffraction coefficient of the medium9,10,11. At short distances, (x , < , {sigma }_{0}^{2}/| D({k}_{0y})|), the diffraction can be neglected. At longer distances, (xgg {sigma }_{0}^{2}/| D({k}_{0y})|) the width of the beam grows linearly with x and the beam approaches a conical shape with the opening angle (Delta theta =arctan (| D| )). 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(k0y) 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 ({{bf{V}}}_{text{g}}=(frac{partial omega }{partial {k}_{x}},frac{partial omega }{partial {k}_{y}})). In an anisotropic medium the group Vg and phase ({{bf{V}}}_{text{ph}}=(omega /k)hat{{bf{k}}}) velocities are not parallel.

Linear diffraction spreading does not occur if D(k0y) = 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 (omega ({bf{k}})=sqrt{{k}_{y}^{2}{c}_{y}^{2}mp {k}_{x}^{2}{c}_{x}^{2}}) the diffraction coefficient is always negative, (D({k}_{y})=-({c}_{y}^{2}/{c}_{x}^{4})({omega }^{2}/{k}_{x}^{3}) , < , 0). It vanishes if ({c}_{x}=sqrt{B/{rho }_{x}}=infty). 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 crystal12. 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 predicted15,16 and observed in the experiment on acoustic transmission through perforated layered elastic structures17. 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(k0y) 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 homogenization18 for a binary structure with high acoustic contrast between the constituents. Generalization of the results18 to a structure with arbitrary contrast is given in Supplementary Note 1:

$${rho }_{x}=frac{{rho }_{text{a}}frac{{k}_{text{b}}}{{k}_{text{a}}} , sin , frac{{k}_{text{a}}a}{2} , cot , frac{{k}_{text{b}}b}{2}+{rho }_{text{b}} , cos , frac{{k}_{text{a}}a}{2}}{1+frac{{k}_{text{b}}}{{k}_{text{a}}} , sin , frac{{k}_{text{a}}a}{2} , cot , frac{{k}_{text{b}}b}{2}-2frac{{rho }_{text{b}}}{{rho }_{text{a}}} , {sin }^{2} , frac{{k}_{text{a}}a}{4}},\ {rho }_{y}={rho }_{z}={rho}_{text{a}}frac{{k}_{text{a}}d}{2}frac{cos , frac{{k}_{text{a}}a}{2},cot , frac{{k}_{text{b}}b}{2}-frac{{rho }_{text{a}}{k_{text{b}}}}{{rho }_{text{b}}k_{text{a}}}sin , frac{{k}_{text{a}}a}{2}}{sin , frac{{k}_{text{a}}a}{2},cot , frac{{k}_{text{b}}b}{2}+frac{{rho }_{text{a}}}{{rho }_{text{b}}}left[frac{{k}_{text{a}}}{{k}_{text{b}}}-2frac{{k}_{text{b}}}{{k}_{text{a}}},{sin }^{2} , frac{{k}_{text{a}}a}{4}right]},\ B={lambda }_{text{a}}frac{{k}_{text{a}}d}{2}frac{cos , frac{{k}_{text{a}}a}{2},cot , frac{{k}_{text{b}}b}{2}-frac{{rho }_{text{a}}{k}_{text{b}}}{{rho }_{text{b}}{k}_{text{a}}},sin , frac{{k}_{text{a}}a}{2}}{sin , frac{{k}_{text{a}}a}{2},cot , frac{{k}_{text{b}}b}{2}+frac{{lambda }_{text{a}}{k}_{text{a}}}{{lambda }_{text{b}}{k}_{text{b}}}-2frac{{k}_{text{b}}{rho }_{text{a}}}{{k}_{text{a}}{rho }_{text{b}}},{sin }^{2} , frac{{k}_{text{a}}a}{4}}.$$

(2)

Here ({k}_{{text{a}}(b)}=sqrt{{(omega /{c}_{{text{a}}(b)})}^{2}-{k}_{y}^{2}}) is the x-component of the wave vector in the layer a(b) and ({c}_{{text{a}}(b)}=sqrt{{lambda }_{{text{a}}(b)}/{rho }_{{text{a}}(b)}}) is the velocity of longitudinal sound. To simplify calculations the transverse mode in the solid layer was ignored in derivation of the EPs Eq. (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 results19, ρx = (aρa + bρb)/d, ({rho }_{y}^{-1}=(a{rho }_{,text{a}}^{-1}+b{rho }_{text{b},}^{-1})/d), and ({B}^{-1}=(a{lambda }_{,text{a}}^{-1}+b{lambda }_{text{b},}^{-1})/d). Frequency-dependent components ρx and ρy of the effective density tensor and the effective compliance are shown in Supplementary Figs. 1 and 2.

### Analysis of isofrequency curves

The dispersion relation

$${omega }^{2}/B={k}_{x}^{2}/{rho }_{x}+{k}_{y}^{2}/{rho }_{y}$$

(3)

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 kx varying within the first half of the Brillouin zone and conserved ky = 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 ∂ω/∂kx < 0. It starts at frequency f = 106.8 kHz. When the angle γ = 0 the band edge is red-shifted to the frequency f0 = 104.8 kHz (see Supplementary Fig. 3).

Near the frequency f0 the dispersion relation, f = f(kxky) 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 f0 give two sets of isofrequency curves topologically similar to hyperbolas. These two sets are separated by the isofrequency contour f(kxky) = f0 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, ky = 0, the band gap occupies the interval f1 < f < f2, where the gap edges f1 = 46.7 kHz and f2 = 77.5 kHz are reached at the edge of the Brillouin zone, kx = π/d = 314 m−1. The isofrequency curves in the first passing band (f < f1) are closed since equation f = f(kxky) has real solutions for any direction of propagation. For the frequencies slightly above f1 a narrow interval of the directions of propagation almost parallel to the axis, ky << kx becomes inaccessible. This interval increases with f. Thus, the topological transition from closed to opened isofrequency curves occurs at f = f1. Evolution of isofrequency curves with frequency is shown in Supplementary Figs. 46.

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 kxd ≈ 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 ky, Eq. (3) becomes a transcendental equation. Unlike canonical hyperbolas, both blue curves contain regions with positive, negative curvatures, and a flex point. Blue solid vertical line in Fig. 2 marks the position of the exact flex point at ky = 0.89 cm−1. It corresponds to the direction of the Bloch vector k that makes an angle (gamma =arctan ({k}_{y}/{k}_{x})approx 6{0}^{circ }) with the negative direction of axis x. Negative x-projections of vector k and phase velocity ({V}_{,text{ph,ref},}^{(x)}) are due to anomalous dispersion. The product ({V}_{,text{ph,ref},}^{(x)}{V}_{,text{g,ref},}^{(x)} , < , 0) where the group velocity is always positive. 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 structures16.

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 ky ≈ 71 m−1 (see Fig. 2). At this point the dynamic density ρx ≈ 32 kg/m3 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 ky-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 Vg << Vph, which is related to strong scattering at the interfaces20.

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 theta /sin gamma ={c}_{text{a}}/{V}_{text{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 ({{bf{V}}}_{text{g,inc}}cdot {{bf{V}}}_{text{ph,inc}}={c}_{text{a},}^{2} , > , 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 xy plane practically does not change along its pass through a 100-period long lattice since for the selected parameters the diffraction coefficient (D={partial }^{2}{k}_{x}/partial {k}_{y}^{2}) vanishes. However, in the 3D map of sound intensity shown in Fig. 4, the width of the beam in the yz plane grows linearly because the corresponding diffraction coefficient ({partial }^{2}{k}_{x}/partial {k}_{z}^{2} , ne , 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. Fig. 3: Nonspreading transmission of acoustic Gaussian beam through 100 periods of the water steel layered structure. Fig. 4: 3D modeling of transmission of acoustic Gaussian beam through water steel layered structure.

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 refraction16

$${n}_{text{g}}=-frac{{k}_{y}{c}_{text{a}}/omega }{partial {k}_{x}/partial {k}_{y}}sqrt{1+{left(frac{partial {k}_{x}}{partial {k}_{y}}right)}^{2}}.$$

(4)

Since in Fig. 2 the derivative ∂kx/∂ky is negative at the flex point (ky = 0.89 cm−1 with ({V}_{,text{g},}^{(x)} , > , 0)), the beam exhibits positive refraction for group velocity, ng > 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 3{5}^{circ }approx 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, kx, ky, 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. Fig. 6: Variation of the phase for forward and backward propagation.

### 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 surfaces22,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 simulations22 show presence of lateral diffraction maxima. Lateral diffraction disappears if the corresponding isofrequency surface has a flat region11,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 crystal25 of steel rods in water, and more recently in perforated Plexiglas plates in air17. Simultaneous collimated propagation of acoustic and electromagnetic waves in phoxonic crystal of silicon rods in air was numerically predicted26. 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 one-dimensional 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.