# Sweeping preconditioners for the iterative solution of quasiperiodic Helmholtz transmission problems in layered media

###### Abstract

We present a sweeping preconditioner for quasi-optimal Domain Decomposition Methods (DDM) applied to Helmholtz transmission problems in periodic layered media. Quasi-optimal DD (QO DD) for Helmholtz equations rely on transmission operators that are approximations of Dirichlet-to-Neumann (DtN) operators. Employing shape perturbation series, we construct approximations of DtN operators corresponding to periodic domains, which we then use as transmission operators in a non-overlapping DD framework. The Robin-to-Robin (RtR) operators that are the building blocks of DDM are expressed via robust boundary integral equation formulations. We use Nyström discretizations of quasi-periodic boundary integral operators to construct high-order approximations of RtR. Based on the premise that the quasi-optimal transmission operators should act like perfect transparent boundary conditions, we construct an approximate LU factorization of the tridiagonal QO DD matrix associated with periodic layered media, which is then used as a double sweep preconditioner. We present a variety of numerical results that showcase the effectiveness of the sweeping preconditioners applied to QO DD for the iterative solution of Helmholtz transmission problems in periodic layered media.

Keywords: Helmholtz transmission problems, domain decomposition methods, periodic layered media, sweeping preconditioners.

AMS subject classifications: 65N38, 35J05, 65T40,65F08

## 1 Introduction

The numerical simulation of interactions between electromagnetic, acoustic, and elastic waves with periodic layered media has numerous applications in the fields of optics, photonics, geophysics [5]. Given the important technological applications of periodic layered media, the simulation of wave propagation in such environments has attracted significant attention [6, 15, 19, 23, 22]. Regardless of the type of discretization (finite elements, finite differences, boundary integral operators), iterative solvers are the preferred method of solution especially for high-frequency layered configurations that involve large numbers of layers that may contain inclusions. The iterative solution of high-frequency Helmholtz and Maxwell equations in complex media is a challenging computational problem [11]. A successful strategy to tackle this problem relies on sweeping preconditioners [9]. We present in this paper a sweeping preconditioner for a Domain Decomposition formulation of Helmholtz transmission problems in two dimensional periodic layered media.

Domain Decomposition Methods are natural candidates for the solution of Helmholtz transmission problems in periodic layered media [23, 22, 19]. Local subdomain solutions (the subdomains may or may not coincide with the periodic layers) are linked iteratively via Robin type transmission conditions defined on inter domain interfaces. Ideally, the transmission operators should act as transparent boundary conditions that allow information to flow out of each subdomain with very little information being reflected back. As such, for a given subdomain, optimal transmission operators on the subdomain interface consist of Dirichlet-to-Neumann (DtN) operators associated with the adjacent subdomain that shares the same interface. In practice, the transmission operators are constructed via various approximations of DtN operators that rely either on Fourier calculus [2, 10] or Perfectly Matched Layers [25, 24]; the ensuing DD are referred to as Quasi-Optimal DD (QO DD) or Optimized Schwartz Methods [11].

The main scope of this paper is QO DD for the solution of Helmholtz transmission problems in periodic layered media separated by grating profiles (i.e. graphs of periodic functions). We present two strategies of subdomain partition whereby (1) the subdomains coincide with the layer subdomains and the subdomain interfaces coincide with the grating profiles of material discontinuity of the layered medium; and (2) the subdomains consist of horizontal strips whose flat boundaries do not intersect any of the grating profiles of material discontinuity. We note that the DD partition strategy (2) is only applicable to layered media configurations where the width of the layers is larger than the roughness of their interfaces. In each subdomain a local Helmholtz quasi-periodic equation with generalized Robin conditions must be solved (the wavenumber may be discontinuous in case (2)), and generalized Robin data on the subdomain boundaries are linked with those corresponding to the adjacent subdomain. The generalized Robin data corresponding to a given subdomain is defined in terms of transmission operators that are approximations of DtN operators corresponding to the adjacent subdomain. Such approximations of periodic DtN operators can be obtained via high-order shape perturbation/deformation series in case (1) [20]. Specifically, using as a small parameter the roughness/elevation height of the grating, the periodic DtN operators are expressed as a perturbation series whose terms can be computed recursively. The zeroth order terms of the perturbation series coincide with DtN of layered domains with flat interfaces, which can be written explicitly in terms of Fourier multipliers. In the case of the subdomain partition (2), since the subdomain interfaces are flat, the transmission operators are chosen to be the aforementioned Fourier multipliers. We establish that the ensuing QO DD corresponding to both subdomain partitions are equivalent to the original transmission problem, with the caveat that the roughness of the grating profiles must be small enough for the subdomain partition (1).

The exchange of Robin data amongst the subdomains in DD is realized via quasi-periodic Robin-to-Robin (RtR) maps that map incoming to outgoing subdomain Robin data. Following the methodology introduced in [22], we express quasi-periodic RtR operators in terms of robust boundary integral equation formulations. The discretization of the RtR maps is realized by extending the high-order Nyström method based on trigonometric interpolation and windowing quasi-periodic Green functions [22] to the case of DtN transmission operators. Since the terms in the shape deformation series expansions of DtN operators are expressed in terms of Fourier multipliers [20], the discretization of the QO transmission operators is straightforward within the framework of trigonometric interpolation. Using Nyström discretization RtR matrices, we discretize the QO DD formulation for layered transmission problems in the form a block tridiagonal matrix which we invert using Krylov subspace iterative methods. However, the numbers of iterations required for the solution of QO DD linear systems grows with the number of layers, especially for high frequencies at high-contrast layer media configurations. In order to alleviate this situation, we construct a double sweep preconditioner based on an approximate LU factorization of the block tridiagonal QO DD matrix. As it was very nicely explained in the recent contribution [11], all of the effective sweeping preconditioners that have been introduced in the last decade [25, 24, 9] can be elegantly described in terms of optimized DD with layered subdomain partitions and approximate LU of the ensuing block tridiagonal DD matrices. The key insight in our construction of the LU factorization is related to the observation that if the transmission operators were to behave as perfect transparent boundary conditions, certain blocks in the QO DD matrix can be approximated by zero [24]. This approximation renders the LU factorization particularly simple as it bypasses altogether the need for inversions of block matrices.

We present a variety of numerical results that highlight the benefits of QO DD formulations for the solution of transmission problems in periodic layered media, as well as the effectiveness of the sweeping preconditioners in the presence of large numbers of layers at high frequencies. With regards to the latter regime, we find that the sweeping preconditioners used in conjunction with QO DD with horizontal strip subdomain partitions is particularly effective. We mention that the quasi-optimal transmission operators based on Fourier square-root principal symbols approximations of DtN operators has been already used in several contributions [2, 24, 13]; we simply extend the square root Fourier calculus to the periodic setting and incorporate it within the high-order shape deformation expansions technology introduced in [20]. Furthermore, the construction of the sweeping preconditioners that we employ in this paper was originally introduced in [24] and further elaborated upon in [11]. The main contributions of this paper are (a) the integration of these two important ideas within a high-order Nyström discretization of robust quasi-periodic boundary integral equation formulations of RtR maps, as well as (b) the analysis of the quasi-periodic QO DD. A fully parallel implementation of the sweeping preconditioner applied to QO DD for transmission Helmholtz equations in layered media is currently being developed. Also, the generalization of the DD with horizontal strip subdomain partitioning is currently under investigation; this would entail careful treatment of cross points (i.e. points on the subdomain boundaries where the wavenumbers are discontinuous), which we plan to pursue along the lines of the contribution [13].

The paper is organized as follows: in Section 2 we present the formulation of Helmholtz transmission problems in periodic layered media; in Section 3 we present QO DD formulations of the periodic Helmholtz transmission problem; we continue in Section 4 with the construction of quasi-optimal transmission operators based on high-order shape perturbation series; we show in Section 5 a means to express the QO DD RtR operators in terms of robust quasi-periodic boundary integral equation formulations, which, in turn, enable us to analyze the equivalence between the QO DD formulations and the original Helmholtz transmission problems; and we conclude in Section 6 with the construction of the sweeping preconditioner and with a presentation of a variety of numerical results that illustrate the effectiveness of these preconditioners in the context considered in this paper.

## 2 Scalar transmission problems

We consider the problem of two dimensional quasi-periodic scattering by penetrable homogeneous periodic layers. We assume that the layers are given by for and and , and all the functions are periodic with principal period , that is for all , and . We assume that the medium occupying the layer is homogeneous and its permitivity is ; the wavenumber in the layer is given by . We assume that a plane wave where impinges on the layered structure, and we are interested in looking for -quasi-periodic fields that satisfy the following system of equations:

(2.1) |

where denote the unit normals to the boundary pointing to the exterior of the subdomain (i.e. for the domain we define on , for the domains we define on and on , and finally, for the domain we define on ) and is the Dirac distribution supported on . We also denote by on for the domain , for the domains , on and on , and finally, for the domain , on ; we point out that for all domains . We note that with this convention on unit normals we have that as well as on . We also assume that and in equations (2.1) are radiative in and respectively. The latter requirement amounts to expressing the solutions and in terms of Rayleigh series

(2.2) |

and

(2.3) |

where , and , where the branches of the square roots in the definition of and in such a way that , and that the branch cut coincides with the negative imaginary axis. We assume that the wavenumbers and the quantities in the subdomains are positive real numbers. In electromagnetic applications, or depending whether the incident radiation is TE or TM. For the sake of simplicity, we consider in this contribution the case ; extensions to general positive are straightforward.

## 3 Domain decomposition approach

The transmission problem (2.1) can be formulated via boundary integral equations (BIEs) [1, 6] or via non-overlapping Domain Decomposition (DD) [22, 19]. Upon discretization, both the BIE and DD amount to solving block tridiagonal linear systems. In the case of large numbers of layers, the ensuing (large) linear systems are solved via direct methods [6, 22] that rely on Schur complements. As such, the applicability of direct solvers for the numerical solution of the transmission problem (2.1) is limited by the size of the Schur complements. Iterative solvers, on the other hand, do not suffer from the aforementioned size limitations, yet are challenged by the presence of significant multiple scattering, especially in high-contrast multi-layer configurations at high frequencies. In the high-frequency regime, relevant to technological applications, efficient preconditioners are needed in order to alleviate multiple scattering. The main scope of this contribution is to present such a preconditioner (referred to as the sweeping preconditioner [25, 24, 9]) in the context of DD formulation of quasiperiodic transmission problems.

The main idea of DD is to divide the computational domain into subdomains, and to match subdomain quasiperiodic solutions of Helmholtz equations via Robin type transmission conditions on the subdomain interfaces. We consider in what follows two strategies of partitioning the computational domain into non-overlapping subdomains: the most natural one in which the DD subdomains coincide with the layer domains , and an alternative one in which the subdomains are horizontal strips. We present in what follows the details of the first subdomain partitioning strategy mentioned above.

### 3.1 DD with subdomains

A natural non-overlapping domain decomposition approach for the solution of equations (2.1) consists of solving subdomain problems in with matching Robin transmission boundary conditions on the common subdomain interfaces for . Indeed, this procedure amounts to computing -quasiperiodic subdomain solutions:

(3.1) | |||||

where , are certain transmission operators for , and denote normal derivatives with respect to the non-unit normals . In addition, we require that and are radiative. We have chosen to double index the interfaces between layer subdomains: the first index refers to the index of the layer , whereas the second index denotes the index of the layer adjacent to the layer so that is the interface between and . Here and in what follows denote Sobolev spaces of -quasiperiodic functions/distributions defined on the periodic interface ; the definition of these spaces is given in terms of Fourier series [22].

Heuristically, in order to give rise to rapidly convergent iterative DD, the transmission operators ought to be good approximations of the restriction to of the DtN operator associated with the -quasiperiodic Helmholtz equation in the domain with wavenumber . This requirement explains why the indices are reversed in the definition of the transmission operators. In addition, the transmission operators and ought to be selected to meet the following two criteria: (1) the subdomain boundary value problems that incorporate these transmission operators in the form of generalized Robin boundary conditions are well-posed for all frequencies, and (2) the DD matching of the generalized Robin data on the interfaces of material discontinuity (which coincide with the layer boundaries) is equivalent to the original transmission conditions (2.1) on the same interfaces.

Specifically, with regards to the issue (1) above, we require that for a given layer domain with , the following -quasiperiodic boundary value problem is well-posed:

(3.2) | |||||

where and are generic -quasiperiodic functions defined on and respectively. The following coercivity properties

(3.3) |

for all in terms of the and duality pairings are sufficient conditions for guaranteeing the well posedness of the boundary value problems (3.2). Indeed, this can be established easily by an application of the Green’s identities in the domain . In the case of the semi-infinite domain , we require that the following -quasiperiodic boundary value problem is well-posed:

(3.4) | |||||

where is a -quasiperiodic function defined on . The coercivity property

(3.5) |

suffices to establish the well posedness of the boundary value (3.4). The latter fact can be established via the same arguments as those in Theorem 3.1 in [22]. A similar coercivity condition imposed on the operator ensures the well posedness of the analogous -quasiperiodic boundary value problem on the semi-infinite domain .

Returning to the requirement (2) above, we ask that the DD matching of the generalized Robin data

is equivalent to the continuity conditions

It can be immediately seen that the equivalence in part (2) is guaranteed provided that is an injective operator. Under the assumption that the coercivity properties (3.5) hold, it follows that

and thus the operators are injective for all . Thus, the coercivity properties (3.5) ensure that both requirements (1) and (2) above are met. We postpone the discussion on the selection of the transmission operators and and we formulate the DD system (3.1) in matrix operator form. To that end, we define certain RtR operators associated with the boundary value problems (3.2). Specifically, we define the RtR map in the following manner:

(3.6) |

Also, associated with the boundary value problem (3.4) posed in the semi-infinite domain we define the RtR map in the form

(3.7) |

The RtR map corresponding to the domain is defined in a similar manner to but for a boundary data defined on .

With these notations in place, the DD formulation (3.1) seeks to find the generalized Robin data associated with each interface

as the solution of the following operator linear system

(3.8) |

where and the right-hand-side vector has zero components with the exception of the first component

and the DD matrix is a tridiagonal block matrix given in explicit form by

(3.9) |

where

(3.10) |

We present in what follows a different strategy of domain decomposition whereby the subdomains are horizontal strips.

### 3.2 DD with stripes subdomains

An alternative DD possibility is to partition the computational domain using horizontal stripes. We restrict ourselves to cases where the layer domains are wide enough so that each periodic interface can be contained in a horizontal strip that does not intersect any other interface . Under this assumption, these horizontal stripes constitute the DD subdomains—see Figure 2 for a depiction of the partitioning in the case of three layers (i.e. ). In general, however, a domain decomposition into horizontal stripes might require that an interface intersect a (flat) boundary of a strip; we leave this challenging scenario for future considerations.

Assuming that there exist real numbers such that for all we have that and , then we can partition into a union of nonoverlapping horizontal strips , where the strip domains are defined as , , and . Using the domain decomposition into layered stripes we seek -quasiperiodic solutions of the following system of PDEs

(3.11) | |||||

where and denotes the jump of the function across the interface . We require that the transmission operators have the following mapping properties and and satisfy coercivity properties similar to those in equations (3.3).

The coercivity properties of the transmission operators and are needed to ensure the well-posedness of the following subdomain PDEs

(3.12) | |||||

for all as well as those posed in the semi-infinite domains and respectively. Associated to the Helmholtz transmission problem (3.12) is the RtR operator defined below

(3.13) |

The DD formulation (3.11) then seeks to find the generalized Robin data associated with each interface

as the solution of the following operator linear system

(3.14) |

where the DD matrix is similar to that defined in equation (3.9), and the right-hand-side vector has zero components with the exception of the first component

Having described two possible DD strategies for the solution of quasi-periodic Helmholtz transmission problems (2.1), we present next a methodology based on Fourier calculus to construct quasi-optimal transmission operators.

## 4 Construction of quasi-optimal transmission operators based on shape perturbation series

We present in what follows a perturbative method to construct quasi-optimal transmission operators and for corresponding to the DD formulation (3.1). To this end, given a generic -periodic profile function we define the periodic interface and the semi-infinite domains and respectively . We assume that the profile function can be expressed in the form , where the -periodic function is analytic (it actually suffices that the profile function is Lipschitz). We employ a perturbative approach [20] to construct approximations of the DtN operator corresponding to the following boundary value problem in the domains :

(4.1) | |||||

where are radiative in the domains and is a -quasiperiodic function defined on , and is the normal to pointing into the domain . Under the assumptions above, the DtN operators are analytic in the shape perturbation variable , and thus we seek the operator in the form of the perturbation series

(4.2) |

where the operators can be computed via explicit recursive formulas [20]. Let us denote by the radius of convergence of the perturbation series (4.2). Following [20], we present next the recursive formulas that lead to closed form expressions of the operators . First, given an -quasiperiodic function which can be represented as

we define the Fourier multiplier operator

(4.3) |

Then, it can be shown that the operators in the perturbation series (4.2) can be computed via the recursion

(4.4) | |||||

where . We note that for all , which is to say that none of the higher-order terms is more regular that the lowest term . Furthermore, since the recursions (4) may lead to significant subtractive cancellations, more stable expressions of the operators were proposed. Indeed, using the commutator

it can be shown that the low-order term corrections can be expressed in the equivalent form

(4.5) |

and

(4.6) |

where . However, the calculation of high-order correction terms via the stable recursions above becomes quite cumbersome. As such, a different strategy based on changes of variables (that straighten out the boundary ) and DtN corresponding to variable coefficient Helmholtz equation in half-planes is advocated in [20] for stable computations of DtN maps. Given that our motivation is to construct readily computable DD transmission operators that are approximations of DtN operators, we will restrict to low-order terms in the perturbation series (4.2), which, as discussed above, can be computed by explicit and stable recursions.

In order to meet the coercivity requirements (3.3), we complexify the wavenumber in the form and we define

(4.7) |

using formulas (4.5) and (4.6) for the definition of the operators in equation (4.7). Indeed, we establish the following

###### Lemma 4.1

Provided that is small enough, the following coercivity property holds

for all .

Proof. By the construction of the Fourier multiplier operator we have that