Computational high frequency scattering from high contrast heterogeneous media

TL;DR

This paper develops a multiscale computational method for simulating high-frequency acoustic wave scattering in strongly heterogeneous media with high contrast, capturing complex physical phenomena efficiently.

Contribution

It introduces a novel multiscale approach with rigorous error analysis for high contrast, high frequency wave propagation in non-periodic heterogeneous structures.

Findings

Method accurately captures wave scattering in high contrast media.

Numerical experiments confirm theoretical error estimates.

Approach effectively models physical phenomena in complex materials.

Abstract

This article considers the computational (acoustic) wave propagation in strongly heterogeneous structures beyond the assumption of periodicity. A high contrast between the constituents of microstructured multiphase materials can lead to unusual wave scattering and absorption, which are interesting and relevant from a physical viewpoint, for instance, in the case of crystals with defects. We present a computational multiscale method in the spirit of the Localized Orthogonal Decomposition and provide its rigorous a priori error analysis for two-phase diffusion coefficients that vary between $1$ and very small values. Special attention is paid to the extreme regimes of high frequency, high contrast, and their previously unexplored coexistence. A series of numerical experiments confirms the theoretical results and demonstrates the ability of the multiscale approach to efficiently capture relevant physical phenomena.

Equations

$$-\mathrm{div}(A\nabla u)-k^{2}u$$

$$\nabla u\cdot n-iku$$

$$\mathcal{B}(u,v):={\int }_{\Omega }A\nabla u\cdot \nabla \overline{v}-k^{2}u\overline{v}dx-ik{\int }_{\partial \Omega }u\overline{v}d\sigma ={\int }_{\Omega }f\overline{v}dx\forall v\in H^1(\Omega ).$$

$$\Vert u{\Vert }_{0,A}$$

$$|u{|}_{1,A}$$

$$\Vert u{\Vert }_{1,A,k}$$

$$k^{1/2}\Vert v{\Vert }_{\partial \Omega }\lesssim \Vert v{\Vert }_{1,A,k}.$$

$$\Vert v{\Vert }_{\partial \Omega }\lesssim {\epsilon }^{-1}\Vert v{\Vert }_{1,A,k}.$$

$$\Vert u{\Vert }_{1,A,k}\lesssim C_{\mathrm{stab}}(k,\epsilon )\Vert f{\Vert }_{L^2(\Omega )}.$$

$$\underset{v\in H^1(\Omega )}{\inf }\underset{\psi \in H^1(\Omega )}{\sup }\frac{\mathrm{Re}\mathcal{B}(v,\psi )}{\Vert v{\Vert }_{1,A,k}\Vert \psi {\Vert }_{1,A,k}}\ge \gamma (\epsilon ,k)$$

$$\text{N}(T)=\bigcup _{K\in {\mathcal{T}}_H,T\cap K\ne \varnothing }K.$$

$$\underset{T\in {\mathcal{T}}_H}{\max }\mathrm{card}\{K\in {\mathcal{T}}_H:K\subset {\text{N}}^m(T)\}\le C_{\mathrm{ol},m}.$$

$$|I_{H}v{|}_{1,A,T}$$

$$\Vert w{\Vert }_{0,A,T}$$

$$\Vert I_{H}v{\Vert }_{1,A,k,T}$$

$$I_H(v)=\sum _{z\in {\mathcal{N}}_{\mathrm{int}}}{P}_{A,{\omega }_z}(v)(z){\lambda }_z$$

$${\int }_{{\omega }_z}A{P}_{A,{\omega }_z}(v){\psi }_{H}dx={\int }_{{\omega }_z}Av{\psi }_{H}dx\text{for all}{\psi }_H\in V_H{|}_{{\omega }_z}.$$

$$I_H^1(v):=\sum _{z\in {\mathcal{N}}_{\mathrm{int}}}{P}_{1,{\omega }_z}(v)(z){\lambda }_z,$$

$${\mathcal{B}}_{{\text{N}}^m(T)}({\mathcal{Q}}_{T,m}{v}_H,w)=-{\mathcal{B}}_T(v_H,w)\forall w\in W({\text{N}}^m(T)).$$

$${\mathcal{Q}}_m=\sum _{T\in {\mathcal{T}}_H}{\mathcal{Q}}_{T,m}.$$

$$\mathcal{B}((\mathrm{id}+{\mathcal{Q}}_m)u_{H,m},(\mathrm{id}+{\mathcal{Q}}_m^{\ast })v_H)=(f,(\mathrm{id}+{\mathcal{Q}}_m^{\ast })v_H)\forall v_H\in V_H.$$

$$u_{\mathrm{LOD},m}:=(\mathrm{id}+{\mathcal{Q}}_m)u_{H,m}.$$

$$\mathcal{B}(\mathcal{Q}v,w)=-\mathcal{B}(v,w)\forall w\in W.$$

$$kH{\epsilon }^{-1}\lesssim 1,$$

$$|w{|}_{1,A}\le \Vert w{\Vert }_{1,A,k}\lesssim |w{|}_{1,A}.$$

$$\mathrm{Re}\mathcal{B}(w,w)\ge \alpha |w{|}_{1,A}^2\gtrsim \alpha \Vert w{\Vert }_{1,A,k}^2\forall w\in W.$$

$$k\Vert w\Vert \lesssim k{\epsilon }^{-1}\Vert w{\Vert }_{0,A}$$

$$|(\mathcal{Q}-{\mathcal{Q}}_m)v_H{|}_{1,A}\lesssim C_{\mathrm{ol},m}^{1/2}{\beta }^m|v_H{|}_{1,A}.$$

$$\Vert (\mathcal{Q}-{\mathcal{Q}}_m)v_H{\Vert }_{0,A}\lesssim C_{\mathrm{ol},m}^{1/2}H{\beta }^m\Vert v_H{\Vert }_{1,A}.$$

$$\mathcal{B}(w,{\mathcal{Q}}^{\ast }v)=-\mathcal{B}(w,v)\forall w\in W$$

$$\mathcal{B}(I_{H}u,(\mathrm{id}+{\mathcal{Q}}^{\ast })v_H)$$

$$=(f,(\mathrm{id}+{\mathcal{Q}}^{\ast })v_H).$$

$$\mathcal{B}(u_{\mathrm{LOD}},(\mathrm{id}+{\mathcal{Q}}^{\ast })v_H)=(f,(\mathrm{id}+{\mathcal{Q}}^{\ast })v_H)\forall v_H\in V_H.$$

Full text

Abstract. This article considers the computational (acoustic) wave propagation in strongly heterogeneous structures beyond the assumption of periodicity. A high contrast between the constituents of microstructured multiphase materials can lead to unusual wave scattering and absorption, which are interesting and relevant from a physical viewpoint, for instance, in the case of crystals with defects. We present a computational multiscale method in the spirit of the Localized Orthogonal Decomposition and provide its rigorous a priori error analysis for two-phase diffusion coefficients that vary between $1$ and very small values. Special attention is paid to the extreme regimes of high frequency, high contrast, and their previously unexplored coexistence. A series of numerical experiments confirms the theoretical results and demonstrates the ability of the multiscale approach to efficiently capture relevant physical phenomena.

AMS subject classifications. 35J05, 65N12, 65N15, 65N30

1 Introduction

Wave propagation in heterogeneous high-contrast materials has received a growing interest in recent years because the combination of microstructures and high contrast can lead to unusual material properties such as opening band gaps, artificial magnetism, or surface plasmons, see, for instance, [19, 34, 47, 52]. In periodic structures, homogenization techniques based on two-scale convergence [1, 36] or Bloch-wave techniques [2, 3, 15] are able to reduce the problem’s complexity and to provide effective models for the overall macroscopic behavior of the wave. In practical applications, however, the periodic structure may easily be destroyed by local defects due to production errors. Perturbation theory for infinite periodic structures indicates that the combination of certain local defects and wave numbers may (suddenly) lead to localized waves, whereas in the perfectly periodic case the corresponding wave cannot propagate, see [4, 34]. While this can be highly favorable, local defects may also destroy the desired tuned (unusual) properties of the periodic structures. However, the above described analytical examination of the phenomena is quite limited to the periodic case with (mild) perturbations, where each case has to be considered separately. Moreover, it is not clear how these asymptotic results can predict the particular behavior for a given (finite) structure.

Therefore, numerical simulations of wave propagation and scattering in heterogeneous high-contrast materials are crucial to predict and also understand more situations and phenomena. Yet, a direct discretization of realistic problems using standard (finite element) methods easily exceeds today’s computational resources because the mesh has to resolve all the material inhomogeneities. Additional challenges and restrictions on the mesh size are encountered in the high-frequency regime by the so-called pollution effect, see [5]. Therefore, multiscale methods need to be employed to reliably simulate wave propagation in high-contrast heterogeneous materials.

By now, the challenges for wave propagation problems on the one side and for high-contrast coefficients on the other side mostly have been studied separately. To reduce the pollution effect for the Helmholtz equation with constant or low-contrast coefficients, a number of approaches has been designed and analyzed, such as (hybridizable) discontinuous Galerkin methods [12, 27], the $hp$-version of the finite element method [21, 40, 41, 42], (plane wave) Trefftz methods [24, 33, 48], or the Localized Orthogonal Decomposition [10, 23, 50]. High-contrast coefficients mainly have been studied for elliptic problems using multiscale finite element methods [13, 18, 45] and the Localized Orthogonal Decomposition [29, 51], just to name a few. Finally, we mention that high-contrast Helmholtz problems have been studied based upon homogenization results in [44]. We shall emphasize that high contrast in the diffusion coefficient can be triggered by the degeneration of the upper or the lower bound of the elliptic operator. In the context of time-harmonic wave propagation, the observed effects are quite different in the two regimes, see, for instance, [34]. Locally, the variation of the diffusion coefficient can be rephrased in terms of a variable wave number. Small diffusion coefficients correspond to large wave numbers which amplifies the aforementioned pollution effect and leads to much more critical conditions on the minimal numerical resolution. That is why we focus in this manuscript on the specific setting where the diffusion coefficient takes a small value $\varepsilon^{2}$ on one part of the domain and $1$ on the complement.

The paper presents a computational multiscale method in the spirit of the Localized Orthogonal Decomposition (LOD) and provides its rigorous analysis for high-contrast Helmholtz problems beyond periodicity and scale separation. We build upon previous ideas for high-frequency wave propagation in [10, 23, 50] and high-contrast elliptic equations in [11, 29, 51]. The combination of these techniques and approaches yields a (quasi-)localized method which does not require structures in the coefficients such as periodicity. Furthermore, a scale separation between the wave length and the size of individual “valleys” (where the coefficient is very small) is not required. The multiscale method can be viewed as a generalized finite element method with special multiscale basis functions. As standard finite element functions on a coarse mesh alone do not yield a faithful approximation space, problem-dependent multiscale functions are added. The latter are defined as solutions of local fine-scale problems. Our numerical analysis shows that the (coarse) mesh, and, hence, the dimension of the space, is coupled linearly to the effective wave number. Moreover, the size of the fine-scale problems to solve scales only logarithmically with this effective wave number. This generalizes the results of [10, 23, 50] for the constant or low-contrast coefficient case. We carefully address the high contrast in our analysis which is crucial to obtain the mildest dependencies on the contrast (encoded in the quantity $\varepsilon$) as possible. We provide a priori error estimates in various norms, which can be motivated from the periodic (high contrast) case and are of the expected order. The error analysis contains the stability constant of the exact solution (see (2.5) and (2.6)). Unfortunately, only little is known in the general setting about bounds for this constant that are explicit in the wave number and the coefficients. Hence, the multiscale method is proved to be contrast- and wave number-robust only under the assumption that the solution depends moderately on these parameters.

The paper is organized as follows. In Section 2 we introduce our model problem and illustrate properties of the solution in the periodic case. We introduce the multiscale method with all necessary notation in Section 3 and analyze its well-posedness and a priori errors in Section 4. The main proofs are detailed in Section 5. Finally, numerical experiments in Section 6 confirm our predicted theoretical results and show the applicability of our method for physically relevant cases concerning band gap simulations and localized defects.

2 Helmholtz problem with high contrast

In this section, we are concerned with the Helmholtz problem with high contrast in general. First, we formulate our model setting in Section 2.1 and also introduce the necessary notation. In Section 2.2, we then explain in detail the considered set-up by motivating it from homogenization theory and Bloch-wave analysis. The discussion moreover sheds a light on what approximation results we can expect for our multiscale method and compares it to low-contrast Helmholtz problems.

2.1 Model problem

Let $\Omega\subset\mathbb{R}^{d}$ be a bounded polyhedral Lipschitz domain. We consider the following model problem: Find $u:\Omega\to\mathbb{C}$ such that

[TABLE]

Here, we assume $f\in L^{2}(\Omega)$ and $k\geq k_{0}>0$. Inhomogeneous (Robin) boundary conditions on $\partial\Omega$ can be treated following the theory of [31], which we omit for simplicity. Other types of boundary conditions [23, 50] and also perfectly matched layers [22] have been studied for homogeneous media and their incorporation in the present work is possible as well. The scalar diffusion coefficient $A$ is piece-wise constant with respect to a quadrilateral background mesh $\mathcal{T}_{\varepsilon}$ with mesh size $O(\varepsilon)$ and $0<\varepsilon\ll 1$. On each quadrilateral, $A$ takes either the value $\varepsilon^{2}$ or $1$. We introduce the non-overlapping partition $\Omega=\Omega_{1}\cup\Omega_{\varepsilon}$, where $\Omega_{1}$ and $\Omega_{\varepsilon}$ are the parts of the domain where $A$ takes the value $1$ or $\varepsilon^{2}$, respectively. We assume that $\Omega_{\varepsilon}\subset\subset\Omega$, so that $A$ does not appear in the Robin boundary condition. This set-up requires some comments. Although we fix the two values of $A$, the present setting still allows much freedom since the geometry of $\Omega_{\varepsilon}$ is relatively flexible. In particular, no periodicity is assumed. We emphasize that the scaling of $A$ is independent of the space dimension. The special choice of $1$ and $\varepsilon^{2}$ induces a high contrast between the two materials, which is important to obtain unusual macroscopic effects in (periodic) homogenization theory, see Section 2.2 below for a detailed discussion. However, we also underline that the results of the present contribution can directly be transferred to other scalings of $A$, see Remark 4.7. We stick with the above choice of $\varepsilon^{2}$ in $\Omega_{\varepsilon}$ for simplicity and since this seems to be the physically (most) relevant set-up. As mentioned in the introduction, we stress that the present choice of high contrast deals with very small values of $A$ in parts of the domain (in contrast to, for instance, [13]). This scaling in particular implies an almost degeneracy of the equation and is therefore harder than the case where $A$ may vary between $1$ and very large values, see Remark 4.7. Finally, note that the assumption of a piece-wise constant $A$ can be relaxed if it is still possible to construct a suitable interpolation operator as described in Section 3.1 for that diffusion coefficient.

Throughout the whole article, we use standard notation on Sobolev spaces. All functions are complex-valued if not mentioned otherwise. The complex conjugate of $v$ is denoted by $\overline{v}$. We write $a\lesssim b$ in short for $a\leq Cb$ with a constant $C$ independent of $k$, $\varepsilon$, the mesh size $H$, and the oversampling parameter $m$. Similarly, $a\approx b$ stands for $a\leq Cb$ and $b\leq C^{\prime}a$ with constants $C,C^{\prime}$ independent of $k$, $\varepsilon$, $H$, and $m$. We then consider the (weak) problem: Find $u\in H^{1}(\Omega)$ such that

[TABLE]

Let $\|\cdot\|$ denote the standard $L^{2}(\Omega)$-norm, we then define the following weighted (semi) norms:

[TABLE]

If necessary, we will indicate by an additional subscript the domain where the norm is taken. Note that $\|\cdot\|_{1,A,k}$ represents the energy norm associated with the sesquilinear form $\mathcal{B}$. If we assume that $\operatorname{dist}(\Omega_{\varepsilon},\partial\Omega)$ is bounded from below by an $\varepsilon$-independent constant, we can deduce the following trace inequality (cf. [9])

[TABLE]

Otherwise, if $\Omega_{\varepsilon}$ is arbitrarily closed to $\partial\Omega$, we would obtain an additional $\varepsilon$-weight in the trace inequality of the following form

[TABLE]

Since the trace inequality is needed for the continuity of $\mathcal{B}$ with respect to the energy norm, we assume that $\operatorname{dist}(\Omega_{\varepsilon},\partial\Omega)$ is bounded from below by an $\varepsilon$-independent constant. Then, the continuity constant of $\mathcal{B}$ is $\varepsilon$-independent, which simplifies the analysis of the multiscale method below, cf. also Remark 4.6.

The unique continuation principle and Fredholm theory provide existence and uniqueness of a (weak) solution to (2.1) in the two-dimensional case or under certain geometrical assumptions on $\Omega_{\varepsilon}$ in general dimensions, see [26, 25] and the discussions therein. Note that if $\Omega_{\varepsilon}$ consists of a finite collection of bounded subdomains as we consider it in the numerical experiments, the unique continuation principle holds in general dimensions, see [25]. The solution depends continuously on the right-hand side in the sense of the stability estimate

[TABLE]

However, the dependencies of the constant $C_{\operatorname{stab}}(k,\varepsilon)$ on $\varepsilon$ and $k$ are not known explicitly in general, even for homogeneous media; see [39] for a positive and [6] for a negative result. Recently, various works have studied stability for the heterogeneous Helmholtz equation, see [10, 25, 43, 53]. We emphasize that the (common) assumptions on the coefficients are not even close to the target setting of this paper, but the recent work [35] indicates that polynomial in $k$ stability results hold for almost all $k$. The stability estimates readily implies an inf-sup condition of the form

[TABLE]

with $\gamma(\varepsilon,k)\approx(1+kC_{\operatorname{stab}}(\varepsilon,k))^{-1}$, see, e.g., [39].

2.2 Photonic crystals and homogenization results

To motivate some of our results and connect and differentiate them to the existing literature, we consider the case of a bounded scatterer with periodic micro-structure. To be more precise, we consider $\Omega_{\varepsilon}=\bigcup_{j\in\mathcal{I}\subset\mathbb{Z}^{2}}\varepsilon(j+\Sigma)$, where $\Sigma$ is an open Lipschitz subset of the unit square and $\mathcal{I}$ is an index set such that $\Omega_{\varepsilon}$ lies completely in a compact subset of $\Omega$. The last condition is necessary to have a uniformly bounded distance between $\Omega_{\varepsilon}$ and $\partial\Omega$.

Two-scale homogenization

In the perfectly periodic case and under the assumption that $\varepsilon$ is much smaller than the wave length $\lambda\sim k^{-1}$, homogenization using two-scale convergence [1, 36] gives information about the limit problem for fixed but arbitrary $k$ and $\varepsilon\to 0$. The surprising result is that one does not only obtain an effective diffusion matrix but also an effective wave number of the form $k^{2}\mu$, where $\mu$ is defined by an additional problem on $\Sigma$. The sign of $\mu$ depends on the choice of $k$, which highly influences the behavior of the macroscopic solution: For positive $\mu$, standard wave propagation is experienced whereas for negative $\mu$ evanescent waves form that decay exponentially inside the scatterer with the periodic microstructure. The above considerations hold for compactly embedded inclusions $\Sigma$, see [7], whereas a layered structure results in a degenerate diffusion matrix and induces so-called surface plasmons, see [8].

Moreover, due to the high contrast inclusions there is no weak $L^{2}$-convergence of the exact solution $u$ (for positive $\varepsilon$) to the macroscopic limit solution (named $u_{0}$), but instead an additional finescale-type contribution in the inclusion has to be added to $u_{0}$, see [1]. In contrast to problems without high contrast as in [14], we therefore cannot hope to find a good macroscopic $L^{2}$-approximation of $u$ – this is only possible in $\Omega_{1}$. However, with an additional correction, we can expect to recover the exact solution in a good $L^{2}$-sense, cf. [44] for corresponding numerical results. This explains parts of the a priori error estimates of our multiscale method in Section 4.2. Roughly speaking, we find a macroscopic approximation which is close to the exact solution in the $A$-weighted $L^{2}$ norm $\|\cdot\|_{0,A}$ (which, in particular, means a good approximation on $\Omega_{1}$.) Adding a suitable correction, we obtain an approximation which is good in the $\|\cdot\|_{1,A,k}$ sense. This in particular implies a good approximation in the standard $L^{2}$ norm.

Bloch wave homogenization and photonic crystals

If the wave length $\lambda\sim k^{-1}$ is of the same order as the periodicity length $\varepsilon$ and hence, also the diameter of one high-contrast inclusion $\varepsilon\Sigma$, other types of resonances occur which can no longer by predicted by two-scale homogenization theory, see [16]. Instead, Bloch wave homogenization methods are employed which are based on the Floquet-Bloch theorem, see [3, 15]. To get an idea of the (optical) properties of the finite material with periodic microstructure (also called a photonic crystal), one can consider its infinite counterpart on the whole space $\mathbb{R}^{d}$. If the wave number $k$ is an eigenvalue of the elliptic operator $-\nabla\cdot(A\nabla\cdot)$ on the whole space, the corresponding waves can propagate inside the (infinite) crystal. Due to the periodicity, the spectrum can be computed as the union of all spectra on one periodicity cell with varying quasi-periodic boundary conditions, the so-called Bloch spectra, see [15]. Band gaps correspond to “forbidden” wave numbers inducing evanescent waves as described above. A high contrast between the material properties allows to open a considerably large gap in the spectrum. However, when the periodic structure is destroyed by (local) defects, perturbation theory indicates that eigenvalues inside the band gaps may appear. The corresponding eigenfunctions are often localized near the defect, see [34] for a general discussion and [4] for the specific example of a line defect. As mentioned, this considers the infinite crystal with the hope that a considerably large finite photonic crystal will have similar properties. The multiscale method presented in this contribution does not rely on two-scale homogenization results and therefore, can also cope with the regime of Bloch wave homogenization. Hence, we can directly study finite photonic crystals with possible local defects and their properties.

3 Multiscale method

In this section, we introduce the multiscale method based on the ideas of the Localized Orthogonal Decomposition (LOD) [38, 32]. First in Section 3.1, we need some more notation on meshes and in particular a suitable interpolation operator which pays attention to the high contrast. After these preliminaries, we can define our method in detail in Section 3.2. We close with some remarks concerning implementation in Section 3.3.

3.1 Meshes, finite element spaces, and interpolation operators

We cover $\Omega$ with a regular mesh $\mathcal{T}_{H}$ consisting either of simplices or of parallelograms/parallelepipeds. The mesh is assumed to be shape regular in the sense that the aspect ratio of the elements of $\mathcal{T}_{H}$ is bounded uniformly from below. We introduce the mesh size $H=\max_{T\in\mathcal{T}_{H}}\operatorname{diam}T$. Although we will assume $H=O(\varepsilon)$ for most of this article, see (4.2), note that $\mathcal{T}_{H}$ does not necessarily resolve the interfaces between $\Omega_{1}$ and $\Omega_{\varepsilon}$. Associated with an element $T\in\mathcal{T}_{H}$ we define its neighborhood as

[TABLE]

Thereby, for any $m\in\mathbb{N}_{0}$, the $m$-layer patches are defined inductively via $\textup{N}^{m+1}(T)=\textup{N}(\textup{N}^{m}(T))$ with $\textup{N}^{0}(T):=T$. The shape regularity implies that there is a bound $C_{\operatorname{ol},m}$ (depending only on $m$) of the number of the elements in the $m$-layer patch, i.e.,

[TABLE]

Throughout this article, we assume that $\mathcal{T}_{H}$ is quasi-uniform, which implies that $C_{\operatorname{ol},m}$ grows at most polynomially with $m$. We discretize the space $H^{1}_{0}(\Omega)$ with the lowest order Lagrange elements over $\mathcal{T}_{H}$, and denote this space by $V_{H}$. This means that $V_{H}=H^{1}(\Omega)\cap\mathcal{S}^{1}(\mathcal{T}_{H})$, where $\mathcal{S}^{1}(\mathcal{T}_{H})$ denotes the space of globally continuous functions which are polynomials of partial degree $\leq 1$ (for quadrilateral elements) or polynomials of total degree $\leq 1$ (for simplicial elements). In the case of quadrilateral meshes, one can easily exploit a possible (periodic) structure or pattern in the coefficient $A$, see Section 3.3.

Let $I_{H}:H^{1}(\Omega)\to V_{H}$ denote a bounded local linear projection operator, i.e., $I_{H}\circ I_{H}=I_{H}$, with the following stability and approximation properties for all $v\in H^{1}(\Omega)$ and all $w\in\ker I_{H}$

[TABLE]

where we recall that the constant hidden in $\lesssim$ is independent of $H$ and $\varepsilon$. We emphasize that we assume stability and approximation properties of $I_{H}$ in $A$-weighted norms, which is the crucial difference between the low-contrast and the high-contrast case. The stability in the energy norm (3.4) can be deduced from (3.2) and (3.3) at least in the following two cases: The mesh satisfies the resolution condition (4.2) (as assumed later on), which restricts the mesh size significantly and may be too pessimistic. The other possibility is that for each $T\in\mathcal{T}_{H}$, the quotient $\frac{|T|}{|T\cap\Omega_{1}|}$ is uniformly bounded from above, which mainly implies that $H$ may not be too fine. Anyhow, the crucial assumptions are (3.2) and (3.3). They have been verified under certain geometric assumptions using special $A$-weighted interpolants in [29, 51]. A possible choice of [51] that we also use in our numerical implementation builds upon $A$-weighted local $L^{2}$ projections in the following sense. Let $\mathcal{N}_{\operatorname{int}}$ denote the interior nodes of $\mathcal{T}_{H}$ and $\lambda_{z}$ the hat function associated with the node $z$. We then define

[TABLE]

with the local $A$-weighted $L^{2}$-projection $P_{A,\omega_{z}}$ on the patch $\omega_{z}:=\{K\in\mathcal{T}_{H}:z\in K\}$. More precisely, for any $v\in H^{1}(\omega_{z})$, $P_{A,\omega_{z}}(v)\in V_{H}|_{\omega_{z}}$ is given by

[TABLE]

Assumptions (3.2) and (3.3) hold for this interpolant if the coefficient $A$ is quasi-monotone. For instance, in the two-dimensional case, this condition is satisfied if – roughly speaking – $\Omega_{\varepsilon}$ consists of small inclusions of diameter $\varepsilon$ that do not touch each other. This setup is considered in our numerical experiments. For more details on quasi-monotonicity we refer to [51] and the originial works in the context of domain decomposition [17, 46]. In [29], a Scott-Zhang-type construction is employed, where the geometrical assumptions on $\Omega_{\varepsilon}$ are not so easy to describe. We emphasize, however, that the periodic structure underlying our numerical experiments in Section 6 is mentioned as an example in [29] (with circle inclusion instead of square inclusions) and is also analyzed in [11]. For this choice, our numerical experiments indicate that for a high contrast, the above described $A$-weighted interpolation operator $I_{H}$ (cf. [51]) is favorable over its unweighted variant

[TABLE]

which utilizes local projections w.r.t. the standard $L^{2}$-norm. Note that this unweighted construction is used for low-contrast problems, e.g. in [49]. Moreover, we stress that especially the construction in [51] can be extended to coefficients which take more than two distinct values (see the discussion in [51]).

3.2 Definition of the method

We formulate a (generalized) Petrov-Galerkin method to discretize (2.1). The method uses special multiscale trial and test functions, see [50]. The idea is to interpret the kernel of the interpolation operator $I_{H}$ as the space of functions with finescale heterogeneities and to incorporate parts of these information by special localized finescale problems.

Let $W(\textup{N}^{m}(T)):=\{w\in\ker I_{H}\;,\,w=0\text{ in }\Omega\setminus\textup{N}^{m}(T)\}$. Note that $\ker I_{H}$ is infinite-dimensional and hence $W(\textup{N}^{m}(T))$ is non-trivial. We introduce the element-wise localized corrector approximation $\mathcal{Q}_{T,m}:V_{H}\to W(\textup{N}^{m}(T))$ via

[TABLE]

Here, $\mathcal{B}_{\omega}$ denotes the restriction of $\mathcal{B}$ to a subset $\omega\subset\Omega$. Note that (3.7) requires only the local computation of finescale problems if $m$ is small. The parameter $m$ is commonly called oversampling parameter and will be coupled to the mesh size $H$ later on. Using the element-wise correctors $\mathcal{Q}_{T,m}$ we define the localized corrector approximation via

[TABLE]

The dual localized corrector is $\mathcal{Q}^{*}_{m}v_{H}:=\overline{\mathcal{Q}_{m}(\overline{v_{H}})}$.

The localized method now uses the pair $((\operatorname{id}+\mathcal{Q}_{m})V_{H},(\operatorname{id}+\mathcal{Q}^{*}_{m})V_{H})$ as the new trial and test spaces in the Petrov-Galerkin formulation: Find $u_{H,m}\in V_{H}$ as solution to

[TABLE]

Note that $u_{H,m}$ lies in the standard finite element space and, hence, is called a macroscopic approximation since it cannot contain finescale information of the exact solution $u$. This information is only recovered in

[TABLE]

The well-posedness of the scheme and the error between the exact and the numerical solution will be analyzed in the next sections.

3.3 Practical aspects

The method presented in the previous section is not yet ready to use since the corrector problems (3.7) are still infinite-dimensional. Moreover, the number of these problems to solve increases with smaller mesh size $H$ so that one has to think about an efficient computation. Both issues are briefly addressed in this section. For general aspects of the implementation, we refer to [20].

Fully discrete version of the multiscale method

To discretize the corrector problems (3.7), we introduce a finescale shape-regular quadrilateral or simplicial mesh $\mathcal{T}_{h}$ of $\Omega$, which resolves all features and discontinuities of $A$. Moreover, we assume that the mesh is fine enough that a direct finite element simulation on $\mathcal{T}_{h}$ would yield a faithful reference solution $u_{h}$ to $u$. We stress, however, that this reference solution is not needed and not computed in our method. Only the local corrector problems are solved on the finescale mesh using the space $W_{h}(\textup{N}^{m}(T)):=\{w\in H^{1}(\Omega)\cap\mathcal{S}^{1}(\mathcal{T}_{h})\;,\,w=0\text{ in }\Omega\setminus\textup{N}^{m}(T)\text{ and }I_{H}w=0\}$. The assumption that $\mathcal{T}_{h}$ is sufficiently fine (as described above) is required for two reasons: (i) It guarantees that the multiscale method (3.9) still possesses a unique solution because for a sufficiently fine mesh, $\mathcal{B}$ fulfills an inf-sup condition over $H^{1}(\Omega)\cap\mathcal{S}^{1}(\mathcal{T}_{h})$ (which is a central argument in the proof of Theorem 4.4 below). (ii) It implies that the error committed by this additional discretization stays negligible. Due to the latter fact, we only analyze (for simplicity) the semi-discrete method as formulated in Section 3.2 and refer to, e.g., [23] for details on the minor technical changes in the proofs for the fully discrete method.

Efficient computation of the correctors

First, we note that all corrector problems are independent and hence, can be easily computed in parallel. What is even more important is that additional structure in the coarse mesh $\mathcal{T}_{H}$ and the coefficient can be exploited to reduce the number of corrector problems. For instance, in the periodic setting the corrector problems are translation invariant if we use a quadratic mesh with a mesh size which is an integer multiple of the periodicity length. Apart from a few problems at the boundary, only one interior corrector problem needs to be solved, cf. [23] for the homogeneous Helmholtz equation. Similar arguments allow to reduce the number of corrector problems for a (periodic) photonic crystal with a few local defects. If one has already computed the correctors for the periodic photonic crystals and then introduces local defects (as perturbations of the periodic structure), one can employ local error indicators as in [30] to automatically detect locations where the correctors need to be recomputed. Numerical experiments in [28] have shown that thereby a large part of the old correctors can be reused for local defects, which also lowers the computational complexity. Here, however, we do not consider this approach in more detail since we assume the geometry of $\Omega_{\varepsilon}$ to be a priori given and not as a deviation from a previous, more structured, set-up.

Efficient computation of the LOD solution

The efficient computation of the LOD solution as well as CPU timings have been addressed in full detail in [20] for diffusion and eigenvalue problems. The algorithms turn over to the Helmholtz case in a very natural way. Without any assumptions on the structure of the coefficient, we need to solve $O(H^{-d})$ cell problems of type (3.8). If these are discretized on a mesh of size $h$ (see above), the number of degrees of freedom for each problem is $O(m^{d}(H/h)^{d})$. Given their coercivity and since $H$ resolves both the wave number and the oscillations of $A$, they can be solved efficiently. As pointed out in the previous paragraph, the number of cell problems can be drastically reduced for structured coefficients. The linear system of the LOD associated with (3.9) is of small dimension, namely of $\operatorname{dim}(V_{H})$. Note that numerical experiments show that $m=2,3$ is sufficient and therefore the sparsity pattern of the matrix is only slightly enlarged by a factor $O(m^{d})$ in comparison to standard FEM. The conditioning is comparable to standard FEM on the coarse mesh $\mathcal{T}_{H}$. The most efficient iterative solution of such systems is an open question in numerical analysis and beyond the scope of this article. The present approach does neither simplify nor overcomplicate this issue of numerical linear algebra. However, the tailored discretization allows the reduction of the degrees of freedom and hence the size of the system to almost optimal in the sense of sampling theory, whereas standard FEM would require much higher costs even in the regime of constant coefficients.

4 Error analysis

In this section, we analyze the multiscale method of Section 3.2 for the high-contrast Helmholtz problem. We first show in Section 4.1 that under a reasonable resolution condition $kH\varepsilon^{-1}\lesssim 1$, the corrector problems are coercive and therefore well-posed. Then, we give the main results of this paper in Section 4.2, namely the inf-sup condition of the method as well as several a priori error estimates. These results are derived under the above resolution condition and the oversampling condition $m\gtrsim|\log(kH\varepsilon^{-1})|$ if the original model problem satisfies a polynomial $k$- and $\varepsilon$-dependence of the stability constant $C_{\operatorname{stab}}(k,\varepsilon)$. The resolution condition $kH\varepsilon^{-1}\lesssim 1$ should be compared to the usual resolution condition encountered for the quasi-optimality in standard finite element methods, which presumably is $k^{2}h\varepsilon^{-2}\lesssim 1$ at best in non-smooth scenarios. Note that the condition for existence of a FEM solution can be relaxed to $h(k\varepsilon^{-1})^{\alpha}\lesssim 1$, but one will always have $\alpha>1$ for standard FEM.

4.1 The corrector problems

We first consider an idealized variant of (3.9) with “$m=\infty$”, i.e., with $\textup{N}^{m}(T)=\Omega$ in the corrector problems (3.7). For notational convenience, we denote by $W$ the kernel of the interpolation operator $I_{H}$. Define the ideal correction operator $\mathcal{Q}:H^{1}(\Omega)\to W$ via

[TABLE]

Note that $\mathcal{Q}=\sum_{T\in\mathcal{T}_{H}}\mathcal{Q}_{T}$, where $\mathcal{Q}_{T}$ solves (4.1) with right-hand side $\mathcal{B}_{T}$. Since the sesquilinear form $\mathcal{B}$ is indefinite, we first have to show the well-posedness of these corrector problems.

Proposition 4.1**.**

If the mesh size $H$, the wave number $k$, and diffusion parameter $\varepsilon$ satisfy the condition

[TABLE]

we have the following equivalence of norms on $W$

[TABLE]

Moreover, the sesquilinear form $\mathcal{B}$ is elliptic over $W$, i.e., there exists $\alpha>0$ (independent of $H$, $\varepsilon$, and $k$) such that

[TABLE]

Proof.

Similar to [50], the crucial observation is that (3.3) implies

[TABLE]

for all $w\in W$. Thus, the energy norm $\|\cdot\|_{1,A,k}$ is equivalent to the weighted $H^{1}$ semi norm $|\cdot|_{1,A}$ if (4.2) is satisfied. Moreover, the $L^{2}$ part of the sesquilinear form $\mathcal{B}$ can be absorbed in the gradient part if (4.2) is satisfied, leading to the ellipticity. ∎

The ellipticity of $\mathcal{B}$ over $W$ implies the unique solvability of (4.1) due to the theory Lax-Milgram-Babuška. The same argument can be applied to obtain the well-posedness of the localized corrector problems (3.8).

Remark 4.2**.**

The resolution condition (4.2) is a natural consequence and generalization of the resolution condition $kH\lesssim 1$ for the homogeneous or low-contrast Helmholtz equation, cf. [23, 50]. In fact, on $\Omega_{\varepsilon}$, $k\varepsilon^{-1}$ is the effective wave number which needs to be resolved. There are certainly geometries of $\Omega_{\varepsilon}$ where this condition is necessary and sharp, but it is in general not clear how “likely” such situations are. As discussed, at least conditions like $\varepsilon^{-2}k^{2}h\lesssim 1$ and $h\ll\varepsilon$ are expected for general non-smooth coefficients for standard finite element methods where possible advantages of high-order methods [41, 42] cannot easily be exploited. (An exact statement or estimate of this condition is not available for this setting as the dependence of the stability constant on $\varepsilon$ and $k$ is not known.)

The error between the idealized corrector and its localized (or truncated) approximation decays exponentially in the number of layers $m$ (also called oversampling parameter) in the following way.

Theorem 4.3**.**

Let $\mathcal{Q}$ be defined by (4.1) and $\mathcal{Q}_{m}$ be defined by (3.8). Let (4.2) be satisfied. Then there exists a constant $0<\beta<1$, independent of $H$, $m$, $k$, and $\varepsilon$, such that for all $v_{H}\in V_{H}$ it holds that

[TABLE]

The proof is given in Section 5.1. Employing that both correctors map into the kernel space and using (3.3), we deduce furthermore that

[TABLE]

We emphasize that neither the constant hidden in $\lesssim$ nor $\beta$ depend on $H$, $m$, $k$, or $\varepsilon$. In particular the independence of $\beta$ from $\varepsilon$ is the so-called contrast-independent localization. It can only be obtained with the $A$-weighted interpolation operator and is the major difference of this work from the previous ones on the low-contrast Helmholtz equation [10, 23, 50].

4.2 Stability of the method and a priori error estimates

The previous section showed that the corrector problems have a unique (and stable) solution. Before analyzing the error of the multiscale method (3.9), we need to show the stability (and well-posedness) of the method by proving an inf-sup condition. Since the arguments utilize the corrector error estimate (4.4), we first examine what happens if we replace $\mathcal{Q}_{m}$ and $\mathcal{Q}^{*}_{m}$ in (3.9) by their idealized counterparts $\mathcal{Q}$ from (4.1) and its adjoint $\mathcal{Q}^{*}$. Hence, $\mathcal{Q}^{*}:H^{1}(\Omega)\to W$ is defined via

[TABLE]

and we note that $\mathcal{Q}^{*}v:=\overline{\mathcal{Q}(\overline{v})}$. By the definition of $\mathcal{Q}^{*}$ and $W$, it holds that $W$ and $(\operatorname{id}+\mathcal{Q}^{*})V_{H}$ are orthogonal with respect to $\mathcal{B}$. More precisely, we have $\mathcal{B}(w,(\operatorname{id}+\mathcal{Q}^{*})v_{H})=0$ for all $w\in W$ and all $v_{H}\in V_{H}$. Note that the same holds true for $(\operatorname{id}+\mathcal{Q})V_{H}$ and $W$. Next we observe that the interpolation $I_{H}u$ of the exact solution satisfies

[TABLE]

Using once more the $\mathcal{B}$-orthogonality of $W$ and $(\operatorname{id}+\mathcal{Q}^{*})V_{H}$, this implies that $u_{\operatorname{LOD}}:=(\operatorname{id}+\mathcal{Q})I_{H}u$ solves

[TABLE]

The well-posedness and the a priori estimates now follow from the properties of the above idealized method and the fact that the ideal and the localized correctors are exponentially close.

Theorem 4.4**.**

Under the resolution condition (4.2) and the oversampling condition

[TABLE]

$\mathcal{B}$* satisfies the following inf-sup condition: there exists $\gamma_{\operatorname{LOD}}\approx\gamma(\varepsilon,k)>0$ such that*

[TABLE]

The proof is similar to [50] and needs the contrast-independent exponential decay of the corrector error, for details see Section 5.2. The oversampling condition depends on the stability constant for the original problem so that the dependence of $m$ on $\varepsilon$ and $k$ is not analytically known. If we assume $\gamma(\varepsilon,k)\approx k^{-q}\varepsilon^{p}$ for non-negative numbers $p$ and $q$ independent of $k$ and $\varepsilon$, the oversampling condition roughly reads $m\gtrsim|\log(k\varepsilon^{-1})|$ (with a multiplicative factor $pq$ hidden in the $\lesssim$ notation). This logarithmic dependence on the effective wave number $k\varepsilon^{-1}$ in $\Omega_{\varepsilon}$ is in agreement with the results for the homogeneous or low-contrast Helmholtz equation in [10, 23, 50]. Numerical experiments moreover indicate that rather small numbers $m=2,3$ are sufficient.

Theorem 4.5**.**

Let $u$ be the solution to (2.1) and let $u_{H,m}$ be the solution to (3.9) and set $u_{\operatorname{LOD},m}:=(\operatorname{id}+\mathcal{Q}_{m})u_{H,m}$. Assume that (4.2) and the oversampling condition (4.5) are satisfied. Then it holds that

[TABLE]

The proof is detailed in Section 5.2.

If $m$ does not only fulfill the oversampling condition (4.5), but is also coupled to $H$ like $m\approx|\log(\gamma^{-1}(\varepsilon,k)C_{\operatorname{ol},m}^{1/2}H)|$, Theorem 4.5 simplifies to the following convergence orders: (i) The error between the exact and the full multiscale solution converges linearly in the energy norm and quadratically in an $A$-weighted $L^{2}$ norm, cf. (4.7) and (4.8). (ii) The error between the exact solution and the FE part of the multiscale solution converges (at least) linearly in the $A$-weighted $L^{2}$ norm, cf. (4.9). For the latter, note that $\inf_{v_{H}\in V_{H}}|u-v_{H}|_{1,A}\lesssim|u|_{1,A}\lesssim C_{\operatorname{stab}}(\varepsilon,k)\|f\|$ independent of the regularity of $u$. If $u$ has more than $H^{1}(\Omega)$ regularity, the rate in (4.9) may improve. We stress that the remaining terms on the right-hand side of (4.9) are of order $H^{2}$ if the coupling $m\approx|\log(\gamma^{-1}(\varepsilon,k)C_{\operatorname{ol},m}^{1/2}H)|$ is satisfied.

These estimates generalize the (expected) approximation results formulated in Section 2.2. We can only find a macroscopic approximation of the exact solution in $\Omega_{1}$, which is reflected in (4.9) by the $A$-weighting in the $L^{2}$ norm. For a good approximation in the energy norm, which also includes the standard $L^{2}$ norm, we need an additional (finescale) corrector, which is also present in (4.7). In the main error estimate (4.7), we furthermore emphasize that the error from the volume term is weighted by $A^{-1/2}$, which induces a factor of $\varepsilon^{-1}$ on the right-hand side of the estimate if the support of $f$ intersects with $\Omega_{\varepsilon}$.

Remark 4.6**.**

Our main results of this paper rely on the $\varepsilon$-independent trace inequality mentioned in Section 2.1, which needs that $\Omega_{\varepsilon}$ is bounded uniformly away from $\partial\Omega$. For settings where this is not the case, one can carefully trace the occurrence of extra $\varepsilon$ powers because of the continuity of $\mathcal{B}$ and we find that this mainly changes the oversampling condition to $m\gtrsim|\log(\gamma^{-1}(\varepsilon,k)\varepsilon C_{\operatorname{ol},m}^{1/2})|/|\log(\beta)|$. This, however, is not critical since we expect a dependence of $m$ on $\varepsilon$ anyway because of $\gamma(\varepsilon,k)$.

Remark 4.7**.**

The proofs of Theorems 4.4 and 4.5 reveal that the $\varepsilon$ dependence in the resolution and the oversampling condition come from the $\varepsilon$ dependence of the lower bound on $A$ (and not the minimal diameter $\varepsilon$ of an individual “valley” of $\Omega_{\varepsilon}$). In fact, one can directly generalize the results to the case where $A$ jumps between the values $1$ and $\alpha^{2}\ll 1$ by replacing all occurrences of $\varepsilon$ by $\alpha$. In particular, this means that our results remain valid for other scalings of $A$ beyond the (physically interesting) one of, for instance, [7].

In the “complementary” high contrast setting, where $A$ jumps between values $1$ and $\alpha\gg 1$, one also needs $A$-weighted interpolation operators to get a contrast-independent decay of the correctors and thereby contrast-independent error estimates for the LOD. Since in this case it trivially holds $\|v\|\lesssim\|v\|_{0,A}$ for any $v\in H^{1}(\Omega)$, the resolution condition (4.2) reduces to the usual $kH\lesssim 1$ (independent of $\alpha$), cf. for instance the proof of Proposition 4.1. For the same reason, one can replace $\|A^{-1/2}f\|$ by $\|f\|$ in (4.7) and (4.9). These improvements of the results reflect and underline the fact that this scaling of $A$ is “easier” as it does not lead to high (effective) wave numbers as discussed in the introduction.

5 Proofs of the main results

5.1 Proof of the corrector error

The goal of this section is to prove (4.4). The main ingredient is the following exponential decay result. We empasize once more that the following decay is contrast-independent, i.e., the rate $\tilde{\beta}$ and also the multiplicative constant hidden in the notation $\lesssim$ both do not depend on $\varepsilon$. This is achieved by using an $A$-weighted interpolation operator and distinguishes the present analysis from previous works on the Helmholtz equation [10, 23, 50].

Proposition 5.1**.**

Let $\mathcal{Q}=\sum_{T\in\mathcal{T}_{H}}\mathcal{Q}_{T}$ be defined by (4.1) and $v_{H}\in V_{H}$. There exists a constant $0<\tilde{\beta}<1$, independent of $H$, $m$, and $\varepsilon$ such that

[TABLE]

Proof.

We define the cutoff function $\eta\in V_{H}$ via

[TABLE]

and set $\mathcal{R}=\operatorname{supp}(\nabla\eta)$. Let $v_{H}\in V_{H}$ and denote $\phi:=\mathcal{Q}_{T}v_{H}$. Elementary estimates lead to

[TABLE]

with

[TABLE]

Since $w:=(\operatorname{id}-I_{H})(\eta\phi)\in W$, the idealized corrector problem (4.1) and the fact that $w$ has only support outside $T$ imply that $\mathcal{B}(\phi,w)=\mathcal{B}_{T}(v_{H},w)=0$. Therefore, we obtain

[TABLE]

Hence, the stability and approximation estimates (3.2) and (3.3) for $I_{H}$, the properties of $\eta$ , and the resolution condition (4.2) give

[TABLE]

where the first term can be hidden on the left-hand side.

Because of $\operatorname{supp}(I_{H}(\eta\phi))\subset\textup{N}(\mathcal{R})$, the properties of $I_{H}$ and the above estimate for $|\eta\phi|_{1,A}$ lead to

[TABLE]

Finally, the properties of $\eta$ and the approximation result (3.3) for $I_{H}$ show that

[TABLE]

All in all, this gives

[TABLE]

with a constant $C$ independent of $H$, $m$, $k$, and $\varepsilon$. We recall that $\textup{N}^{2}(\mathcal{R})=\textup{N}^{m}(T)\setminus\textup{N}^{m-5}(T)$. Because of

[TABLE]

we obtain

[TABLE]

Note that $(1+2C^{-1})^{-1}<1$, so that a repeated application of the above argument and algebraic manipulations finish the proof. Note in particular that $\tilde{\beta}$ is independent of $H$, $m$, $k$, and $\varepsilon$ because $C$ is independent of these quantities. ∎

The idea for the proof of Theorem 4.3 is now that the ideal corrector is only truncated in $\Omega\setminus\textup{N}^{m}(T)$, where we know by the previous Proposition that the contribution is exponentially small.

Proof of Theorem 4.3.

We define the cutoff function $\eta\in\mathcal{S}^{1}(\mathcal{T}_{H})$ via

[TABLE]

With the ellipticity and continuity of $\mathcal{B}$ over $W$, we deduce Céa’s Lemma

[TABLE]

Inserting $w_{T,m}=(\operatorname{id}-I_{H})(\eta\mathcal{Q}_{T}v_{H})$ and using stability and approximation properties of $I_{H}$ yields – similar to the proof of Proposition 5.1 –

[TABLE]

Applying (5.1) results in an exponential decay result for the error between the element correctors $\mathcal{Q}_{T}$ and $\mathcal{Q}_{T,m}$.

To obtain the global estimate, we set $z:=(\mathcal{Q}-\mathcal{Q}_{m})v_{H}$ and $z_{T}:=(\mathcal{Q}_{T}-\mathcal{Q}_{T,m})v_{H}$ and define the cutoff function $\eta\in\mathcal{S}^{1}(\mathcal{T}_{H})$ via

[TABLE]

The ellipticity of $\mathcal{B}$ yields

[TABLE]

The second term vanishes because $(\operatorname{id}-I_{H})(\eta z)\in W$ with support outside $\textup{N}^{m}(T)$. The function $(1-\eta)z$ vanishes on $S:=\{\eta=1\}$ and we obtain with the scaling properties of $\eta$, (3.3), (3.2), and the resolution condition (4.2) that

[TABLE]

$I_{H}(\eta z)$ vanishes on $\Omega\setminus\textup{N}(\operatorname{supp}(1-\eta))$. Hence, we infer with the properties of $\eta$, (3.3), (3.2), and the resolution condition (4.2) that

[TABLE]

All in all, summation over all $T$ and Cauchy inequality yield with the finite overlap of patches

[TABLE]

Combination with (5.2) finishes the proof. ∎

5.2 Proof of the well-posedness of the method and of the error estimates

We first prove the discrete inf-sup condition which implies that the solution of the multiscale method is indeed well-defined.

Proof of Theorem 4.4.

Let $v_{H}\in V_{H}$. From the inf-sup condition for the model problem we infer that there exists $\psi\in H^{1}(\Omega)$ with $\|\psi\|_{1,A,k}=1$ such that

[TABLE]

Define now $\psi_{H}:=I_{H}\psi$. We have

[TABLE]

Since $\mathcal{Q}^{*}$ is a projection onto $W$, we have $(\operatorname{id}+\mathcal{Q}^{*})I_{H}\psi=(\operatorname{id}+\mathcal{Q}^{*})\psi$. The solution properties of $\mathcal{Q}^{*}$ and (4.1) of $\mathcal{Q}$ imply $\mathcal{B}((\operatorname{id}+\mathcal{Q})v_{H},(\operatorname{id}+\mathcal{Q}^{*})\psi_{H})=\mathcal{B}((\operatorname{id}+\mathcal{Q})v_{H},\psi)$. Hence, the continuity of $\mathcal{B}$ and (4.4) imply

[TABLE]

The stability of $I_{H}$ yields $\|\psi_{H}\|_{1,A,k}\lesssim\|\psi\|_{1,A,k}=1$ and

[TABLE]

Moreover, we have the following stability for $\mathcal{Q}_{m}$:

[TABLE]

This finally gives

[TABLE]

which together with the oversampling condition (4.5) finishes the proof. ∎

We can now prove the a priori error estimates which shed alight onto the approximation properties of the multiscale method.

Proof of Theorem 4.5.

Proof of (4.7): Denote $e:=u-u_{\operatorname{LOD},m}=u-(\operatorname{id}+\mathcal{Q}_{m})u_{H,m}$ and define $e_{H,m}:=(\operatorname{id}+\mathcal{Q}_{m})I_{H}e=(\operatorname{id}+\mathcal{Q}_{m})(I_{H}u-u_{H,m})$. The triangle inequality gives

[TABLE]

We will show that under the oversampling condition (4.5) the second term can be bounded by the first term.

Let $z_{H}\in V_{H}$ be the solution of the adjoint problem

[TABLE]

Then we obtain with the orthogonality of $W$ and $(\operatorname{id}+\mathcal{Q})V_{H}$ and the Galerkin orthogonality $\mathcal{B}(e,(\operatorname{id}+\mathcal{Q}_{m}^{*})\psi_{H})=0$ for all $\psi_{H}\in V_{H}$ that

[TABLE]

The solution $z_{H}$ of the adjoint problem fulfills the following stability

[TABLE]

Using this stability and the corrector error (4.4) results in

[TABLE]

Dividing by $\|e_{H,m}\|_{1,A,k}$ on both sides and using the oversampling condition (4.5) gives the bound $\|e_{H,m}\|_{1,Ak}\lesssim\|e-e_{H,m}\|_{1,A,k}$. Second, we estimate the error $e-e_{H,m}=u-(\operatorname{id}+\mathcal{Q}_{m})I_{H}u\in W$. With the ellipticity of $\mathcal{B}$ over $W$, the $\mathcal{B}$-orthogonality between $W$ and $(\operatorname{id}+\mathcal{Q})V_{H}$ and Galerkin orthogonality we obtain

[TABLE]

The last term can be estimated employing (4.4) as

[TABLE]

where we also used the continuity of $\mathcal{B}$ and the stability of $I_{H}$. $\|u\|_{1,A,k}$ can be estimated against the data with the inf-sup constant of the model problem. For the first term, we can insert $I_{H}(e-e_{H,m})$ and use the approximation properties of $I_{H}$ to obtain

[TABLE]

Proof of (4.8): We again abbreviate $e:=u-u_{\operatorname{LOD},m}$ and consider the following two adjoint problems: (i) Find $z\in H^{1}(\Omega)$ such that

[TABLE]

(ii) Find $z_{H,m}\in V_{H}$ such that

[TABLE]

Set $z_{\operatorname{LOD},m}:=(\operatorname{id}+\mathcal{Q}_{m}^{*})z_{H,m}$. Then we deduce with Galerkin orthogonality

[TABLE]

Note that $z$ and $z_{H,m}$ are solutions to adjoint problems with a volume term on the right-hand side of the form $Ae$. Hence, we deduce similar to the above estimate for $\|e\|_{1,A,k}$ that

[TABLE]

Combination of the foregoing estimates finishes the proof.

Proof of (4.9): Let $u_{\operatorname{LOD}}$ be the solution to the idealized variant of (3.9), i.e., with $\textup{N}^{m}(T)=\Omega$ for all elements $T$, as introduced at the beginning of Section 4.2. Because of (3.3) and the definition of the norms in (2.2) and (2.4), we obtain

[TABLE]

The projection property of $I_{H}$ and its stability (3.2) imply for the first term

[TABLE]

To estimate the second term, we abbreviate $e_{H}:=I_{H}u-u_{H,m}\in V_{H}$ and note that by the definition of $\mathcal{Q}$ and the stability (3.4), we have

[TABLE]

Because of the discrete inf-sup condition (4.6) (which also holds for the idealized version of the LOD), there exist $z_{H}\in V_{H}$ with $\|z_{H}\|_{1,A,k}=1$ such that

[TABLE]

Proceeding in a similar way as for the estimate of $e_{H,m}$ above, we obtain with the definition of $\mathcal{Q}$ and Galerkin orthogonality that

[TABLE]

where we used (4.4) in the last estimate. The remaining error $\|u-u_{\operatorname{LOD},m}\|_{1,A,k}$ has been estimated in (4.7). ∎

6 Numerical experiments

We investigate the method and the a priori estimates in several experiments. The convergence history plots display the relative error in the indicated norm versus the (coarse) mesh size $H$. We consider the LOD solution $u_{\operatorname{LOD},m}$ defined in (3.10) as well as its FE part $u_{H,m}$, cf. (3.9). For comparison, we also compute a standard FE solution (denoted as P1FEM in the convergence plots) on $V_{H}$ as well as the best approximation in $V_{H}$ with respect to the indicated norm (denoted as P1-best in the convergence plots). The solution plots display the real part with a color map truncated to the interval $[-2,2]$ to visualize the wave behavior outside of the scatterer, which is our main interest. We use a uniform simplicial mesh on the unit square $\Omega=(0,1)^{2}$ in all our experiments and compute all correctors, i.e., we do not exploit the periodic structure as discussed in Section 3.3. As right-hand side we use

[TABLE]

with varying center $x_{0}$. We compare results for the $A$-weighted interpolation operator $I_{H}$ (3.5) as well as its unweighted variant $I_{H}^{1}$ (3.6).

6.1 Periodic structure with Mie resonances

We consider a periodic scatterer in the domain $(0.25,0.75)^{2}$ with

[TABLE]

We choose the center for $f$ as $x_{0}=(0.125,0.5)$ and the wave number $k=9$. This choice of $k$ is connected to a negative-valued effective wave number $k^{2}\mu$ in homogenization theory and this unusual behavior is induced by (Mie) resonances in the small inclusions, cf. [44]. Therefore, we expect that an approximation of the exact solution might be hard in that case and consider it as a good test for the multiscale method. Note that this setting fulfills our assumptions, in particular $\Omega_{\varepsilon}$ is uniformly bounded away from $\partial\Omega$. A reference solution (still denoted $u$) is computed with standard finite elements with mesh size $h=2^{-9}$ which is also the fine mesh used for the discretization of the corrector problems as explained in Section 3.3. The convergence history in $H$ for the errors of $u_{\operatorname{LOD},m}$ and $u_{H,m}$ with respect to the reference solution $u$ is studied for patch sizes $m=1,2,3$.

For $\varepsilon=2^{-3}$, the results obtained with $I_{H}$ and $I_{H}^{1}$ are very similar so that we restrict the confirmation of the convergence rates to $I_{H}^{1}$. Figure 6.1 shows the convergence histories for $u-u_{\operatorname{LOD},m}$ and we verify the linear rate in the energy norm (and the standard $L^{2}$ norm) as well as the quadratic rate in the weighted $L^{2}$ norm. A patch size of $m=2,3$ seems to be sufficient. We observe a clear superiority of the multiscale method over a standard finite element method on the coarse mesh, which fails to yield good approximations for these values of $H$. The standard FEM even deviates significantly from the FE best approximations with respect to the energy norm and weighted $L^{2}$-norm (computed by projecting the reference solution onto $V_{H}$). The FE best approximations are outperformed by the multiscale method which takes fine-scale features into account.

In this context it is very interesting to analyze the error $u-u_{H,m}$ of the FE part $u_{H,m}$ of the multiscale method in Figure 6.2. In the standard as well as the weighted $L^{2}$ norm, we observe a convergence which closely follows the behavior of the FE best approximations in the space $V_{H}$ with respect to the corresponding norms. This clearly underlines that there exists a macroscopic approximation to the exact solution which is good in an $L^{2}$ sense. While the multiscale method is able to (almost) find this best approximation, the FEM (using the same approximation space) fails completely. Note that we observe a numerical convergence rate of about $2$ for the error $u-u_{H,m}$ in this example. This is explained by higher regularity of the exact solution, which results in an order $H$ for the best-approximation error $\inf_{v_{H}\in V_{H}}|u-v_{H}|_{1,A,k}$ as discussed after Theorem 4.5.

For $\varepsilon=2^{-4}$, we compare $\|u-u_{\operatorname{LOD},m}\|_{1,A,k}$ and $\|u-u_{H,m}\|_{0,A}$ for the unweighted interpolation operator $I_{H}^{1}$ (3.6) and its $A$-weighted variant $I_{H}$ (3.5). As Figure 6.3 shows, only the weighted interpolant succeeds to yield good approximations (for patch size $m=3$ at least). In the $\|\cdot\|_{0,A}$-norm, the error even follows the best approximation error for the weighted operator (bottom right), which is the best we can hope for. We moreover observe a larger pre-asymptotic range for $\varepsilon=2^{-4}$ than before for $\varepsilon=2^{-3}$ which indicates that the resolution condition between $k$, $\varepsilon$ and $H$ is sharp in this setting.

6.2 Periodic structure with local defects

We now place the periodic structure described in the previous section into a slab-like scatterer, i.e., we set

[TABLE]

We fix the periodicity to $\varepsilon=2^{-3}$, see Figure 6.4 for a representation of $A$. Since the contrast is rather moderate, the results using $I_{H}$ and $I_{H}^{1}$ are very similar and we only depict computations with $I_{H}^{1}$. We choose the same function $f$ as data term, this time at $x_{0}=(0.25,0.5)$, i.e., closer to the periodic structure. Moreover, we change the wave number and study $k=28$. The reference solution is computed on a fine mesh with $h=2^{-9}$. We emphasize that this a regime where the wave comes into (geometric) resonance with the periodic structure and where homogenization arguments can no longer be applied, see also the discussion in [16].

In this section, we consider the periodic structure and two perturbations (a point and a line defect, see below), which yield physically interesting effects. As we mainly focus on the phenomena here, we do not study convergence histories, but show how (and for which meshes) the upscaled approximation $u_{\operatorname{LOD},m}$ and its macroscopic part $u_{H,m}$ are able to faithfully represent the (qualitative) behavior of the exact solution. Compared to the wave number $k=28$, a rather coarse mesh of $H=2^{-6}$ and a moderate patch size of $m=2$ is sufficient to produce a qualitatively good approximation with our multiscale method in all experiments. Already in the previous subsection, we observed that the periodic setting with $\varepsilon=2^{-3}$ is not such a great challenge from the high contrast point of view so that a rough coupling of $kH\lesssim 1$ seems to be sufficient. We also observe that we need $H=2^{-7}$ and $m=2$ to have a faithful approximation by the macroscopic part $u_{H,m}$. All in all, the experiments of this section show that the multiscale method is able to capture physically relevant settings and is applicable also in the regime where no scale separation between wave length and fine-scale geometry exists.

For the perfectly periodic setting, we observe a sort of lensing effect: The wave pattern behind the scatterer looks (qualitatively) as if a source were located also behind the scatterer, see Figure 6.4. This effect is closely related to negative refraction of the material, see [16, 37]. Next, we introduce a point defect into the periodic structure: We eliminate one inclusion and leave the setting otherwise unchanged, see Figure 6.5. This local defect, however, although rather far away from the point source has a tremendous effect on the behavior of the wave. The wave pattern behind the scatterer is different from the periodic case (for instance, it is no longer symmetric around the axis $y=0.5$), see Figure 6.5. Hence, one may conclude that the lensing effect is destroyed. From the point defect we now move to a line defect: We widen the area where $A=1$ in the middle of the scatterer (around the line $y=0.5$) from $0.5\varepsilon$ to $\varepsilon$, see Figure 6.6. The new (thin) channel now acts as a wave guide (cf. [4]), i.e., the wave issued form the source is mainly traveling through the new channel, see Figure 6.6.

Conclusion

We analyzed the Localized Orthogonal Decomposition for high-contrast high-frequency Helmholtz problems beyond periodicity and scale separation. It yields faithful approximations of the exact solution if the (coarse) mesh size is of the order of the effective wave length. The method relies on the solution of local finescale problems where the size of the local patches only grows logarithmically with the effective wave number. By effective wave number or wave length we denote the actual wave number or wave length divided by square root of the lowest value of the diffusion coefficient. We proved optimal convergence rates in the energy, the $L^{2}$, and a weighted $L^{2}$ norm and linked the results to expectations from the periodic case. The numerical experiments confirmed the theoretical convergence rates and showed the applicability of the method to the Bloch wave regime and periodic structures perturbed by local (point or line) defects. We thereby simulated interesting physical effects such as lensing and wave guides. Since the method is also applicable for totally unstructured or disordered coefficients, unusual and interesting phenomena in random media can be studied in future work.

Acknowledgments

We thank S. Sauter for his comment on the unique continuation principle for the Helmholtz equation and E. A. Spence for the discussion on stability estimates. We are grateful for the remarks of the anonymous reviewers which helped to improve the paper.