Convergence of parallel overlapping domain decomposition methods for the Helmholtz equation

We analyse parallel overlapping Schwarz domain decomposition methods for the Helmholtz equation, where the exchange of information between subdomains is achieved using first-order absorbing (impedance) transmission conditions, together with a partition of unity. We provide a novel analysis of this method at the PDE level (without discretization). First, we formulate the method as a fixed point iteration, and show (in dimensions 1, 2, 3) that it is well-defined in a tensor product of appropriate local function spaces, each with L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document} impedance boundary data. We then obtain a bound on the norm of the fixed point operator in terms of the local norms of certain impedance-to-impedance maps arising from local interactions between subdomains. These bounds provide conditions under which (some power of) the fixed point operator is a contraction. In 2-d, for rectangular domains and strip-wise domain decompositions (with each subdomain only overlapping its immediate neighbours), we present two techniques for verifying the assumptions on the impedance-to-impedance maps that ensure power contractivity of the fixed point operator. The first is through semiclassical analysis, which gives rigorous estimates valid as the frequency tends to infinity. At least for a model case with two subdomains, these results verify the required assumptions for sufficiently large overlap. For more realistic domain decompositions, we directly compute the norms of the impedance-to-impedance maps by solving certain canonical (local) eigenvalue problems. We give numerical experiments that illustrate the theory. These also show that the iterative method remains convergent and/or provides a good preconditioner in cases not covered by the theory, including for general domain decompositions, such as those obtained via automatic graph-partitioning software.

boundary condition, geometry and variable k, our theory is restricted here to the homogeneous interior impedance problem with k constant, and the boundary condition ∂u ∂n − iku = g on ∂Ω, (1.2) where ∂u/∂n is the normal derivative, outward from Ω.

Parallel domain decomposition method
To solve (1.1), (1.2), we use a parallel overlapping Schwarz method with impedance transmission condition, based on a set of Lipschitz polyhedral subdomains {Ω j } N j=1 , forming an overlapping cover of Ω. To derive this, note that if u solves (1.1), (1.2), then, the restriction of u to Ω j : where ∂/∂n j denotes the outward normal derivative on ∂Ω j . To iteratively solve (1.4) -(1.6), we introduce a partition of unity (POU), {χ j } N j=1 , with the properties for each j : suppχ j ⊂ Ω j , 0 ≤ χ j (x) ≤ 1 when x ∈ Ω j , and j χ j (x) = 1 for all x ∈ Ω.
   (1.7) Then, the parallel Schwarz method reads as follows: given an iterate u n defined on Ω, we solve each local problem on Ω j for u n+1 j , (∆ + k 2 )u n+1 ∂ ∂n j − ik u n+1 j = g on ∂Ω j ∩ ∂Ω, (1.10) and the new iterate is the weighted sum of the local solutions: This well-known method can be thought of as a generalization of the classical algorithm of Després [16], [4] to the case of overlapping subdomains. The form of the algorithm stated above can be found in [18, §2.3]. The novel contribution of this paper is the convergence analysis of the method.

The main results and structure of the paper
The main results of this paper are as follows. 1. A proof that the iterative method (1.8)-(1.11) is well-defined for general Lipschitz subdomains (Theorem 2.12) using well-posedness results for the Helmholtz equation on Lipschitz domains and the harmonic-analysis results of [42].
2. The formulation of the fixed-point iteration (3.9) for the error vector e n , where e n j = u j −u n j , j = 1, . . . , N , and the expression of powers of the fixed-point operator T in terms of "impedance-toimpedance maps" linking pairs of subdomains with non-trivial intersection (Theorem 3.9).
3. For 2-d rectangular domains covered by overlapping strips, with each subdomain only overlapping its immediate neighbours, sufficient conditions for T N to be a contraction, where N is the number of subdomains (Theorem 4.13 and its corollaries). These conditions are formulated in terms of norms of impedance-to-impedance maps and compositions of such maps. 4. A summary and explanation of the results from the companion paper [44] that bound impedanceto-impedance maps using rigorous high-frequency asymptotic analysis (a.k.a., semiclassical analysis). In particular, these results indicate that Theorem 4.13 implies contractivity of T N for a model case with two subdomains and provided that the overlap is sufficiently large (see §4.4.4).

5.
Finite element experiments ( §6) that both back up the theory and investigate scenarios out of the theory's current reach. Since the theory in Points 1-4 is at the continuous level without discretization, §5 first gives a description of the finite element algorithms used, along with justification that the results in §6 are reliable.
The experiments related to the theory illustrate how the good behaviour of the relevant impedanceto-impedance maps induces good convergence of the iterative method. Those beyond the theory show, for square domains and general domain decompositions, that the fixed point operator still enjoys a power contractivity property ( §6).
Regarding 2 and 3: to our knowledge this is the first time that overlapping DD methods for Helmholtz have been analysed in terms of impedance-to-impedance maps. This analysis therefore gives a route to analyse overlapping DD methods for Helmholtz using the PDE theory of the Helmholtz equation, which will then be the focus of [44]. Interest in impedance-to-impedance (a.k.a., Robin-to-Robin) maps can be widely found in the non-overlapping literature -see the references in the literature review below. These maps also arise in the formulation of fast direct methods (e.g. [34], [55]) and the recent work [3] analyses these maps in this setting (using complementary techniques to those in [44]). Previous work of some of the authors (e.g., [38,36]) also used PDE theory to analyse overlapping DD preconditioners; while this work was able to cover very general geometries, it was limited to the case when k > 0, corresponding to media with some absorptive properties. In the present paper we consider only the propagative case k = 0.
Regarding 3: In 1-d we recover the classical result (a special case of [54]) that, with N subdomains, the N th power of the fixed point operator is zero, and so the iterative method converges in N iterations.
The structure of the paper is as follows. §2 contains the necessary well-posedness and regularity theory for the Helmholtz equation needed to prove the results described in Point 1 above. §3 formulates the fixed-point iteration as described in Point 2. §4 focuses on 2-d strip decompositions as described in Point 3. §5 describes the set-up for the finite-element computations used to illustrate the results, with the results of these computations given in §6.

Literature review
In the last decade there has been an explosion of progress in the construction and analysis of solvers for frequency-domain wave problems. Here we highlight those parts of the literature most related to the present paper; more substantial recent reviews can be found, for example in [33] and in the introductions to [38] and [60].
The method (1.8)-(1.11) can be thought of as a straightforward extension of the parallel nonoverlapping method of Després [16], [4] to the case of overlapping subdomains. In [16], [4] the coupling between subdomains at each iteration is achieved by feeding to each subdomain impedance data from its neighbours at the previous iteration. In the overlapping case considered here, the boundary impedance data for the next iterate on a given subdomain is a weighted average of data coming from all subdomains overlapping it (see (1.9) and (1.11)).
here has the more modest aim of proving power contractivity of the fixed point operator. This turns out to be not only provable for a model problem in simple-enough geometry but also observable in more general situations, giving hope that the present theory can be generalised. Power contractivity of the fixed point operator also ensures convergence rates for preconditioned GMRES.
The estimates ensuring power contractivity in §4, involving the norms of impedance-to-impedance maps (which can be computed by solving canonical eigenvalue problems) are in some sense analogous to condition number estimates in the positive definite case: both estimates provide upper bounds that can be controlled in certain parameter ranges, but the actual value of the bound is rarely computed when solving a particular problem (so the bounds are "descriptive" rather than "prescriptive").

1.6
The discrete analogue of the results of this paper In [35] it was shown (see also §5.1) that a natural finite element counterpart of the parallel iterative method considered in the current paper is the Restricted Additive Schwarz method with impedance transmission condition (often called the Optimized Restricted Additive Schwarz (ORAS) method).
Since this paper was written, three of the current authors developed a convergence theory for ORAS, thus providing a discrete version of the theory given here. These results are presented in [37].
The theory in [37] applies to discrete Helmholtz systems arising from conforming nodal finite elements of any polynomial order and a general theory for discrete fixed point iteration analogous to § §2, 3 on general Lipschitz domains and partitions of unity is presented. For domain decompositions in strips in 2-d, we show that, when the mesh size is small enough, ORAS inherits the convergence properties of the parallel iterative method at the PDE level which are proved here, independent of the polynomial order of the elements. The proof relies on characterising the ORAS iteration in terms of discrete 'impedance-to-impedance maps' on local discrete Helmholtz-harmonic spaces, which we prove (via a novel weighted finite-element error analysis) converge as h → 0 in the operator norm to their non-discrete counterparts (i.e., the operators analysed here). This discrete theory thus justifies the use of the finite element method to illustrate the properties of the iterative method at the PDE level, as we have done in §6 of this paper.
Sesquilinear forms. We define the global and local sesquilinear forms by Prolongation and restriction. We build a prolongation R : H 1 (Ω ) → H 1 (Ω) by multiplication by the POU, i.e., for each v ∈ H 1 (Ω ), Recalling that u denotes the restriction of u to Ω (see (1.3)), we have the important property The main purpose of this section is to justify step (1.9) in the domain decomposition algorithm, by showing that for each n, the impedance trace of u n is in L 2 (∂Ω ), for each . This then ensures that, for each , u n+1 is well-defined in the space U (Ω ) defined below and hence u n+1 ∈ U (Ω), so that u n+1 in turn provides suitable L 2 impedance traces on each ∂Ω for the next iteration. The main result of this section is Theorem 2.12. To prove it we need to analyse (1.8)-(1.11) in a space of higher regularity than H 1 . In what follows we need the following further property of the partition of unity {χ }.
Such a partition of unity exists by, e.g., [48,Theorem 3.21 and Corollary 3.22]. We note for later that Assumption 2.1 implies that ∂χ /∂n = 0 on ∂Ω \ ∂Ω, and thus, for any v ∈ H 1 (Ω ), Notation 2.2. Where possible, we explicitly indicate the dependence of our estimates on the wavenumber k. In this context, we always assume k ≥ k 0 where k 0 > 0 is chosen a priori and we use the notation A B to mean A ≤ CB with a constant C independent of k (but possibly depending on k 0 ) and A ∼ B to mean A B and B A.

The Helmholtz problem in spaces U (D) and U 0 (D)
In this subsection D denotes a general Lipschitz domain, with boundary ∂D.
Since the trace operator maps H 1 (D) to H 1/2 (∂D) ⊂ L 2 (∂D), an equivalent definition of U (D) is in the rest of the paper we refer to U 0 (D) as the space of Helmholtz-harmonic functions on D.
Lemma 2.4 (Well-posedness of the Interior Impedance Problem in U (D)). Suppose D is Lipschitz and consider the problem ∆u + k 2 u = −f in D and ∂u/∂n − iku = g on ∂D, (2.8) with f ∈ L 2 (D) and g ∈ L 2 (∂D). Then there exists a unique solution u ∈ U (D) and C j = C j (k), j = 1, 2, such that Proof of Lemma 2.4. By the standard result about well-posedness of the interior impedance problem for Lipschitz D (see, e.g., [57, § §6.1.3, 6.1.6]), a unique solution u exists and there exist C j = C j (k), j = 1, 2, such that Without loss of generality we can assume that C j (k) 1, for j = 1, 2. Then, multiplying the PDE in (2.8) by u and using Green's identity (see, e.g., [48,Lemma 4.3]), we find that Inserting the boundary condition from (2.8), taking the imaginary part, and using the Cauchy-Schwarz inequality gives k u 2 Now, multiplying (2.11) through by k and using the inequality 2ab ≤ εa 2 + ε −1 b 2 , for any a, b, ε > 0 to estimate both terms on the right-hand side, we have Combining this with (2.10), we obtain which is in the required form (2.9).
To complete the bound on u U (D) , we therefore only need to bound k −1 ∆u L 2 (D) and ∂u/∂n L 2 (∂D) by the right-hand side of (2.9); a bound on k −1 ∆u L 2 (D) follows from the PDE in (2.8) together with (2.12). The required bound on ∂u/∂n L 2 (∂D) follows from the boundary condition in (2.8) together with (2.12). If D is either Lipschitz star-shaped or C ∞ , then (2.9) holds with C 1 (k) ∼ 1 and C 2 (k) ∼ 1 by [52, Equation 3.12] and [2, Theorem 1.8 and Corollary 1.9] respectively.
We now use results from the harmonic-analysis literature about the Laplacian on Lipschitz domains to give an alternative characterisation of the space U (D). Let The following theorem is a consequence of [42,Corollary 5.7]; see [15,Lemma 2].
Since the H 3/2 (D; ∆) norm is characterised only through norms on D (as opposed to the norm on U (D), which is characterised through norms on both D and ∂D), the following corollary holds. By Theorem 2.6 and the definition of U (D) in §2.2, the norms · H 3/2 (D;∆) and ||| · ||| U (D) defined by are equivalent. Moreover the equivalence constants are k−independent (since k does not feature in the definition of either norm.). We therefore have the following corollary.
(i.e., the omitted constant is independent of k).
The following lemma studies the behaviour of the impedance trace of a function u ∈ U (D) on an interface interior to D. This plays a key role in the analysis of the iterative method (1.8)-(1.11).
i.e., the impedance-to-impedance map (defined more precisely in §3.2) is bounded as an operator from L 2 (∂D) to L 2 (∂D ).
We make two remarks: (i) Lemma 2.10 makes no assumptions about the geometries of D and D, other than that they are both Lipschitz and D ⊂ D.
(ii) The powers of k in (2.19) and (2.20) are almost-certainly not sharp; this is because the righthand side of the trace result (2.15) involves a norm on H 3/2 (D; ∆) that does not weight derivatives with the appropriate powers of k (in contrast to, e.g., (2.7)). As far as we are aware, the result analogous to (2.15) with an H 3/2 (D; ∆) norm weighted with k and a potentially-k-dependent constant has not yet been proved for general Lipschitz domains.

Framework for the convergence analysis
In this section we develop the tools needed to analyse the algorithm (1.8)-(1.11) in the space U 0 .

The error propagation operator T
This motivates the introduction of the operator-valued matrix T = (T j, ) N j, =1 , defined as follows. For v ∈ U (Ω ), and any j, Therefore, e n+1 j = T j, e n , and so e n+1 = T e n . (3.9) (ii) It is convenient here to introduce the notation so that (3.7) holds on Γ j, and (3.8) holds on ∂Ω j \Γ j, .
Since e n j is Helmholtz-harmonic in Ω j for each j, we aim to analyse convergence of (3.9) in the space U 0 defined in (2.22). For the rest of this section we restrict to 2-d and 3-d, using the following norm, previously introduced in [4, Equation 12]. (The 1-d case is discussed brielfy in §4.3, where the norm on the boundary data is trivially the modulus on C.) where ∂/∂n denotes the outward normal derivative on ∂D.

12)
and so · 1,k,∂D is a norm on U 0 (D). Furthermore, if D is either Lipschitz star-shaped or C ∞ , then · 1,k,∂D is equivalent to · U (D) , with equivalence constants independent of k.
Proof. If v ∈ U 0 (D), then by Green's first identity (see, e.g., [48,Lemma 4.3]), Taking the imaginary part, we have Thus, To obtain the norm equivalence, observe that, for v ∈ U 0 (D), Moreover since ∂v/∂n and v both belong to L 2 (∂D), Lemma 2.4 implies that The stated k−independence follows from Remark 2.5.
Using (3.11), we define the norm on U 0 : For simplicity, we now assume that each Ω is star-shaped Lipschitz, so that the norm equivalence in Lemma 3.3 holds with constants independent of k. Analogues of the following results for general Lipschitz Ω hold, but with different k-dependence.
Furthermore, to simplify the notation, we define the operator The next theorem summarises some basic properties of the operator T j, on the space U 0 (Ω ). Also,

16)
and T j, which, recalling (3.10), vanishes on ∂Ω j \Γ j, . Differentiating the product on the right-hand side of (3.17) yields (3.15). Then, taking norms of both sides of (3.15) and using Assumption 2.1 and the fact that 0 ≤ χ ≤ 1, we obtain Then, using the fact that Γ j, ⊂ Ω is an interface in Ω , using the multiplicative trace theorem and then Lemma 3.3, we obtain Combining this with (3.18) yields (3.16). Finally, to obtain the boundedness of T j, : U 0 (Ω ) → U 0 (Ω j ), we use Lemma 2.10 (ii), to obtain and we then combine this with (3.16).
Remark 3.6. The same arguments show that T j, : In the following section we are interested in proving power contractivity of the error propagation operator T . This motivates us to study the composition T j, T ,j ; indeed, where the sum is over all ∈ {1, 2, . . . , N }\{j, j }, with Γ j, = ∅ = Γ ,j . A useful expression for the action of (3.20) can be obtained by inserting v = T ,j z j , with z j ∈ U (Ω j ), into (3.15), to obtain The first term on the right-hand side of (3.21) is of key interest in this paper. We note that its value is obtained by (i) finding T ,j z j , i.e., the unique function in U 0 (Ω ) with impedance data on Γ ,j given by imp (χ j z j ) ; (ii) evaluating imp j (T ,j z j ) on Γ j, and (iii) then multiplying the result by χ . Combining steps (i) and (ii) leads us to the following key definition.

24)
and Proof. Since T ,j z j ∈ U 0 (Ω ) and imp (T ,j z j ) vanishes on ∂Ω \Γ ,j , we have Substituting this in (3.21) gives (3.24). The estimate (3.25) is obtained by following the proof of (3.16), with v = T ,j z j .
We now see from (3.25) that (at least for sufficiently large k and/or sufficiently small ∇χ L ∞ (Γ j, ) ) the right-hand side of (3.24) is dominated by the first term. The norm of the impedance-to-impedance map lies at the heart of the convergence theory in §4.

Convergence of the iterative method for strip decompositions
In this section we obtain a convergence theory for the iterative method (1.8) -(1.11) when the domain Ω is either an interval in 1-d or a rectangle in 2-d. In 1-d the subdomains are intervals and in 2-d the subdomains are sub-rectangles.    To illustrate these definitions, given g ∈ L 2 (Γ − ), let u denote the Helmholtz-harmonic function on Ω with left-facing impedance data g on Γ − and zero impedance data, elsewhere. Then I Γ − →Γ + −1 g is the right-facing impedance data of u on Γ + −1 and I Γ − →Γ − +1 g is the left-facing impedance data of u on Γ − +1 . Note that Γ + −1 and Γ − +1 are both interior interfaces in Ω (see Figure 4.2). Recall from Remark 3.1 that T , = 0 and T j, = 0 if Ω j ∩ Ω = ∅. Therefore, T takes the block tridiagonal form

Notation common to both 1-d and 2-d
where L and U are the lower and upper triangular components of T . We record for later that, by the Cayley-Hamilton theorem, In what follows a crucial role is played by the products: We remark that in 2-d the structure (4.3) remains the same if the vertical interfaces in Γ ± are replaced by non-intersecting polygonal pieces; however, we do not pursue this generalisation here. The diagonal entries in (4.5) can be estimated in terms of impedance-to-impedance maps -see Theorem 3.9.

The results of §3 specialised to strip decompositions
Since strip decompositions have the property (4.2), Theorem 3.9 simplifies to the following.
and thus Proof. To obtain (4.6), without loss of generality, consider the case j = − 1 = j . Then, for any v ∈ U 0 (Ω ), In (4.8), the first equality comes from the definition of T −1, , the second comes from the fact that (by Assumption 2.1), (∂χ /∂n −1 ) = 0 on Γ + and the final equality comes from the fact (see (4.1) and Using this instead of (3.15) and propagating this simplification through the arguments using to prove Theorems 3.5 and 3.9 gives the result.

One dimension
The following result is known from [54, Propositions 2.5 and 2.6] (restricted to 1-d), but we state it here in our notation, because it helps motivate our approach in the 2-d case.

Two dimensions
In the rest of this section our goal is to estimate T n , where T = L + U is given by (4.3). In 1-d, T n takes the simple form given in Proposition 4.5, however this is not the case in 2-d. Our bounds in 2-d on T n are therefore based on the following elementary algebraic result.

An elementary algebraic result and its consequences
For integers n ≥ 1 and 0 ≤ j ≤ n − 1, let P(n, j) denote the set of monomials of order n in the two variables x, y that take the form p(x, y) = x s 0 y s 1 x s 2 . . . x s j (j even) or x s 0 y s 1 x s 2 . . . y s j (j odd) (4.17) or p(x, y) = y s 0 x s 1 y s 2 . . . x s j (j odd) or y s 0 x s 1 y s 2 . . . y s j (j even),   Proof. The formula (4.21) follows directly from Proposition 4.6. To obtain the final statement, note that P(n, 0) = {x n , y n }, so by (4.4), when n ≥ N , p(L, U) = 0, for p ∈ P(n, 0).  i.e., Γ δ is an interior interface; see Figure 4.3. Let u be the solution to Then define the canonical left-to-right and left-to-left impedance-to-impedance maps by

24)
and define the following norms of these maps ρ(k, δ, L) = sup (4.25) By Part (ii) of Lemma 2.10, ρ and γ are well-defined. We record the following simple relationship between γ and ρ.  Proof. Let u ∈ U 0 ( Ω) be the solution to (4.23). By Lemma 3.3 ∂ Using the boundary conditions in (4.23) together with (4.27), we obtain Now let Ω − denote the subdomain of Ω with height 1 and vertical sides Γ δ and Γ + (see Figure 4.3).
Since u ∈ U 0 (Ω − ), repeating the argument above gives Since ∂Ω − \Γ δ ⊂ ∂Ω, we can use the definition (4.25) of ρ(k, δ, L) to estimate the first term on the right-hand side of (4.29) and use (4.28) to estimate the second term on the right-hand side of (4.29): The result follows from the definition of γ(k, δ, L) (4.25).  Up to now we have developed the theory with general L and δ to emphasise that these can vary with . To reduce technicalities in the remainder of the theory, we introduce the simplifying notation.  We make this slight abuse of notation to avoid introducing additional symbols for the maxima above.
Proof. We use Theorem 4.7 with n = N , so the j = 0 term in (4.21) vanishes. We now claim that, for each p ∈ P(N, j) with j ∈ {1, . . . , N − 1}, and for any v ∈ U 0 , We prove (4.42) in the case j = 1, where p(x, y) = x s 0 y s 1 with s 0 + s 1 = N ; the case of higher j is obtained by induction. Then using Lemma 4.12 freely, Hence, combining (4.42) and (4.22), and an application of the Binomial Theorem gives (4.41).
Corollary 4.14 (Estimate of T N , useful for ρ small). Assume ρ ≤ ρ 0 ≤ γ and N ≥ 3. For any v ∈ U 0 , Thus, if ρ is small (relative to γ and N ), then T N is a contraction.
Proof of Corollary 4.14. By Theorem 4.13, Taylor's theorem, and the fact that ρ ≤ γ, where we have used the fact that the function x → (γ + x) N −3 is increasing on [0, ρ 0 ] to bound the Taylor-theorem remainder. Corollary 4.14 provides an estimate for T N 1,k,∂ that is first order in ρ. An estimate with a higher order in ρ can be obtained by considering higher powers of T .
Proof. The proof uses estimate (4.22) with n = sN . Consider any monomial p ∈ P(sN, j) with j ≤ s − 1. This is a monomial of order sN with j ≤ s − 1 transitions from x to y or from y to x. Thus it must contain at least one string of length ≥ N without jumps. (For example, any p ∈ P(2N, 1) must contain one string of length ≥ N without a jump.) Hence, using (4.4), p(L, U ) = 0 for all p ∈ P(sN, j) when j ≤ s − 1.
and thus the result follows in a similar way to that in the proof of Theorem 4.13.

Estimating the canonical map I −+ using semiclassical analysis
Theorem 4.13 and Corollaries 4.14 and 4.15 show that convergence of the iterative method improves as ρ gets smaller. We now describe results from [44] that give sharp bounds on the large-k limit of ρ (with other parameters, such as δ, fixed).
In the canonical domain ( For simplicity, [44] considers the case when the (outer) boundary condition (1.2) is replaced by a condition that the solution is "outgoing" (in a sense made precise by the notion of the wavefront set); i.e., that no outgoing rays hitting ∂ Ω are reflected. Studying this situation therefore focuses on the effect of the impedance boundary conditions coming from the domain decomposition, and ignores the effect of any high-frequency reflections from absorbing boundary conditions on ∂ Ω (see [23] for a precise description of these reflection effects). The outgoing condition replacing (1.2) is, in some sense, the ideal absorbing boundary condition at high frequency on ∂ Ω. Since perfect matched layers approximate the outgoing condition exponentially well at high frequency [24], we expect that the results of [44] will also hold when the outgoing condition is replaced by perfectly matched layers (and this is work in progress).
The paper [44] considers the following two model problems. Model Problem 1: the canonical problem specified in Definition 4.9 with outgoing conditions on the top and bottom, impedance data posed on the left, and zero impedance data on the right (i.e., that discussed above), and Model Problem 2: the canonical problem with outgoing conditions on the top, bottom, and right sides, and impedance data posed on the left.
Model Problem 2 is the canonical problem for the strip-decomposition algorithm applied with two subdomains when the global problem is (1.1) with outgoing boundary conditions. The reason for considering this further simplification is that in Model Problem 1 a ray moving from left to right can still be reflected an infinite number of times, and the reflection coefficient on Γ − depends on the data; thus an upper bound for general impedance data in this situation is more challenging to prove.
Upper and lower bounds on I −+ for Model Problem 2. Let Ω be the canonical domain of where the outgoing condition near Γ out is defined in terms of the wavefront set, as will be explained in [44]. In analogy with (4.24), I −+ : L 2 (Γ − ) → L 2 (Γ δ ) is defined by Then, for any > 0, there exists k 0 ( ) > 0 such that, for all k ≥ k 0 , Furthermore, for any 0 < < θ max , Observe that there exist C 1 , C 2 > 0 such that and thus Theorem 4.16 shows that lim k→∞ I −+ L 2 (Γ − )→L 2 (Γ δ ) is bounded above and below by multiples of (θ max ) 2 , and hence multiples of δ −2 , where we recall that δ is the distance of Γ δ from Γ − . The idea behind Theorem 4.16. The tools of semiclassical/microlocal analysis decompose solutions of PDEs in both frequency and space variables. These tools show that, at high-frequency, Helmholtz solutions propagate along the rays of geometric optics, in the sense that the wavefronts are perpendicular to the ray direction. The ideas behind Theorem 4.16 can therefore be understood by first looking at the impedance-to-impedance map for plane-wave solutions of the Helmholtz equation (since these are simple Helmholtz solutions travelling along rays), ignoring the fact that these do not satisfy the outgoing condition on all of Γ out , and thus are not solutions of Model Problem 2.
(∂ x − ik)u| Γ δ = ik(cos θ − 1) exp ik(δ cos θ + y sin θ) , so that, for this class of u, We now use (4.50) as a heuristic for the behaviour of the impedance-to-impedance map on solutions of Model Problem 2 travelling on rays at angle θ to the horizontal. Since the solution of Model Problem 2 is outgoing on Γ out , anything reaching Γ δ must arrive on a ray emanating from Γ − and hitting Γ δ , and the maximum angle such rays can have with the horizontal satisfies tan θ max = δ −1 ; see Figure  4.3. The right-hand side of (4.50) is largest when θ = θ max , with this expression then (modulo the presence of and ) the right-hand sides of (4.47) and (4.48). The arguments in [44] use these ideas in a rigorous way; for example, to prove the lower bound (4.48), we take a sequence of data (g(k)) k>0 where the Helmholtz solutions it creates are concentrated at high frequency in a beam coming from one point of Γ − and traveling in one direction (cos θ, sin θ), and we take θ to be arbitrarily close to θ max . The notion of concentration at high frequency is understood in a rigorous way using so-called semiclassical defect measures; see [45, §9.1] for an informal overview of these, and [64,Chapter 5], [9], [50], [25], [24].
Finally, we highlight that these ray arguments and angle considerations are similar to those in [26, §5] used to optimise boundary conditions in domain decomposition for the wave equation.

Estimating higher order products of L and U
The estimates in Theorem 4.13 and Corollaries 4.14 and 4.15 use Lemma 4.12 repeatedly to bound T n 1,k,∂ in terms of powers of ρ and γ. For example, to bound the term L s U for an integer s > 0, the argument in Theorem 4.13 uses (4.36)-(4.38) to obtain Thus if ρ is controllably small, Corollary 4.14 implies power contractivity for T . The use of Corollary 4.14 is illustrated in Experiment 6.1 below, which shows that the convergence rate of the domain decomposition method improves as ρ decreases. However, we expect that estimates like that in Corollary 4.14 are not in general sharp. In particular, looking at the case k = 80 in Figure 6.2 and Table 3 of §6 we see a case when ρ ≈ 0.15, but the method converges effectively for N = 4, 8, 16, even though (4.43) grows linearly in N . Thus, we expect that sharper results may be possible by bounding composite maps such as L s U directly, rather than estimating each of their factors, as in (4.51). In fact, in [44], ray arguments are used to give insight into the behaviour of these composite maps in the k → ∞ limit, and these arguments do indeed indicate that the compositions of the maps behave better than the products of the norms of the individual components.
To illustrate the use of composite maps, we consider the dominant (j = 1) term in (4.22): One of the terms appearing inside the maximum corresponds to p(L, U ) = L N −1 U . This operator is blockwise very sparse; for N ≥ 2 all its nonzero blocks lie on the (N − 2)th diagonal below the main diagonal (see (4.5) for the case N = 2). The (N, 2)th element of L N −1 U is where the operator product is understood as concatenated on the left as the counting index j increases. Rewriting (4.6) using the notation (4.1), we see that, for any s, A straightforward induction argument then shows that In Experiment 6.3 we use (4.54) to compute the norm of the composite operator L N −1 U directly.

Finite-element approximations
In this section we describe how we use finite-element computations to illustrate our theoretical results. Due to space considerations, we restrict here to a description of algorithms and brief remarks on finite-element convergence; more details are in [37]. For any domain Ω, let T h be a family of shape-regular meshes on Ω with mesh diameter h → 0. We assume each mesh resolves the boundaries of all subdomains. Let V h be an H 1 -conforming nodal finite-element space of polynomial degree p defined with respect to T h . For any subset (domain or surface) Λ that is resolved by T h , we define V h (Λ) = {w h | Λ : w h ∈ V h }. We let N (Λ) denote the set of nodes of the space V h that lie in Λ.

The iterative method
Here we describe the computation of finite-element approximations of the iterates u n defined in (1.8) -(1.11). With a as in (2.2), and for any F ∈ H 1 (Ω) , we consider finding u ∈ H 1 (Ω) satisfying a(u, v) = F (v) for all v ∈ H 1 (Ω); (5.1) this includes the weak form of (1.1), (1.2) as a special case. To discretize (5.1), we define To formulate the discrete version of (1.8) -(1.11) on each subdomain Ω , we introduce the local space V h := V h (Ω ), and define the local operators A h, : Note that the extension R h, v h, ∈ V h is defined nodewise: it coincides with v h, at nodes in Ω and vanishes at nodes in Ω\Ω . Thus R h, v h, ∈ H 1 (Ω) is a finite-element approximation of the operator of extension by zero, even though the (true) extension by zero does not, in general, map H 1 (Ω ) to H 1 (Ω). We define the restriction operator R h, : It is shown in [35] that a natural discrete analogue of (1.8)-(1.11) is The algorithm (5.4), (5.5) is derived in [35] as a finite-element approximation of (1.8)-(1.11). In fact (5.4), (5.5) is equivalent to the well-known Restricted Additive Schwarz method with impedance transmission condition (also known as WRAS-H [43] and ORAS [18,Definition 2.4] and [58]).
Moreover, since u h is the exact solution of (5.1), we have, trivially, The error is then e n h := (e n h,1 , e n h,2 , · · · , e n h,N ) , where e n h, = u h | Ω − u n h, . Subtracting (5.4) from (5.6), we obtain the error equation The two expressions in (5.7) can be combined and written in the operator matrix form: providing a finite element analogue of (3.9). The matrix form of T h is discussed in [35, §5].
In §6 we plot error histories for this method. To do this, we need to choose a suitable norm in which to measure the error. Since e n h, ≈ e n ∈ U 0 (Ω ) (defined in Definition 2.3), it is natural to try to analyse e n h, in a finite-element analogue of U 0 (Ω ). In fact, one can show that, for each n, e n h, ∈ V h ,0 := {v h, ∈ V h : a (v h, , w h, ) = 0 for any w h, ∈ V h ∩ H 1 0 (Ω )}, which indicates that the error is 'discrete Helmholtz harmonic'. Therefore we define the norm: This is a norm for h sufficiently small because the sesquilinear form a satisfies a discrete inf-sup condition on V h (Ω ) × V h (Ω ) (with h-independent constant) [49,Theorem 4.2]. The norm of the error vector e n h is then given by

The impedance-to-impedance maps
We now describe the computation of the canonical impedance-to-impedance maps I s,t , defined on the canonical domain Ω in Figure 4.3, for any s, t ∈ {−, +}. We emphasise that this computation is used only to verify the theory of this paper, and is not needed in the implementation of the domain decomposition solver.
To construct finite-element approximations of these maps, we first derive a variational problem satisfied by them. To do this, we introduce the space V ( Ω), defined as the completion of C ∞ ( Ω) in . Then we define the sesquilinear form where a denotes the sesquilinear form (2.2) defined on Ω. With t ∈ {+, −} and v t ∈ H 1 (Ω t ), let R t v t ∈ V ( Ω) denote the function that coincides with v t on Ω t and is zero elsewhere on Ω. Another application of Green's first identity yields the following result.
Proposition 5.1 (Variational formulation of impedance-to-impedance map). For s, t ∈ {−, +}, let g ∈ L 2 (Γ s ), and let u s ∈ U 0 ( Ω) be the Helmholtz-harmonic function with impedance data g on Γ s and zero elsewhere. Then Motivated by (5.11) and (5.12), we define a finite-element approximation I h,s,t : L 2 (Γ s ) → V h (Γ δ ) to the map I s,t as follows. Analogously to (5.3), for any v h ∈ V h (Γ δ ), we define its node-wise zero extension to all of V h ( Ω) by 0 for all x ∈ N ( Ω\Γ δ ).
Note that R Γ δ ,h v h ∈ V h ⊂ H 1 ( Ω) but is supported only on the union of all elements of the mesh T h that touch Γ δ . Using this, we define I h,s,t by the variational problem where u h,s ∈ V h ( Ω) is the standard finite-element approximation of the function u s (from Proposition 5.1), obtained by solving the homogeneous Helmholtz problem on Ω with impedance data g on Γ s and zero elsewhere. Note that several approximations have been made here. First, in going from (5.12) to (5.13), the test function v t ∈ H 1 (Ω t ) has been replaced by v h ∈ V h (Γ δ ) on the left-hand side and R Γ δ ,h v h on the right-hand side. Moreover the formula (5.11), which requires u ∈ U ( Ω), has been formally applied here with u replaced by u s,h ∈ V h ⊂ U ( Ω). Despite these 'non-conforming' approximations, it can be shown (with details in [37]) that, with · denoting the operator norm from L 2 (Γ s ) → L 2 (Γ δ ), the following convergence result holds.
Thus, the computations of I h,s,t , given in §6, are reliable approximations of I l,j .
A key point in the computation is the realisation that, for any g ∈ L 2 (Γ s ), I h,s,t g = I h,s,t g h , with g h denoting the L 2 -orthogonal projection of g onto V h (Γ s ). (This is because the finite-element solution of the Helmholtz problem only 'sees' the impedance data through its L 2 moments against the finiteelement basis functions.) The operator I h,s,t thus acts only on finite-dimensional spaces, and its norm can be computed by solving an appropriate matrix eigenvalue problem. In §6 this is done using the code SLEPc, within the finite-element package FreeFEM++.

Numerical experiments
In this section, we verify the theoretical results in Theorems 4.13 and 4.16 and Corollaries 4.8, 4.14, 4.15 using the finite-element approximations described in §5. We also perform some extra experiments that provide insight into the performance of the iterative method in situations not covered by the theory. All experiments are on rectangles, the domain is discretized using a uniform triangular mesh with diameter h, and we use the Lagrange conforming element of degree 2. We use mesh diameter h ∼ k −5/4 , which is sufficient to ensure a bounded relative error as k increases [19,Corollary 5.2]. The experiments are implemented using the package FreeFEM++ [41].

Numerical illustration of our theory
In this subsection we consider the 2-d strip domain as in Notation 4.1. The global domain Ω has height H = 1 and length L Ω . For the domain decomposition, we divide Ω into N equal non-overlapping rectangular domains and then extend each subdomain by adding to it neighbouring elements of distance ≤ rL Ω /N away, where r > 0 is a parameter. Thus the interior subdomains have length L = (1 + 2r)L Ω /N , while the end subdomains have length (1 + r)L Ω /N . The global overlap size is δ = 2rL Ω /N . In the first two experiments we examine how the convergence rate depends on the parameters ρ, γ, defined in (4.34), (4.35) and (4.25). Experiment 6.1. (Computation of ρ and γ and convergence of the iterative method as ρ decreases.) Corollaries 4.14 and 4.15 suggest that the convergence rate should improve as ρ decreases, and Theorem 4.16 suggests that the large-k limit of ρ should decrease as δ increases. Table 1 gives values of ρ and γ as functions of k, δ and L as defined in (4.25). These are computed using the method outlined in §5.2. Here r is chosen so that δ = L/3. The top part of Table 1 shows that ρ decreases as L increases, as suggested by Theorem 4.16.
For fixed k, the observed decay rate of ρ is slightly faster than O(δ −1 ). The bottom part of this table shows the corresponding values of γ. Here γ ≤ 1, somewhat smaller than the upper bound predicted by Lemma 4.10. There is a very modest growth of the values of ρ and γ as k increases, for each fixed L; given the lower bound in Theorem 4.16, we expect that the values of ρ in Table 1 are in the preasymptotic regime for k → ∞. To obtain the relative error in the iterative method, we solve the problem (5.2) with right-hand side F = 0, so that the finite-element solution is u h = 0 and the relative error is simply The nodal values of the starting guess u 0 h ∈ V h were chosen to be uniformly distributed in the unit disc in the complex plane. The relative error (6.1) was computed with respect to the first iterate u 1 h ∈ V h 0 , because the initial guess u 0 h is not in this space.    .1 shows that the convergence rate improves as L and hence L Ω increases. This is consistent with Corollary 4.14, which shows that with N fixed and γ bounded, the iterative method is power contractive for small enough ρ. The convergence rate is apparently unaffected by increasing k, a bit better than expected from the k-dependence of ρ in Table 1.
The next experiment investigates the effect of letting the number of subdomains N grow. In this case, Corollary 4.14 guarantees contractivity of T N only for small enough N . However we see that in fact the iterative method continues to work well as N grows. The explanation for this is that, as discussed in §4.4.5, the composite impedance-to-impedance maps are better behaved than the individual ones; this is illustrated in Experiment 6.3 below. Experiment 6.2 (Dependence on N ). We repeat the experiments in Figure 6.1 but instead of N = 3 (i.e, 3 subdomains) we use N = 4, 8, 16. For each N , we choose L Ω so that the sizes of the subdomains and overlaps do not depend on N , and thus ρ and γ remain fixed as N grows. The subdomain length is L = 2 and the overlap is δ = L/3. In Figure 6.2 we plot the relative error histories for k = 20, 40, 80.
The relative error histories show a sudden reduction of the error after each batch of N steps, and, after each such reduction, the convergence rate appears to be higher than before. This can be partially explained by Corollary 4.15; indeed, as the number of iterations n passes through sN for s = 2, 3, . . ., the order of the estimate for the norm of T n increases from O(ρ s−1 ) to O(ρ s ). However this explanation can not be completely rigorous because the coefficient of the powers of ρ in Corollary 4.15 also grows with N . To understand the behaviour of the iterative method better we need to consider composite maps, which is the purpose of Experiment 6.3. Before that, Table 2 gives the average number of iterations needed to reach a relative error of 10 −6 for each of the scenarios depicted in Figure 6.2, computed over 50 random starting guesses. This table clearly indicates that the number of iterations needed to obtain a fixed error tolerance is roughly O(N ) as N grows. We also observe modest improvement in the iteration numbers as k increases; similar results were seen in [38, Table 2: Average number of iterations to reach a relative error of 10 −6 in Figure 6.2 Experiment 6.3 (Robustness to N explained via composite maps). As discussed in §4.4.5, the dominant term in (4.22) with n = N is the j = 1 term (4.52) The goal of this experiment is to show that the behaviour of (4.52) is better than that predicted by estimating its norm by the product of the norms of its components (as in (4.51)). Following (4.54), for N = 4, 8, 16, L = 2 and δ = L/3, we compute , (6.2) and use this as a proxy for (4.52), with this replacement justified by (4.54), and the fact that L N −1 U is a representative element of {p(L, U ) : p ∈ P(N, 1)}. The results in Table 3 show that ζ N remains small and bounded as N increases. Although we have here computed only one term in (4.52), this gives some explanation why the convergence rate of the iterative method remains stable as N increases, as observed in Figure 6.2 and Table 2.  Table 3: Norm of the composite impedance-to-impedance map (6.2), δ = L/3 For the most efficient parallel implementations, the overlap δ should be as small as possible. In our final experiment for the strip domain we therefore study the dependence of the convergence of the iterative method on the overlap parameter. Here, the length of the global domain L Ω = N (L − δ) is chosen so that L = 2, i.e., the subdomains have length 2. In Figure 6.3 we plot the relative error histories. These histories indicate that for small N there is quite a big difference in performance between δ = 2h and the other two choices of δ. However, as N increases the difference between the three choices of overlap becomes less pronounced. With N = 16 we again see clearly the 'staircase' form of the error decay, as in Experiment 6.2. To give some heuristic explanation for Figure 6.3, Tables 4 and 5 provide the analogous results to  Table 3 for the new choices of overlap. As N and k increase, the different choices of overlap all give similar values of ζ N , thus explaining the competitiveness of the small overlap method in this case.   As discussed in §5.2, the parameters ρ, γ are computed above by the finite-element method, and Corollary 5.2 ensures that these approximations converge to the true values of ρ, γ as h → 0. In practice, we compute ρ, γ with h ∼ k −5/4 , which is sufficient for ensuring a bounded error for the Helmholtz problem as k increases by [19,Corollary 5.2]. The following experiment shows that this choice of h also leads to accurate computation of the impedance-to-impedance maps. Experiment 6.5 (Accuracy of the impedance-to-impedance map computation). We compute ρ on the canonical domain Ω depicted in Figure 4.3, with L = 1 and δ = 1/3 for increasing k. In Table 6, we list computed values for ρ(k, 1/3, 1) (i.e., the norm of the left-to-right impedanceto-impedance map -see (4.25)), using mesh sizes h, chosen as decreasing multiples of k −5/4 . A 'brute force' computation of the impedance-to-impedance map by numerical differentiation of the finiteelement solution gave almost identical results to those given in Table 6; we therefore conclude that the computation of ρ is sufficiently accurate when h = k −5/4 .  6.2 Domain decompositions that are not of strip type Experiment 6.6 (Unit square with uniform checkerboard decomposition). We partition Ω := (0, 1) 2 into N 2 non-overlapping equal subsquares each with side length H = 1/N , and then extend to an overlapping cover by adding to each subdomain neighbouring elements that have distance ≤ δ from its boundary (so the actual overlap is 2δ). Tables 7-9 give the iteration counts for the method (1.8)-(1.11) required to achieve a reduction of 10 −6 in the Euclidean norm of the relative residual (with zero right-hand side and starting from a random initial guess), with overlap parameter δ = H/4, H/10, and h, respectively. We also give (in brackets in each table) the number of iterations needed by the corresponding GMRES-accelerated iteration (that is GMRES using the 'ORAS' preconditioner implicit in (5.4), (5.5) -see also [35] ) to obtain a relative residual of 10 −6 . The number of iterations of the iterative method grows somewhere between O(N ) and O(N 2 ), where N 2 is the number of subdomains. In contrast, the number of GMRES iterations seems to grow somewhat more slowly -close to O(N ). It appears that while the iterative method is not effective as a solver when the subdomains do not contain enough wavelengths, the GMRES method continues to function acceptably and is robust or even improving as k increases. (Related theory and observations are given around [38, Table 3].) In contrast to the strip domain case, we also observe that reducing the overlap does significantly affect the performance of both the iterative method and GMRES.    . We consider the same set-up as in Experiment 6.6, but instead of using the checkerboard domain decomposition, we use METIS to generate a non-overlapping domain decomposition -see Figure 6.4 for plots for N = 4, 16, 64 subdomains -and then we extend to an overlapping cover in the same way as before. Tables 10-12 give iteration counts for the iterative method (and GMRES in brackets) for each choice of δ. The iteration counts behave similarlly to those given in Experiment 6.6. Again, we notice the almost-robustmess of the GMRES method as k increases for all choices of δ, and particularly for generous overlap.      To finish the paper we include some remarks on how the fundamental general theory developed here could be extended to get convergence estimates for more general decompositions. Here we only (and rather tentatively) discuss the case of a 2×2 checkerboard decomposition depicted in Figure 6.5. More general decompositions could be approached by extending this reasoning, although, we also admit such reasoning will become increasingly complex for more general domain decompositions. Consider the unit square divided into four quarters, yielding non-overlapping subdomains labelled 1,2,3,4 in Figure 6.5. These are extended to overlapping subdomains Ω , = 1, 2, 3, 4. We depict only Ω 1 (in red) and Ω 2 (in blue) in the figure.  from Ω j to Ω , with m 1 = m 2 , m 2 = m 3 , m 3 = m 4 , and m 4 = m 5 . The first thing to note is that many of the component products in (6.3) have already been investigated earlier in this paper. For example, consider the product T 1,2 T 2,1 ; by Theorem 3.9 we see that, as k gets larger then the norm of this product is well-estimated by the norm of the impedance-toimpedance map I Γ 2,1 →Γ 1.2 . This map operates only in the x−direction and was analysed in §4.4.2. In fact it corresponds to the left-to-right map defined by (a) finding the Helmholtz-harmonic function on Ω 2 with data on ∂Ω 2 ∩ Ω 1 and (b) evaluating the right-facing impedance data of the solution on ∂Ω 1 ∩ Ω 2 . This map was analysed in §4 and shown to have norm of order O(ρ). Similar remarks apply to T 2,1 T 1,2 , T 1,4 T 4,1 , etc. However the 'diagonal switches' T 1,3 T 3,1 , T 3,1 T 1,3 , T 1,4 T 4,1 and T 4,1 T 1,4 are genuinely two-dimensional and the properties of the corresponding impedance maps would need to be studied numerically.
More generally we could conjecture that if the sequence in (6.4) visits the same subdomain twice then the norm of the corresponding term in (6.3) is small in size. This conjecture would have to be examined numerically, but would then give estimates for all off-diagonal terms in (6.3).
Finally, one would have to analyse the cycles which appear in the diagonal terms in (6.3) e.g. the term corresponding to 1 → 2 → 3 → 4 → 1.
This further analysis is outside the scope of the present paper. However, while it is clear from experiments that higher powers of T are still contractive in the checkerboard (and more general) cases, the convergence profile (at least in certain norms [35]) does not exhibit the same jumps when one passes from n = sN to n = (s + 1)N as are present in the strip domain case (see Figures 6.1, 6.2).