New preconditioners for the Laplace and Helmholtz integral equations on open curves: analytical framework and numerical results

Helmholtz wave scattering by open screens in 2D can be formulated as first-kind integral equations which lead to ill-conditioned linear systems after discretization. We introduce two new preconditioners in the form of square-roots of on-curve differential operators both for the Dirichlet and Neumann boundary conditions on the screen. They generalize the so-called “analytical” preconditioners available for Lipschitz scatterers. We introduce a functional setting adapted to the singularity of the problem and enabling the analysis of those preconditioners. The efficiency of the method is demonstrated on several numerical examples.


Introduction
For the resolution of wave scattering problems, a well-established approach is the boundary element method (BEM), which involves the discretization of integral equations and leads to smaller linear systems than finite element methods, but to dense ones. However, they can be solved efficiently by combining iterative methods such as GMRES [43], with matrix compression methods, e.g. the fast multipole method, see [21] and references therein. The number of GMRES iterations can be very large for highly refined meshes and large wave numbers. Reducing this number of iterations, either by using well-conditioned integral equations or by finding preconditioners for the linear systems, has attracted a lot of attention since more than two decades. In the case of an obstacle with a C ∞ smooth boundary Γ , general methods to build efficient preconditioners are known [1,2,4,14,24,47], which in most cases can also be applied in theory and/or in practice if the scatterer is only assumed to have Lipschitz regularity. Among them, one efficient approach is the Generalized Combined Source Integral Equation (GCSIE) method [2] invented by Levadoux in [30]. It involves the inversion of the discretized GCSIE operator where S k and D k stand for the classical single-and double-layer potentials on Γ , I d is the identity operator andΛ k is an approximation of the exterior Dirichlet-to-Neumann (DtN) map Λ k for the Helmholtz equation. In the ideal caseΛ k = Λ k , G k is the identity and, consequently, whenΛ k is a compact perturbation of Λ k , G k is a compact perturbation of the identity. The Galerkin discretization of such operators leads to well-conditioned linear systems for which the GMRES method converges super-linearly [9,23,47]. To build candidates forΛ k , a generic tool is the theory of pseudo-differential operators (see e.g. [26]). This has been succesfully applied by Antoine and Darbas in [4] who proposed the choicẽ based on low-order expansions of the symbol of Λ k . The authors also introduce an efficient scheme to approximate the square root operator in (1) relying on Padé approximants. Their numerical results clearly demonstrate that the preconditioning performances are independent of the discretization parameters and quite robust in k.
The full theoretical analysis of this last feature still remains an open problem. When the scattering object is a "screen", that is a curve in 2D and an open surface in 3D with boundary (therefore not the boundary of a Lipschitz domain), the robustness of the preconditioners in the discretization parameters and in k is lost in practice. In this work, we propose a generalized version of the square-root preconditioners of Antoine and Darbas on screens in 2D that overcomes this limitation. To fix notations, let us write the classical first-kind integral equations corresponding to the Dirichlet and Neumann problems as where S k and N k are the classical single-and hypersingular layer potentials on the screen Γ and u D and u N are smooth right-hand sides. We define the operators where ∂ τ is the tangential derivative on Γ and ω(x) is defined for x ∈ Γ as the square root of the product of the distances from x to the end-points of Γ . Those are optimal preconditioners (i.e. the number of iterations are bounded uniformly with respect to the discretization parameters for fixed k), which additionally perform robustly with respect to the wavelength parameter k. We present numerical evidence supporting this claim. The design of P k and Q k , and the proof that they provide parametrices of S k and N k relies on two new classes of pseudo-differential operators on open curves. This is treated in detail in [6]. To our knowledge, recent literature mainly contains two important directions for the problem of preconditioning for Helmholtz screen scattering problems, and the present work is connected to both of them. The first one, see e.g. [24,25,27,39,50], builds on recently found closed-form expressions for the inverses of the Laplace layer potentials (k = 0) on straight screens in 2D and 3D (recovered by a unified approach in [20]). Those explicit formulas allow to construct "optimal preconditioners", i.e. that lead to provably well-conditioned systems independently of the mesh size, solving the so called "duality mismatch". By usual domain parametrization and compact perturbations this also gives optimal preconditioners for the Helmholtz (k > 0) layer potentials and for more general screens. However, the preconditioners do not perform well uniformly in k, and in those works, no attempt is made at obtaining robustness with respect to k. In this work, we show that when Γ is a segment, P 0 and Q 0 coincide with those exact inverses (up to minor modifications), and that for large k, P k and Q k greatly outperform P 0 and Q 0 .
The second direction, pursued by Bruno and Lintner [13] is a generalization of the Calderón preconditioners available for Lipschitz scatterers. In their work, some weighted versions of the Helmholtz layer potentials, S k,ω and N k,ω , are considered and it is shown that S k,ω N k,ω is a second-kind operator leading to another optimal preconditioning strategy that is also very efficient for large values of k in practice. The weighted layer potentials S k,ω and N k,ω also play a key role in the present study. Furthermore, we compare numerically our preconditioning method to that of Bruno and Lintner. In our implementation, we find that they lead to similar numbers of iterations, but our preconditioners are significantly faster to evaluate due to their quasilocal form.
The outline of the paper is as follows. We use the first section to fix notations and recall some important results concerning integral equations on screens. In the second section, we introduce the preconditioners for the Laplace problem on a straight screen. The formulas are generalized to the Helmholtz equation and non-flat arcs in the third section. In Sect. 4, we describe the weighted Galerkin setup that we use to discretize the integral equations. In Sect. 5, we detail the discretization of the preconditioners introduced in Sects. 2 and 3. Finally, in Sect. 6, we show the performance of our preconditioners in a variety of cases and compare them to other ideas of the literature.
where γ : [−1, 1] → R 2 is an injective C ∞ function. Let k ≥ 0 be the wavenumber and G k be the Green's kernel where H (1) 0 is the Hankel function of the first kind and order zero and |x| is the Euclidean norm of x. The classical single-layer boundary integral operator, denoted by S k , is defined for φ is the space D(Γ ) of smooth compactly supported functions on Γ by where dσ is the uniform measure on Γ . The Helmholtz hypersingular boundary integral operator N k : where φ and ψ denote the arclength derivatives of φ and ψ, n is a smooth unit normal vector on Γ and D (Γ ) is the set of distributions on Γ .
Sobolev spaces on the screen are denoted by H s (Γ ) andH s (Γ ) for s ∈ R and defined as in [33, chap. 3]. Specifically, considering a smooth domain Ω such that Γ ⊂ ∂Ω, a distribution u on Γ is said to be in H s (Γ ) if there exists an U ∈ H s (∂Ω) such that u = U |Γ . Furthermore, u ∈H s (Γ ) if u is in H s (∂Ω) and supp u ⊂ Γ .
Proposition 1 (see [48,Thm 1.8] and [49,Thm 1.4]) The operator S k has a continuous extension which is bicontinuous when k = 0. When k = 0, it is bicontinuous if and only if the logarithmic capacity of Γ (see e.g. [40]) is not equal to 1, and positive definite when the logarithmic capacity is strictly less than 1. Similarly, the operator N k can be extended continuously as and this extension is bicontinuous, and positive definite when k = 0.
For convenience, we assume that the logarithmic capacity of Γ is strictly less than 1, bearing in mind that this condition can be enforced by rescaling, and that our main focus is the case k = 0 where the capacity does not play a role.
In this work, we are concerned with the resolution of the integral equations where u D ∈ H 1/2 (Γ ) and u N ∈ H −1/2 (Γ ). They are related to the Helmholtz equation in R 2 \ Γ with prescribed Dirichlet boundary data u D ("sound-soft") or Neumann boundary data u N ("sound-hard") respectively. Due to the singularity of the manifold Γ , λ and μ have edge singularities even if u D and u N are arbitrarily smooth. More precisely, Costabel and where E 1 = γ (−1) and E 2 = γ (1) are the two end-points of Γ .
In addition to the singular nature of the problem, the Galerkin discretization of first-kind integral equations is known to produce ill-conditioned linear systems, see [45,Sec. 4.5]. The aim of this paper is to introduce a formalism that enables to resolve the singularity and provide efficient preconditioners for the linear systems.

Laplace equation on a segment
We start with the particular case where Γ is the segment Γ = [−1, 1] × {0} and furthermore k = 0. The associated integral equations, with logarithmic kernels, are the object of a considerable number of works, for instance [5,19,29,36,44,51].
Recently, exact iverses of those operators were exhibited [27], through explicit variational forms (Prop 3.1 and 3.3) and identities involving tangential square-roots of differential operators with boundary conditions at the end points (Prop 3.10 and subsequent remark). In Gimperlein et al. [20], those results are recovered with a different method, and it is shown that S 0 and N 0 are related to fractional powers of the Laplace-Beltrami operator on the line {(x, 0) | x ∈ R }, by considering the natural extensions of functions by 0 outside Γ .
Here we give a new expression for the inverses of S 0 and N 0 (Theorem 1 and Theorem 2), in the form of square roots of suitably weighted tangential differential operators. Those formulas do not involve boundary conditions at the end points nor extensions by 0 outside Γ . The proofs of the identities are very simple; nevertheless, their consequences in terms of preconditioning have not been exploited so far.

Analytical setting
Like many works on the topic e.g. [28,46], we introduce a simple functional setting based on Cheybyshev polynomials The Chebyshev polynomials of first and second kind [32] are respectively given by We denote by ω the operator ω : 1 − x 2 , and by ∂ x the derivation operator. The Chebyshev polynomials satisfy the ordinary differential equations which can be rewritten in divergence form as where we emphasize that ∂ x ω f should be understood as the composition of the operators ∂ x and ω applied to f , that is ). Both T n and U n are polynomials of degree n, and form complete orthogonal families respectively in the Hilbert spaces see e.g. [32] and in particular Thm. 5.2.
As a consequence, any u ∈ L 2 1 ω and v ∈ L 2 ω can be expanded in where the sums converge respectivey in L 2 1 ω and L 2 ω . The "Fourier-Chebyshev" coefficientsû n andv n are given bŷ with the inner products Those properties can be used to define Sobolev-like spaces.
Definition 1 For all s ≥ 0, we define Endowed with the scalar product T s is a Hilbert space for all s ≥ 0. Similarly, we set which is a Hilbert space for the scalar product (1 + n 2 ) sǔ nvn .
We write u 2 T s := (u, u) T s and u 2 U s := (u, u) U s . One can extend the definitions of T s ad U s for s ∈ R, in which case they form interpolating scales of Hilbert spaces, see [ For s = ± 1 2 , those spaces were analyzed (with different notation) e.g. in [27] and verify with equivalent norms, that is to say, densely defined respectively on T s and U s , are self-adjoint with compact resolvent and purely discrete and real spectrum.
Proof For k = 0, the proof is elementary, and relies on the families of eigenfunctions T n and U n . We omit this part. For k = 0, let us show that −k 2 ω 2 : T s−2 → T s−2 (resp. U s−2 → U s−2 ) is symmetric, closed, and compact from T s to T s−2 (resp. U s → U s−2 ). For all n ∈ N, one has the formula xT n = T |n−1| + T n+1 2 from which the symmetry and boundedness (hence closedness) of φ → xφ in T s−2 follows. By composition and linearity, this extends to ω 2 . By similar arguments, ω 2 is also closed and symmetric in U s−2 . Since the inclusions T s ⊂ T s−2 and U s ⊂ U s−2 are compact, we can apply Weyl's theorem to finish the proof.

Remark 1
Let √ · stand for the operator defined as usual for x > 0 and by √ x = i √ −x for x < 0. This coincides for example with the square root defined on the complex plane with branch cut on the negative imaginary axis. By what precedes, the operators can be meaningfully defined via functional calculus. Using the explicit eigenvalues of −(ω∂ x ) 2 and −(∂ x ω) 2 and the relative compactness of ω 2 , it is easy to show that they are continuous respectively from T s to T s−1 and U s to U s−1 for all s ∈ R.

Lemma 3
There exists an increasing sequence k D n (resp. k N n ) tending to +∞ of wavenumbers such that the operator −(ω∂ has a bounded inverse if and only k is not one of the k D n (resp. k N n ).
The proof also relies on the spectral theorem.

Laplace single-layer equation
We start with the single-layer integral equation The following result is fundamental: Using the decomposition of g on the basis (T n ) n , we see at once that the solution λ to equation (8) admits the expansion As a corollary, we obtain by an alternative proof the result of Costabel et al. cited in Proposition 2 in this particular case: , the solution λ to the equation , then g ∈ T ∞ , and by equation (9), from which we deduce that α also belongs to Following [13], we introduce the weighted single layer operator as the operator that appears in Lemma 4.

Definition 2
Let S 0,ω be the weighted single layer operator defined by Clearly, S 0,ω extends to an isomorphism from T s to T s+1 for all s ∈ R.
Theorem 1 There holds in the space of linear operators T s → T s−1 for all s. As a consequence, as announced. In particular, the equality holds from T −1/2 to T 1/2 . Writing and exploiting the equalities T 1/2 = H 1/2 and T −1/2 = ωH −1/2 (Γ ), we deduce the second assertion.

Laplace hypersingular equation
We now turn our attention to the hypersingular equation Similarly to the previous section and again following [13], we consider the weighted version of the hypersingular operator N 0,ω := N 0 ω. We can get the solution to equation (11) by solving and letting μ = ωβ. We now show that N 0,ω can be analyzed using this time the spaces U s .

Lemma 5 For any β, β , one has
Proof We use the well-known integration by parts formula valid when u and v are regular enough and vanish at the endpoints of the segment. For a smooth β, we thus have which implies the claimed identity.
Lemma 6 For all n ∈ N, there holds Proof From the identity ∂ x T n+1 = (n + 1)U n and Eq. (4), one has Therefore, by Lemma 5 which implies the result.

Corollary 2
The operator N 0,ω is bijective from U s to U s−1 for all s ∈ R.

Helmholtz equation
In this section, we aim at generalizing the preceding analysis to the case of Helmholtz equation (k > 0), still for a flat screen Γ = [−1, 1] × {0}. Let S k,ω := S k 1 ω and N k,ω := N k ω. The following commuting relationship holds: where we use the notation ω y and ∂ y to emphasize the dependence on the variable y. Thus, Since G k is a solution of the Helmholtz equation, we have for all Note that y 2 − x 2 = ω 2 x − ω 2 y so the first term vanishes and we find All the previous computations can be justified rigorously by breaking the integral into y < x and y > x. A careful analysis reveals that no Dirac mass appears in the previous formula. We conclude by density of T ∞ in T s .
There also holds the following identity: The proof can be found in [8, Chap. 2, Thm 2.2]. Those commuting relationships imply that the operators S k,ω and N k,ω share the same eigenvectors as, respectively, The eigenfunctions of the operator −(ω∂ x ) 2 − k 2 ω 2 thus provide us with a diagonal basis for S k,ω . They are the solutions of another Sturm-Liouville problem There exists a discrete set of values a 2n (q) for which this equation possesses even and 2π periodic solutions, which are known as the Mathieu cosine functions, and usually denoted by ce n . Here, we use the notation ce k n to emphasize the dependency in the parameter k = √ 2q of those functions. The normalization is taken as The Mathieu cosine functions are L 2 orthogonal: so that any even 2π periodic function in L 2 (−π, π) can be expanded into the functions ce n , with the coefficients obtained by orthonormal projection. Setting in analogy to the zero-frequency case, we have For large n, using the general results from the theory of Hill's equations (see e.g. [37, eqs. (21), (28) and (29)]), we have the following asymptotic formula for λ n,k : The first commutation established in Theorem 3 implies that the Mathieu cosine functions are also the eigenfunctions of the single-layer operator. (An equivalent statement is given in [12,Thm 4.2], if we allow the degenerate case μ = 0.) A similar analysis can be applied to the hypersingular operator. The eigenfunctions of −(∂ x ω) 2 − k 2 ω 2 are given by where se k n are the so-called Mathieu sine functions, which also satisfy the Mathieu differential equation (13), but with the condition that they are 2π periodic and odd functions.
From all the previous considerations, one could expect that the operators provide compact perturbations of the inverses of S k,ω and N k,ω respectively. This is true but the proof is not as simple as in the Laplace case, due to the fact that no explicit knowledge of the eigenvalues of the operators is available. To state the result precisely, let us introduce the following terminology: Theorem 4 ([6, Thm. 7]) The operators where the equalities are understood in the set of endomorphisms of T s and U s for each s, I d is the identity operator, and K 1 and K 2 are of order −4 respectively in the scales T s and U s .
Since K 1 and K 2 are of negative order, they are in particular compact endomorphisms of L 2 1 ω and L 2 ω respectively, due to the compact inclusions that hold for Therefore, the operators appearing in (14) are of second-kind in those spaces. The extent to which the dependence in k in P k and Q k is optimal is reflected by the next theorem.
Theorem 5 ([6, Thm. 5 ]) Let K be an operator of order 0 in the scale T s and let Δ K be defined by Then where L is an operator of order −2 in the scale T s .
Taking K = 0, we see that P 0 and Q 0 are also compact equivalent inverses of S k,ω and N k,ω , however up to a less regularizing remainder than P k and Q k . A clear link between the order of the remainder and the numerical performance of the preconditioner remains to be elucidated. We conjecture that a smoother remainder leads to better performances and especially, more robustness with respect to the parameter k. This is strongly supported by our numerical results presented in Sect. 6.
We note that all the previous analysis (except the commuting relationships) carries over to the more general case of a C ∞ non-intersecting open curve Γ . The weight ω is replaced by where Z 1 and Z 2 are the end-points of the curve. P k and Q k are replaced by where ∂ τ is the tangential derivative.
The previous theoretical analysis suggests that we use P k and Q k as operator preconditioners for S k,ω and N k,ω . The remainder of this paper is dedicated to testing this idea in practice.

Galerkin method
It is known that the naive piecewise polynomial Galerkin discretization of the nonweighted integral equations with a uniform mesh leads to very slow convergence in terms of the mesh size (the error in energy norm converges in O(  [16,38], Galerkin method with special singular functions [15,48], spectral collocation [46] or Galerkin [28] methods, cosine change of variables [51], and high-order Nyström methods [13]. In this work, we refrain from using spectral methods relying on the explicit knowledge of eigenvectors (here T n and U n ) despite their exponential order of convergence. Instead, we prefer to use a Galerkin space that only uses intrinsic information from the curve, hoping that this will render the method more amenable to generalization.
For this reason, our Galerkin space is the set of piecewise linear functions. They are defined on a non-uniform mesh, which is refined towards the end-points as follows. Let X : [0, |Γ |] be the parametrization of Γ by the arclength, where |Γ | is the length of the curve. We choose the nodes (X i ) 1≤i≤N as where (s i ) 1≤i≤N are such that the value is (almost) independent of i. Such a mesh turns out to be analogous to an algebraically graded mesh with a grading parameter β = 2. That is to say, near an edge of the curve, the width of the i − th interval is approximately (ih) 2 for some constant parameter h. Notice that this modification alone, i.e. using the h-BEM with a polynomial order p = 1, is not sufficient to get an optimal rate of convergence. Indeed, it only leads to a convergence rate in O(h) for the energy norm (cf. [38, Theorem 1.3]) instead of the expected O(h 5/2 ) behavior (to reach such an order of convergence would require β = 5).
The key ingredient to recover optimal convergence, beside the graded mesh, is to use a weighted L 2 scalar product (with weight 1 ω or ω depending on the considered equation), in order to assemble the operators in their natural spaces. We state here the orders of convergence that one gets with this new method, and refer the reader to [8, Chap. 3, Sec. 2.3] for the proofs. To keep the exposition simple, we also restrict our presentation to the case where Γ = [−1, 1] × {0} and k = 0.
In the following, we take the notation h = max i h i where h i are defined in Eq. (15) and use a h subscript to highlight the dependence of the Galerkin matrices and solutions with respect to the discretization. Dirichlet problem. For the resolution of the single-layer equation (8) we use a variational formulation of (10) to compute an approximation α h of α. Namely, let V d h be the Galerkin space of discontinuous piecewise affine functions on the mesh (x i ) 0≤i≤N defined above, and α h the unique solution in V d h of We then compute λ h = α h ω .
where C > 0 is a constant independent of h.
In particular, when u D is smooth, the solution α = ωλ belongs to T ∞ , and we get the optimal rate of convergence of the error in O(h 5/2 ).

Remark 2
From numerical results it seems that this result holds also when replacing V d h by V h the set of continuous piecewise linear functions. All numerical simulations are performed in V h .
Neumann problem. For the numerical resolution of (11), we use a variational form for equation (12) to compute an approximation β h of β, and solve it using a Galerkin method in V h . We denote by β h the unique solution in V h to the variational equation: Then, the proposed approximation for μ, given by μ h = ωβ h , satisfies the following error estimate.
Numerical validation. In fact, estimates in the local weighted L 2 norms can be derived from the previous reasoning, namely: for the Dirichlet problem and for the Neumann problem. We verify those rates numerically. For the Dirichlet problem, we solve two test cases S 0,ω α 1 = u 1 and S 0,ω α 2 = u 2 having the explicit solutions α 1 (x) = ω(x) and α 2 = ω(x) 3 , for adequately chosen right hand sides (rhs) u 1 and u 2 . One can check that α 1 ∈ T s for s < 3 2 and α 1 / ∈ T 3/2 , while α 2 ∈ T 2 . The L 2 1 ω error is plotted in Fig. 1 in each case as a function of the mesh size h. We find that the expected rates O(h 3/2 ) and O(h 2 ) predicted by the theory are precisely observed in practice. Fig. 1 Effective order of convergence of the approximation of the solution α of (10) by the weighted Galerkin method in V h . Two cases are considered where α ∈ T s for all s < 3/2 but α / ∈ T 3/2 (solid line and circles) and α ∈ T 2 (dashed line, crosses). The approximate slopes are displayed above each curve, and correspond to the theoretical convergence rates proved for V d h Similarly, for the Neumann case, we solve a a test case N 0,ω β = u N where the solution β is explicit. We take u N = U 2 the second Chebyshev polynomial of the second kind. The corresponding solution β is proportional to U 2 and thus belongs to U ∞ . The theory therefore predicts a convergence rate of the error in the L 2 ω and U 1 norms respectively in O(h 2 ) and O(h). This behavior is again confirmed by our numerical results, presented in Fig. 2.

Discretization of the preconditioners
Let (φ i ) 1≤i≤n dof be the usual "tent functions" basis of V h . For operators A : T s → T −s and B : U s → U −s , we denote by [A] 1 ω and [B] ω their Galerkin matrices defined by

The operators
introduced in Sect. 2 for k = 0 (Theorem 1 and Theorem 2) and Sect. 3 for k > 0 are at the base of our preconditioning strategy, as P k S k,ω and Q k N k,ω are compact perturbations of the identity (Theorem 4). To define preconditioners for the linear systems, we would ideally need to compute the Galerkin matrices of P k and Q k .

Fig. 2
Effective order of convergence of the approximation of the solution β of (12) by the weighted Galerkin method in V h . Here, β ∈ U 2 and the error is measured in the norms L 2 ω (solid line, circles) and U 1 (dashed line, crosses). The approximate slopes are displayed above each curve and correspond to the theoretical convergence rates Then, it would be possible, following [23], to prove that the condition numbers of the matrices are independent of the mesh size as soon as P k and Q k are invertible. We will take k = k = k in what follows, and assume that k is not in the countable set of problematic frequencies (cf. Lemma 3). Otherwise, it would suffice to replace P k by −(ω∂ τ ) 2 − (k 2 + iη)ω 2 1/2 for some non-zero real number η to ensure that it is invertible for all k (and perform a similar modification for Q k ). It is not clear how to compute the required Galerkin matrices exactly, so we describe a way to approximate them numerically. Let us defineD 1 One can verify that ] ω . Note thatD 1 andD 2 are linear maps in a finite dimensional vector space, so they are be represented in the boundary element basis by square matrices MD 1 and MD 2 , explicitly given by The matrices MD 1 It remains to propose a method to compute the square root (or its inverse) of an invertible matrix with real eigenvalues. When k = 0, the eigenvalues are positive, and in this case we use the method of [22], relying on the discretization of contour integrals in the complex plane. When the frequency is non-zero, the previous method fails since the spectrum of the matrix contains negative values. We then follow Antoine and Darbas [4] by using a Padé approximation of the square root with regularization and rotation of the branch cut. We consider the classical Padé approximation where the coefficients c 0 , a i and b i are given by the formulas We use a "rotated" version of this approximation to avoid the singularity related to the branch cut for X < −1: This yields the new approximation This provides a good approximation of X → √ 1 + X in any region of the real line away from X = −1, as described by the next result.
Lemma 7 [31,Thm. 3.1] Let θ ∈ R, X ∈ R and r = |X + 1|. One has As a consequence, it is not difficult to check that, if θ ∈ (−π, π), then when N p → ∞, the Padé approximants converge exponentially in the uniform error in any region of the form δ < |X + 1| < R where 0 < δ < R . The previous scheme can be exploited to approximate D 0 − k 2 I d where D 0 is a matrix with real eigenvalues. Writing where I is the identity matrix. Using Lemma 7, one can conclude that, when N p is large enough, this yields a good approximation in the eigenspaces that are associated to eigenvalues λ of D 0 such that In our context, the eigenvalues λ ≈ k 2 correspond to the so-called "grazing modes".
To deal with them, in her thesis [18], Marion Darbas introduced a regularization recipe, which consists in adding a small imaginary part to the wavenumber. Namely, choosing some small ε > 0, the approximation is replaced by The following lemma gives sufficient conditions for the rhs in (19) to exist.

Lemma 8
Let D 0 be a symmetric non-negative matrix. Let θ ∈ (0, π) and ε ∈ (0, k). Then, for all j = 0..N p , the matrix is not an eigenvalue of D 0 . But under the conditions of the theorem, the following elementary computations reveal that the number z j has a non-zero imaginary part, while the eigenvalues of D 0 are real. The imaginary part of (k+iε) 2

B j vanishes if and only if
The two solutions of this equation are given by Remark that Im(B j ) has the same sign as sin θ (b 2 j − b j ) which is negative, while Re(B j ) has the same sign as (1 − cos θ)b 2 j + cos θ b j which is positive. Thus, u 1 ≤ 0 while u 2 = |B j |+|Re(B j )| |Im(B j )| ≥ 1. Therefore, the only two possible values of ε for which z j is real are outside the interval (0, k).
where D 0 = MD 2 + c 2 k 2 I . Then in both cases we apply formula (19).

Remark 3
The use of Padé approximation for Dirichlet-to-Neumann maps is common in wave problems [4,35]. Here we do not have to impose artificial boundary conditions on the end-points of the screen, an approach that has been adopted recently in the literature [35] for the approximation of a DtN map on parts of a polygonal boundary.

Remark 4
Although we have shown that the preconditioners are well-defined, it would also be desirable to guarantee some spectral equivalence between the approximate and the exact Galerkin matrices of P k and Q k . In particular, this would ensure that the preconditioners that we construct are invertible. We have to leave this question open for the moment.
In our implementation, we use the parameters N p = 15, θ = π 3 and ε = 0.05k 1/3 . Our numerical results do not depend crucially on those choices, see Sect. 6.4. The choice of N p is dictated by the ratio N dof k where N dof is the dimension of the Galerkin space. In our tests, this ratio will be held fixed unless stated otherwise, and thus, a fixed value of N p yields good results for all k. An informal way to explain this fact is that for our choices of operator X , the largest eigenvalues of the Galerkin matrix behave as λ M ≤ C 1 N 2 dof . Using a fixed number of points per wavelength gives in turn N dof = C 2 k and thus Choosing N p such that (18) holds with R = C 1 C 2 , and forgetting about the grazing modes, we thus ensure a small error in Eq. (17) independently of k. On the other hand, when k is fixed and the mesh is refined, it is necessary to increase N p to maintain accuracy.

Numerical results
In this section, we present some numerical results concerning the efficiency of the preconditioners defined above for solving the linear systems arising from the Galerkin method detailed in Sect. 4. The linear systems are solved with the GMRES method [43] with no restart (restarting GMRES does not improve the performance in our implementation) and a relative residual tolerance of ε = 10 −8 . Execution times are reported to show in which case the preconditioned linear system is solved faster than the non-preconditioned system. We report the best timing of three successive runs when the time is less than 5 seconds. If the number of iterations is greater than 500, we stop the calculations and report the time to reach the 500th iteration. The computations are also interrupted when they last more than 15 minutes and in this case the iteration number reached after 15 minutes is reported. All the simulations are run on a personal laptop with an eight-core intel i7 processor, a clock rate of 2.8GHz, and 16GB of RAM. The method is implemented in the language Matlab R2018a, and uses the environment Gypsilab 1 developed by Matthieu Aussal and François Alouges [3]. The full code of our method is freely available online 2 . In what follows, k stands for the wave number, N for the number of mesh points and |Γ | for the length of a curve Γ . Moreover, ω Γ is the weight defined on any curve Γ by where Z 1 and Z 2 are the end-points of the curve. This gives ω Γ (X ) = 1 − X 2 1 in the particular case Γ = [−1, 1]. For problems such that N > 5 × 10 3 , we compress the matrices using the Efficient Bessel Decomposition [7], in association with the "fiNUFFT" code from the Flat Iron institute [10,11].
As suggested by one of our referees, it is possible to use the mass matrix to precondition the Galerkin matrices. This does improve significantly the convergence rate of GRMES (probably owing to the non-uniformity of the mesh). However, this is far from begin competitive, both in time and in number of iterations, with the preconditioners analyzed here. For the sake of conciseness, numerical results for mass matrix preconditioning are not included.

Laplace equation on the segment
We start by testing the performance of the exact inverses, characterized in Sect. 2, as preconditioners. This is rather a validation stage. Segment, Laplace-Dirichlet problem. In Table 1, we report the timings and number of GMRES iterations for the Laplace weighted single-layer equation Two cases are considered, first without any preconditioner, and then with a preconditioner given by the exact inverse 2 −(ω∂ x ) 2 + 2 ln(2) π 0 (see Theorem 1 and the previous section for the detailed construction of the preconditioner). The rhs is chosen as A graph of the history of the GMRES relative residual is given in Fig. 3 for a mesh with N = 2000 elements.  we also report in Table 2 the timings and number of iterations of the GMRES method first without preconditioner, and then with the preconditioner obtained from the operator [−(∂ x ω) 2 ] −1/2 by the method described in the previous section. The rhs is chosen as A graph of the history of the GMRES relative residual is given in Fig. 3 for a mesh with N = 2000 elements. We observe in both Dirichlet and Neumann cases a low and stable number of iterations which is the expected behavior. In the Neumann problem, the presence of the preconditioner leads to huge speedups.

Helmholtz equation on the segment
We now turn our attention to the Helmholtz equation (k > 0). From now on, unless stated otherwise, the number of segments in the discretization is set to N ≈ 5k |Γ | (5 elements per wavelength). Segment, Helmholtz-Dirichlet problem. In Table 3 we report the number of GMRES iterations for the numerical resolution of the weighted single-layer integral equation , when using a preconditioner based on the opretor as compared to the case where no preconditioner is used. We take, for the Dirichlet data, the plane wave u D (x) = e ikx . We also provide, in Fig. 4a, the history of the GMRES relative residual with and without preconditioner, for a problem with k |Γ | = 200π . When the preconditioner is used, only a small increase of the number of iterations occurs as k increases. Segment, Helmholtz-Neumann problem. We run the same numerical comparisons, this time for the Helmholtz weighted hypersingular integral equation and taking the preconditioner based on the operator Results are given in Table 4 for different meshes and in Fig. 4b    Huge differences, both in time and number of iterations are shown in favor of the preconditioned system.

Helmholtz equation on non-flat arcs
Spiral-shaped arc We first consider a spiral-shaped arc of equation 1]. The curve has a length |Γ | ≈ 3.87. We report in Tables 5 and 6 the number of iterations and computing times respectively for the Dirichlet and Neumann problems with the rhs given by where for all (x, y), u inc (x, y) = e ikx . The results show that the preconditioning strategy is also efficient in presence of non-zero curvature. To illustrate the problem,  the scattering pattern with Neumann boundary conditions is illustrated in Fig. 5 for this geometry. Non-smooth arc We consider now the case where Γ is the V-shaped arc given by the parametric equations where θ is a parameter. When θ = π , Γ is the segment. When θ = π , this arc has a corner in the middle of angle θ . Since it is not a smooth arc, the theory presented in this work does not apply. For example, the solution α to the weighted single-layer integral equation where u D is a smooth function, has a singularity at the corner of Γ . This is illustrated in Fig. 6 where we plot the solution α(s) as a function of the arclength s for a rhs u D = e ikx , when θ = π 2 . Despite this singularity, the number of GMRES iterations remains independent of the mesh size for a fixed frequency. This result is reported in Table 7, where we compare, for k |Γ | = 50 fixed, the number of GMRES iterations for the resolution of the Helmholtz weighted single-layer integral equation respectively on the segment and on the V-shaped curve, for different values of N .
We show in Table 8 the influence of the frequency on the preconditioning performances for the Helmholtz weighted single-layer integral equation on the V-shaped arc of angle θ = π 2 . The results are qualitatively the same as in the case of a smooth curve. To illustrate the problem, the scattering pattern for a sound-hard V-shaped arc of angle π 2 (Neumann conditions) is shown in Fig. 7 for k |Γ | = 50. for a V-shaped sound-hard (Neumann boundary conditions) screen with θ = π 2 and k |Γ | = 50π . On the left, notice that the energy is deflected in the orthogonal directions. On the right, notice the resonating aspect of the solution

Influence of the number of Padé approximants
The method is not very sensitive to the number of Padé approximantes, i.e. the parameter N p in formula (19). We show this in Fig. 8 in the case of the Dirichlet problem for the spiral-shaped arc with k |Γ | = 200π and u D = e ikx . The parameter ε and the angle of the branch rotation θ remain fixed (see Sect. 5).

Importance of the correction
It is crucial to include the correct dependence in k in the preconditioners. We report here some numerical results in several situations where this dependence is not respected. Laplace preconditioning First, we precondition the Helmholtz weighted single-layer integral equation on the segment with the operator An identity is added under the square root to P 0 to make it invertible. It is easy to check that 2P 0 is spectrally equivalent to the inverse of S 0,ω with Since the Laplace preconditioner is also a compact perturbation of the inverse of S k,ω , the usual theory predicts that the number of iterations remains bounded when the frequency is fixed and the mesh is refined. This result is confirmed numerically in Fig.  9. We see indeed in practice that no matter how the mesh is refined, the number of iterations remains constant. However, we see in the previous example that this number of iterations is not independent of k and in fact grows with k. Including the dependence in k in the preconditioner reduces a lot this behavior, as illustrated in Fig. 10.
No singularity correction Second, we test the preconditioner without singularity correction, which is the method obtained when we take ω ≡ 1. That is, we solve the Fig. 9 History of the GMRES relative residual for the resolution of the Helmholtz weighted single-layer integral equation on the segment for k = 15 (left) and k = 50 (right). The mesh is progressibely refined, the level of refinement begin indicated by the color of the curves. We compare the situation where the linear system is preconditioned (circles) as opposed to the case where no preconditioner is used (crosses). The preconditioner is based on the operator P 0 = −(ω∂ x ) 2 + I d . Notice that the curves with circles are almost superimposed. We thus verify in practice that the number of iterations for the preconditioned system is independent of the discretization (though not of k) Fig. 10 History of the GMRES relative residual for solving the Helmholtz weighted single-layer integral equation preconditioned either with P 0 (circles) or P k (stars). Different colors represent different wave numbers. In each case, we keep the proportionality N ≈ 5k |Γ | non-weighted integral equation on the segment with a standard P 1 Galerkin method, and build a preconditioner based on the operator This is the direct application of the method of Antoine and Darbas [4] to the context of an open curve. As stated at the beginning of Sect. 4, if a uniform mesh is used, the Galerkin approximation converges at the rate O( √ h) only. To better capture the Table 9 Number of GMRES iterations and H −1/2 error for the resolution of the singe-layer integral equation S k λ = u D on the segment, with k = 10π , N = 80. In the first five lines, the equation is solved using a standard Galerkin method on a graded mesh with parameter β ranging from 1 (uniform mesh) to 5. In all those 5 cases, we precondition the linear system using a discrete version of P k Eq. (20). The last line refers to the weighted Galerkin method described in this work, with the new square root preconditioner. Notice that even for this small problem, the Galerkin matrices associated to graded meshes become very close to singular, and the preconditioners are no longer accurately evaluated using the Padé approximation as in Sect. 5. The square root is instead computed directly using an eigenvalue decomposition of the matrices.
To compute a good reference solution this problem, we had to use a 20 times finer mesh, preventing us to do the same test for a larger wave number singularity, one simple idea is to use a mesh graded towards the end-points. A graded mesh of parameter β is a mesh such that near the edge, the width of the i-th interval is approximately (ih) β . The parameter β = 2 corresponds to the mesh in our Galerkin method defined in Sect. 4. The parameter β = 5 is the one that theoretically leads to the same rate of convergence as in our method [38]. In Table 9, we report the number of iterations in the GMRES method for the resolution of the Helmholtz (non-weighted) single-layer integral equation preconditioned by P k for k = 10π and several graded meshes corresponding to different choices of β, ranging from 1 (uniform mesh) to 5. We also report the number of iterations for the resolution of the weighted Galerkin method with the weighted square-root preconditioners. In each case, we indicate the H − 1 2 error. One can see that mesh-refinement allows to decrease the error at the price of losing the performance of the preconditioner. This highlights the need for the method introduced in this paper.

Comparison with the generalized Calderón preconditioners
We finally adapt to our context the idea of Bruno and Lintner [13], namely to use S k,ω and N k,ω as mutual preconditioners. Notice that the way we discretize the problem is different from [13] where a spectral method is used. In our setting, using the notation of Sect. 5, we define the preconditioners respectively for the weighted single-layer and weighted hypersingular integral equations. We report the number of iterations and computing times respectively for the Dirichlet and Neumann problems on the segment respectively in Table 10 and Table 11. The performance is compared to that of our new preconditioners. The rhs are respectively u D = u inc and u N = ∂u inc ∂n where u inc is a plan wave of angle of incidence π 4 . Our results confirm the efficiency of the Generalized Calderón preconditioners, whose iteration counts remain very stable with respect to k. Despite a slightly larger increase of the iterations for our method, the resolution remains faster for the tests presented here, particularly for the Dirichlet problem. This is due to the fact that our preconditioners are evaluated faster.

Conclusion
We have presented a new approach for the preconditioning of integral equations coming from the discretization of wave scattering problems in 2D by open arcs. The methodology is very effective and proven to be optimal for Laplace problems on straight segments. It generalizes the formulas mainly proposed in [4] for regular domains. They are only modified by a suitable weight. We deeply believe that the methodology opens new perspectives for such problems. First, it is possible to generalize the approach in 3D for the diffraction by a disk (see [8,Chap. 4]). Second, the strategy that we used here can probably be extended to the half line and, hopefully, to 2D sectors, giving, on the one hand a new pseudo-differential analysis more suitable than classical ones (see e.g. [34,41,42]) for handling Helmholtz-like problems on singular domains, and, on the other hand, a completely new preconditioning technique adapted to the treatment of BEM operators on domains with corners or wedges in 3D. Eventually, the weighted square root operators that appeared in the present context might well be generalized to give suitable approximation of the exterior Dirichlet to Neumann map for the Helmholtz equation which is of particular importance in e.g. domain decomposition methods. Having such approximations might therefore lead to better methods in that context too.