# A complex-scaled boundary integral equation for time-harmonic water waves

Anne-Sophie Bonnet-Ben Dhia, Luiz Faria, Carlos Pérez-Arancibia

► To cite this version:

Anne-Sophie Bonnet-Ben Dhia, Luiz Faria, Carlos Pérez-Arancibia. A complex-scaled boundary integral equation for time-harmonic water waves. SIAM Journal on Applied Mathematics, 2024, 84 (4), pp.1532-1556. 10.48550/arXiv.2310.04127. [hal-04286461]

HAL Id: hal-04286461

https://hal.science/hal-04286461v1

Submitted on 15 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
A COMPLEX-SCALED BOUNDARY INTEGRAL EQUATION FOR
TIME-HARMONIC WATER WAVES

ANNE-SOPHIE BONNET-BEN DHIA*, LUIZ M. FARIA*, AND
CARLOS PÉREZ-ARANCIBIA†

**Abstract.** This paper presents a novel boundary integral equation (BIE) formulation for the two-dimensional time-harmonic water-waves problem. It utilizes a complex-scaled Laplace's free-space Green's function, resulting in a BIE posed on the infinite boundaries of the domain. The perfectly matched layer (PML) coordinate stretching that is used to render propagating waves exponentially decaying, allows for the effective truncation and discretization of the BIE unbounded domain. We show through a variety of numerical examples that, despite the logarithmic growth of the complex-scaled Laplace's free-space Green's function, the truncation errors are exponentially small with respect to the truncation length. Our formulation uses only simple function evaluations (e.g. complex logarithms and square roots), hence avoiding the need to compute the involved water-wave Green's function. Finally, we show that the proposed approach can also be used to find complex resonances through a *linear* eigenvalue problem since the Green's function is frequency-independent.

**Key words.** Water waves, finite-depth, perfectly matched layers, boundary integral equations, Nyström method.

**1. Introduction.** Wave propagation and scattering problems in physics and engineering often involve unbounded domains that include infinite boundaries. These types of models are particularly important in various scenarios, such as seismic waves interacting with the free surface, outdoor sound propagation, guided (acoustic, electromagnetic, elastic) wave phenomena, and water waves propagating at the ocean's free surface, to mention a few. This paper specifically focuses on the linear (finite-depth) time-harmonic water wave problem, which features two infinite boundaries modeling the bottom topography and the free surface. This classical problem can be traced back to the work of Cauchy [10] and originated from the study of the hydrodynamics of floating and submerged bodies. Extending beyond water waves, the simple model considered in this paper could serve as a prototypical example for investigating a wide range of simplified wave problems where the oscillations of volumetric waves can be neglected at the scale of surface waves. This is typically the case in plasmonics, where the relevant physical phenomenon is the presence of (subwavelength) surface plasmon polaritons [27].

While boundary integral equation (BIE) methods are well known for their ef-
fectiveness in handling unbounded domains, they encounter notable challenges when
applied to problems involving infinite boundaries. BIEs naively formulated using free-
space Green's functions over the entire domain inherently result in equations defined
on the domain's infinite boundaries. Therefore, suitable truncation techniques become
essential to enable their numerical solution, provided, of course, sufficient decay of the
integrands exists. The simple approach of abruptly truncating the computational do-
main at some distance far from the wave sources, offers limited success in the case of
the Helmholtz equation due to the slow decay of the free-space Green's function, and
becomes completely ineffective when applied to two-dimensional water wave problems
due to the logarithmic growth of Laplace's free-space Green's function. (We mention
in passing, that a domain decomposition approach that uses the free-space Green's
function together with a BIE formulation over a bounded domain, has been pursued

<footnote>*POEMS, CNRS, Inria, ENSTA Paris, Institut Polytechnique de Paris, France.</footnote>

<footnote>†Department of Applied Mathematics, University of Twente, Enschede, The Netherlands</footnote>
in [31, 44]. In those works, a Dirichlet-to-Neumann map based on a modal series expansion of the scattered velocity potential is used to effectively reduce the problem to a bounded domain, an idea that later seemed to have influenced the development of finite element methods for water wave problems [21].)

Leveraging the geometric simplicity of the infinite boundaries involved, many BIE formulations for water wave problems rely instead on the use of problem-specific Green's functions that exactly satisfy the relevant boundary conditions over infinite flat portions of the domain's boundaries [2, 3, 6, 11, 15, 26]. While these formulations yield BIEs posed on a bounded portion of the domain's boundary, the numerical evaluation of the associated Green's function often becomes the major bottleneck in the solution process since simple expressions for the free surface Green's function are in general not available. Intricate numerical procedures have been developed in order to compute the free surface Green's function (see e.g. [43] and [26] for a comparison of different techniques and recent developments on this subject). This approach typically entails Fourier-like integrals, non-standard special functions, and/or power series expansions, which in view of the large number of Green's function evaluations involved in the BIE solution process, can be computationally demanding [32]. Furthermore, off-the-shelf acceleration techniques for the Laplace equation, such as the fast multipole method (FMM) [16], cannot be directly applied to the free surface Green's function. Interestingly, an explicit closed-form expression (in terms of exponential integrals) for the infinite depth problem in two dimensions, is presented in [19]. Even in this relatively simple case, developing a multipole expansion is non-trivial [34]. Unfortunately, as far as the authors are aware, there are no closed-form expressions (nor multipole expansions) available for the two-dimensional finite depth problem as well as for the three-dimensional finite and infinite depth problems. We refer the reader to [28, 18, 22] where these topics are discussed in great detail. In addition to the aforementioned efficiency issues, it is important to mention that problem-specific Green's functions always incorporate image sources which are essential for enforcing the free surface and Neumann boundary conditions at flat boundaries. Their presence, however, imposes significant constraints on the types of domains to which the associated BIE formulation can be applied. In fact, domains suitable for the application of such BIE formulations must be strictly confined within the strip defined by the flat free surface and the flat bottom topography.

Over the last decade, there have been notable advancements in the development of novel BIE formulations for Helmholtz wave problems involving unbounded interfaces [8, 25]. These formulations revisit the idea of utilizing the simple efficient-to-evaluate free-space Green's function and truncating the resulting BIE posed on infinite boundaries, but they incorporate mechanisms to effectively account for the correct "outgoing" solution behavior outside the bounded region of interest. In detail, [8] put forth the so-called Windowed Green's function (WGF) method, which features a real-valued suitably-scaled smooth window function that multiplies the free-space Green's function kernel to effectively reduce the BIE to a bounded domain. In essence, this method exploits the fact that windowed integrands exhibiting decaying and oscillatory behavior can be neglected from the tails of the integrals over the infinite boundaries, introducing errors that decay super-algebraically fast as the size of the truncated region increases. Around the same time, a technique was introduced by Lu et al. in [25] that utilizes a BIE formulation of the PDE boundary value problem transformed by a perfectly matched layer (PML) complex stretching and subsequently truncated to a bounded domain. As the WGF method, the resulting PML-BIE formulation is posed on a bounded portion of the infinite boundary and features a modified free-space
Green's function that encodes the outgoing oscillatory character of the PDE solution.

Although the WGF proves effective at solving Helmholtz problems in the context of layered and periodic media as well as waveguides [8, 9, 7, 39], its direct application to the water wave problem under consideration is not feasible. While it may be tempting to attribute the failure of the WGF method to the logarithmic growth of the free space Green's function, an argument can be made that by combining an easy-to-evaluate harmonic function with the free-space Green's function, a decaying Green's function can be formed (in a manner akin to the construction used in [12]). However, even in that scenario, it is easy to see that the WGF method would fail to approximate the correct radiative water wave problem solution. Indeed, when considering real boundary data, the WGF method produces real-valued (stationary-wave) solutions that fail to approximate the correct radiative water-wave problem solutions. We suspect that the failure of the WGF method in the context of the water waves problem is due to the fact that, unlike in the Helmholtz case, the Laplace free space Green's function contains no information about the "correct" radiation condition of the problem.

Building upon the PML-BIE method for the Helmholtz equation [25], this paper proposes a complex-scaled BIE method for the time-harmonic water wave problem in two dimensions. Our BIE formulation is derived from a direct integral representation of the analytically extended scattered velocity potential, obtained by suitably stretching the unbounded domain in the complex plane. This representation involves integrals over the infinite boundaries corresponding to the free surface and the bottom, and is obtained by applying Green's third identity to the velocity potential and the logarithmic free-space Green's function over the complex stretched domain. Interestingly, despite the unboundedness of Laplace's Green's function along the complex stretched infinite boundaries, the resulting integral representation formula remains meaningful due to the provable exponential decay of the complexified velocity potential. Indeed, this decay is achieved by devising complex stretching coordinates on which the modal decomposition of the far-field scattered velocity potential, which encompasses a (constant-amplitude) propagating surface-wave mode and infinitely many evanescent modes, decays exponentially fast. Upon enforcing the boundary conditions and truncating the infinite boundaries, we then obtain a BIE that we can solve numerically via a high-order Nyström method which suitably handles the complexified kernels. Through our numerical experiments, we demonstrate that the error in the computed solutions using this methodology decays exponentially fast as the length of the PML region increases.

A noteworthy difference with the PML-BIE technique applied to the Helmholtz equation [25] stems from the fact that, due to the lack of wave-like behavior in the fundamental solution of Laplace's equation, its complex extension still exhibits a logarithmic growth at infinity (a more detailed discussion is given in Section 3). This makes it not immediately obvious that truncation errors stemming from the method will converge to zero as the length of the PML region is increased. While for the exact solution it is straightforward to show that the truncation of the involved integral operators to a bounded domain yields exponentially small errors (for the analytic extension of the exact solution decays exponentially at infinity), proving a stability result on the proposed truncation technique remains an open question. The main difficulty lies in the fact that the kernels in the proposed integral equation do not decay exponentially at infinity (in fact they may even grow logarithmically in some cases), and thus the integral operators over the unbounded curves are only well-defined for densities decaying sufficiently fast. The numerical results of Section 5
suggest, however, that the lack of decay in the Green's function is not a fundamental problem in the method, and much like in the Helmholtz case we observe truncation errors that decay exponentially fast as the length of the PML layer is increased.

To further showcase the capabilities and possible advantages of the proposed methodology, we numerically investigate the challenging problem of the scattering of a surface wave by an infinite step bottom topography, where a direct application of a BIE formulation that makes use of the constant depth water-wave Green's function (satisfying free surface and Neumann boundary conditions at the flat infinite boundaries), is not feasible [4, 5, 35]. Finally, we apply our methodology to approximate the resonant frequencies that arise in floating body problems [17, 29]. In contrast to alternative approaches for computing complex resonances based on the use of the water-waves Green's function, which yields a (non-linear) eigenvalue problem where the frequency appears in the kernels of the integral operators [17], our PML-transformed formulation gives rise to a linear (generalized) eigenvalue problem that can be efficiently solved using standard matrix eigenvalue solvers.

The paper is organized as follows. In Section 2 we provide the setup of the problem under consideration, recalling some standard results of water-waves theory. Section 3 focuses on the construction of the complex stretching based on the modal decomposition of the scattered velocity potential, and the associated complex-scaled BIE formulation. In Section 4, we introduce the Nyström method employed in our approach. Next, Section 5 presents a variety of numerical experiments designed to validate the proposed methodology and also demonstrate its practical applicability (the method’s implementation, as well as the scripts used to generate the figures, are publicly available in an accompanying GitHub repository[^1]; all of the numerical results of Section 5 can be reproduced by following the instructions there). In Section 6, we briefly discuss some open problems and theoretical difficulties. Finally, we present in Section 7 some concluding remarks and possible future directions.

**2. Problem setup.** This work focuses on linear and time-harmonic water waves in a bath of finite depth $d$. For the sake of presentation simplicity, in what follows we derive dimensionless equations, using the depth $d$ as a characteristic length scale (this means, in particular, that the dimensionless depth is one).

Throughout this paper, we let $\Omega_f \subset \mathbb{R}^2$ denote the infinite fluid domain bounded above by the free surface $\Gamma_s$ and below by a bottom topography $\Gamma_b$. We also let $\Omega_o \subset \mathbb{R}^2$ denote the domain occupied by bounded solid obstacles (either fully or partially) immersed in the fluid, with the immersed part of their boundary denoted by $\Gamma_o := \partial\Omega_o \cap \bar{\Omega}_f$. More precisely, the unperturbed free surface is assumed to be placed at $\mathbb{R} \times \{0\}$ so that $\Gamma_f := \mathbb{R} \times \{0\} \setminus \bar{\Omega}_o$. The bottom topography, on the other hand, is assumed to be given by a constant depth $1$ outside of a bounded domain, i.e., there exists an $L_0 > 0$ such that $\{|x_1| > L_0\} \times \{-1\} \subset \Gamma_b$. In the derivations that follow we use $L_1 > \max\{L_0, \text{diam } \Omega_o\}$. A schematic of the geometry is shown in Figure 1.

<footnote n="1">¹https://github.com/maltezfaria/HarmonicWaterWaves/tree/v0.1</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1: Schematics of the domain and main variables involved in the setup of water wave problem (2.1) under consideration.</figcaption></figure>

Here we seek to solve the following boundary value problem:

$$ (2.1a) \qquad \Delta\varphi(\mathbf{x}) = 0, \qquad \mathbf{x} := (x_1, x_2) \in \Omega := \Omega_f \setminus \bar{\Omega}_o, $$

$$ (2.1b) \qquad \nabla\varphi(\mathbf{x}) \cdot \mathbf{n}(\mathbf{x}) + \nu\varphi(\mathbf{x}) = f_1(\mathbf{x}), \qquad \mathbf{x} \in \Gamma_f := \mathbb{R} \times \{0\} \setminus \bar{\Omega}_o, $$

$$ (2.1c) \qquad \nabla\varphi(\mathbf{x}) \cdot \mathbf{n}(\mathbf{x}) = f_2(\mathbf{x}), \qquad \mathbf{x} \in \Gamma_b \cup \Gamma_o, $$

$$ (2.1d) \qquad \lim_{x_1 \to \pm\infty} \partial_{x_1} \varphi(\mathbf{x}) \mp ik\varphi(\mathbf{x}) = 0, $$

where letting $\omega > 0$ denote the angular frequency and $g$ the gravitational acceleration, the term $\nu = (\omega^2 d)/g$ in the Robin boundary condition (2.1b) is a dimensionless parameter, while the constant $k > 0$ in the radiation condition (2.1d) is uniquely determined by the dispersion relation

$$ (2.2) \qquad k \tanh(k) = \nu. $$

As usual, $\mathbf{n} \in \mathbb{R}^2$ denotes the unit normal which points toward the interior of the fluid domain $\Omega$ (see Figure 1). It is worth noting that the unknown velocity potential $\varphi$ is related to the fluid velocity $\mathbf{u}$ through $\mathbf{u}(\mathbf{x}, t) = \text{Re}(e^{-i\omega t}\nabla\varphi)$, where $t$ is the time variable.

We here assume that the boundary data $f_1$ and $f_2$ are (compactly) supported in $\Gamma_f$ and $\Gamma_b$ vanishing over $\{|x_1| > L_2\} \times \{0\}$ and $\{|x_1| > L_2\} \times \{-1\}$, respectively. Note that the step-like bottom topography problem considered in Section 5 is not covered by the previous description; nevertheless, a slight modification of the assumptions made here allows us to extend the proposed PML-BIE approach to this challenging problem.

In order to construct a suitable complex scaling for problem (2.1), we rely on the analytic extension of the PDE solution $\varphi$. To justify the extension and gain insight into the construction of the complex stretching, in what follows in this section we recall some classical results about the modal decomposition of $\varphi$ within semi-infinite strips of constant unit depth [20, 23]. Indeed, letting $L = \max\{L_1, L_2\}$, the exact solution $\varphi$ of (2.1) admits a modal decomposition within the domain $\{|x_1| > L\} \times (-1, 0)$. To define the modal decomposition we first introduce the sequence $\{\gamma_n\}_{n \ge 1}$, such that

$$ (2.3) \qquad \gamma_n \tan(\gamma_n) = -\nu, \qquad n\pi - \frac{\pi}{2} < \gamma_n < n\pi. $$

Then one can prove that, with a suitable choice of the normalization constants ($a_n$),
the functions {$u_n$}$_{n\ge0}$ defined by

$$
(2.4) \qquad u_0(x_2) = a_0 \cosh(k(x_2 + 1)),
$$

$$
(2.5) \qquad u_n(x_2) = a_n \cos(\gamma_n(x_2 + 1)), \quad n \ge 1,
$$

with $k > 0$ being the solution of (2.2) and $\gamma_n$ the solution of (2.3), form an orthonormal basis of $L^2([-1, 0])$. Invoking the radiation condition (2.1d) to select the correct behavior at $\pm\infty$, the following modal expansion for the solution $\varphi$ of (2.1) can be readily obtained:

$$
(2.6) \qquad \varphi(\mathbf{x}) =
\begin{cases}
A_0^+ u_0(x_2) e^{ik(x_1-L)} + \sum_{n \ge 1} A_n^+ u_n(x_2) e^{-\gamma_n(x_1-L)}, & \mathbf{x} \in (L, \infty) \times [-1, 0], \\
A_0^- u_0(x_2) e^{-ik(x_1+L)} + \sum_{n \ge 1} A_n^- u_n(x_2) e^{\gamma_n(x_1+L)}, & \mathbf{x} \in (-\infty, -L) \times [-1, 0],
\end{cases}
$$

where

$$
(2.7) \qquad A_n^{\pm} = \int_{-1}^{0} \varphi(\pm L, x_2) \overline{u_n(x_2)} \, dx_2, \quad n \ge 0.
$$

The well-posedness of (2.1) can be studied by reducing the problem to a bounded
domain via a DtN map approach derived from the modal expansion (2.6). Indeed, it
can be established that the problem (2.1) is well posed for all frequencies except for
a countable set (see e.g. [40]).

Clearly, on each connected component of the exterior domain {|x₁| > L} × [-1, 0]
the velocity potential ϕ is a superposition of the outgoing surface wave and of infinitely
many evanescent modes.

**3. Complex-scaled boundary integral equation.** In this section, we derive a complex scaled integral equation for solving (2.1) with the desirable property that its solutions are exponentially decreasing as |x₁| → ∞. Using the modal decomposition (2.6), we show in Section 3.1 that the velocity potential ϕ admits an analytic extension in x₁ into the complex plane for |x₁| > L (i.e. outside the support of the perturbations). Using the analytic extension of the modal decomposition, we then derive in Section 3.2 a modified water waves problem, as well as its boundary integral representation.

**3.1. Complex stretching.** The modal decomposition (2.6) can be used to study analytic properties of $\varphi$ with respect to the coordinate $x_1$. More precisely, for $x_1 > L$ (resp. $x_1 < -L$) the solution $\varphi$ of (2.1) has an analytic extension to complex values of its first variable with a real part larger than $L$ (resp. smaller than $-L$). To see why this is so, let us consider for instance the side $x_1 > L$. It is clear that for any $N \in \mathbb{N}$, the finite sum

$$
(3.1) \qquad A_0^+ u_0(x_2) e^{ik(z_1-L)} + \sum_{n=1}^{N} A_n^+ u_n(x_2) e^{-\gamma_n(z_1-L)}
$$

is well-defined for complex values of $z_1$, and is a holomorphic function of $z_1$. Moreover, it converges when $N \to +\infty$, uniformly in all compact subsets of $\text{Re}(z_1) > L$, to a function that we denote (with a slight abuse of notation) by $\varphi(z_1, x_2)$. Finally, thanks
to Morera's theorem, a uniform limit of a sequence of holomorphic functions is itself holomorphic, so we conclude that $z_1 \mapsto \varphi(z_1, x_2)$ is analytic within $\text{Re}(z_1) > L$. The exact same argument applies to the analytic extension $z_1 \mapsto \varphi(z_1, x_2)$ for $\text{Re}(z_1) < -L$.

It follows then by the analyticity of $\varphi$ that

$$ (3.2) \qquad (\partial_{z_1}^2 + \partial_{x_2}^2) \varphi(z_1, x_2) = 0, \quad (z_1, x_2) \in \{|\text{Re } z_1| > L\} \times [-1, 0]. $$

We now introduce a (continuous and piece-wise differentiable) stretching path $\tau: \mathbb{R} \to \mathbb{C}$ along which the transformed function

$$ (3.3) \qquad \tilde{\varphi}(x_1, x_2) := \varphi(\tau(x_1), x_2) $$

decays exponentially fast as $|x_1| \to \infty$. In view of the modal decomposition (2.6), the following conditions guarantee the exponential decay of $\tilde{\varphi}(\cdot, x_2)$:

$$ (3.4a) \qquad \operatorname{Im}(\tau(x_1)) > 0 \quad \text{for} \quad x_1 > L $$

$$ (3.4b) \qquad \operatorname{Im}(\tau(x_1)) < 0 \quad \text{for} \quad x_1 < -L $$

Furthermore, in order to recover the original solution to (2.1) for $|x_1| < L$, we require

$$ (3.4c) \qquad \tau(x_1) = x_1 \quad \text{for} \quad |x_1| < L. $$

Finally, we suppose that

$$ (3.4d) \qquad 0 < c < \text{Re}(\tau'(x_1)) < C < \infty \quad \text{for} \quad x_1 \in \mathbb{R}, $$

which implies that $\tau$ is injective.

A simple application of the chain rule then yields the identity $\partial_{z_1}\varphi(\tau(x_1), x_2) = \partial_{x_1}\tilde{\varphi}(x_1, x_2)/\tau'(x_1)$, therefore, in view of (3.2), $\tilde{\varphi}$ defined in (3.3) satisfies

$$ (3.5a) \qquad \nabla \cdot (A(\mathbf{x})\nabla\tilde{\varphi}) = 0 \quad \text{with} \quad A(\mathbf{x}) := \begin{bmatrix} \tau'(x_1)^{-1} & 0 \\ 0 & \tau'(x_1) \end{bmatrix}, \quad \mathbf{x} \in \Omega. $$

(We have chosen a matrix formalism to simplify the discussion in the subsequent section.) Note that by (3.4c), $\tau' = 1$ for $|x_1| < L$, so (3.5a) becomes the Laplace equation in that case. Finally, we note that because $f_1$ and $f_2$, as well as the perturbations of the geometry vanish for $|x_1| > L$, it can be easily shown that $\tilde{\varphi}$ satisfies the same boundary conditions as $\varphi$, i.e.,

$$ (3.5b) \qquad \nabla\tilde{\varphi}(\mathbf{x}) \cdot \mathbf{n}(\mathbf{x}) + \nu\tilde{\varphi}(\mathbf{x}) = f_1(\mathbf{x}), \quad \mathbf{x} \in \Gamma_f, $$

$$ (3.5c) \qquad \nabla\tilde{\varphi}(\mathbf{x}) \cdot \mathbf{n}(\mathbf{x}) = f_2(\mathbf{x}), \quad \mathbf{x} \in \Gamma_b \cup \Gamma_o. $$

It is worth mentioning at this point that the radiation condition (2.1d) is simply replaced by the exponential decay of $\tilde{\varphi}$. Equations (3.5) comprise the PML-transformed water waves system that we seek to solve in lieu of the original system (2.1).

**3.2. Boundary integral equation.** In order to recast the problem as a BIE by leveraging fundamental potential theory results, we first show that our PDE (3.5a) is strongly elliptic in the sense of [30]. Indeed, given the form of the matrix function $A$ introduced in (3.5a) and the stretching path property (3.4d), it follows that

$$ (3.6) \qquad \text{Re}(\xi^* A \xi) = \text{Re} \left( \frac{1}{\tau'} |\xi_1|^2 + \tau' |\xi_2|^2 \right) > \min \left( c, \frac{1}{C} \right) |\xi|^2, \quad \xi \in \mathbb{C}^2, $$
which shows the desired property.

A crucial aspect of the subsequent BIE formulation is the explicit relationship

$$ (3.7) \qquad \tilde{G}(\mathbf{x}, \mathbf{y}) := G(\tau(\mathbf{x}), \tau(\mathbf{y})) = -\frac{1}{4\pi} \log ((\tau(x_1) - \tau(y_1))^2 + (x_2 - y_2)^2), $$

between the Green's function $\tilde{G}$ associated with (3.5a), and the Laplace's fundamental solution $G$, along with the transformation $\tau$ (see Appendix A). Here, the principal logarithm is used with a branch cut on the negative real line and, for the sake of compactness of notation, we have introduced the vector-valued transformation

$$ \tau(\mathbf{x}) := (\tau(x_1), x_2) : \mathbb{R}^2 \to \mathbb{C} \times \mathbb{R}. $$

With these definitions at hand, we are in position to derive the boundary integral representation for $\tilde{\varphi}$, defined in (3.3), over the (unbounded) domain $\Omega$. We then begin by considering the truncated domain $\Omega_M := \{\mathbf{x} \in \Omega : |x_1| < M\}, M > 0$. On $\Omega_M$, the following integral representation formula holds [30]:

$$ (3.8) \qquad \tilde{\varphi}(\mathbf{x}) = \mathcal{D}_{\partial\Omega_M}[\tilde{\varphi}](\mathbf{x}) - \mathcal{S}_{\partial\Omega_M}[\nabla\tilde{\varphi} \cdot A\mathbf{n}](\mathbf{x}), \quad \mathbf{x} \in \Omega_M, $$

where $\mathcal{S}_\Sigma$ and $\mathcal{D}_\Sigma$ denote respectively the single- and double-layer potentials, defined as

$$ (3.9a) \qquad \mathcal{S}_\Sigma[\sigma](\mathbf{x}) := \int_{\Sigma} \tilde{G}(\mathbf{x}, \mathbf{y})\sigma(\mathbf{y}) \, \mathrm{d}\mathbf{s}(\mathbf{y}) \quad \text{and} $$

$$ (3.9b) \qquad \mathcal{D}_\Sigma[\sigma](\mathbf{x}) := \int_{\Sigma} \left( \nabla_y \tilde{G}(\mathbf{x}, \mathbf{y}) \cdot A(\mathbf{y})\mathbf{n}(\mathbf{y}) \right) \sigma(\mathbf{y}) \, \mathrm{d}\mathbf{s}(\mathbf{y}), \quad \mathbf{x} \in \mathbb{R}^2 \setminus \Sigma. $$

In (3.9), $\Sigma$ stands for a generic sufficiently regular (Lipschitz) curve.

Using the fact that $\tilde{\varphi}$ decays exponentially fast as $|x_1| \to \infty$, while the integral kernels in (3.10) grow at most logarithmically, the contribution from the lateral curves $\{\pm M\} \times [-1, 0] \subset \partial\Omega_M$ in (3.8) decays exponentially fast as $M \to \infty$. As such, we may replace $\Omega_M$ by its limit $\Omega$ in (3.8), yielding the following integral representation:

$$ (3.10) \qquad \tilde{\varphi}(\mathbf{x}) = \mathcal{D}_{\Gamma}[\tilde{\varphi}](\mathbf{x}) - \mathcal{S}_{\Gamma}[\nabla\tilde{\varphi} \cdot A\mathbf{n}](\mathbf{x}), \quad \mathbf{x} \in \Omega, $$

where $\Gamma := \partial\Omega = \Gamma_b \cup \Gamma_o \cup \Gamma_f$.

It is worth mentioning that an integral representation on the unbounded curve $\Gamma$, such as (3.10), can only be obtained for $\tilde{\varphi}$, and not for $\varphi$. In fact, in the absence of a complex-stretching equation (3.10) ceases to make sense because $\varphi$ does not decay, and the single-layer kernel exhibits a logarithmic growth, meaning that the integrals are not convergent (even in the conditional sense). This is a fundamental difference from the Helmholtz equation, where the decay of the Green's function and of the scattered field suffices to render integrals over unbounded interfaces conditionally convergent (but not absolutely so).

Next, to obtain an integral equation for $\tilde{\varphi}$ over $\Gamma$, we evaluate the limit $\Omega \ni \mathbf{x} \to \mathbf{x}' \in \Gamma$ in (3.10). Taking into account the jump condition of the the double-layer potential (3.9b) across $\Gamma$ [30], we arrive at

$$ (3.11) \qquad \frac{\tilde{\varphi}(\mathbf{x})}{2} = \mathcal{D}_{\Gamma}[\tilde{\varphi}](\mathbf{x}) - \mathcal{S}_{\Gamma}[\nabla\tilde{\varphi} \cdot A\mathbf{n}](\mathbf{x}), \quad \mathbf{x} \in \Gamma, $$
where $S_Γ$ and $D_Γ$ are respectively the single- and double-layer operators, defined as

$$ (3.12a) \qquad S_Γ[\sigma](\mathbf{x}) := \int_Γ \tilde{G}(\mathbf{x}, \mathbf{y})\sigma(\mathbf{y}) \, ds(\mathbf{y}) \quad \text{and} $$

$$ (3.12b) \qquad D_Γ[\sigma](\mathbf{x}) := \text{p.v.} \int_Γ (\nabla_y \tilde{G}(\mathbf{x}, \mathbf{y}) \cdot A(\mathbf{y})\mathbf{n}(\mathbf{y})) \sigma(\mathbf{y}) \, ds(\mathbf{y}), \quad \mathbf{x} \in Γ. $$

The letters p.v. in front of the integral sign in (3.12b) refer to the fact that the integral must be interpreted in the sense of the Cauchy principal value.

Note that the double-layer kernel is given in explicit terms by

$$ (3.13) \qquad \nabla_y \tilde{G}(\mathbf{x}, \mathbf{y}) \cdot A(\mathbf{y})\mathbf{n}(\mathbf{y}) = \tau'(y_1) \frac{\partial G}{\partial \mathbf{n}_y}(\tau(\mathbf{x}), \tau(\mathbf{y})) = \frac{\tau'(y_1)(\tau(\mathbf{x}) - \tau(\mathbf{y})) \cdot \mathbf{n}(\mathbf{y})}{2\pi(\tau(\mathbf{x}) - \tau(\mathbf{y})) \cdot (\tau(\mathbf{x}) - \tau(\mathbf{y}))}, $$

which is simply a (complex) multiple of Laplace's free-space double-layer kernel evaluated at the complex points $\tau(\mathbf{x})$ and $\tau(\mathbf{y})$. Furthermore, under the additional assumption that $Γ$ and $\tau$ are smooth at $\mathbf{x}$, it holds that

$$ (3.14) \qquad \tau(\mathbf{y}) = \tau(\mathbf{x}) + J(\mathbf{x})(\mathbf{y} - \mathbf{x}) + \mathcal{O}(|\mathbf{x} - \mathbf{y}|^2), $$

$$ (3.15) \qquad \mathbf{n}(\mathbf{y}) = \mathbf{n}(\mathbf{x}) + \mathcal{O}(|\mathbf{x} - \mathbf{y}|) $$

as $\mathbf{y} \to \mathbf{x}$, where

$$ (3.16) \qquad J(\mathbf{x}) = \begin{bmatrix} \tau'(x_1) & 0 \\ 0 & 1 \end{bmatrix} $$

denotes the Jacobian matrix of $\tau$. Noting that $J(\mathbf{x})\mathbf{n}(\mathbf{x}) = \mathbf{n}(\mathbf{x})$ for any $\mathbf{x} \in Γ$, where we have used the fact that $J$ is either the identity (outside the PML) or $\mathbf{n}$ points in the $x_2$ direction (inside the PML), we obtain that

$$ (3.17) \qquad \nabla_y \tilde{G}(\mathbf{x}, \mathbf{y}) \cdot A(\mathbf{y})\mathbf{n}(\mathbf{y}) = \frac{\tau'(x_1)(\mathbf{x} - \mathbf{y}) \cdot \mathbf{n}(\mathbf{x})}{J(\mathbf{x})(\mathbf{x} - \mathbf{y}) \cdot J(\mathbf{x})(\mathbf{x} - \mathbf{y})} + \mathcal{O}(1) \quad \text{as } \mathbf{y} \to \mathbf{x}. $$

Therefore, since $(\mathbf{x} - \mathbf{y}) \cdot \mathbf{n}(\mathbf{x}) = \mathcal{O}(|\mathbf{x} - \mathbf{y}|^2)$ as $Γ \ni \mathbf{y} \to \mathbf{x} \in Γ$ for a smooth curve, we conclude that the double-layer kernel (3.13) is a bounded function at $\mathbf{y} = \mathbf{x}$, and hence the principal value in (3.12b) is not necessary in such cases.

With the validity of the representation formula (3.10) established, and the availability of the free-space fundamental solution (3.7), we can now derive a BIE for $\tilde{\varphi}$. Indeed, using the boundary conditions (3.5b) and (3.5c) to replace $\nabla\tilde{\varphi} \cdot A\mathbf{n} = \tau'\tilde{\varphi} \cdot \mathbf{n}$ in favor of $\tilde{\varphi}$ in (3.11), we arrive at

$$ (3.18) \qquad -\frac{\tilde{\varphi}(\mathbf{x})}{2} + D_Γ[\tilde{\varphi}](\mathbf{x}) + \nu S_{Γ_f} [\tau'\tilde{\varphi}](\mathbf{x}) = S_{Γ_f}[f_1](\mathbf{x}) + S_{Γ_b \cup Γ_o}[f_2](\mathbf{x}), \quad \mathbf{x} \in Γ. $$

To further simplify the notation, we introduce a “global” Robin coefficient $\alpha : Γ \to \mathbb{R}$, as well as a “global” source $f : Γ \to \mathbb{R}$ given by:

$$ (3.19) \qquad \alpha(\mathbf{x}) = \nu \mathbb{1}_{Γ_f}, \quad f(\mathbf{x}) = f_1(\mathbf{x}) \mathbb{1}_{Γ_f} + f_2(\mathbf{x}) \mathbb{1}_{Γ_o \cup Γ_b}, $$

where $\mathbb{1}_Σ$ denotes the indicator function over a domain $Σ$. Equation (3.18) is then rewritten in a more compact form as

$$ (3.20) \qquad -\frac{\tilde{\varphi}(\mathbf{x})}{2} + D_Γ[\tilde{\varphi}](\mathbf{x}) + S_Γ[\alpha\tau'\tilde{\varphi}](\mathbf{x}) = S_Γ[f](\mathbf{x}), \quad \mathbf{x} \in Γ. $$

In the subsequent sections, we will present the development of a Nyström method for the numerical approximation of solutions to (3.20).
**4. Numerical method.** In this section, we provide the details of the discretization scheme employed, as well as the concrete choice of the (linear) PML function used along with the parameters involved. We also present an argument for the consistency of the truncated BIE which includes an integral equation for the error, that exhibits a right-hand-side that vanishes exponentially as $M \to \infty$. The numerical examples of Section 5 demonstrate that, as expected, the truncation errors decay exponentially with respect to the length of the PML layer.

In order to discretize (3.20), we proceed in two steps. First, after a choice of PML function $\tau$ in (3.4) has been made, the infinite boundary $\Gamma$ is truncated to a finite curve $\Gamma_M = \{\mathbf{x} \in \Gamma : |x_1| < M\}$. Second, the resulting integral equation over $\Gamma_M$ is discretized using a Nyström scheme, wherein a composite P-point Gauss-Legendre quadrature rule is used to approximate the integrals over $\Gamma_M$. A correction, which accounts for the singular and nearly-singular interactions arising in the numerical discretization of the underlying integral operators is obtained by means of a Gauss-Kronrod adaptive quadrature; more details are given in Section 4.2.

The particular choice of the BIE discretization here developed overcomes certain limitations of some existing Nyström methods. For instance, the well-known spectrally accurate Martensen-Kussmaul approach [13, section 3.6] becomes numerically unstable in the current setting due to the underlying kernel splitting involving complex arguments [24]. Other approaches, such as the Density Interpolation Method [33, 14], developed by some of the authors and collaborators, does not handle open curves such as the ones that arise in the current BIE formulation. The hybrid Gauss-trapezoidal quadrature rule of Alpert [1], on the other hand, can in principle be used in the current setting and was the choice of spatial discretization used by Lu et al. in [25].

Throughout this section, we will use $\tilde{\varphi}_{M,h}$ to denote the approximate solution over $\Gamma_M$ using a P-point quadrature rule on elements of approximate size $h$, where the dependency on $P$ is omitted for notational simplicity.

**4.1. Domain truncation.** To ensure a clear and concise presentation, we will focus exclusively on what we call “linear PMLs”, wherein $\tau$ is given by

$$ (4.1) \qquad \tau(x_1) = x_1 + i \begin{cases} c(x_1 + a), & x_1 \le -a, \\ 0, & |x_1| < a, \\ c(x_1 - a) & x_1 \ge a, \end{cases} $$

where $a > 0$ is a parameter controlling the start of the complex stretching, while the parameter $c > 0$ controls its slope. We know from (2.6) that $\tilde{\varphi}(x_1, x_2) = \varphi(\tau(\mathbf{x}), x_2)$ admits the following modal decomposition for $x_1 > L$:

$$ (4.2) \qquad \tilde{\varphi}(x_1, x_2) = A_0^+ u_0(x_2) e^{-ck(x_1-a)} e^{ik(x_1-L)} + \sum_{n \ge 1} A_n^+ u_n(x_2) e^{-\gamma_n(x_1-L)} e^{-ic\gamma_n(x_1-a)} $$

with the coefficients $A_n^+$ given in (2.7). This shows that

$$ (4.3) \qquad \tilde{\varphi}(x_1, x_2) \sim \mathcal{O}(e^{-\min(ck, \gamma_1)x_1}) \quad \text{as} \quad x_1 \to \infty. $$

A similar argument holds for $x_1 \to -\infty$, and thus $\tilde{\varphi}$ decays exponentially as $|x_1| \to \infty$.

Given the exponential decay of $\tilde{\varphi}$, we replace the integral equation (3.20) by its truncated version

$$ (4.4) \qquad -\frac{\tilde{\varphi}_M(\mathbf{x})}{2} + D_{\Gamma_M}[\tilde{\varphi}_M](\mathbf{x}) + S_{\Gamma_M}[\alpha \tau' \tilde{\varphi}_M](\mathbf{x}) = S_{\Gamma_M}[f](\mathbf{x}), \quad \mathbf{x} \in \Gamma_M. $$
Subtracting (4.4) from (3.20), it is straightforward to see that the truncation error $\mathcal{E}_M(\mathbf{x}) := \tilde{\varphi}(\mathbf{x}) - \tilde{\varphi}_M(\mathbf{x})$ satisfies the following integral equation:

$$ (4.5) \qquad -\frac{\mathcal{E}_M(\mathbf{x})}{2} + D_{\Gamma_M}[\mathcal{E}_M](\mathbf{x}) + S_{\Gamma_M}[\alpha\tau'\mathcal{E}_M](\mathbf{x}) = -D_{\Gamma_M^c}[\tilde{\varphi}](\mathbf{x}) - S_{\Gamma_M^c}[\alpha\tau'\tilde{\varphi}](\mathbf{x}) $$

for $\mathbf{x} \in \Gamma_M$, where $\Gamma_M^c = \Gamma \setminus \Gamma_M$. Because of the exponential decay of $\tilde{\varphi}$, it follows that the right-hand side is exponentially small in $M$ as $M \to \infty$, which shows the consistency of the truncation technique. To further conclude that $\mathcal{E}_M$ is also exponentially small, however, a stability result of the type

$$ (4.6) \qquad \|\mathcal{L}_M^{-1}\|_{L^2(\Gamma_M) \to L^2(\Gamma_M)} \lesssim M^q, \quad \mathcal{L}[\sigma] := -\frac{\sigma}{2} + D_{\Gamma_M}[\sigma] + S_{\Gamma_M}[\alpha\tau'\sigma] $$

is needed as $M \to \infty$, with $q > 0$ independent of $M$. Although such a result is lacking in our analysis, the numerical tests of Section 5 suggest that the resulting truncation errors do indeed decay exponentially fast as $M$ increases.

**REMARK 4.1 (PML parameters).** While it may be tempting to take the parameter $c$ large in (4.1) so as to improve the exponential decay of $\tilde{\varphi}$, consideration of discretization errors ($\tilde{\varphi}$ has to be approximated on a grid) will typically put bounds on the size of $c$. Even if we assume, as is the case in this paper, that the same grid size is used inside and outside the PML layer, the question of determining the “optimal” PML parameters is subtle and depends closely on the discretization scheme employed, for one seeks to balance truncation and discretization errors. In this paper, we simply take $c = 1$, and avoid questions related to an optimal choice of PML parameters.

**4.2. Singular and nearly-singular integration.** With a truncated curve $\Gamma_M$ at hand, we are now in position to discretize the integral operators in (4.4). To this end, we assume that the curve $\Gamma_M$ is covered by $N$ non-overlapping possibly curved elements $\{\mu_n\}_{n=1}^N$ of approximate size $h$; i.e.

$$ (4.7) \qquad \Gamma_M = \bigcup_{n=1}^N \bar{\mu}_n. $$

For each element $\mu_n$, a $C^\infty$ bijection $\chi_n : [0, 1] \to \bar{\mu}_n$ is assumed known.

On each $\mu_n$, a P-point Gauss-Legendre quadrature rule $Q_n = \{(\mathbf{x}_{p,m}, w_{p,m})\}_{p=1}^P$ is constructed by mapping the quadrature on the reference segment $[0, 1]$ to $\bar{\mu}_m$ using $\chi_m$. A composite quadrature $Q$ for $\Gamma_M$ is then obtained by collecting the element quadratures:

$$ (4.8) \qquad Q = \{Q_n\}_{n=1}^N = \left\{ \left\{ (\mathbf{x}_{p,n}, w_{p,n}) \right\}_{p=1}^P \right\}_{n=1}^N. $$

Given an element index $p$ and a local quadrature index $n$, it will be sometimes convenient to use a linear index $i = n(P-1)+p$, $1 \le i \le NP$, to refer to the node $\mathbf{x}_{p,n}$ as $\mathbf{x}_i$. With this linear index, the boundary quadrature can be alternatively written as $Q = \{(\mathbf{x}_i, w_i)\}_{i=1}^{NP}$.

Using a Nyström discretization, the integral operators in (3.20) are approximated using the quadrature $Q$, and equality is enforced (strongly) at the quadrature nodes. This yields the following discrete linear system for $\tilde{\varphi}_i \approx \tilde{\varphi}_M(\mathbf{x}_i)$:

$$ (4.9) \qquad -\frac{\tilde{\varphi}_i}{2} + \sum_{j=1}^{NP} K_{i,j} \tilde{\varphi}_j = \sum_{j=1}^{NP} G(\mathbf{x}_i, \mathbf{x}_j) f(\mathbf{x}_j) w_j, \quad 1 \le i \le NP, $$
where $K_{i,j} := K(\mathbf{x}_i, \mathbf{x}_j)w_j$ with

$$ (4.10) \qquad K(\mathbf{x}, \mathbf{y}) := \nabla_{\mathbf{y}} \tilde{G}(\mathbf{x}, \mathbf{y}) \cdot A(\mathbf{y})\mathbf{n}(\mathbf{y}) + \alpha(\mathbf{y})\tau'(\mathbf{y})\tilde{G}(\mathbf{x}, \mathbf{y}). $$

It is well-known that, for a given element $\mu_n \subset \Gamma_M$, the quadrature rule $Q_n$ need to be corrected for points $\mathbf{x}_i$ close to (or on) $\mu_n$ due to the singular nature of the kernels in the integral operators. In other words, the discrete problem (4.9), needs to be modified to properly account for the singular and nearly-singular interactions, where the Gauss-Legendre quadrature becomes inaccurate due the lack of smoothness of the integrand. In this paper, we opt for a simple correction method whereas $\tilde{\varphi}$ is replaced by its Lagrange interpolant on the element $\mu_n$, and the integration against the Lagrange interpolant is performed using an adaptive Gauss-Kronrod quadrature. In more detail, for a target point $\mathbf{x}_i$ “close” to the element $\mu_n$, where close is measured relative to the local mesh size, the matrix entries $K_{i,j_p}$ obtained from the regular quadrature are replaced by the following integrals

$$ (4.11) \qquad \tilde{K}_{i,j_p} = \int_{\mu_n} K(\mathbf{x}_i, \mathbf{y}) \ell_{n,p}(\mathbf{y}) \, d\mathbf{s}_{\mathbf{y}}, \quad 1 \le p \le P, $$

where $j_p = n(P-1) + p$ and $\ell_{n,p}$ is defined as $\ell_{n,p}(\chi_n(s)) = \hat{\ell}_p(s)$, $0 \le s \le 1$ with $\hat{\ell}_p$ being the $p$-th Lagrange basis over the reference $P$-point Gauss-Legendre quadrature $\{s_p, w_p\}_{p=1}^P$; i.e., $\hat{\ell}$ is the $P-1$ degree polynomial satisfying $\hat{\ell}_p(s_q) = \delta_{pq}$ for $1 \le q \le P$. In practice, the integrals in (4.11) are computed using an adaptive quadrature rule to a prescribed absolute tolerance; in all the examples presented in this paper, we set the integration tolerance to $10^{-10}$. The same correction approach is taken to accurately evaluate the right-hand side of (4.9) which involves the single-layer operator.

The errors introduced by the corrected quadrature approach for the singular and nearly-singular entries can be estimated as follows. First, using the definition of $\tilde{K}_{i,j_p}$ in (4.11) and assuming that $\tilde{\varphi}(\mathbf{x}_{j_p}) = \tilde{\varphi}_{j_p}$, the quadrature error for the element $\mu_n$ with a target point $\mathbf{x}_i$ can be written as

$$ (4.12) \qquad \int_{\mu_n} K(\mathbf{x}_i, \mathbf{y}) \tilde{\varphi}(\mathbf{y}) \, d\mathbf{s}_{\mathbf{y}} - \sum_{p=1}^{P} \tilde{\varphi}_{j_p} \tilde{K}_{i,j_p} = \int_{\mu_n} K(\mathbf{x}_i, \mathbf{y}) \left( \tilde{\varphi}(\mathbf{y}) - \sum_{p=1}^{P} \tilde{\varphi}_{j_p} \ell_{n,p}(\mathbf{y}) \right) \, d\mathbf{s}_{\mathbf{y}}. $$

From standard interpolation results (e.g. [38, Theorem 1]), we have that

$$ (4.13) \qquad \max_{\mathbf{y} \in \tau_m} \left| \tilde{\varphi}(\mathbf{y}) - \sum_{j=1}^{P} \tilde{\varphi}_{j_p} \ell_{n,j}(\mathbf{y}) \right| < C h^P, $$

for some constant $C$ which depends on the derivatives up to order $P$ of $\tilde{\varphi}$ and the element parametrization $\chi_n$. Provided $K(\mathbf{x}_i, \cdot)$ is absolutely integrable over $\mu_n$, as in all examples of this paper, we may then write the following bound:

$$ (4.14) \qquad \left| \int_{\mu_n} K(\mathbf{x}_i, \mathbf{y}) \left( \tilde{\varphi}(\mathbf{y}) - \sum_{p=1}^{P} \tilde{\varphi}_{j_p} \ell_{m,j}(\mathbf{y}) \right) \, d\mathbf{s}_{\mathbf{y}} \right| \le C h^P \int_{\mu_n} |K(\mathbf{x}_i, \mathbf{y})| \, d\mathbf{s}_{\mathbf{y}}. $$

It is now important to distinguish between the single and double-layer kernels making up $K$, as their singularities as $\mathbf{y} \to \mathbf{x}$ differ. The single-layer kernel has a
singularity of logarithmic type; this means, in particular, that its contribution to the integral on the right-hand side of (4.14) can be bounded by $h|\log(h)|$. The behaviour of the double-layer kernel is slightly more subtle. If $\mathbf{x}_i \in \mu_m$, the double-layer kernel is in fact bounded, and its contribution to the integral on the right-hand side of (4.14) can be bounded by $h$, yielding a logarithmic improvement over the single-layer error. If $\Gamma \ni \mathbf{x}_i \notin \mu_m$ with $\text{dist}(\mathbf{x}_i, \mu_m) \sim h$, however, two cases arise depending on the regularity of the curve $\Gamma_M$ (which is assumed to be at least Lipschitz). These are

$$ (4.15) \qquad \int_{\tau_m} \frac{(\tau(\mathbf{x}_i) - \tau(\mathbf{y})) \cdot \mathbf{n}(\mathbf{y})}{|\mathbf{x}_i - \mathbf{y}|^2} d\mathbf{s}_{\mathbf{y}} \sim \begin{cases} \mathcal{O}(h) & \text{if } \Gamma \text{ is } \mathcal{C}^2 \\ \mathcal{O}(1) & \text{otherwise} \end{cases} $$

Since $K$ in (4.10) involves both the single- and double-layer kernels, we expect the error in (4.14) to be proportional to $h^{P+1}|\log(h)|$ for smooth curves, or $h^P$ for curves featuring corners. The numerical experiments presented in next section, closely match these estimates.

**5. Numerical results.** In this section, we present a variety of numerical results designed to showcase the capabilities of the proposed methodology, as well as to test the effect of the parameters involved in the domain truncation and discretization. All linear systems are solved iteratively by means of GMRES [37], with a relative residual tolerance of $10^{-10}$.

In the examples that follow, the reported errors are measured at a finite set of test points $\mathbf{X}$ located either on the boundary $\Gamma$ (in which case they coincide with the quadrature nodes) or in the volume $\Omega$. In detail, errors are measured via

$$ (5.1) \qquad E_{M,h} = \max_{\mathbf{x} \in \mathbf{X}} |\varphi_{M,h}(\mathbf{x}) - \varphi_{\text{ref}}(\mathbf{x})|, $$

where $\varphi_{\text{ref}}$ denotes a reference solution which is the exact PDE solution when it is available or a solution obtained using a fine discretization of $\Gamma_M$ with a large truncation parameter $M$. When $\mathbf{x} = \mathbf{x}_j \in \mathbf{X} \cap \Gamma_M$ is a quadrature node, the value $\varphi_{M,h}(\mathbf{x})$ is taken from $\tilde{\varphi}_j$ which is in turn obtained by solving the underlying linear system. Otherwise, when $\mathbf{x} \in \mathbf{X} \cap \Omega$, we use the domain quadrature $\mathcal{Q}$ to evaluate an approximation to (3.10):

$$ (5.2) \qquad \varphi_{M,h}(\mathbf{x}) = \sum_{j=1}^{NP} K(\mathbf{x}, \mathbf{x}_j) \tilde{\varphi}_j - \tilde{G}(\mathbf{x}, \mathbf{x}_j) f(\mathbf{x}_j) w_j. $$

**5.1. Wave-maker.** The first example problem that we consider is that of a wavemaker whose geometry consists of the semi-infinite strip $[0, \infty) \times [-1, 0]$. For validation purposes, we construct an exact solution by summing the right-going propagating mode and the first evanescent mode in (2.6):

$$ (5.3) \qquad \varphi_{\text{ref}}(x_1, x_2) = \frac{\cosh(k(x_2 + 1))}{\cosh(k)} e^{ikx_1} + \frac{\cos(\gamma_1(x_2 + 1))}{\cos(\gamma_1)} e^{-\gamma_1 x_1}, $$

where $k$ and $\gamma_1$ are implicitly defined by (2.2) and (2.3), respectively.

The Neumann trace of $\varphi_{\text{ref}}$ is given as a source on the left boundary $\{0\} \times [-1, 0]$, and a numerical approximation $\tilde{\varphi}_{M,h}$ is found by solving the linear system (4.9). For this example, the error is measured at all the quadrature nodes of $\Gamma_M$; for points inside the PML, we use $\tilde{\varphi}_{\text{ref}}(x_1, x_2) = \varphi_{\text{ref}}(\tau(x_1, x_2), x_2)$ as the exact solution. The PML
function is given by (4.1), with the PML strength $c = 1$. The value of $a$, denoting the starting point of the PML, is indicated in each example.

Figure 2a displays the solution of the aforementioned problem on the free surface for the $\nu = d\omega^2/g = 1$, which yields $k \approx 1.2$, and $\gamma_1 \approx 2.8$. The PML layer starts at $a = 2\lambda$ (indicated by a dashed vertical line), where $\lambda = 2\pi/k$ is the surface wavelength. Since we are interested in assessing the effect of the truncation on the error in this example, we fix $h = \lambda/10$ and $P = 10$; i.e., the resulting domain discretization consists of 100 points-per-wavelength, which in this example suffices to achieve 10 digits of accuracy in the approximate BIE solution. A visual comparison between the exact solution, its analytical extension, and the numerical approximation $\tilde{\varphi}_{M,h}$ is shown in Figure 2a, where it is apparent that $\tilde{\varphi}_{M,h}$ converges to $\tilde{\varphi}$ on all of $\Gamma_M$, and therefore to $\varphi$ for $x_1 < a = 2\lambda$ (we only show here the values at $\Gamma_f$). For Figure 2a, the domain is truncated at $M = 3\lambda$. To study the rate of convergence as $M$ increases, we then present Figure 2b which displays the maximum error $E_{M,h}$, measured again on all the quadrature nodes, as a function of the truncation parameter $M$ for various values of $\nu$. (As in the previous example, we use a fine domain discretization to suppress the effect of BIE discretization errors.) We vary $M$ from $a$ to $a+5\lambda$ and display $E_{M,h}$ as a function of the “normalized” PML length $\ell = (M-a)/\lambda$. Given the exponential decay rate of $\varphi$ in (4.3), we expect a truncation error of $O(e^{-\min(k,\gamma_1)M})$ which can be more conveniently expressed in terms of the dimensionless PML length as $O(e^{-\min(2\pi,\gamma_1\lambda)\ell})$. The expected decay rate is clearly observed in Figure 2b, where for small $\nu$ the error is dominated by the propagating mode (and thus exhibits a $O(e^{-kM})$ rate), while for large values of $\nu$ we observe a slower exponential decay due to the evanescent mode (which decays as $O(e^{-\gamma_1 M})$). The plateau observed around $10^{-10}$ is an indication of errors stemming from the BIE solution (which include, for instance, the iterative solver residual and the quadrature rule errors) becoming dominant.

<figure><img src="image_3.png" /><figcaption>(a) Exact solution $\varphi_{\text{ref}}$, its analytic extension $\tilde{\varphi}_{\text{ref}}$, and the numerical approximation $\tilde{\varphi}_{\ell,h}$. The vertical dashed line indicates the start of the PML layer.</figcaption></figure>

<caption>(b) Convergence to exact solution as the truncation length $M$ is increased for different values of the parameter $\nu$. The dashed lines are the expected reference slopes.</caption>

Fig. 2: Solution to the wave-maker problem (left) and convergence for various values of $\nu$ (right).

As observed in Figure 2b, the exponential decay rate of truncation errors appears to deteriorate as $\nu$ increases. This is easy to understand if we consider $k$ and $\gamma_1$ as functions of $\nu$; both functions are plotted in Figure 3. As can be seen (and it easily follows from (2.2) and (2.3)) we have that $k \to \nu$, $\gamma_1 \to \pi/2$ in the limit $\nu \to \infty$. As
such, the truncation of the evanescent mode becomes dominant for sufficiently large $\nu$, and a large PML layer may be needed to obtain a small truncation error. While not a major drawback of the method for reasonably small $\nu$ values — the method still converges exponentially fast for any fixed $\nu$ — we propose next a possible modification to the PML function (4.1) to improve the truncation errors in the presence of slowly-decaying evanescent modes.

<figure><img src="image_4.png" /><figcaption>Fig. 3: Wavenumber $k$ versus $\gamma_1$ as a function of $\nu$.</figcaption></figure>

Inspired by [42], we improve our complex stretching to account not only for the propagating surface-wave mode but also for the evanescent modes, hence ensuring the effectiveness of the PML for a larger range of $\nu$ values. The main idea is to employ two different absorbing layers; a complex scaled layer, as before, to absorb the propagating mode, followed by a real stretching to accelerate the convergence of the evanescent modes. In detail, we replace the PML function (4.1) by:

$$ (5.4) \qquad \tau_{\text{two-layer}}(x_1) = \begin{cases} \tau(x_1), & |x_1| \le b, \\ \tau(\nu x_1) - \tau(\nu b) + \tau(b) & |x_1| > b, \end{cases} $$

where $b > a$ is a new parameter controlling the start of the real-stretching layer. The two-layered PML is depicted in Figure 4a. It is important to note that the real stretching increases the frequency of the oscillatory mode, and therefore should only be applied after the complex stretching layer.

Given the start of the PML layer $x_1 = a$, and the truncation parameter $M$, we chose $b$ by the following simple rationale. Supposing the propagating mode decays by $e^{-ck(b-a)}$, while the slowest decaying evanescent mode decays by $e^{-\gamma_1\nu(M-b)}$, requiring both errors to be equal yields

$$ (5.5) \qquad b = \frac{cka + \gamma_1\nu M}{ck + \gamma_1\nu}. $$

Equation (5.5) provides a way of choosing $b$ given the other problem parameters. Inserting this expression into the decay rate of the propagating mode, $e^{-ck(b-a)}$, and expressing the result in terms of the truncation parameter $M$, we obtain:

$$ (5.6) \qquad e^{-ck(b-a)} \sim \mathcal{O}(e^{-ck\mu M}), \quad \text{where} \quad \mu = \frac{\gamma_1\nu}{ck + \gamma_1\nu}. $$

Equation (5.6) provides the expected decay rate of $\tilde{\varphi}$ when the two-layer PML function (5.4) is used. The constant $\mu < 1$ in (5.6) provides an idea of the deterioration
compared to (4.3) for the cases where $ck < \gamma_1$ (i.e., cases for which the decay is dominated by the propagating mode). Note that, as $\nu \to \infty$, we have $\mu \to \pi/(2c+\pi)$, meaning the decay rate is at worst a constant factor slower than $e^{-ck\mu}$; for $c=1$, for example, we have $\mu \approx 0.61$ as $\nu \to \infty$. As can be seen in Figure 4b, the proposed idea significantly improves the truncation error, even for values of $\nu$ as large as 32. Furthermore, we observe that as $\nu$ increases the decay appears to match the expected $O(e^{-k\pi/(\pi+2)M})$ asymptote, shown as a dashed line in Figure 4b.

<figure><img src="image_5.png" /><figcaption>(a) Two-layered PML with $a=1, b=2, c=1$ (b) Convergence for large values of $\nu$ with $1, \nu=2$ (c) Convergence for large values of $\nu$ with $1, \nu=2$</figcaption></figure>

<caption>Fig. 4: Two-layered PML function (left) and its convergence properties (right)</caption>

Finally, although not the main focus of this paper, we present in Figure 5 convergence results with respect to mesh refinement. Following the discussion in Section 4.2, we expect a convergence of order $P$ for the wavemaker, which is what is observed Figure 5. (For this example we took a PML of length $4\lambda$, where $\lambda$ is the wavelength, which was enough to ensure sufficiently small truncation errors.)

**5.2. Scattering problem.** Having in the previous sections validated the proposed algorithm and its implementation, we here focus on the classical scattering problem. In the examples that follow we consider an incident field, corresponding to a right-going propagating mode, given by

$$ (5.7) \qquad \varphi^{\text{inc}}(x_1, x_2) = \frac{\cosh(k(x_2+1))}{\cosh(k)} e^{ikx_1}. $$

Expressing the total field as

$$ (5.8) \qquad \varphi^{\text{tot}} = \varphi^{\text{inc}} + \varphi, $$

we seek to find the scattered field $\varphi$ satisfying the radiation condition (2.1d). It follows from the problem linearity that $\varphi$ solves (2.1) with $f_1 = 0$ and $f_2 = -\nabla\varphi^{\text{inc}} \cdot \mathbf{n}$; i.e. $\varphi$ satisfies the following boundary condition:

$$ (5.9) \qquad \varphi = -\nabla\varphi^{\text{inc}} \cdot \mathbf{n}, \quad \mathbf{x} \in \Gamma_o \cup \Gamma_b. $$

We next consider three different settings for the scattering problem consisting of (i) a finite number of fully submerged jellyfish-like obstacles over a variable topography, (ii) a constant-depth topography with two piercing obstacles, and (iii) a
<figure><img src="image_6.png" /><figcaption>Fig. 5: Convergence to the modal solution for $\nu = 2$, and $P = 2, 5, 10$, as the number of points per wavelength is increased. The expected convergence rate of order $P$ described in Section 4.2, shown as dashed lines of the same color, is achieved.</figcaption></figure>

topography containing a step profile. In all of the examples considered in this section we use a mesh of size $\lambda/30$ and the parameter values $P = 5$, and $\nu = 4$ (which yields a wavelength $\lambda \approx 1.57$); the PML parameters are fixed as $a = 3\lambda$, $c = 1$, and the domain is truncated at $M = 4\lambda$.

**5.2.1. Submerged obstacles.** In our first scattering problem, we consider a bottom topography given by the following height function:

$$ (5.10) \qquad \text{topography}(s) = (s, -1 + e^{-s^2/4} \sin(2\pi s/\lambda)), \quad s \in \mathbb{R}, $$

which gives rise to a “bump-like” perturbation to a constant depth $^2$. In addition to the bump-like topography, we consider three jellyfish-like obstacles whose boundaries are given by the following parametrization (after appropriate rotation, translation, and scaling):

$$ (5.11) \qquad \text{jellyfish}(s) = r(s)(\cos(s), \sin(s)), \quad s \in \mathbb{R}, $$

where $r(s) = 1 + 0.3 \cos(4s + 2 \sin(s))$.

The resulting geometry is displayed in Figure 6 along with the incident, total, and scattered, fields. The damping of the scattered field inside the PML-layer is clearly visible in Figure 6 (bottom), especially for $x_1 > 2\lambda$ (i.e., passed the dashed lines). The main effect of the topography and obstacles is to cause a phase shift in the transmitted wave, as can be observed by comparing the incident field (top) to the total field (center) in Figure 6.

Because no exact solution is known for this scattering problem, we verified convergence of our Nyström method with respect to the mesh size by performing a self-convergence test, where we evaluate $\varphi$ at a uniform mesh of the rectangle $[1.5\lambda, 2.5\lambda] \times [-0.75 \times -0.25] \subset \Omega$, and compare the results to a reference solution obtained using 256 points-per-wavelength BIE domain discretization. We expect from the discussion in Section 4.2 the error to decrease as $O(h^6)$ as $h \to 0$ (this is a smooth geometry), which appears to agree with the observed convergence rate in Figure 8a.

<footnote>$^2$Although this is not a compact perturbation of a constant depth, the perturbations decay exponentially fast, and are numerically indistinguishable from a compact perturbation while being easier to construct.</footnote>
<figure><img src="image_7.png" /><figcaption>Fig. 6: Incident (top), total (middle), and scattered (bottom) fields for three jellyfish-like obstacles over a variable topography. The dashed lines indicate the start of the PML layer, after which we observe a decay of the scattered field (bottom) due to the complex stretching.</figcaption></figure>

**5.2.2. Piercing obstacles.** While we have insofar considered only fully submerged obstacles (i.e., the free surface coincides with $\mathbb{R} \times \{0\}$), piercing objects can be handled in a similar fashion. From an application viewpoint, piercing obstacles are particularly interesting since many relevant scatterers (ships, docks, buoys) lie precisely on the free surface. For simplicity, we take here a geometry given by two piercing semi-circles of radius $r = 1/4$ centered at $x_1 = \pm x_c$, with $x_c = r + 1/2$. This geometry creates a free surface with three disconnected components: $(-\infty, -1] \times \{0\}$, $[-0.5, 0.5] \times \{0\}$, and $[0.5, \infty) \times \{0\}$. As we will see later in Section 5.3, the part of the surface disconnected from infinity can act as a resonator for certain values of the Robin parameter $\nu$.

As for the other geometries considered above, we display in Figure 7 the incident, total, and scattered fields. Unlike the submerged obstacle, we observe a large scattered field both on the left and on the right of the obstacles. Again, as expected, the scattered field decays within the PML (bottom figure), both on the right and on the left of the obstacles.

Finally, we perform the same self-convergence test described in Section 5.2.1, where we evaluate the maximum error of the solution on a test domain inside of $\Omega$. The convergence with respect to the mesh refinement is shown in Figure 8b, where the reference line has the slope corresponding to $O(h^{-5})$. Since the piercing geometry features corners, we expect such behavior according to Section 4.2, although a detailed analysis of the discretization errors in such cases is likely to be significantly more challenging due to the abrupt switch in the boundary condition arising at the junction points connecting the piercing obstacle and the free surface.

**5.2.3. Step topography.** The final scattering geometry considered is that of a topography with a different depth as $x_1 \to \pm\infty$. In particular, we use here the following smooth step-like topography function:

$$ (5.12) \qquad \text{step}(s) = \left( s, -1 + \frac{\Delta d}{2} (1 + \tanh(5s)) \right), \quad s \in \mathbb{R}, $$

where, for this particular example, we take $\Delta d = 0.75$ so that the shallow region on the right-hand side of the bottom has a depth of 0.25. It is worth mentioning that this
<figure><img src="image_8.png" /><figcaption>Fig. 7: Incident (top), total (middle), and scattered (bottom) fields in the presence of two piercing bodies. The dashed vertical lines indicate the start of the PML layer.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8: Self-convergence of scattering problem for the geometries shown in Figure 6(a) and Figure 7(b). A reference solution with 256 points-per-wavelength was used when computing the error.</figcaption></figure>

is a challenging problem for methods based on the use of the standard water waves Green's function because the Neumann trace of the Green's function fails to vanish on an unbounded curve, making it necessary to introduce suitable approximations to achieve a BIE that can be discretized and numerically solved (see e.g. [35]).

This type of topography does not enter the framework introduced in Section 2 as it does not give rise to a computational domain that is a compact perturbation of a unit depth waveguide. Incidentally, the source term $f_2$ is not compactly supported, since the Neumann trace of the incident right-going wave does not vanish at $x_2 = -0.25$ for large values of $x_1$ (i.e., it does not satisfy the boundary condition on the right-hand side of the bottom boundary). Upon closer inspection, however, it is straightforward to see that both $\tilde{\varphi}^{\text{inc}}(\mathbf{x}) = \varphi^{\text{inc}}(\tau(\mathbf{x}))$ and its Neumann trace $\mathbf{n} \cdot \nabla \tilde{\varphi}^{\text{inc}}$, decay exponentially fast as $x_1 \to \infty$ (recall the definition of $\varphi^{\text{inc}}$ in (5.7)), implying that the analytic extension of $f_2$, although not compactly supported, decays exponentially fast as $x_1 \to$
∞. Although only a formal argument, the exponential decay of the source term means that we can expect our numerical method to behave in this case similarly to the compactly supported case.

As it turns out, this is indeed what we observe in Figure 9, where we display the incident, total, and scattered field for the step-like geometry. The decay of the scattered field $\tilde{\varphi}$ inside the PML layer becomes apparent (bottom image), and we see that for this value of the Robin parameter $\nu$ ($\nu = 4$ in this case), the scattered field is mostly concentrated on the shallow region. The change in the wavelength in the total field (middle image) is clearly visible on the right of the step, where the effective frequency is increased in the shallower region. While further theoretical and numerical investigations are needed to fully validate the effectiveness of the proposed technique for step-like topographies, we believe the results obtained are promising, in particular given the simplicity of the method.

<figure><img src="image_10.png" /><figcaption>Fig. 9: Incident (top), total (middle) and scattered (bottom) field for step-like bottom topography given by the curve parametrization in (5.12).</figcaption></figure>

**5.3. Resonant frequencies.** In this final set of examples, we consider the problem of finding resonant frequencies to the water-waves problem. More precisely, we seek to find the complex values of the Robin parameter $\nu$ such that (2.1) has a non-trivial solution for $f_1 = f_2 = 0$. Using the boundary integral formulation (3.18), this corresponds to finding eigenvalues $\nu$ and eigenfunctions $\tilde{\varphi}_\nu$ solving the following (generalized) linear eigenvalue problem:

$$ (5.13) \qquad -\frac{\tilde{\varphi}_\nu(\mathbf{x})}{2} + D_\Gamma[\tilde{\varphi}_\nu](\mathbf{x}) = -\nu S_{\Gamma_f}[\tau'\tilde{\varphi}_\nu](\mathbf{x}), \quad \mathbf{x} \in \Gamma. $$

Note that we do not allow for $\tau$ to depend on $\nu$, which would make for a nonlinear eigenvalue problem. In what follows, we fix the PML parameters to be $a = 2.5$, $c = 1$, and take $M = 5$; for the discretization parameters, we set $P = 5$ and take a mesh size of $h = 0.2$.

Figure 10 displays the spectrum corresponding to the piercing geometry of Section 5.2.2, as well as the eigenfunction corresponding to an eigenvalue with a small imaginary part. The eigenvalues on a large region of the complex plane enclosing $\{z \in \mathbb{C} : 0 \le \operatorname{Re} z \le 20, -20 \le \operatorname{Im} z \le 0\}$ are shown on the top-left figure. There we can clearly observe that a part of the spectrum follows a straight line of approximate slope $-1$; these eigenvalues are typically PML-dependent and correspond to
the complex stretching of the continuous spectrum $\mathbb{R}^+$ in the original problem. The bifurcation at $Re \nu \approx 10$ has also been observed in [41, e.g. Figures 5,6]. The more interesting and physically relevant eigenvalues appear at nearly regular intervals close to the real axis. These correspond to eigenfunctions that are localized near the obstacles, and are not noticeable affected by variations in the PML parameters. We show in the top-right image a zoomed-in version of the spectrum, and mark in red one of the aforementioned eigenvalues with a small imaginary part. Finally, the bottom image displays the (normalized) eigenfunction corresponding to the highlighted red eigenvalue. We see that portion of the free surface between the two-obstacles acts as a trapping region, where localized eigenfunctions can exist.

<figure><img src="image_11.png" /><figcaption>Fig. 10: Complex resonances of a piercing geometry with two semi-circles. The top-left image shows the spectrum on a large portion of the complex plane, while the top-right image provides a zoomed-in version where we can clearly see the presence of an “almost real” eigenvalue, depicted in red. The bottom image displays the eigenfunction corresponding to the red eigenvalue in the top-right image.</figcaption></figure>

**6. Open theoretical problems.** This work has opened up several challenging theoretical questions that were not addressed in this paper. In what follows we mention a few of them.

Firstly, the logarithmic growth of complex-scaled Green function (3.7) makes it difficult to pose a functional framework for analyzing the well-posedness of (3.20) in the usual Sobolev spaces. A possible way to circumvent such theoretical difficulties is to consider instead a Green’s function $\tilde{G}^D$ satisfying a (homogeneous) Dirichlet boundary condition at $x_2 = H$ for some artificial height $H > 0$ [36]. By the method of images, we have that $\tilde{G}^D$ is given by

$$ (6.1) \qquad \tilde{G}^D(\mathbf{x}, \mathbf{y}) = -\frac{1}{4\pi} \log ((\tau(x_1) - \tau(y_1))^2 + (x_2 - y_2)^2) \\ \qquad \qquad + \frac{1}{4\pi} \log ((\tau(x_1) - \tau(y_1))^2 + (x_2 - y_2^r)^2) $$

where $y_2^r = 2H - y_2$ is the image point. It is possible to verify that, as in the usual Laplace case, this complexified Dirichlet Green’s function decays linearly as $|y| \to \infty$, which in turn may help define an appropriate functional framework for equations like (3.20).

Secondly, the equivalence between the boundary integral equation (3.20) and the complex-scaled water-waves problem (3.5) was not addressed. To briefly indicate the source of difficulty, let us consider a density $\sigma : \Gamma \to \mathbb{C}$ solution of (3.20), and define
the function $\tilde{\varphi}: \Omega \to \mathbb{C}$ by

$$ (6.2) \qquad \tilde{\varphi}(\mathbf{x}) = \mathcal{D}_{\Gamma}[\sigma](\mathbf{x}) - S_{\Gamma}[f - \alpha\tau'\sigma](\mathbf{x}). $$

One would like to show that (6.2) is a solution to the original problem (3.5). The fact that (6.2) satisfies the PDE (3.5a) follows directly from the choice of the integral representation. This means that we may write

$$ (6.3) \qquad \tilde{\varphi}(\mathbf{x}) = \mathcal{D}_{\Gamma}[\tilde{\varphi}](\mathbf{x}) - S_{\Gamma}[\tau' \nabla \tilde{\varphi} \cdot \mathbf{n}](\mathbf{x}). $$

Taking the Dirichlet trace of (6.2) and using the fact that $\sigma$ solves (3.20) we obtain that $\tilde{\varphi} = \sigma$ on $\Gamma$. Subtracting (6.3) from (6.2) yields

$$ (6.4) \qquad S_{\Gamma}[\tau' \nabla \tilde{\varphi} \cdot \mathbf{n} - f + \tau' \alpha \tilde{\varphi}](\mathbf{x}) = 0, \quad \mathbf{x} \in \Omega. $$

Equality (6.4) can be extended to $\Gamma$ to give:

$$ (6.5) \qquad S_{\Gamma}[\tau' \nabla \tilde{\varphi} \cdot \mathbf{n} - f + \tau' \alpha \tilde{\varphi}](\mathbf{x}) = 0, \quad \mathbf{x} \in \Gamma. $$

As a consequence, $\tilde{\varphi}$ as given by (6.2) solves the water-waves system (3.5) if the single-layer operator $S_{\Gamma}$ is “injective” in a certain sense. Due to the current limited understanding of the mapping properties of the single-layer operator in this scenario, we are unable to provide a more precise analysis at this time.

Finally, while the use of a Dirichlet Green's function was motivated above for theoretical reasons, the use of a Green's function satisfying a homogenous Neumann condition at the flat bottom $x_2 = -1$ is more natural from a numerical point of view since it eliminates many degrees of freedom at the bottom (i.e. only the perturbations to the topography need to be discretized when $\Omega$ lies above $x_2 = -1$). Such a Green's function is easy to construct once again by the method of images:

$$ (6.6) \qquad \tilde{G}^N(\mathbf{x}, \mathbf{y}) = -\frac{1}{4\pi} \log \left( (\tau(x_1) - \tau(y_1))^2 + (x_2 - y_2)^2 \right) \\ -\frac{1}{4\pi} \log \left( (\tau(x_1) - \tau(y_1))^2 + (x_2 - y_2^d)^2 \right) $$

where $y_2^d = -2 - y_1$. Unlike the Dirichlet Green's function, however, $\tilde{G}^N$ still possesses a logarithmic growth at infinity. Because of the aforementioned computational advantage of such a Neumann Green's function, it would be interesting to better understand the properties of the integral equations derived from it.

**7. Conclusions and future directions.** In this paper, we presented a novel boundary integral formulation for solving the two-dimensional finite depth water waves problem. We showed that the underlying truncation methodology leads to exponentially small truncation errors, and it is applicable to many challenging problems including scattering in the presence of a step topography. We believe the ability to compute resonances (both real and complex) in the form of a linear eigenvalue problem, while enjoying the dimensional reduction of boundary integral equations, is a very interesting feature of the method worthy of further exploration. Future directions of research include extensions to infinite depth and the three-dimensional problem.

**Appendix A. Complex-scaled fundamental solution.** Let $x_1 \in \mathbb{R} \mapsto \tau(x_1) \in \mathbb{C}$ be as in section 3.1 a continuous and piecewise continuously differentiable function satisfying (3.4d). In this section, we provide a direct proof that

$$ \tilde{G}(\mathbf{x}, \mathbf{y}) = G(\tau(\mathbf{x}), \tau(\mathbf{y})) \text{ where } \tau(\mathbf{x}) = (\tau(x_1), x_2) $$
is the free-space fundamental solution of the complex-scaled Laplace equation, i.e. that $\tilde{G}$ satisfies

$$ (A.1) \qquad \int_{\mathbb{R}^2} \nabla_y \tilde{G}(x, y) \cdot A(y) \nabla v(y) \, dy = v(x), \quad \forall x \in \mathbb{R}^2, $$

for all test functions $v \in \mathcal{D}(\mathbb{R}^2)$, where $A$ is given by (3.5).

By continuity arguments, it suffices to prove (A.1) for points $x$ where $\tau$ is differentiable. Note that $\tilde{G}(x, y)$ is well-defined for all $x, y$ in $\mathbb{R}^2$, as soon as $x \neq y$. Indeed, suppose that

$$ (\tau(x_1) - \tau(y_1))^2 + (x_2 - y_2)^2 = 0, $$

then

$$ \mathrm{Re}(\tau(x_1) - \tau(y_1)) = 0, $$

which implies by (3.4d) that $x_1 = y_1$, and finally that $x = y$. Now by the properties of the Lebesgue integral and the symmetry of $A(y)$, it holds that, by letting $I(x)$ denote the left hand-side of (A.1),

$$ I(x) = \lim_{\epsilon \to 0} \int_{|y-x|>\epsilon} A(y) \nabla_y \tilde{G}(x, y) \cdot \nabla v(y) \, dy. $$

Integrating by parts and using that

$$ \nabla_y \cdot (A(y) \nabla_y \tilde{G}(x, y)) = 0 \quad \text{if } y \neq x, $$

we obtain that

$$ I(x) = \lim_{\epsilon \to 0} \frac{1}{\epsilon} \int_{|y-x|=\epsilon} A(y) \nabla_y \tilde{G}(x, y) \cdot (x-y) v(y) \, ds_y. $$

On the other hand, the chain rule and the definition of the function $G$ lead to

$$ \nabla_y \tilde{G}(x, y) = \frac{1}{2\pi} \frac{1}{(\tau(x_1) - \tau(y_1))^2 + (x_2 - y_2)^2} \left[ \frac{\tau'(y_1)(\tau(x_1) - \tau(y_1))}{x_2 - y_2} \right], $$

so that by definition of $A(y)$:

$$ I(x) = \lim_{\epsilon \to 0} \frac{1}{2\pi\epsilon} \int_{|y-x|=\epsilon} \tau'(y_1) \frac{\tau'(y_1)^{-1}(\tau(x_1) - \tau(y_1))(x_1 - y_1) + (x_2 - y_2)^2}{(\tau(x_1) - \tau(y_1))^2 + (x_2 - y_2)^2} v(y) \, ds_y. $$

Now parametrizing the circle by $y = x + \epsilon r(\theta)$, with $r(\theta) = (\cos(\theta), \sin(\theta))$ and using straightforward Taylor expansions, this leads to

$$ I(x) = \frac{v(x)}{2\pi} \int_{0}^{2\pi} \frac{\tau'(x_1)}{\tau'(x_1)^2 \cos(\theta)^2 + \sin(\theta)^2} d\theta. $$

We conclude by using the following classical result of complex analysis:

$$ \int_{0}^{2\pi} \frac{1}{\alpha^2 \cos(\theta)^2 + \sin(\theta)^2} d\theta = \frac{2\pi}{\alpha} $$

which holds for all $\alpha \in \mathbb{C}$ such that $\mathrm{Re}(\alpha) > 0$. The result can be proved first for real $\alpha$, rewriting the identity as

$$ \mathrm{Im} \left( \int_C \frac{dz}{z} \right) = 2\pi $$

where $C$ is an ellipse defined by $C = \{\alpha \cos(\theta) + i \sin(\theta); 0 \le \theta \le 2\pi\}$. Then one concludes that the identity still holds for complex $\alpha$ by analyticity.
[1] B. K. ALPERT, *Hybrid Gauss-trapezoidal quadrature rules*, SIAM J. Sci. Comput., 20 (1999), pp. 1551-1584.

[2] T. ANGELL, G. HISAO, AND R. KLEINMAN, *Mathematical Approaches in Hydrodynamics*, SIAM, 1991, ch. Recent developments in floating body problems.

[3] T. ANGELL, G. HISAO, AND R. KLEINMAN, *An integral equation for the floating-body problem*, J. Fluid Mech., 166 (1986), pp. 161-171.

[4] G. ATHANASSOULIS AND K. BELIBASSAKIS, *A consistent coupled-mode theory for the propagation of small-amplitude water waves over variable bathymetry regions*, J. Fluid Mech., 389 (1999), pp. 275-301.

[5] K. BELIBASSAKIS AND G. ATHANASSOULIS, *Three-dimensional Green's function for harmonic water waves over a bottom topography with different depths at infinity*, J. Fluid Mech., 510 (2004), pp. 267-302.

[6] J. L. BLACK, *Wave forces on vertical axisymmetric bodies*, J. Fluid Mech., 67 (1975), pp. 369-376.

[7] O. P. BRUNO, E. GARZA, AND C. PÉREZ-ARANCIBIA, *Windowed Green function method for nonuniform open-waveguide problems*, IEEE Trans. Antennas Propag., 65 (2017), pp. 4684-4692.

[8] O. P. BRUNO, M. LYON, C. PÉREZ-ARANCIBIA, AND C. TURC, *Windowed green function method for layered-media scattering*, SIAM J. Appl. Math., 76 (2016), pp. 1871-1898.

[9] O. P. BRUNO AND C. PÉREZ-ARANCIBIA, *Windowed Green function method for the helmholtz equation in the presence of multiply layered media*, Proc. R. Soc. A: Math. Phys. Eng. Sci., 473 (2017), p. 20170161.

[10] A. L. B. CAUCHY, *Theory of wave propagation à la surface of a heavy fluid of indefinite depth*, Royal Academy of Sciences, 1816.

[11] S. K. CHAKRABARTI, *Application and verification of deepwater green function for water waves*, J. Sh. Res., 45 (2001), pp. 187-196.

[12] S. N. CHANDLER-WILDE AND B. ZHANG, *A uniqueness result for scattering by infinite rough surfaces*, SIAM J. Appl. Math., 58 (1998), pp. 1774-1790.

[13] D. L. COLTON, R. KRESS, AND R. KRESS, *Inverse Acoustic and Electromagnetic Scattering Theory*, vol. 93, Springer, 1998.

[14] L. M. FARIA, C. PÉREZ-ARANCIBIA, AND M. BONNET, *General-purpose kernel regularization of boundary integral equations via density interpolation*, Comput. Methods Appl. Mech. Eng., 378 (2021), p. 113703.

[15] J. FENTON, *Wave forces on vertical bodies of revolution*, J. Fluid Mech., 85 (1978), pp. 241-255.

[16] L. GREENGARD AND V. ROKHLIN, *A fast algorithm for particle simulations*, J. Comput. Phys., 73 (1987), pp. 325-348.

[17] C. HAZARD AND M. LENOIR, *Determination of scattering frequencies for an elastic floating body*, SIAM J. Math. Anal., 24 (1993), pp. 1458-1514.

[18] R. HEIN, *Green's Functions and Integral Equations for the Laplace and Helmholtz Operators in Impedance Half-Spaces*, PhD thesis, Ecole Polytechnique, 2010.

[19] R. HEIN, M. DURÁN, AND J.-C. NÉDÉLEC, *Explicit representation for the infinite-depth two-dimensional free-surface Green's function in linear water-wave theory*, SIAM J. Appl. Math., 70 (2010), pp. 2353-2372.

[20] N. KUZNETSOV, V. MAZ'Å, AND B. VAINBERG, *Linear Water Waves: A Mathematical Approach*, Cambridge University Press, 2002.

[21] M. LENOIR AND A. TOUNSI, *The localized finite element method and its application to the two-dimensional sea-keeping problem*, SIAM J. Numer. Anal., 25 (1988), pp. 729-752.

[22] C. LINTON AND P. MCIVER, *Mathematical techniques for wave/structure interactions*, 2001.

[23] C. M. LINTON AND P. MCIVER, *Handbook of Mathematical Techniques for Wave/Structure Interactions*, Chapman and Hall/CRC, 2001.

[24] W. LU AND Y. Y. LU, *Efficient high order waveguide mode solvers based on boundary integral equations*, J. Comput. Phys., 272 (2014), pp. 507-525.

[25] W. LU, Y. Y. LU, AND J. QIAN, *Perfectly matched layer boundary integral equation method for wave scattering in a layered medium*, SIAM J. Appl. Math., 78 (2018), pp. 246-265.

[26] E. MACKAY, *The green function for diffraction and radiation of regular waves by two-dimensional structures*, European Journal of Mechanics-B/Fluids, 87 (2021), pp. 151-160.

[27] S. A. MAIER, *Plasmonics: Fundamentals and Applications*, vol. 1, Springer, 2007.

[28] P. A. MARTIN, *Multiple scattering: interaction of time-harmonic waves with N obstacles*, Cambridge University Press, 2006.

[29] M. MCIVER, *An example of non-uniqueness in the two-dimensional linear water wave problem*,
J. Fluid Mech., 315 (1996), pp. 257–266.

[30] W. MCLEAN AND W. C. H. MCLEAN, *Strongly elliptic systems and boundary integral equations*, Cambridge university press, 2000.

[31] C. C. MEI, *Numerical methods in water-wave diffraction and radiation*, Annu. Rev. Fluid Mech., 10 (1978), pp. 393–416.

[32] J. N. NEWMAN, *Algorithms for the free-surface green function*, J. Eng. Math., 19 (1985), pp. 57–67.

[33] C. PÉREZ-ARANCIBIA, L. M. FARIA, AND C. TURC, *Harmonic density interpolation methods for high-order evaluation of Laplace layer potentials in 2D and 3D*, J. Comput. Phys., 376 (2019), pp. 411–434.

[34] C. PÉREZ-ARANCIBIA, P. RAMACIOTTI, R. HEIN, AND M. DURÁN, *Fast multipole boundary element method for the Laplace equation in a locally perturbed half-plane with a Robin boundary condition*, Comput. Methods Appl. Mech. Eng., 233 (2012), pp. 152–163.

[35] R. PORTER AND D. PORTER, *Water wave scattering by a step of arbitrary profile*, J. Fluid Mech., 411 (2000), pp. 131–164.

[36] M. D. PRESTON, P. G. CHAMBERLAIN, AND S. N. CHANDLER-WILDE, *An integral equation method for a boundary value problem arising in unsteady water wave problems*, J. Integral Equ. Appl., (2008), pp. 121–152.

[37] Y. SAAD AND M. H. SCHULTZ, *Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems*, SIAM J Sci Comput, 7 (1986), pp. 856–869.

[38] G. STRANG, *Approximation in the finite element method*, Numer. Math., 19 (1972), pp. 81–98.

[39] T. STRAUSZER-CAUSSADE, L. M. FARIA, A. FERNANDEZ-LADO, AND C. PÉREZ-ARANCIBIA, *Win-dowed Green function method for wave scattering by periodic arrays of 2D obstacles*, Stud. Appl. Math., 150 (2023), pp. 277–315.

[40] J. THOMAS BEALE, *Eigenfunction expansions for objects floating in an open sea*, Commun. Pure Appl. Math., 30 (1977), pp. 283–313.

[41] F. TREYSSEDE, K.-L. NGUYEN, A.-S. BONNET-BENDHIA, AND C. HAZARD, *Finite element com-putation of trapped and leaky elastic waves in open stratified waveguides*, Wave Motion, 51 (2014), pp. 1093–1107.

[42] B. VIAL, F. ZOLLA, A. NICOLET, M. COMMANDRÉ, AND S. TISSERAND, *Adaptive perfectly matched layer for wood's anomalies in diffraction gratings*, Opt. Express, 20 (2012), pp. 28094–28105.

[43] C. XIE, Y. CHOI, F. RONGÈRE, A. H. CLÉMENT, G. DELHOMMEAU, AND A. BABARIT, *Comparison of existing methods for the calculation of the infinite water depth free-surface green function for the wave-structure interaction problem*, Appl. Ocean Res., 81 (2018), pp. 150–163.

[44] R. W. YEUNG, *Numerical methods in free-surface flows*, Annu. Rev. Fluid Mech., 14 (1982), pp. 395–442.