Willis coupling in water waves

For mechanical waves, Willis coupling means a cross-coupling between stress and velocity or between momentum and strain. In contrary to its realization in acoustic and elastic waves, whether Willis coupling exists in water waves, as another kind of mechanical wave, is still unknown. Here, we propose and establish the concept of Willis coupling in water waves as the cross-coupling between the horizontal velocity at the free surface and the acceleration potential or between the vertical displacement of the water free surface and the flux. Thanks to the surface wave feature of water waves, the proposed metamaterial’s resonating conditions can be tuned by using the wave shoaling effect. With a proper three-dimensional design, Willis coupling can still have significant effects with resonance in the long-wavelength regime. Furthermore, by adding loss in the Willis metamaterial, asymmetric reflectance and absorption can be achieved, which are useful for applications such as seismic isolation, coastal protection, water-wave energy-harnessing, and also for constructing non-Hermitian exceptional points.

The concept of metamaterials has also been successfully implemented for water waves to achieve various unique material parameters and phenomena, including large anisotropy, negative index, negative gravity, and invisibility cloaking [18][19][20][21][22][23][24]. With the background of the recent demonstrations of Willis coupling in propagates along x-direction on the water surface with height h 0 from a rigid bottom. It consists of mushroom-shaped pillars (green color, radius r o and thickness d 1 ), and crescent-shaped walls (gray color, with water depth d 2 below water surface) supported by narrow poles of radius r p and a height of h1. The crescent shape is generated by subtracting a circular disk of radius r i from a larger disk of radius r o with displacement Δ. The metamaterial is periodic in both x-and y-directions with lattice constants a = 0.6 m and w = 0.6 m, respectively. acoustics and elastic flexural waves, it becomes natural to ask whether Willis coupling also exists for water waves, which are indeed a kind of mechanical waves, like acoustic waves and elastic waves in a solid, for completeness.
In this paper, we propose a three-dimensional (3D) Willis metamaterial with unit cell composed of a crescent pedal supported by a thin pillar, whose resonating conditions in different directions can be tuned by using the wave shoaling effect of water waves [53]. By breaking mirror symmetry, the metamaterial can induce Willis coupling and provide different wave impedances for forward and backward incidence. The existence of Willis coupling in water waves is theoretically demonstrated by an effective medium model and numerical simulations of the 3D metamaterial. Here, the Willis coupling can occur even at very small frequencies by resonance effect from wave shoaling, which may not be straightforward in acoustic or elastic waves [50][51][52]. Although metamaterials for water waves are far less developed at the moment than their electromagnetic, acoustic, and mechanical counterparts, their investigations as we propose here may find potential applications in coastal protection, damping, and also energy harvesting [54,55].

Willis coupling in water waves from wave shoaling
Here, we consider a 3D water-wave metamaterial consisting of a square lattice of unit cells (a × w) with a water depth of h 0 above a flat bottom, as shown in figure 1. Each unit cell consists of a mushroom-shaped pillar (green color) with radius r 0 , thickness d 1 and a water depth of d 2 below the water surface and is supported by a narrow pole of radius r p and a height of h 1 . Each pillar has a crescent-shaped wall (gray/white color below/above the water surface) on top of it, which is geometrically generated by the difference between two circular disks of radii r o and r i respectively with a relative shift of Δ. The water wave propagates on the water free surface (z = 0) along the x-direction, passing through several layers of these unit cells. The wave shoaling in the region above each pillar has an abrupt decrease in water depth, causing a decrease of phase velocity and induces resonance. Here, we set h 0 = 0. The water wave, described using a potential flow approach [21,56] for an incompressible and irrotational fluid with viscosity being neglected in this section [57], has its velocity potential Φ (velocity v = −∇Φ) in the bulk region and satisfies the Laplace equation: together with a boundary condition on the water free surface at z = 0: and zero normal gradient of Φ (i.e. no flow condition) on all the rigid walls, including the bottom at z = −h 0 and the surfaces of the mushroom-shaped pillars. The standard acceleration g caused by gravity is taken as 9.8 m s −2 .
In the background region without the structures (labeled with subscript '0'), equations (1) and (2) can be reduced by using separation of variables between x and z as a wave equation for water waves propagating along the x-direction (where the time variable has been suppressed by time-dependence exp (−iωt) with ω the angular frequency in units of rad s −1 ): The acceleration potentialφ on the water surface and the flux q in the wave equation are defined withh 0 being the ratio between the flux and the horizontal surface velocity in the background region.h 0 varies with the wavenumber k 0 bỹ h 0 = tanh (k 0 h 0 ) /k 0 , giving rise to the frequency dispersion ω 2 /k 2 0 = c 2 0 = gh 0 with wave speed c 0 decreases with an increasing frequency. Now, we consider the region with the structures. The crescent-shape walls are designed to break the mirror symmetry along the propagation direction and induce a Willis coupling. Then, the water wave equation, the homogenized version of equation (3), has to be written in a matrix form as: where ζ is the vertical displacement of the free water surface, v x is the horizontal velocity at the water surface, g eff is the effective gravity and h eff is the effective height of the metamaterial's section. A non-zero τ eff is unavoidably introduced in the 2 × 2 constitutive matrix to represent the Willis coupling arising from a mirror symmetry breaking along x. Equation (4) approaches to equation (3) when g eff = g, h eff =h 0 and, τ eff = 0 and the metamaterial's parameters go back to those of the background medium. We have assumed reciprocity to enforce a zero sum for the two off-diagonal elements [58]. If the system is lossless, the constitutive matrix is Hermitian, and τ eff becomes purely real. On the other hand, if the system does not break mirror symmetry, we have equation (4) preserved by q → −q and x → −x so that τ eff should be zero as expected and the surface displacement ζ of the free surface is solely linked to the acceleration potentialφ through gravity (Bernoulli equation on elevated water surface), while the horizontal velocity v x at the free surface is solely linked to the flux q through an effective height of water. The Willis coupling for water waves means that the vertical surface displacement ζ is related not only to acceleration potential but also to the flux while the horizontal surface velocity v x is now not only related to the flux but also to the acceleration potential.
A direct consequence of a non-zero Willis coupling is considered by solving equation (4) for the eigenmodes to obtain the dispersion relationship of the metamaterial: which now contains the Willis coupling parameter τ eff that modifies the dispersion relation of the background medium. The sign of the τ eff depends on the direction of the crescent-shaped wall (and is revealed by the difference between the forward and backward reflectance in the later part). The wavenumber k, although modified, is the same for both the forward and backward modes when mirror symmetry is broken. This comes from the time-reversal symmetry of the lossless system instead. On the other hand, the sign of τ eff and its magnitude will affect the impedances defined here by Z ± ±φ/q, for the forward (+) and backward (−) propagating or decaying modes. They become different from each other for a non-zero Willis coupling [46] and are governed by: The Willis coupling parameter τ eff affects the expression of k and introduces an additional term in Z ± . The dispersion relationship together with the impedances enable us to identify the three constitutive parameters (effective medium) as These formulas generalize the conventional case 1/g eff = k/ (ωZ) and 1/h eff = kZ/ω when there is no Willis coupling (Z + = Z − = Z and τ eff = 0). On the other hand, we can perform 3D full-wave simulations, consisting of the microstructures (COMSOL multiphysics in solving equations (1) and (2), with two distinct excitations indexed by j = 1, 2. Deep inside the metamaterial and across a single unit cell from x = 0 to x = a, we can obtain theφ and q fields on the two sides of a unit cell. Then, a discrete-difference approximation of equation (4), i.e. in the local limit |k| a π, becomes where the discrete difference Δφ (also a similar expression for Δq) can be obtained by and the macroscopic average φ can be intuitively obtained by averaging over the boundary of the unit cell as φ = 1 2w dy φ a, y +φ 0, y , with a and w being the periodicities of the metamaterial in the x-and y-direction. To get higher accuracy in compensating the numerical dispersion arising from the discrete-difference approach, the field averaging formula (10) can be modified as which is used for the effective medium extractions in this work. Z 0 is the impedance in the background, with its expression introduced later. It is derived by assuming the secondary radiation comes from point dipole moments at the center of each unit cell (see supplementary note 1 for derivation (https://stacks.iop. org/NJP/23/073004/mmedia)). We note that equation (11) returns to equation (10) in the limit of small k 0 a, in which the numerical dispersion becomes smaller. Similar treatment holds for the flux q in obtaining the macroscopic field q . As a result, we have developed a scheme based on an integral approach in obtaining these macroscopic fields in order to achieve a better accuracy at a finite frequency. Equation (8) then allows us to extract the constitutive matrix numerically by where subscript j (=1, 2) stands for the jth excitation. Here, we excite the metamaterial with 3 unit cells (beyond which the results converge in our current case) and extract the constitutive matrix using equation (12) from the fields of the central unit cell by generating an incident wave from the left-(j = 1) and right-(j = 2) hand sides of the whole metamaterial. We note that for our design, the resonance is actually very localized and the number of unit cells is not critical for the effective medium extraction. The numerically extracted effective medium is shown in dimensionless units in figures 2(a)-(c) for g/g eff ,h 0 /h eff and τ eff c 0 respectively. In the long-wavelength limit, the effective medium parameters start from g/g eff ∼ = 1,h 0 /h eff ∼ = 1 (because the flow is nearly unaltered without blocking) and τ eff = 0. Their magnitude increases significantly when near the resonance frequency around 0.28 Hz, which is caused by the wave shoaling of water waves. The Willis coupling significantly deviates from zero (we will see its effect in reflectance later). The dimensionless impedances Z ± /Z 0 for the forward and backward modes evaluated from the effective medium perspective are shown in figure 2(d), where Z 0 is the impedance of the background medium defined as Z 0 = g/h 0 . The Willis coupling now introduces an imaginary part of the two impedances with the same magnitude but opposite signs in the passbands (equation (6) with real k, g eff and τ eff ). In the bandgap, k becomes purely imaginary and the real part of the two impedances becomes zero. The band structure evaluated from the effective medium parameters (g eff , h eff and τ eff ) with equation (5) is shown as solid lines in figure 2(e), which clearly exhibits the signature of a resonance gap caused by the wave shoaling of water waves at around 0.28 Hz (with gap lower bound at Brillouin zone edge and gap upper bound at the centre). For comparison, the simulated band structure from 3D-full wave simulations is also plotted as red symbols in the same figure, and it can be observed that the band structures from effective medium approximate with great accuracy the simulated ones. This demonstrates the validity of this effective medium approach not only for the acoustic band in the low-frequency regime, but also for the band above the stopband, which is reminiscent of high-frequency homogenization theories that go beyond the usual quasistatic limit [59,60]. It is worthy of note that from the solved eigenmodes for the Willis metamaterial in the forward (with Bloch factor exp (ink 0 a)) and backward (with Bloch factor exp (−ink 0 a)) cases, the impedance on one side of the unit cell for each eigenmode can also be evaluated by Z ± = ± dyφ / dy q , with the solved impedance and band structure, the effective medium parameters can also be evaluated with equation (7) [13]. Compared to the excitation method we adopt in this work, the eigenmode method gives an accurate k but can only be used within a passband.  It can be seen that the transmission and reflection spectra reproduced from the same effective medium (lines) agree well with the simulated results (symbols) with different number of unit cells. It shows the validity of the effective medium description as the wavelength in the background medium (e.g. 3.7 m at 0.5 Hz) is larger than six times of the lattice constant.

Asymmetric reflection phase from Willis coupling
Here, |t f | = |t b | and |r f | = |r b | (subscript 'f '/'b' for forward/backward incidence) for a reciprocal and lossless system. The non-zero Willis coupling induces a difference in reflection phases in the two directions, as shown in figure 3(c). The relative phase arg r f /r b for both cases with 1 and 4 unit cells from full-wave simulations are shown in figure 3(c) as symbols while the result for 1 unit cell generated from the effective medium is shown as the blue curve. They all overlap with each other except at individual wavenumbers where |r| becomes nearly zero with ambiguous phase, showing that this asymmetry of reflection phase is the property of individual unit cells. For the current system, we have some flexibility to define the boundary of the unit cell by a small shift ξ in the x-direction, which will induce a phase shift on arg(r f /r b ) by −4ξk 0 . Figure 3(c) shows the linear fit (black dashed line) of the simulated relative phase, but it will not be able to renormalize the resonating asymmetry of reflection phases by such a shift of unit cell boundaries. It indicates the Willis coupling is genuine and cannot be eliminated by the convention in defining the boundaries of a unit cell. The deviation is largest around the resonance, where the Willis coupling reaches its maximum. We note that arg(r f /r b ) flips sign when we apply a mirror operation in the x-direction for the structure, which is also equivalent to a sign flipping on τ eff .

Asymmetric reflection amplitudes from Willis coupling in lossy cases
When a material loss exists in the Willis metamaterial, both the reflection amplitude and phase can become asymmetric for the two incident directions. Here, we would like the loss to be realized by viscosity to induce friction from the boundary layer above the mushroom-shaped pillars. In the original design, the boundary layer is much thinner than water height above the pillars and the viscous effect is negligible. In this section, we scale down the metamaterial to have significant effect of viscosity. Figure 4(a) plots the transmission (green dotted line), forward reflection (blue solid line), and backward reflection (red dashed line) amplitudes evaluated from 3D simulations using the linearized Navier-Stokes equations (COMSOL multiphysics). To include the effect of viscosity, we scale down the structures by 10 times. In this case, the viscous loss from the boundary layers leads to the asymmetric reflection amplitudes. |r f | is generally lower than |r b | as the water wave enters the top of the mushroom-shaped pillars with significant viscous loss in the forward direction while the water wave is blocked (by crescent walls) in the backward direction. The Willis coupling term becomes imaginary around resonance (shown in supplementary note 2). In the simulations, the boundary condition on the water surface is a prescribed normal stress of −ρ 0 gζẑ at z = 0, no-slip boundary condition is applied on the rigid bottom and all the surfaces of the metamaterial, while continuity periodic boundary condition is applied along the y-direction. The material parameters for water is taken as density ρ 0 = 1000 kg m −3 , sound speed 1500 m s −1 , dynamic viscosity 0.001 Pa s and bulk viscosity 0.0028 Pa s at room temperature. We note that the previous potential flow approach approximates well the simulations in linearized Navier-Stokes equations at zero viscosity. The introduction of no-slip boundary condition and a finite sound speed (comparing to the water wave velocity) does not alter the numerical spectra significantly at zero viscosity.
To enhance the viscous effect, we then scale down the metamaterial in figure 3 by 28 times, the corresponding results are shown in figure 4(b). The asymmetry in reflection amplitudes becomes larger while |r f | drops to nearly zero at 0.8 Hz. It indicates the phenomenon of unidirectional zero reflection, equivalently a non-Hermitian exceptional point (EP) with both eigenvalues and eigenvectors degenerate for the scattering matrix, and also for the constitutive matrix [61]. Such a Willis coupling induced EP occurs in the deep subwavelength regime with λ/a ∼ 21. The low frequency EP is possible due to several factors. First, the resonance can occur at low frequency due to the wave shoaling effect. The local wave speed is very small above the pillars, comparing to the background medium, equivalent to resonance from high refractive indices. Second, the normalized impedance of our metamaterial near zero frequency is not far from value one because the water flow is nearly unaltered without blocking as we are using narrow supporting poles. A low frequency resonating Willis coupling then splits the impedances (in forward and backward directions) and easily match one of them to have its real part being 1. Finally, with an appropriate absorption (from scaling in the current case), the imaginary part of the impedance can also be tuned to zero so that the impedance in the forward direction is now exactly one to get |r f | = 0. Such a low-frequency EP may not be straightforward for other kinds of Willis coupling designs in acoustic and elastic waves [50][51][52]. Interestingly, the asymmetric reflection from Willis coupling also means asymmetric absorption, which has been demonstrated for the case of the airborne acoustic sound [62,63].

Conclusion
In this work, we establish the concept of Willis coupling and demonstrate the existence of Willis coupling in linear surface water waves with both effective medium theory and designed 3D structures. The existence of Willis coupling in water waves is revealed through asymmetry in reflection phases and amplitudes (in the lossy case) in the forward and backward incidence. We further note Willis coupling, in terms of direction-dependent impedances, might be enhanced in nonlinear water wave configurations such as for ocean waves and may have important ramifications and applications such as in maritime engineering, seismic isolation, energy-harnessing.