Mehul P Makwana

Nicolas Laforge

Richard V Craster

Guillaume Dupont

Sébastien Guenneau

Vincent Laude

Muamer Kadic

Experimental observations of topologically guided water waves within non-hexagonal structures

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Considerable recent activity in wave phenomena is motivated through topological effects and focused on identifying situations where topological protection occurs that can enhance, or create, robust wave guidance along edges or interfaces. Remarkably, the core concepts that gave rise to topological insulators, originating within quantum mechanics 2 carry across, in part, to classical wave systems [START_REF] Khanikaev | Two-dimensional topological photonics[END_REF][START_REF] Ozawa | Topological Photonics[END_REF] . Topological insulators can be divided into two broad categories: those that preserve time-reversal symmetry (TRS), and those which break it. We concentrate upon the former due to the simplicity of their construction that solely requires passive elements. By leveraging the discrete valley degrees of freedom, arising from degenerate extrema in Fourier space, we are able to create robust symmetry-protected waveguides. These valley states are connected to the quantum valley-Hall effect and hence this research area has been named valleytronics [START_REF] Behnia | Polarized light boosts valleytronics[END_REF][START_REF] Schaibley | Valleytronics in 2D materials[END_REF] .

Hexagonal structures are the prime candidates for valleytronic devices as they exhibit symmetry induced Dirac cones at the high-symmetry points of the Brillouin zone (BZ); when perturbed these Dirac points can be gapped, leading to well-defined KK valleys distinguished, from each other, by their opposite chirality or pseudospin.

This pseudospin has been used in a wide variety of dispersionless wave settings to design valleytronic devices [START_REF] Dong | Valley photonic crystals for control of spin and topology[END_REF][START_REF] Xiao | Valley-Contrasting Physics in Graphene: Magnetic Moment and Topological Transport[END_REF] . Here we extend the earlier research by examining a highlydispersive physical system, i.e. water waves and move away from hexagonal structures. The topological protection afforded by these valley states is attributed to, both, the orthogonality of the pseudospins as well as the Fourier separation between the two valleys [START_REF] Makwana | Designing multidirectional energy splitters and topological valley supernetworks[END_REF] . The vast majority of the valleytrona) Electronic mail: mm107@ic.ac.uk b) Electronic mail: nicolas.laforge@femto-st.fr ics literature, inspired by graphene, opts to use hexagonal structures [START_REF] He | Acoustic topological insulator and robust one-way transport[END_REF][START_REF] Schomerus | Helical scattering and valleytronics in bilayer graphene[END_REF][START_REF] Ye | Observation of acoustic valley vortex states and valley-chirality locked beam splitting[END_REF][START_REF] Cheng | Robust reconfigurable electromagnetic pathways within a photonic topological insulator[END_REF][START_REF] Wu | Direct observation of valley-polarized topological edge states in designer surface plasmon crystals[END_REF][START_REF] Xia | Topological phononic insulator with robust pseudospin-dependent transport[END_REF][START_REF] Qiao | Electronic highways in bilayer graphene[END_REF][START_REF] Makwana | Geometrically navigating topological plate modes around gentle and sharp bends[END_REF][START_REF] Tang | Observations of symmetry induced topological mode steering in a reconfigurable elastic plate[END_REF][START_REF] Proctor | Manipulating Topological Valley Modes in Plasmonic Metasurfaces[END_REF] . However a negative that emerges with these, especially when dealing with complex topological domains [START_REF] Makwana | Designing multidirectional energy splitters and topological valley supernetworks[END_REF] , is that certain propagation directions are restricted due to mismatches in chirality between incoming and outgoing modes. Notably, this has led to hexagonal structures being prohibited from partitioning energy in more than two-directions [START_REF] Cheng | Robust reconfigurable electromagnetic pathways within a photonic topological insulator[END_REF][START_REF] Wu | Direct observation of valley-polarized topological edge states in designer surface plasmon crystals[END_REF][START_REF] Qiao | Electronic highways in bilayer graphene[END_REF] .

In this Letter, we demonstrate experimentally how a strategically designed square structure also allows for the emergence of valley-Hall edge states as well as allowing for the excitation of modes that are not ordinarily ignited within hexagonal valley-Hall structures. Additionally, the system chosen differs from the vast majority of the earlier literature [START_REF] He | Acoustic topological insulator and robust one-way transport[END_REF][START_REF] Schomerus | Helical scattering and valleytronics in bilayer graphene[END_REF][START_REF] Ye | Observation of acoustic valley vortex states and valley-chirality locked beam splitting[END_REF][START_REF] Cheng | Robust reconfigurable electromagnetic pathways within a photonic topological insulator[END_REF][START_REF] Wu | Direct observation of valley-polarized topological edge states in designer surface plasmon crystals[END_REF][START_REF] Xia | Topological phononic insulator with robust pseudospin-dependent transport[END_REF][START_REF] Qiao | Electronic highways in bilayer graphene[END_REF][START_REF] Makwana | Geometrically navigating topological plate modes around gentle and sharp bends[END_REF][START_REF] Tang | Observations of symmetry induced topological mode steering in a reconfigurable elastic plate[END_REF][START_REF] Proctor | Manipulating Topological Valley Modes in Plasmonic Metasurfaces[END_REF] that has focused on an idealised situation in which the dispersion of the host medium has been avoided. This assumption restricts the applicability of the earlier studies to a small subset of, potentially useful, physical platforms that could host topological effects. Most notably, this assumption does not hold for water wave systems, which generally support highly dispersive surface waves [START_REF] Lighthill | Waves in Fluids[END_REF] . The combination of topological physics applied to water waves is a relatively unexplored area 2,[START_REF] Wu | Topological helical edge states in water waves over a topographical bottom[END_REF] ; those that have conducted experiments have either focused on 1D systems [START_REF] Wu | Topological helical edge states in water waves over a topographical bottom[END_REF] or the hexagonal valley-Hall structure 2 . Potential applications of this budding area include controlling ocean wave energy [START_REF] Bennetts | Graded resonator arrays for spatial frequency separation and amplification of water waves[END_REF] , in a non-intrusive manner, for energy-harvesting or erosion mitigation [START_REF] Dupont | Type of dike using C-shaped vertical cylinders[END_REF] .

The fluid within our domain has a constant depth of h = 4 cm and contains a periodic array of rigid, vertical and bottom mounted, square objects (2 cm of side length in a 4 cm square array) that perforate the free surface of the liquid (see the experimental set-up in Fig. 1). The planar coordinates are denoted by x = (x 1 , x 2 ) whilst the vertical upward direction has the coordinate x 3 ; the origin is prescribed to be at the mean free surface. Under the usual assumptions of linear wa-arXiv:1912.03232v2 [physics.flu-dyn] 16 May 2020 FIG. 1: Experimental setup: A crystal is assembled using square shaped aluminum tubes 7 cm in height arranged in a square array with different orientations using a plastic positioning frame at the bottom of the tank (80 × 80 cm 2 with 60-degree oblique edges made of soft polystyrene to mimic PMLs). A mechanical straight paddle holding a small plastic cylinder is used to generate water waves. The tank is continuously illuminated and images of water waves are recorded with a high speed camera placed on the top. A black and white random pattern is placed under the tank to provide water elevation measurement using an image cross-correlation algorithm. The experimental setup was inspired by the work of Moisy et al; [START_REF] Moisy | A synthetic Schlieren method for the measurement of the topography of a liquid interface[END_REF] .

ter wave theory [START_REF] Lighthill | Waves in Fluids[END_REF] , where the fluid is assumed to be inviscid, incompressible with irrotational flow, there exists a velocity potential Φ 25 such that

Φ(x,t) = ℜe φ (x) cosh(k(x 3 + h)) cosh(kh) exp(-iωt) , (1)
where ω denotes the angular frequency. The wavenumber, k the real positive solution of the dispersion relation

gk + σ ρ k 3 tanh(kh) = ω 2 , (2)
is used as a proxy for the frequency; 2 in Eq. ( 2) g = 9.81 m s -2 is gravitational acceleration, σ = 0.07 N m -1 is the surface tension between air and water and ρ = 10 3 kg m 3 the water density. Then φ , the reduced potential, satisfies the Helmholtz equation,

∇ 2 x + k 2 φ (x) = 0, (3)
where this equation holds at the mean free surface and the subscript x indicates differentiation with respect to x. and no-flow boundary conditions on the vertical rigid cylinders: taking n = (n 1 , n 2 ) as the unit outward normal to the square tubes' surface, ∂ φ /∂ n = 0 on each of them.

When the problem is posed in terms of the reduced potential, φ , as the Helmholtz equation, with periodically arranged inclusions (the tubes), this directly maps across to the phononic crystal literature. Recognising the periodicity guides us to setting φ (x 1 , x 2 ) = φ κ (x 1 , x 2 ) exp(iκ • x) with κ as the Bloch wavenumber and φ κ as the periodic piece of the Bloch solution 2 . A key ingredient, that guides the experiments, is an understanding of the dispersion relation, relating k to the Bloch-wavenumber, κ spanning the BZ κ ∈ [0; π/L] × [0; π/L] (see Fig 2(a)) for an infinite perfectly periodic square lattice system; we determine this relationship numerically. The geometry and band structures are shown in Fig. 2; for a square array, lattice constant L, the irreducible Brillouin zone (IBZ) is an eighth of the BZ. Despite this, we opt to plot around a quadrant of the BZ as this will incorporate the two distinct Dirac cones that are essential for our valley-Hall states. The desired quadrant has the following vertices:

X = (π/L, π/L), N = (π/L, 0), Γ = (0, 0), M = (0, π/L).
Results -The unrotated cellular structure chosen, Fig. 2(a), contains, both, horizontal and vertical mirror symmetries along with four-fold rotational symmetry. Hence, in its entirety the structure has C 4v point group symmetries. Notably, it is the presence of these mirror symmetries that yield Dirac cones along the outer edges of the BZ, Fig. 2(b,c) [START_REF] He | The Emergence of Dirac points in Photonic Crystals with Mirror Symmetry[END_REF][START_REF] Makwana | Topological beam-splitting in photonic crystals[END_REF][START_REF] Makwana | Tunable three-way topological energysplitter[END_REF][START_REF] Xia | Observation of valleylike edge states of sound at a momentum away from the high-symmetry points[END_REF][START_REF] Manzanares-Martinez | Dirac cone switching in two-dimensional photonic crystals by rotation of the crystal basis[END_REF] . Note, that rectangular structures (wallpaper group P2mm) which possess these mirror symmetries will also yield these nonsymmetry repelled Dirac cones [START_REF] Makwana | Topological beam-splitting in photonic crystals[END_REF][START_REF] Makwana | Tunable three-way topological energysplitter[END_REF] . In contrast to hexagonal structures the position of these degeneracies can be tuned by varying the geometrical or material parameters of the system [START_REF] Xia | Observation of valleylike edge states of sound at a momentum away from the high-symmetry points[END_REF] . By rotating the internal square inclusion, both mirror symmetries are broken thereby yielding the band-gap shown for the dispersion curves in Fig. 2(d). The residual valleys, that demarcate the band-gap, are locally imbued with a nonzero valley-Chern number [START_REF] Ochiai | Photonic realization of the (2+1)-dimensional parity anomaly[END_REF] ,

C v = i 2π S ∇ κ ×φ * κ (x)•∇ x φ κ (x)dκ = i 2π γ φ * κ (x)∇ x φ κ (x)•dl (4)
where the path integrated around (γ) encircles a particular valley and the superscript * denotes the complex conjugate. Despite the calculation (and name) of the C v resembling that of its TRS breaking counterpart, namely the Chern number, there is an important difference: the former is not a quantized quantity whilst the latter is. The surface associated with γ is not on a closed manifold hence the Gauss-Bonnet theorem [START_REF] Berger | Differential Geometry: Manifolds, Curves, and Surfaces[END_REF] does not hold. Despite this, the opposite pseudospin modes have a bijective relationship with sgn(C v ) = ±1 which itself can be classed as a topological integer [START_REF] Qian | Theory and Experimental Investigation of the Quantum Valley Hall Effect[END_REF] . From this we can apply the bulk-boundary correspondence for certain edge terminations thereby guaranteeing the existence of valley-Hall edge states.

Motivated by this, we place a perturbed cellular structure, that contains a positively or negatively rotated inclusion, above its reflectional twin. This results in a pair of gapless edge modes that almost span the entirety of the band-gap, see Fig. 3. Here we use "gapless" to refer to the crossing of the concave and convex (opposite parity) modes. This distinguishes valley-Hall systems, that are topological, with those that are not and have coupled edge states; for example, the armchair termination within hexagonal structures produces gapped edge states that are, in turn, less robust [START_REF] Fefferman | Bifurcations of edge states-topologically protected and non-protected-in continuous 2d honeycomb structures[END_REF] .

The gapless nature of the states, and in turn the applicability of the Gauss-Bonnet theorem, is contingent upon the termination chosen containing projections of valleys with iden- tical sgn(C v ). Unique to this specific square structure, the different parity eigenmodes belong to the same interface (see Fig. 3), rather than different interfaces. Despite this both, concave and convex states, have opposite parity and hence remain orthogonal. The relationship between the interfaces arises due to the mirror-symmetry relationship between the media either side of the interface in Fig. 3. This also implies that a right-propagating mode along one of the interfaces is a left-propagating mode on the other. A numerical illustration of this phenomenon is found in [START_REF] Makwana | Tunable three-way topological energysplitter[END_REF] where decaying Hermite polynomials were used to oust a specific parity edge state, along both interfaces, non-simultaneously. This phenomenon does not occur for hexagonal structures where the different parity eigenmodes belong to different interfaces. This relationship between the two interfaces allows for propagation, within our square structure, that is ordinarily forbidden within graphene-like structures. Coupling between modes, that are hosted along different interfaces, is crucial for energy navigation around sharp corners [START_REF] Makwana | Geometrically navigating topological plate modes around gentle and sharp bends[END_REF] and within complex topological domains [START_REF] Makwana | Designing multidirectional energy splitters and topological valley supernetworks[END_REF][START_REF] Cheng | Robust reconfigurable electromagnetic pathways within a photonic topological insulator[END_REF][START_REF] Wu | Direct observation of valley-polarized topological edge states in designer surface plasmon crystals[END_REF][START_REF] Qiao | Electronic highways in bilayer graphene[END_REF] . Further explanation for this phenomenon can be found in [START_REF] Makwana | Topological beam-splitting in photonic crystals[END_REF][START_REF] Makwana | Tunable three-way topological energysplitter[END_REF] . The propagation of water waves is imaged at the surface of the water tank of Fig. 1. A mechanical paddle holding a circular cylinder is shaken at a controlable frequency. Cylindrical waves originating from the monopolar source are observed numerically and experimentally in Fig. 4(a,e). The experimental setups for a topologically nontrivial interface, with two different lengths, is shown in Figs. 4(b,f); the upper/lower halves have square inclusions rotated clockwise/ anti-clockwise in order to break the mirror symmetries and generate the valley-edge states required. Images were acquired by a high speed camera and post processed using a cross-correlation algorithm 1 ; each image was discretized into 360 areas each composed of 16 pixels.

Full-wave numerical simulations, performed using COM-SOL Multiphysics (a commercial finite element scheme), for tightly confined valley-Hall edge states, Figs. 4(c,d), show excellent agreement with the experiments, Figs. 4(g,h), despite our model not taking into account contact-line effects that occur between the water and the solid pillars, viscosity or nonlinearity. These square structure valley-Hall edge states have longer-wavelengths than their hexagonal counterparts and hence the distance between the pillars is subwavelength. A frequency modulated monopolar source is generated that ignites the even-parity valley-Hall edge state. The observed patterns are associated with the surface curvature where the coloured regions are indicative of the vertical elevation of the water level. The localisation of the topological In the absence of any rods the energy would radiate isotropically away from the source. The broadbandedness of this effect is demonstrated via the experimental results shown in the supplementary material. The tight-confinement of these dispersive water waves, within a strategically designed square structure, is a highly nontrivial and unique observation.

We now strategically extend our earlier design, Figs. 4(b, f), to engineer four structured quadrants that results in a threeway topological energy-splitter, Fig. 5. We rotate the bottomright and top-right inclusion sets anti-clockwise and clockwise, respectively, thereby creating four distinct domain walls upon which the valley-Hall states reside. The monopolar source triggers a wave, from the leftmost interface, into upward and downward modes along with continued rightward propagation. Incidentally, the most pronounced displacement pattern is along the two geometrically distinct horizontal interfaces. This continued rightward propagation is forbidden for hexagonal systems [START_REF] Makwana | Designing multidirectional energy splitters and topological valley supernetworks[END_REF][START_REF] Cheng | Robust reconfigurable electromagnetic pathways within a photonic topological insulator[END_REF][START_REF] Wu | Direct observation of valley-polarized topological edge states in designer surface plasmon crystals[END_REF][START_REF] Qiao | Electronic highways in bilayer graphene[END_REF] . For coupling between the incident mode and the right-sided mode the chirality's must match and this does not happen for hexagonal structures. Contrastingly, this mismatch in chirality is overcome for the square structure as the right-sided interface is the reflectional partner of the left-sided interface. Hence, the incident mode need only to couple to itself in order to continue it's rightward propagation. This subtle relationship between the mirror-symmetry generated Dirac cones, and the subsequent, mirror-symmetry related interfaces allows for propagative behaviour not readily found within the valleytronics literature, Fig. 5.

In the experiments there are multiple loss mechanisms on the lengthscales we are operating at: viscous attenuation [START_REF] Padrino | Correction of Lamb's dissipation calculation for the effects of viscosity on capillary-gravity waves[END_REF] , contact line losses associated with the frictional drag of the meniscus moving up and down the rigid pillars i.e. contact angle hysteresis [START_REF] Joanny | A model for contact angle hysteresis[END_REF] , Marangoni effects due to surface tension variations [START_REF] Craster | Dynamics and stability of thin liquid films[END_REF] and their effect on capillary-gravity waves, and then the nonlinear inertial effects that are ignored through linearisation of the Navier-Stokes equations [START_REF] Lighthill | Waves in Fluids[END_REF] . Inspecting Fig. 5, which is a simulation, we note that the amplitude of the interface modes propagating up and down are about 1/10 of that propagating left to right. In Fig. 5(b) we insert losses, by lumping them into a complex wave velocity, of just 2% (we choose this to be simply illustrative and to demonstrate that losses can in an experiment obscure the subtle effects often sought in topological systems) and this reduces the signal quite dramatically in the up and down interfaces; in experiments the amplitudes were too small to be accurately measured and we attribute this to the loss mechanisms we describe above.

We have experimentally shown the existence of topological valley-Hall transport for gravity-capillary water waves within a non-hexagonal structure. We have also simulated a threeway topological multiplexer for the same highly-dispersive system and cautioned that losses may lead to low amplitudes. These demonstrations open up a useful way for design in energy transport: the conventional symmetry constraints asso-ciated with hexagonal structures can be relaxed leading to richer designs of waveguides and multiplexers within highlydispersive systems. We note that the measured vertical displacement of the liquid surface η in Figure S1 is related to the potential φ computed numerically via

η(x 1 , x 2 ,t) = ℜ{-i ω g φ (x 1 , x 2 ) exp(-iωt)} , ( S8
)
where t is the time variable, g and ω denote gravitational acceleration and angular frequency, respectively, as defined in the main article.

BROADBAND DOMAIN WALL STATES

The monopolar source was used to ignite the edge states in Figs.

FIG. 2 :

 2 FIG. 2: Geometry, band structure and topological features: (a) Periodic cell (physical space) of the square lattice with sidelength L showing a square inclusion of sidelength l s inside it. Mirror-symmetry breaking rotation (arrows) and lines (dashed) are also shown. (b) In reciprocal space, the points ΓNX denote the extrema of the IBZ that we extend to ΓNXM to show the two topologically inequivalent regions; the two distinct sgn(C v ) values are indicated by ± signs around the perimeter of the BZ and these are associated with the + perturbation in panel (a) (theperturbation would result in opposite sgn(C v )'s, see 29 ). The sgn(C v ) positions resemble those in 27,28 . (c) Band diagram for the configuration in (a), with two circles marking the position of the strategically engineered Dirac cones(d) Band diagram when the inclusion is rotated through an angle of 20 • . A band-gap highlighted in green emerges from the symmetry breaking perturbation at Dirac points.

FIG. 3 :

 3 FIG. 3: Valley-Hall edge states: Band diagram for a ribbon with the upper/lower inclusions rotated clockwise/anti-clockwise. The real parts of the even and odd eigenmodes within the band-gap are shown (in red) as are several close-by ribbon modes (also in red); the latter's values range between their max and min and a.u stands for arbitrary units (η is the vertical displacement of the surface).The blue curves are from Fig.2(d), i.e. bulk modes along ΓM. Numerically, using finite elements, we take a long ribbon of N inclusions, apply Dirichlet boundary conditions top and bottom and extract the modes decaying away from the interface.

FIG. 4 :

 4 FIG. 4: Valley-Hall edge states: Experiments, and simulations. (a,d) Experimental set-up showing the top view of the water surface perforated by 4 × 8 and 8 × 8 square rigid inclusions respectively. (b,e) Corresponding numerical calculations; (c,f) Experimentally observed wavepatterns. These valley-Hall states are generated by a monopolar source operating at a frequency of 7.3 Hz and placed 6 cm from the domain wall edge and propagate between inclusions rotated by ±20 • .

FIG. 5 :

 5 FIG. 5: Numerical illustration for a three-way splitter at 7.3 Hz: Four quadrants of alternately squares rotated clockwise and counter-clockwise. The valley-Hall phases, associated with the propagating modes, are shown within the yellow inset. (a) has no losses whilst (b) incorporates attenuation introduced by considering a complex wave velocity with 2% imaginary part relative to the real part.

FIG. S1 :

 S1 FIG. S1: (a) Numerical and (b) experimental propagation of cylindrical water waves from a monopolar source.

4 ,

 4 FIG. S2: Valley-Hall edge states: Experiments, and simulations. The water surface is perforated by 4 × 8 square rigid inclusions. Numerics (top) and experiments (bottom) are shown at frequencies (a, f) 6.7 Hz, (b, g) 7 Hz, (c, h) 7.3 Hz, (d, i) 7.6 Hz and (e, j) 7.9 Hz.

Supplementary Material-The supplementary material contains additional experimental images for valley-Hall edge states in crystals of different size. Also shown are images and simulations in the absence of the crystal. It also contains basic information about the finite element models used.

We acknowledge V. Pêcheur, E. Carry, T. Daugey and F. Chollet for the experimental part. This work has been supported by the EIPHI Graduate School (contract ANR-17-EURE-0002) and by the French Investissements d'Avenir program, project ISITE-BFC (contract ANR-15-IDEX-03). R.V.C. and M.P.M. thank the UK EPSRC for their support through grants EP/L024926/1 and EP/T002654/1.

Supplementary Material

SIMULATING WATER WAVES USING FINITE ELEMENTS

The implementation of the water wave equation within a finite element package (such as Comsol Multiphysics) is fairly straightforward. The way to proceed is as follows: We first multiply the main equation valid at the mean free surface

(with the subscript x indicating differentiation with respect to x) by a smooth test function Ψ with a compact support in the periodic cell Y = [0; L] × [0; L] before integrating over Y . Using Green's formula, and the fact that the Neumann boundary condition on the inclusion in Y is a natural condition we obtain the so-called weak form of the Helmholtz equation

Solving the Floquet-Bloch spectral problem amounts to looking for the countable set of real positive eigenvalues k 2 and associated non-zero eigensolutions φ κ (x) for a given κ such that φ (

Here, κ is the Bloch wavenumber and φ κ (x) is the periodic piece of the Floquet-Bloch solution. When we let κ vary within the Brillouin zone BZ = [0; π/L] × [0; π/L], we obtain the band spectrum. It is customary in condensed matter theory to only consider the edges of the BZ, and this enables us in our water wave problem to represent the band spectrum with dispersion curves, rather than dispersion surfaces, in a way similar to what was done in S1,S2 .

The discrete formulation is set up with nodal elements (first-order triangular elements where the unknowns are the values of the scalar field at the vertices of the triangles and the interpolated field is piecewise linear in the triangles). In order to find Floquet-Bloch modes, with the finite element method, we must discretise our grid and apply discrete Bloch conditions. A scalar discrete field on the square cell Y with Floquet-Bloch conditions relates the left and the right sides, and the bottom and the top sides of Y , respectively. The details of this algorithm can be found e.g. in S3 . We implement the Floquet-Bloch water wave problem in comsol multiphysics using the global PDE equation form and built-in Floquet-Bloch conditions. We model the scattering problem using the weak form of the Helmholtz equation (S1) where the right hand side now has a forcing term (a point source). Perfectly matched layers (PMLs) provide a reflectionless interface between the region of interest (the waterwave tank of square shape when viewed from above) and the PMLs (four rectangles and four small squares surrounding the square domain) at all incident angles. PMLs were originally introduced by Berenger S4 in electromagnetism, where the regions of PMLs consisted of anisotropic absorptive dielectric media. In our case, we consider some 'liquid PMLs' defined by the Helmholtz equation, ) which holds at the mean free surface of the liquid. Moreover, the matrix µ is given by µ = s y s z s x 0 0

and the scalar parameter

where s x ,s y and s z represent complex stretched coordinates (as introduced by Berenger S4 ) defined as,

for the absorption of water waves along the x 1 -direction, with the roles of x and y interchanged for absorption of water waves along the x 2 direction. Moreover, for absorption of water waves along the x 1 = x 2 direction one considers s x = s y = a + ib , s z = 1 . (S7)