Interface transmittance and interface waves in acoustic Willis media

Acoustics Willis media, known as bianisotropic acoustic media, incorporate additional coupling between pressure and velocity and between momentum and volumetric strain in their constitutive equation. The extra coupling terms have a significant influence on acoustic wave behavior. In this paper, the unusual wave phenomena relevant to interfaces between homogeneous acoustic Willis media are theoretically studied. We show that Willis media offer more flexible control in wave front and energy flow when waves are transmitted through an interface. Different from traditional acoustic fluid, Willis acoustic media support edge and interface waves, for which the existence conditions and corresponding wave features are systematically investigated. The study unveils more possibilities for manipulating acoustic waves and may inspire new functional designs with acoustic Willis metamaterials.


Introduction
In the past 20 years, with the emergence and development of metamaterials, the design space of wave devices and other functionality structures has been enlarged unprecedentedly. Metamaterials often exhibit abnormal material properties that natural materials usually do not have, which can lead to many novel wave phenomena, such as negative refraction [1,2], super lens [3,4], and wave cloaking [5][6][7], providing broad application prospects and meanwhile appealing more sophisticated homogenization for the characterization of the dynamic effective properties. In this background, the theory of Willis materials, initially proposed by Willis [8] in the 1980s for the dynamic behavior of solid composite materials, has regained much attention [9][10][11][12].
Acoustic Willis media (known as acoustic bianisotropic media) incorporate coupling between pressure and velocity and between momentum and volumetric strain. It has been found that the local Willis coupling is directly related to the local asymmetry of unit cells [12]. Accordingly, different designs of acoustic Willis meta-atoms have been proposed, such as the membrane unit [13,14], folded channel [15], and Helmholtz resonators [16,17]. The extra degree of design freedom offered by the Willis coupling is utilized to realize various novel wave functionalities. Several studies have observed asymmetric reflection [11, [18][19][20] when waves are incident from different directions, based on which the unidirectional absorber may be realized [21][22][23]. When Willis meta-atoms are used in metasurfaces [13,[24][25][26] or metagratings [15,27,28] for anomalous refraction or reflections, independent and more efficient control of transmission and reflection can be realized. In addition, active mechanisms can be introduced to enhance the significance and flexibility of the coupling effect [29,30], with which many non-reciprocal phenomena are demonstrated [30,31]. Although remarkable progress has been made in recent years, most of the research on acoustic Willis coupling concentrates on the physical origin and design of Willis meta-atoms. Extended wave functionalities are usually demonstrated in reduced dimensionality, such as metasurfaces. Relatively, systematic theoretical study of wave behaviors in continuous Willis media has not received much attention. In recent works, some phenomena, such as sound scattering [32], sound focusing [33], and the topological phase transition [34] in Willis acoustic media, have been investigated.
In this paper, the unusual wave phenomena relevant to interfaces between homogeneous acoustic Willis media are theoretically studied. Section 2 discusses the general properties of bulk waves in Willis media, such as slowness curves, wave modes, and impedance. Section 3 explains the interface transmittance when a wave is incident into a Willis medium, which exhibits more flexible control in wave front and energy flow through material parameters. Section 4 considers the edge and interface waves in Willis media, partly demonstrated in [35]. Here, we present a more systematical examination considering other possibilities along with the corresponding parameter conditions and wave modes. Finally, conclusions are drawn in Section 5.

General bulk wave properties
Assuming harmonic motion with circular frequency ω and time convention e −iωt , the momentum and continuity equations of Willis media are, respectively, as follows: where p is the acoustic pressure, ε the volumetric strain, v the particle velocity, and μ the momentum density. Distinct from traditional acoustic fluid, the constitutive relations of Willis media are characterized by where κ and ρ are the bulk modulus and mass density, respectively. The density is in tensorial form and can be anisotropic. The vector S represents the acoustic Willis coupling term, whose fundamental physics originates from the locally monopolar-dipolar coupling and non-local phase effects of acoustic scatterers [12]. Here, we assume that the non-local effects could be ignored; hence, S is purely imaginary. Combining Eqs 1, 2 gives the wave equation of Willis acoustic media: The wave equation is quite complex, and to simplify the analysis and highlight the Willis coupling effect, we consider in this section the isotropic density ρ, i.e., ρ ρI. According to [32], a dimensionless real vector W iS/ κρ √ is introduced for the Willis coupling for the sake of convenience. For a plane wave p p exp(ik · r) withp and k being the complex amplitude and the wave vector, respectively, the dispersion equation is where W = |W| and k = |k|. Eq. 4 has been formulated by [32]. In a two-dimensional (2D) case, the scenario of involved directions regarding the bulk wave propagation is depicted in Figure 1A, where ϕ is the azimuthal angle of W and ϕ′ is the angle between W and k. Eq. 4 clearly shows that the slowness curve is an ellipse because of W. This is a natural consequence because the coupling vector brings directionality. Without loss of generality, align the x-axis with W, then Eq. 4 is simplified as ( 5 ) Figure 1C shows the slowness curves of Willis media with different values of W. The major axis of the ellipse is collinear with the direction of W, and the ellipse tends to be flatter as W becomes larger. Compared with ordinary medium (W = 0), the phase velocity parallel to W remains unchanged, whereas the phase velocity perpendicular to W increases.
To further characterize the particle movements in Willis media, we investigate the velocity field of plane waves. Expressing the velocity field in the medium as v v exp(ik · r), the complex amplitudev can be derived aŝ where the velocity components are decomposed in directions parallel to and perpendicular to the wave vector, represented by subscripts " " and "⊥," respectively. Eq. 6 reveals that waves in the Willis acoustic medium are generally elliptically polarized, and because |v ⊥ | ≥ |v |, the long axis of the ellipse is collinear with the wave vector direction. This is a typical feature different from the traditional acoustic media that only supports the longitudinal wave. A pure longitudinal wave (v ⊥ 0) only happens when the coupling vector is collinear with the wave vector (ϕ′ 0 or π). A typical wave picture is shown in Figure 1B, where the grid points of solid lines represent the real-time positions of particles, different colors represent their phases, and the dotted grid points represent their initial positions. For a wave propagating to the right, particles rotate clockwise around their equilibrium position in an elliptical orbit.
In the aforementioned discussion, ρ and κ are assumed to have positive values. As the Willis medium is often used to characterize metamaterials, it is reasonable to allow the negative values as well. If one of ρ and κ is negative, the definition of W should be modified as W iS/ −κρ √ . As a result, Eq. 5 changes to In Eq. 7, if W ≤ 1, the medium does not support any traveling waves as no real solution of k exists. However, if W > 1, it is possible Frontiers in Physics frontiersin.org to find a real solution of k. In that case, the media have hyperbolic slowness curves, as shown in Figure 1D. The real axis of the hyperbola is perpendicular to W, and its eccentricity increases with W. The group velocity can be determined by calculating the timeaveraged intensity of power flow I Re(p † v)/2. For different sign combinations of ρ and κ, from Eqs 5, 7, respectively, they are expressed as The corresponding directions of group velocity are also plotted in Figures 1C,D for elliptic and hyperbolic slowness curves. For a medium with ρ < 0 and κ > 0, the group velocity points to the outer normal of the hyperbola, whereas for a medium with ρ > 0 and κ < 0, the group velocity points to the inner normal, as depicted in Figure 1D. If both ρ and κ are negative, Eq. 5 remains unchanged, so the slowness curve is elliptic. However, v changes in the opposite direction, as well as I. Hence, the group velocity points to the inner normal of the ellipse, as depicted in Figure 1C. The directions of group velocity can also be derived from the gradient of Eq. 5 or Eq. 7, but the causality constraint must be considered as in [36]. Further analysis shows that allowing the density to be anisotropic only changes the eccentricity of the ellipse or hyperbola (see Supplementary Material S1).

Interface transmittance and abnormal refraction
Having been acquainted with the wave properties, we consider in this section the transmittance of acoustic waves through the interface between an ordinary acoustic medium (Medium I) and a Willis medium (Medium II), as presented in Figure 2A. A plane wave is incident from the left side, and the incident, reflection, and refraction angles are θ i , θ r , and θ t , respectively. The parameters on both sides are marked in the figure.
The reflected wave is on the side of the ordinary medium, and the reflection angle follows θ r θ i . On the other hand, the refraction angle in the Willis medium can be determined by the continuity of the tangential wave vector on the interface (k II sin θ t k I sin θ i ) as well as Eq. 4, where k I and k II are wave numbers on both sides, respectively. The refraction angle θ t is ruled by As it is a transcendental equation, in general, θ t can only be calculated numerically. However, for normal incidence (θ i = 0), it is obvious that θ t = 0, regardless of the magnitude and direction of the Willis coupling vector. In this case, the impedance of the Willis medium in the normal direction of the interface is Z n p/v . Substituting Eq. 6 and noting that ϕ′ = ϕ, we get Therein, "±" represents the different signs when propagating to the positive or negative direction of the x-axis. In comparison to traditional acoustic fluid, an extra factor ( 1 + W 2 sin 2 ϕ ± iW cos ϕ) is added to the impedance. As the absolute value of this factor is greater than 1, the presence of Willis coupling will always increase the impedance, regardless of its azimuthal angle. The impedance matching condition at the interface is which indicates that the direction of the coupling vector must be perpendicular to the wave vector if the full transmission is required. We learn in Section 2 that the direction of energy transmission may differ from the wave vector in the Willis medium. Even for the normal incidence, the energy flow direction may still deviate. Using Eq. 4, the reflection power (R |I r |/|I i |) and the two components of the transmitted power flux (T x I tx /|I i |; T y I ty /|I i |) are calculated. Figure 2B presents the variation of transmission and reflection power versus the azimuthal angle ϕ under normal incidence. Other parameters are set as ρ = ρ 0 , κ = 0.5κ 0 (ρ 0 = 1.2 kg/m 3 , κ 0 = 1.4 × 10 5 Pa, here and after), and W = 1 so that when ϕ = 90°or ϕ = 270°, the impedance matching conditions are satisfied. The calculated frequency is 7,000 Hz. Figure 2B shows that at these two points, the full transmission happens. When W points to other directions, there will be energy flow in the y direction and reflection. In Figure 2C, the deflected energy propagation for the case of Willis coupling with ϕ = 60°( as indicated by the dot in Figure 2B) is verified by finite element method (FEM) simulation considering the normal Gaussian beam incidence. All FEM simulations in this paper study are carried out via COMSOL software. The energy flux component in the y-direction corresponding to ϕ = 60°is negative. Correspondingly, it is observed that the wave beam deflects downward after passing through the interface.
For oblique incidence and considering negative parameters, more interesting phenomena can be obtained, as exemplified in Figure 3, wherein panels in the first row are the refractive patterns drawn from the slowness curve analysis, whereas the second row presents the results of FEM simulations of wave beams with 7,000 Hz. As shown in Figures 3A,D, we use a Willis medium with ρ > 0 and κ < 0, whereas in Figures 3B,E, we use a Willis medium with ρ < 0 and κ > 0. The used material parameters are given in Figures 3D,E. Because of the hyperbolic slowness curve, the wave number k y of the incident wave cannot be less than a certain critical value to get a real wave number for the refraction wave, which is opposite to the case of positive parameters. Taking the 60°incident angle as an example, as indicated by the k i arrow, the refractive k r is determined from Figures 3A,B. The group velocity must have a positive x component to ensure energy always goes forward. For κ < 0 cases, as the group velocity points to the inner normal of the hyperbolic, negative refraction of phase velocity and positive refraction of group velocity are predicted, as shown in Figure 3A. Conversely, for the ρ < 0 case, as shown in Figure 3B, group velocity points to the outer normal of the hyperbolic. Thus, positive refraction of phase velocity and negative refraction of energy flow are predicted. From the FEM simulations in Figures 3D,E, the aforementioned analysis is confirmed by observing the refracted wave beam and wave front. As shown in Figures  3C,F, we use a Willis medium with ρ < 0 and κ < 0. The typical negative refraction is realized like the double ordinary doubly negative medium.

Edge and interface waves
When boundary condition is considered in solving the wave equation, it is possible to find surface modes, for example, the Frontiers in Physics frontiersin.org well-known Rayleigh surface wave for solids. For an ordinary acoustic fluid, surface modes are not supported. However, in the realm of acoustic metamaterials, surface modes can be found in acoustic media with negative parameters [3,37] or with gyrotropic mass [38]. In [35], we have partly demonstrated the existence of interface waves at the interface of two Willis media. Here, we present a more systematic examination of the edge and interface waves of Willis media.

Edge waves with sound hard boundary
As sketched in Figure 4A, we consider a semi-infinite of Willis acoustic medium, and the Cartesian coordinate system is established to make the open edge along the y-axis. Considering acoustic field explicitly expressed by p p exp(ik x x + ik y y), for a possible edge mode, the field is traveling along the y-direction and attenuated away from the surface. Thus, k y must be real, and k x must have a positive (negative) imaginary part when the surface is on the left (right).
For sound soft (free) boundary, the requisite boundary condition is p (x = 0) = 0, which holds only forp = 0. Then, the velocityv vanishes at the same time, and it is concluded that no edge mode is supported on the sound soft boundary. For sound hard boundary, it is required that the normal velocity on the boundary vanishes, that is, v x (x = 0) = 0. Combining the boundary condition and the requirements of the wave vectors, we can conclude that surface waves may exist on the hard boundaries, and the parameter conditions for their existence are (Supplementary Material S2) ρ yy κ > 0 and i S x ρ yy − S y ρ xy > 0, when the surf ace is on the lef t, ρ yy κ > 0 and i S x ρ yy − S y ρ xy < 0, when the surf ace is on the right.
(12)  The wave number along the surface is k y ± ω ρ yy /κ . For the case of isotropic density and ρ > 0, κ > 0, the aforementioned condition reduces to be simply W x > 0, when the surf ace is on the lef t, W x < 0, when the surf ace is on the right. (13) Eq. 13 has an obvious geometric meaning that the vector W has to point away from the edge to find the edge wave. In this case, the wave mode can be analytically expressed as In Eq. 14, k x is purely imaginary only when W is perpendicular to the boundary (W y = 0), and the wave is exponentially attenuated away from the boundary. In other cases, k x has both real and imaginary parts, which means that the edge wave is oscillatory attenuated away from the boundary (see Supplementary Material S3). In addition, non-zero W y will make the imaginary part of k x smaller. Hence, the attenuation will slow down further. It is also noticed that, unlike the bulk wave, the edge wave is linearly polarized, and particle velocity possesses only components along the interface. The wave pattern of the edge mode is shown in Figure 4B. The grid points on the solid lines represent the real-time position of particles, and the colors represent their phases. The dotted lines represent their initial positions.
FEM simulations are performed to verify the surface wave on the hard boundaries, as shown in Figure 5. Material parameters of the Willis medium are ρ = ρ 0 , κ = κ 0 , and W = 1. Two directions of the coupling vector, ϕ = 0°( Figure 5A) and ϕ = 180°( Figure 5B), are used here. The calculated frequency is 5,000 Hz. A pair of point sources with a half wavelength distance and opposite phases is set on the hard boundary to form a dipole, which can stimulate the surface mode more efficiently. In Figure 5A, the condition that W points away from the surface is satisfied. Correspondingly, the edge wave is observed along the y-axis. When W points to the surface, as in Figure 5B, no surface mode exists, and only the bulk mode is excited. Figure 5A shows that, for large W perpendicular to the edge, the wave vectors of edge mode and bulk mode differ much, so they are not easily coupled with each other. For edge waveguides having corners not so sharp, the edge wave can pass through without obvious scattering into the bulk, showing some robustness, as depicted in Figure 5C.
To design a waveguide with more robustness, we can utilize a medium supporting edge wave that does not allow bulk waves. A simple choice is to use a single negative medium with isotropic Frontiers in Physics frontiersin.org density. However, such a medium does not possess edge mode because Eq. 12 cannot be satisfied, so we seek the possibility in the Willis medium with anisotropic density. The derived condition for the Willis medium that does not support bulk waves is (see If a set of material parameters can satisfy Eq. 12 and 15 at the same time, the edge wave transmission would be very stable. In the example in Figure 5D, we choose a set of parameters as ρ xx = -ρ 0 /2, ρ yy = ρ 0 , ρ xy = 0, κ = κ 0 , S x −i κ 0 ρ 0 √ , and S y = 0 so that Eqs 12, 15 are simultaneously met for the edge on the left. All energy is concentrated at the boundary without any bulk waves. Moreover, in this exclusively edge-mode medium, a single monopole source is enough to excite the edge wave without matching its mode. In Figure 5E, a waveguide containing corners with the right angle is established using two different Willis media to meet the corresponding parameter requirements on the boundary of different orientations. Therein, Medium I is the same as that in Figure 5D, and Medium II is simply obtained by rotation from Medium I. The edge wave transmits through the designed route without energy leaking into the bulk.

Interface waves between two media
In practice, if the impedance of the Willis and adjacent media significantly differs, their interface can be treated as an ideal soft or hard boundary. For other cases, further analysis is required to estimate whether interface modes exist.
Consider the problem shown in Figure 6A. The ordinary acoustic medium (Medium I) on the left side and the Willis medium (Medium II) on the right side form an interface along the y-axis. In order to simplify the problem, the discussion is limited to Willis media with isotropic density and ρ > 0, κ > 0. The density and bulk modulus of the ordinary medium are ρ 0 and κ 0 , respectively. The parameters of the Willis medium are ρ, κ, and W. For an interface wave, the pressure fields on the two sides can be written as p I p exp(ik I x x + ik s y) and p II p exp(ik II x x + ik s y), respectively, where k s is the tangential wave vector along the interface. It must be continuous and real on both sides. k I x and k II x are normal wave vectors, which should satisfy Im(k I x ) < 0 and Im(k II x ) > 0 to ensure attenuation in the direction away from the interface. The conditions for the existence of interfacial waves are as follows (see Supplementary Material S5): where ρ ρ/ρ 0 and κ κ/κ 0 . The interface mode is difficult to be expressed analytically, so we demonstrate it with a specific numerical example in Figure 6B. The parameters are set as ρ 0.2, κ 0.1, W = 1, and ϕ = 0°. The grid points in Figure 6B represent the real-time positions of the particle. Due to the non-zero normal velocity, particles near the interface rotate in an elliptical trajectory around their equilibrium positions, which is very different from the edge mode. At the interface, only the normal velocity maintains continuity. The two media have a relative slip in the tangential direction. The FEM simulation of the interface wave is shown in Figure 7A with the same parameters. The influence of the material parameters on the robustness of interface waves is analyzed in Supplementary Materials S6.
In the problem setup shown in Figure 6A, if the two domains are both Willis media, the analysis of the interface wave is more complex, so the general situation is not studied. However, if the two media meet some specific conditions, the problem will become intuitive. In Section 4.1, we use the boundary condition that the normal velocity is zero to derive the existence condition Eq. 12 for edge waves and the corresponding edge wave number k y ± ω ρ yy /κ . The derivation is reversible; that is, if k y ± ω ρ yy /κ and Eq. 12 holds, the normal velocity on the boundary must be zero. Thus, we conclude that, in the current problem of two Willis domains, if the tangential wave numbers on both sides of the interface are the same and Eq. 12 holds for each side, the continuity condition of the interface will always be met. In this situation, the media on both sides are mutually sound hard FIGURE 6 (A) Schematic representation of the interface mode between an ordinary medium (Medium I) and a Willis medium (Medium II). The interface mode is a traveling mode along the surface and evanescent away from the interface. (B) Wave mode of the interface wave, where the interface expands at a finite width to clearly display the particle motion on both sides. The insets show typical particle orbits on each side.
Frontiers in Physics frontiersin.org boundaries to each other, so the interface wave can be supported. On the premise that ρ I yy /κ I ρ II yy /κ II , an existence condition for the interface wave reads An FEM example is shown in Figure 7B. For convenience, the parameters of Medium II are obtained by the mirror symmetry operation of Medium I about the y-axis, so if Eq. 17 holds on one side, it also holds on the other side. The parameters are ρ I = ρ II (components are marked in Figure 7B), κ I = κ II = κ 0 , S I  Figure 7B satisfy Eq. 15, so a pure interface wave is excited without bulk leaking.
When the density is isotropic, ρ > 0, and κ > 0, Eq. 17 reduces to Similarly, Eq. 18 implies that vectors W on both sides point away from the interface. An example is shown in Figure 7C. The parameters in Medium I are ρ I = 5ρ 0 , κ I = 5κ 0 , W I = 1, and ϕ I = 180°, and those in Medium II are ρ II = ρ 0 , κ II = κ 0 , W II = 1, and ϕ II = 0°. The density and bulk modulus of the Willis medium on the left side are both five times those on the right side, so the tangential wave vector is continuous on the interface. The coupling vector point away from the interface on both sides. An interface wave is observed as Eq. 18 is satisfied.

Conclusion
In this study, unusual wave phenomena relevant to interfaces in the homogeneous acoustic Willis media are studied theoretically. We show that the media exhibit anisotropic features due to the Willis coupling vector terms, and the slowness curve can be tuned between elliptical and hyperbolic shapes. The interface transmittance can be adjusted by the magnitude and direction of the coupling vector, which offers more flexible control in the compared with the traditional acoustic fluids. The Willis acoustic media support edge waves at acoustic hard boundaries and interface waves at interfaces between an ordinary acoustic fluid and a Willis medium or between two Willis media. Particularly, the edge modes may also exist in certain Willis media that do not support bulk modes, in which case they can achieve high transmission.
The study unveils more possibilities for manipulating acoustic waves and may inspire new functional designs with acoustic Willis metamaterials. It should be noted that this theoretical study assumes continuous acoustic media already with Willis coupling. On the experimental side, designing metamaterials with the wanted coupling effect is still an ongoing and challenging task. Especially for those wave phenomena calling for a strong S vector and negative density or modulus, the experimental demonstration may necessitate a more sophisticated design.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.