Elastic transmission eigenvalues and their computation via the method of fundamental solutions

ABSTRACT

A stabilized version of the fundamental solution method to catch ill-conditioning effects is investigated with focus on the computation of complex-valued elastic interior transmission eigenvalues in two dimensions for homogeneous and isotropic media. Its algorithm can be implemented very shortly and adopts to many similar partial differential equation-based eigenproblems as long as the underlying fundamental solution function can be easily generated. We develop a corroborative approximation analysis which also implicates new basic results for transmission eigenfunctions and present some numerical examples which together prove successful feasibility of our eigenvalue recovery approach.

Keywords:

• Interior transmission eigenvalues
• method of fundamental solutions
• elastic scattering

2010 Mathematics Subject Classifications:

• 65N80
• 65N35
• 35P25

1. Introduction

Interior transmission eigenvalues (ITEs) arise primarily in the study of inverse scattering theory, cf. [1–3]. Their major role comes along with conventional reconstruction algorithms for which incident wave frequencies of that specific order need to be excluded to fully justify the feasibility of qualitative methods for the target recovery process. While the acoustic and electromagnetic cases have been thoroughly investigated in the last years, see [4–13], there are comparably only a few papers covering the numerical study of elastic ITEs, e.g. [14–17], especially when focusing on the overall complex-valued spectrum. It is the purpose of this paper to introduce the approved algorithm from [9,10] also for their computation. As for the latter, we confine to ideally isotropic, homogeneous and planar-symmetric scatterers to resemble the 2D case but allow for less regular shapes within our theoretical foundation now which are assumed to be known in advance.

Numerical challenges in computing ITEs origin from the fact that the underlying interior transmission eigenvalue problem (ITP) is non-self-adjoint, non-elliptic and non-linear in the eigenvalue parameter. Its straightforward discretization would therefore result in non-Hermitian matrices whose pseudospectra are generally harder to capture, especially when the matrix size becomes large. In this context, the advantage of Ritz-type methods such as the method of fundamental solution (MFS) whose prototype goes back to Kupradze in the 1960s for approximating solutions of well-posed boundary value problems, see [18], is that surprisingly accurate results can already be obtained for many regular domains while only a relatively small number of trial functions is used. Due to this observation our MFS-based eigenvalue algorithm similar to [19] will be efficient at remarkably low numerical costs and is still easy-to-implement as being mesh- and integration-free unlike the usual competitive methods, cf. [14].

In the course of substantiating our approach mathematically, we prove some novel findings concerning norm relations of ITP eigenfunction pairs which surprisingly differ for purely real and complex-valued ITEs with non-vanishing imaginary part. Since the field of eigenfunction properties is quite undiscovered unlike their eigenvalues themselves in the ITP context, see [20], we want to point out that most of our results presented here should be adaptable for the acoustic and electromagnetic case, too.

The remainder of this paper is structured as follows: In Section 2, we will setup the modeling framework of the ITP. Section 3 will then consider the eigenvalue approximation problem from an abstract and thus more general perspective that also fits other boundary control techniques such as the boundary element method. Definitions and prerequisite results that are specific for our MFS ansatz will be given in Section 4 and are followed by numerical examples in Section 5. A final conclusion will summarize the merits of our proposed method at the end.

2. Setup of the elastic interior transmission eigenvalue problem

We model time harmonic vibrations of a bounded and elastically penetrable solid D with homogeneous mass density $\rho =\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}>0$ which is assumed to differ from its normalized background material through $\rho \ne 1$. Mathematically, we think of D as a domain $D\subset {\mathbb{R}}^2$ of class $C^{1,1}$ whose regional displacement, when absent or present in comparison with its surrounding, will be modeled by vector fields $u,v:D\to {\mathbb{C}}^2$, respectively. For both the generalized but isotropic version of Hooke's law shall relate their gradients in a linear way to the internal stress $\sigma :{\mathbb{C}}^{2\times 2}\to {\mathbb{C}}^{2\times 2}$ via $\sigma \left(z\right)=2\mu ϵ\left(z\right)+\lambda \mathrm{t}\mathrm{r}\left(ϵ\right(z\left)\right)\mathrm{I}.$ Here, $\mathrm{I}\in {\mathbb{R}}^{2\times 2}$ is the identity matrix and $ϵ\left(z\right)=(\mathrm{\nabla }z+(\mathrm{\nabla }z{)}^{\mathrm{T}})/2$ is the symmetric part of the gradient for any displacement field z measuring its corresponding strain. Further, $\mu ,\lambda$ are Lamé parameters which coincide for simplicity for both the scatterer and the background and which are constrained to $\mu >0$ as well as $2\mu +\lambda >0$ to guarantee strong ellipticity of the governing Navier system, see [21] (1) $\mathrm{d}\mathrm{i}\mathrm{v}\left(\sigma \right(z\left)\right)+\varrho {\omega }^{2}z=(\mu \mathrm{\Delta }z+(\lambda +\mu )\mathrm{\nabla }\left(\mathrm{d}\mathrm{i}\mathrm{v}z\right))+\varrho {\omega }^{2}z=:{\mathrm{\Delta }}^{\ast }z+\varrho {\omega }^{2}z=0.$(1) These equations describe from a physical perspective the spatial part of the elastic wave propagation, where ϱ represents any mass density under consideration. With this notation at hand we may now formulate the elastic interior transmission problem (ITP): (2) $\begin{array}{rl}& {\mathrm{\Delta }}^{\ast }u+{\omega }^{2}u=0\mathrm{i}\mathrm{n}D,\\ & {\mathrm{\Delta }}^{\ast }v+\rho {\omega }^{2}v=0\mathrm{i}\mathrm{n}D,\\ & u=v\mathrm{o}\mathrm{n}\mathrm{\partial }D,\\ & \sigma \left(u\right)\nu =\sigma \left(v\right)\nu \mathrm{o}\mathrm{n}\mathrm{\partial }D,\end{array}$(2) where $u,v\in L^2(D,{\mathbb{C}}^2)$ and the scattered part $(u-v)$ is supposed to be an element of the more regular Sobolev space $H_0^2\left(D\right)=\{\phi \in H^2(D,{\mathbb{C}}^2):\phi =0\mathrm{a}\mathrm{n}\mathrm{d}\sigma (\phi )\nu =0\mathrm{o}\mathrm{n}\mathrm{\partial }D\}$ with outer normal ν along the boundary $\mathrm{\partial }D$. Frequencies $\omega \in \mathbb{C}\mathrm{\setminus }\left\{0\right\}$ which admit non-trivial solutions $(u,v)$ to the above ITP will be called interior transmission eigenvalues (ITEs). Apparently, they refer to those pairs of harmonic waves in D whose behavior along the scattering boundary coincides and which thus lock the possibility of detecting the scatterer on the basis of close-by data.

Note that the co-normal derivative appearing in the Neumann boundary condition of (2(2) $\begin{array}{rl}& {\mathrm{\Delta }}^{\ast }u+{\omega }^{2}u=0\mathrm{i}\mathrm{n}D,\\ & {\mathrm{\Delta }}^{\ast }v+\rho {\omega }^{2}v=0\mathrm{i}\mathrm{n}D,\\ & u=v\mathrm{o}\mathrm{n}\mathrm{\partial }D,\\ & \sigma \left(u\right)\nu =\sigma \left(v\right)\nu \mathrm{o}\mathrm{n}\mathrm{\partial }D,\end{array}$(2) ) is connected naturally to the highest order term of their partial differential equations (PDEs) via integration by parts, also known as Betti's first formula, see [22] (3) $\begin{array}{rl}{\int }_D{\mathrm{\Delta }}^{\ast }\phi \cdot \psi \mathrm{d}x& =-{\int }_D\sigma \left(\phi \right):\mathrm{\nabla }\psi \mathrm{d}x+{\int }_{\mathrm{\partial }D}\left(\sigma \right(\phi \left)\nu \right)\cdot \psi \mathrm{d}s\\ & =-{\int }_{D}2\mu ϵ\left(\phi \right):ϵ\left(\psi \right)+\lambda \mathrm{d}\mathrm{i}\mathrm{v}\phi \mathrm{d}\mathrm{i}\mathrm{v}\psi \mathrm{d}x+{\int }_{\mathrm{\partial }D}\left(\sigma \right(\phi \left)\nu \right)\cdot \psi \mathrm{d}s\end{array}$(3) which holds by duality for any $\psi \in L^2(D,{\mathbb{C}}^2)$ and $\phi \in H^2(D,{\mathbb{C}}^2)$. Here, the colon symbol denotes the Frobenius inner product given by $A:B=\mathrm{t}\mathrm{r}\left(A{B}^{\mathrm{H}}\right)$, whereas the single dot refers to the scalar-product-like operation $a\cdot b=a^{\mathrm{H}}b$ for $a,b\in {\mathbb{C}}^2$. Having thus set the mathematical framework of our eigenproblem, in the sequel we try to recover ITEs ω by solving (2(2) $\begin{array}{rl}& {\mathrm{\Delta }}^{\ast }u+{\omega }^{2}u=0\mathrm{i}\mathrm{n}D,\\ & {\mathrm{\Delta }}^{\ast }v+\rho {\omega }^{2}v=0\mathrm{i}\mathrm{n}D,\\ & u=v\mathrm{o}\mathrm{n}\mathrm{\partial }D,\\ & \sigma \left(u\right)\nu =\sigma \left(v\right)\nu \mathrm{o}\mathrm{n}\mathrm{\partial }D,\end{array}$(2) ) approximately in the sense that we allow for small deviations within the boundary data that are assumed to vanish in some limiting procedure.

3. Approximation analysis of elastic interior transmission eigenvalues via boundary control

We define our relaxed space of trial functions for approximating solutions of (2(2) $\begin{array}{rl}& {\mathrm{\Delta }}^{\ast }u+{\omega }^{2}u=0\mathrm{i}\mathrm{n}D,\\ & {\mathrm{\Delta }}^{\ast }v+\rho {\omega }^{2}v=0\mathrm{i}\mathrm{n}D,\\ & u=v\mathrm{o}\mathrm{n}\mathrm{\partial }D,\\ & \sigma \left(u\right)\nu =\sigma \left(v\right)\nu \mathrm{o}\mathrm{n}\mathrm{\partial }D,\end{array}$(2) ) subject to boundary optimization by (4) $\mathcal{H}:=\bigcup _{0\le \arg \left(\omega \right)<\frac{\pi }{4}}\mathcal{H}\left(\omega \right),$(4) where for any $\omega \in \mathbb{C}\mathrm{\setminus }\left\{0\right\}$ $\mathcal{H}\left(\omega \right):=\{(u,v)\in C^{\mathrm{\infty }}\left(\overline{D}\right)\times C^{\mathrm{\infty }}\left(\overline{D}\right):{\mathrm{\Delta }}^{\ast }u+{\omega }^{2}u=0,{\mathrm{\Delta }}^{\ast }v+{\omega }^2\rho v=0\}.$ By definition, any pair in $\mathcal{H}\left(\omega \right)$ fulfills the required PDE conditions from (2(2) $\begin{array}{rl}& {\mathrm{\Delta }}^{\ast }u+{\omega }^{2}u=0\mathrm{i}\mathrm{n}D,\\ & {\mathrm{\Delta }}^{\ast }v+\rho {\omega }^{2}v=0\mathrm{i}\mathrm{n}D,\\ & u=v\mathrm{o}\mathrm{n}\mathrm{\partial }D,\\ & \sigma \left(u\right)\nu =\sigma \left(v\right)\nu \mathrm{o}\mathrm{n}\mathrm{\partial }D,\end{array}$(2) ) automatically and it is the choice of ω that determines in how far their boundary data are compatible in the sense of the ITP. Note that the above restrictions on ω in $\mathcal{H}$ are due to the following theorem on the overall locations of complex ITEs and the fact that both $\overline{\omega }$, $-\overline{\omega }$ and $-\omega$ are each ITEs if and only if some ω from the first quadrant in the complex plane is. We drop the proof of the former since it would exactly follow the lines from its acoustic analogon in [23] with the obvious operator adaptions.

Theorem 3.1

Let $\omega ={\omega }_1+\mathrm{i}{\omega }_2\in \mathbb{C}\mathrm{\setminus }\left\{0\right\}$ be an ITE with ${\omega }_1,{\omega }_2\in \mathbb{R}$ for the scatterer D. Then it holds that ${\omega }_1^2>{\omega }_2^2\mathrm{a}\mathrm{n}\mathrm{d}{\omega }_1^4+{\omega }_2^4+\frac{2\rho +6}{\delta -1}{\omega }_1{\omega }_2-\frac{\lambda \left(D\right)}{\rho }({\omega }_1^2-{\omega }_2^2)>0,$ where $\lambda \left(D\right)$ is the smallest interior Dirichlet eigenvalue for the Navier problem (1(1) $\mathrm{d}\mathrm{i}\mathrm{v}\left(\sigma \right(z\left)\right)+\varrho {\omega }^{2}z=(\mu \mathrm{\Delta }z+(\lambda +\mu )\mathrm{\nabla }\left(\mathrm{d}\mathrm{i}\mathrm{v}z\right))+\varrho {\omega }^{2}z=:{\mathrm{\Delta }}^{\ast }z+\varrho {\omega }^{2}z=0.$(1) ) with $\varrho =\rho$. In particular, if ω lies in the first quadrant of the complex plane, then $0\le \arg \left(\omega \right)<\pi /4$.

We aim to extract those ω which allow for approximate ITP eigenfunctions in $\mathcal{H}$ with a relatively small ratio of boundary misfit to interior norm. The next theorem states that if these residual quotients can be made arbitrarily small while the corresponding ω accumulate, their limit is indeed an ITE. Its proof slightly refines the technique from [9, Theorem 2] to encompass also complex-valued eigenvalues now.

Theorem 3.2

Assume that $\left\{\right(u_m,v_m,{\omega }_m){\}}_{m\in \mathbb{N}}\subset \mathcal{H}\times \mathbb{C}$ fulfill for some $1\le C<\mathrm{\infty }$ the following conditions:

1. eigenvalue convergence: ${\omega }_m\to \omega$ such that $\mathrm{a}\mathrm{r}\mathrm{g}\left(\omega \right)<1/4,$

2. uniform interior bound: $1/C\le (\parallel u_m{\parallel }_{L^2(D,{\mathbb{C}}^2)}^2+\parallel v_m{\parallel }_{L^2(D,{\mathbb{C}}^2)}^2)\le C$ for all m large enough,

3. vanishing boundary misfit: $(\parallel u_m-v_m{\parallel }_{H^{\frac{3}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}+\parallel \sigma (u_m-v_m\left)\nu {\parallel }_{H^{\frac{1}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}\right)\to 0$ when $m\to \mathrm{\infty }$.

Then, ω is an ITE and a subsequence of $(u_m,v_m)$ converges weakly to some eigenfunction pair.

Proof.

By rescaling and redefining $(u_m,v_m)$, if necessary, we can assume without loss of generality that C = 1 and aim to apply weak compactness in order to construct a solution candidate which will indeed meet all the required ITE criteria. By assumption $\left(ii\right)$ we know (modulo the extraction of subsequences which we will not relabel in m) that $u_m\rightharpoonup u$ and $v_m\rightharpoonup v$ in $L^2(D,{\mathbb{C}}^2)$ which implies that $(u,v)\in L^2(D,{\mathbb{C}}^2)\times L^2(D,{\mathbb{C}}^2)$ fulfills the interior conditions of the ITP, or equivalently the PDE system (1(1) $\mathrm{d}\mathrm{i}\mathrm{v}\left(\sigma \right(z\left)\right)+\varrho {\omega }^{2}z=(\mu \mathrm{\Delta }z+(\lambda +\mu )\mathrm{\nabla }\left(\mathrm{d}\mathrm{i}\mathrm{v}z\right))+\varrho {\omega }^{2}z=:{\mathrm{\Delta }}^{\ast }z+\varrho {\omega }^{2}z=0.$(1) ), according to $\begin{array}{rl}{\int }_{D}u\cdot ({\mathrm{\Delta }}^{\ast }\phi +{\omega }^2\phi )\mathrm{d}x& =\underset{m\to \mathrm{\infty }}{lim}{\int }_D{u}_m\cdot ({\mathrm{\Delta }}^{\ast }\phi +{\omega }_m^2\phi )\mathrm{d}x\\ & =\underset{m\to \mathrm{\infty }}{lim}{\int }_D({\mathrm{\Delta }}^{\ast }{u}_m+{\omega }^2{u}_m)\cdot \phi \mathrm{d}x=0\end{array}$ for any bump function $\phi \in C_c^{\mathrm{\infty }}(D,{\mathbb{C}}^2)$ and with a similar calculus for v. In order to prove that $(u-v)$ also has the correct ITP boundary data, it suffices to prove that the differences $(u_m-v_m)$ are bounded in $H^2(D,{\mathbb{C}}^2)$ in combination with assumption $\left(iii\right)$ because $w\mapsto w_{|\mathrm{\partial }D}$ and $w\mapsto \sigma (w{)}_{|\mathrm{\partial }D}\nu$ are continuous as operators $H^2(D,{\mathbb{C}}^2)\to H^{\frac{3}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)$ and $H^2(D,{\mathbb{C}}^2)\to H^{\frac{1}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)$, respectively, due to our regularity assumptions on $\mathrm{\partial }D$.

Noting that $(u_m-v_m)$ are solutions to the inhomogeneous Navier system ${\mathrm{\Delta }}^{\ast }(u_m-v_m)+{\omega }_m^2(u_m-v_m)=(1-\rho )v_m\mathrm{i}\mathrm{n}D,$ elliptic estimates like in [21] tell us that (5) $\parallel u_m-v_m{\parallel }_{H^2(D,{\mathbb{C}}^2)}\le C\left(\parallel u_m-v_m{\parallel }_{H^{\frac{3}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}+\parallel u_m{\parallel }_{L^2(D,{\mathbb{C}}^2)}+\parallel v_m{\parallel }_{L^2(D,{\mathbb{C}}^2)}\right)$(5) which therefore gives the desired uniform bound with respect to m.

It remains to show that $(u,v)$ is non-trivial. For this we recall that the embedding $H^2(D,{\mathbb{C}}^2)\hookrightarrow L^2(D,{\mathbb{C}}^2)$ is compact which implies $(u_m-v_m)\to (u-v)$ strongly in $L^2(D,{\mathbb{C}}^2)$. Apparently, $(u,v)\ne 0$ if $\parallel u-v{\parallel }_{L^2(D,{\mathbb{C}}^2)}>0$ so we will assume contrarily that $(u_m-v_m)\to 0$ in $L^2(D,{\mathbb{C}}^2)$. Then, on the one hand, the bounded sequence $a_m:={\int }_D{u}_m\cdot {\overline{v}}_m\mathrm{d}x$ may be singled out to converge to some $a\in \mathbb{C}$ from which we know by $\left(ii\right)$ with C = 1 that (6) $\frac{1}{2}\ge |a|\ge \mathrm{R}\mathrm{e}a=\underset{m\to \mathrm{\infty }}{lim}\frac{\parallel u_m{\parallel }_{L^2(D,{\mathbb{C}}^2)}^2+\parallel v_m{\parallel }_{L^2(D,{\mathbb{C}}^2)}^2-\parallel u_m-v_m{\parallel }_{L^2(D,{\mathbb{C}}^2)}^2}{2}\ge \frac{1}{2},$(6) i.e. $a=1/2$. Since $\mathrm{a}\mathrm{r}\mathrm{g}\left(\omega \right)<1/4$ and $\rho \ne 1$, we can even conclude $|\mathrm{R}\mathrm{e}a({\omega }^2-\rho {\overline{\omega }}^2)|>0$.

On the other hand, (3(3) $\begin{array}{rl}{\int }_D{\mathrm{\Delta }}^{\ast }\phi \cdot \psi \mathrm{d}x& =-{\int }_D\sigma \left(\phi \right):\mathrm{\nabla }\psi \mathrm{d}x+{\int }_{\mathrm{\partial }D}\left(\sigma \right(\phi \left)\nu \right)\cdot \psi \mathrm{d}s\\ & =-{\int }_{D}2\mu ϵ\left(\phi \right):ϵ\left(\psi \right)+\lambda \mathrm{d}\mathrm{i}\mathrm{v}\phi \mathrm{d}\mathrm{i}\mathrm{v}\psi \mathrm{d}x+{\int }_{\mathrm{\partial }D}\left(\sigma \right(\phi \left)\nu \right)\cdot \psi \mathrm{d}s\end{array}$(3) ) with exchanged roles of its test functions $\phi ,\psi \in C^{\mathrm{\infty }}\left(\overline{D}\right)$ yields the analogon of Green's second identity ${\int }_D\psi \cdot {\mathrm{\Delta }}^{\ast }\phi -\phi \cdot {\mathrm{\Delta }}^{\ast }\psi \mathrm{d}x={\int }_{\mathrm{\partial }D}\psi \cdot \left(\sigma \right(\phi \left)\nu \right)-\phi \cdot \left(\sigma \right(\psi \left)\nu \right)\mathrm{d}s.$ Duality shows that the separated boundary contributions of the approximate eigenfunctions $\parallel u_m{\parallel }_{H^{-\frac{1}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}$,$\parallel v_m{\parallel }_{H^{-\frac{1}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}$,$\parallel \sigma \left(u_m\right)\nu {\parallel }_{H^{-\frac{3}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}$ and $\parallel \sigma \left(v_m\right)\nu {\parallel }_{H^{-\frac{3}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}$ are even uniformly bounded in m thanks to $\left(i\right)$–$\left(iii\right)$. Therefore we may further compute $\begin{array}{rl}0& <|\mathrm{R}\mathrm{e}a({\omega }^2-\rho {\overline{\omega }}^2)|\\ & =\underset{m\to \mathrm{\infty }}{lim}\left|\mathrm{R}\mathrm{e}{\int }_D({\omega }^2-\rho {\overline{\omega }}^2)u_m\cdot {\overline{v}}_m\mathrm{d}x\right|\\ & =\underset{m\to \mathrm{\infty }}{lim}\left|\mathrm{R}\mathrm{e}{\int }_D\left({\omega }_m^2{u}_m\right)\cdot {\overline{v}}_m-u_m\cdot \left(\rho {\overline{\omega }}_m^2{\overline{v}}_m\right)\mathrm{d}x\right|\\ & =\underset{m\to \mathrm{\infty }}{lim}\left|\mathrm{R}\mathrm{e}{\int }_D{\overline{v}}_m\cdot {\mathrm{\Delta }}^{\ast }{u}_m-u_m\cdot {\mathrm{\Delta }}^{\ast }{\overline{v}}_m\mathrm{d}x\right|\\ & =\underset{m\to \mathrm{\infty }}{lim}\left|\mathrm{R}\mathrm{e}{\int }_{\mathrm{\partial }D}{\overline{v}}_m\cdot \left(\sigma \right(u_m\left)\nu \right)-u_m\cdot \left(\sigma \right({\overline{v}}_m\left)\nu \right)\mathrm{d}s\right|\\ & =\underset{m\to \mathrm{\infty }}{lim}\left|\mathrm{R}\mathrm{e}{\int }_{\mathrm{\partial }D}{\overline{v}}_m\cdot \left(\sigma \right(u_m\left)\nu \right)-{\overline{u}}_m\cdot \left(\sigma \right(v_m\left)\nu \right)\mathrm{d}s\right|\\ & =\underset{m\to \mathrm{\infty }}{lim}\left|\mathrm{R}\mathrm{e}{\int }_{\mathrm{\partial }D}({\overline{v}}_m-{\overline{u}}_m)\cdot \left(\sigma \right(u_m\left)\nu \right)+{\overline{u}}_m\cdot \left(\sigma \right(u_m-v_m\left)\nu \right)\mathrm{d}s\right|\\ & \le \underset{m\to \mathrm{\infty }}{lim}\parallel v_m-u_m{\parallel }_{H^{\frac{3}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}\parallel \sigma \left(u_m\right)\nu {\parallel }_{H^{-\frac{3}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}\\ & +\underset{m\to \mathrm{\infty }}{lim}\parallel \sigma (v_m-u_m)\nu {\parallel }_{H^{\frac{1}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}\parallel u_m{\parallel }_{H^{-\frac{1}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}\\ & =0,\end{array}$ which is a contradiction.

In order to establish for fixed m, some more qualitative relation between the boundary misfit of some sufficiently approximate eigenfunction pair and its absolute eigenvalue deviation from the actual one, we derive an interconnecting estimate in the following which can also be seen as practical guiding principle for numerical calculations. Similar to the acoustic case, see [9, Lemma 5], but now improved for even less regular shapes, its rigorous validity relies on a critical integral expression which must not vanish.

Lemma 3.3

Let $(u,v)$ be an ITP eigenfunction pair with ITE ω and assume that $(\tilde{u},\tilde{v})\in \mathcal{H}\left(\tilde{\omega }\right)$ for some arbitrary frequency $\tilde{\omega }$. If the integral constraint (7) $|{\int }_{D}u\cdot \tilde{u}-\rho v\cdot \tilde{v}\mathrm{d}x|\ge \tilde{\epsilon }>0$(7) is fulfilled, then there exists a constant C>0 determined only by the boundary data of u (or equivalently $v)$ such that (8) $|{\omega }^2-{\tilde{\omega }}^2|\le \frac{C}{\tilde{\epsilon }}\sqrt{\parallel \tilde{u}-\tilde{v}{\parallel }_{H^{\frac{3}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}^2+\parallel \sigma (\tilde{u}-\tilde{v})\nu {\parallel }_{H^{\frac{1}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}^2}.$(8)

Proof.

Applying Betti's formula twice, using the identical ITP boundary data for u and v, yields $\begin{array}{rl}& {\omega }^2{\int }_D\tilde{u}\cdot u-\rho \tilde{v}\cdot v\mathrm{d}x\\ & ={\int }_D-\tilde{u}\cdot {\mathrm{\Delta }}^{\ast }u+\tilde{v}\cdot {\mathrm{\Delta }}^{\ast }v\mathrm{d}x\\ & ={\int }_{D}2\mu ϵ\left(\tilde{u}\right):ϵ\left(u\right)+\lambda \mathrm{d}\mathrm{i}\mathrm{v}\tilde{u}\mathrm{d}\mathrm{i}\mathrm{v}u\mathrm{d}x-{\int }_{D}2\mu ϵ\left(\tilde{v}\right):ϵ\left(v\right)+\lambda \mathrm{d}\mathrm{i}\mathrm{v}\tilde{v}\mathrm{d}\mathrm{i}\mathrm{v}v\mathrm{d}x\\ & -{\int }_{\mathrm{\partial }D}(\tilde{u}-\tilde{v})\cdot \left(\sigma \right(u\left)\nu \right)\mathrm{d}s\\ & ={\int }_D({\mathrm{\Delta }}^{\ast }\tilde{u}\cdot u+{\mathrm{\Delta }}^{\ast }\tilde{v}\cdot v)\mathrm{d}x+{\int }_{\mathrm{\partial }D}\left(\sigma \right(\tilde{u}-\tilde{v}\left)\nu \right)\cdot u\mathrm{d}s-{\int }_{\mathrm{\partial }D}(\tilde{u}-\tilde{v})\cdot \left(\sigma \right(u\left)\nu \right)\mathrm{d}s\\ & ={\tilde{\omega }}^2{\int }_D\tilde{u}\cdot u-\rho \tilde{v}\cdot v\mathrm{d}x+{\int }_{\mathrm{\partial }D}\left(\sigma \right(\tilde{u}-\tilde{v}\left)\nu \right)\cdot u\mathrm{d}s-{\int }_{\mathrm{\partial }D}(\tilde{u}-\tilde{v})\cdot \left(\sigma \right(u\left)\nu \right)\mathrm{d}s\end{array}$ and after rearranging we obtain $({\omega }^2-{\tilde{\omega }}^2){\int }_D\tilde{u}\cdot u-\rho \tilde{v}\cdot v\mathrm{d}x={\int }_{\mathrm{\partial }D}\left(\sigma \right(\tilde{u}-\tilde{v}\left)\nu \right)\cdot u\mathrm{d}s-{\int }_{\mathrm{\partial }D}(\tilde{u}-\tilde{v})\cdot \left(\sigma \right(u\left)\nu \right)\mathrm{d}s.$ Taking absolute values gives $\begin{array}{rl}|{\omega }^2-{\tilde{\omega }}^2|& \le \frac{1}{\tilde{\epsilon }}\left({\int }_{\mathrm{\partial }D}|\left(\sigma \right(\tilde{u}-\tilde{v}\left)\nu \right)\cdot u|\mathrm{d}s+{\int }_{\mathrm{\partial }D}|(\tilde{u}-\tilde{v})\cdot \left(\sigma \right(u\left)\nu \right)|\mathrm{d}s\right)\\ & \le \frac{1}{\tilde{\epsilon }}\left(\parallel \sigma (\tilde{u}-\tilde{v})\nu {\parallel }_{H^{\frac{1}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}\parallel u{\parallel }_{H^{-\frac{1}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}+\parallel \tilde{u}-\tilde{v}{\parallel }_{H^{\frac{3}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}\parallel \sigma \left(u\right)\nu {\parallel }_{H^{-\frac{3}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}\right)\\ & \le \frac{C}{\tilde{\epsilon }}\sqrt{\parallel \tilde{u}-\tilde{v}{\parallel }_{H^{\frac{3}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}^2+\parallel \sigma (\tilde{u}-\tilde{v})\nu {\parallel }_{H^{\frac{1}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}{\parallel }^2},\end{array}$ where $C:=\sqrt{\parallel u{\parallel }_{H^{-\frac{1}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}^2+\parallel \sigma \left(u\right)\nu {\parallel }_{H^{-\frac{3}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}^2}.$ Duality shows that $C<\mathrm{\infty }$ since u solves the Navier equation.

We state a direct consequence with respect to the $L^2(D,{\mathbb{C}}^2)$-norm of eigenfunctions for frequencies $\omega \in \mathbb{C}\mathrm{\setminus }\mathbb{R}$.

Corollary 3.4

If $(u,v)$ is an ITP eigenfunction pair with $\omega \in \mathbb{C}\mathrm{\setminus }\mathbb{R},$ then we have that $\parallel u{\parallel }_{L^2(D,{\mathbb{C}}^2)}^2=\rho \parallel v{\parallel }_{L^2(D,{\mathbb{C}}^2)}^2$.

Proof.

Let $\left\{\right(u_m,v_m){\}}_{m\in \mathbb{N}}$ be a sequence in $\mathcal{H}\left(\omega \right)$ such that $w_m:=(u_m-v_m)\to (u-v)$ in $H^2(D,{\mathbb{C}}^2)$, see for example Theorem 4.2 later. In particular, $({\overline{u}}_m,{\overline{v}}_m)\in \mathcal{H}\left(\overline{\omega }\right)$, $\parallel {\overline{w}}_m{\parallel }_{H^{\frac{3}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}^2+\parallel \sigma \left({\overline{w}}_m\right){\parallel }_{H^{\frac{1}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}^2\to 0$ and the left hand side of (7(7) $|{\int }_{D}u\cdot \tilde{u}-\rho v\cdot \tilde{v}\mathrm{d}x|\ge \tilde{\epsilon }>0$(7) ) converges because of ${\overline{u}}_m=({\mathrm{\Delta }}^{\ast }{\overline{w}}_m+\rho {\overline{\omega }}^2{\overline{w}}_m)/\left(\overline{\omega }\right(1-\rho \left)\right)$ and ${\overline{v}}_m=({\mathrm{\Delta }}^{\ast }{\overline{w}}_m+{\overline{\omega }}^2{\overline{w}}_m)/\left(\overline{\omega }\right(1-\rho \left)\right)$ to $\parallel u{\parallel }_{L^2(D,{\mathbb{C}}^2)}^2-\rho \parallel v{\parallel }_{L^2(D,{\mathbb{C}}^2)}^2$ when setting $\tilde{u}:={\overline{u}}_m$ and $\tilde{v}:={\overline{v}}_m$. Since $|{\omega }^2-{\overline{\omega }}^2|>0$ by assumption, the right hand side of (8(8) $|{\omega }^2-{\tilde{\omega }}^2|\le \frac{C}{\tilde{\epsilon }}\sqrt{\parallel \tilde{u}-\tilde{v}{\parallel }_{H^{\frac{3}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}^2+\parallel \sigma (\tilde{u}-\tilde{v})\nu {\parallel }_{H^{\frac{1}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}^2}.$(8) ) forces any uniform bound on $\tilde{\epsilon }$ to vanish in the limit $m\to \mathrm{\infty }$.

The situation is different for $\omega \in \mathbb{R}$ and especially when it is the ITE of smallest magnitude. In this case, with the additional assumption that $\rho >1$ or $0<\rho <1$ is sufficiently large or small to be made precise next, respectively, we can show that $\parallel u{\parallel }_{L^2(D,{\mathbb{C}}^2)}^2-\rho \parallel v{\parallel }_{L^2(D,{\mathbb{C}}^2)}^2\ne 0$ which also guarantees a uniform bound $\tilde{\epsilon }>0$ for approximate eigenfunction pairs $(\tilde{u},\tilde{v})$ in the vicinity of $(u,v)$ (cf. Corollary 3.6 later). To avoid misleading confusion at that point, we want to emphasize that Corollary 3.4 does not imply the non-existence of some positive constant $\tilde{\epsilon }$ when dealing with complex-valued ITEs ω admitting a non-vanishing imaginary part.

Theorem 3.5

Let ω be the smallest real-valued ITE for the scatterer D with constant density $\rho \ne 1$. If $\rho >1$ is large or $0<\rho <1$ is small enough, where the corresponding thresholds depend only on D and the Lamé parameters $\mu ,\lambda ,$ we have the relations $\parallel u{\parallel }_{L^2(D,{\mathbb{C}}^2)}^2-\rho \parallel v{\parallel }_{L^2(D,{\mathbb{C}}^2)}^2<0\mathrm{o}\mathrm{r}\parallel u{\parallel }_{L^2(D,{\mathbb{C}}^2)}^2-\rho \parallel v{\parallel }_{L^2(D,{\mathbb{C}}^2)}^2>0,$ respectively.

Proof.

For the sake of presentation, we will assume that $\rho >1$ since the case $0<\rho <1$ works structurally similar. Because u and v can be expressed each in terms of their difference w: = u−v by $u=({\mathrm{\Delta }}^{\ast }w+\rho {\omega }^{2}w)/\left({\omega }^2\right(\rho -1\left)\right)$ and $v=({\mathrm{\Delta }}^{\ast }w+{\omega }^{2}w)/\left({\omega }^2\right(\rho -1\left)\right)$, the basic idea of our proof will be to exploit isometry of the Fourier transform with respect to the single field w to obtain a well-behaved algebraic integrand in terms of ρ. The fact $w\in H_0^2\left(D\right)$ then shows that u, v, w extend naturally by zero outside of D so Plancherel's identity gives $\begin{array}{rl}{\int }_D|u{|}^2-\rho |v{|}^2\mathrm{d}x& =\frac{1}{(2\pi {)}^2}{\int }_{{\mathbb{R}}^2}|\mathcal{F}u{|}^2-\rho |\mathcal{F}v{|}^2\mathrm{d}\xi \\ & =\frac{1}{(2\pi {)}^2}{\int }_{{\mathbb{R}}^2}\frac{|\mathcal{F}({\mathrm{\Delta }}^{\ast }w+\rho {\omega }^{2}w){|}^2-\rho |\mathcal{F}({\mathrm{\Delta }}^{\ast }w+{\omega }^{2}w){|}^2}{(\rho -1{)}^2{\omega }^4}\mathrm{d}\xi \\ & =\frac{1}{(2\pi {)}^2}{\int }_{{\mathbb{R}}^2}\frac{-|\mathcal{F}{\mathrm{\Delta }}^{\ast }w{|}^2+\rho {\omega }^4|\mathcal{F}w{|}^2}{(\rho -1){\omega }^4}\mathrm{d}\xi \\ & \le \frac{1}{(2\pi {)}^2}{\int }_{{\mathbb{R}}^2}\frac{-c|\xi {|}^4+\rho {\omega }^4}{(\rho -1){\omega }^4}|\mathcal{F}w{|}^2\mathrm{d}\xi .\end{array}$ In the last step, we employed the pointwise estimate $|\mathcal{F}{\mathrm{\Delta }}^{\ast }w\left(\xi \right){|}^2\ge c|\xi {|}^4|\mathcal{F}w\left(\xi \right){|}^2$ with $c:=min\{\mu ,\lambda +2\mu \}>0$ inherited from strong ellipticity. As will be shown later, we have for the smallest ITE $\omega =\omega \left(\rho \right)$ with density ρ that $\rho {\omega }^4\to 0$ and likewise $R\left(\rho \right)\to 0$ for $\rho \to \mathrm{\infty }$, where $R\left(\rho \right):=\omega \sqrt[4]{(2\rho -1)/c}$ is the solution of $p\left(t\right):=(-c{t}^4+\rho {\omega }^4)/(\rho -1){\omega }^4=-1$ in t. With $B_{R\left(\rho \right)}$ being the disc centered at the origin of radius $R\left(\rho \right)$, we split the integral above and exploit monotonicity of p to deduce $\begin{array}{rl}& {\int }_{B_{R\left(\rho \right)}}\frac{-c|\xi {|}^4+\rho {\omega }^4}{(\rho -1){\omega }^4}|\mathcal{F}w{|}^2\mathrm{d}\xi +{\int }_{{\mathbb{R}}^2\mathrm{\setminus }{B}_{R\left(\rho \right)}}\frac{-c|\xi {|}^4+\rho {\omega }^4}{(\rho -1){\omega }^4}|\mathcal{F}w{|}^2\mathrm{d}\xi \\ & \le p\left(0\right){\int }_{B_{R\left(\rho \right)}}|\mathcal{F}w{|}^2\mathrm{d}\xi +p\left(R\right(\rho \left)\right){\int }_{{\mathbb{R}}^2\mathrm{\setminus }{B}_{R\left(\rho \right)}}|\mathcal{F}w{|}^2\mathrm{d}\xi \\ & =\frac{\rho }{\rho -1}\parallel \mathcal{F}w{\parallel }_{L^2(B_{R\left(\rho \right)},{\mathbb{C}}^2)}^2-\left(\parallel \mathcal{F}w{\parallel }_{L^2({\mathbb{R}}^2,{\mathbb{C}}^2)}^2-\parallel \mathcal{F}w{\parallel }_{L^2(B_{R\left(\rho \right)},{\mathbb{C}}^2)}^2\right)\\ & \le \left(\frac{\rho }{\rho -1}+1\right)\parallel \mathcal{F}w{\parallel }_{L^2(B_{R\left(\rho \right)},{\mathbb{C}}^2)}^2-(2\pi {)}^2\parallel w{\parallel }_{L^2(D,{\mathbb{C}}^2)}^2.\end{array}$ The first summand can be made arbitrarily small in terms of ρ because of $\begin{array}{rl}\parallel \mathcal{F}w{\parallel }_{L^2(B_{R\left(\rho \right)},{\mathbb{C}}^2)}^2& \le (max|\mathcal{F}w|{)}^2\pi R(\rho {)}^2\\ & \le \parallel w{\parallel }_{L^1(D,{\mathbb{C}}^2)}^2\pi R(\rho {)}^2\\ & \le \parallel w{\parallel }_{L^2(D,{\mathbb{C}}^2)}^{2}R\left(\rho {)}^2{\mathcal{L}}^2\right(D),\end{array}$ where ${\mathcal{L}}^2\left(D\right)$ denotes the two-dimensional Lebesgue measure of D. Putting everything together, we finally obtain ${\int }_D|u{|}^2-\rho |v{|}^2\mathrm{d}x=\parallel w{\parallel }_{L^2(D,{\mathbb{C}}^2)}^2\left(\left(\frac{\rho }{\rho -1}+1\right){\mathcal{L}}^2\left(D\right)R(\rho {)}^2-(2\pi {)}^2\right)<0$ for ρ large enough due to the decay of $R\left(\rho \right)$.

As announced for the latter, it remains to show that $\rho {\omega }^4\to 0$ as $\rho \to \mathrm{\infty }$, where ω is the smallest real-valued ITE of D with density ρ. According to Corollary 1 from [1], the magnitude of ω can be bounded from above by the smallest ITE from any included disc $B_r\subset D$ which thus amounts to show the asymptotics assertion just for the unit disc as scatterer. In this case, as was derived in the appendix of [14] with an ansatz for a purely compressional wave instead of its orthogonal shear part that we will consider now for completion, ω is some ITE corresponding to a radial symmetric eigenfunction if and only if (9) $\mathrm{d}\mathrm{e}\mathrm{t}\left(\begin{array}{cc}{J}_1\left(\omega \sqrt{\frac{1}{\mu }}\right)& J_1\left(\omega \sqrt{\frac{\rho }{\mu }}\right)\\ \omega \sqrt{\frac{1}{\mu }}{J}_1^{\prime }\left(\omega \sqrt{\frac{1}{\mu }}\right)& \omega \sqrt{\frac{\rho }{\mu }}{J}_1^{\prime }\left(\omega \sqrt{\frac{\rho }{\mu }}\right)\end{array}\right)=0,$(9) where $J_1$ is the first Bessel function of order one. This condition can be restated as finding roots ω of the piecewise continuous function $g\left(\omega \right):=f\left(\omega \right)-f\left(\sqrt{\rho }\omega \right)$ which can be recursively decomposed into $f\left(x\right)=h\left(\frac{x}{\sqrt{\mu }}\right)\mathrm{a}\mathrm{n}\mathrm{d}h\left(y\right):=\frac{y{J}_1^{\prime }\left(y\right)}{J_1\left(y\right)}.$ Let $j_1<j_2$ be the two smallest positive roots of $J_1$ and choose $\rho >1$ large enough to have $j_2<j_1\sqrt{\rho }$. Then set ${\omega }_1:=j_1\sqrt{\mu /\rho }$ as well as ${\omega }_2:=j_2\sqrt{\mu /\rho }$ and observe that g is singular at those points, but continuous in between. Also, those poles have different signs according to $\underset{\omega \searrow {\omega }_1}{lim}g\left(\omega \right)=-\underset{\omega \searrow {\omega }_1}{lim}f\left(\sqrt{\rho }\omega \right)=-\mathrm{\infty }\mathrm{a}\mathrm{n}\mathrm{d}\underset{\omega \nearrow {\omega }_2}{lim}g\left(\omega \right)=-\underset{\omega \nearrow {\omega }_2}{lim}f\left(\sqrt{\rho }\omega \right)=\mathrm{\infty },$ which follows from the basic facts that $J_1<0$ in $(j_1,j_2)$, $J_1^{\prime }\left(j_1\right)<0$ but $J_1^{\prime }\left(j_2\right)>0$ and that both $f\left({\omega }_1\right)$, $f\left({\omega }_2\right)$ are finite. Therefore, we can make use of the intermediate value theorem which guarantees for any large ρ a root ω of g fulfilling ${\omega }_1\le \omega \le {\omega }_2$ or equivalently the uniform bound $j_1^2\mu \le \rho {\omega }^2\le j_2^2\mu$ as $\rho \to \mathrm{\infty }$. In particular, $\rho {\omega }^4\to 0$ for the same limit procedure which finally proves our lemma.

Finally, we cite a result from [9, Corollary 6] that summarizes the previous findings for our conceptual ITE recovery approach.

Corollary 3.6

Let the conditions of Theorem 3.2 hold for $\left\{\right(u_m,v_m,{\omega }_m){\}}_{m\in \mathbb{N}}\subset \mathcal{H}\times \mathbb{C}$ with ITE ω and eigenfunction pair $(u,v)\in L^2(D,{\mathbb{C}}^2)\times L^2(D,{\mathbb{C}}^2)$. Assume additionally that (10) ${\int }_D(u^2-\rho v^2)\mathrm{d}x\ne 0.$(10) Then, for sufficiently large $m\in \mathbb{N},$ we have (modulo the relabeling of the weakly convergent subsequence) $|{\omega }^2-{\omega }_m^2|\le C\sqrt{\parallel u_m-v_m{\parallel }_{H^{\frac{3}{2}}\left(\mathrm{\partial }D\right)}^2+\parallel \sigma (u_m-v_m)\nu {\parallel }_{H^{\frac{1}{2}}\left(\mathrm{\partial }D\right)}^2},$ where C>0 depends on the boundary data of u (or equivalently $v)$ and the magnitude of (10(10) ${\int }_D(u^2-\rho v^2)\mathrm{d}x\ne 0.$(10) ).

Proof.

$(u_m,v_m)\rightharpoonup (u,v)$ in $L^2(D,{\mathbb{C}}^2)\times L^2(D,{\mathbb{C}}^2)$ implies with (10(10) ${\int }_D(u^2-\rho v^2)\mathrm{d}x\ne 0.$(10) ) that $\tilde{\epsilon }>0$ in Lemma 3.3 uniformly for m large enough and thus the existence of C>0.

Altogether, we have proved that the conditioned process from Theorem 3.2 for detecting ITEs is spurious free as limiting procedure and provided a constrained a posteriori estimate for the eigenvalue approximation accuracy at each step m. To meet the conditions therein, we will construct ‘simple’ trial functions in the next section that are unspecified in (4(4) $\mathcal{H}:=\bigcup _{0\le \arg \left(\omega \right)<\frac{\pi }{4}}\mathcal{H}\left(\omega \right),$(4) ) so far.

4. The method of fundamental solutions

4.1. The abstract setting

In this section, we want to focus on how to generate interior solutions to a given PDE on the basis of its fundamental solution. A fundamental solution for the free space Navier system (1(1) $\mathrm{d}\mathrm{i}\mathrm{v}\left(\sigma \right(z\left)\right)+\varrho {\omega }^{2}z=(\mu \mathrm{\Delta }z+(\lambda +\mu )\mathrm{\nabla }\left(\mathrm{d}\mathrm{i}\mathrm{v}z\right))+\varrho {\omega }^{2}z=:{\mathrm{\Delta }}^{\ast }z+\varrho {\omega }^{2}z=0.$(1) ) with arguments $x\ne y\in {\mathbb{R}}^2$, constant coefficients $\mu ,\lambda ,\varrho$ and unknown eigenfrequency parameter ω is given by (11) $\begin{array}{rl}{\mathrm{\Phi }}_{\varrho {\omega }^2}(x,y)& :={\mathrm{\Phi }}_{\varrho {\omega }^2}^{\mu ,\lambda }(|x-y|)\\ & :=\frac{\mathrm{i}}{4\mu }{H}_0^{\left(1\right)}(k_s|x-y|)\cdot \mathrm{I}+\frac{\mathrm{i}}{4{\omega }^2}{\mathrm{\nabla }}^{\mathrm{\top }}\mathrm{\nabla }\left(H_0^{\left(1\right)}(k_p|x-y|)-H_0^{\left(1\right)}(k_s|x-y|)\right),\end{array}$(11) with wave numbers $k_s^2:=\frac{\varrho {\omega }^2}{\mu },k_p^2:=\frac{\varrho {\omega }^2}{\lambda +2\mu }$ and $H_0^{\left(1\right)}$ being the first Hankel function of order zero. The indices s and p originate from the well-known Helmholtz decomposition that divides any properly decaying solution $e=e_p+e_s$ of the exterior Navier problem into its compressional and shear wave field. More precisely, $e_p=(-1/k_p^2)\mathrm{\nabla }\left(\mathrm{d}\mathrm{i}\mathrm{v}e\right)$ and $e_s=(-1/k_s^2)(-{\mathrm{\partial }}_2({\mathrm{\partial }}_1{e}_2-{\mathrm{\partial }}_2{e}_1),{\mathrm{\partial }}_1({\mathrm{\partial }}_1{e}_2-{\mathrm{\partial }}_2{e}_1){)}^{\mathrm{\top }}$ solve $\mathrm{\Delta }{e}_s+k_s^2{e}_s=0\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{\Delta }{e}_p+k_p^2{e}_p=0\mathrm{i}\mathrm{n}{\mathbb{R}}^2\mathrm{\setminus }D$ and are constrained to fulfill Kupradze's (outgoing) radiation condition in two dimensions $\underset{r\to \mathrm{\infty }}{lim}\sqrt{r}({\mathrm{\partial }}_r{e}_s-\mathrm{i}{k}_s{e}_s)=0,\underset{r\to \mathrm{\infty }}{lim}\sqrt{r}({\mathrm{\partial }}_r{e}_p-\mathrm{i}{k}_p{e}_p)=0$ uniformly in all directions of the radial distance r from the origin. In particular, freezing one of the arguments in (11(11) $\begin{array}{rl}{\mathrm{\Phi }}_{\varrho {\omega }^2}(x,y)& :={\mathrm{\Phi }}_{\varrho {\omega }^2}^{\mu ,\lambda }(|x-y|)\\ & :=\frac{\mathrm{i}}{4\mu }{H}_0^{\left(1\right)}(k_s|x-y|)\cdot \mathrm{I}+\frac{\mathrm{i}}{4{\omega }^2}{\mathrm{\nabla }}^{\mathrm{\top }}\mathrm{\nabla }\left(H_0^{\left(1\right)}(k_p|x-y|)-H_0^{\left(1\right)}(k_s|x-y|)\right),\end{array}$(11) ) such as y without loss of generality, these properties apply column-wise (and by symmetry of the fundamental matrix also row-wise) to $x\mapsto {\mathrm{\Phi }}_{\varrho {\omega }^2}(x,y)$ in ${\mathbb{R}}^2\mathrm{\setminus }\left\{y\right\}$ for any source point $y\in {\mathbb{R}}^2$ and density $\varrho =\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}$. Choosing some simply closed and sufficiently smooth contour $\mathrm{\Gamma }\in {\mathbb{R}}^2\mathrm{\setminus }D$, called the artificial or source boundary, we easily see that any coefficient function $c:\mathrm{\Gamma }\to {\mathbb{C}}^2$ would generate a smooth solution of (1(1) $\mathrm{d}\mathrm{i}\mathrm{v}\left(\sigma \right(z\left)\right)+\varrho {\omega }^{2}z=(\mu \mathrm{\Delta }z+(\lambda +\mu )\mathrm{\nabla }\left(\mathrm{d}\mathrm{i}\mathrm{v}z\right))+\varrho {\omega }^{2}z=:{\mathrm{\Delta }}^{\ast }z+\varrho {\omega }^{2}z=0.$(1) ) in D by the continuous superposition (12) $x\mapsto {\int }_{\mathrm{\Gamma }}{\mathrm{\Phi }}_{\varrho {\omega }^2}(x,y)c\left(y\right)\mathrm{d}s\left(y\right)=:\left({\mathrm{\Phi }}_{\varrho {\omega }^2}{\ast }_{{|}_{\mathrm{\Gamma }}}c\right)\left(x\right),$(12) which is the starting point for the MFS. We now refine (13) $\mathcal{H}\left(\omega \right):=\left\{u={\mathrm{\Phi }}_{{\omega }^2}{\ast }_{{|}_{\mathrm{\Gamma }}}{c}_u,v={\mathrm{\Phi }}_{\rho {\omega }^2}{\ast }_{{|}_{\mathrm{\Gamma }}}{c}_v:(c_u,c_v)\in L^2\left(\mathrm{\Gamma }\right)\times L^2\left(\mathrm{\Gamma }\right)\right\}$(13) from (4(4) $\mathcal{H}:=\bigcup _{0\le \arg \left(\omega \right)<\frac{\pi }{4}}\mathcal{H}\left(\omega \right),$(4) ) and prove in the following that the resulting approximation space is still sufficiently dense to finally recover exact eigenfunctions via boundary control. For this, we need some prerequisites concerning the construction of fundamental solutions for higher order PDEs. It also shows that the MFS ansatz in (12(12) $x\mapsto {\int }_{\mathrm{\Gamma }}{\mathrm{\Phi }}_{\varrho {\omega }^2}(x,y)c\left(y\right)\mathrm{d}s\left(y\right)=:\left({\mathrm{\Phi }}_{\varrho {\omega }^2}{\ast }_{{|}_{\mathrm{\Gamma }}}c\right)\left(x\right),$(12) ) for approximating u and v separately is equivalent to recovering their difference u−v solving (14) $({\mathrm{\Delta }}^{\ast }+{\omega }^2)({\mathrm{\Delta }}^{\ast }+\rho {\omega }^2)(u-v)=0$(14) via its synthesized fundamental solution.

Lemma 4.1

If ${\mathrm{\Phi }}_{\rho {\omega }^2}$ and ${\mathrm{\Phi }}_{{\omega }^2}$ are fundamental solutions for the free space Navier equation (1(1) $\mathrm{d}\mathrm{i}\mathrm{v}\left(\sigma \right(z\left)\right)+\varrho {\omega }^{2}z=(\mu \mathrm{\Delta }z+(\lambda +\mu )\mathrm{\nabla }\left(\mathrm{d}\mathrm{i}\mathrm{v}z\right))+\varrho {\omega }^{2}z=:{\mathrm{\Delta }}^{\ast }z+\varrho {\omega }^{2}z=0.$(1) ) with densities ρ and 1, respectively, then the function ${\mathrm{\Phi }}_{\rho {\omega }^2,{\omega }^2}:=({\mathrm{\Phi }}_{\rho {\omega }^2}-{\mathrm{\Phi }}_{{\omega }^2})/\left(\right(1-\rho \left){\omega }^2\right)$ is a fundamental solution for the fourth order operator $({\mathrm{\Delta }}^{\ast }+{\omega }^2)({\mathrm{\Delta }}^{\ast }+\rho {\omega }^2)$.

Proof.

Let $\phi \in C_c^{\mathrm{\infty }}(D,{\mathbb{C}}^2)$ be an arbitrary bump function and set ${\mathrm{\Phi }}_{\rho {\omega }^2,{\omega }^2}^y\left(x\right):={\mathrm{\Phi }}_{\rho {\omega }^2,{\omega }^2}(y-x)$ for some fixed $y\in {\mathbb{R}}^2$ (and likewise for its generating fundamental solutions). Then we can easily check $\begin{array}{rl}& {\int }_D{\mathrm{\Phi }}_{\rho {\omega }^2,{\omega }^2}^y({\mathrm{\Delta }}^{\ast }+{\omega }^2)({\mathrm{\Delta }}^{\ast }+\rho {\omega }^2)\cdot \phi \mathrm{d}x\\ & ={\int }_D\frac{{\mathrm{\Phi }}_{\rho {\omega }^2}^y-{\mathrm{\Phi }}_{{\omega }^2}^y}{(1-\rho ){\omega }^2}({\mathrm{\Delta }}^{\ast }+{\omega }^2)({\mathrm{\Delta }}^{\ast }+\rho {\omega }^2)\cdot \phi \mathrm{d}x\\ & =\frac{1}{(1-\rho ){\omega }^2}{\int }_D{\mathrm{\Phi }}_{\rho {\omega }^2}^y({\mathrm{\Delta }}^{\ast }+\rho {\omega }^2)({\mathrm{\Delta }}^{\ast }+{\omega }^2)\cdot \phi \mathrm{d}x\\ & -\frac{1}{(1-\rho ){\omega }^2}{\int }_D{\mathrm{\Phi }}_{{\omega }^2}^y({\mathrm{\Delta }}^{\ast }+{\omega }^2)({\mathrm{\Delta }}^{\ast }+\rho {\omega }^2)\cdot \phi \mathrm{d}x\\ & =\frac{{\mathrm{\Delta }}^{\ast }\phi \left(y\right)+{\omega }^2\phi \left(y\right)}{(1-\rho ){\omega }^2}-\frac{{\mathrm{\Delta }}^{\ast }\phi \left(y\right)+\rho {\omega }^2\phi \left(y\right)}{(1-\rho ){\omega }^2}\\ & =\phi \left(y\right),\end{array}$ which proves the lemma.

We are now ready to prove that for any eigenfunction pair $(u,v)$ of (2(2) $\begin{array}{rl}& {\mathrm{\Delta }}^{\ast }u+{\omega }^{2}u=0\mathrm{i}\mathrm{n}D,\\ & {\mathrm{\Delta }}^{\ast }v+\rho {\omega }^{2}v=0\mathrm{i}\mathrm{n}D,\\ & u=v\mathrm{o}\mathrm{n}\mathrm{\partial }D,\\ & \sigma \left(u\right)\nu =\sigma \left(v\right)\nu \mathrm{o}\mathrm{n}\mathrm{\partial }D,\end{array}$(2) ) with eigenfrequency ω lying in the first complex quadrant we can find approximations in $\mathcal{H}\left(\omega \right)$ with arbitrarily small boundary misfits. In particular, all the recovery criteria from Theorem 3.2 can be satisfied for corresponding ITEs.

Theorem 4.2

Let $w\in H^2(D,{\mathbb{C}}^2)$ be any distributional solution to the fourth order equation (14(14) $({\mathrm{\Delta }}^{\ast }+{\omega }^2)({\mathrm{\Delta }}^{\ast }+\rho {\omega }^2)(u-v)=0$(14) ) for some frequency ω such that $0\le \arg \left(\omega \right)<\pi /4$ (holds even more generally for $\mathrm{I}\mathrm{m}\omega \ge 0)$. Then there exists a sequence of elements $(u_m,v_m{)}_{m\in \mathbb{N}}\subset \mathcal{H}(\omega )$ from (13(13) $\mathcal{H}\left(\omega \right):=\left\{u={\mathrm{\Phi }}_{{\omega }^2}{\ast }_{{|}_{\mathrm{\Gamma }}}{c}_u,v={\mathrm{\Phi }}_{\rho {\omega }^2}{\ast }_{{|}_{\mathrm{\Gamma }}}{c}_v:(c_u,c_v)\in L^2\left(\mathrm{\Gamma }\right)\times L^2\left(\mathrm{\Gamma }\right)\right\}$(13) ) such that $(u_m-v_m)=:w_m\to w$ in $H^2(D,{\mathbb{C}}^2)$. If $\mathrm{\partial }D$ is of class $C^{1,1},$ then in particular $\parallel w_m-w{\parallel }_{H^{\frac{3}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}\to 0$ and $\parallel \sigma (w_m-w)\nu {\parallel }_{H^{\frac{1}{2}}(\mathrm{\partial }D,{\mathbb{C}}^2)}\to 0$.

Proof.

Assume $\tilde{w}\in {\tilde{H}}^{-2}(D,{\mathbb{C}}^2)$, where the latter denotes the negative Sobolev space with compact support in D identifying the dual space of $H^2(D,{\mathbb{C}}^2)$, is chosen such that (15) ${\int }_D\tilde{w}\cdot ({\mathrm{\Phi }}_{\rho {\omega }^2,{\omega }^2}{\ast }_{{|}_{\mathrm{\Gamma }}}{c}_1+{\mathrm{\Phi }}_{\rho {\omega }^2}{\ast }_{{|}_{\mathrm{\Gamma }}}{c}_2)\mathrm{d}x=0$(15) for all $(c_1,c_2)\in L^2(\mathrm{\Gamma },{\mathbb{C}}^2)\times L^2(\mathrm{\Gamma },{\mathbb{C}}^2)$. The integral expression is hereby overloaded in notation with the corresponding duality pairing and by definition of ${\mathrm{\Phi }}_{\rho {\omega }^2,{\omega }^2}$ from the previous lemma, we see that the kernel of $\tilde{w}$ contains functions of the form $(u-v)$ with $(u,v)\in \mathcal{H}\left(\omega \right)$. Therefore, if we can show that (15(15) ${\int }_D\tilde{w}\cdot ({\mathrm{\Phi }}_{\rho {\omega }^2,{\omega }^2}{\ast }_{{|}_{\mathrm{\Gamma }}}{c}_1+{\mathrm{\Phi }}_{\rho {\omega }^2}{\ast }_{{|}_{\mathrm{\Gamma }}}{c}_2)\mathrm{d}x=0$(15) ) implies (16) ${\int }_D\tilde{w}\cdot w^{\ast }\mathrm{d}x=0$(16) for every distributional solution $w^{\ast }\in H^2(D,{\mathbb{C}}^2)$ of $({\mathrm{\Delta }}^{\ast }+{\omega }^2)({\mathrm{\Delta }}^{\ast }+\rho {\omega }^2)w^{\ast }=0$, the Hahn–Banach theorem would yield the desired density claim since no other extension is possible.

For this, we define the auxiliary functions $w:={\mathrm{\Phi }}_{\rho {\omega }^2,{\omega }^2}{\ast }_{{|}_D}\tilde{w}\in L^2(D,{\mathbb{C}}^2)\cap C^{\mathrm{\infty }}({\mathbb{R}}^2\mathrm{\setminus }D,{\mathbb{C}}^2)$ and $v:={\mathrm{\Phi }}_{\rho {\omega }^2}{\ast }_{{|}_D}\tilde{w}\in L^2(D,{\mathbb{C}}^2)\cap C^{\mathrm{\infty }}({\mathbb{R}}^2\mathrm{\setminus }D,{\mathbb{C}}^2)$, where the gain of regularity in D results from the fact that the convolution with ${\mathrm{\Phi }}_{\rho {\omega }^2,{\omega }^2}$ or ${\mathrm{\Phi }}_{\rho {\omega }^2}$ are pseudo-differential operators of order $-2$ due to their logarithmic singularity type. Rewriting (15(15) ${\int }_D\tilde{w}\cdot ({\mathrm{\Phi }}_{\rho {\omega }^2,{\omega }^2}{\ast }_{{|}_{\mathrm{\Gamma }}}{c}_1+{\mathrm{\Phi }}_{\rho {\omega }^2}{\ast }_{{|}_{\mathrm{\Gamma }}}{c}_2)\mathrm{d}x=0$(15) ), we obtain ${\int }_{\mathrm{\Gamma }}{c}_1\left(y\right)\cdot w\left(y\right)+c_2\left(y\right)\cdot v\left(y\right)\mathrm{d}s\left(y\right)=0$ for all $(c_1,c_2)\in L^2(\mathrm{\Gamma },{\mathbb{C}}^2)\times L^2(\mathrm{\Gamma },{\mathbb{C}}^2)$ and setting one of the coefficient functions to zero, respectively, we may conclude that $w_{|\mathrm{\Gamma }}=v_{|\mathrm{\Gamma }}=0$. Using the pointwise estimate $|\sqrt{r}({\mathrm{\partial }}_r{v}_p-\mathrm{i}(\rho {\omega }^2\left)v_p\right)\left(x\right)|\le \sqrt{r}\parallel {\mathrm{\partial }}_r({\mathrm{\Phi }}_{\rho {\omega }^2}{)}_p-\mathrm{i}\rho {\omega }^2({\mathrm{\Phi }}_{\rho {\omega }^2}{)}_p{\parallel }_{H^2(x-D,{\mathbb{C}}^{2\times 2})}\parallel \tilde{w}{\parallel }_{H^{-2}(D,{\mathbb{C}}^2)}$, and similarly for $v_s$, while employing standard differentiation properties and decay estimates for the resulting Hankel functions expansion within the first norm on the right, we deduce that Kupradze's radiation conditions are completely inherited by v as convolution of the correspondingly radiating fundamental solution and some compactly supported distribution. By uniqueness of the exterior Navier problem for $\mathrm{I}\mathrm{m}\omega \ge 0$ and $({\mathrm{\Delta }}^{\ast }+\rho {\omega }^2)v=\tilde{w}$ in the sense of distributions, we may then conclude that v = 0 outside of Γ. Due to analyticity, v even needs to vanish completely in ${\mathbb{R}}^2\mathrm{\setminus }\overline{D}$ because the right hand side of the Navier equation is zero here by assumption. Similarly, we want to prove in the following that $w\in H_0^2\left(D\right)$ for justifying its role as a valid test function later:

Using the definition of ${\mathrm{\Phi }}_{\rho {\omega }^2,{\omega }^2}$, direct calculations show the distributional relations $({\mathrm{\Delta }}^{\ast }+{\omega }^2)w=v=({\mathrm{\Delta }}^{\ast }+{\omega }^2)\left({\mathrm{\Phi }}_{{\omega }^2}{\ast }_{{|}_D}v\right)$ as well as $({\mathrm{\Delta }}^{\ast }+\rho {\omega }^2)w={\mathrm{\Phi }}_{{\omega }^2}{\ast }_{{|}_D}\tilde{w}=({\mathrm{\Delta }}^{\ast }+\rho {\omega }^2)\left({\mathrm{\Phi }}_{{\omega }^2}{\ast }_{{|}_D}v\right)$. This combines to $0=\left(\right({\mathrm{\Delta }}^{\ast }+\rho {\omega }^2)-({\mathrm{\Delta }}^{\ast }+{\omega }^2\left)\right)(w-{\mathrm{\Phi }}_{{\omega }^2}{\ast }_{{|}_D}v)/(\rho -1)=w-{\mathrm{\Phi }}_{{\omega }^2}{\ast }_{{|}_D}v$ and implies the additional representation $w={\mathrm{\Phi }}_{{\omega }^2}{\ast }_{{|}_D}v$. The same uniqueness and analyticity reasoning as for v above but with frequency ${\omega }^2$ now yields w = 0 in ${\mathbb{R}}^2\mathrm{\setminus }\overline{D}$ again and by a bootstrap argument due to $v\in L^2(D,{\mathbb{C}}^2)$ we then conclude $w\in H_{loc}^2({\mathbb{R}}^2,{\mathbb{C}}^2)$ or equivalently $w\in H_0^2\left(D\right)$. Therefore, we can find a sequence of bump functions $\{{\phi }_k{\}}_{k\in \mathbb{N}}$ such that $({\mathrm{\Delta }}^{\ast }+\rho {\omega }^2)({\mathrm{\Delta }}^{\ast }+{\omega }^2)(w-{\phi }_k)\to 0$ in $H^{-2}(D,{\mathbb{C}}^2)$. Taking then any distributional solution $w^{\ast }\in H^2(D,{\mathbb{C}}^2)$ of $({\mathrm{\Delta }}^{\ast }+{\omega }^2)({\mathrm{\Delta }}^{\ast }+\rho {\omega }^2)w^{\ast }=0$ as mentioned in (16(16) ${\int }_D\tilde{w}\cdot w^{\ast }\mathrm{d}x=0$(16) ), we may finally compute $\begin{array}{rl}{\int }_D\tilde{w}\cdot w^{\ast }\mathrm{d}x& ={\int }_D\left(({\mathrm{\Delta }}^{\ast }+\rho {\omega }^2)({\mathrm{\Delta }}^{\ast }+{\omega }^2)w\right)\cdot w^{\ast }\mathrm{d}x\\ & =\underset{k\to \mathrm{\infty }}{lim}{\int }_D({\mathrm{\Delta }}^{\ast }+\rho {\omega }^2)({\mathrm{\Delta }}^{\ast }+{\omega }^2){\phi }_k\cdot w^{\ast }\mathrm{d}x\\ & =0,\end{array}$ where we especially incorporated in the first step that $({\mathrm{\Delta }}^{\ast }+\rho {\omega }^2)({\mathrm{\Delta }}^{\ast }+{\omega }^2)w=\tilde{w}$ in the sense of distributions according to our original definition of w. Since $w^{\ast }$ was an arbitrary homogeneous solution, the desired density result for the interior domain is thereby proven. An application of standard trace theorems eventually takes the approximation result over to the boundary of D in corresponding norms.

In the next section, we will focus on the numerical implications of these abstract findings and especially present its associated recipe of how to approximate ITEs in practice.

4.2. The numerical setting

The MFS aims to discretize (12(12) $x\mapsto {\int }_{\mathrm{\Gamma }}{\mathrm{\Phi }}_{\varrho {\omega }^2}(x,y)c\left(y\right)\mathrm{d}s\left(y\right)=:\left({\mathrm{\Phi }}_{\varrho {\omega }^2}{\ast }_{{|}_{\mathrm{\Gamma }}}c\right)\left(x\right),$(12) ) in terms of a Riemann sum with unknown integration weights determined in some optimization process which we derive in the following. For arbitrary frequencies ω and some approximation order $m\in \mathbb{N}$, we focus on trial functions defined on $\overline{D}$ of the form $u_m\left(x\right)=\sum _{j=1}^m{\mathrm{\Phi }}_{{\omega }^2}(x,y_j)({\tilde{c}}_u{)}_{1\le i\le 2,j},v_m(x)=\sum _{j=1}^m{\mathrm{\Phi }}_{\rho {\omega }^2}(x,y_j\left)\right({\tilde{c}}_v{)}_{1\le i\le 2,j},$ where $\{y_1,\dots ,y_m\}\subset \mathrm{\Gamma }$ are preselected source points and ${\tilde{c}}_u,{\tilde{c}}_v\in {\mathbb{C}}^{2\times m}$ are up to now unspecified block coefficient vectors. Note $(u_m,v_m)$ are actually pseudo-elements of $\mathcal{H}\left(\omega \right)$ by using delta-distribution-type coefficient functions $c_u$ and $c_v$ along Γ, respectively. However, it can be readily seen that they are dense in $\mathcal{H}\left(\omega \right)$ for $m\to \mathrm{\infty }$ in any interior Sobolev norm, cf. [9] where this was shown for the acoustic framework. In order to characterize the optimal pair $({\tilde{c}}_u,{\tilde{c}}_v)\in {\mathbb{C}}^{2\times m}\times {\mathbb{C}}^{2\times m}$, we will now develop an MFS-based scheme that is supposed to meet the criteria listed in Theorem 3.2 from a numerical perspective. It uses the collocation concept as a discretized representative for the Sobolev norm-conditioned quantities involved.

First, we treat the vanishing boundary misfit condition $\left(iii\right)$. Therefore we pick m collocation points $\{x_1,\dots ,x_m\}\subset \mathrm{\partial }D$ and consider for fixed ω the constrained minimization of $({\tilde{c}}_u,{\tilde{c}}_v)\mapsto \sum _{i=1}^m|u_m\left(x_i\right)-v_m\left(x_i\right){|}^2.$ This optimization must be subject to criterion $\left(ii\right)$ in a coherently discretized way for which we will take quadrature-like sample points $\{{\hat{x}}_1,\dots ,{\hat{x}}_{m_I}\}\subset D$ with $m_I$ being large but fixed and demand $({\tilde{c}}_u,{\tilde{c}}_v)\mapsto \sum _{i=1}^{m_I}|u_m\left({\hat{x}}_i\right){|}^2+|v_m\left({\hat{x}}_i\right){|}^2\approx 1.$ That is, we scale all eigenfunction candidate pairs to some almost constant but non-vanishing discretized interior norm to guarantee $(u_m,v_m)\ne 0$ for all ω. Both sums above can be restated in compact matrix form, for which we define (17) $\begin{array}{r}M\left(\omega \right):=\left(\begin{array}{cc}{\tilde{\mathrm{\Phi }}}_{{\omega }^2}& {\tilde{\mathrm{\Phi }}}_{\rho {\omega }^2}\\ \sigma \left({\tilde{\mathrm{\Phi }}}_{{\omega }^2}\right)\nu & \sigma \left({\tilde{\mathrm{\Phi }}}_{\rho {\omega }^2}\right)\nu \\ {\hat{\mathrm{\Phi }}}_{{\omega }^2}& 0\\ 0& {\hat{\mathrm{\Phi }}}_{\rho {\omega }^2}\end{array}\right)\in {\mathbb{C}}^{(4m+4m_I)\times 4m},\end{array}$(17) where the upper dense rows correspond to the boundary control matrices $\begin{array}{rl}{\left({\tilde{\mathrm{\Phi }}}_{\varrho {\omega }^2}\right)}_{2i-1,2j-1}& =({\mathrm{\Phi }}_{\varrho {\omega }^2}(x_i,y_j){)}_{1,1},{\left({\tilde{\mathrm{\Phi }}}_{\varrho {\omega }^2}\right)}_{2i-1,2j}=({\mathrm{\Phi }}_{\varrho {\omega }^2}(x_i,y_j){)}_{1,2},\\ {\left({\tilde{\mathrm{\Phi }}}_{\varrho {\omega }^2}\right)}_{2i,2j-1}& =({\mathrm{\Phi }}_{\varrho {\omega }^2}(x_i,y_j){)}_{2,1},{\left({\tilde{\mathrm{\Phi }}}_{\varrho {\omega }^2}\right)}_{2i,2j}=({\mathrm{\Phi }}_{\varrho {\omega }^2}(x_i,y_j){)}_{2,2}\end{array}$ with $1\le i,j\le m$ and co-normal derivatives affecting only the first arguments on the right hand side of involved definitions. A similar pattern applies for $\sigma \left({\tilde{\mathrm{\Phi }}}_{\varrho {\omega }^2}\right)\nu$. The diagonal lower part of (17(17) $\begin{array}{r}M\left(\omega \right):=\left(\begin{array}{cc}{\tilde{\mathrm{\Phi }}}_{{\omega }^2}& {\tilde{\mathrm{\Phi }}}_{\rho {\omega }^2}\\ \sigma \left({\tilde{\mathrm{\Phi }}}_{{\omega }^2}\right)\nu & \sigma \left({\tilde{\mathrm{\Phi }}}_{\rho {\omega }^2}\right)\nu \\ {\hat{\mathrm{\Phi }}}_{{\omega }^2}& 0\\ 0& {\hat{\mathrm{\Phi }}}_{\rho {\omega }^2}\end{array}\right)\in {\mathbb{C}}^{(4m+4m_I)\times 4m},\end{array}$(17) ) then similarly embodies the interior samples $\begin{array}{rl}{\left({\hat{\mathrm{\Phi }}}_{\varrho {\omega }^2}\right)}_{2i-1,2j-1}& =({\mathrm{\Phi }}_{\varrho {\omega }^2}({\hat{x}}_i,y_j){)}_{1,1},{\left({\hat{\mathrm{\Phi }}}_{\varrho {\omega }^2}\right)}_{2i-1,2j}=({\mathrm{\Phi }}_{\varrho {\omega }^2}({\hat{x}}_i,y_j){)}_{1,2},\\ {\left({\hat{\mathrm{\Phi }}}_{\varrho {\omega }^2}\right)}_{2i,2j-1}& =({\mathrm{\Phi }}_{\varrho {\omega }^2}({\hat{x}}_i,y_j){)}_{2,1},{\left({\hat{\mathrm{\Phi }}}_{\varrho {\omega }^2}\right)}_{2i,2j}=({\mathrm{\Phi }}_{\varrho {\omega }^2}({\hat{x}}_i,y_j){)}_{2,2},\end{array}$ where now $1\le i\le m_I$ and again $1\le j\le m$. The benefit of this matrix reformulation is that we can now perform a QR-factorization $M\left(\omega \right)=Q_M\left(\omega \right)R_M\left(\omega \right)=\left(\begin{array}{c}Q\left(\omega \right)\\ Q_I\left(\omega \right)\end{array}\right)R_M\left(\omega \right)$ with $Q\in {\mathbb{C}}^{4m\times 4m}$ and $Q_I\in {\mathbb{C}}^{4m_I\times 4m}$, which conveniently reflects the above minimal-boundary-to-maximal-interior coupling through the unitary property of $Q_M$. This observation goes back to Betcke and Trefethen, see [19], who analyzed Dirichlet eigenvalues with a slightly different ansatz but still via boundary control. Our common ingredient left is then to find local minima of (18) $\omega \mapsto \underset{r\in {\mathbb{C}}^{4m},|r|=1}{min}|Q\left(\omega \right)r|=:{\sigma }_1\left(\omega \right).$(18) Those solutions that are (almost) roots for the smallest singular value ${\sigma }_1\left(\omega \right)$, see Figure 1 (left), will be denoted by ${\omega }_m$ to trace back to the underlying trial space dimension of $(u_m,v_m)$. Our described solution algorithm will be called modified MFS and its output, ${\omega }_m$, approximate ITEs. In the spirit of $\left(i\right)$ from Theorem 3.2 we hope them to converge for $m\to \mathrm{\infty }$ also in practice and Lemma 3.3 would then provide a measure for their convergence speed in terms of the continuous analogon of ${\sigma }_1\left({\omega }_m\right)$. The next section will demonstrate this for some exemplary scatterers.

Figure 1. Exemplary plots of (18(18) $\omega \mapsto \underset{r\in {\mathbb{C}}^{4m},|r|=1}{min}|Q\left(\omega \right)r|=:{\sigma }_1\left(\omega \right).$(18) ) along the real axis for the disc of radius 0.5 as scatterer generated for m = 50 (left) and m = 80 (right). Both cases confirm 3 approximate ITEs ${\omega }_{50}$ within the interval $[0,2]$. However, for critically large m the graph gets polluted by ill-conditioning artifacts.

5. Numerical results

We use the modified MFS from the previous section to compute some ITEs for a disc-, ellipse-, kite- and square-shaped scatterer D whose (collocation) boundaries are given by $\begin{array}{rl}\mathrm{\partial }{D}_d& :=\left(\begin{array}{c}0.5\cos \left(t\right)\\ 0.5\sin \left(t\right)\end{array}\right),t\in [0,2\pi ),\\ \mathrm{\partial }{D}_e& :=\left(\begin{array}{c}0.5\cos \left(t\right)\\ \sin \left(t\right)\end{array}\right),t\in [0,2\pi ),\\ \mathrm{\partial }{D}_k& :=\left(\begin{array}{c}0.75\cos \left(t\right)+0.3\cos \left(2t\right)\\ \sin \left(t\right)\end{array}\right),t\in [0,2\pi ),\\ \mathrm{\partial }{D}_s& :=\mathrm{\partial }\left([-0.5,0.5]\times [-0.5,0.5]\right),\end{array}$ respectively. Note that $D_s$ is actually not covered by our ITE approximation analysis due to its corners but still a very typical scatterer feasible for our algorithm since collocation is invisible with respect to the boundary regularity except for singular points that need to be excluded. Throughout, we fix the constitutive parameters $\mu =\frac{1}{16},\lambda =\frac{1}{4},\rho =4$ which were also used in the context of [14] and therefore serve as independent reference values for our exemplary MFS-based findings listed in Figure 2. Those approximate the first four real-valued ITEs from each scatterer and were individually obtained via the source point boundaries $\begin{array}{rl}{\mathrm{\Gamma }}_d& :=\left(\begin{array}{c}\cos \left(t\right)\\ \sin \left(t\right)\end{array}\right)=2\cdot \mathrm{\partial }{D}_d,t\in [0,2\pi ),\\ {\mathrm{\Gamma }}_e& :=\left(\begin{array}{c}0.95\cos \left(t\right)\\ 1.9\sin \left(t\right)\end{array}\right)=1.9\cdot \mathrm{\partial }{D}_e,t\in [0,2\pi ),\\ {\mathrm{\Gamma }}_k& :=\left(\begin{array}{c}\left(1.2\cos \left(t\right)+0.48\cos \left(2t\right)\right)\\ 1.6\sin \left(t\right)\end{array}\right)=1.6\cdot \mathrm{\partial }{D}_k,t\in [0,2\pi ),\\ {\mathrm{\Gamma }}_s& :=\mathrm{\partial }\left([-0.65,0.65]\times [-0.65,0.65]\right)=1.3\cdot \mathrm{\partial }{D}_s.\end{array}$

Figure 2. Approximations of the first four real-valued ITEs (counting without multiplicity) for some exemplary scatterers obtained via the modified MFS with material parameter $\mu =1/16,\lambda =1/4,\rho =4$: One clearly sees that the more advanced the shape of the boundary becomes, the less ITE digits can be effectively recovered.

The extracted scaling factors $\{2,1.9,1.6,1.3\}$ were preselected to follow the conclusions from [9] and approach unity the more scattering shapes seem to deviate from the disc that is considered as the most promising state. Having thus set all the necessary computational contours for the modified MFS, both the m-dependent collocation points $\{x_1,\dots ,x_m\}$ as well as the auxiliary sources $\{y_1,\dots ,y_{2m}\}$ were distributed equiangular on their corresponding boundaries, e.g. equidistant with respect to t if our representation above allows. However, it is well known that the individual optimization of Γ itself and its source point distribution do have a noticeable impact on the MFS approximation quality as was shown in [24], but resulting in a more advanced non-linear problem for each scatterer in total. In contrast, concerning the interior points, $m_I=10$ of them were fixed randomly in a centered disc with radius 0.5 as their overall location and number turns out not to affect the ITE output significantly.

With these input arrangements and the focus on real-valued ITEs first, the modified MFS can be successfully exploited by incrementing m within $40\le m\le 80$, where the lower bound just provided an averaged value for when to expect good results. However, exceeding the given D-specific threshold for our setup, the graph of (18(18) $\omega \mapsto \underset{r\in {\mathbb{C}}^{4m},|r|=1}{min}|Q\left(\omega \right)r|=:{\sigma }_1\left(\omega \right).$(18) ) always began to suffer drastically from the intrinsic ill-conditioning effects of the discretization matrix $M\left(\omega \right)$ in (17(17) $\begin{array}{r}M\left(\omega \right):=\left(\begin{array}{cc}{\tilde{\mathrm{\Phi }}}_{{\omega }^2}& {\tilde{\mathrm{\Phi }}}_{\rho {\omega }^2}\\ \sigma \left({\tilde{\mathrm{\Phi }}}_{{\omega }^2}\right)\nu & \sigma \left({\tilde{\mathrm{\Phi }}}_{\rho {\omega }^2}\right)\nu \\ {\hat{\mathrm{\Phi }}}_{{\omega }^2}& 0\\ 0& {\hat{\mathrm{\Phi }}}_{\rho {\omega }^2}\end{array}\right)\in {\mathbb{C}}^{(4m+4m_I)\times 4m},\end{array}$(17) ) via impeding oscillations and finally lead to unreliable approximations apart, see Figure 1 (right). Obtained via the regime in between, we believe our given cut-off results from Figure 2 to be correct up to that point modulo round-off errors compared to the exact ITE mantissa since they correspond to the nonfluctuating digits within the modified MFS output ${\omega }_m$ when $m\nearrow 80$. Some exact reference values for $D_d$ corresponding to rotational-symmetric eigenfunctions can be obtained by computing roots of (9(9) $\mathrm{d}\mathrm{e}\mathrm{t}\left(\begin{array}{cc}{J}_1\left(\omega \sqrt{\frac{1}{\mu }}\right)& J_1\left(\omega \sqrt{\frac{\rho }{\mu }}\right)\\ \omega \sqrt{\frac{1}{\mu }}{J}_1^{\prime }\left(\omega \sqrt{\frac{1}{\mu }}\right)& \omega \sqrt{\frac{\rho }{\mu }}{J}_1^{\prime }\left(\omega \sqrt{\frac{\rho }{\mu }}\right)\end{array}\right)=0,$(9) ) which even confirms our approximation for its smallest ITE to be exact up to machine precision, see Figure 3 for the convergence history. In general, the recoverable accuracy strongly correlates with the scattering shape and deteriorates with the existence of corners or concave parts. These technical observations mostly agree with the ones from the acoustic case analyzed in [9,10] and only differ in the retarded yet D-specific convergence regime in m for starting the actual ITE approximations which was about $m\ge 20$ before. It manifests the fact that the eigenfunctions from elasticity are vector-valued and thus numerically slightly more expensive.As in the particular scope of this paper, we also want to present some approximations from the complex-valued eigenvalue spectrum. Therefore, Figure 4 plots (18(18) $\omega \mapsto \underset{r\in {\mathbb{C}}^{4m},|r|=1}{min}|Q\left(\omega \right)r|=:{\sigma }_1\left(\omega \right).$(18) ) in the complex range $0\le \mathrm{R}\mathrm{e}\omega \le 2.5,-2\le \mathrm{I}\mathrm{m}\omega \le 2$ for the disc $D_d$ and the square $D_s$ since only their distribution of ITEs seems not that dense among our analyzed scatterers according to Figure 2 and is thus better suited for broad contour visualizations. We can immediately see that ITEs are correctly displayed in conjugated pairs and as an example the closest of them with respect to the positive imaginary axis are computed for the two scatterers to be $1.8624+0.3104\mathrm{i}$ and $1.987178187576699+0.283125784408650\mathrm{i}$, respectively. The correctness of all given ITE digits for $D_d$ can even be confirmed again by comparing with the correspondingly coinciding root of (9(9) $\mathrm{d}\mathrm{e}\mathrm{t}\left(\begin{array}{cc}{J}_1\left(\omega \sqrt{\frac{1}{\mu }}\right)& J_1\left(\omega \sqrt{\frac{\rho }{\mu }}\right)\\ \omega \sqrt{\frac{1}{\mu }}{J}_1^{\prime }\left(\omega \sqrt{\frac{1}{\mu }}\right)& \omega \sqrt{\frac{\rho }{\mu }}{J}_1^{\prime }\left(\omega \sqrt{\frac{\rho }{\mu }}\right)\end{array}\right)=0,$(9) ). Regarding accuracy, the modified MFS does not show any remarkable difference for its extension to the upper half space of the complex plane. However, although the plot indicates a certain axial symmetry with below, one clearly sees that inaccuracies due to large condition numbers of the underlying matrix (17(17) $\begin{array}{r}M\left(\omega \right):=\left(\begin{array}{cc}{\tilde{\mathrm{\Phi }}}_{{\omega }^2}& {\tilde{\mathrm{\Phi }}}_{\rho {\omega }^2}\\ \sigma \left({\tilde{\mathrm{\Phi }}}_{{\omega }^2}\right)\nu & \sigma \left({\tilde{\mathrm{\Phi }}}_{\rho {\omega }^2}\right)\nu \\ {\hat{\mathrm{\Phi }}}_{{\omega }^2}& 0\\ 0& {\hat{\mathrm{\Phi }}}_{\rho {\omega }^2}\end{array}\right)\in {\mathbb{C}}^{(4m+4m_I)\times 4m},\end{array}$(17) ) dominantly propagate from the lower half space upwards for increasing m, giving a first hint where the oscillations in Figure 1 (right) originate from. Fortunately, ITEs always arise in conjugated pairs and because of Theorem 4.2 we even know that restricting to complex numbers with non-negative imaginary part within our investigations, as already implemented in our initial ansatz (4(4) $\mathcal{H}:=\bigcup _{0\le \arg \left(\omega \right)<\frac{\pi }{4}}\mathcal{H}\left(\omega \right),$(4) ) and (13(13) $\mathcal{H}\left(\omega \right):=\left\{u={\mathrm{\Phi }}_{{\omega }^2}{\ast }_{{|}_{\mathrm{\Gamma }}}{c}_u,v={\mathrm{\Phi }}_{\rho {\omega }^2}{\ast }_{{|}_{\mathrm{\Gamma }}}{c}_v:(c_u,c_v)\in L^2\left(\mathrm{\Gamma }\right)\times L^2\left(\mathrm{\Gamma }\right)\right\}$(13) ), is not limiting at all.

Figure 3. Exponential decay of the modified MFS with respect to the absolute ITE deviation from 1.451304027606383 being the smallest real-valued one of the disc with radius 0.5 as scatterer. A similar behavior for the boundary data misfit of the corresponding approximate eigenfunctions in terms of the smallest singular value is shown. Their linear correlation in the plot confirms the continuous analogon from Lemma 3.3.

To get a rough impression about how competitive our method is compared to others, we use the approximation results for real-valued $k^2$ from the paper [14] in which the authors also analyzed the scatterers $D_d$ and $D_s$ but via special finite element methods. With a mesh size of $h\le 0.0125$, they obtained the approximations $k_d^2\approx 2.110723$ and $k_s^2\approx 1.9428781$ for the smallest ITE squared of the disc and the square, respectively. Allowing for possible round-off deviations in the last digit within our corresponding cut-off values, we obtain the tolerance ranges $k_d^2\in [2.106283380546506,2.106283380546512]$ and $k_s^2\in [1.9423,1.9429]$ by setting at most m = 80 which may be converted via a boundary partitioning to a collocation point distance of the same order as the alternative mesh size. While the results regarding $D_s$ are of the same order for both solution approaches, the modified MFS clearly dominates in accuracy for $D_d$ and is therefore believed to do so for slightly perturbed scatterers D, too.

Figure 4. Contour extract from the complex plane evaluating (18(18) $\omega \mapsto \underset{r\in {\mathbb{C}}^{4m},|r|=1}{min}|Q\left(\omega \right)r|=:{\sigma }_1\left(\omega \right).$(18) ) for the unit square (below) and the disc of radius 0.5 (above) as scatterers with m = 45 (left) and m = 55 (right), respectively. Approximate ITEs are denoted by the centers of emerging concentric circles and spoilt primarily from the lower half space through the ill-conditioning artifacts already encountered in Figure 1.

6. Conclusion

We analyzed the method of fundamental solution in a stabilized version for the computation of complex-valued elastic transmission eigenvalues in two dimensions. Our theoretical studies show that the short algorithm shall produce spurious-free results in the limit whose approximation error per step we could quantify in terms of some discretized residual output value. Our numerical experiments confirm these expectations in practice within some scatterer-specific collocation point regime for which the method's feasibility is unaffected from ill-conditioning pollutions. This was tested for a collection of scatterers whose boundaries are easily-parametrizable. In accordance with the conclusions from previous works that analyzed our algorithm in the context of related transmission problems but restricted to real-valued eigenvalues so far, the best results, including also the complex spectrum from now, can still be obtained for the unit disc. The more the actual scattering shape then deviates from the disc, the less accurate approximations are finally achievable while even more collocation points are needed in total. However, depending on the choice of fundamental solution for generating the trial functions, different areas of the complex spectrum (in our case the lower half space) should be avoided as initial guess input for the algorithm due to ill-conditioning effects which are, however, not that restrictive because all eigenvalues arise in conjugated pairs. Generally, since the source boundary necessary for the MFS setup was individually preselected by intuition for simplicity, even more promising approximations can be achieved by its optimization which was, however, not covered in the scope of this paper. As a conclusion, the modified MFS is especially effective for regular convex domains and dominates here over many competitive methods in the general context of transmission eigenvalue problems.

Disclosure statement

No potential conflict of interest was reported by the author(s).