### scattering problem in pseudo-chiral media and a practical reconstruction method

Tiexiang Li

School of Mathematics, Southeast University, Nanjing 211189, People’s Republic of China.

E-mail: txli@seu.edu.cn

Tsung-Ming Huang

Department of Mathematics, National Taiwan Normal University, Taipei 116, Taiwan.

E-mail: min@ntnu.edu.tw

Wen-Wei Lin

Department of Applied Mathematics, National Chiao Tung University, Hsinchu 300, Taiwan.

E-mail: wwlin@math.nctu.edu.tw

Jenn-Nan Wang

Institute of Applied Mathematical Sciences, NCTS, National Taiwan University, Taipei, 106, Taiwan.

E-mail: jnwang@math.ntu.edu.tw

Abstract. In this paper, we consider the two-dimensional Maxwell’s equations with
the TM mode in pseudo-chiral media. The system can be reduced to the acoustic
equation with a negative index of refraction. We first study the transmission eigenvalue
problem (TEP) for this equation. By the continuous finite element method, we
discretize the reduced equation and transform the study of TEP to a quadratic
eigenvalue problem by deflating all nonphysical zeros. We then estimate half of the
eigenvalues are negative with order of O(1) and the other half of eigenvalues are
positive with order of O(10^{2}). In the second part of the paper, we present a practical
numerical method to reconstruct the support of the inhomogeneity by the near-field
measurements, i.e., Cauchy data. Based on the linear sampling method, we propose
the truncated singular value decomposition to solve the ill-posed near-field integral
equation, at one wave number which is not a transmission eigenvalue. By carefully
chosen an indicator function, this method produce different jumps for the sampling
points inside and outside the support. Numerical results show that our method is able
to reconstruct the support reliably.

Keywords: Two-dimensional transmission eigenvalue problem, pseudo-chiral model,

transverse magnetic mode, linear sampling method, singular value decomposition

1. Introduction

The transmission eigenvalue problem (TEP) has attracted a lot of attention recently in the study of direct/inverse scattering problems in inhomogeneous media [3, 4, 5, 6, 9, 11,12, 20, 25]. The existence of transmission eigenvalues is intimately connected to the

”bijectivity” of the far-field operator, which is crucial in some reconstruction methods such as the linear sampling method (LSM) [1, 8] and the factorization method [21].

Conversely, the transmission eigenvalues (also eigenvalues) carry the information of the
scatterer and can be estimated by the far-field data [2] or the near-field data (Cauchy
data) [27]. This observation leads to the development of some reconstruction methods
using the transmission eigenvalues or eigenvalues, see for example, [10, 29]. In [29],
an eigenvalue method using multiple frequency near-field data (EM^{2}F) was proposed
to detect Dirichlet or transmission eigenvalues and to reconstruct the support of the
scatterer.

In this paper, we would like to propose a reliable numerical method to investigate the distribution of transmission eigenvalues for the 2d acoustic equation with an inhomogeneous index of refraction. The model is derived from the Maxwell’s equations with the TM mode in pseudo-chiral media. It turns out the index of refraction of the reduced acoustic equation decreases to negative infinity when we increase the charity parameter. Precisely, we consider the TEP

∆u + λε(x)u = 0, in D, (1a)

∆v + λv = 0, in D, (1b)

u − v = 0, on ∂D, (1c)

∂u

∂ν − ∂v

∂ν = 0, on ∂D (1d)

for the scattering of acoustic wave on a bounded and simply connected inhomogeneous
domain D ⊆ R^{2}, where ν is the outer normal to the smooth boundary ∂D, u, v ∈ L^{2}(D)
with u−v ∈ H_{0}^{2}(D) = {w ∈ H^{2}|w = 0,^{∂w}_{∂ν} = 0}, and ε(x) is the index of refraction. Any
λ ∈ C such that (1) has nontrivial solutions u and v is called a transmission eigenvalue
and u, v are called the corresponding transmission eigenfunctions for D .

Equation (1a) can be considered as a reduced Maxwell’s equations with transverse magnetic (TM) mode. For the standard Maxwell’s equations, ε(x) is the electric permittivity corresponding to the product of the dielectric constant and the free-space dielectric constant. In practice, ε can not be taken as large as we wish. In [23]

where the Maxwell’s equations for Tellegen media is studied, ε(x) is the sum of the electric permittivity and the square of Tellegen parameter. By choosing large Tellegen parameter, we can enlarge the parameter ε(x) as we want. On the contrary, here ε(x) is the sum of the electric permittivity minus the square of pseudo-chiral media parameters. Therefore, by selecting the pseudo-chiral parameters sufficiently large, ε(x) becomes negative (see Section 2below).

In recent years, several papers have been devoted to developing efficient algorithms

for computing transmission eigenvalues of 2d/3d TEP [3, 11, 14, 16, 17,20, 21, 22, 23, 24, 26, 27]. Three finite element methods (FEMs) and a coupled boundary element method were developed to solve the 2d/3d TEP in [11, 14, 17] (see the book [30] for more details). Two iterative methods and the corresponding convergence analysis were given in [28]. In [17], a mixed finite element method for 2d TEP was proposed, which leads to a non-Hermitian quadratic eigenvalue problem (QEP) that was solved by an adaptive Arnoldi method. Furthermore, a multilevel correction method was used to reduce the solution of TEP into some linear boundary value problems which could be solved by the multigrid method [18]. For results related to our current work, in [23, 24], the TEP for general inhomogeneous media are discretized into QEP with symmetric coefficient matrices. For such QEP, a secant-type iteration for computing several smallest positive transmission eigenvalues accurately was proposed in [24], and a quadratic Jacobi-Davidson method with nonequivalent deflation technique for computing a large number positive transmission eigenvalues was developed in [23].

Remind that the index of refraction ε(x) increases to infinity as we increase the Tellegen parameter. Numerical simulations demonstrate that the positive transmission eigenvalues are densely distributed in an interval near the origin. Similar to the method in [23], the TEP (1) is discretized into a QEP and a quadratic Jacobi-Davidson method with nonequivalence deflation is applied to compute a large number of positive eigenvalues. It turns out in the case here there exists an eigenvalue-free interval near the origin. The existence of the eigenvalue-free interval motivates us to study the reconstruction of the support of ε(x) from the near-field data in the spirit of LSM.

Corresponding to the TEP(1), the scattering problem is described by

∆u + k^{2}ε(x)u = 0, in R^{2}\ {x_{0}}, (2a)

u = u^{i}+ u^{s}, (2b)

r→∞lim

√r ∂u^{s}

∂r − iku^{s}

= 0, (2c)

where D := supp(ε(x) − 1) and u^{i} is the incident field due to a point source at x_{0}, i.e.,
u^{i}(x, x_{0}) := Φ(x, x_{0}), x_{0} ∈ C, (3)
where Φ(x, x_{0}) = _{4}^{i}H_{0}^{(0)}(k|x − x_{0}|), the fundamental solution of the Helmholtz equation
in R^{2}.

Our reconstruction method is based on the set up in [29]. We assume that the target
D is inside some domain Ω which itself is inside a curve C (see Figure1). Suppose that
u^{s}is measured on Γ = ∂Ω for all point sources x on C. We define the near-field operator
N : L^{2}(Γ) → L^{2}(C) for v ∈ L^{2}(Γ) by

(N v)(x) = Z

Γ

u^{s}(y, x)v(y)ds(y), x ∈ C. (4)
The LSM with near-field data relies on solving the ill-posed integral equations

(N v)(x) = Φ(x, z) for all x ∈ C, (5)

Figure 1. The target D is inside some domain Ω (Γ := ∂Ω) which itself is surrounded
by a curve C. The scattered field u^{s} is due to the scattering of the incident field u^{i}
having a point source at x_{0}∈ C.

where z ∈ T is a sampling point and T is a sampling domain inside Ω containing the target D (see Figure 1). In general, the above ill-posed equation do not have a solution.

The following theorem serves as the backbone of the LSM.

Theorem 1.1. [8, 7, 29] Assume that λ = k^{2} is not a transmission eigenvalue of (1).

Let N be the near-field operator defined by (4).

• If z ∈ D, then there exists a convergent sequence v_{n} in L^{2}(D), such that

n→∞lim N v_{n}= Φ(·, z).

• If z ∈ Ω \ D, then for every sequence vn satisfying

n→∞lim N v_{n}= Φ(·, z),
we have

n→∞lim ||v_{n}||_{L}^{2}_{(D)}= ∞.

Our reconstruction method is a direct application of Theorem1.1. By this theorem,
we can see that if k^{2} is not a transmission eigenvalue then we can determine whether
z ∈ D or z ∈ Ω \ D from the behaviors of the solutions to (4). Since there exists
an eigenvalue-free interval near the origin in TEP (1), we choose a k near the origin
so that k^{2} is not a transmission eigenvalue, we then discretize the integral equation
(5). The strategy here is to choose m points x1, · · · , xm ∈ C forming a m-vector
b = [Φ(x_{1}, z), · · · , Φ(x_{m}, z)]^{>} ∈ C^{m} and n unknowns [v(y_{1}), · · · , v(y_{n})]^{>} := v satisfying
the linear system

Av = b. (6)

The idea is to take m ≥ n, i.e., (6) is overdetermined. To determine whether z ∈ D or z ∈ Ω \ D, we look at the distance between the vector b and the space A = span(A), the subspace spanned by the columns of A, namely,

dist(b, A) ≡ min

v∈C^{n}||Av − b||_{2}.

If dist(b, A) > 0, then, intuitively, any ”solution” v satisfying (6) must contain some sufficiently large components. Mimicing the second case of Theorem1.1, we thus assign z ∈ Ω \ D. On the other hand, if dist(b, A) = 0, i.e., b ∈ A, then the norm of v is clearly finite. We thus say that z ∈ D in view of the first case of Theorem 1.1. An easy application of truncated SVD will quickly determine the value of dist(b, A). Our criterion is very simple and easily to be implemented.

This paper is organized as follows. In Section 2, we demonstrate the derivation of 2d TEP (1) from the Maxwell’s equations with TM mode in pseudo-chiral media.

A corresponding discretized QEP and its spectral analysis are given in Section 3. In Section 4, a practical numerical method based on the LSM and truncated SVD for the reconstruction of the the target D is developed. Related numerical results are presented in Section 5. A concluding remark is given in Section 6. Finally, in Appendix A, we present a framework of solving the direct scattering problem using FEM. The purpose of solving the direct problem is to obtain synthetic data for numerical simulations.

2. Maxwell’s equations with the TM mode in pseudo-chiral media

Physically, the governing equations for the propagation of electromagnetic wave in bi- isotropic materials with complex media is modeled by the 3d source-free frequency domain Maxwell’s equations

∇ × E = ik (µH + ζE) , (7a)

∇ × H = −ik (ε_{r}E + ξH) , (7b)

where E and H are the electronic field and magnetic field respectively, k is the frequency,
ε_{r} and µ are electric permittivity and magnetic permeability respectively, ξ and ζ are 3-
by-3 magnetoelectric parameter matrices in various forms for describing different types
of complex media.

We consider E and H in (7) in the transversal magnetic (TM) mode as

E = [0, 0, E3(x)]^{>}, H = [H1(x), H2(x), 0]^{>} (8)
with x = (x_{1}, x_{2})^{>} ∈ R^{2}, ζ and ξ are pseudo-chiral media of the forms

ζ =

0 0 ζ_{1}

0 0 ζ_{2}

−ζ1 −ζ2 0

, ξ =

0 0 ξ_{1}

0 0 ξ_{2}

−ξ1 −ξ2 0

(9)

with ξ_{1} = ζ_{1} = iγ_{1}, ξ_{2} = ζ_{2} = iγ_{2}, and γ_{1}, γ_{2} being real numbers. We assume in this
paper µ = 1. Then the equation (7a) can be simplified to

∂x2E3

−∂_{x}_{1}E_{3}
0

= ik

H1

H_{2}
0

+

ζ1E3

ζ_{2}E_{3}
0

. (10)

Substituting (10) into (7b) yields

(ik)^{−1}

0 0

−

∂^{2}

∂x^{2}_{1} + _{∂x}^{∂}^{2}2
2

E_{3}

−

0 0

∂

∂x1 (ζ_{2}E_{3}) − _{∂x}^{∂}

2 (ζ_{1}E_{3})

= − ik

ε_{r}

0
0
E_{3}

−

0
0
ξ_{1}H_{1} + ξ_{2}H_{2}

, which implies that

− ∆E_{3}

= k^{2}

ε_{r}E_{3}− ξ_{1}

(ik)^{−1} ∂

∂x_{2}E_{3}− ζ_{1}E_{3}

+ ξ_{2}

(ik)^{−1} ∂

∂x_{1}E_{3}+ ζ_{2}E_{3}

+ ik

∂

∂x_{1} (ζ_{2}E_{3}) − ∂

∂x_{2}(ζ_{1}E_{3})

= k^{2}(ε_{r}+ ξ_{1}ζ_{1}+ ξ_{2}ζ_{2}) E_{3}+ ik

∂

∂x_{1} (ζ_{2}E_{3}) − ∂

∂x_{2} (ζ_{1}E_{3}) + ξ_{1} ∂

∂x_{2}E_{3}− ξ_{2} ∂

∂x_{1}E_{3}

=k^{2}(ε_{r}− γ_{1}^{2}− γ_{2}^{2})E_{3}

=k^{2}εE3

with ε = ε_{r}− (γ_{1}^{2}+ γ_{2}^{2}).

Let E_{0} = [0, 0, E_{0,3}(x)]^{>} and H_{0} = [H_{0,1}(x), H_{0,2}(x), 0]^{>} be, respectively, the
electronic and magnetic plane waves with TM mode in vacuum which are governed
by the free Maxwell’s equations

∇ × E0 = ikµ0H0, (11a)

∇ × H_{0} = −ikε_{0}E_{0} (11b)

with ε0 = µ0 = 1. We now suppose the Maxwell’s equations (7) and (11) are defined on a cylindrical set D × R satisfying the boundary conditions

E ×ν = Ee _{0}×ν,e (12a)

(H + ζE) ×ν = He _{0}×ν,e (12b)

whereν = [νe _{1}, ν_{2}, 0]^{>} = [ν, 0]^{>}is the outer unit normal to the smooth boundary ∂D×R.

It is easily seen from (12a) that the Dirchlet boundary condition

E3 = E0,3 on ∂D (13a)

holds. Multiplying (12b) by ik and using (11a) we get the Neumann boundary condition

∇ × E ×ν = ∇ × Ee _{0}×ν,e
i.e.,

∂E_{3}

∂ν = ∂E_{0,3}

∂ν on ∂D. (13b)

In view of the boundary conditions (13a), (13b), choosing u = E_{3} and v = E_{0,3}, we
arrive at the TEP (1) with λ = k^{2} and ε(x) = ε_{r}− (γ_{1}^{2}+ γ_{2}^{2}). The index of refraction
will become negative as long as either |γ1| or |γ2| are sufficiently large.

3. Spectral analysis of discretized TEP

In this section, we first briefly introduce the discretized TEP and give its spectral
analysis. Let {φ_{i}}^{n}_{i=1} and {ψ_{i}}^{m}_{i=1} be standard nodal bases for spaces of continuous
piecewise linear functions on D ⊂ R^{2}, that have vanishing DoF on ∂D and D,
respectively, where DoF denotes the degree of freedom. Applying the standard piecewise
linear FEM (we refer to [11] for details) to (1) with

u =

n

X

i=1

u_{i}φ_{i}+

m

X

i=1

w_{i}ψ_{i}, (14a)

v =

n

X

i=1

v_{i}φ_{i}+

m

X

i=1

w_{i}ψ_{i}. (14b)

Hereafter, we denote 0 and I as a zero matrix/submatrix and the identity matrix,
respectively, with appropriate dimensions if it is clear in the context. Now, setting
u = [u_{1}, · · · , u_{n}]^{>}, v = [v_{1}, · · · , v_{n}]^{>} and w = [w_{1}, · · · , w_{m}]^{>} appeared in (14), we can
derive a generalized eigenvalue problem (GEP)

L(λ)z = (A − λB)z = 0, (15)

in which λ = k^{2}, where

A =

K 0 E

0 K E

E^{>} −E^{>} 0

, B =

−M_{ε} 0 −F_{ε}

0 M_{1} F_{1}

−F_{ε}^{>} −F_{1}^{>} −G_{ε}− G_{1}

, z =

u v w

(16)

with K, E, M_{1}, M_{ε}, F_{1}, F_{ε}, G_{1} and G_{ε} being given in Table1. Here, A 0 means A
is symmetric positive definite.

We now define

M = M_{ε}+ M_{1} 0, G = G_{ε}+ G_{1} 0, F = F_{ε}+ F_{1}, (17a)
K = K − EGb ^{−1}F^{>}, Mc1 = M1− F1G^{−1}F^{>}, M = M − F Gc ^{−1}F^{>}, (17b)
and

S =h

K E

i

, T_{1} =h

M1 F1

i

, M =

"

M F

F^{>} G

#

. (17c)

stiffness matrix for interior meshes K = [(∇φi, ∇φj)] 0 ∈ R^{n×n}
stiffness matrix for

E = [(∇φ_{i}, ∇ψ_{j})] ∈ R^{n×m}
interior/boundary meshes

mass matrices for interior meshes M1 = [(φi, φj)] 0 ∈ R^{n×n}
M_{ε} = [−(εφ_{i}, φ_{j})] 0 ∈ R^{n×n}
mass matrices for interior/boundary meshes F_{1} = [(φ_{i}, ψ_{j})] ∈ R^{n×m}

F_{ε}= [−(εφ_{i}, ψ_{j})] ∈ R^{n×m}
mass matrices for boundary meshes G_{1} = [(ψ_{i}, ψ_{j})] 0 ∈ R^{m×m}

G_{ε} = [−(εψ_{i}, ψ_{j})] 0 ∈ R^{m×m}

Table 1. Stiffness and mass matrices with ε(x) < 0 for x ∈ ¯D.

Lemma 3.1. Let M and cM be defined in (17a), (17b). Then M, cM 0.

Proof. Because

"

M_{1} F_{1}
F_{1}^{>} G_{1}

#

0,

"

M_{ε} F_{ε}
F_{ε}^{>} G_{ε}

#

0, we see that

M =

"

M F

F^{>} G

#

=

"

M_{1} F_{1}
F_{1}^{>} G_{1}

# +

"

M_{ε} F_{ε}
F_{ε}^{>} G_{ε}

#

0.

On the other hand, observe that

0 ≺

"

I −F G^{−1}

0 I

# "

M F

F^{>} G

# "

I 0

−G^{−1}F^{>} I

#

=

"

Mc 0

0 G

# ,

which implies cM 0.

Theorem 3.2. For λ 6= 0, the GEP (15) can be reduced to the following QEP

Q(λ)p = (λ^{2}A_{2}− λA_{1}− A_{0})p = 0, (18)
where p = u − v, A_{2}, A_{1} and A_{0} are all n × n symmetric matrices given by

A_{2} = M_{1}− cM_{1}Mc^{−1}Mc_{1}^{>}− F_{1}G^{−1}F_{1}^{>} (19a)

= M_{1}− T_{1}M^{−1}T_{1}^{>},

A1 = K − bK cM^{−1}Mc_{1}^{>}− cM1Mc^{−1}Kb^{>}− EG^{−1}F_{1}^{>}− F1G^{−1}E^{>} (19b)

= K − SM^{−1}T_{1}^{>}− T_{1}M^{−1}S^{>},

A_{0} = bK cM^{−1}Kb^{>}+ EG^{−1}E^{>} (19c)

= SM^{−1}S^{>}.

Proof. By (16), subtracting the 1st equation by the 2nd equation in (15), and rewriting the 3rd equation in (15) we have

"

K
E^{>}

# p − λ

"

M_{1}
F_{1}^{>}

#

p = −λ

"

M_{ε}+ M_{1} F_{ε}+ F_{1}
F_{ε}^{>}+ F_{1}^{>} G_{ε}+ G_{1}

# "

u v

#

. (20)

On the other hand, from the 1st equation of (15) and the 1st equation of (20), we have h

K Ei

"

u w

#

− λh

M_{1} F_{1}i

"

u w

#

+ λM_{1}p = Kp. (21)

From (17a) and (17c), equations (20) and (21) can be rewritten as

(S − λT_{1})^{>}p = −λM

"

u w

#

, (22a)

(S − λT1)

"

u w

#

+ λM1p = Kp. (22b)

Expressing

"

u w

#

in terms of p in (22a) and plugging it into (22b), we end up with a QEP

λ^{2}M_{1}− (S − λT_{1})M^{−1}(S − λT_{1})^{>} p = λKp. (23)
Rewriting (23) as

λ^{2}(M_{1}− T_{1}M^{−1}T_{1}^{>}) + λ(−K + SM^{−1}T_{1}^{>}+ T_{1}M^{−1}S^{>}) − SM^{−1}S^{>} p = 0
and using the fact that

M^{−1} =

"

M F

F^{>} G

#−1

=

"

Mc^{−1} 0

−G^{−1}F^{>}Mc^{−1} G^{−1}

# "

I −F G^{−1}

0 I

# .

We can show by careful calculation that the coefficient matrices in (23) satisfy those in (19).

Corollary 3.3. [15] Let L(λ) and Q(λ) be defined in (15) and (18), respectively. Then σ(L(λ)) = σ(Q(λ)) ∪ {0, · · · , 0

| {z }

m

},

where σ(·) denotes the spectra of the associated matrix pencil.

Theorem 3.4. Let Q(λ) in (18) with matrices A_{2}, A_{1} and A_{0} defined in (19). Then A_{2}
and A0 are positive definite. Furthermore there are n negative and n positive eigenvalues
for Q(λ).

Proof. From Lemma 3.1 and S^{>} being of full column rank, we get that the matrix A_{0}
in (19c) is positive definite. On the other hand, from

"

M_{1} F_{1}
F_{1}^{>} G_{1}

#

0,

"

M_{ε} F_{ε}
F_{ε}^{>} G_{ε}

#

0

it follows

C ≡

M_{1} 0 F_{1}

0 0 0

F_{1}^{>} 0 G_{1}

+

0 0 0

0 M_{ε} F_{ε}
0 F_{ε}^{>} G_{ε}

=

M_{1} 0 F_{1}
0 M_{ε} F_{ε}
F_{1}^{>} F_{ε}^{>} G

0.

Let

L_{1} =

I_{n} 0 0
I_{n} I_{n} 0
0 0 I_{m}

, L_{2} =

I_{n} 0 0

0 I_{n} −F G^{−1}

0 0 I_{m}

,

L_{3} =

I_{n} − cM_{1}Mc^{−1} 0

0 I_{n} 0

0 0 Im

, L_{4} =

I_{n} 0 −F_{1}G^{−1}

0 I_{n} 0

0 0 I_{m}

. Then

0 ≺ L_{4}L_{3}L_{2}L_{1}CL^{>}_{1}L^{>}_{2}L^{>}_{3}L^{>}_{4}

= L_{4}L_{3}L_{2}

M_{1} 0 F_{1}
M_{1} M_{ε} F
F_{1}^{>} F_{ε}^{>} G

I_{n} I_{n} 0
0 I_{n} 0
0 0 I_{m}

L^{>}_{2}L^{>}_{3}L^{>}_{4}

= L_{4}L_{3}

M_{1} M_{1} F_{1}
Mc_{1}^{>} Mc 0

F_{1}^{>} F^{>} 0

I_{n} 0 0

0 In 0

0 −G^{−1}F^{>} Im

L^{>}_{3}L^{>}_{4}

= L4

M_{1}− cM_{1}Mc^{−1}Mc_{1}^{>} 0 F_{1}
Mc_{1}^{>} Mc 0

F_{1}^{>} 0 G

I_{n} 0 0

− cM^{−1}Mc_{1}^{>} In 0

0 0 I_{m}

L^{>}_{4}

=

A_{2} 0 0
0 Mc 0

0 0 G

, which implies that A2 0.

Let (λ, p) be an eigenpair of (18), then

λ^{2}(p^{H}A_{2}p) − λ(p^{H}A_{1}p) − (p^{H}A_{0}p) = 0. (24)
Here and hereafter, the superscript “H” in (24) denotes the conjugate transpose.

Because A_{1} is symmetric and A_{2}, A_{0} 0, we have

a2 ≡ p^{H}A2p > 0, a1 ≡ p^{H}A1p ∈ R, a0 ≡ p^{H}A0p > 0, (25)

which implies that the roots of the quadratic equation (24) are

λ_{+}= a_{1}+pa^{2}_{1}+ 4a_{2}a_{0}

2a_{2} > 0, λ_{−} = a_{1}−pa^{2}_{1}+ 4a_{2}a_{0}

2a_{2} < 0. (26)

Hence, there are n negative and n positive eigenvalues for (18) and the associated eigenvectors are real vectors.

Theorem 3.5. Let

W_{0} =

"

M F

F^{>} G

#−1/2"

K
E^{>}

#

, W_{1} =

"

M F

F^{>} G

#−1/2"

M_{1}
F_{1}^{>}

#

, (27)

and

d_{0} = kW_{0}k_{2}, d_{1} = kW_{1}k_{2}.
Suppose that

a_{1} = λ_{min}(K) − d_{0}d_{1} > 0, (28a)
a_{0} = λ_{min}(A_{0}), ¯a_{0} = λ_{max}(A_{0}) = d^{2}_{0}, (28b)
a_{2} = λmin(A2) , ¯a2 = λmax(A2). (28c)
Then the n negative and n positive eigenvalues of (18) are, respectively, in the intervals
(−β∗, 0) and (β^{∗}, ∞), where

β_{∗} = 2¯a_{0}

pa^{2}_{1}+ 4a_{2}a_{0}+ a_{1}, β^{∗} = a_{1}+pa^{2}_{1}+ 4a_{2}a_{0}
2¯a2

. (29)

Proof. By the definitions of W_{0} and W_{1}, A_{1} in (19b) can be represented as

A_{1} = K − W_{0}^{>}W_{1}− W_{1}^{>}W_{0}. (30)
Given an orthogonal vector p, from (27) and (28a), equation (30) implies that

p^{H}A1p = p^{H}Kp − p^{H}W_{0}^{>}W1p − p^{H}W_{1}^{>}W0p ≥ λmin(K) − d0d1 = a_{1} > 0. (31)
From (25)-(27) and (31), it follows that

−λ− = 2a_{0}

pa^{2}_{1}+ 4a_{2}a_{0}+ a_{1} ≤ 2¯a_{0}

pa^{2}_{1}+ 4a_{2}a_{0}+ a_{1} = β∗,
and

λ+= a_{1}+pa^{2}_{1}+ 4a_{2}a_{0}

2a_{2} ≥ a_{1}+pa^{2}_{1}+ 4a_{2}a_{0}
2¯a_{2} = β^{∗}.

4. Numerical method for reconstructing the unknown domain

In this section, we will propose a practical numerical method for the reconstruction of
the support of the target D based on the LSM and the truncated SVD technique. As
described in the Introduction, let z be a sampling point in the sampling domain T (see
Figure 1) and k^{2} is not a transmission eigenvalue for TEP (1). According to Theorem
1.1, from (4) and (5), we can determine whether z is inside D or not by solving v ∈ L^{2}(Γ)
which satisfies the near field integral equation

Z

Γ

u^{s}(x, y)v(y)ds(y) = Φ(x, z) for all x ∈ C, (32)
where Φ(x, z) is given in (3).

To solve (32), we first discretize (32) as follows. For each point source x_{i} ∈ C,
i = 1, · · · , m, suppose we have measured the scattered fields on the discrete points
y_{1}, · · · , y_{n} on Γ. With these scattered fields, the discretized equation of (32) can be
formulated as the overdetermined linear system

Av = b, (33)

where b = [Φ(x_{1}, z), · · · , Φ(x_{m}, z)]^{>} ∈ C^{m}, A ∈ C^{m×n} depends on x_{i} and y_{j} for
i = 1, · · · , m and j = 1, · · · , n, and the measured scattered fileds. v ∈ C^{n} is the
approximation to the unknown values [v(y_{1}), · · · , v(y_{n})]^{>}. In general, the linear system
(33) is very sensitive due to the ill-posedness of (32). For the case of m < n, the problem
(33) can be solved by the Tikhonov regularization method which gives

v∈Cmin^{n}{||Av − b||^{2}+ δ||v||^{2}} (34)
with δ > 0 be a regularization parameter. The proper choice of regularization
parameters plays an important role in (34) in order to achieve accurate and stable
numerical results. Actually, over the last four decades, many different methods for
selecting regularization parameters have been proposed [13, 19, 31]. Even so, the
selection of the optimal regularization parameter remains a difficult question in solving
the discretized ill-posed linear system (33).

In this paper, we will propose a novel numerical method to reconstruct the domain D by solving the linear system (33) for the case of m ≥ n. Let A = span(A) be the subspace spanned by the columns of A. Recall from Section1 that

dist(b, A) = min

v∈C^{n}||Av − b||_{2}.

From the perspective of matrix analysis, we reinterpret Theorem 1.1 in the following
way. Assume that k^{2} is not an eigenvalue of (18), for instance, we choose 0 < k^{2} < β^{∗}.
If dist(b, A) > 0, then v satisfying (33) is in the sense that some components of v are
infinity. In other words, we have kvk = ∞. We then assign z ∈ Ω \ D. On the other
hand, if dist(b, A) = 0, then we can find a regular vector v with kvk < ∞ satisfying

(33). In this case, we say that z ∈ D. Our reconstruction method is based on this interpretation.

We will use a truncated SVD to determine the size of v. Let A = U

"

Σ 0

#
V^{H}

be the SVD of A, where U ∈ C^{m×m} and V ∈ C^{n×n} are unitary, Σ = diag(σ_{1}, · · · , σ_{n})
is diagonal matrix of singular values with σ_{1} ≥ · · · ≥ σ_{n} ≥ 0. Usually, the discretized
linear system (33) is sensitive that is inherited from the ill-posedness of (32), and the
coefficient matrix A in (33) is not of full column rank. Suppose rank(A) = r < n, which
means that σ_{r} > σ_{r+1} = · · · = σ_{n} = 0. Denote Σ_{r} = diag(σ_{1}, · · · , σ_{r}) and partition
U and V with respect to Σ_{r} as U = [U_{1}, U_{2}] and V = [V_{1}, V_{2}] with U_{1} ∈ C^{m×r} and
V_{1} ∈ C^{n×r}, respectively. Obviously, the following statements hold true.

A = U_{1}Σ_{r}V_{1}^{H}, A = span(U_{1}) ≡ the subspace spanned by columns of U_{1}, (35)
and

dist(b, A)^{2} = min

v∈C^{n}||Av − b||^{2}_{2} = min

v∈C^{r}||Σrv − U_{1}^{H}b||^{2}_{2} + ||U_{2}^{H}b||^{2}_{2} = ||U_{2}^{H}b||^{2}_{2},
where U_{1}Σ_{r}V_{1}^{H} is called the truncated SVD of A. Let v be a vector satisfying (33).

Then we have

Av = b ⇔

"

Σ_{r} 0
0 0

# "

V_{1}^{H}v
V_{2}^{H}v

#

=

"

U_{1}^{H}b
U_{2}^{H}b

#

≡

"

Σ_{r} 0
0 0

# "

ˆ
v_{1}
ˆ
v_{2}

#

=

" ˆb_{1}
bˆ2

#

. (36)

We interpret that v = ∞ is the solution of the equation 0 v = c if c 6= 0. Thus, from (36) follows that

Σ_{r}vˆ_{1} = ˆb_{1} ⇒ ˆv_{1} = Σ^{−1}_{r} bˆ_{1} ⇒ ||ˆv_{1}|| < ∞,
and

0ˆv_{1}+ 0ˆv_{2} = ˆb_{2} ⇒ if ||ˆb_{2}|| 6= 0, then ||ˆv|| = ∞.

This implies that

if b ∈ span{U_{1}} ⇔ rank(A) = rank([A, b]), then ||v|| < ∞; (37a)
if b /∈ span{U1} ⇔ rank(A) 6= rank([A, b]), then ||v|| = ∞. (37b)
On the other hand, let bA = [A, b], then

Ab^{H}A =b

"

V1Σr 0
bˆ^{H}_{1} bˆ^{H}_{2}

# "

Σ_{r}V_{1}^{H} bˆ_{1}
0 bˆ_{2}

#

=

"

V_{1}Σ^{2}_{r}V_{1}^{H} V_{1}Σ_{r}bˆ_{1}
bˆ^{H}_{1} Σ_{r}V_{1}^{H} ||b||^{2}_{2}

#

. (38)

It is easy to see that bA^{H}A in (38) is similar tob

Σ^{2}_{r} Σ_{r}bˆ_{1} 0
bˆ^{H}_{1} Σ_{r} ||b||^{2}_{2} 0

0 0 0

≡

"

Σb^{2} 0
0 0

#

(39)

by similarity transformations V ⊕ I and 0 ⊕

"

0 I I 0

#

. Here “ ⊕ ” denotes the direct sum
of matrices. Let ˆσ_{1} ≥ · · · ≥ ˆσ_{r+1} ≥ 0 be eigenvalues of bΣ. It follows immediately that
the singular values of bA are ˆσ_{1} ≥ · · · ≥ ˆσ_{r+1} ≥ ˆσ_{r+2} = · · · = ˆσ_{n+1} = 0.

In view of the special structure of matrix bΣ in (39), we apply the interlacing theorem for singular values of Σ and bΣ which leads to

ˆ

σ_{1} ≥ σ_{1} ≥ ˆσ_{2} ≥ σ_{2} ≥ · · · ˆσ_{r} ≥ σ_{r} ≥ ˆσ_{r+1} ≥ σ_{r+1} = 0. (40)
For convenience, we set 0/0 = 1. Let σn+1= 0, then we have

1 ≤ ˆσj/σj < ∞ for j = 1, · · · , r, and ˆσj/σj = 1 for j = r + 2, · · · , n + 1. (41) Combing (40) and (41), the statements (37) can be represented as

b ∈ span{U_{1}} ⇔ ˆσ_{r+1} = 0, thus ˆσ_{r+1}/σ_{r+1} = 0/0 = 1; (42a)
b /∈ span{U1} ⇔ ˆσr+1 6= 0, thus ˆσr+1/σr+1 = ˆσr+1/0 = ∞. (42b)
In conclusion, in our method, to decide whether z ∈ D or not relies heavily on
effectively determining the ranks of A and bA. However, since the integral equation (33)
is ill-posed, the matrices A and bA are normally ill-conditioned and there is often no
gap in the spectrums of singular values. From numerical point of view, it is difficult to
compute the truncated SVD for A. Consequently, we usually set r = n.

We now propose a practical criterion to numerically realize (42). As above, let the
singular values of A be σ_{1} ≥ · · · ≥ σ_{n} ≥ 0 and those of bA be ˆσ_{1} ≥ · · · ≥ ˆσ_{n+1} ≥ 0. We
first compute

j(z) ≡ arg max_{1≤j≤n}σˆ_{j}
σ_{j}
and define

I_{z}= ˆσ_{j(z)}.

Here I_{z}is used as an indicator to determine whether z is in D or not. Precisely, we pick
a small threshold parameter 0 < M. Then we choose the following dichotomy:

if Iz< M, then we set z ∈ D,
if I_{z}≥ M, then we set z /∈ D.

The algorithm of computing the indicator I_{z} is summarized in Algorithm 1 below.

5. Numerical experiments

In what follows, we will demonstrate the efficiency of our numerical method for four different domains shown in Figure 2: (a) a disk centered at (0, 0) with radius 0.5; (b) an ellipse region centered at the origin with axes 0.6 and 0.4; (c) a peanut-like region

Algorithm 1 Practical numerical reconstruction method based on the singular values Input: Discrete point sources xi ∈ C, detection points yj ∈ Γ for i = 1, · · · , m and

j = 1, · · · , n, respectively, with m > n. The wave number k ∈ R and sampling points z ∈ T .

Output: The indicator I_{z}.

1: For each point sources x_{i}, collect all the measured scattered field u^{s}_{ij} for i =
1, · · · , m, j = 1, · · · , n.

2: Using the measured scattered field u^{s}_{ij} to discretize the integral of (32), then
construct the coefficient matrix A ∈ C^{m×n} of (33).

3: Compute the SVD of A = U_{1}Σ_{n}V_{1}^{H} as in (35) with Σ_{n} = diag(σ_{1}, · · · , σ_{n}).

4: For the sampling point z, compute the vector on the right hand side of (33),
b = [Φ(x_{1}, z), · · · , Φ(x_{m}, z)]^{T}, where Φ(·, ·) is given in (3).

5: Let bΣ ≡

"

Σ^{2}_{n} Σ_{n}bˆ_{1}
bˆ^{H}_{1} Σ_{n} ||b||^{2}_{2}

#

with ˆb_{1} = U_{1}^{H}b.

6: Calculate the singular values of bΣ as ˆσ_{1} ≥ · · · ≥ ˆσ_{n+1}.

7: Determine the index j(z) = arg max_{1≤j≤n}σˆ_{j}/σ_{j}.

8: Set I_{z}= ˆσ_{j(z)}.

(a) Disk (b) Ellipse (c) Peanut (d) Heart

Figure 2. Four model domains that represent the region D.

enclosed by the equation ((x_{1}− 0.7)^{2} + x^{2}_{2})((x_{1}+ 0.7)^{2}+ x^{2}_{2}) = 0.72^{4}; (d) a heart-like
region enclosed by the equation x^{2}_{1}+ x^{2}_{2}+ 0.5x_{2} = 0.5px^{2}_{1}+ x^{2}_{2}.

All computations for numerical test examples are carried out in MATLAB 2017a.

For the hardware configuration, we use an HP server equipped with the RedHat Linux operating system, two Intel Quad-Core Xeon E5-2643 3.33 GHz CPUs and 96 GB of main memory.

5.1. Distribution of the transmission eigenvalues

We carry out the FEM method described in Section 3to the TEP (1) on four different
domains as in Figure2with the regular mesh size h ≈ 0.004 for triangles of each domain
D. We set ε_{r} = 16, γ_{1} = 10, γ_{2} = 4, and thus ε(x) = −100. The associated dimensions
n and m of matrices given in Table1 are shown in Table2. For each domain, we derive
a corresponding QEP as given in (18).

Because of the two negative signs in (18), we should modify the quadratic Jacobi-

Table 2. Dimensions n, m (K ∈ R^{n×n}, E ∈ R^{n×m}) of matrices for the benchmark
problems with the mesh size h ≈ 0.004.

Domain Disk Ellipse Peanut Heart

(n, m) (124631, 1150) (71546, 976) (149051, 1871) (168548, 1492)

-80 -40 0 40 80 120 160 200 240

Disk Ellipse Peanut Heart

17.7 30.3

28.3 100.3

Transmission eigenvalues

Figure 3. The eigenvalues λ of the QEP (18) in the intervel [−80, 250] for the four domains in Figure2with ε(x) = −100. The arrows point to the first positive eigenvalue of the QEP corresponding to each domain.

Davidson method in [23] to solve these QEPs. Of course, the partial locking and partial deflation schemes of Algorithm 2 and Algorithm 3 in [23] can be employed as well. The numerical results are shown in Figure 3.

Owing to the negative inhomogeneous medium ε(x) = ε_{r}− (γ_{1}^{2}+ γ_{2}^{2}) in TEP (1)
derived from the pesudo-chiral model. From Figure3, we immediately observe that there
indeed exists an eigenvalue-free interval (0, β^{∗}) = (0, 100.3) (disk), (0, 28.3) (ellipse),
(0, 30.3) (peanut) and (0, 17.3) (heart), which verifies Theorem 3.5. Consequently, it is
legitimate to say that k^{2} is not a transmission eigenvalue if k^{2} ∈ (0, β^{∗}). In comparison,
we note that the TEP (1a) with the Tellegen model in [23] is reduced to

−∆E_{3} = ω^{2}(ε_{r}+ (γ_{1}^{2}+ γ_{2}^{2}))E_{3} ≡ ω^{2}ε(x)E_{3} (43)
with a large positive inhomogeneous medium ε(x), and its transmission eigenvalues are
densely distributed on the interval (0, O(1)).

5.2. Reconstruction of the unknown domain from the near-field measurements

In this subsection, we apply Algorithm1to reconstruct the domain D from the near-field measurements for four distinct shapes as in Figure 2. As described above, our method is based on the LSM and the SVD technique.

For LSM, as shown in Figure 1, we make the following preparations for numerical experiments. Consider each domain in Figure 2 to be the target D, and let the circles

with radii 3 and 6 be Γ and C, respectively. Choose the rectangle domain [−1, 1]×[−1, 1]

to be the sampling domain T , which contains all possible targets D. From Figure 3, we choose a k ∈ (0,√

β^{∗}), in other words, k^{2} is not a transmission eigenvalue. Different
k ∈ (0,√

β^{∗})’s are chosen for testing in our numerical simulations.

Divid C and Γ into m and n segments uniformly, and denote the nodes as x_{1}, · · · , x_{m}
and y_{1}, · · · , y_{n}, respectively. Place a uniform grid on T by drawing vertical and
horizontal lines through the points with coordinates (x_{1i}, x_{2j}), where x_{1i} = −1 + ih_{1}
for i = 1, · · · , p and x_{2j} = −1 + jh_{2} for j = 1, · · · , q, respectively. Let each mesh point
(x_{1i}, x_{2j}) to be the sampling point z_{ij}. In our experiment, for these parameters, we
choose m = 1269, n = 693 and p = q = 201, respectively.

To construct the discretized near-field integral equation, we need to collect all
the scattered fields u^{s} on the detection points y_{1}, · · · , y_{n} for each point source x_{i},
i = 1, · · · , m. For the purpose of numerical experiments, for a given point source x_{i},
we can obtain the scattered fields in row [u^{s}_{i1}, · · · , u^{s}_{in}] by solving the direct scattering
problem (2) using the FEM in Appendix A. Here, the FEM is applied to the domain
enclosed by C, and the corresponding dimensions of E in Table 1 are n = 143834 and
m = 1269, respectively. Adding 3% noise to the computed scattered fields u^{s} to produce
the detected waves, with which, we discretize (32) into an ill-posed overdetermined linear
system (33) with A = [δ_{Γ}u^{s}_{ij}(x_{i}, y_{j})]_{m×n} ∈ C^{m×n}, where δ_{Γ} = 6π/n is the arc length of
the uniform segment of Γ .

Applying Algorithm1to each sampling point zij ∈ T with all the above parameters to obtain the corresponding indicator Izij for i = 1, · · · , p and j = 1, · · · , q. From numerous numerical experiments by tuning k ∈ (0,√

β^{∗}), we find that the reconstruction
results will be heavily affected by the wave number k. In Figure 4, we show the 3d
surface figures and the corresponding 2d contour figures by plotting I_{z}_{ij} with respect
to (x_{i}, y_{j}) ∈ T , and the actual targets in Figure 2 are plotted by red lines. In
fact, Figure 4 shows the best results for reconstructing targets listed in Figure 2 by
choosing appropriate wave numbers k^{∗} ∈ (0,√

β^{∗}). In Figure 5, we show the 2d
contour figures for other wave numbers k ∈ (0,√

β^{∗}). According to our experimental
experience, the areas of the reconstructed domains are less than the areas of the actual
domains by approximately 10% − 20%. Our numerical experiments also suggest that
the reconstructions of convex domains are more stable with respect to the choice of the
wave number k and noise.

5.3. Tellegen model – positive index of refraction

In comparison, we demonstrate some numerical results of reconstructing the unknown target D for the inverse scattering problem with positive index of refraction. In [24], we considered the Maxwell’s equations with the TM mode in Tellegen media and derived the acoustic equation (1a) with a positive ε(x) increasing with respect to the Tellegen parameters. When ε(x) is large enough, we have shown that the transmission eigenvalues are densely distributed near the origin. In this subsection, we show some numerical

(a) Disk

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

(b) Disk

(c) Ellipse

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

(d) Ellipse

(e) Peanut

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

(f) Peanut

(g) Heart

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

(h) Heart

Figure 4. Reconstruction results of four targets in 3d surface figures and in 2d contour
figures with ε(x) = −100 and k^{∗}∈ (0,√

β^{∗}). The domains enclosed by red curves are
the exact targets D. The scattered fields used have 3% noise.

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -1

-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

(a) Disk

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

(b) Ellipse

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

(c) Peanut

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

(d) Heart

Figure 5. Reconstruction results of four targets in 2d contour figures with ε(x) =

−100 and the different k ∈ (0,√

β^{∗}). The domains enclosed by red curves are the exact
targets D. The scattered fields have 3% noise.

reconstruction results based on our method when k^{2} is an eigenvalue of the associated
QEP (presumably a transmission eigenvalue). We test two different kinds of index of
refraction on the disk domain, one is a normal index of refraction ε(x) = 16 [24] and the
other is ε(x) = 500 (Tellegen model) [23]. As shown in [24], when ε(x) = 16, the lowest
transmission eigenvalue is approximately 1.988. In Figure 6(a), 6(b), we see that the
reconstruction results are very sensitive to the choice of the wave number k. For the case
of ε(x) = 500, our numerical result in [23] indicates that the transmission eigenvalues
are distribute densely in the interval (0, 50). Figure 6(c) (choosing k = 5) shows that
the disk is completely unrecognizable. The numerical results in this subsection verify
the need of avoiding the transmission eigenvalues in the LSM.

6. Conclusion

In this paper, we propose a practical numerical method based on the LSM and the truncated SVD to reconstruct the support of the inhomogeneity in the acoustic equation with negative index of refraction, resulting from in pseudo-chiral media. It turns out

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -1

-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

Disk domain with (x)=16 and k=2

(a) Reconstruction result of a disk (k = 2)

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

Disk domain with (x)=16 and k=2.1

(b) Reconstruction result of a disk (k = 2.1)

(c) Reconstruction result of a disk (large index of refraction)

Figure 6. Numerical reconstruction results of a disk when k^{2} is a transmission
eigenvalue. The scattered fields have 3% noise.

the index of refraction is in the form ε(x) = ε_{r} − (γ_{1}^{2}+ γ_{2}^{2}). We are also interested in
the distribution of transmission eigenvalues for the corresponding TEP. The associated
discretized TEP can be reduced into a QEP with symmetric coefficient matrices whose
all the nonphysical zero eigenvalues are deflated. We prove that the corresponding QEP
has half of negative eigenvalues and half of positive eigenvalues in (−β∗, 0) and (β^{∗}, ∞),
respectively, and there exists an eigenvalue-free interval (0, β^{∗}). We then apply the LSM
using the near-field data to reconstruct the domain D with ε(x) 6= 1 and propose an