1. Introduction
In this short paper, we consider the classical problem of gravity waves propagating at the (non-overturning) free surface of a homogeneous non-viscous fluid in irrotational motion over an impermeable (uneven but non-overturning) seabed. Mathematically, in two dimensions without obstacles (i.e. for a simply connected fluid domain extending to infinity in all horizontal directions), this leads to the system of equations (for 
$x\in \mathbb {R}$, 
$t\in \mathbb {R}$ or 
$t\geqslant t_0$) (Wehausen & Laitone Reference Wehausen and Laitone1960)
\begin{gather} \partial_x^{2}\phi+\partial_y^{2}\phi=0 \quad\mathrm{for}\ {-}d(x)\leqslant y\leqslant\eta(x,t), \end{gather}
\begin{gather}\partial_y\phi+ (\partial_xd)(\partial_x\phi)=0 \quad\mathrm{at}\ y={-}d(x), \end{gather}
\begin{gather}\partial_y\phi-\partial_t\eta-(\partial_x\eta)(\partial_{x}\phi)=0 \quad\mathrm{at}\ y= \eta(x,t), \end{gather}
\begin{gather}\partial_t\phi+g\eta+{{\textstyle{1\over 2}}}(\partial_x\phi)^{2}+ {{\textstyle{1\over 2}}}(\partial_y\phi)^{2}=0\quad \mathrm{at}\ y=\eta(x,t), \end{gather}
where 
$\phi (x,y,t)$ is a velocity potential such that 
$u\stackrel {\textrm def}{=}\partial _x\phi$ is the horizontal velocity and 
$v\stackrel {\textrm def}{=}\partial _y\phi$ is the vertical one, 
$g>0$ is the acceleration due to gravity (directed downwards), with 
$(x,y)$ respectively the horizontal and upward vertical Cartesian coordinates, and 
$t$ is the time. Here 
$y=\eta (x,t)$, 
$y=0$ and 
$y=-d(x)$ are, respectively, the equations of the free surface, of the still-water level and of the bottom; and 
$h(x,t)\stackrel {\textrm def}{=}\eta (x,t)+d(x)$ is the total water depth. Physically, (1.1) means that the motion is irrotational and isochoric, (1.2) and (1.3) characterise the impermeability of the bottom and of the free surface, while (1.4) expresses that the pressure at the free surface equals the constant atmospheric pressure (set to zero without loss of generality). Capillarity and other surface effects can be considered but they do not affect the analysis below, so they are of no interest here. Also, extensions of (1.1)–(1.4) in higher dimensions and/or moving bottoms are straightforward; these generalisations are considered at the end of the present paper. However, further generalisations (e.g. overturning surface and/or bottom, submerged obstacles, floating bodies, lateral solid boundaries, rough bottom) are beyond the scope of the present study; they require ad hoc investigations.
A Dirichlet–Neumann (or Dirichlet-to-Neumann) operator (DNO) takes as input a function expressed at a point of the domain boundary and outputs its (outward) normal derivative at the same point. Here, the DNO producing the (non-unitary outgoing) normal derivative at the free surface is 
$G({{\phi _{s}}})\stackrel {\textrm def}{=}[\partial _y\phi -(\partial _x\eta ) (\partial _x\phi )]_{y=\eta }$, where 
${{\phi _{s}}}(x,t)\stackrel {\textrm def}{=}\phi (x,\eta,t)$ denotes the velocity potential at the free surface. Fulfilling the Laplace equation (1.1) and the bottom impermeability condition (1.2), the DNO is a homogeneous linear function of 
${{\phi _{s}}}$, i.e. 
$G({{\phi _{s}}})=\mathscr {G}{{\phi _{s}}}$, where 
$\mathscr {G}$ is a self-adjoint positive-definite pseudo-differential operator depending nonlinearly on 
$\eta$ and 
$d$ (Craig & Sulem Reference Craig and Sulem1993; Craig et al. Reference Craig, Guyenne, Nicholls and Sulem2005). The operator 
$\mathscr {G}$ is a fundamental mathematical object because it ‘encodes’ the domain geometry, the kinematics of the fluid motion and the bottom impermeability; moreover, it appears explicitly in the Hamiltonian formulation (Zakharov Reference Zakharov1968) of the equations (1.1)–(1.4). Understandably, 
$\mathscr {G}$ has been the subject of many mathematical studies – see Lannes (Reference Lannes2013) and Nicholls & Reitich (Reference Nicholls and Reitich2001) for details – and it is at the heart of several rigorous investigations on water waves (e.g. Alazard, Burq & Zuily Reference Alazard, Burq and Zuily2012; Alazard & Baldi Reference Alazard and Baldi2015). Knowledge of the mathematical features of the DNO is certainly important, but its explicit construction is at least as important, in particular for practical applications.
For a flat horizontal free surface and bottom, the fluid domain is a strip and the DNO is easily obtained analytically, e.g. via a Fourier transform. For a wavy surface and bottom, the DNO can be constructed as a perturbation of the strip, assuming small amplitudes. This is the route followed in two dimensions by Craig & Sulem (Reference Craig and Sulem1993) and in three dimensions by Craig & Groves (Reference Craig and Groves1994) for flat seabeds, then extended to varying bottoms (Craig et al. Reference Craig, Guyenne, Nicholls and Sulem2005); these authors provided recurrence relations for computing the DNO to an arbitrary order of their perturbative expansion. For small perturbations of the flat surface and seabed, other series representations of the DNO are available in the literature (Dommermuth & Yue Reference Dommermuth and Yue1987; West et al. Reference West, Brueckner, Janda, Milder and Milton1987). Although all these series are formally equivalent, this is not necessarily the case with their truncations at the same order, as outlined by Schäffer (Reference Schäffer2008). Moreover, such expansions are badly conditioned, so prone to large numerical errors and instabilities (Wilkening & Vasan Reference Wilkening and Vasan2015). An explicit formulation of the DNO in expected to facilitate various reformulations for more efficient computations, for example, but this is not the scope of the present paper.
The main purpose of this paper is to show how explicit DNOs can be derived and, via a few examples, to show their interest for analytic manipulations. Although some indications on potential issues and remedies with numerical computations are briefly discussed, it is not the purpose here to derive the most effective way to compute numerically a DNO.
The paper is organised as follow. In § 2, an explicit DNO is derived in two dimensions via rather elementary algebra. This DNO being in complex form, a real reformulation is introduced in § 3 in order to facilitate analytical approximations. Some approximations for small amplitudes in finite depth and for finite amplitudes in shallow water are then derived in § 4. The DNO is extended to higher dimensions in § 5, and its generalisation for moving bottoms is provided in § 6. Finally, a summary and perspectives are briefly drawn in § 7.
2. Two-dimensional Dirichlet–Neumann operator
Let 
$\psi$ be the streamfunction harmonic conjugate of the velocity potential 
$\phi$ (Milne-Thomson Reference Milne-Thomson2011). These two functions are related by the Cauchy–Riemann relations 
$\phi _x =\psi _y=u$ and 
$\phi _y=-\psi _x=v$. Thus, the complex potential 
$f\stackrel {\textrm def}{=}\phi +\mathrm {i}\psi$ is a holomorphic function of 
$z\stackrel {\textrm def}{=} x+\mathrm {i} y$, with 
$z={{z_{s}}}\stackrel {\textrm def}{=} x+\mathrm {i}\eta$ at the free surface and 
$z={{z_{b}}}\stackrel {\textrm def}{=} x-\mathrm {i} d$ at the bottom. (As general notation, subscripts ‘
$s$’ and ‘
$b$’ denote quantities written, respectively, at the free surface and at the bottom.) The seabed being impermeable and static, it is a streamline where 
$\psi ={\psi _{b}}$ is constant. Without loss of generality, we then choose 
${\psi _{b}}=0$ for simplicity.
For any complex abscissa 
$z_0$, the Taylor expansion around 
$z_0=0$ is (omitting temporal dependences for brevity)
\begin{equation} f(z-z_0)=\exp[{-}z_0\partial_z]f(z)\ \stackrel{\textrm def}{=} \sum_{n=0}^{\infty}\,\frac{({-}1)^{n}z_0^{n}}{n!}\,\frac{\partial^{n}f(z)}{\partial z^{n}}. \end{equation}
For instance, taking 
$z_0=\mathrm {i} h=\mathrm {i}(d+\eta )$, the relation (2.1) written at the free surface becomes
\begin{equation} f({{z_{s}}}-\mathrm{i} h)=\exp[-\mathrm{i} h\partial_{{{z_{s}}}}]f({{z_{s}}}), \end{equation}
with the formal operator 
$\exp [-\mathrm {i} h \partial _{{{z_{s}}}}]\stackrel {\textrm def}{=}\sum _{n=0}^{\infty } (n!)^{-1}(-\mathrm {i} h)^{n}\partial _{{{z_{s}}}}^{n}$ together with 
$\partial _{{{z_{s}}}} \stackrel {\textrm def}{=}(1+\mathrm {i}\eta _x)^{-1}\partial _x$, 
$\partial _{{{z_{s}}}}^{2}=(1+\mathrm {i}\eta _x)^{-1} \partial _x(1+\mathrm {i}\eta _x)^{-1}\partial _x$, etc. (Throughout this paper, we use the classical convention that any operator acts on everything it multiplies on its right, unless parentheses enforce otherwise.)
It should be noticed that exponents denote differential compositions, so 
$h^{n}$ is the 
$n$th power of the function 
$h=d+\eta$, while 
$\partial _{{{z_{s}}}}^{n}$ is the 
$n$th iteration of the differential operator 
$\partial _{{{z_{s}}}}$. Therefore, for example, if 
$h$ is not constant then 
$h^{2}=h(x)^{2}\neq h(h(x))$, 
$h^{2}\partial _{{{z_{s}}}}^{2} \neq (h\partial _{{{z_{s}}}})^{2}=h\partial _{{{z_{s}}}}h\partial _{{{z_{s}}}}$, 
$\exp [-\mathrm {i} h\partial _{{{z_{s}}}}]\neq \sum _{n=0}^{\infty }(n!)^{-1}(-\mathrm {i} h\partial _{ {{z_{s}}}})^{n}$ and the operator inverse of 
$\exp [-\mathrm {i} h\partial _{{{z_{s}}}}]$ is not 
$\exp [\mathrm {i} h\partial _{{{z_{s}}}}]$ (but 
$\exp [\mathrm {i} h\partial _{{{z_{b}}}} ]$ as shown in § 3).
Since 
${{z_{s}}}-\mathrm {i} h=x-\mathrm {i} d={{z_{b}}}$ then 
$f({{z_{s}}}-\mathrm {i} h)=f({{z_{b}}}) ={{\phi _{b}}}$ is real (recall that 
${\psi _{b}}=\mathrm {Im}\,{f_{b}}=0$ by definition), while 
$f({{z_{s}}})={{\phi _{s}}}+\mathrm {i}{{\psi _{s}}}$ is complex. Therefore, the imaginary part of (2.2), i.e.
\begin{equation} 0 = \mathrm{Re} \{\exp [-\mathrm{i} h\partial_{{{z_{s}}}}] \} {{\psi_{s}}} + \mathrm{Im} \{\exp [-\mathrm{i} h\partial_{{{z_{s}}}}] \} {{\phi_{s}}}, \end{equation}yields at once
\begin{equation} {{\psi_{s}}} ={-} \left(\mathrm{Re} \{ \exp [-\mathrm{i} h\partial_{{{z_{s}}}}] \}\right)^{{-}1} \mathrm{Im} \{\exp [-\mathrm{i} h\partial_{{{z_{s}}}}] \} {{\phi_{s}}}. \end{equation}
The equation for the free-surface impermeability being 
$\partial _t\eta =\mathscr {G}{{\phi _{s}}}=-\partial _x{{\psi _{s}}}={{v_{s}}}-{{u_{s}}}\partial _x\eta$, an explicit definition of the DNO is obtained directly from (2.4) as
\begin{equation} \mathscr{G} = \partial_x\left(\mathrm{Re} \{ \exp [-\mathrm{i} h \partial_{{{z_{s}}}} ] \}\right)^{{-}1} \mathrm{Im} \{\exp [-\mathrm{i} h\partial_{{{z_{s}}}} ] \}. \end{equation}
Formula (2.5) provides an explicit expression for the DNO, i.e. 
$\mathscr {G}$ appears only on the left-hand side. It is the main result of this paper, which can be generalised to higher dimensions and for moving bottoms (see below). It is also suitable to derive various approximations, in particular high-order shallow-water approximations without assuming small amplitudes (see § 4.3 below; in fact, this goal was the original motivation for deriving (2.5)).
For applications, it is convenient to introduce an operator 
$\mathscr {J}$ such that 
$\mathscr {G}=-\partial _x\mathscr {J}\partial _x$, so
\begin{equation} \mathscr{J} ={-}\left(\mathrm{Re}\{\exp[- \mathrm{i} h \partial_{{{z_{s}}}}]\}\right)^{{-}1} \mathrm{Im}\{ \exp[-\mathrm{i} h\partial_{{{z_{s}}}}]\}\partial_x^{{-}1}. \end{equation}
Since 
$\mathscr {G}$ is a self-adjoint positive-definite operator (Lannes Reference Lannes2013), so is 
$\mathscr {J}$. Further, it is also convenient to introduce the operators 
$\mathscr {R}$ and 
$\mathscr {I}$ defined by
\begin{equation} \mathscr{R}\stackrel{\textrm def}{=} \mathrm{Re} \{\exp[-\mathrm{i} h \partial_{{{z_{s}}}}]\}, \quad \mathscr{I}\stackrel{\textrm def}{=}-\mathrm{Im} \{\exp[-\mathrm{i} h \partial_{{{z_{s}}}}]\}\partial_x^{{-}1}, \end{equation}
so 
$\mathscr {J}=\mathscr {R}^{-1}\mathscr {I}$.
Although explicit, the formulae (2.5) and (2.6) are not quite in closed form since they involve series (via the definition of the exponential operator) and operator inversion. Additional relations, suitable for practical applications, are thus derived below.
3. Auxiliary relations
With different choices of 
$z$ and 
$z_0$, the Taylor expansion (2.1) provides various relations of practical interest. Several variants of (2.5) can then be derived, their convenience depending on the problem at hand.
With the choice 
$z={{z_{b}}}$ and 
$z_0=-\mathrm {i} h=-\mathrm {i}(d+\eta )$, the relation (2.1) becomes
\begin{equation} f({{z_{b}}}+\mathrm{i} h) = f({{z_{s}}}) = \exp [\mathrm{i} h \partial_{{{z_{b}}}}]f({{z_{b}}}), \end{equation}so a comparison with (2.2) yields at once
\begin{equation} \left(\exp[-\mathrm{i} h \partial_{{{z_{s}}}}]\right)^{{-}1} = \exp[\mathrm{i} h\partial_{{{z_{b}}}}] \quad \Longleftrightarrow \quad \left(\exp[\mathrm{i} h\partial_{{{z_{b}}}}]\right)^{{-}1} = \exp[-\mathrm{i} h\partial_{{{z_{s}}}}]. \end{equation}
With this relation, the operator involving 
$\partial _{{{z_{s}}}}$ in the DNO (2.5) can be replaced by one involving 
$\partial _{{{z_{b}}}}$. This is somewhat convenient in constant depth because, then, 
$\partial _{{{z_{b}}}}=\partial _x$. However, two operators then need to be inverted instead of one with (2.5), so further simplifications are desirable.
Taking 
$z=x$ together with 
$z_0=-\mathrm {i}\eta$ and 
$z_0=\mathrm {i} d$, (2.1) yields
\begin{equation} f(x+\mathrm{i}\eta) = f({{z_{s}}}) = \exp\!\left[\mathrm{i}\eta\partial_x\right]f(x), \quad f(x-\mathrm{i} d) = f({{z_{b}}}) = \exp\!\left[-\mathrm{i} d\partial_x\right]f(x), \end{equation}
and the elimination of 
$f(x)$ between these two relations, together with (2.2), yields
\begin{equation} \exp[-\mathrm{i} h\partial_{{{z_{s}}}}]=\exp\!\left[-\mathrm{i} d\partial_x\right] \left(\exp\!\left[\mathrm{i}\eta \partial_x\right]\right)^{{-}1} =\exp\!\left[-\mathrm{i} d\partial_x\right] \left(\exp\!\left[-\mathrm{i}\eta\partial_x\right]\,\right)^{{\dagger}}\left(1+ \mathrm{i}\eta_x\right), \end{equation}
where a 
${\dagger}$ denotes the adjoint operator. (For any complex function 
$\gamma$ of a single real variable 
$x$, the operator 
$\exp [\gamma \partial _x]\stackrel {\textrm def}{=}\sum _{n=0} ^{\infty }(n!)^{-1}\gamma ^{n}\partial _x^{\,n}$ has for Hermitian adjoint 
$(\exp [\gamma \partial _x])^{{\dagger}}\stackrel {\textrm def}{=}\sum _{n=0}^{\infty }(n!)^{-1}\partial _x^{n} (-\gamma ^{\ast })^{n}$, a star denoting the complex conjugate. We then have 
$(\exp [\gamma \partial _x])^{-1}=(\exp [\gamma ^{\ast }\partial _x])^{{\dagger}}(1+\gamma _x)$.) We thus have relations that allow us to avoid the computation of the 
$\partial _z$ operators, moreover without inversions. The operator 
$\mathscr {R}$ remains to be inverted, however.
For a real or complex function 
$\gamma$ depending on a single real variable 
$x$, let the operators and their Hermitian adjoints be
\begin{gather} \mathscr{C}_\gamma \stackrel{\textrm def}{=} \sum_{n=0}^{\infty}\frac{({-}1)^{n}}{(2n)!}\,\gamma^{2n} \partial_x^{2n},\quad \mathscr{S}_\gamma \stackrel{\textrm def}{=} \sum_{n=0}^{\infty}\frac{({-}1)^{n}}{(2n+1)!}\, \gamma^{2n+1}\,\partial_x^{2n+1}, \end{gather}
\begin{gather}\mathscr{C}_\gamma^{{\dagger}}=\sum_{n=0}^{\infty}\frac{({-}1)^{n}}{(2n)!}\,\partial_x^{2n} {\gamma^{{\ast}}}^{2n},\quad \mathscr{S}_\gamma^{{\dagger}}=\sum_{n=0}^{\infty}\frac{({-}1)^{n+1}}{(2n+1)!}\, \partial_x^{2n+1}\,{\gamma^{{\ast}}}^{2n+1}. \end{gather}We then have
\begin{equation} \exp\!\left[-\mathrm{i} d\partial_x\right]=\mathscr{C}_d - \mathrm{i}\,\mathscr{S}_d, \quad \left(\exp\!\left[-\mathrm{i}\eta\,\partial_x\right]\,\right)^{{\dagger}}= \mathscr{C}_\eta^{{\dagger}}+\mathrm{i}\,\mathscr{S}_\eta^{{\dagger}}, \end{equation}and the relation (3.4) is split into real and imaginary parts as
\begin{gather} \mathrm{Re}\{\exp[-\mathrm{i} h\partial_{{{z_{s}}}}]\}= \mathscr{C}_d\,\mathscr{C}_\eta^{{\dagger}}+\mathscr{S}_d\,\mathscr{S}_\eta^{{\dagger}}- \mathscr{C}_d\,\mathscr{S}_\eta^{{\dagger}}\eta_x + \mathscr{S}_d\,\mathscr{C}_\eta^{{\dagger}} \eta_x, \end{gather}
\begin{gather}\mathrm{Im}\{\exp[-\mathrm{i} h\partial_{{{z_{s}}}}]\}= \mathscr{C}_d\,\mathscr{S}_\eta^{{\dagger}}-\mathscr{S}_d\,\mathscr{C}_\eta^{{\dagger}}+ \mathscr{C}_d\,\mathscr{C}_\eta^{{\dagger}}\eta_x + \mathscr{S}_d\,\mathscr{S}_\eta^{{\dagger}} \eta_x. \end{gather}
With the operator relation 
$\eta _x=\partial _x\eta -\eta \partial _x$ (resulting from the Leibniz rule), we have
\begin{equation} \mathscr{C}_\eta^{{\dagger}}\eta_x=\partial_x^{{-}1}\mathscr{S}_\eta^{{\dagger}}\partial_x - \mathscr{S}_\eta^{{\dagger}}, \quad \mathscr{S}_\eta^{{\dagger}}\eta_x= \mathscr{C}_\eta^{{\dagger}} -\partial_x^{{-}1}\,\mathscr{C}_\eta^{{\dagger}}\partial_x, \end{equation}so the relations (3.8) and (3.9) yield
\begin{equation} \mathscr{R}=\mathscr{C}_d\partial_x^{{-}1}\,\mathscr{C}_\eta^{{\dagger}}\partial_x + \mathscr{S}_d\partial_x^{{-}1}\,\mathscr{S}_\eta^{{\dagger}}\partial_x,\quad \mathscr{I}=\mathscr{S}_d\partial_x^{{-}1}\,\mathscr{C}_\eta^{{\dagger}}- \mathscr{C}_d\partial_x^{{-}1}\,\mathscr{S}_\eta^{{\dagger}}. \end{equation}The latter relations are particularly convenient to derive analytic approximations and to extrapolate the DNO to higher dimensions, as shown below.
4. Approximate Dirichlet–Neumann operators
From the explicit DNO (2.5) and the relations derived in the previous section, several approximations of practical interest can be easily obtained. We consider here only two special cases.
4.1. Infinitesimal waves in arbitrary depth
Assuming that the free surface 
$\eta$ remains close to zero, one can formally expand the DNO in increasing order of nonlinearities in 
$\eta$ (Craig & Sulem Reference Craig and Sulem1993; Craig et al. Reference Craig, Guyenne, Nicholls and Sulem2005). Thus, writing 
$\mathscr {G}=\mathscr {G}_0+\mathscr {G}_1+\mathscr {G}_2+\cdots$ and similarly for 
$\mathscr {R}$, 
$\mathscr {I}$ and 
$\mathscr {J}$, one obtains at once from (3.11) the following:
\begin{gather} \mathscr{R}_0 = \mathscr{C}_d, \quad \mathscr{R}_1 ={-}\mathscr{S}_d\eta\partial_x, \quad \mathscr{R}_2 ={-}{{\textstyle{1\over 2}}}\,\mathscr{C}_d\partial_x\eta^{2}\partial_x, \quad \mathrm{etc.}, \end{gather}
\begin{gather}\mathscr{I}_0 = \mathscr{S}_d\partial_x^{{-}1}, \quad \mathscr{I}_1 = \mathscr{C}_d\eta, \quad \mathscr{I}_2 ={-}{{\textstyle{1\over 2}}}\,\mathscr{S}_d\partial_x\eta^{2}, \quad \mathrm{etc}. \end{gather}The relation
\begin{align} \mathscr{R}^{{-}1}&=[\mathscr{R}_0+\mathscr{R}_1+\mathscr{R}_2+\cdots]^{{-}1} =[1+\mathscr{R}_0^{{-}1}\,\mathscr{R}_1+\mathscr{R}_0^{{-}1}\,\mathscr{R}_2+ \cdots]^{{-}1}\mathscr{R}_0^{{-}1}\nonumber\\ &=[1-\mathscr{R}_0^{{-}1}\,\mathscr{R}_1- \mathscr{R}_0^{{-}1}\,\mathscr{R}_2+\mathscr{R}_0^{{-}1}\,\mathscr{R}_1\, \mathscr{R}_0^{{-}1}\,\mathscr{R}_1+\cdots]\mathscr{R}_0^{{-}1} \end{align}then yields, after some algebra,
\begin{equation} \left.\begin{gathered} \mathscr{J}_0=\mathscr{C}_d^{{-}1}\,\mathscr{S}_d\,\partial_x^{{-}1}, \quad \mathscr{J}_1 = \eta + \mathscr{J}_0\partial_x\eta\partial_x\mathscr{J}_0, \\ \mathscr{J}_2={{\textstyle{1\over 2}}}\,\partial_x\eta^{2}\partial_x\mathscr{J}_0 + \mathscr{J}_0\partial_x\eta\partial_x\,\mathscr{J}_1-{{\textstyle{1\over 2}}}\,\mathscr{J}_0 \partial_x^{2}\eta^{2}, \quad \mathrm{etc.}, \end{gathered}\right\} \end{equation}hence
\begin{equation} \left.\begin{gathered} \mathscr{G}_0={-}\partial_x\,\mathscr{C}_d^{{-}1}\,\mathscr{S}_d, \quad \mathscr{G}_1={-}\partial_x\eta\partial_x-\mathscr{G}_0\eta\,\mathscr{G}_0, \\ \mathscr{G}_2={{\textstyle{1\over 2}}}\,\partial_x^{2}\eta^{2}\,\mathscr{G}_0-\mathscr{G}_0\eta\,\mathscr{G}_1 -{{\textstyle{1\over 2}}}\,\mathscr{G}_0\partial_x\eta^{2}\partial_x, \quad\mathrm{etc}. \end{gathered}\right\}\end{equation} In constant depth, the expansion of Craig & Sulem (Reference Craig and Sulem1993) is, as expected, recovered by introducing the operator 
$\mathscr {D}\stackrel {\textrm def}{=}\mathrm {i}\partial _x$, i.e. replacing 
$\partial _x$ by 
$-\mathrm {i}\mathscr {D}$. With a variable bottom, the expansion of Craig et al. (Reference Craig, Guyenne, Nicholls and Sulem2005) is also recovered, except for the definition of 
$\mathscr {G}_0$. Indeed, Craig et al. (Reference Craig, Guyenne, Nicholls and Sulem2005) define 
$\mathscr {G}_0$ with an expansion for small amplitudes of the bottom corrugation (i.e. 
${\max }|d(x)-\bar {d}|$ is small, where 
$\bar {d}$ is the mean depth), and they provide a recursion formula for computing this series. In (4.5), 
$\mathscr {G}_0$ is defined explicitly for arbitrary (non-overturning) bottom and no additional expansions are required.
4.2. Remarks
For higher-order approximations, the recursion formula of Craig & Sulem (Reference Craig and Sulem1993) can be used verbatim with 
$\mathscr {G}_0$ defined here in (4.5). This approach is convenient for the derivation of (rather low-order) analytical approximations. However, with numerical computations, this recursion is prone to cancellation errors leading to large numerical errors and instabilities (Wilkening & Vasan Reference Wilkening and Vasan2015). This problem is more pronounced in higher dimensions (W. Craig, personal communication 2005).
These difficulties come mostly from the expansion of the inverse operator 
$\mathscr {R}^{-1}$. For numerical computations, this expansion should be avoided to obtain 
$B=\mathscr {R }^{-1}A$ (for some functions 
$A$ and 
$B$). It is generally more efficient to solve 
$\mathscr {R}B=A$ via an iterative procedure. This is a similar problem to the resolution of linear systems of equations, for which iterative methods are often more efficient (Isaacson & Keller Reference Isaacson and Keller1994). For the DNO, the relation (3.11a) (see also relation (5.11)) shows that 
$\mathscr {R}$ behaves (roughly) like a 
$\cosh$ function, so 
$\mathscr {R}^{-1}$ behaves like a sech function. The Maclaurin series of 
$\cosh (z)$ having an infinite radius of convergence, while that of 
$\operatorname {sech}(z)$ converges only for 
$|z|<{\rm \pi} /2$, this provides an informal/heuristic argument showing why 
$B=\mathscr {R}^{-1}A$ should not be computed but 
$\mathscr {R}B=A$ should be solved instead. With other representations (than truncated Taylor series) of 
$\mathscr {R}^{-1}$, the computation of 
$B=\mathscr {R}^{-1}A$ may be efficient, however.
For linear waves in the context of a highly variable bathymetry, the improvements of the DNO expansion proposed by Andrade & Nachbin (Reference Andrade and Nachbin2018) could be exploited to reformulate the explicit DNO in a more effective form for numerical computations. However, when speed and high numerical accuracy are required, the DNO perturbation expansions are not competitive (specially for steep waves) and boundary integral formulations should be preferred (Clamond & Grue Reference Clamond and Grue2001; Fructus et al. Reference Fructus, Clamond, Grue and Kristiansen2005; Fructus & Grue Reference Fructus and Grue2007).
4.3. Long waves in shallow water
For long waves in shallow water, the characteristic wavelength 
$L_c$ is much larger than the characteristic depth 
$d_c$, so 
$\sigma \stackrel {\textrm def}{=} d_c/L_c\ll 1$ is a ‘shallowness’ dimensionless small parameter. The horizontal derivative 
$\partial _x$ is then of first order in shallowness and the DNO can be expanded is power series of 
$\sigma$, without assuming small amplitude for the waves and/or for the bottom corrugation. Thus, we do not need to explicitly introduce scalings to assess the order of terms – it is sufficient to count the number of derivatives. For instance, 
$\partial _x^{3} \eta$, 
$(\partial _x^{2}\eta )(\partial _x\eta )$ and 
$(\partial _x\eta )^{3}$ are all of third order in shallowness, as well as 
$\partial _x^{3}d$, 
$(\partial _x^{2} d)(\partial _xd)$ and 
$(\partial _xd)^{3}$.
We then have the shallow-water even-terms expansions 
$\mathscr {R}=\mathscr {R}_0+\mathscr {R}_2 +\mathscr {R}_4+\cdots$ (and similarly for 
$\mathscr {I}$ and 
$\mathscr {J}$), so, from (3.11),
\begin{gather} \mathscr{R}_0=1,\quad \mathscr{R}_2={-}{{\textstyle{1\over 2}}}\,d^{2}\partial_x^{2}- d\partial_x\eta\partial_x-{{\textstyle{1\over 2}}}\,\partial_x\eta^{2}\partial_x,\quad \mathrm{etc.}, \end{gather}
\begin{gather}\mathscr{I}_0=h,\quad \mathscr{I}_2={-}{{\textstyle{1\over 6}}}\,d^{3}\partial_x^{2}- {{\textstyle{1\over 2}}}\,d^{2}\partial_x^{2}\eta-{{\textstyle{1\over 2}}}\,d\partial_x^{2}\eta^{2}- {{\textstyle{1\over 6}}}\,\partial_x^{2}\eta^{3},\quad \mathrm{etc.}, \end{gather}hence, after some algebra,
\begin{equation} \mathscr{J}_0=h,\quad \mathscr{J}_2={{\textstyle{1\over 2}}}\,h^{2}d_{xx}+ hh_xd_x-hd_x^{2}+{{\textstyle{1\over 3}}}\,\partial_xh^{3}\partial_x, \quad \mathrm{etc.} \end{equation}
Note that 
$\mathscr {J}_0$ and 
$\mathscr {J}_2$ are obviously self-adjoint, as they should be.
It should be emphasised that these approximations were obtained directly from the explicit DNO, considering weak variations in 
$x$ (i.e. long waves in shallow water) but without assuming small amplitudes of the free surface and of the seabed (i.e. there are no restrictions on the magnitude of 
$|\eta |$ and 
$|d(x)-\bar {d}|$, 
$\bar {d}$ being the mean depth).
5. Dirichlet–Neumann operator in higher dimensions
It is rather straightforward to extrapolate the DNO given by (2.5) to three (and more) spatial dimensions. In higher dimensions, the holomorphic functions cannot be used but series representations remain. This feature is exploited here to obtain an explicit expression for the DNO in an arbitrary number of dimensions.
With 
$\boldsymbol {x}=(x_1,x_2,\ldots,x_N)\in \mathbb {R}^{N}$ referring to the ‘horizontal’ coordinates, the mathematical problem is then posed in the 
$(N+1)$-dimensional Cartesian 
$(\boldsymbol {x},y)$-space, with 
$y$ the ‘upward vertical’ coordinate. Obviously, only the two-dimensional (i.e. 
$N=1$) and three-dimensional (i.e. 
$N=2$) cases are of physical interest for water waves. Let 
$\boldsymbol {\nabla }\stackrel {\textrm def}{=}(\partial _{x_1},\ldots,\partial _{x_N})$, 
$\varDelta \stackrel {\textrm def}{=} \boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\nabla }$ and 
$\mathscr {D}\stackrel {\textrm def}{=}(-\varDelta )^{1/2}$ denote, respectively, the horizontal gradient, Laplacian and semi-Laplacian operators.
The DNO is naturally extended in higher dimensions by extrapolating the relation 
$\mathscr {G}=-\partial _x\mathscr {R}^{-1}\mathscr {I}\partial _x$, the operators 
$\mathscr {R}$ and 
$\mathscr {I}$ having to be redefined. In the two-dimensional case, these operators are defined via complex expressions in § 2. In order to extend these operators to higher dimensions, one must consider their real form (3.11), so their extrapolation is natural.
One-dimensional operators involving only even-order derivatives have straightforward extensions in higher dimensions by replacing the second-order horizontal derivative 
$\partial _x^{2}$ by the horizontal Laplacian 
$\varDelta$. For instance,
\begin{gather} \mathscr{C}_d \mapsto \sum_{n=0}^{\infty}\frac{({-}1)^{n}}{(2n)!}\,d^{2n} \varDelta^{n} \stackrel{\textrm def}{=} \cosh (d\,\mathscr{D}), \end{gather}
\begin{gather}\mathscr{C}_\eta^{{\dagger}} \mapsto \sum_{n=0}^{\infty}\frac{({-}1)^{n}}{(2n)!}\, \varDelta^{n}\eta^{2n}\stackrel{\textrm def}{=} \cosh(\mathscr{D}\,\eta), \end{gather}
\begin{gather}\mathscr{S}_d\partial_x^{{-}1}\mapsto \sum_{n=0}^{\infty}\frac{({-}1)^{n}}{(2n+1)!} \,d^{2n+1}\varDelta^{n} \stackrel{\textrm def}{=} \sinh(d\,\mathscr{D})\mathscr{D}^{{-}1}, \end{gather}
\begin{gather}\partial_x^{{-}1}\mathscr{S}_\eta^{{\dagger}} \mapsto \sum_{n=0}^{\infty}\frac{({-}1)^{n+1}}{(2n+1)!}\, \varDelta^{n}\eta^{2n+1} \stackrel{\textrm def}{=} -\mathscr{D}^{{-}1}\sinh(\mathscr{D}\,\eta). \end{gather}
It should be emphasised that, as in the one-dimensional case, the operators do not commute, so, for example, 
$\cosh (d\mathscr {D})\neq \cosh ( \mathscr {D}d)$ and 
$\cosh (d\mathscr {D})^{-1}\neq \operatorname {sech}(d\mathscr {D})$, the equalities holding only in constant depth because then 
$d\mathscr {D}=\mathscr {D}d$.
A natural extension of 
$\mathscr {I}\partial _x$ is thus 
$\mathscr {I}\boldsymbol {\nabla }$ with
\begin{equation} \mathscr{I} \mapsto \sinh(d\,\mathscr{D})\mathscr{D}^{{-}1} \cosh(\mathscr{D}\,\eta)+\cosh(d\,\mathscr{D}) \mathscr{D}^{{-}1}\sinh(\mathscr{D}\,\eta). \end{equation}
In order to find the extension of 
$\partial _x\mathscr {R}^{-1}$, the operator 
$\mathscr {R}$ given by (3.11a) is rewritten as
\begin{equation} \mathscr{R}=\mathscr{C}_d\partial_x{}^{{-}1}[\mathscr{C}_\eta^{{\dagger}} +\partial_x\,\mathscr{C}_d^{{-}1}(\mathscr{S}_d\,\partial_x^{{-}1}) \partial_x (\partial_x^{{-}1}\,\mathscr{S}_\eta^{{\dagger}})]\partial_x. \end{equation}Thus, we have the natural extension
\begin{equation} \partial_x\,\mathscr{R}^{{-}1}\mapsto[\cosh(\mathscr{D}\eta) +\mathscr{G}_0\,\mathscr{D}^{{-}1}\sinh(\mathscr{D}\eta)]^{{-}1} \boldsymbol{\nabla}\boldsymbol{\cdot}\cosh(d\mathscr{D})^{{-}1}, \end{equation}
where 
$\cosh (d\mathscr {D})^{-1}$ is the inverse operator of 
$\cosh (d\mathscr {D})$, and
\begin{equation} \mathscr{G}_0 \stackrel{\textrm def}{=} -\boldsymbol{\nabla}\boldsymbol{\cdot}\cosh(d\mathscr{D})^{{-}1}\sinh (d\mathscr{D})\mathscr{D}^{{-}1}\,\boldsymbol{\nabla}.\end{equation}Therefore, the DNO becomes at once
\begin{equation} \mathscr{G}={-}[\cosh(\mathscr{D}\eta)+\mathscr{G}_0 \,\mathscr{D}^{{-}1}\sinh(\mathscr{D}\eta)]^{{-}1} \boldsymbol{\nabla}\boldsymbol{\cdot}\cosh(d\mathscr{D})^{{-}1}\mathscr{I}\,\boldsymbol{\nabla}. \end{equation}
In order to avoid misinterpretations of the formula (5.9), it is worth re-emphasising here that: (i) any operator acts on everything it multiplies on its right, so (5.9) should be applied successively leftward starting from the furthest right; and (ii) exponents denote operator compositions, so an exponent 
$-1$ means an operator inversion.
Note that 
$\mathscr {G}\to \mathscr {G}_0$ as 
$\eta \to 0$. Moreover, processing as in § 4.1 for infinitesimal waves, one finds the expansion of Craig et al. (Reference Craig, Guyenne, Nicholls and Sulem2005), except for 
$\mathscr {G}_0$ that is defined implicitly by Craig et al. (Reference Craig, Guyenne, Nicholls and Sulem2005) but explicitly here.
5.1. Constant depth
In constant depth, 
$d$ commuting then with both 
$\mathscr {D}$ and 
$\boldsymbol {\nabla }$, we have the simplified relation 
$\mathscr {G}_0=\mathscr {D}\tanh (d \mathscr {D})$, while 
$\mathscr {I}$ and 
$\mathscr {G}$ become (see Appendix A for details)
\begin{align} \mathscr{I}&=\mathscr{D}{}^{{-}1}[\sinh(d\mathscr{D}) \cosh(\mathscr{D}\eta)+\cosh(d\mathscr{D})\sinh(\mathscr{D}\eta)] =\mathscr{D}^{{-}1}\sinh(\mathscr{D}h),\end{align}
\begin{align} \mathscr{G}&={-}[\cosh(d\mathscr{D})\cosh(\mathscr{D}\eta) +\sinh(d\mathscr{D})\sinh(\mathscr{D}\eta)]^{{-}1} \boldsymbol{\nabla}\boldsymbol{\cdot}\mathscr{I}\,\boldsymbol{\nabla}\nonumber\\ &={-}\cosh(\mathscr{D}h)^{{-}1}\,\mathscr{D}^{{-}1}\,\boldsymbol{\nabla}\boldsymbol{\cdot} \sinh(\mathscr{D}h)\boldsymbol{\nabla}. \end{align}
A better-conditioned formulation, avoiding the computation of 
$\mathscr {D}$, is
\begin{equation} \mathscr{G}={-}[\operatorname{sech}(d\mathscr{D}) \cosh(\mathscr{D}h)]^{{-}1}\,\boldsymbol{\nabla}\boldsymbol{\cdot}[ \operatorname{sech}(d\mathscr{D})\operatorname{sinhc} (\mathscr{D}h)h]\boldsymbol{\nabla}, \end{equation}with
\begin{equation} \operatorname{sech}(d\mathscr{D})\stackrel{\textrm def}{=}\sum_{n=0}^{\infty} \frac{E_{2n}}{(2n)!}\,d^{2n}\,\mathscr{D}^{2n}, \quad \operatorname{sinhc}(\mathscr{D}h)\stackrel{\textrm def}{=}\sum_{n=0}^{\infty} \frac{1}{(2n+1)!}\,\mathscr{D}^{2n}h^{2n}, \end{equation}
where 
$E_n$ are the Euler numbers (Abramowitz & Stegun Reference Abramowitz and Stegun1965) (since 
$d$ and 
$\mathscr {D}$ commute, we have 
$\operatorname {sech}(d\mathscr {D})=\cosh (d \mathscr {D})^{-1}$). The relation (5.12) involving only even powers of 
$\mathscr {D}$, only Laplacian and gradient operators need to be evaluated, i.e. the computation of the non-local operator 
$\mathscr {D}$ can be avoided.
As the DNO appears in Hamiltonian formulations of water waves, its functional variations are crucial to derive the equations of motion and to investigate stability (Fazioli & Nicholls Reference Fazioli and Nicholls2010). Thanks to the explicit DNO (5.11), these variations can be obtained quite effortlessly. Indeed, with the relations (cf. Appendix A)
\begin{gather} \cosh(\mathscr{D}(h+\delta h))= \cosh(\mathscr{D}h)+ \mathscr{D}\sinh(\mathscr{D} h)\delta h + O(\delta h^{2}), \end{gather}
\begin{gather}\sinh(\mathscr{D}(h+\delta h))= \sinh(\mathscr{D}h)+\mathscr{D}\cosh(\mathscr{D} h)\delta h+O(\delta h^{2}), \end{gather}the first variation of the DNO is obtained at once as
\begin{align} \mathscr{G}(h+\delta h)&=\mathscr{G}(h)-\cosh(\mathscr{D}h )^{{-}1}\,\boldsymbol{\nabla}\boldsymbol{\cdot}\cosh(\mathscr{D}h)\delta h\,\boldsymbol{\nabla}\nonumber\\ &\quad+\cosh(\mathscr{D}h)^{{-}1}\mathscr{D}\sinh(\mathscr{D} h)\delta h\,\mathscr{G}(h)+O(\delta h^{2}), \end{align}or
\begin{align} \mathscr{G}(h+\delta h)&=\mathscr{G}(h) -\cosh(\mathscr{D}h )^{{-}1}\,\boldsymbol{\nabla}\boldsymbol{\cdot}\cosh(\mathscr{D} h)\delta h\boldsymbol{\nabla} + \mathscr{G}(h)\delta h\,\mathscr{G}(h)\nonumber\\ &\quad- \cosh(\mathscr{D}h)^{{-}1}\boldsymbol{\nabla}\boldsymbol{\cdot} \cosh(\mathscr{D}h)(\boldsymbol{\nabla} h)\delta h\,\mathscr{G}(h) + O(\delta h^{2}). \end{align}
Similarly, higher-order functional variations of 
$\mathscr {G}$ can be easily obtained. This is one illustration of the advantage of dealing with an explicit DNO.
5.2. Remarks
Since the multi-dimensional DNO was derived by extrapolating from the two-dimensional case, one can then naturally ask if (5.9) is a correct expression.
First, we note that the explicit expression of the DNO is not unique. For instance, as in two dimensions and as suggested by the Taylor expansion around 
$\eta =0$, the DNO could also be written as 
$\mathscr {G}=-\boldsymbol {\nabla }\boldsymbol {\cdot }\mathscr {J}\boldsymbol {\nabla }$ for some operator 
$\mathscr {J}(d,\eta, \mathscr {D},\boldsymbol {\nabla })$ to be specified.
In two dimensions (i.e. for 
$N=1$), one can exploit the theory of holomorphic functions to directly check that the explicit DNO (5.9) is a correct one. This procedure is simply the reverse of the derivations made in §§ 2 and 3. This is not possible in higher dimensions (i.e. 
$N>1$) because holomorphic functions cannot be used. The validity of (5.9) was then checked by expanding it à la Craig & Sulem, checking that the two expansions match. (This is detailed in Appendix A for constant depth.)
6. Moving bottom
We consider finally the generalisation of a moving bottom, i.e. 
$d=d(x,t)$. Of course, for simplicity, we begin with the two-dimensional case, the generalisation to higher dimensions being straightforward.
When 
$\partial _td\neq 0$, the bottom is no longer a streamline, so the streamfunction is not zero at the seabed, i.e. 
${\psi _{b}}={\psi _{b}}(x,t)\neq 0$. The lower boundary condition (1.2) becomes 
$\partial _td=\partial _x{{\psi _{b}}}=-{{v_{b}}}- {{u_{b}}}\partial _xd$. With a moving bottom, the relations (2.1), (2.2) and (3.1) still hold, but (2.3) becomes
\begin{equation} {{\psi_{b}}} = \mathrm{Re}\{\exp[-\mathrm{i} h\partial_{{{z_{s}}}}]\}{{\psi_{s}}} + \mathrm{Im}\{\exp[-\mathrm{i} h\partial_{{{z_{s}}}}]\}{{\phi_{s}}}. \end{equation}
The condition for the bottom impermeability yielding 
${{\psi _{b}}}=\partial _x^{-1} \partial _td$, the relation (6.1) gives
\begin{equation} {{\psi_{s}}}=\mathrm{Re}\{\exp[-\mathrm{i} h\partial_{{{z_{s}}}}] \}{}^{{-}1}(\partial_x^{{-}1}\partial_td-\mathrm{Im}\{\exp[-\mathrm{i} h\partial_{{{z_{s}}}}]\}{{\phi_{s}}}). \end{equation}
The relation (6.2) shows that the Dirichlet-to-Neumann transformation at the free surface is no longer a homogeneous linear function of 
${{\phi _{s}}}$. The impermeability of the free surface is then 
$\partial _t\eta =G({{\phi _{s}}})$, the generalised DNO 
$G$ being
\begin{equation} G\!\left({{\phi_{s}}}\right) = \mathscr{G}\,{{\phi_{s}}} - \partial_x\, \mathrm{Re}\{\exp[-\mathrm{i} h\,\partial_{{{z_{s}}}}]\}^{{-}1}\partial_x^{{-}1}\,\partial_td, \end{equation}
where 
$\mathscr {G}$ is given by (2.5). Note that 
$\partial _x^{-1}\partial _td$ is not uniquely defined due to the antiderivative, uniqueness being enforced by the definitions of the mean water level and of the frame of reference.
In higher dimensions, with 
$\partial _x^{-1}=\partial _x\partial _x^{-2}\mapsto \boldsymbol {\nabla }\varDelta ^{-1}=-\mathscr {D}^{-2}\boldsymbol {\nabla }$, the DNO obviously becomes
\begin{equation} G\!\left({{\phi_{s}}}\right) =[\cosh(\mathscr{D}\,\eta)+\mathscr{G}_0 \,\mathscr{D}^{{-}1}\sinh(\mathscr{D}\,\eta)]^{{-}1} \boldsymbol{\nabla}\boldsymbol{\cdot}\cosh(d\,\mathscr{D})^{{-}1}(\mathscr{D}^{{-}2}\boldsymbol{\nabla} \partial_td-\mathscr{I}\boldsymbol{\nabla}{{\phi_{s}}}), \end{equation}
with 
$\mathscr {I}$ and 
$\mathscr {G}_0$ being defined, respectively, by (5.5) and (5.8).
7. Discussion
Using elementary algebra, we obtained explicit formulae for the DNOs involved in water wave problems. We first derived the DNO for two-dimensional waves over a static (uneven) bottom. We then extrapolated the formula to higher dimensions and generalised the formula for moving bottoms. The latter generalisation is interesting for its applications, such as tsunami generation (Iguchi Reference Iguchi2011), but also because it shows that extensions to fluids stratified into several homogeneous layers is possible (Craig, Guyenne & Kalisch Reference Craig, Guyenne and Kalisch2005; Constantin & Ivanov Reference Constantin and Ivanov2019). The DNO is also used in some water wave problems with vorticity (Constantin, Ivanov & Martin Reference Constantin, Ivanov and Martin2016; Groves & Horn Reference Groves and Horn2020), and the derivation of an explicit DNO for rotational waves is conceivable.
In this short paper, the focus is on the DNO at the free surface assuming a given bottom shape and motion. Obviously, one can as easily obtain the DNO at the bottom from an assumed free surface, which, in particular, should find applications in bottom detection from free-surface measurements (Fontelos et al. Reference Fontelos, Lecaros, López-Ríos and Ortega2017).
The explicit DNOs derived here are expressed with pseudo-differential operators formally defined in terms of series. Such definition supposes sufficient regularity of the free surface and the bottom, regularity yet to be specified by rigorous mathematical analysis. When these regularity conditions are not met, other more general representations of the operators should be used instead, such as integral formulations. Once these operators are properly defined, the explicit DNO should then be usable verbatim, allowing the investigation of rough bottoms and waves with angular crests, for example.
The main purpose of this paper is to show how explicit DNOs can be derived and, via examples, to show their interest for analytic manipulations. Although some indications on potential issues and remedies with numerical computations are briefly discussed, it is not the purpose here to derive the most effective way to compute the DNO numerically. For special functions, their definitions via power series are often not suitable for accurate fast computations, at least not in every case and without extra knowledge (e.g. periodicity, symmetries, locations of singularities). The situation is similar with DNOs defined via series, with the substantial extra difficulty that they involve non-commutative algebra.
DNOs appear in many fields of research in physics (acoustics, elasticity, electromagnetism, etc.) and, more generally, in the theory of partial differential equations. The use of a DNO is not restricted to problems involving the Laplace equation; it is also commonly employed in close relatives, such as the Helmholtz equation. The elementary formal approach presented here could thus be adapted in these contexts.
Declaration of interests
The author reports no conflict of interest.
Appendix A. Some operator relations in constant depth
In constant depth, the algebra is significantly simplified because 
$d$ commutes with both 
$\mathscr {D}$ and 
$\boldsymbol {\nabla }$. As mentioned at the end of § 5, we then have 
$\mathscr {G}_0=\mathscr {D}\tanh (d\mathscr {D})$. We also have, from the definition of the operators,
\begin{align} \cosh(d\mathscr{D})\cosh(\mathscr{D}\eta)&= \left(\sum_{i=0}^{\infty}\,\frac{d^{2i}\,\mathscr{D}^{2i}}{(2i)!}\right) \left(\sum_{j=0}^{\infty}\,\frac{\mathscr{D}^{2j}\eta^{2j}}{(2j)!}\right) =\sum_{i,j=0}^{\infty}\,\frac{d^{2i}\,\mathscr{D}^{2i+2j}\eta^{2j}}{(2i)!(2j)!}\nonumber\\ &=\sum_{i=0}^{\infty}\,\sum_{j=0}^{i}\,\frac{\mathscr{D}^{2i}\,d^{2i-2j}\,\eta^{2j}}{(2j)! (2i-2j)!}=\sum_{i=0}^{\infty}\,\sum_{j\ \text{even}}^{2i} \,\frac{\mathscr{D}^{2i}\,d^{2i-j}\eta^{j}}{j!(2i-j)!}, \end{align}
\begin{align} \sinh(d\mathscr{D})\sinh(\mathscr{D}\eta)&= \sum_{i,j=0}^{\infty}\,\frac{d^{2i+1}\,\mathscr{D}^{2i+2j+2}\eta^{2j+1}} {(2i+1)!(2j+1)!}=\sum_{i=1}^{\infty}\,\sum_{j=0}^{i-1}\,\frac{\mathscr{D}^{2i} d^{2i-2j-1}\eta^{2j+1}}{(2j+1)!(2i-2j-1)!}\nonumber\\ &=\sum_{i=1}^{\infty}\,\sum_{j\,\text{odd}}^{2i-1} \frac{\mathscr{D}^{2i}d^{2i-j}\eta^{j}}{j!(2i-j)!}. \end{align}
Since 
$|n!|=\infty$ for all negative integers 
$n$, the summation 
$\sum _{i=1}^{\infty }$ in (A2) can be replaced by 
$\sum _{i=0}^{\infty }$. Thus, we have
\begin{align} \cosh(d\mathscr{D})\cosh(\mathscr{D}\eta)+ \sinh(d\mathscr{D})\sinh(\mathscr{D}\eta)&= \sum_{i=0}^{\infty}\frac{\mathscr{D}^{2i}}{(2i)!}\sum_{j=0}^{2i}\, \frac{(2i)!\,d^{2i-j}\eta^{j}}{j!(2i-j)!}\nonumber\\ &= \sum_{i=0}^{\infty}\,\frac{\mathscr{D}^{2i}\,(d+\eta)^{2i}}{(2i)!}= \cosh(\mathscr{D}h). \end{align}Similarly, one can easily derive the relations
\begin{gather} \sinh(d\mathscr{D})\cosh(\mathscr{D}\eta)+ \cosh(d\mathscr{D})\sinh(\mathscr{D}\eta)= \sinh(\mathscr{D}h), \end{gather}
\begin{gather}\cosh(\eta\mathscr{D})\sinh(\mathscr{D}d)+ \sinh(\eta\mathscr{D})\cosh(\mathscr{D}d)= \sinh(h\mathscr{D}), \end{gather}
\begin{gather}\cosh(\eta\mathscr{D})\cosh(\mathscr{D}d)+ \sinh(\eta\mathscr{D})\sinh(\mathscr{D}d)= \cosh(h\mathscr{D}), \end{gather}and, obviously,
\begin{equation} \cosh(\mathscr{D}h)\pm\sinh(\mathscr{D}h)=\exp(\pm\mathscr{D}h). \end{equation}
Note that, as 
$\mathscr {D}$ commutes with 
$d$, but not with 
$\eta$ and 
$h$, these relations are not valid for uneven bottoms and, in constant depth, 
$\sinh (\mathscr {D}h )\neq \sinh (h\mathscr {D})$ for example. However, for varying bottoms, similar relations can be easily obtained if 
$\eta$ is constant.
Taylor expansions around 
$\eta =0$ yield
\begin{gather} \sinh(\mathscr{D}h)=\sinh(d\mathscr{D}) [1+{{\textstyle{1\over 2}}}\mathscr{D}^{2}\eta^{2}+\cdots] +\cosh(d\mathscr{D})[\mathscr{D}\eta+{{\textstyle{1\over 6}}}\mathscr{D}^{3}\eta^{3}+\cdots], \end{gather}
\begin{gather}\cosh(\mathscr{D}h)=\cosh(d\mathscr{D}) [1+{{\textstyle{1\over 2}}}\mathscr{D}^{2}\eta^{2}+\cdots] +\sinh(d\mathscr{D})[\mathscr{D}\eta+{{\textstyle{1\over 6}}}\mathscr{D}^{3}\eta^{3}+\cdots], \end{gather}
hence, with 
$\mathscr {G}_0\stackrel {\textrm def}{=}\mathscr {D}\tanh (d\mathscr {D})$,
\begin{align} (\cosh(\mathscr{D}h))^{{-}1} &=[1+\mathscr{G}_0 \eta+{{\textstyle{1\over 2}}}\,\mathscr{D}^{2}\eta^{2}+{{\textstyle{1\over 6}}}\,\mathscr{G}_0\,\mathscr{D}^{2}\eta^{3}+ \cdots]^{{-}1}\operatorname{sech}(d\mathscr{D})\nonumber\\ &=[1-\mathscr{G}_0\eta-{{\textstyle{1\over 2}}}\,\mathscr{D}^{2}\eta^{2}+\mathscr{G}_0 \eta\,\mathscr{G}_0\eta+\cdots]\operatorname{sech}(d\mathscr{D}), \end{align}
\begin{align} \boldsymbol{\nabla}\boldsymbol{\cdot}\mathscr{D}^{{-}1}\sinh\left(\mathscr{D}h\right)\boldsymbol{\nabla}&={-} \cosh(d\mathscr{D})[\mathscr{G}_0-\boldsymbol{\nabla}\boldsymbol{\cdot}\eta\boldsymbol{\nabla}- {{\textstyle{1\over 2}}}\,\mathscr{G}_0\boldsymbol{\nabla}\boldsymbol{\cdot}\eta^{2}\boldsymbol{\nabla}+\cdots]. \end{align}
Thus, with 
$\mathscr {G}$ defined in (5.12), one gets
\begin{equation} \mathscr{G}=[1+\mathscr{G}_0\eta+{{\textstyle{1\over 2}}}\,\mathscr{D}^{2}\eta^{2}+\cdots]{}^{{-}1} [\mathscr{G}_0-\boldsymbol{\nabla}\boldsymbol{\cdot}\eta\boldsymbol{\nabla} -{{\textstyle{1\over 2}}}\,\mathscr{G}_0\boldsymbol{\nabla}\boldsymbol{\cdot}\eta^{2}\boldsymbol{\nabla}+\cdots]. \end{equation}
Expanding the DNO as 
$\mathscr {G}=\mathscr {G}_0+\mathscr {G}_1+\mathscr {G}_2+\cdots$, with 
$\mathscr {G}$ defined in (5.11), one obtains
\begin{equation} \mathscr{G}_1 ={-}\mathscr{G}_0\eta\,\mathscr{G}_0 - \boldsymbol{\nabla}\boldsymbol{\cdot}\eta\boldsymbol{\nabla}, \quad \mathscr{G}_2={-}{{\textstyle{1\over 2}}}\,\mathscr{D}^{2}\eta^{2}\,\mathscr{G}_0- \mathscr{G}_0\eta\,\mathscr{G}_1- {{\textstyle{1\over 2}}}\,\mathscr{G}_0\boldsymbol{\nabla}\boldsymbol{\cdot}\eta^{2}\boldsymbol{\nabla}, \quad\text{etc.}, \end{equation}so the expansion of Craig & Sulem (Reference Craig and Sulem1993) is recovered.
Substituting 
$h+\delta h$ for 
$h$, for some small 
$\delta h$, we have the first-order Taylor expansions
\begin{align} \cosh(\mathscr{D}(h+\delta h))&\stackrel{\textrm def}{=} \sum_{n=0}^{\infty}\frac{\mathscr{D}^{2n}\,(h+\delta h)^{2n}}{(2n)!}\nonumber\\ &=\sum_{n=0}^{\infty}\frac{\mathscr{D}^{2n}h^{2n}}{(2n)!} + \sum_{n=1}^{\infty} \frac{\mathscr{D}^{2n}h^{2n-1}\delta h}{(2n-1)!}+O(\delta h^{2}) \nonumber\\ &=\cosh(\mathscr{D}h)+\mathscr{D}\sinh(\mathscr{D}h)\delta h + O(\delta h^{2}) \end{align}and
\begin{align} \sinh(\mathscr{D}(h+\delta h))&= \sum_{n=0}^{\infty}\frac{\mathscr{D}^{2n+1}h^{2n+1}}{(2n+1)!}+\sum_{n=0}^{\infty} \frac{\mathscr{D}^{2n+1}h^{2n}\delta h}{(2n)!}+O(\delta h^{2})\nonumber\\ &=\sinh(\mathscr{D}h)+\mathscr{D}\cosh(\mathscr{D} h)\delta h + O(\delta h^{2}), \end{align}hence
\begin{equation} \cosh(\mathscr{D}(h+\delta h))^{{-}1}= [1-\cosh(\mathscr{D}h)^{{-}1}\mathscr{D}\sinh(\mathscr{D} h)\delta h] \cosh(\mathscr{D}h)^{{-}1} + O(\delta h^{2}). \end{equation}We are then in position to compute explicitly the functional variations of the DNO.
