An efficient spectral-Galerkin approximation and error analysis for Maxwell transmission eigenvalue problems in spherical geometries

TL;DR

This paper develops an efficient spectral-Galerkin method for Maxwell transmission eigenvalue problems in spherical geometries, utilizing vector spherical harmonics and specialized algorithms for TE and TM modes, with proven error estimates and validated numerical results.

Contribution

It introduces a novel spectral-Galerkin approach with error analysis for Maxwell transmission eigenvalues in spherical geometries, including algorithms for TE and TM modes.

Findings

Spectral-Galerkin method achieves high accuracy for eigenvalues.

Parallelizable solution process for TE and TM modes.

Numerical experiments confirm theoretical error estimates.

Abstract

We propose and analyze an efficient spectral-Galerkin approximation for the Maxwell transmission eigenvalue problem in spherical geometry. Using a vector spherical harmonic expansion, we reduce the problem to a sequence of equivalent one-dimensional TE and TM modes that can be solved individually in parallel. For the TE mode, we derive associated generalized eigenvalue problems and corresponding pole conditions. Then we introduce weighted Sobolev spaces based on the pole condition and prove error estimates for the generalized eigenvalue problem. The TM mode is a coupled system with four unknown functions, which is challenging for numerical calculation. To handle it, we design an effective algorithm using Legendre-type vector basis functions. Finally, we provide some numerical experiments to validate our theoretical results and demonstrate the efficiency of the algorithms.

Equations

$$H(\mathbf{curl},D):=\{\mathbf{u}\in {L^2(D)}^3:\mathbf{curl}\mathbf{u}\in L^2(D{)}^3\},$$

$$H_0(\mathbf{curl},D):=\{\mathbf{u}\in H(\mathbf{curl},D):\mathbf{u}\times \mathbf{\nu }=0\text{on}\partial D\}$$

$$\mathcal{U}(D):=\{\mathbf{u}\in H(\mathbf{curl},D):\mathbf{curl}\mathbf{u}\in H(\mathbf{curl},D)\},$$

$${\mathcal{U}}_0(D):=\{\mathbf{u}\in H_0(\mathbf{curl},D):\mathbf{curl}\mathbf{u}\in H_0(\mathbf{curl},D)\}$$

$$\bar{\xi }\cdot \mathbf{N}\xi \ge \gamma |\xi {|}^2,\forall \xi \in {\mathbb{C}}^3\text{ a.e in }D.$$

$$\mathbf{curl}\mathbf{curl}\mathbf{E}-k^2\mathbf{N}\mathbf{E}=0$$

$$\mathbf{curl}\mathbf{curl}{\mathbf{E}}_0-k^2{\mathbf{E}}_0=0$$

$$\mathbf{E}\times \mathbf{\nu }={\mathbf{E}}_0\times \mathbf{\nu }$$

$$\mathbf{curl}\mathbf{E}\times \mathbf{\nu }=\mathbf{curl}{\mathbf{E}}_0\times \mathbf{\nu }$$

$$x_1=r\sin \theta \cos \varphi ,x_2=r\sin \theta \sin \varphi ,x_3=r\cos \theta .$$

$${\mathbf{e}}_r=\frac{\mathbf{x}}{r},{\mathbf{e}}_{\theta }=(\sin \theta \cos \varphi ,\sin \theta \sin \varphi ,-\sin \theta ),{\mathbf{e}}_{\varphi }=(-\sin \varphi ,\cos \varphi ,0).$$

$${△}_{S}u=\frac{1}{\sin \theta }\frac{\partial }{\partial \theta }(\sin \theta \frac{\partial u}{\partial \theta })+\frac{1}{{\sin }^2\theta }\frac{{\partial }^{2}u}{\partial {\varphi }^2},$$

$${\nabla }_{S}u=\frac{\partial u}{\partial \theta }{\mathbf{e}}_{\theta }+\frac{1}{\sin \theta }\frac{\partial u}{\partial \varphi }{\mathbf{e}}_{\varphi }.$$

$${\Delta }_S{Y}_l^m=-l(l+1)Y_l^m,l\ge 0,|m|\le l;$$

$${\int }_S{Y}_l^m{Y}_{l^{\prime }}^{m^{\prime }}dS={\delta }_{l{l}^{\prime }}{\delta }_{m{m}^{\prime }}.$$

$${\mathbf{T}}_l^m={\nabla }_S{Y}_l^m\times {\mathbf{e}}_r=\frac{1}{\sin \theta }\frac{\partial Y_l^m}{\partial \varphi }{\mathbf{e}}_{\theta }-\frac{\partial Y_l^m}{\partial \theta }{\mathbf{e}}_{\varphi },$$

$${\mathbf{V}}_l^m=(l+1)Y_l^m{\mathbf{e}}_r-{\nabla }_S{Y}_l^m,$$

$${\mathbf{W}}_l^m=l{Y}_l^m{\mathbf{e}}_r+{\nabla }_S{Y}_l^m,$$

$${\int }_S{\mathbf{T}}_l^m\cdot \overline{{\mathbf{V}}_{l^{\prime }}^{m^{\prime }}}dS={\int }_S{\mathbf{T}}_l^m\cdot \overline{{\mathbf{W}}_{l^{\prime }}^{m^{\prime }}}dS={\int }_S{\mathbf{V}}_l^m\cdot \overline{{\mathbf{W}}_{l^{\prime }}^{m^{\prime }}}dS=0,$$

$${\int }_S{\mathbf{V}}_l^m\cdot \overline{{\mathbf{V}}_{l^{\prime }}^{m^{\prime }}}dS=(l+1)(2l+1){\delta }_{l{l}^{\prime }}{\delta }_{m{m}^{\prime }},{\int }_S{\mathbf{T}}_l^m\cdot \overline{{\mathbf{T}}_{l^{\prime }}^{m^{\prime }}}dS=l(l+1){\delta }_{l{l}^{\prime }}{\delta }_{m{m}^{\prime }},$$

$${\int }_S{\mathbf{W}}_l^m\cdot \overline{{\mathbf{W}}_{l^{\prime }}^{m^{\prime }}}dS=l(2l+1){\delta }_{l{l}^{\prime }}{\delta }_{m{m}^{\prime }},{\int }_S{\mathbf{T}}_l^m\cdot {\nabla }_S{Y}_l^{m}dS=0,$$

$${\int }_S{\nabla }_S{Y}_l^m\cdot {\nabla }_S{Y}_{l^{\prime }}^{m^{\prime }}dS=l(l+1){\delta }_{l{l}^{\prime }}{\delta }_{m{m}^{\prime }}.$$

$$d_l^{+}=\frac{d}{dr}+\frac{l}{r},d_l^{-}=\frac{d}{dr}-\frac{l}{r}.$$

$$\mathbf{curl}(f{\mathbf{V}}_l^m)=(d_{l+2}^{+}f){\mathbf{T}}_l^m,\mathbf{curl}(f{\mathbf{W}}_l^m)=-(d_{l-1}^{-}f){\mathbf{T}}_l^m,$$

$$(2l+1)\mathbf{curl}(f{\mathbf{T}}_l^m)=(l+1)(d_{l+1}^{+}f){\mathbf{W}}_l^m-l(d_l^{-}f){\mathbf{V}}_l^m.$$

$$\mathbf{E}=\sum _{l=0}^{\infty }\sum _{|m|=0}^l(v_l^m(r){\mathbf{V}}_l^m(\theta ,\varphi )+t_l^m(r){\mathbf{T}}_l^m(\theta ,\varphi )+w_l^m(r){\mathbf{W}}_l^m(\theta ,\varphi )),$$

$${\mathbf{E}}_0=\sum _{l=0}^{\infty }\sum _{|m|=0}^l({\bar{v}}_l^m(r){\mathbf{V}}_l^m(\theta ,\varphi )+{\bar{t}}_l^m(r){\mathbf{T}}_l^m(\theta ,\varphi )+{\bar{w}}_l^m(r){\mathbf{W}}_l^m(\theta ,\varphi )).$$

$$\mathbf{curl}\mathbf{curl}(v_l^m(r){\mathbf{V}}_l^m(\theta ,\varphi ))=\mathbf{curl}((d_{l+2}^{+}{v}_l^m(r)){\mathbf{T}}_l^m(\theta ,\varphi ))$$

$$=\frac{l+1}{2l+1}(d_{l+1}^{+}{d}_{l+2}^{+}{v}_l^m(r)){\mathbf{W}}_l^m(\theta ,\varphi )-\frac{l}{2l+1}(d_l^{-}{d}_{l+2}^{+}{v}_l^m(r)){\mathbf{V}}_l^m(\theta ,\varphi ),$$

$$\mathbf{curl}\mathbf{curl}(t_l^m(r){\mathbf{T}}_l^m(\theta ,\varphi ))=-(\frac{l}{2l+1}{d}_{l+2}^{+}{d}_l^{-}+\frac{l+1}{2l+1}{d}_{l-1}^{-}{d}_{l+1}^{+})t_l^m(r){\mathbf{T}}_l^m(\theta ,\varphi ),$$

$$\mathbf{curl}\mathbf{curl}(w_l^m(r){\mathbf{W}}_l^m(\theta ,\varphi ))=\frac{l}{2l+1}(d_l^{-}{d}_{l-1}^{-}{w}_l^m(r)){\mathbf{V}}_l^m(\theta ,\varphi )$$

$$-\frac{l+1}{2l+1}(d_{l+1}^{+}{d}_{l-1}^{-}{w}_l^m(r)){\mathbf{W}}_l^m(\theta ,\varphi ).$$

$$\frac{l}{2l+1}(d_l^{-}{d}_{l-1}^{-}{w}_l^m(r)-d_l^{-}{d}_{l+2}^{+}{v}_l^m(r))-k_l^{2}n{v}_l^m(r)=0,$$

$$\frac{l+1}{2l+1}(d_{l+1}^{+}{d}_{l+2}^{+}{v}_l^m(r)-d_{l+1}^{+}{d}_{l-1}^{-}{w}_l^m(r))-k_l^{2}n{w}_l^m(r)=0,$$

$$\frac{-1}{2l+1}(l{d}_{l+2}^{+}{d}_l^{-}+(l+1)d_{l-1}^{-}{d}_{l+1}^{+})t_l^m(r)-k_l^{2}n{t}_l^m(r)=0,$$

$$\frac{l}{2l+1}(d_l^{-}{d}_{l-1}^{-}{\bar{w}}_l^m(r)-d_l^{-}{d}_{l+2}^{+}{\bar{v}}_l^m(r))-k_l^2{\bar{v}}_l^m(r)=0,$$

$$\frac{l+1}{2l+1}(d_{l+1}^{+}{d}_{l+2}^{+}{\bar{v}}_l^m(r)-d_{l+1}^{+}{d}_{l-1}^{-}{\bar{w}}_l^m(r))-k_l^2{\bar{w}}_l^m(r)=0,$$

$$\frac{-1}{2l+1}(l{d}_{l+2}^{+}{d}_l^{-}+(l+1)d_{l-1}^{-}{d}_{l+1}^{+}){\bar{t}}_l^m(r)-k_l^2{\bar{t}}_l^m(r)=0,$$

$$\mathbf{E}\times \mathbf{\nu }=\mathbf{E}\times {\mathbf{e}}_r=\sum _{l=0}^{\infty }\sum _{|m|=0}^l((w_l^m(r)-v_l^m(r)){\mathbf{T}}_l^m-t_l^m(r){\nabla }_S{Y}_l^m).$$

$${\mathbf{E}}_0\times \mathbf{\nu }={\mathbf{E}}_0\times {\mathbf{e}}_r=\sum _{l=0}^{\infty }\sum _{|m|=0}^l(({\bar{w}}_l^m(r)-{\bar{v}}_l^m(r)){\mathbf{T}}_l^m-{\bar{t}}_l^m(r){\nabla }_S{Y}_l^m).$$

$$w_l^m(R)-v_l^m(R)={\bar{w}}_l^m(R)-{\bar{v}}_l^m(R),t_l^m(R)={\bar{t}}_l^m(R).$$

Full text

Keywords: Maxwell transmission eigenvalue problems, spherical geometry, TE and TM modes, spectral-Galerkin approximation, error analysis

1 Introduction

The Maxwell transmission eigenvalue problem is a boundary value problem in a bounded domain which arises in inverse scattering theory for inhomogeneous media. It is well known that the transmission eigenvalues play a critical role in the reconstruction of inhomogeneous non-absorbing media [1, 2, 3]. In addition, the transmission eigenvalues are also used to estimate the index of refraction of a non-absorbing inhomogeneous medium in recent years[4, 5, 6]. The method consists of several steps. First of all, the support of the scattering obstacle can be recovered by using the measured scattering data and the linear sampling method [7] and the transmission eigenvalues can be identified from either the far field or near field data [8, 4]. Then, the bounds for smallest and largest eigenvalues of the (matrix) index of refraction can be obtained in terms of the support of the scattering obstacle and the first transmission eigenvalue of the anisotropic media [9]. Finally, reconstructions of the electric permittivity (if it is a scalar constant) or an estimate of the eigenvalues of the matrix in the case of anisotropic permittivity can be obtained [7].

However, the effectiveness of the above method rests on having an efficient and robust algorithm for computing Maxwell transmission eigenvalues. Although the Maxwell transmission eigenvalue problem is stated in a simple form, its solution is complex since it is nonstandard so that the classical theory can not be applied directly. There are only a few papers dealing with the numerical computation of Maxwell transmission eigenvalues due to this notorious difficulty. In [10], Monk et al. propose two finite element methods in computing a few lowest Maxwell’s transmission eigenvalues, curl-conforming finite element and mixed finite element methods. In [3], Sun et al. presented an iterative method to compute the Maxwell’s transmission eigenvalue, where the transmission eigenvalue problem is written as a quad-curl eigenvalue problem. Then the real transmission eigenvalues are shown to be the roots of a nonlinear function whose value is the generalized eigenvalue of a related self-adjoint quad-curl eigenvalue problem which is computed by using a mixed finite element method. In [11], Huang et al. also proposed a numerical algorithm for computing Maxwell transmission eigenvalue problems. However, all these methods are based on low-order finite element methods, so they become very expensive if high accuracy is needed. To the best of our knowledge, no error analysis is given yet for any existing numerical method, even for three-dimensional spherical domains.

In practical applications, we often need to solve Maxwell transmission eigenvalue problems on the special domain of spherical geometries. We know of only a few reports on spectral-Galerkin approximation for the Maxwell transmission eigenvalue problems on the special domain of spherical geometries. Thus, we will present in this paper an efficient spectral-Galerkin approximation for the Maxwell transmission eigenvalue problem in spherical geometry. Using a vector spherical harmonic expansion, we reduce the problem to a sequence of equivalent one-dimensional TE and TM modes that can be solved individually in parallel. For the TE mode, we derive associated generalized eigenvalue problems and corresponding pole conditions. Then we introduce weighted Sobolev spaces corresponding to the pole condition and prove error estimates for the generalized eigenvalue problem. The TM mode is a coupled system with four unknown functions, which is challenging for numerical calculation. To handle it, we design an effective algorithm using Legendre-type vector basis functions. Finally, we provide some numerical experiments to validate our theoretical results and demonstrate the efficiency of the algorithms.

The rest of this paper is organized as follows. In the next section, we introduce the Maxwell transmission eigenvalue problem. In §3, we derive a dimension reduction scheme under spherical geometries. In §4, we derive the weak formulation and error estimate for the TE mode. In §5, we describe in detail an efficient implementation of the algorithm. We present several numerical experiments in §6 to demonstrate the accuracy and efficiency of our method. Finally, in §7 we give some concluding remarks.

2 Maxwell transmission eigenvalue problem

Let $D\subset\mathbb{R}^{3}$ be a bounded simply connected region with piecewise smooth boundary $\partial D$ and denote by $\boldsymbol{\nu}$ the outward normal vector to $\partial D$. Let $(.,.)_{D}$ denote the scalar product in $L^{2}(D)^{3}$ and define the Hilbert spaces

[TABLE]

equipped with the scalar product $(\mathbf{u},\mathbf{v})_{\mathbf{curl}}=(\mathbf{u},\mathbf{v})_{D}+(\mathbf{curl}\;\mathbf{u},\mathbf{curl}\;\mathbf{v})_{D}$ and the corresponding norm $\|\cdot\|_{\mathbf{curl}}$. Following [12], we also define

[TABLE]

equipped with the scalar product $(\mathbf{u},\mathbf{v})_{\mathcal{U}}=(\mathbf{u},\mathbf{v})_{\mathbf{curl}}+(\mathbf{curl}\;\mathbf{u},\mathbf{curl}\;\mathbf{v})_{\mathbf{curl}}$ and the corresponding norm $\|\cdot\|_{\mathcal{U}}$.

Let $\boldsymbol{N}$ be a $3\times 3$ matrix valued function defined on $D$ such that $\boldsymbol{N}\in L^{\infty}(D,\mathbb{R}^{3\times 3})$. Definition 2.1 A real matrix field $\boldsymbol{N}$ is said to be bounded positive definite on $D$ if $\boldsymbol{N}\in L^{\infty}(D,\mathbb{R}^{3\times 3})$ and there exists a constant $\gamma>0$ such that

[TABLE]

We further assume that $\boldsymbol{N},\boldsymbol{N}^{-1}$ and $(\boldsymbol{N}-\boldsymbol{I})^{-1}$ or $(\boldsymbol{I}-\boldsymbol{N})^{-1}$ are bounded positive definite real matrix fields on $D$. The interior transmission eigenvalue problem for the anisotropic Maxwell’s equations in terms of electric fields is formulated as the problem of finding two vector valued functions $\mathbf{E}\in{L^{2}(D)}^{3}$ and $\mathbf{E_{0}}\in{L^{2}(D)}^{3}$ such that $\mathbf{E}-\mathbf{E_{0}}\in\mathcal{U}_{0}(D)$ satisfies [3, 10]:

[TABLE]

3 Dimension reduction scheme under spherical geometries

We shall restrict our attention to the case where $D$ is a ball of radius $R$ and $\boldsymbol{N}=n\boldsymbol{I}$ with $n$ being a function along the radial direction. In this case, by using vector spherical harmonic expansion we can reduce the problem to a sequence of equivalent one-dimensional TE and TM modes that can be solved individually in parallel.

3.1 Vector spherical harmonics

There are several versions with different notation and properties for vector spherical harmonics (see e.g.,[13, 14, 15, 16]). We adopt in this paper the family of vector spherical harmonics in [13, 14].

The spherical coordinates $(r,\theta,\phi)$ are related to the Cartesian coordinates $\mathbf{x}=(x_{1},x_{2},x_{3})$ by

[TABLE]

Let

[TABLE]

Then $\{\mathbf{e}_{r},\mathbf{e}_{\theta},\mathbf{e}_{\phi}\}$ forms a moving orthonormal coordinate basis in $\mathbb{R}^{3}$.

Let $S$ be the unit spherical surface, and denote by $\Delta_{S}$ and $\nabla_{S}$ the Laplace-Beltrami and tangential gradient operators on $S$, namely,

[TABLE]

The spherical harmonics $\{Y_{l}^{m}\}$ (as normalized in [13]) are eigenfunctions of $\Delta_{S}$, i.e.,

[TABLE]

and form an orthonormal basis for $L^{2}(S)$:

[TABLE]

The family of vector spherical harmonics is defined by [13]:

[TABLE]

where $\mathbf{V}_{0}^{0}=\frac{1}{\sqrt{4\pi}}\mathbf{e}_{r}$. With the understanding of $\mathbf{T}_{0}^{0}=\mathbf{W}_{0}^{0}=\mathbf{0},$ the indexes $\{l,m\}$ run over $\{(l,m):l\geq 0,0\leq|m|\leq l\}.$

The following lemma was proved in [13] ($(A.1)$ and $(A.2)$)

Lemma 3.1

[TABLE]

Let $f$ is a function along the radial direction and define the differentiation operators:

[TABLE]

The following lemma is again proved in [13] ($(A.9)$):

Lemma 3.2

[TABLE]

3.2 Dimension reduction scheme

Let us write

[TABLE]

From Lemma 3.2 we can derive that

[TABLE]

From (3.19)-(3.21) and the nonzero of eigenfunction, together with Lemma 3.1, (2.1) can be reduced to

[TABLE]

where $l\geq 1,0\leq|m|\leq l.$ Similarly, (2.2) can be reduced to

[TABLE]

where $l\geq 1,0\leq|m|\leq l.$ From (3.7)-(3.9) and (3.17), we can derive that

[TABLE]

Similarly, we can derive that

[TABLE]

From Lemma 3.1, the boundary condition (2.3) can be reduced to

[TABLE]

From Lemma 3.2 and (3.7)-(3.9) we can derive that

[TABLE]

Similarly, we can derive that

[TABLE]

From Lemma 3.1, the boundary condition (2.4) can be reduced to

[TABLE]

Note that the modes $t_{l}^{m}$ (coefficients of $\mathbf{T}_{l}^{m}$) are decoupled from the modes $v_{l}^{m}$ and $w_{l}^{m}$.

In summary, we only have to solve the following sequence ($l\geq 1$ and $|m|\leq l$) of one-dimensional eigenvalue problems, i.e., the so-called TE mode:

[TABLE]

and TM mode:

[TABLE]

4 Weak formulation and error estimation of the TE mode

For brief, we shall only give the error analysis in detail for the TE mode. For the TM mode, it is a coupled system with four unknown functions which is a challenging problem for numerical calculation. We propose an efficient numerical algorithm by using Legendre approximation based on vector basis functions.

4.1 Weak formulation and discrete formulation

By simplification, the TE mode can be rewritten as follows:

[TABLE]

Next, we formulate (4.1) - (4.4) as a equivalent fourth-order eigenvalue problem. Let $\tilde{u}_{l}=t_{l}^{m}(r)-\bar{t}_{l}^{m}(r)$, $\tilde{L}_{l}\tilde{u}_{l}=\frac{1}{r^{2}}\partial_{r}(r^{2}\partial_{r}\tilde{u}_{l})-\frac{l(l+1)}{r^{2}}\tilde{u}_{l}$. Subtracting (4.2) from (4.1), we obtain

[TABLE]

Divding $n-1$ and applying $\tilde{L}_{l}+k_{l}^{2}n$ to both sides of above equation, we obtain

[TABLE]

Then the (4.1)-(4.4) is equivalent to the following fourth order eigenvalue problem:

[TABLE]

Let $r=\frac{t+1}{2}R$, $\tilde{n}(t)=n(\frac{t+1}{2}R)$, $u_{l}=\tilde{u}_{l}(\frac{t+1}{2}R),$ $L_{l}u_{l}=\frac{1}{(t+1)^{2}}\partial_{t}((t+1)^{2}\partial_{t}u_{l})-\frac{l(l+1)}{(t+1)^{2}}u_{l}$. Then (4.7)-(4.8) can be restated as follows:

[TABLE]

We now define the usual weighted Sobolev space:

[TABLE]

equipped with the inner product and norm

[TABLE]

where $I=(-1,1)$ and $\omega=1+t$ is a weight function. Next, we introduce the following non-uniformly weighted Sobolev space:

[TABLE]

equipped with the corresponding inner product and norm

[TABLE]

Then the weak form of (4.9) -(4.10) is: Find $(k^{2}_{l},u_{l})\in\mathbb{C}\times H_{0,\omega,l}^{2}(I)$, such that

[TABLE]

Note the condition $u_{l}(-1)=0$ in $H_{0,\omega,l}^{2}(I)$ is a essential polar condition which should be imposed for the well-posedness of the weak form (4.17) (and the same type of pole condition for $v$).

We now introduce an associated generalized eigenvalue problems. First we define

[TABLE]

for $\tilde{n}>1$, and

[TABLE]

for $\tilde{n}<1$,

[TABLE]

where $\tau_{l}=k_{l}^{2}$. Then (4.17) can be written as either

[TABLE]

or

[TABLE]

The associated generalized eigenvalue problem is: Find $(\lambda(\tau_{l}),u_{l})\in\mathbb{C}\times H_{0,\omega,l}^{2}(I)$, such that

[TABLE]

or

[TABLE]

where $\lambda(\tau_{l})$ is a continuous function of $\tau_{l}$. From (4.21)-(4.22), we know that a transmission eigenvalue is the root of $f(\tau_{l})=\lambda(\tau_{l})-\tau_{l}$.

Let $P_{N}$ be the space of polynomials of degree less than or equal to $N$, and setting $X_{N}=P_{N}\cap H_{0,\omega,l}^{2}$. Then the discrete formulation of (4.23) is: Find $(\lambda_{N}(\tau_{l}),u_{lN})\in\mathbb{C}\times X_{N}$, such that

[TABLE]

The discrete formulation of (4.24) is: Find $(\lambda_{N}(\tau_{l}),u_{lN})\in\mathbb{C}\times X_{N}$, such that

[TABLE]

4.2 Error estimation of approximate eigenvalues

Hereafter, we shall use the expression $a\lesssim b$ to mean that there exists a positive constant $C$ such that $a\leq Cb$.

     

Lemma 4.1 It holds:

[TABLE]

Proof. Since

[TABLE]

then

[TABLE]

From Hardy inequality (cf. B 8.6 in [17]) we have

[TABLE]

Thus, we can obtain

[TABLE]

Theorem 1

Let $\tilde{n}\in L^{\infty}(I)$ satisfy

[TABLE]

for some $\alpha>0$ and $\beta>0$ positive constants. Then $A_{\tau_{l}}$ or $\tilde{A}_{\tau_{l}}$ is a continuous and coercive sesquilinear form on $H_{0,\omega,l}^{2}(I)\times H_{0,\omega,l}^{2}(I)$, i.e.,

[TABLE]

where $\tilde{n}_{*}=\inf_{I}(\tilde{n})$ and $\tilde{n}^{*}=\sup_{I}(\tilde{n})$.

Proof. We shall only give the proof in the case of $1+\alpha\leq\tilde{n}_{*}\leq\tilde{n}\leq\tilde{n}^{*}<\infty$. It can be similarly derived for the case of $~{}~{}0<\tilde{n}_{*}\leq\tilde{n}\leq\tilde{n}^{*}<1-\beta$.

Since

[TABLE]

for $\frac{1}{\tilde{n}^{*}-1}<\varepsilon<\frac{1}{\tilde{n}^{*}-1}+1$, then from Lemma 4.1 we can obtain (4.27). Due to

[TABLE]

then from Lemma 4.1 we can obtain (4.28).           $\Box$

Similar to the proof of Theorem 1, we can derive the following Theorem:

Theorem 2

$B(u,v)$* is a continuous and coercive bilinear form on $H_{0,\omega,l}^{1}(I)\times H_{0,\omega,l}^{1}(I)$, i.e.,*

[TABLE]

To give the error analysis, we will use extensively the minimax principle.

Lemma 4.2 Let $\lambda^{m}(\tau_{l})$ denote the $m$-th eigenvalues of (4.23) and $V_{m}$ be any $m$-dimensional subspace of $H_{0,\omega,l}^{2}(I)$. Then, for $\lambda^{1}(\tau_{l})\leq\lambda^{2}(\tau_{l})\leq\cdots\leq\lambda^{m}(\tau_{l})\leq\cdots$, there holds

[TABLE]

Proof. See Theorem 3.1 in [18].          $\Box$

Lemma 4.3 Let $\lambda^{m}(\tau_{l})$ denote the $m$-th eigenvalues of (4.23) and be arranged in an ascending order, and define

[TABLE]

where $u_{l}^{i}$ is the eigenfunction corresponding to the eigenvalue $\lambda^{i}(\tau_{l})$. Then we have

[TABLE]

Proof. See Lemma 3.2 in [18].           $\Box$

It is true that the minimax principle is also valid for the discrete formulation (4.25) (see [18]).

Lemma 4.4 Let $\lambda_{N}^{m}(\tau_{l})$ denote the $m$-th eigenvalues of (4.25), and $V_{m}$ be any $m$-dimensional subspace of $X_{N}$. Then, for $\lambda_{N}^{1}(\tau_{l})\leq\lambda_{N}^{2}(\tau_{l})\leq\cdots\leq\lambda_{N}^{m}(\tau_{l})\leq\cdots$, there holds

[TABLE]

Let $\prod_{N}^{2,l}:H_{0,\omega,l}^{2}(I)\rightarrow X_{N}$ be an orthogonal projection, defined by

[TABLE]

Theorem 3

Let $\lambda_{N}^{m}(\tau_{l})$ be obtained by solving (4.25) and be an approximation of $\lambda^{m}(\tau_{l})$, an eigenvalue of (4.23). Then, we have

[TABLE]

Proof. According to the coerciveness of $A_{\tau_{l}}(u,v)$ and $B(u,v)$ we can get $\lambda^{m}(\tau_{l})>0$. From $X_{N}\subset H_{0,\omega,l}^{2}(I)$, together with (4.29) and (4.32) we can obtain $\lambda^{m}(\tau_{l})\leq\lambda_{N}^{m}(\tau_{l})$. Let $\Pi_{N}^{2,l}E_{1,m}$ denote the space spanned by $\Pi_{N}^{2,l}u_{l}^{1},\Pi_{N}^{2,l}u_{l}^{2},\cdots,\Pi_{N}^{2,l}u_{l}^{m}$. It is obvious that $\Pi_{N}^{2,l}E_{1,m}$ is a $m$-dimensional subspace of $X_{N}$. From the minimax principle, we have

[TABLE]

Since $A_{\tau_{l}}(v,v)=A_{\tau_{l}}(\Pi_{N}^{2,l}v,\Pi_{N}^{2,l}v)+2A_{\tau_{l}}(v-\Pi_{N}^{2,l}v,\Pi_{N}^{2,l}v)+A_{\tau_{l}}(v-\Pi_{N}^{2,l}v,v-\Pi_{N}^{2,l}v),$ from $A_{\tau_{l}}(v-\Pi_{N}^{2,l}v,\Pi_{N}^{2,l}v)=0$ and the non-negativity of $A_{\tau_{l}}(v-\Pi_{N}^{2,l}v,v-\Pi_{N}^{2,l}v)$, we have

[TABLE]

Thus, we have

[TABLE]

The proof of Theorem 3 is complete.           $\Box$

We denote the Jacobi weight function of index $(\alpha,\beta)$ by $\omega^{\alpha,\beta}(x)=(1-t)^{\alpha}(1+t)^{\beta}$ and introduce the following non-uniformly weighted Sobole spaces:

[TABLE]

equipped with the inner product and norm

[TABLE]

and

[TABLE]

equipped with the inner product and norm

[TABLE]

Define the orthogonal projection $\pi_{N,\omega^{-2,-2}}:L_{\omega^{-2,-2}}^{2}(I)\rightarrow Q_{N}^{-2,-2}$ by

[TABLE]

where $Q_{N}^{-2,-2}=\{\phi\in P_{N}:\phi(\pm 1)=\phi^{\prime}(\pm 1)=0\}.$ From the Theorem 1.8.2 in [19] we have the following Lemma:

Lemma 4.5. For any $u\in H_{\omega^{-2,-2},*}^{s}(I)$, the following inequality holds:

[TABLE]

Theorem 4

There exists an operator $\pi_{N}^{2,l}:H_{0,\omega,l}^{2}(I)\rightarrow P_{N}^{0,l}$ such that $\pi_{N}^{2,l}u(\pm 1)=u(\pm 1)=0$, and $\partial_{t}\pi_{N}^{2,l}u(-1)=\partial_{t}u(-1)$, $\partial_{t}\pi_{N}^{2,l}u(1)=\partial_{t}u(1)=0$ and for $u\in H_{\omega^{-2,-2},l}^{s}(I)$ with $s\geq 2$, there holds

[TABLE]

where $P_{N}^{0,l}=\{\phi\in P_{N}:\phi(\pm 1)=\phi^{\prime}(1)=0\}$.

Proof. Let $u_{*}(t)=\frac{1}{4}(t+1)(1-t)^{2}\partial_{t}u(-1)$ for $\forall u\in H_{0,\omega,l}^{2}(I)$. By construction, we have $\partial_{t}^{k}(u-u_{*})(\pm 1)=0,(k=0,1).$ If $u\in H_{\omega^{-2,-2},l}^{s}(I)$ , then we have $u-u_{*}\in H_{\omega^{-2,-2},*}^{s}(I)$. In fact, from Hardy inequality (cf. B 8.8 in [17]) we have

[TABLE]

Thus, we can derive that

[TABLE]

Since

[TABLE]

then we have

[TABLE]

Similarly, we can derive that $\int_{I}\omega^{1,1}(\partial_{t}^{3}u_{*})^{2}dt=3(\partial_{t}u(-1))^{2}\lesssim\int_{I}(\partial^{2}_{t}u)^{2}dt.$ Thus we have

[TABLE]

For $k>3$, we have

[TABLE]

Thus, $u-u_{*}\in H_{\omega^{-2,-2},*}^{s}(I)$ and we can define

[TABLE]

Then by Lemma 4.5 we can obtain

[TABLE]

Together with (4.35), (4.36) and (4.37) we can get desired results. $\Box$

Theorem 5

*Let $\lambda^{m}_{N}(\tau_{l})$ is the $m$-th approximate eigenvalue of $\lambda^{m}(\tau_{l})$.

If $\{u_{l}^{i}\}_{i=1}^{m}\subset H_{\omega^{-2,-2},l}^{s}(I)$ with $s\geq 2$, then we have*

[TABLE]

where $C(m)$ is a constant independent of $N$.

Proof. For any $v\in E_{1,m}$ it can be represented by $v=\sum_{i=1}^{m}\mu_{i}u_{l}^{i}$, we then have

[TABLE]

From Cauchy-Schwarz inequality we have

[TABLE]

From Hardy inequality (cf. B 8.6 in [17]) we have

[TABLE]

Then from Poincar$\acute{e}$ inequality we can obtain

[TABLE]

From the property of orthogonal project (4.34) and continuity of $A_{\tau_{l}}(u,v)$ in Theorem 1 we can derive that

[TABLE]

Since

[TABLE]

then from Theorem 3 and Theorem 4 we can get desired results. $\Box$

5 Efficient implementation of the algorithm

We describe in this section how to solve the TE mode (3.35)-(3.38) and TM mode (3.39)-(3.44) efficiently.

5.1 Efficient implementation of the algorithm for TE mode

We start by constructing a set of basis functions for $X_{N}$. Let

[TABLE]

where $d_{i}=\frac{1}{\sqrt{2(2i+3)^{2}(2i+5)}}$ and $L_{i}$ is the Legendre polynomial of degree $i$.

It is clear that

[TABLE]

where $\phi_{N-3}(t)=\frac{1}{4}(t+1)(t-1)^{2}$. Setting

[TABLE]

We shall look for

[TABLE]

Now, plugging the expression of (5.1) in (4.25), and taking $v_{N}$ through all the basis functions in $X_{N}$, we will arrive at the following linear eigenvalue system:

[TABLE]

where

[TABLE]

5.2 Efficient implementation of the algorithm for TM mode

Let $r=\frac{t+1}{2}R$, $W(t)=w_{l}^{m}(\frac{t+1}{2}R)$,$V(t)=v_{l}^{m}(\frac{t+1}{2}R)$,$\tilde{n}(t)=n(\frac{t+1}{2}R)$, $w(t)=\bar{w}_{l}^{m}(\frac{t+1}{2}R)$,$v(t)=\bar{v}_{l}^{m}(\frac{t+1}{2}R)$, $D_{l}^{\pm}=\partial_{t}\pm\frac{l}{l+1}$. Then (3.39)-(3.44) can be restated as follows:

[TABLE]

Let us define

[TABLE]

Then the weak form of (5.3)-(5.8) is: Find $(k_{l}^{2},(W(t),V(t),w(t),v(t))\neq\mathbf{0})\in\mathbb{C}\times\mathbf{H}(I),$ such that for all $(h_{1}(t),h_{2}(t),h_{3}(t),h_{4}(t))\in\mathbf{H}(I)$,

[TABLE]

Let $\mathbf{H}_{N}=\mathbf{H}(I)\cap P_{N}^{4}$, then the discrete form of (5.9) is: Find $(k_{Nl}^{2},(W_{N},V_{N},w_{N},v_{N})\neq\mathbf{0})\in\mathbb{C}\times\mathbf{H}_{N},$ such that for all $(h_{1N}(t),h_{2N}(t),h_{3N}(t),h_{4N}(t))\in\mathbf{H}_{N}$,

[TABLE]

We now construct a set of basis functions for $\mathbf{H}_{N}$. Let

[TABLE]

Setting

[TABLE]

where $j=0,1,\cdots,N-1$. It is clear that

[TABLE]

We shall look for

[TABLE]

Now, plugging the expression of (5.11) in (5.10), and taking $(h_{1N}(t),h_{2N}(t),h_{3N}(t),h_{4N}(t))$ through all the basis functions in $\mathbf{H}_{N}$, we will arrive at the following linear eigenvalue system:

[TABLE]

where $\mathcal{U}=(\alpha_{0},\alpha_{1},\cdots,\alpha_{4N+2})^{T},$ $\mathcal{A}$ and $\mathcal{B}$ are the corresponding stiff matrix and mass matrix, respectively .

Note that the stiff matrix $A$ (resp. $\mathcal{A}$) and the mass matrix $B$ (resp. $\mathcal{B}$) in (5.2) (resp. (5.12)) are all sparse for the constant $n$. Thus either direct or iterative eigen solvers can be efficiently applied in parallel. For the case of variable $n$, since one is mostly interested in a few smallest transmission eigenvalues, it is most efficient to solve (5.2) (resp. (5.12)) by using a shifted inverse power method (cf., for instance, [20]) which requires solving, repeatedly for different righthand side $f$ (resp. $\tilde{f}$),

[TABLE]

where $\lambda_{a}(\tau_{l})$ (resp. $k_{al}$) is some approximate value for the transmission eigenvalue $\lambda_{N}(\tau_{l})$ (resp. $k_{Nl}$). The above system can be efficiently solved by the Schur-complement approach, we refer to [21] for a detailed description on a related problem. In summary, the approximate transmission eigenvalue problem (5.2) (resp. (5.12)) can be solved very efficiently.

6 Numerical experiments

We now perform a sequence of numerical tests to study the convergence behavior and show the effectiveness of our algorithm. We operate our programs in MATLAB 2015b.

6.1 Homogeneous medium $n$

We take $R=1$, $n=16$, and $l=1,2,3$ in our examples. Numerical results for the first four eigenvalues with different $l$ and $N$ are listed in Tables 6.1-6.3 for the TE mode and in Tables 6.4-6.6 for the TM mode, respectively.

Table 6.1 The first four eigenvalues of the TE mode for $l=1$ and different $N$ in unit ball

N $1st$ $2nd$ $3rd$ $4th$

10 1.460855902327352 2.309268980991891 3.142098536003481 4.004131427018431

15 1.460842273223355 2.309270674683650 3.141592652865679 4.028313168694923

20 1.460855902076009 2.309270674683547 3.141592653589792 4.028312376370235

25 1.460855902076010 2.309270674683545 3.141592653589794 4.028312376370704

30 1.460855902076010 2.309270674683548 3.141592653589792 4.028312376370695

**Table 6.2 The first four eigenvalues of the TE mode for $l=2$ and different $N$ in unit ball **

N $1st$ $2nd$ $3rd$ $4th$

10 1.764042417090797 2.631690641913292 3.465152456664587 4.304027448259460

15 1.764042422029338 2.631678257808169 3.465236228216971 4.293583051636833

20 1.764042422029338 2.631678257809425 3.465236224179554 4.293582919867683

25 1.764042422029339 2.631678257809425 3.465236224179551 4.293582919866948

30 1.764042422029338 2.631678257809422 3.465236224179552 4.293582919866945

Table 6.3 The first four eigenvalues of the TE mode for $l=3$ and different $N$ in unit ball

N $1st$ $2nd$ $3rd$ $4th$

10 2.061050417316723 2.949614759058087 3.789381768848592 4.729312472707909

15 2.061050433015994 2.949488215613814 3.792296568087934 4.619875856328957

20 2.061050433015994 2.949488215659479 3.792296458205384 4.619887058309071

25 2.061050433015994 2.949488215659483 3.792296458205412 4.619887058253892

30 2.061050433015993 2.949488215659481 3.792296458205412 4.619887058253892

Table 6.4 The first four eigenvalues of the TM mode for $l=1$ and different $N$ in unit ball

N $1st$ $2nd$ $3rd$ $4th$

10 1.165407223825574 2.045867252670490 3.158017632244895 3.442289242670909

15 1.165407223827102 2.045867782103358 3.418097647742193 4.292461561664239

20 1.165407223827106 2.045867782103363 3.418097651533288 4.292488875027842

25 1.165407223827098 2.045867782103357 3.418097651533263 4.292488875029357

30 1.165407223827108 2.045867782103361 3.418097651533326 4.292488875029386

**Table 6.5 The first four eigenvalues of the TM mode for $l=2$ and different $N$ in unit ball **

N $1st$ $2nd$ $3rd$ $4th$

10 1.475116524235235 2.340653677153939 3.231276889331592 3.354682219133095

15 1.475116524493845 2.340657592735246 3.233313705482311 4.557043549625480

20 1.475116524493846 2.340657592735365 3.233313708702767 4.557097304644416

25 1.475116524493844 2.340657592735367 3.233313708702756 4.557097304725269

30 1.475116524493843 2.340657592735366 3.233313708702761 4.557097304725263

Table 6.6 The first four eigenvalues of the TM mode for $l=3$ and different $N$ in unit ball

N $1st$ $2nd$ $3rd$ $4th$

10 1.777410985980160 2.656258796160595 3.455344578108596 3.582237584322505

15 1.777410996101287 2.656264636190942 3.512014047361030 4.421816274456907

20 1.777410996101286 2.656264636197187 3.512014051598614 4.421843661628902

25 1.777410996101286 2.656264636197187 3.512014051598621 4.421843661635361

30 1.777410996101284 2.656264636197187 3.512014051598614 4.421843661635358

We see from Tables 6.1-6.6 that numerical eigenvalues achieve at least fourteen-digit accuracy with $N\geq 25$. As a comparison, we list in Table 6.7 the results in [10] computed by the finite element method with mesh size $h\thickapprox 0.2$. We observe that numerical eigenvalues listed in Table 6.7 have only two- to three-digit accuracy, while our method captured at least fourteen-digit accuracy with much less computational efforts.

Note that the multiplicity of the eigenvalue can be predicted by $|m|\leq l$. For $l=1$, we have $m=0,\pm 1$ and hence the first eigenvalue has multiplicity 3; for $l=2$, we have $m=0,\pm 1,\pm 2$ and therefore the second eigenvalue has multiplicity 5. In general, the $j$-th eigenvalue has multiplicity $2j+1$. The first three numbers (read in horizontal) in Table 6.7 are for the TM mode with $l=1$, which compared with the first column in Table 6.4; the second three numbers are for the TE mode with $l=1$, which compared with the first column in Table 6.1; the next five numbers are for the TM mode with $l=2$, which compared with the first column in Table 6.5; the following five numbers are for the TM mode with $l=2$, which compared with the first column in Table 6.2; the last seven numbers in Table 6.7 are for the TE mode with $l=3$, which compared with the first column in Table 6.6.

Table 6.7 Computed Maxwell’s transmission eigenvalues for the unit ball

1.1741 1.1717 1.1721 1.4665 1.4667 1.4671 1.4824 1.4828

1.4828 1.4830 1.4836 1.7690 1.7690 1.7698 1.7700 1.7705

1.7857 1.7859 1.7862 1.7865 1.7867 1.7868 1.7872

In order to further demonstrate accuracy and efficiency of our algorithm, we use numerical solutions with $N=60$ as reference solutions and plot errors of approximated eigenvalues with different $N$ in Figures 1-6. We observe exponential decay in all cases.

6.2 Inhomogeneous medium $n$

We take $R=1$, $n=8+4r^{2}$, and $l=1$ in our examples. Numerical results of the first four eigenvalues with different $N$ are listed in Table 6.8 for the TE mode and in Table 6.9 for the TM mode, respectively.

Table 6.8 The first four eigenvalues of the TE mode for $l=1$ and different $N$ in unit ball

N $1st$ $2nd$ $3rd$ $4th$

10 1.924760241224309 3.066320509754762 4.953581565194202 6.219217568046522

15 1.924760240239596 3.066318451356159 4.944962746965171 6.162042706007884

20 1.924760240239595 3.066318451356097 4.944962719618220 6.162013704025296

25 1.924760240239596 3.066318451356097 4.944962719618172 6.162013703949516

30 1.924760240239597 3.066318451356096 4.944962719618174 6.162013703949522

Table 6.9 The first four eigenvalues of the TM mode for $l=1$ and different $N$ in unit ball

N $1st$ $2nd$ $3rd$ $4th$

10 1.546722576754737 3.418052693123285 4.366864963075772 4.640915773722154

15 1.546722576768443 3.418109299464758 4.616102608978945 6.423176009030875

20 1.546722576768448 3.418109299467635 4.616102624493460 6.425723290604534

25 1.546722576768440 3.418109299467642 4.616102624493454 6.425723292013815

30 1.546722576768443 3.418109299467622 4.616102624493481 6.425723292013920

We see from Tables 6.8-6.9 that numerical eigenvalues achieve at least thirteen-digit accuracy with $N\geq 25$. Here again, we use numerical solutions with $N=60$ as reference solutions and plot the errors in Figure 7-8.

7 Conclusions

In summary, the method developed in this paper is a first but important step towards a robust and accurate algorithm for more general transmission eigenvalue problems which will be the subject of our future work.