## Abstract

The scattering of a wave from a periodic rough surface in two dimensions is considered. The proposed method is based on the use of a conformal mapping and of the multimodal admittance method. The use of the conformal mappings induces a spatially varying refractive index and transforms the boundary condition on the rough surface into a boundary condition on a flat surface. Then, the multimodal admittance method reduces this problem to a Riccati equation for the modal admittance matrix, which is solved numerically with a Magnus–Möbius scheme. The method is shown to converge exponentially with the number of Fourier modes. The method also allows to find geometries having trapped modes, or quasi-trapped modes, at a given frequency, by looking at singularities, or quasi-singularities, of this Riccati equation. Besides, a simple perturbation expansion, based on a small roughness approximation, is developed.

## 1. Introduction

The scattering of a wave by a rough surface is a subject which has been intensively studied in acoustics but also in optics and electromagnetism. As a result, a very large literature exists on the subject, including several books which deal with this topic (for instance [1–3]). Introduced by Nevière & Cadilhac [4], and then used by other authors (for instance [5–7]), one approach consists in using a conformal mapping. Despite being an old mathematical tool, conformal mappings still remain the object of recent and active research works [8–10] and are now used in different domains of physics. In our case, a conformal mapping is used to first simplify the boundary value problem of the Helmholtz equation: the rough surface is transformed into a flat surface, whereas the bulk Helmholtz equation is not altered drastically. Indeed, the latter equation becomes another Helmholtz equation but with a spatially varying refractive index, corresponding to an equivalent problem. The challenging difficulty is then moved from applying the boundary condition on a complex geometry (the rough surface), to solving this new equation. The family of conformal mappings used for rough surface has the particularity to tend to the identity far from the boundary. As a consequence, the effect of the mapping is localized near the boundary, so that, far from the boundary, we recover a free space Helmholtz equation (with no variation of the index). Thus, the scattered field can be computed in a band adjacent to the boundary, with a simple boundary condition far from the boundary.

In this paper, we propose a multimodal admittance method based on the use of a conformal mapping to compute the scattering of a wave by a periodic rough boundary in the harmonic regime. The Helmholtz problem with varying index, obtained by applying the conformal mapping, is treated as a waveguide problem in the *y*-direction, orthogonal to the boundary. This problem is solved with the multimodal admittance method which was first developed to treat the wave propagation problem in a waveguide with varying cross section [11]. It relies on the projection of the governing equations on the basis of the local transverse Fourier modes. Thanks to this projection, a set of coupled mode equations is obtained and the computed field can be matched to the radiation condition far from the boundary as in preceding works [4,12]. To avoid exponential divergence of the error, due to the presence of evanescent modes, an admittance matrix is introduced [13]. Its variation in the *y*-direction is governed by a Riccati equation of the form Y′(*y*)=*f*(Y(*y*),*y*). This equation is known to be subject to movable singularities in the complex *y* plane [13,14]. These singularities are responsible for very high peaks of coefficients of the admittance matrix (quasi-singularities) on the real *y*-axis. Despite these quasi-singularities, the Riccati equation can be integrated very efficiently using a Magnus–Möbius scheme [14–16]. These singularities present great interest in order to find quasi-trapped modes [17,18]. Indeed, the existence of a trapped mode is linked to a singularity of the admittance/impedance matrix [13] (the admittance matrix cannot be defined at the point where such a trapped mode exists). In this paper, we adapt this method in order to compute the field due to the scattering from a rough surface of an incident plane wave and to find quasi-trapped modes, or quasi-edge waves, of such a rough surface. The particularity of this method is that it allows us to find, among a given family of geometries, the geometries having quasi-trapped modes at a given frequency.

In §2, the governing equations are presented, and the conformal mapping is applied to them. The multimodal method is detailed in §3. Section 4 presents the different conformal mappings that will be used in the paper. Particular attention is given to a generic mapping, able to describe an arbitrary geometry. Applications of our method are given in §5, and some numerical considerations are treated in §6. The method to find trapped modes is presented in §7. Incidentally, a simple perturbative approach to solve the problem at the first order is developed in §8.

## 2. The governing equations

We aim at solving the problem of the scattering of an incident plane wave by a rigid periodic rough boundary with Neumann boundary condition. Note that the case of Dirichlet boundary condition might be easily treated also. We consider the harmonic regime with implicit time dependence *e*^{−iωt}, that will be omitted in the following. The space coordinates (*X*,*Y*) are supposed to be non-dimensionalized by *D*/2*π* and the time by *D*/2*πc*, with *D* being the spatial period of the boundary and *c* the sound velocity supposed constant. Incidentally, the boundary is 2*π*—periodic in the *X*-direction. The boundary is represented by a curve *X*,*Y*) plane, as sketched in figure 1*a*. The domain considered is the upper infinite half plane above this curve *Φ* satisfies the Helmholtz equation with rigid boundary condition on *a*
and
*b*
where *k*=*ωD*/2*πc* is the dimensionless wavenumber and where Δ_{XY} and **∇**_{XY} are the Laplacian and the gradient operators with respect to the (*X*,*Y*) variables, and ** n** is the unit normal vector to

The most challenging issue to solve equations (2.1) comes from applying the boundary condition on *b*. More precisely, the conformal mapping is an analytic function of the complex variable *z*=*x*+i*y* to the variable *Z*=*X*+i*Y*
*Ω* is a periodic function of the form
*f* satisfies *f*(0)=0, so that the mapping tends to the identity for *y*=0. The mapping function *Ω* maps the *computational plane* (*x*,*y*), in which the problem will be solved (figure 1*b*), to the *physical plane* (*X*,*Y*), in which the problem is initially stated (figure 1*a*). It transforms the Laplacian operator in equation (2.1*a*) into another Laplacian operator
_{xy} is the Laplacian operator with respect to the (*x*,*y*) variables. Defining
*a*) in terms of the (*x*,*y*) variables results in another Helmholtz equation
*n*(*x*,*y*) acts as an equivalent spatially varying index. Importantly, as a conformal mapping locally conserves angles, the boundary condition (2.1*b*) simply reduces to
*n* that we obtain here is well suited for a multimodal method and for a simple asymptotic approach.

The total field *Φ* is now split into the sum of an incident field *φ*_{in} and of a scattered field *φ*
*a*), it also satisfies equation (2.6). Owing to the 2*π*—periodicity of the mapping function in the *x*-direction, and to the form of the incident field, we assume that the total field, and in particular the scattered field, are quasi-periodic in the *x*-direction
*X*,*Y*), that is, the scattered field represents outgoing waves (this condition will be defined more precisely in §3). The problem (2.11) can be viewed as a waveguide in the *y*-direction with periodic boundary condition and is solved with a multimodal method [11,13].

## 3. Multimodal resolution

In this section, a method to solve equations (2.11), based on the multimodal method, is presented. The scattered field *φ* is first decomposed as a Fourier series
*α*_{n}=*α*_{0}+*n*. If we define the vector *a*) on the basis of the *g*_{n} functions yields the second-order equation
*Γ*_{n,m} coefficients are due to the varying index and are responsible for the coupling between the Fourier modes. To rewrite equation (3.3) in the form of a first-order ordinary differential equation, we introduce the derivative of *φ* in the direction of the waveguide
*φ* in equation (3.1), we expand it as a Fourier series

This equation cannot be directly integrated numerically because of the presence of evanescent modes which are responsible for exponential divergence of the error [11,13]. Besides, it cannot be integrated as an initial value problem because the boundary conditions are imposed, respectively, at *y*=0 and at *φ* and as *Γ* is zero at infinity, the initial condition on Y is simply

Once the admittance matrix is known, to compute the scattered field *φ*, equation (2.11*c*), the boundary condition at the wall, has also to be projected on the modes *g*_{n}. This projection in the computational plane yields
** φ**(0) can directly be obtained, and the field can be computed everywhere.

### (a) Numerical integration

In practice, the infinite series has to be truncated. The question of the convergence depending on the number of modes retained in the truncation is discussed later in §6. The initial condition on the admittance is imposed at a point *y*=*L*, where *L* is sufficiently large so that *Γ*(*L*) can be neglected in comparison with K. Since we use a mapping that tends exponentially to the identity, *L* can have a reasonable value.

The admittance matrix is computed numerically with a Magnus–Möebius scheme [13]. To keep the paper self-contained, we remind here the principles of the numerical integration presented in [13]. The *y*-axis is discretized backward into *M* points *y*=*L* to the boundary at *y*=0, so that *y*=*L* so that *Ω* depends on the order retained for the scheme. For the order 2 scheme, its expression is
_{j} are the elements of the matrix exponential

Once the admittance matrix Y is computed, we move on the initial condition on the scattered field (3.11) that is computed numerically as it may be very difficult to compute it analytically. The *Chebfun* package [20], based on Chebyshev expansions [21], is used here and it provides a rapid and accurate computational tool to evaluate the integral in equation (3.11). With the computed value of Y(0) and the boundary condition (3.11), a complete set of initial values for ** φ**(0) and

**(0) at**

*ψ**y*=0 is obtained. Eventually, the field is calculated by again using scheme (3.12) but in the reverse direction from

*y*=0 to

*y*=

*L*. Note that the expression of the exponential matrices in (3.12) has already been computed during the integration of Y and, as said before, corresponds in fact to the resolvent operator of the problem (3.7). It is therefore very convenient to store in memory this resolvent operator, as it allows us to compute easily the solution to any other initial condition.

For *y*>*L*, as there is no more coupling between the modes, the scattered field becomes a sum of plane waves
*y*=0 and *y*=*L*.

## 4. Conformal mappings

Consider, as an example, the simple mapping
*X*=*π* when *a* is real, as in figure 2*a* for which the value *a*=0.9 has been used. The parameter *a* must fulfil |*a*|<1 as otherwise the mapping ceases to be bijective. As explained in the preceding section, the conformal mapping is a tool to transform a complicated geometry into a simpler geometry. In this case, it maps the straight grid in the computational (*x*,*y*) plane in figure 2*b* into the deformed grid in the physical (*X*,*Y*) plane in figure 2*a*. The rough boundary from figure 2*a* is the image of the straight line *y*=0 from figure 2*b*. The price we have to pay for this simplification is the introduction of a space-dependent refractive index *n* in the Helmholtz equation (2.11*a*), that leads to equation (2.11*a*). This refractive index *n* is represented in figure 2*b*. It can be observed that higher values of *n* are localized near the points on the boundary which are farther, near *X*=0 and *X*=2*π*.

The use of conformal mappings allows us to describe a large variety of geometries. For instance, figure 3 shows some geometries that can be obtained with simple conformal mappings. The method described in the preceding section can then be applied directly using these mappings.

However, given a particular boundary geometry, the major difficulty of the method is to find the appropriate mapping to model the boundary. Some methods exist to determine the mapping depending on the type of the boundary. Floryan proposed a method which can handle boundaries made of successive straight segments [22]. In this paper, we focus on a particular mapping based on a Fourier series
*Z*=*X*+i*h*(*X*). Some methods have been developed to determine the *μ*_{p} coefficients. One of these methods [23] relies on three steps: mapping the flat boundary into the unit circle, then mapping this unit circle into a rough circle that is eventually mapped into the rough surface. In the following, we use the iterative method proposed more recently by Vandembroucq & Roux [24,25], which we reproduce hereafter. This method has the advantages of being very simple yet efficient. The problem to solve is to determine the *μ*_{p} coefficients such that
*a*
and
*b*
An order 0 component *μ*_{0}, responsible for a constant shift of the mapping, could also be considered, but we do not take it into account here without loss of generality for the sake of simplicity. Equation (4.4*b*) shows that the *μ*_{p} can be obtained by computing the Fourier transform of *h*(*X*(*x*)). The problem is that we know the function *X*↦*h*(*X*) but we do not know *x*↦*h*(*X*(*x*)). However, for moderate *μ*_{p} values, equation (4.4*a*) can be first approximated by *X*≃*x*, and the *μ*_{p} can be computed from the Fourier transform of *h*(*x*). With these coefficients, a better approximation of equation (4.4*a*) can be used to compute iteratively higher order approximations. The algorithm proposed by Vandembroucq relies on the discretization of the segment [0 2*π*] in 2*N* sampling points *x*_{j}=*jπ*/*N* with 0≤*j*≤2*N*–1 and takes advantage of the fast Fourier transform to compute successive approximations of the images (*X*(*x*_{j},0),*Y* (*x*_{j},0)) of the flat boundary points (*x*_{j},0). This algorithm can be rephrased very simply in terms of the hilbert () function present in toolboxes for computing languages such as Matlab or Octave. This function actually consists in performing a Fourier transform, cancelling half of the signal, and then performing the inverse Fourier transform. Hereafter is given the Matlab/Octave code corresponding to the algorithm of the iterative method:

x = ( 0 : 2* N−1 ) ' *

**pi**/ N ;X = x ;

**for**m = 1 : NB_ITERATIONSY = h( X ) ;

w = hilbert ( Y ) ;

X = x −

**imag**( w ) ;**end**c =

**fft**( w ) ;mu = c (1:N) * 1 i / 2 / N ;

Eventually, the vector mu contains the sought coefficients *μ*_{0}⋯*μ*_{N−1}. To converge, this algorithm requires that the slope of the boundary remains less than 1
*Ω* into two (or more) mappings *Ω*=*Ω*_{1}°*Ω*_{2} where one of these mappings fulfils condition (4.5) so that this method can be applied.

Note that the crest mapping (4.1), of which an example is represented in figure 2, is actually the simplest form of the general mapping (4.2). Indeed, it corresponds to the case where all the *μ*_{p} are null except *μ*_{1}=−i*a*. Figure 4 shows two other examples of geometries modelled with this mapping and that are used in the following section. The first one, in figure 4*a*, is a randomly generated rough surface, whose coefficients *μ*_{p} are given in table 1. The second one, in figure 4*b*, is an almost parabolic reflector which is composed of a parabola of equation *h*(*X*)=(*X*−*π*)^{2}/4*f*_{0}−(2*π*−*d*)/4*f*_{0}+*d*(*π*−*d*)/4*f*_{0} for *X*∈[*d*;2*π*−*d*] and of a polynomial prolongation to obtain a smooth surface. The distance *d* represents the width of the parabolic section in the *X* direction, and *f*_{0} is the focal length. Coefficients *μ*_{p} are given in table 2.

In the following, we will focus on four particular geometries, a geometry being the given set of a conformal mapping and of its parameters. Three of them are obtained with the mapping (4.2) with different parameters. We will refer to the geometry in figure 2 with *a*=0.9 (or *μ*_{1}=−0.9*i*) as to the *crest geometry*. Geometries in figure 4*a*,*b* will be named *rough geometry* and *parabolic geometry*, respectively. The geometry from figure 3*a*, which is obtained with the mapping
*a*=0.02 and *b*=1.05, will be called the *cavity geometry*. The expression of the coupling matrices corresponding to these mappings are given in the appendix.

## 5. Applications

As a first example of application, we consider the simple case of the crest geometry. The modulus of the total field *Φ* is represented on figure 5*a* for *k*=20.3 and for an incident field given by (2.9) with *α*_{0}=20 (*α*_{0}/*k*=0.99). This geometry corresponds to a high-frequency wave with near grazing angle. It can be noted that this method efficiently predicts the shadowing effect at the bottom right of the crest. The case of the parabolic reflector from figure 4*a* is shown in figure 5*b* at *k*=40.3 for an incident wave with normal incidence (*α*_{0}=0). It can be seen that the scattered wave focuses at the point (*π*/2,*f*_{0}). Finally, the case of the rough surface proposed in figure 4*b* is treated in figure 6 with *k*=5.3 and *α*_{0}=4 (*α*_{0}/*k*=0.76), in a case where the incident wave and the surface roughness have comparable typical length scales.

## 6. Numerical considerations

The accuracy of the result given by the proposed method depends on different parameters and is subject to several error sources. The first source of error comes from the initial condition (3.10) on the admittance matrix Y, where this impedance matrix is approximated by the characteristic admittance Y_{c} at *y*=*L*. The error due to this approximation decreases exponentially as *L* increases. In all the figures shown before the value *L*=4*π* has been used and gives very good results. The second source of error comes from the Möbius Magnus integration scheme and depends on the strategy retained for the discretization of the *y*-axis. A strategy with an adaptive step is detailed in [15]. The third source of error is due to the numerical integration of the initial condition on the scattered field (3.11). In the present work, the Chebfun package has been used. Its precision has been estimated by comparing *ψ*_{n}(0) and *ψ*_{−n}(0) in cases of symmetric geometries with normal incident wave, for which, theoretically, *ψ*_{n}=*ψ*_{−n}. It appears that the precision is therefore of the order 10^{−15}, which is the computer precision.

Eventually, the last and most important source of error is due to the mode truncation. The number of Fourier modes which has to be retained depends on both the frequency and on the mapping used. In all cases, all the propagative modes have to be kept. But for the calculation to converge, a certain number of evanescent modes has also to be considered. A pragmatic way to see if such a method has converged is to compare results obtained with different number of modes in mode truncation [11,26].

To study the effect of the truncation on the convergence, we consider as a test case the crest geometry. We consider here a frequency *k*=2.3, and a vertical incident wave *α*_{0}=0, with *L*=8*π* and *n*_{y}=8500 points on the *y*-axis. At this frequency, there are five propagative modes (from *n*=−2 to +2). We performed four computations with different modal truncations, first retaining only these five propagative modes, and then adding 2, 4 and 56 evanescent modes. The contours of the modulus of the total field are represented in figure 7*a*–*d* for each modal truncation, respectively. No particular difference can be noted between figure 7*c*,*d*, meaning that the method converged with four evanescent modes.

To provide a more careful analysis, we propose here an estimation of the rate of convergence depending on the number of modes. The contribution of a mode relies on the *φ*_{n} coefficients, which are by definition
*φ* is the exact solution. Integrating by parts repeatedly equation (6.1), and using the definition of the modes and the pseudo-periodicity boundary condition (2.11b) show that
*φ* and all its derivatives are quadratically integrable functions, then
*x* variable. This means that the *φ*_{n} coefficients converge more rapidly than every power of *n*, that is, they converge exponentially. Figure 8 displays the modulus of the modal amplitudes at the origin |*φ*_{n}(0)|, showing the exponential convergence (note the logarithmic scale for the *φ*_{n} on the figure). From mode 15, the amplitudes reach the level of the error on the initial condition and become dominated by roundoff errors (here around 10^{−16}).

### (a) Energy flux

A last tool to evaluate the accuracy of the method is to check the conservation of the energy flux across a surface *y*=const. [27]. Thanks to the modal formulation and to the conformal property of the mapping, this energy flux can be very easily obtained numerically
*W*
*Γ* matrix is self adjoint, as can be seen from equation (3.4) defining its coefficients. Furthermore, K^{2} is a real diagonal matrix so that *W*′(*y*)=0. It shows that the multimodal resolution preserves the conservation of the energy flux *W*. It is true even when the expansion is truncated at a finite number of modes as can be seen in equation (6.6).

## 7. Quasi-trapped modes

As mentioned earlier, the Riccati equation (3.9) is known to possess movable (quasi-) singularities. As pointed out in [13], these singularities turn out to be an efficient way to find trapped modes as is explained hereafter. Such quasi-singularities can be found for instance in the cavity geometry (figure 3*a*) with *k*=2.9 at *y*_{c}=0.257, *y*_{c}=0.0050 and *y*_{c}=0.0016, as can be seen from figure 9*a*, which shows the real part of *Y* _{10,10}, one of the coefficients of the impedance matrix that is particularly affected by the quasi-singularity. At such points, a peak of the modulus of one eigenvalue *λ*_{c} occurs. Considering *ϕ*_{c} be the eigenvector associated to that singular eigenvalue *λ*_{c}, then it satisfies by definition
*b*, showing the modulus of the highest eigenvalue ** φ** to equation (3.3) but with the initial condition

*y*=

*y*

_{c}(and

**and ∂**

*φ*_{y}

**are finite), we obtain that**

*φ**y*=

*y*

_{c}, while not being the trivial identically null solution (because ∂

_{y}

**is not null). This means that this field (quasi-) fulfils a Dirichlet boundary condition on**

*φ**y*=

*y*

_{c}. As this field is by construction a solution to equation (2.6) only composed of outgoing waves, it is a trapped mode of the flat surface

*y*=

*y*

_{c}(or surface wave). This flat surface

*y*=

*y*

_{c}in the computational plane corresponds to a rough surface (or at least a non-flat surface) in the physical plane ((

*X*,

*Y*) variables). Therefore, the corresponding solution in the physical plane is a trapped mode of that rough surface. Eventually, we obtain that for each (quasi-) singularity, it is possible to exhibit a trapped mode of the rough surface that is the image of the line

*y*=

*y*

_{c}by the mapping function

*Ω*. The three quasi-trapped modes associated to the quasi-singularities from figure 9

*a,b*are represented in figure 10. We can see that, indeed, the energy is localized inside the cavity.

Scanning the *y*-axis looking for singularities constitutes one method to find such geometries having trapped modes. Note that to apply this method, the frequency *k* has been chosen, so that this method also permits us to find these geometries as a function of a given frequency, rather than finding the frequencies at which trapped modes exist for a given geometry.

In the same manner, the (quasi-) zeros of the eigenvalues of Y can be used to find geometries having trapped modes with respect to the Neumann boundary condition (rigid wall). In that case, we consider this time an eigenvalue *λ*_{0} which becomes (quasi-) zero at *y*=*y*_{0}
*ϕ*_{0} be its associated eigenvector, then, the solution ** φ** with initial condition

*y*=

*y*

_{0}. As these (quasi-) zeros of Y corresponds to (quasi-) singularities of the impedance matrix Z=Y

^{−1}, the same method can be applied by considering the highest eigenvalue

*ν*

_{max}(

*y*) of all eigenvalues

*ν*

_{m}(

*y*) of the Z matrix. Figure 11

*a*–

*d*shows three examples of these trapped modes obtained by using the cavity geometry, the crest geometry and the rough geometry, respectively.

To have a better understanding of how to find these trapped modes, it is interesting to follow the associated singularities of the Y and Z matrices depending on *k* and *y*. Figure 12*a* shows the evolution in the (*k*,*y*) plane for the cavity geometry of the modulus of the highest eigenvalues *λ*_{max}(*y*) of the admittance matrix Y. Quasi-singularities correspond to dark areas and form structures resembling curved lines. They also evolve towards lower values of *k* as *y* is close to zero. The integration of the Riccati equation (3.9) for *k*=2.9 and *α*_{0}=0 consists in fact in examining the (*k*,*y*) plane along the vertical line *k*=2.9 (dashed line). Doing this, we found the three quasi-singularities from figure 9*a*,*b*, here represented by cross marks. Note that a third dimension could be added by considering the variation of *α*_{0}. Similar phenomena happen when considering the highest eigenvalue *ν*_{max}(*y*) of the impedance matrix Z in figure 12*b*. The cross mark represents the quasi-singularity corresponding to the quasi-trapped mode from figure 3*a*. The singularities of the impedance matrix Z also form curved line structures. They however seem to have a particular behaviour as these lines connect to the singularities due to the cut-off frequencies of the transverse modes of the waveguide for *y*=const. to curved lines so as to match the rough boundary at *y*=0. These horizontal straight lines, which are 2*π* long, have a lower length than their image. The more *y* is close to zero, the higher is the difference of length between the straight lines and their images. As a result, the cut-off frequencies, and the associated singularities, is supposed to move to lower frequencies when *y* tends to zero. This is when one of these singularity can cross the vertical path of integration *k*=const. The fact that the frequencies of the trapped mode tend to the cut-off frequencies might be not surprising. Indeed, a similar phenomenon is known in water channels. For *α*_{0}*d*=*π* (*d* being the width of the channel), an obstacle at the centre allows the existence of water waves trapped modes. For a cylindrical obstacle, it has been shown that as the diameter of the obstacle tends to zero, the frequency of the trapped mode tends to the cut-off frequency of one of the mode of the channel [28].

## 8. Asymptotic approach

We propose in this section an asymptotic approach based on a small roughness approximation, to obtain the first order of the scattered field. For this, we consider the case of the mapping (4.2). Its particular form based on a Fourier series makes that the coupling matrix has a relatively simple form, that will play a crucial role in the derivation of the asymptotic solution. The mapping (4.2) is first rewritten under the form
*ε*>0 represents the amplitude of the roughness. Considering a small roughness for which *ε*≪1, the coupling coefficients *Γ*_{n,m} (see equation (A 3) in appendix A(a)) are at first order in *ε*
*φ*_{n} together with the initial condition (3.11), becomes at leading order
*λ*_{p}=*β*_{0}+i*p* and *γ*_{0}=*α*_{0}−i*β*_{0}.

We therefore look for a solution under the form of an *ε* power series
*A*_{n} are obtained from the initial condition (8.4) retaining again terms of order 1. This yields

Now, in order to obtain the order *ε* component *ε*. This yields
*B*_{n} are determined by the initial condition. To do this, we identify the derivative of equation (8.10) with the component of order *ε* in equation (8.4). This yields

Thus, the final form of the scattered field is given by
*λ*_{n} have a positive imaginary part, the two exponentials involving the *λ*_{n} vanish in the far field.

Figure 13 shows a comparison of the total field computed by this perturbative method and by the numerical method developed in the first part of this paper. We considered here a crest geometry *Ω*(*z*)=*z*−*ε*i e^{iz} with *ε*=0.1. We can see that the perturbative method is able to give a good result in comparison to the numerical method. Figure 14 corresponds to the *ε*=0.3, where some discrepancies appear.

Eventually, the asymptotic solution in equation (8.12) has a rather simple form. The fact that the mapping considered and the modal projection are in fact Fourier series, led to a number of simplifications in the derivations of this asymptotic solution. Eventually, the scattered field *φ* is itself a Fourier series in the *x*-direction, whose coefficients have a very simple dependence on the *η*_{p}. Therefore these *η*_{p} coefficients can be simply extracted from the scattered field by computing its Fourier transform along *x*. This could constitute a sequel to this work.

## 9. Concluding remarks

A method to compute the scattering by an infinite periodic rough boundary has been presented. The method we proposed in this paper requires no assumption on the relative orders of the characteristic lengthscale of the boundary geometry and of the wavelength. The central idea of the proposed method is to use a conformal mapping to move the mathematical difficulty from the boundary condition into the bulk equations. Eventually, the problem can be seen as a straight waveguide with periodic boundary condition, leading to a Helmholtz equation with a spatially varying index. This new problem is solved with the multimodal impedance matrix method, resulting in a Riccati equation governing the variation of the impedance matrix along the waveguide. This formulation allows an exponential convergence with the number of modes and conserves the energy flux along the waveguide. Three illustrations of our method have been given. Incidentally, the formulation presented in this paper also allowed to develop a simple perturbative approach, that gives the first order of the scattered field under a small roughness assumption.

The integration of the admittance matrix, which is the core of the multimodal method, is subject to quasi-singularities of the admittance matrix. These quasi-singularities turn to be very useful to find geometries having quasi-trapped modes. They form structures resembling curved lines in the (*k*,*y*) plane. A quasi-trapped mode can then be found when one of these curved lines crosses the line *k*=const. corresponding to the considered frequency. These quasi-singularities move towards lower values of *k* near the boundary due to the strong effect of the mapping in this region. The particularity of this method is that it provides the geometry having a quasi-trapped mode as a function of the frequency, instead of finding the frequencies at which a quasi-trapped mode exist inside a particular geometry.

## Authors' contributions

G.F. carried out the calculations and simulations and drafted the manuscript. V.P. proposed the main ideas, supervised the study, and helped draft and correct the manuscript. All authors gave final approval for publication.

## Funding statement

V.P. gratefully acknowledges the Agence Nationale de la Recherche for financial support (PLATON project, ANR-12-BS09-0003-03). G.F. acknowledges Ministère de la Recherche for funding his thesis through the AMX fellowship.

## Conflict of interests

We have no competing interests.

## Appendix A. Coupling matrices

**(a) General mapping based on Fourier series**

Mapping function:

**(b) Cavity mapping**

See figure 3*a*.

Mapping function:

- Received October 10, 2014.
- Accepted December 22, 2014.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.