meals I

ng7 .UBRARY
MIC Sig-:15. State
University

 

 

 

 

This is to certify that the
dissertation entitled

STUDIES OF NONLINEAR PROBLEMS FOR MAXWELL’S
EQUATIONS

presented by

YING LI

has been accepted towards fulfillment
of the requirements for the

Ph.D degree in Mathematics

 

 

l‘ \(“" é 6gp

 

Major Professor’ 8 Signature
07/: 4 / 2p 0 I

Date

MSU is an affinnative-action, equal-opportunity employer

 

..--—--o-----.—---n-o-a-n-o-n-u---I-u--I--o---l-I--c-I-o_-|-o----u-u-n-u-n-I-o-u-n-u-u—u-n-n-I-n-u-----n-oon—n-

l PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DAIEDUE

DAIEDUE

DAIEDUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/07 p:/CIRC/DateDue.indd-p.1

 

STUDIES OF NONLINEAR PROBLEMS FOR
MAXWELL’S EQUATIONS

By

Ying Li

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Mathematics

2007

ABSTRACT

STUDIES OF NONLINEAR PROBLEMS FOR
MAXWELL’S EQUATIONS

By

Ying Li

Consider the electromagnetic field scattered by a nonlinear optical medium. Be-
cause of inhomogeneity of the medium, the governing equations are Maxwell’s equa-
tion with jump coefficients and a source term. By using the Sommerfeld radiation
condition, the model scattering problem may be truncated into a bounded domain.
In this paper, LP estimates for Maxwell’s equation are established. The solution of
Maxwell’s equation is represented by spherical harmonics. LP estimate is for the
Maxwell equations with jump coefficients. An application of our LP estimates gives
rise to the wellposedness of a linearized model.

In part two, an adaptive finite element method is developed for solving Maxwell’s
equations in a nonlinear periodic structure. The medium or computational domain
is truncated by a perfect matched layer (PML) technique. Error estimates are
established. Numerical examples are provided, which illustrate the efficiency of the
method.

In part three, an inverse scattering problem is formulated for breast cancer de-
tection. A recursive linearization algorithm is used to solve the inverse scattering
problem. We employed the idea of finite element boundary integral method and
added suitable boundary conditions on the surface of the breast and an artificial
boundary which encloses the tumor. Finite element method is used for the inte-
rior domain containing inhomogeneity. Nystr'o'm method is used for the integral

equations and exterior domain. Numerical examples are presented.

To my beloved grandma and grandpa.

iii

ACKNOWLEDGMENTS

I would like to convey my gratitude to all those who gave me the possibility to
complete this thesis.

My deepest thanks go to my supervisor, Professor Gang Bao, who guided me into
the field of research, and shepherded me through the bulk of the work. His scientific
insight has had a remarkable influence on my entire career. I am very grateful for
his enthusiastic support and helpful advice throughout the time it took to compile
the work, and his unwavering belief in my ability to successfully complete the task.

I am sincerely thankful to Professor Zhengfang Zhou, for his patience, motiva-
tion, enthusiasm, and immense knowledge in PDE theories and other aspects in
applied mathematics, and for his untiring help during my difficult moments. His
detailed and constructive comments were vital to the development of this thesis.

Special thanks go to Professor Haijun Wu from Nanjing University of China,
for his patient help during the time he stayed in MSU, which led to my better
understanding of the theory and algorithm used in this study. Many thanks are
also due to Dr. Peijun Li from University of Michigan, who gave me invaluable help
including answering my numerous questions.

I would like to thank the members of my thesis committee, Professor Chichia
Chiu, Professor Tien-Yen Li, Professor Di Liu, and Professor Promislow, for the time
they spent on serving on my committee. I wish to express my special thanks to my
friend, Xiaoyue Luo, for her friendship, emotional support and helpful discussions.

I am forever indebted to my dearest grandparents, who brought me up and were
always proud of me, and my Aunt, Huichun Liu, who gave me help whenever I
needed. Finally, my heartfelt thanks go to my husband, Gang Tang, for his love,
support and inspiring ideas during my PhD study, and my daughter, Mingming, for

showering me with her unconditional love through all my moods.

iv

TABLE OF CONTENTS

LIST OF FIGURES
LIST OF TABLES
Introduction

1 Lp estimate Of Maxwell’s equations in a bounded domain
1.1 Introduction ................................
1.2 Spherical Harmonics ...........................
1.3 Boundary Estimate ............................
1.4 Existence and Uniqueness ........................

2 Numerical Solution of Nonlinear Diffraction Problems
2.1 Introduction ................................
2.2 Modeling of the Nonlinear Grating Problem ..............
2.3 Variational Formulation .........................
2.4 PML Formulation .............................
2.5 Discrete Problem .............................
2.6 Numerical Examples ...........................

3 Inverse Medium Scattering in Breast Cancer Detection

3.1 Introduction ................................
3.2 Analysis of the Scattering Map .....................
3.3 Inverse Medium Scattering ........................

3.3.1 Born Approximation .......................

3.3.2 Recursive Linearization ......................
3.4 Implementation ..............................
3.5 Numerical Experiments ..........................

BIBLIOGRAPHY

vi

vii

13
20

23
23
25
29
32
4O
47

50
50
56
63
64
66
70
72

81

1.1
2.1
2.2
2.3
2.4
2.5
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
3.10
3.11
3.12
3.13

3.14
3.15

LIST OF FIGURES

Geometry of the scattering problem ................... 8
Geometry of the grating problem .................... 27
Geometry of the PML problem ..................... 33
ZnS overcoated binary silver gratings .................. 47
Groove depth and enhancement ..................... 48
Second harmonic enhancement ...................... 49
Geometry of the inverse scattering problem ............... 55
Dielectric properties at frequencies described by Debye model [42] . . 56
Geometry of the inverse scattering problem ............... 63
Real part of smooth scatterer function ................. 73
Imaginary part of smooth scatterer function .............. 74
Born Approximation of the real part of smooth scatterer function . . 74
Born Approximation of the imaginary part of smooth scatterer function 75

Final construction of the real part of smooth scatterer function . . . . 75

Final construction Of the imaginary part of smooth scatterer function 76

Real part of piecewise scatterer function ................ 76
Imaginary part of piecewise scatterer function ............. 77
Born Approximation of the real part of piecewise scatterer function . 77
Born Approximation of the imaginary part of piecewise scatterer func-

tion .................................... 78
Final construction of the real part of piecewise scatterer function . . . 79
Final construction of the imaginary part of piecewise scatterer function 79

vi

LIST OF TABLES

3.1 Typical dielectric properties of various tissues in the breast [27] . . . 51

vii

Introduction

This research focuses on mathematical and computational studies of second har-
monic generation in electromagnetism and optics. Second harmonic generation
arises from the nonlinearity of optical materials. When a plane wave with frequency
w is incident on a nonlinear structure, the nonlinearity of the structure gives rise
to the scattered waves at both frequency w and 2w. This important phenomenon
is known as Second Harmonic Generation (SHG) in nonlinear optics. A significant
application of SHG is to obtain coherent beams of light in parts of the spectrum
at which lasers cannot be made and to construct optoelectronic devices based on
nonlinear effects in waveguides and Optical fibers.

In chapter one, we established the uniqueness and existence of the solution of
Maxwell’s equations in a bounded domain containing a nonlinear medium. Consider
the electromagnetic field scattered by the nonlinear medium. Because of inhomo-
geneity of the medium, the governing equations are Maxwell’s equation with jump
coefficients and a source term. By using the Sommerfeld radiation condition, the
model scattering problem may be truncated into a bounded domain. The solution
of Maxwell’s equation is represented by spherical harmonics. LP estimate for the
equations are established, which gives rise to the wellposedness of a linearized model.

In practice, the SHG optical effects are often too weak to be observed. Therefore,
modeling and enhancement of SHG are of great interest to potential real applica-

tions. It is pointed out that the SHG can be greatly enhanced in periodic structures.

In chapter two, questions on the existence and uniqueness have been studied. W’e
have developed an adaptive finite element method for solving the model scattering
problem. The medium or computational domain is truncated by a perfect matched
layer (PML) technique. Error estimates are established. Numerical examples are
provided, which illustrate the efficiency of the method. Numerical solution of the
nonlinear model problem in three dimensions is completely Open, which will be one
of my future projects.

In chapter three, we consider the inverse scattering problem arising from breast
cancer detection. The problem is to determine the dielectric property of the tissues
from the measurements of electromagnetic field on the surface, given the incident
field. In addition to the ill-posedness and nonlinearity of the inverse scattering
problem, one major difficulty lies in the multiple scales of the problem. The tumor
is comparably small in the computational domain, which makes the computation
challenging. Another difficulty is due to the dispersive nature of the human body.
Maxwell’s equations in dispersive media must be studied. A continuation method
is developed for this problem. The algorithm needs multi-frequency Dirichlet and
Neumann data on the surface. The initial guess comes from Born approximation.

The dielectric constant is updated by using higher and higher wavenumber k.

CHAPTER 1

LP estimate Of Maxwell’s

equations in a bounded domain

1.1 Introduction

Second harmonic generation (SHG) is a well known nonlinear optical effect. It was
first demonstrated by P. A. Franken, A. E. Hill, C. W. Peters, and G. Weinreich [33]
in 1961. The demonstration was made possible by the invention of laser in 1960,
which created the required high intensity monochromatic light. In the experiment,
they focused a ruby laser with a wavelength of 694 nm into a quartz sample. They
sent the output light through a spectrometer, recording the spectrum on photo-
graphic paper, which indicated the production of light at 347 nm. The physical
mechanism behind SHG can be understood as follows. Due to the nonlinearity, the
incident (pump) wave generates a nonlinear polarization which oscillates with twice
the fundamental frequency. According to Maxwell’s equations, this nonlinear po-
larization radiates an electromagnetic field with this doubled frequency. The latter
also interacts with the fundamental wave, so that the pump wave can be attenuated

(pump depletion) when the second harmonic intensity develops. Energy is trans-

ferred from the pump wave to the second harmonic wave. It is a very important
nonlinear optical effect because theoretically coherent beams of light can be obtained
in parts of the spectrum at which lasers cannot be made and optoelectronic devices
based on nonlinear effects in waveguides and Optical fibres can be constructed. The
reader is referred to [52] for detailed descriptions of nonlinear optics.

A PDE model was introduced in [48], [49] and [50] to describe nonlinear SHG in
periodic structures. In Bao, Minut and Zhou [18], the regularity is studied for the
solutions of Maxwell’s equations with source term in a domain with jump dielectric
coefficients. These Lp estimates are further employed to solve the linearized SHG
problem in a periodic structure. In this paper, we study the regularity of Maxwell’s
equations for the scattering by a bounded domain of a nonlinear medium. Although
the interior estimate is the same as in [18], the boundary estimate requires a new
technique. A striking difference is due to the decay rate of the fields away from the
medium. In the periodic case [18], the fields decay exponentially, which makes the
L()0 norm estimate easier. For the scattering from a bounded medium, because of
the slow decay of the fields away from the medium, we must use fine properties of
the Hankel functions and spherical harmonics.

For simplicity we assume the medium is nonmagnetic (p. E [1.0) and no external

current or charge is present in the field. The following Maxwell’s equations hold:

as -
__: E

at vx ,
013 ..
‘37—”va

where E is the electric field, H is the magnetic field, B is the magnetic induction,

and D is the electric induction. The constitutive equations are:

[3'

Hoff,

U1
ll

CUE-F f3,

where #0 is the constant magnetic permeability, 60 is the dielectric permittivity in
vacuum. P is the polarization.
The time harmonic solutions of Maxwell equations, also called plane waves,

are complex-valued fields:

that satisfy the system of time-harmonic Maxwell equations:

—o

V XE = —z'w,uf_f,

where w stands for the frequency of the electromagnetic waves.

In the linear case, the polarization is induced linearly by the electric field:

where XU) is called the linear susceptibility tensor. Thus, if a beam of angular fre-
quency w is passing through the medium, a polarization oscillating at w is produced
which in turn serves as a source for the further propagation of the original wave at
w.

If the light intensity is high (as lasers), the nonlinear effect will play a role and

the polarization will depend nonlinearly on the electric field [52]:
13 = X(1)E + x(2)E2 + X(3>E3 +

where Xm is the i-th order susceptibility tensor. Throughout, we restrict our at-

tention to the 2nd order susceptibility by ignoring the higher order terms. Hence

"_ (1) (2) .
P“ 2 X-j Ej+ Z X-jkEjEk’
3.213233 jak:1a2a3

where 132' is the i-th component of E It is clear that the second term in the

polarization may generate fields of frequency 2w. Let
Em) = Eéjw)(x)e—-ijwt + 56(jw)(x)ezjwt_
We can write the total field as
E‘ = EM + [35(2‘”).

Then, by omitting eiiwt terms,

EjEk = (E83?)($)e—zwt +E(()j)($ ) ezwt+E$w)($ )e—i2wt+E(()2W)(x ) ei2wt).

(Eé:)($)e—iwt+E(()k)(x ) eth+E(()i (”(13% )e —i2wt+E(giw)(x )6 iZwt)

w —zw w 2w to) w 2w
_E(()j)E(()k) 2t E<j>E P3681.) t+EélEék) 2t

(“1 ) (zwle —zwt+ (2w) (w ) —iwt (2w) (w )z'wt
+E0j EOke +E0j EOke +E0j EOke .

It follows that

Pl‘ 2X2)“; (fie—taut +E(()j)e iwt+E83w)e _z'2wt+E(()2w)e i2wt)

(2) (WE (301:) —i2wt (w ) (2W) eiwt (w ) (w) ei2wt
+JZX Xl]k( (E0 6 +on EOke +on EOke

{(91 ) (2w) —twt (2w)—(w) —iwt (20’) (w ) iwt
+5703} EOk e +E0j EOk e +E0j EOk e )

by ignoring the 3w and higher order terms. Evidently, the combination of two fields
of frequency w generates the second harmonic field.
The time-harmonic Maxwell’s equations become (by dropping subscript O to

simplify the notation)

V XEW) = —'iw,u.0ff:w),

v xI-I(w)= zchEW )+ iwjzk(2 x_ jl)c(E Elia”) + E§2w)E(w)),

“(2&1)

V XE H(2u1)

= —'i2tu;10 ,

~ 9 . -' . . 2 w w
v wa) = 22..)(0EW) +22kaxlng§ )EIE ).
J,

The nonlinear polarization may be treated as a source term in lV'laxwell’s equations.

Rewrite the equations as
V XE = —iwrtfi,
v xH=iweE+§, (1.1.1)

where g’ is the source term.

The rest of this paper is organized as follows. In the next section, spherical
harmonics are introduced to represent the magnetic field. The boundary condition
is derived on the artificial boundary S R' Section 3 is devoted to establishing the
LP estimate on S R' The wellposedness of the model problem is proved in Section 4

and Section 5.

1 .2 Spherical Harmonics

Consider a bounded nonlinear medium enclosed by a boundary surface S. Assume
that the dielectric coefficient is 61 inside S; 50 outside of S. It is assumed to be
vacuum outside the medium.

N ow let S R be the sphere of radius R such that S R encloses the whole nonlinear
medium. Outside the medium, the electromagnetic fields satisfy (1.1.1) with g = 0.
For simplicity, the arrows are omitted throughout this section. Taking curl of the

second equation and eliminating the electric field E gives:
v x (v x H) = —w2ueH.
By employing the vector identity:

VX<VxA)=v(v-A)—(v-V)A.

Figure 1.1. Geometry of the scattering problem

Nonlinear Medium 8:81

 

and

vlvxm=m
we obtain the Helmholtz equation:

AH+kfir=a (ran

where k = w2

err is the wave number in vacuum.
To ensure the uniqueness of the solution, the following Sommerfeld’s radiation

condition is imposed:

H
lim r(%— — ikH) = 0.

r —> 00 r
We will now use spherical harmonics to represent the solution of this equation.

Spherical harmonics are the angular portion of an orthogonal set of solutions to

Laplace’s equation represented in a system of spherical coordinates. The readers
are referred to [45] and [1] for detailed discussions of spherical harmonics.
In spherical coordinates, equation (1.2.1) becomes

18(7‘2%—H

2
7287‘ )+ TizASUH+kH= 0,

where

1 a2 1 a 8
Aso _ 5,395,? + EM "1600)

is the Laplace-Beltrami operator on the unit sphere SU. The Hermitian product in

L2(SU) is given by

27r
/ uvdo— — / [flaw t(,6 (1)) sin 6d6d¢.
SU

Let H1(SU) be the Hilbert space

H1<s> = {u 6 flag), vsUu e (L2(SU)>2}

with its Hermitian product

u,v =— uvdo+ V av odo.
( )H1(SU) 4 S S( S S )

The Laplace—Beltrami operator is self-adjoint in the space L2(S). It admits a family
of eigenfunctions which constitutes an orthogonal basis for the space L2(S). This
basis is also orthogonal for the scalar product in H1(S). These eigenfunctions are
called spherical harmonics.

Denote by 371 the space of homogeneous polynomials of degree I that are har-
monic, with restrictions to the unit sphere. We list the following theorem from [45]

without proof:

Theorem 1.2.1. Let Ylm, —l g m S 1, denote an orthonormal basis of 371 for the
hermitian product of L2(S). The functions Ylm, forl Z O and —l S m g l, consti-

tute an orthogonal basis in L2(S), which is also orthogonal in H1(S). Moreover, yl

coincides with the subspace spanned by the eigenfunctions of the Laplace-Beltrami

operator associated with the eigenvalue —l(l + 1), z e
ASYlm +l(l+1)Ylm = 0,
and the eigenvalue —l(l + 1) has multiplicity 2l + 1.

The spherical harmonics of order l are the 21 + 1 functions of the form:
1
l +— ( l _ m)!

1
imgb m cos
27r2 (l+m)l '128 Pl ( (0))

 

where an(cos(6’)) are associated Legendre functions. They have the following

properties:

—m mll‘m! m
“”1 2H) (l+ml!Pl’

WW, ¢> = sums—me, a

 

where the superscript * denotes complex conjugation. The spherical harmonics form
a complete set of orthonormal functions and thus form a vector space analogous to
unit basis vectors. On the unit sphere, any square integrable function can thus be

expanded as a linear combination of these:

=2 2 flnfi’W»

l=0m=—l

The expansion coefficients can be obtained by:

27r 7r
17” = (2f(9.¢mm*(9.¢)do= f0 ems/0 d6sin6f(6.¢mm*(6.¢>.

We now expand the magnetic field by spherical harmonics. Let

Z

Hosts): 2191(7) Z lilelm (0,95).
120 mz—l

10

Suppose that on the boundary SR:

we»: zzu

 

 

l— — Om— — —
for some constants ugn. Substituting this solution into the Helmholtz equation, we
have
d2Hl 2d_H, 2 1(1+ 1)
k — H 2
dr2 7‘ (17‘ ( r2 ) l 0,

i.e., the spherical Bessel equation. It can be transformed by rescaling to

1(1+1)

d2H 2dH
___l+__l+(1_ Tml =0, (1.2.2)

(172 r dr

A useful Lemma from [45] about the spherical Bessel equation is stated here without

proof.

Lemma 1.2.1. The spherical Bessel equation (1.2.2) admits two families of solu-

tions, known as spherical Hankel functions, which satisfy the recursion formulas

 

d 1 1+1
(1,111: ;:Hz H1+1=——H1+Hl—
2z+1

Hz+1+Hz—1= ,. Hz-

The spherical Hankel functions are given by the expressions

 

1 1d ei
h§ )m = <—r>‘<;;,;>’<7>,
2 1d 6‘”
h) )(r>=<—r>’(;;,;>’< 7, );

more specifically by

ir

h1)(7)=(—i)le——T (”15+ U31: +-- -l+()lfil(‘1‘ )1)

l
(m + l)!
m!(l — m)!2m'

 

(

l
,,§2)( ):h<1>( )
flirt =

11

The function

satisfies the recursion formula
(zl _ 1 — (z — 1))(zl +z+ 1) = —r2.
Moreover,

1 1 1
1+ oar—2 + + 05m = r2|h§ )(r)|2, (1.2.3)

l

am = flinfliii

In addition

Writing the solution as:

Hl(r) = 711h§1)(kr) + ’yl‘ghl(2)(kr).

By looking at the explicit forms of spherical Hankel functions:

 

hl(1)(r):(_z tr)“:— ([30+ifl1:+ +(lllfill(‘1‘)la)
hi

we obtain:

(:)_H_ —ik.=H Z Z [[k7l1h1( —h(1)(kr)—ihl(1)(kr)),

l=0m=—

+ k7'2h;"‘(%hl(2) (kr) — ili§2)(kr))]Ylm(6, e).

12

From Theorem 1.2.1, we know that

 

 

 

 

%hl(1)(k7") — z"1(1)(kv) N 1'3()((—z illmelkggrg eikr — i(—z')lei:cr)
= fi6(-(- )1 :71“;

219%?th _ MPH”) N fié((i)lflik€_i:;:2‘ e—z'kr _ i(i)le—:kr)
= Btu-flit + 13—71? — 11%).

But by the Sommerfeld’s radiation condition, we need Mfg—H — ikH) —+ 0 as r ——> 00,
r

which concludes that 712 = O.

Plugging in the boundary condition, the solution becomes

H(r,9,¢) = 20 m2— umYl’” (6 95) (1.2.4)
1: = —lhl(1)(kR)
On S R1 we have
00 l
giFISR =2 Zifi (kR) )ul MY; (9 as). (1.2.5)

=0m:—

It follows that the problem may be truncated into a bounded domain with the

boundary S R and a boundary condition (1.2.5) on S R'

1.3 Boundary Estimate

Our goal is to establish global estimates for the solutions of the scattering problem.
We first present a local estimate from [18] which provides the LP estimate inside
the dielectric and on the interface. Throughout the paper, C stands for a positive
generic constant whose value may vary step by step but should always be clear from
the context.

Let 1< p < 00, let B be an open ball in R3 and let 97 E Lp(B). Let E E Lp(B)
and H E W17p(B) be a solution of (1.1.1).

13

Let S be a C2 surface embedded in R3 such that S divides B into two connected

components B+ and B‘. Assume the electric permittivity c is defined by
6+ in B+,
6" in 8‘.

Theorem 1.3.1. (Local estimate) For any 3’ with B' C B,

IIEIILP(BI) + llHHW1,p(BI) S C(llfillLP(B) + llglllfvug) +“EHW_1’p(B))’

where C is a constant depending only on p, B’ and B.

Therefore, in order to establish global estimates, it suffices now to obtain LP
estimates for the solutions of Maxwell’s equations on SR. Our main result in this

paper is:

Theorem 1.3.2. Let Q = {:cIR — (5 < |:L‘| g R}, 1 < p < oo. Assume that
H e whim), E e we) in e satisfy:

V XE = —iwufi,
V XI? = iweE.
Then
I? < C H
n IIW1,p(Q,)_ n 1me)
for any/9’: {IBIR—6<R'< M < R} CO.

Remark 1.3.3. In fact since no forcing term is present, Theorem 2.5.2 holds also for

p=1,oo.

Let Q,- 29 and 96 = {.rlR < |:1:| < R+6}, and 02 = {:1:|R—6 < |:c| < R+6}.
Thus in 92, the Helmholtz equation (1.2.1) holds. We have

H . < C H

14

for any 91 C 02 from standard interior elliptic estimates [34]. Let Q’ = Q,- 091.
Then

HW < C H < C H

In order to prove the theorem, we need to estimate ”HUI/pale) by HHHLP(QZ-)'

Claim: anmee) g C||H|| 121% SR)'

If the claim is true, then we have

IIHIIW1,p(Q,1)_ C(IIHHme HIIHIIW ))

The trace theorem and Nirenberg-Gagliardo inequality [2] imply that

< < a 1—0,
IIHIILp(SR) _ CIIHII 1+ ,p _CIIHIIW1,,,(Q,,1)IIHII,,D(Q,1)
W2 ( '1)

< H C H I ,

1
where 0 < C < 2’ 0 < a < 1, n is any positive constant, and C77 depends on the

choice of 77. Thus

H < C H .
Proof of the claim. We prove the result by examining each of the three possible
cases.
Case 1: p = 2. Note that
27r 00
IIHHiQ R22] ammo
0 l: 0 —
m—

00

[:0 mZ—zl agni/[W (6 ¢)R251110d6d¢

Z: le um]2R2

220m:—

by the orthonormality of spherical harmonics. It follows from (1.2.3) that

15

l R+6 h(1)(kr)
um2 l
g Ill/R |———

h]” (kR)

00 l R+6 h(1) hr 2
= Z Z IuPIQRZ/ I—l ( )IQT—dr

[2r2dr

 

 

l R+6 2
= Z Z luml2R2 7“— k27‘2 d7.

 

 

 

k2R2
00 l R 51+al—1——+ +al—1—
= Z Z lainIQRZ/ + 11.272 ”“2321 dr
R z_1___ l__
l=0m=—l 1+“1k2R2+°"+alk2zR2l
00 l
s 6 Z Z MFR?
1: 0m 2 —l
2
C H
_ n was )
It is worthwhile to note that from the computation above
H S H for all r 2 R.
n ”mgr, H 1,26,?)

Case 2: 2 < p 3 00.

We only need to consider the case p = 00, the rest follows from Riesz convexity
theorem. The idea here is that we show first that for small enough 6 > 0, there
is a positive function f (r) on [R, R + 6] such that the maximum of f (r)|H | on Qe

is on 096, i.e., we have the maximum principle for f (r)H on the domain 98. Re-

16

member that this f (r) is needed since H itself can not have the maximum principle.

After this modified maximum principle is established, we only need to show that

H < C H .
” ”LOO(SR+6) _ ll ”LOO(SR)
To construct the function f (r), set V = f (r)H . We have

Av: (1(§(2§»+-} —A§V

2
_-:T—I2{f(r )+2-a£f(r )+f"(r )H+g%—I:-f(r )+§Hf’(r)+;12-A5Hf(7‘)

=a4¥HRa+00%greraflewr+;Hfle>

8V ,
=—k2V+2ff+TV+-12:Vf7
V’f__f—2__ VfI+ +fTNV+2Vf7l
ff)? +—f’—’+§£]V+2;2flv

In order to use maximum principle of elliptic differential equations, we need to find

= —k2V+2f’

=I—k2—u

 

a function f such that

f/2
f

f” 2 fl<

)2_____

k2——+2( f Tf_0

fl

Let —f— — y. Then the equation becomes
2

Now set y = k tan(k(r—R)), (5 < 217%. Then y’ = k2(1+tan2(k(r—R))) = k2+y2,

2.
and y 2 O on [R, R + (5]. Consequently, k2 + y2 — y' — 7y 3 0 on [R, R + 6].

It is easy to verify that

 

l
f : cosk(r — R)
is a solution for y = k tan(k(r — R)). Therefore 1 g f S x/2 on [R, R + (5] for any
6 S ZR. W'hen r e [R, R + (5], the function [V] achieves its maximum value on the
boundary.

17

Next, we will show that
H < C H .
I] ”LOO(SR + 6) _ H “LOO(SR)
From the arguments for p = 2, we know that

CllHlle s CIIHII 190(33)-

“H “L2(R<r<R+26)< L2(SR)

The standard elliptical theory concludes that

H <C H .
H “W2’2(R+6/2<r<R+36/2)_ H llL2(R<T<R+26)

The Sobolev embedding theorem implies that
H < C H .
H “LOO(SR + 5) — H ”W230? + 6/2 < r < R + 36/2)
Combing these estimates yields
H < C H .
[l ”LOO(SR+ 6) _ ll ”LOO(SR)

Case3:1§p<2.

Consider a sequence of smooth function {Hn} on S R1 such that

[[Hn — H|| L1(S:)) —> 0 as n —> 00. Write Hn in spherical har-
monics: Hn(R =2 2— lmh Define THn(R,0,¢) =
l— — Om— —
[20771: —(l———1)(k)hnlYlm(6 Q). There exists g(r,6’,q§) E LOO(Qe) and
= = ‘lhl )(kR)
N
g(r,6, (1‘)) 2219;1(6, ¢>X1k= where X1}: is the characteristic function on 1k,
MAT

U 116: U1] [Tkfl‘k +1] : [R, R + (5], and [[gk(6,¢)[]Loo(QA) S 2, where
k— -1 k— — '
9k = {(r, 6, ¢)[r E 1A.}, such that

18

llHnll = llTHn(R,9,<f>)ll

L1(Qe)
= | / THn(R,0,¢)9(T19,¢)dl/l
fie

L1(Qe)

N
=l 2: / THn(R,6,¢)gk(6,¢)xldeI

N
=I 2 f9 THn<R,0,¢>gk(6.¢)d1/I=I
= 11:

Write gk(6, (1)) in expansion of spherical harmonics:

277 7r 00 1
f0 TH(R,0,¢) Z Z77 17,2}1/ m,(6 ¢)r2 sin0d0d¢dr|
m:

k=1rk l=0 —I
N r 27r 7r
k+1
=| 23/ / H(R,6,¢)

———)—nZ]Ylm (6, (1’))?“2 sinddddodrl

From Case 2, we know that ||Tgk(6,¢ (15)||Loo(Q e) g
t
Cl] 2 Z nEI/lm( (6, ¢)||Loo(S R.) Therefore
l— — Om—— - —l R
I < H T I- < C H .
_ C” ”L1 (2:1 [I QkHLOOm ell 1.]— I] ”L1(SR)

By the density of COO(SR) in L1(SR), we proved the conclusion for p = 1. Again

from the Reisz Convexity Theorem, we get

“H“LPme) S CllHllLP(SR)

forlgpg2. El

19

1.4 Existence and Uniqueness

In this section, we establish the existence and uniqueness of solutions for the lin-

earized SHG problem. By linearization, we mean the following equations:

.7 x50») = _.-..,.Og<w>,

H

V xH(w) = iweOE(w)
and

V XE<2W> = —i2quI-f(2w),

- . ~ . 2
V XH(2W) = z2weOE(2w) + 22w 2k X(j;)cE('w)E}(cw)-
J,

 

This approximation assumes that the electric field of the incident and diffracted
waves at the initial frequency to inside the nonlinear medium acts as a source for
field generation at 2w. In addition, the SHG is assumed to be so weak that its
influence on the field at the initial frequency is negligible.

From the regularity result proved in this paper, we obtain the following well-

posedness result.

Theorem 1.4.1. Let SR 2 {:1:| |:c| = R}, BR 2 {:I:| |:1:| < R} and Q C BR.
Suppose Xe]: E LOO(Q), x2226 = 0 in BR \ Q, supp(m) C Br \ Q for some 1" < R,

and m E LOO(Q). Let

 

60 in BR\Q,
f:
61 in Q.

Then the linearized SHG model problem

V xE<wl = —2'w,LLH(w),

V XH<W> = -iweE(w) + 773.,

20

V XEQw) = —2°2w11H(2w),

v xHQW): i2weE(2w) + 121.229 x- k )Ejl‘” ”Eff” in 83,
j. k

8H ~

0n _=|SR TRH on SR

has a unique solution (H(w),fj(2w)) 6 W111? and (E(w),§(2w)) 6 LP for any

1<p<oo.

Proof. By [1], we know that there is a unique solution H(w) E W11 2(9) and El“) €
L2(a).
From the regularity result Theorem 2.5.2, H(w) E W11p(Q) and 13(0)) 6 Lp(Q).

_. . 2
Hence 9 = z2kaxlngJ<w)El(€w) E Lp/2(Q), for 1 < p < oo.
.7,

It follows that there exists a unique solution H (2w) 6 W11p(Q) and consequently

E9”) 6 [19(9) by a similar argument. El

Remark 1.4.2. If in addition, xii): 6 Cam) and m 6 C062), we have H(2w) E
W1’p(Q) {1011062) and HQ“) 6 LPG?) 000(9) by the following argument.
According to the standard elliptic regularity theory in [34], H E C5 (0’) for some

0<fi<1and||H|| stithQ’CQ.

03(9')
The standard elliptic regularity results indicate the C 1’ a regularity of H away
from a tubular neighborhood of S and near the boundary S R'

For any fixed :130 E S and r > 0, denote
QT = {.17| ICL‘1 — 1(1)] < r, Ircg — 1:9] < r, [$3 — (cgl < r}.

One may choose R such that Q R C 9. Using a transformation, Q R is mapped into
a set containing a neighborhood Q R0 of the origin. Without loss of generality, the
preimage of Q R0 is assumed to contain Q R /2. For simplicity, we shall omit the

primes and set (in the new coordinate system)

21

Q1210 = QRO flag > 0},Q;,0 = QRo flag, < 0}.
Consider a more general model problem in QRO

ox, Wax]- ox, J” Jar]-

+du+f =0. (1.4.1)

Suppose that u E 07(QR0) flH1(QR0). Suppose also that the coefficients aij,

7 J)

part of the operator is elliptic in Q R 0, i.e., there is a constant CO such that

bj E 07(62350) and cj, d, f E C(QEO) have ajump at $3 = 0 b- and the principle

Ease-5,12 cola?

The 01,01 regularity of H near the boundary S can be obtained by the following

theorem. See [8] for a proof.

Theorem 1.4.3. Under the above assumptions, the solution of (1.4.1) satisfies

u e Cl»a(QIi,O/4)

forsomeO<01 <1.

22

CHAPTER 2

Numerical Solution of Nonlinear

Diffraction Problems

2. 1 Introduction

Consider a plane wave with frequency w which is incident on a nonlinear periodic
structure (grating). The nonlinearity of the structure gives rise to the diffracted
waves at both frequency w and 2a). This important phenomenon is known as Sec-
ond Harmonic Generation (SHG) in nonlinear optics. A significant application of
SHG is to obtain coherent beams of light in parts of the spectrum at which lasers
cannot be made and to construct optoelectronic devices based on nonlinear effects in
waveguides and Optical fibers. We refer to [52] for a detailed description of nonlinear
optics.

In practice, however, the SHG optical effects are often too weak to be observed.
Therefore, modeling and enhancement of SHG are of great interest to potential real
applications. Recently, a PDE model was introduced in [48], [49] and [50] to describe
nonlinear SHG in diffractive gratings. It is also pointed out in [49] and [50] that

the SHG can be greatly enhanced in periodic structures. Questions on the existence

23

and uniqueness have been studied in [8]. We also refer to [12] for more recent results
on the optimal design of nonlinear gratings.

To solve the model scattering problem, the first difficulty is to truncate the
domain into a bounded computational domain. In [7] and [11] the authors used finite
element method based on variational formulation in the bounded domain containing
the medium, with periodic condition in (1:1 direction and transparent boundary
condition on the top and bottom boundaries. The derived transparent boundary
condition is represented by a quasi-differential operator and is nonlocal. In practical
computations, the infinite series in the definition of the quasi-differential operator
have to be truncated. Here we apply the perfectly matched layer (PML) technique
to truncate the unbounded domain. PML was first introduced by Berenger in [21]. It
provides a reflectionless interface between the region of interest and the PML layers
at all incident angles. The layers themselves are lossy, so that after a few layers the
wave is significantly attenuated. The main advantage of a perfectly matched layer
as a. boundary condition is that it provides a reflectionless interface for the outgoing
wave at all incident angles. Another advantage is that it preserves the sparse nature
of the FEM matrix, so that the matrix system can be solved easily. We refer to [54]
for a review on PML methods. In practical applications involving the PML method,
there is a judicial compromise between a thin layer, which requires a rapid variation
of the artificial material property, and a thick layer, which requires more grid points
and hence more computer time and more storage (See [25]). In this paper, we use
an a posteriori error estimate to determine the PML parameters. Moreover, the
derived a posteriori error estimate shows exponential decay in terms of the distance
to the computational domain. This property leads to coarse mesh size away from
the computational domain and thus makes the total computational cost insensitive
to the thickness of the PML absorbing layer.

Moreover, since the grating surface is usually piecewise smooth, and across the

24

surface the dielectric coefficient is discontinuous, the solution of the scattering prob-
lem will have singularities which slow down the finite element convergence when
using uniform mesh refinements. The a posteriori estimate adaptively determines
the finite element mesh size which overcomes this difficulty.

We refer the readers to [10],]11], and [47] for a general review of the modeling
and computation of the grating problem. An introduction of PML and adaptive

finite element method applied to linear grating problems may be found in [22].

2.2 Modeling of the Nonlinear Grating Problem

Assume the medium is nonmagnetic (u S no) and no external current or charge is
present in the field. For convenience, the magnetic permeability is assumed to be

unity everywhere. The following time harmonic Maxwell’s equations (time depen-

dence e-iwt) hold:

VxE'zzi‘J-H V-H=0, (2.2.1)
C
V x H = ——D V - D = 0, (2.2.2)
C

where E is the electric field, H is the magnetic field, E is the electric induction, and

c is speed of the light. The constitutive equation is:

_.—-.

[j = 6E + 47rX(2)(:1:,w) : EE,

where e is the dielectric permittivity, w is angular frequency, and X0) is the second

2)

order nonlinear susceptibility tensor of third rank, i.e., X( : RE is a vector whose

3
jth component is Z xfi] EkEl, j = 1, 2, 3. The medium is said to be linear if

k,l=1

13 = 65 or X02) vanishes. In principle, essentially all optical media are nonlinear,

i.e., D is a nonlinear function of E.

25

In this paper, we only consider the 1-D grating problem by assuming that all
fields are constant in the $3 direction. The medium is determined by the dielectric
coefficient e(:r,w) = 6(1‘1,:1:2,w). Assume that the dielectric coefficient is periodic

in 51:1 direction with period L:
e(:r1+ nL,a:2,w) = €(:1:1,:1:2,w), Vzrl,:1:2 E R, n integer.
Assume that the nonlinear medium is contained in the region
Q = {($1,132): 0 < 51:1< L and b2 < x2 < b1}

for some positive constants b1 and b2.
Assume that e is constant away from a region Q, i.e., there exist constants 61

and 62, such that:

6(331,;1:2,w) = 61(w) in 91 ={(:1:1,:1:2):a:2 2 b1},

6($1»$2aw)= E2(a)) in 92 = {($61,952) = 5132 S b2},

and e is piecewise constant in Q with jumps at certain interfaces. Assume further
that (21, {22 are linear media.

The assumption on the piecewise linear medium is technical which is needed
to assure proper regularity of the solution. The main theoretical results (Theorem
2.4.1, Theorem 2.5.2) remain valid in the case of an inhomogeneous medium with
sufficient smooth 6 and interfaces by the regularity result in [8]

The electric field of the incident and diffracted waves at the fundamental fre-
quency c121 inside the nonlinear medium acts as a source for the field generation at
the second harmonic frequency 21121, and it is assumed that the SHG is so weak that
its influence on the field at the fundamental frequency is negligible. This is the well

known undepleted pump approximation in the literature. See [49] and [50]. Under

26

 

 

 

5C2 A 91
b1 F1
a
b2 {‘2
Q2
5’31

 

 

Figure 2.1. Geometry of the grating problem

this assumption, the electric induction [3 may be written as:

final) = 5(I,w1)E($,w1),

fi($,w2) = c(:c,w2)E(:1:,w2) + 47rx(2)(:1:,w2) : E(:1:,w1)E(:r,w1),

where 1122 = 21111.

In the linear case, TE polarization means the electric field is transversal to (3:1,
332) plane. TM polarization means the magnetic field is transversal to ($1, 2:2)
plane. In the nonlinear case, however, the polarization is determined by group sym-
metry properties of X(2)- Here we assume that the electromagnetic fields are TM
polarized at the frequency “’1 and TE polarized at the frequency L122. This polariza-
tion assumption is known to support a large class of nonlinear optical materials, for

example, crystals with cubic symmetry structures. See [12] for detailed information.

Therefore

H($,w1) = H(rr,w1)1“73,

1

E(:1:,w2) = E(:r,w2):r:'3.

27

Define for convenience
Lu’j .
kj(:1:) = 7 e(:1:,wj), 111 Q,
= — 61(wj), j,l = 1,2.

From equation (2.2.2), we deduce that

8H(:1:, 1111) _ 8H(:1:,w1)
8172 i (9171

 

( ,0) = —iw716(x,w1)-(E1(x,w1), E2(:1:,w1), E3(ac,w1)).

So

8H(:1:,w1) _ iwl
W - C €($,w1lE1(l‘7w1),
8H(a:,w1) _ ital
——————a$1 — c c(:1:,w1)E2(2:,w1).

Also from equation (2.2.1),
ital

8E3 6E2 6E1 8E3 8E2 8E1
i a _ ) = (03 O, H):
82:2 823 (9233 (91:1 8:131 8:132 c

 

where for simplicity, we omitted the variables. Hence

ltd—1H—aEz_6E1

 

c _ 81:1 022'
It follows that
, 2 2 2 2
1 0 6H 8H c a H 8 H
V- —VH=V- —————,—-——,O =— + ,
(k? ) (w%e(3$1 6371 )) w¥e( 6113? (92:3 )
0 BBQ 8E1 c2 02H 82H
H:T_(—_— =_T—2_+——2)'
zwl 6x1 @5132 “116 8331 0:132

Therefore we get the equation at frequency col:

1
"1
By a similar derivation, the equation at frequency L122 can be obtained:

47w. 2
(A + k3)E : ——;2—2— Z ng](;r,w2)(E(27,w1))J-(E(13,w1))l
I j7l:132,3
= Z pj,la.x,Hax,H.
3,1 = 1,2

28

- l 167r (2)

(—1)]+ ___—M --

(emu)? 33’
Throughout, we assume that klj > 0, 9?ij > 0, S‘s/£2]- 2 O, 3161(2) > 0,

where pjl =

(311:1(1‘) Z 0.

2.3 Variational Formulation

2041.131

Let u I = e — i5 151:2 be the incoming incident plane wave upon the grating

surface from the top, where 011 = kllsin 6, HI = kllcosl9, and —g < 0 < Z;- is

the angle of the incident. We are interested in the “quasiperiodic” solutions, i.e.,
solutions (H, E), such that u = He_ialxl and v = Ee—i0‘2xl (012 = 1:21 sin 0)
are periodic in 2:1 with period L.

n)

2 n
For each integer n, let 01( = %. Since 21 and v are periodic in 11:1, they have

the following Fourier expansions:

n
u(,:131 11:2) 2: u(n ia( )xl,

nEZ

Tl
v()$17$2 =2 ,1an )fl,

nEZ

L . L .
where u(n)(:rg) = %‘/0 ue_7’a(n)x1d$1 and u(n)(:1:2) = %/0 116—201(71)“de-

Hence we have the expansion for H and E:

H ___ ueialxl = Z u(n 01(705171 + ialxl
n E Z

E : 116202171 2 Z 7J,(n)eior(n):r.1 +ic125131.
n E Z

Define Fj = {($1,222): 0 < 2:1 < L,:1:2 = bj},j=1,2. We wish to reduce
the problem to the bounded domain 0. The radiation condition for the diffraction
problem insists that (H, E) is composed of bounded outgoing plane waves in {21 and

$22, plus the incident wave u I in {21.

29

Since H satisfies the Helmholtz equation AH + [£ng = 0 in {21, we have

2 . n
2 [111+ oil)2 + 5—,,(n)(,2,,6.(a< ) + “1)“ = o,
n E Z $2
2
k21— —(a (n) + (11)2 + $11(n)(22) = O for (152 2 ()1. (2.3.1)
I2

For any integer n, and j,l- —— 1, 2, define fill] that satisfies (fl? [)2- — [€31- (01 (n)+

a )2 and u(finl) > 0. One can easily verify that ,611— — El. Solution of (2.3.1) can

J

then be written as
- n ~ - n

with complex constants U f”) and 0]”) . The radiation condition implies H 101) = 0

in 91 and gives:

(71)
H=u1+ Z 11]" +all$1+lfl11$21$691
nEZ

Similarly, we can deduce the following equations:

2: U2(n)e (n)+ (11).??1 — zfl12$2 :L‘ E Q2,
n E Z

:2: Vln( (71).]. 012).??1-1- 2321232 :1: E 91,
n E Z
n E Z

- n
For any quasiperiodic function at frequency 1.121, f = E: f (n)ei(a( ) + (311)31: 1
n E Z

or at frequency 1122, f = 2 f (n)eZ(Q( ) + (.12).’1:1, define respectively the Dirichlet

n E Z
to Neumann operator, which is introduced in [11],

Tut: Z islnlfmleimmuf’j)“, 0<I1<L,j,l=1,2.

30

The 1-D grating problem can then be formulated as follows:

V-(—1§VH)+H=O inQ,
’“1
.2 _ .
(A+k2)E— Z ijoxjHale in Q,
331: 1,2

8(H—u1)

T—T11(H—u])=0 our],
%—5— T12H— — 0 on F2,
BE

5; —T21E=0 on F1,
@5—5— T22E— — 0 on F2,

where u is the unit outer normal to 652. See [8] for the details.
To find the variational form of this problem, we need to use the above transparent

boundary conditions and introduce the following subspace of H1(Q):

XJ-(Q) = {w E H1(Q): wa = we mjxl is periodic in 2:1 with period L}.
Define Bj: Xj(f2) x Xj(Q) H C:

811.2 w=> f(kwwmczz-i/jkl—51T13991wdx
j2=1

=/ Wow-[13w— 22: / flay-«122.6.
Q Q j=1 Fj

Bu

Note that —a—1 —T11u1 = —2i1‘31u[; the weak formulation of the nonlinear 1-D grat-
1/

ing problem then reads as follows: Given incoming plane wave u I = 62013: 1 _ zfil 332,

find H E X1(Q) and E E X2(Q) such that:

 

B H,1/) = —/ u 'wdr, VLOEX (Q,
1( ) F1k21 I 1 )

8__H_ 8H—
j,l= 12p

31

Assume in the following that the variational problem has a unique solution. Then
the general theory in [6] implies that there exists constants 71 > O, 72 > 0, such

that:
l8j(80,1/))l
sup ———-—

”I” VCpEXj(Q).
o¢weHHm VHVm

2.4 PML Formulation

Now we introduce the PML layers. We surround our computational domain 0 with
two PML layers of thickness 61 and 62 in 01 and 92, respectively. The specially
designed model medium in the PML layers should be chosen such that either the
wave never reaches the outside boundary or the reflected wave is so small that it
essentially does not affect the solution in Q. Let 3(x2) = 81(1‘2) + i32(a:2) be
the model medium property which satisfies 81,82 6 C(R), 31 2 0, 32 _>_ 0, and

s(:z:2) = 1 for b2 3 2:2 S ()1. Introduce the PML regions:

95A“; ={(:z:1,:132)2 0 < 11:1 < L and b1< 332 < b14451},

92”“ = {($1,352) 0 < x1< L and b2 —62 < £2 < b2},

PfML ={(:1:1,:c2): O < 2:1< L and $2 = b1+61},

Pgfl’ll’ ={($1,1‘2)Z O<:1:1< L and 11:2 =b2—62}

and the PML differential operators:

£—~34 lsei34+ a 1 1 ai+ari
1 (91:1 kflli) .2 0:171 8:192 1%(33) s(:1:2)(9:z:2 '2 ’

a a a 1 a 2
= — k. .. . .
011(5(12)011)+ 01:2 3(1‘2) 0:1‘2) + 2(T)S(T2)

 

 

 

 

£2

Let D = {(171,172) :0 < 11:1 < L,b2 — 62 < 172 < b1+ 61}. The PML model can

be formulated as follows:

32

$2“

 

 

 

 

 

b1 + 61 QPML FiDML
1
b1 P1
9
b2 Q£3.le F2
b2 — 62 Pf M L
E£1

 

 

Figure 2.2. Geometry of the PML problem

C H = —g
1 1 in D,
32E = —92
where
£111.] in QfA'IL,
91 =
0 elsewhere,
92 Z

— Z pj lamj H 8171 H elsewhere
with boundary conditions:

A

H(O, 1172) = e_m1LH(L,:I:2) for ()2 — 52 < 1:2 < b1+ 51,

A

E(0,$2) = 8—i02LE(L,fB2) for ()2 — 52 < 1:2 < b1+ (51,

33

H2211 on FfML,
H = 0 on rgML,
E: 0 on prL,
5:0 0.. rgML.

Let G be any open set in D, introduce the subspace of H1(G):

Xj(G) = {w E H1(G) : we, 2 "we—0.7.221 is periodic in 1:1 with period L}

and the sesquilinear form AjG: X j(G) x X j(G) +—+ C as follows:

1 (flaac’MJr 1 16¢aa/3_
S 2 8:131 82:1 k%($) 8(232) 83:2 8:132
8d) 615 1 do 62;

_ __ __ .2 f
Agawm — /G(S($2lax1 an + W) 8352 BIZ k2<x>s(x2>¢w)dx.

 

 

 

Alcoa) = j ( seam/w.

G k%(:z:)

 

Define X§(D) = {w E Xj(D), w = 0 on FfML UI‘gML}. Then the weak formu-
lation of the PML model reads as follows: Find H 6 X1 (D) and 1.3} 6 X20 (D), such
that qul on FfML, H=Oon FgML, and

A1 009,1» = [D 911;de w e x‘fw), (2.4.1)
«@me = [D 921/de W e X§(D). (2.4.2)
1

Let A711]. 2 “5% —(a(n) +01)2|2 and U1j={n:k¥j > (a(n)+al)2},j = 1,2.
And let

[Xi] = min{A7fj Z n E Ulj}, A17: min{ATllj I n ¢ Ulj}'

Introduce the following notations:

b1 + (51 b2
01 Z/b .3(T)d7', 02 = /b 6 3(T)dT.
1 2 _ ‘2

34

— +
2AM 2A11

 

 

A411 : max(
I — ’ R + ’
€201A11_1 6201 A11 _1
2A" 2A+
max( 12 , 12 if 362(011) = 0,
201 a- 20RA+
2|k12|

 

'f 36 (w )> 0
R 1 2 1 v
3202 |C3k12| _1

R I

where a]. and a]. are the real and imaginary parts of a -, respectively. Define

 

C = \/1+(b1— b2)_1.
It is proved in [22] that the problem (2.4.1) has a unique solution H, if (M11 +
A112)C'2 s 71, and the following estimate holds:
|31(H - HAM
uwu 1,16,,

+(

lllH—Hlllm == SUP
0¢¢6Hkm

 

M126
2
kn

 

 

)IIHIIL2(F2)-

Next we prove the existence and uniqueness of (2.4.2) and derive an error

estimate between E and E. We first find an equivalent form of (2.4.2) in do-

main 9. Similar to the previous argument, we write [:3 in the expansion E‘ =

- 72
Z 0(")(12)e2(0( )+a2)$1 and deduce that

n E Z
n [T2 n $2
2 S(7’)d7’ . —ifl / 3(T)dT
E = W 21 b wink 21 > M”) +a2>r1
n E Z
in QffllL’
b b
2332/ 28(T)dT . #11332] 28(T)d7’ ( )
E : (V2006 1‘2 + {/2006 ”1:2 ) 81(0 n + a2):r1
n E Z
1.. 951%.

35

Then the constants V101), V101) , V201) and V201) can be uniquely determined by the

boundary conditions E = 0 on FfML and FéDML:

1
:15 11:2
21331 5(T)dT . . —ifigl 8(T)d7‘
175% f + fink. / 0,
“201) + 12/2“) = 13(")(b2),
b b2
zfln 3(T)dT . —i/3n S(’T)d7'
Vénle 22[132 + Vén)e 22 [$2 _ 0

Thus we have:

 

 

( )
E = (32(x2)0(n)(b )ei(a(n) + 0‘2)“ in QgML,

where

b1 + 51 ()1 + 51
—ifl31/ 8(T)d7' #331] S(T)dT

“n 532 n £132
—2 22/1) 6 S(T)dT 2322/!) 6 S(T)d7'
2 — 2 — e 2 _ 2 ,

- n
For any quasiperiodic function f = Z f(n)ez(a( ) + 0(2):”, define

. . . . .- (n) . ,. .

n E Z
Then
81? . .
‘5; —T2P1AILE= O OI] F1,
053 .

36

Introduce the sesquilinear form 35%le X2(Q) x X2(Q) H C as follows:

2
85’ ML(¢,212)= A)(V¢Vd3—k3($)¢i5)dr- Z) Aug/"Loam.
.7

j = 1
We get the following variational problem: Find 19 E X262), such that:
PM L 6H 6H2) _
B2 19 — 2.4.3
( ,w— 1231,22,]an 0,le < )
j) —

The following lemma establishes the relation of this variational problem to the PML

model problem (2.4.2). The proof is straight-forward from the above constructions.

Lemma 2.4.1. Any solution E of (2.4.2) is a solution of (2.4.3). Conversely, any
solution 29 of (2.4.3) can be uniquely extended to the whole domain D to be a solution

of (2.4.2).

1
Let A21]- : |k2j - — ((101) +ag)2|§ and U2]: 2 {n : kgj > (04") +012)2}, j = 1,2,

then E2? = A2j for n E U2]: and fi2j— — iA2J- for n E U2j. Let
A2]. = min{Agj,n E Ugj}, A2j— — min{A2j,n E Ugj}.

From Lemma 2.2 and Section 5 of [22], we have the following lemma which plays an

important role in the subsequent analysis.

Lemma 2.4.2. For any qb, Lb E X262), the following estimate holds:

 

 

A 7
I [F .(sza5 - T21} ”man s hay-Haw ,qung ,). (244)
J J J
where
25 2M

[”21 = max( 82 IA?1 7 RAEI a

0121-1301 21-1

2A“ 2N
max( 2 [Ag2 , 2 RA? ) if (3552022) = O,

A122: €02 22—1802 22—1

2lk22|

 

. O. , 0
R .. Zf\562(a12) > .
€202 Ignaz! _ 1

37

Lemma 2.4.3. For any 4/) E X2(Q),

llwll 2

L (rj) (2.4.5)

<I' <6
_ ”w”11%(r )._ Hwflfflal)

j

 

with e = \/1+(b1— (2)-1, if

/)(b.’131,])=:Zi/102)(b])e (0(7)) + 02):!)1 on Fj’
l
“4” H. (L 2 ((7‘1+|a)2+a2| 444513124.
H2< J nEZ

Theorem 2.4.1. Let 72 > 0 be the constant in the inf-sup condition, and (4121 +
M22)C'2 < '72. Then the problem (2.4.3) attains a unique solution E. Moreover,

the following estimate holds:

|32(E - E. 4)!
“4421(0)
+ C‘M22IIEH

 

mE—mmg: mp
044eHHm

S (3M21HEII

(2.4.6)

L2(P1) L203)

+ our? — HIIH1m) [IIHIIH1 + 5(2) + HHHH1+ 4(2)||l

for some constant C which depends on the data of the original grating problem and

1

.5).

some constant 6 E (0

Proof. It follows from Lemma 2.4.2 and Lemma 2.4.3 that

IBPML(4 4)I>IB2(44) )I— Z/( T24— TRAIL4>4441
j=1

2 (32(444 — (M21 + M2462 “t“ H1(o)"'*"’”H1(Q)'

From the assumption (21121 + A122)C'2 g 72, it is obvious that the bilinear form
1323]” L

satisfies the inf—sup condition, and hence the problem (2.4.3) has a unique

solution. It remains to prove the estimate (2.4.6).

38

Clearly,

62(4) — E, 4) = 3514142, 4) — 8203‘. 4) + 82(5, 4) — 85“”(13, 4)
2 fr (T21 —TPML)E44141+ [1‘2 (T22—T21‘32A4L)H4441
1
3H aH 6H 8H _
4 z 44/ ___—.444
j,l—1 2 3 Q 871758—331 ij 8ch

By Holder’s inequality,

0H 3H 6H 8H
I Z pjl/Q( 5—4 0—4, “Jamal

 

 

j,l=1,2
(6H8(H— 1H) 0(H—H)8H_
j,l=1,2
aHaH— H (H— H)aH_
3 Z ijzlll 81:] ( >441 .4|+|/‘9((9 ___—WI]
- r2
],l=1,2
‘ 3H 2—2 1
E Z ijlllllH—H” 1M(/(——)4d4)4
j,l=1,2 H (a) 9‘9“]
~ 8H _ l
+ ”H — Hll 141(9)‘ [9(5x—l)24244)4]
< Z l-IIIH—Hu [ (51124,, ]2—,,( 942%);
j,l=1,2

aH _ _
+ ”9(a)“ )21‘9414] 2” (fa 42444) 2"}

1 1
for some 1<p,q < 00, such that — + — =1,
P q

362 Z lPfllHH—HIIHI 9 “W”

21:1 2 ( ) WI’QWQ)

0,1). Take

for some constant 62. By [12], [28], H,H E H1+ 6(0) for some 6 E ( 2

and using Sobolev imbedding theorem, we have:

P=1—_64q=34

39

8H 8H 8H 8H -
Z

I 4-, (—— — ———)444(
j,l=1,2 J Q ij 82:1 ij 811
< (3 H — H H + H - ,
_ II n H1 (9) [n IIH1+ 45(9) ll IIH1+ 5<mlllllwl|H1(m
where 5' depends on the data of the original grating problem. [:1

2.5 Discrete Problem

In this section we introduce the finite element approximation of the PML prob-
lems. Let M h be a regular triangulation of the domain D. Recall that any trian-

gle T E M h is considered closed. We assume that any element T must be com-

 

pletely included in QfML, Q§ML or {-2. We also require that if (0, z) is a node
on the left boundary, then (L, z) is also a node on the right boundary, and vice
versa. Let th(D) E X j(D) be the conforming linear finite element spaces and
gph(D) = th(D)flX3-)(D). Denote by 1h: C(D) +—+ V1h(D) the standard finite
element interpolation operator.
The finite element approximation to the PML problem reads as follows: Find
Hh E V1h(D) and Eh E V2h(D), such that Hh = [huI on FfML, Hh = 0 on
FgML, Eh = 0 on Film/1135'}, = O on FgML, and

AlD(Hh,'i/Jh) = /D glib—hdxfltbh E Vloh(D), (2.5.1)
A2 0(137444494.) = /D 92h12’7hd‘chl/Jh e V2330), (2.5.2)
where
0 in QfML U QQPAIL,

()2 = ~ ~
2h _ Z leanHhaleh elsewhere.

j4l=142

4O

Assume the problem (2.5.1) and (2.5.2) has a unique solution (Hh, Eh) E V1h(D) x

V2}1(D). Let

Then

 

 

£1= V°(A1(1‘)V)+ 31(4).
£2 = v'(A2(113)V)+ 32(4).

411,44, 4) = /D(Al(4)V4V4 — 3444444,

2420(4. 4) = /D(Az(4)V4v4 — 82(4)44)d2

For any T E M h’ denote by hT its diameter. Let Bh denote the set of all sides

that do not lie on I}, j = 1,2. For any 6 E Bh, he stands for its length. For any

T E 1121],, introduce the residuals:

RlT == filfihlr + 91|T (2.5.3)

RQT 1: £2EHIT + 92th (2.5.4)

For any interior side 6 E B h which is the common side of T1 and T2 E 1%,, define

the jump residuals across 6 as:

J16 =(A1VthT1 — AIVHh|T2) . 46, (2.5.5)

J26 = (.42\7H,,|T1 — AQVEhITQ) - 24.. (2.5.6)

41

where V6 is the unit normal vector to e pointing from T2 to T1. If e is the side on
the left boundary and e’ is the corresponding side on the right boundary which is

also a side of some element of T’, we define the jump residuals as:

Jle = Anna—25144.14) - 6’221L52—1(HHIT2)14 (2-5-7)
J1.4 = AlllleialLBZ—1(thT)— 5212442)), (248)
424 = 2121115214th) — 4—222L5:_1(2,,(T,)1, (2.5.9)
12.2 = 4211122022— 2 (E414)— ,: (E4144) (2.5.10)

For any T E Alh, denote by an the local error estimator, which is defined as

follows:
1
an= maX.pj($2)[thlR jmTllL2 (562T hellJ je||L2 )2 l4 i=1424
:cET (e)

where T is the union of all elements having nonempty intersection with T and

l8($2)l€—Rjk($2) 23 E QAfIML’k 21,2, .
Pj($2) = J =14?
1 :r E Q,

with Rjk(2:2) defined as:

 

 

I2
Rj1(a:2) =min(]A.1/b:232 T),dr A1 1./b 81(T)d7'),
1
min(A fl)? 82(7' T)dT, A+ fb2 81(7' T)dT)) if‘352(w-)=0,
R4242 = 2 22 2 22 2
ng‘jQI £52 81(7')d7’ if 362(0).” > 0.
Define also
— I + R
2Aile A1121 2A1L1e_A1101
A! =max( , ),
13 —2A“ a] —2A+ R
1—4 11 1 1.. 1101

42

IOIH

l
- 612 2(1 + 261(A;1+ kj1))

 

 

2 1
0'1 = Cmax( J , ,
J —2A—. 0{ —2A,+. of?
1—e 31 1—e 11
l 1
- 2k; 62 2(1+2<52(A2 + k'2))2
Cmax( J22A2‘ 1’ 222+ I]? ) 1f 4620123) = O,
Cj2: 1—8 j20 1‘8 j202

1
62lmaX(14|k_7°2|)(1 + 252(l‘3kj2I + |kj2|))l2

 

R
for j = 1,2. The following is the a posterior estimate of the magnetic field at

frequency to which is proved in [22]:

Theorem 2.5.1. There exists a constant C > 0, depending only on the minimum
angle of the mesh 111),, such that the following a posterior error estimate is valid:

(311111 61412

 

 

|H—H| g( ”H —u H +(— H H (25.11

H hllm £321 ) h 1 L2(F1) [:22 )H )1 L204,» )
@4413

+ I u —u 2.5.12

( ki1 )H h 1 IHL2(PIIDML) ( )

1
+ C(1+ 011+ C12)( Z "i222-
T 6 Mb.

We next establish the corresponding estimate for the electric field at frequency

2421, which is the main error estimate of this paper.

Theorem 2.5.2. The following estimate holds:

 

HIE-13’ Ill 3((3M )IIE‘ ll +(C‘M W?! 2513)
h 29 21 14142041) 22 hlL2(P2) (
l
+C(1+C11+C12)( Z 713T)?
6611111 A -
+ H + H. H —-
( (221 )(H HHH‘W) II 14IIH1+5(Q))II h, u1||L2(I.1)
0041,,

+(

 

( || +HH H )IIH

43

661413

+ , H + H I u —u
( 1221 )(ll I|H1+2(Q) ll 4|IH1+4(Q))141 [(52,442)

 

~ - l
+CC(1+021+022)())H1 14.; +1th 11.4 )( Z 42442.
H (9) H (9) TeM
'h

In order to prove this result, we first establish some lemmas. For any 7,0 E X 2(9),

we extend it to be in X2(D) denoted by (b as follows:

' ’ 444:2 424014) 0‘2 2

i(a(n)+a2):c1 in QPML 3.21 2
.7 1 3 '

Lemma 2.5.1. ([22]) Let uj be unit outer normal to Dill/IL. Then for any (I),
(,1) E X262), the following identity holds:

PML ,7 _ 6472
F] P] J

In what follows, for the sake of simplicity, whenever no confusion of the notation

is incurred, we shall not distinguish 1b from w in 9;) M L .

Lemma 2.5.2. ( error representation formula) For any (1) E X261), let 21) be extended

to the whole domain D as above, and 1,49}, E V20 (D),

32(E — 13,144) = /D 9242(4) - 4),)42 - 1120(5)“? - film)

+A‘ (T21 —T£‘ML)EthId$1
1

+ A. (T22 — T2I22'IL)EHW$1
2

 

6H,, 6H,, 8H 6H -
+ p- f — ————)4d.4. 2.5.15
jlgl 2 fl 9(5xj 8201 at]- 8:121 ( )

Proof. From the definition,

82(E _ 1231,44?) = 32(E — E,v,’r) + 82(E‘ _ Eh???”

44

2 fr (T21 — TszliAIL)EA"I,l-Jd;r1 +/I‘ (T22 — TQISPPIL)E1Z‘dI1
1 2

BHBH BHBH — .
+ Z: pjl/Q( B—xj 67—23 Basj 6:1: —)z,l'd1:
jl=1,2 l l

+BPML(E— Eh,to)— Z/rj (ng— {J’MLxE— Ehwdxl
j=1

.—_ f1“ (Tm—T21?! L)Ehi/2dat1 +/F (TQQ—TQISML)E,,¢7dx1
1 2

0H 8H 6H 8H PML
+. Z p]l_/§‘2(8:L'j 8:131 —j6:z: —Bz—l-)w wda: +28 (E_ Eh '2)

Also

BPML(E— Em»): Ang— Eh,rb)— i/jr TPML(E— Ehwdxl

]— — 1
2 . 82b
_ $1ng — Eh, 11)) + 2 (E — Myers]
. F V
and from £21/2 —— 0 and Green formula
A (E E ) /(E—E)a—de
29PML h? Pj h Buj 1

Hence

which completes the proof. El

Lemma 2.5.3. ([22]) For any 1b E X2(§2), let {b be extended to D as before. Then

the following estimate holds:

—1RQJ'V <0. ,, ":12 2516
H5 6 ¢||L2(QPML)_ 23||LIIIH1(Q), J 3° (-- )

.7

We are now ready to prove Theorem 2.5.2.

45

Proof. Proof of Theorem 2.5.2 From the definition,
82093 -— Eh, w) = [Um—w — was: — ~420(Eh,2/1 — 1%)
+ / (T21 — TilMLlPh’chl
I‘1
+ fr (T22 — T2IEML)Eh1/3d$1
2

thaE, 8H 8H -
+ p' f(— ‘ -—-—-)t/de
j, €12 fl 9 6333' ‘9‘”! 3333' 0951

 

=III+IV+V+VI+VII.

By integration by parts and using (2.5.3)-(2.5.10),

III+IV= Z (AR2T(2/J-1/)h)d$+ Z éLJ2€(Ib—I/Jh)da§).

Tth eCBT

Using the interpolation estimate in [51] and Lemma 2.5.3, we get

II” + IVI S 0 Z U2T||p§1V¢IIL2
T E Mh

l
5 00+ 021+ 022x 2 n3T>2H¢nH1m>
T E Mh

(T)

It follows from Lemma 2.4.2 and Lemma 2.4.3 that

v VI< CM E C'M E " .
I + |_( 21“ hl|L2(I‘1)+ 22” hHL2(I‘2))”P”H1(Q)

By an argument similar to that used in the proof of Theorem 2.4.1, we conclude

that

 

31%,, 31%,, 8H 6H _
IVIII =| Z p-l f < — —)¢d:r|
3.712172 J Q 6173' 31:1 322]“ 8271

s C(IIHII H1 + a<m+||HhHH1+ 6(9))H'Ph ‘ ”H'm'l‘fl'Hlm)’

where 6' depends on the data of the original grating problem by Theorem 2.4.1.

The proof is now complete. Cl

46

 

 

Air

 

 

 

ZnS

 

 

 

 

 

Figure 2.3. ZnS overcoated binary silver gratings

2.6 Numerical Examples

Our first example is from [12]. See also [22] for the implementation of the algorithm.
This example is concerned with the grating enhancement of the SHG effects for ZnS
overcoated binary silver gratings, see Figure 2.3. The enhancement of the field at 2w
is computed and shown in Figure 2.4 with respect to the associated flat structure.
Here, the period of the grating is L = 0.4 pm, the incident angle is 28.920, and the
wavelength A = 1.06 am. The results were obtained for a thickness t = 0.33 pm of
the coating layer, the fill factor 0.43 with respect to the groove depth d.

Another example comes from [43]. The structure is a subwavelength square
grating of the period L = 0.65/\, the fill factor F = 0.09, and the depth varying
from d = 0.01/\ to d = 1.0M. The refractive index of the material is taken to be
n1 2 3.346 at the fundamental frequency and n2 2 3.539 at the second harmonic
frequency, corresponding to the material properties of GaAs at A1 = 1.907pm and
A2 : 0.954/tm, respectively. The nonlinear coefficient is taken to be 240nm / V. The
incident angle is 300. Figure 2.5 shows the enhancement of the grating structure

with respect to the bulk material.

47

Enhancement

70

Figure 2.4. Groove depth and enhancement

 

10-

 

k—L

i

 

 

8.05

0.1 0.15
Groove Depth

0.2 0.25
pm

48

0.3

0.35

0.4

Enhancement

Figure 2.5. Second harmonic enhancement

 

l

2.5

1.5-

 

 

 

 

0.5 - L

 

 

 

 

'\ JL J
0’ [JL 1 -1.“ 1 _1_ _LA_,

0 0.5 1 1.5 2 2.5 3 3.5
material thinkness (pm)

49

CHAPTER 3

Inverse Medium Scattering in

Breast Cancer Detection

3. 1 Introduction

Breast cancer is a big threat to women’s health. Early detection is the best way
for protection. Mammography is the most effective technology presently available
for breast cancer screening. According to the report from the US. Institute of
Medicine (IOM), limitations of mammography include missing up to 15% of breast
cancers, difficulty in imaging women with dense breasts, and inconclusive results.
The limitations of X-ray mammography provide clear motivation for the develop—
ment of a complementary breast-imaging tool to assist in detection and diagnosis.
Studies show that the dielectric properties of normal breast tissues are significantly
different from those of malignant breast tissues in the microwave frequency range.
Other tissues in the breast, like the glandular tissue and blood vessels, also have
dielectric properties different from the surrounding fatty tissue, but the difference
is not nearly as significant as for the cancerous tissue. Typical dielectric properties

of various tissues in the breast are listed in Table 3.1. Readers are also referred

Table 3.1. Typical dielectric properties of various tissues in the breast [27]

 

 

Media 63 03 67* (6GHZ) (7(6GHZ)
skin 37.00 1.10 34.72 3.89
tumor 54.00 0.70 50.74 4.82
fatty tissue 10.00 0.15 9.80 0.40
average tissue 16.29 0.23 15.66 1.03
fibroglandular tissue 21.57 0.31 21.5 1.7

 

to [19] and [41] for the measurement of dielectric properties at specific radiowave
and microwave frequency bands. It is apparent that there is a wide range of di—
electric properties between the tumor tissue and other normal tissues in the human
breast. Under microwave illumination, tissues with different dielectric properties
will generate different responses. Significant dielectric contrast will lead to a high
imaging contrast.

Microwave imaging for breast cancer detection has gained intense attention.
Ultra-wideband (UWB) confocal microwave imaging (CMI) approach provides qual-
itative high resolution images of backscattered energy distributions of the interior
of the breast. See [30] for a review of different approaches on CMI. See, for ex-
ample, [53], [31], [44], for analysis and numerical examples on different approaches.
The nonlinear inverse scattering approach is able to provide quantitative images of
dielectric properties of objects with high contrasts. An iterative algorithm needs
to be used. In each iteration, an equation describing electric field distribution in
heterogeneous media is solved. Then the dielectric properties are adjusted by min-
imizing the errors between measured and calculated electric fields. An iterative
reconstruction algorithm based on the Levenberg-Marquardt method is presented
in [32]. In [24], the authors use Newton-type reconstruction combined with Mar-
quardt and Tikhonov regularizations to update an initial dielectric property distri-
bution iteratively in order to minimize the squared difference between computed and

measured data. Recently the group at Rensselaer has designed and built electrical

51

impedance tomography process that applies currents through electrodes attached to
the surface of the body and measures the resulting voltages. This system uses the
electrical measurements to reconstruct and display approximate pictures of the elec-
tric conductivity and permittivity inside the body. The mammography geometry is
modeled as a rectangular box with electrode arrays on the top and bottom planes.
The reconstruction algorithm is based on linearizing the conductivity about a con-
stant value. In the case that the target conductivity is extremely high compared
to that of the background solution, there is a large discrepancy between the true
and the reconstructed conductivity values. A review of EIT techniques for breast
cancer detection can be found in [56]. See also [23] for numerical implementation.
By considering the tumor tissue as small inhomogeneity in the surrounding normal
tissues, an algorithm using small volume asymptotics has been used to reconstruct
conductivity distributions. In [4] and [5], the authors give theoretical derivation
of the asymptotic formula. In [3] some numerical examples are given. Since the
approach is perturbative, the significant contrast makes the reconstruction more
challenging.

In this paper, we formulate the problem as an inverse scattering problem, which
is to determine the dielectric property of the tissues from the measurements of
electromagnetic field on the breast surface, given the incident field. Our approach
follows the general idea of [20] and employ the recursive linearization algorithm
from [16] and [15].

In two dimensional cases, the electromagnetic intensity satisfies the Helmholtz

equation:
Au+k2(.r)(1+q(k0,:r))u=0, (3.1.1)

where u is the total field; k(:t) = k0t/e(k0,a:) is the wavenumber; k0 is the

wavenumber in vacuum; q(k0, at) is the scatterer which has a compact support and

52

6(k0, 3:) = 1 + q(k0,:r) is the dielectric permittivity in dispersive medium, which is

assumed to satisfy the Debye model [35]

_o k ,1:
“£70,110: €7~(k0,:1:)— l—(LE—l
k0 —0
#0
20—60052?)—

 

 

(3.1.2)

600(513') )‘1‘6 I
1(+ikO‘/;—% zko‘fé—g

in microwave range frequencies, where 67- is the relative permittivity; o is the con-

ductivity; 63 = lim e; 600 = lim 6; as = lim 0; 60 is the permittivity in
0 —> 0 [co —-> 00 k0 -—> 0

vacuum. In the following, we assume that the material 18 nonmagnetic, i.e., :“0 = 1.

The scatterer is illuminated by a one-parameter family of plane waves
ui = eiko ' 5. (3.1.3)
Evidently, such incident waves satisfy the homogeneous equation
Aui + k311i 2 0.
The total electric field 11 consists of the incident field iii and the scattered field as:
u = ui + as.
It follows from the equations (3.1.1) and (3.1.3) that the scattered field satisfies
Ans + k2(:1:)(1+ (1)118 = (—k2(:1:)(1 + q) + 13m (3.1.4)

In free space, the scattered field is required to satisfy the following Sommerfeld

radiation condition

011

r l-lanO\/:(—é—7: —ik0u 8) =0, r=|:1:|,

uniformly along all directions—— .In practice, it is convenient to reduce the problem

to a bounded domain. For the sake of simplicity, we employ the first order absorbing

53

boundary condition [36] on the surface of the breast:
flu: — ikous = 0. (31-5)

Given the incident field iii, the direct problem is to determine the scattered
field as for the known scatterer q(k0,:z:). Using the Lax-Milgram lemma and the
Hedholm alternative, the direct problem is shown in [16] to have a unique solution
for all k0 > 0. An energy estimate for the scattered field is given in this paper,
which provides a criterion for the weak scattering. Furthermore, properties on the
continuity and the Fréchet differentiability of the nonlinear scattering map are ex-
amined. For the regularity analysis of the scattering map in an open domain, the
reader is referred to [9], [39] and [26]. The inverse medium scattering problem is
to determine the scatterer q(k0,;1:) from the measurements on the surface of the
breast, aslpb, given the incident field ui. Two major difficulties for solving the
inverse problem by optimization methods are the ill-posedness and the presence of
many local minima. In this paper we developed a continuation method based on
the approach introduced in [16]. The algorithm requires multi-frequency scattering
data. Using an initial guess from the Born approximation, each update is obtained
via recursive linearization on the wavenumber k0 by solving one forward problem
and one adjoint problem of the Helmholtz equations.

In addition to the ill-posedness and nonlinearity of the inverse scattering prob-
lem, one major difficulty lies in the multiple scales of the problem. The tumor is
comparably small in the computational domain, which makes the computation chal-
lenging. The strategy is to map the boundary data to the artificial boundary of a
fairly small domain that encloses the tumor. The idea of mapping follows from [55].
The problem may be then reduced to a smaller domain which can be solved by finite
element method. Suitable boundary conditions and jump conditions must then be

added on the boundary of the smaller domain and the surface of the breast. Nystro'm

Figure 3.1. Geometry of the inverse scattering problem

method is used on the integral equation in the annulus region between the surface
and the smaller domain. Another difficulty is due to the dispersive nature of the

human body. We employed the simplified Debye model:

.0803)

etkow) = estx) _. 1.0 «26’

which is the approximation when the frequency is below the range of the Debye

 

model. The Debye model (3.1.2) is suitable for dielectrics exhibiting resonance
effects at the frequency 1011 ~ 1012Hz. The frequency used in our experiments
range from 109 ~ 1010Hz. With this model, the reconstruction of q can be done
separately for the real and imaginary part.

The plan of this paper is as follows. The analysis of the variational problem
for direct scattering is presented in Section 3.2. The Fréchet differentiability of the
scattering map is also given. In Section 3.3, an initial guess of the reconstruction
from the Born approximation is derived in the case of weak scattering. Section
3.4 is devoted to numerical study of a regularized iterative linearization algorithm.

Numerical examples are presented in Section 3.5.

55

 

Figure 3.2. Dielectric properties at frequencies described by Debye model [42]

 

 

 

 

 

"

 

 

3.2 Analysis of the Scattering Map

In this section, the direct scattering problem is studied to provide some criterion
for the weak scattering, which plays an important role in the inversion method.
The F réchet differentiability of the scattering map for the problem (3.1.4), (3.1.5) is
examined.

To state our boundary value problem, we introduce the bilinear form a :

H1(Qb) >< H1(Qb) .3 cc

3(3, 21) = (W, W) — 12((1 + (1115.11) 41113.11),

and the linear functional on H1(Qb)
, __ ”2 ,2 i 1
b0?) — ((1. (1+ (1) - 1.0)11 .112)-

Here, we have used the standard inner products

(MT/93PM and <¢,1>>=/Fb¢-Eds.

where the overline denotes the complex conjugate.

56

Then, we have the weak form of the boundary value problem (3.1.4) and (3.1.5):
find us E H1(Qb) such that

u(u3,g) = 1(5), vs 6 H1(nb). (3.2.1)

Throughout the paper, the constant C stands for a positive generic constant
whose value may change step by step, but should always be clear from the contexts.

For a given scatterer q and an incident field ui, we define the map S (q,ui) by
us = S(q,ui), where us is the solution of the problem (3.1.4) and (3.1.5) or the
variational problem (3.2.1). It is easily seen that the map S (q, uz) is linear with
respect to ui but is nonlinear with respect to q. Hence, we may denote S (q, uz) by
S (q)uz

Concerning the map S (q), a continuity result for the map S (q) is presented in

Lemma 3.2.3.

Lemma 3.2.1. Given the scatterer q E LOO(Qb), the direct scattering prob-

lem (3.1.4) and (3.1.5) has at most one solution.
Please see [16] for the proof.

Lemma 3.2.2. If the wavenumber k0 is suficiently small, the variational prob-
lem (3.2.1) admits a unique weak solution in H1(Qb) and S (q) is a bounded linear
map from L2(Qb) to H1(Qb). Furthermore, there is a constant C dependent of 9b,

such that

k2—1172 k2

2' 1 2'
||S(q)u ”(116211) _ T + 130 + EBWQHLOOWbVHU ||L2mb>~ (3-2-2)

Proof. Decompose the bilinear form a into a = a1 + k2a2, where

Quasi) = (Vu8,V€) - ik0<u3,€>,

agate = —<<1 + auto.

57

We conclude that al is coercive from

,s s us2 u
Ialtu,u11_>.0111v ”21:11.1”0” 3211 1/21 $1

2 s 2
2 Cko(|qu3|l +1111 11 1
21111.1 III/2111.1
2
2 0101121511 .
H 1(91))

where the last inequality may be obtained by applying standard elliptic esti-
mates [34]. Next, we prove the compactness of a2. Define an operator .A : L2(Qb) —>
H1015) by

6116411316) = agate, vg e 11119,),
which gives
(17.4113, vg) — 110(211132) = —((1 + q)u3,g), vg e H1(Qb).
Using the Lax—Milgram Lemma, it follows that

113411511011” S||L2 (32-3)

1111111,< 1 - 11—0 1.2111151
where the constant C is independent of k0. Thus A is bounded from L2(Qb) to
H1(Qb) and H1(Qb) is compactly imbedded into L2(Qb). Hence A : L2(Qb) -—>
L2(Qb) is a compact Operator.

Define a function 45 E L2(Qb) by requiring (15 E H1(Qb) and satisfying

01015.6): 11151. vs 6 H1191).

It follows from the Lax—Milgrarn Lemma again that

“P - k8| 11:12

1911 1 :01 1.0 + [011+ 1$11q11Loomw111u11Lg 1.21,) (32.41

H (12b)

Using the operator A, we can see that the problem (3.2.1) is equivalent to find

 

us E L2(Qb) such that

(I + 12/0118 = a. (3.2.5)

58

When the wavenumber k0 is small enough, the operator T + k2A has a uniformly

bounded inverse. We then have the estimate

3 011111,; (32.61

llu ”12(le (9b),

where the constant C is independent of k0. Rearranging (3.2.5), we have us =
d) — kgAus, so as E H1(Qb) and, by the estimate (3.2.3) for the operator A, we

have

112131 ӣ12131 2 .
k0 L (9b)

<
-11111H19 )+C

( b

The proof is complete by combining the estimates (3.2.6) and (3.2.4) and observing

H1 (91))

that as = S(q)uz. E]

For a general wavenumber k0 > 0, from the equation (3.2.5), the existence follows
from the Fredholm alternative and the uniqueness result. However, the constant C

in the estimate (3.2.2) depends on the wavenumber.

Remark 3.2.1. It follows from the explicit form of the incident field (3.1.3) and the

estimate (3.2.2) that

llusll

Hlmb) 5 1011/2101 + 0211q11Lm(Qb)1,

where Q is the compact support of the scatterer q and the constant C1, C2 depends

on 130, 90'
Lemma 3.2.3. Assume that q1,qQ E LOO(Qb). Then
s 2'—s Vi <0 — 1' , 3.2.7
h th ' t tC d d kt ,Q , d .
u1 ere e cons an epen s on 0 b an [[qQIILoo(Qb)
Proof. Let vii = S(q1)ui and 11.23 = S(q2)ui. It follows that for j = 1, 2

All; + kt2(1+ qflui = (—k2(1 + (1]) + haul.

59

By setting to = u‘i — ug, we have
Aw + k2(1+ q1)'w = —k2(q1 — q2)(ui + 113).

The function 111 also satisfies the boundary condition (3.1.5).

We repeat the procedure in the proof of Lemma 3.2.2 to obtain

.. (C _ i S .
11w11H1(9b)_ 1111 quILoo(Qb)11u+u211L2mb)

Using Lemma 3.2.2 again for u; yields

us < C ui ,

which gives

115311116 — S<q21ui11H1(D) .<. 01111 — an m 19,1231 L2 (911’

where the constant C depends on 9b, kg, and [[qQIILoo (9b). [:1

Let 7 be the restriction (trace) operator to the boundary I‘b. By the trace
theorem, 7 is a bounded linear operator from H1(Qb) onto H1/2(I‘b). We can now
define the scattering map M (q) = ’73 (q)

Next, consider the Fréchet differentiability of the scattering map. Recall the map
S (q) is nonlinear with respect to q. Formally, by using the first order perturbation
theory, we obtain the linearized scattering problem of (3.1.4), (3.1.5) with respect

to a reference scatterer q,

Av + k2(1 + q)v = —k26q(ui + us), (3.2.8)
3% _ 11.01 = 0, (32.9)

where us 2 S (q)v.i.

Define the formal linearzation T(q) of the map S (q) by v = T((1)(6q,ui), where
v is the solution of the problem (3.2.8), (3.2.9). The following is a boundedness
result for the map T(q). A proof may be given by following step by step the proofs

of Lemma 3.2.2. Hence we omit it here.

60

Lemma 3.2.4. Assume that q,6q E LOO(Qb) and ui is the incident field. Then
v = T(q)(6q,ui) E H1(Qb) with the estimate

T 6 , i < C 6 i , 3.2.10
11 (q11 q u ”1111911 _ 11 qIILoo(Qb)11u 11 L2mb) 1 1
where the constant C depends on 160,91), and [IqllLoo(Qb).
The next lemma is concerned with the continuity property of the map.

Lemma 3.2.5. For any q1,q2 E Loomb) and an incident field ui, the following

estimate holds

T 6 i — T 6 i < C — -
- 6 i , 3.2.11
11 q11 Loombfllu 11 11211211 1 1
where the constant C depends on 190,121), and “qQHLOO(Qb)'
Proof. Let vj = T(qj)(6q,ui), for j = 1,2. It is easy to see that
Am - 2121+ 1911 + (1,111. v21:
- k26qtuff - 113) — 162011 - Q2102,

where “g S(qJ)u.

Similar to the proof of Lemma 3.2.2, we get

1 — 1 <C 6 3— S
1111 1215,1me (1 q11Loo(Qb)11u1 251,916,!»
+ — l .
”‘11 (12[[L00(Qb) “P2HH1(Qb))
From Lemma 3.2.2, we obtain
11211 v21 H1 (011 _ 1111 (1211 Looms“ 111Loo(9b)11u “1.211111

which completes the proof. Cl

61

The following result concerns the differentiability property of S (q)

Lemma 3.2.6. Assume that q, 6g E LOO(Qb). Then there is a constant C dependent

of k, Qb, and [IqHLOO(Qb)’ for which the following estimate holds

“S(q + 6(1)'ui- S(qiui - T(q)(6q1ui)llH1 S Cll‘SQIlioo

(9b) ”u2 “12(be
(3.2.12)

(91))

Proof. By setting it; = S(q)ui, u; = S(q + bq)ui, and v = T(q)(6q,ui), we have

Aus + k2(1+ q)u‘i = (—k2(1 + q) + (C8)“,
A115 + k2(1+ q + 6q)u§ = (—k2(1 + q + 5(1) + k811i,

Av + k2(1 + q)v = —lc26qu‘i — k26qui.

In addition, u'i, ug, and v satisfy the boundary condition (3.1.5).

Denote U = ug — u? — v. Then
AU + 162(1 + q)U = —326q(u§ — of).
Similar arguments as in the proof of Lemma 3.2.2 give
U <C6 3—3 .
From Lemma 3.2.2, we obtain further that

3 0111111336

“UHH1(Qb) (9b) [In2 ”L2(Qb).

Finally, by combining the above lemmas, we arrive at

Theorem 3.2.2.. The scattering map 111(q) is Fréchet diflerentiable with respect to

q and its Fréchet derivative is

D111(q) = ')'T(q). (3.2.13)

62

Q0

Figure 3.3. Geometry of the inverse scattering problem

3.3 Inverse Medium Scattering

In this section, a regularized recursive linearization method for solving the inverse
scattering problem of the Helmholtz equation in two dimensions is proposed. The
algorithm requires multi-frequency Dirichlet and Neumann scattering data, and the
recursive linearization is obtained by a continuation method on the wavenumber
190. It first solves a linear equation (Born approximation) at the lowest 1130, which
gives the initial guess of q(k0,a:). Updates are subsequently obtained by using
a sequence of increasing wavenumbers. For each iteration, one forward and one
adjoint equation are solved. Since in this specific problem, the tumor is very small
compared to the breast, finite element method is time-consuming. Our strategy is
to map the boundary data to the artificial boundary of a fairly small domain that
encloses the tumor.

Let 9b be the circle that contains the breast; let P b = 69b be the surface of the
breast; let 90 = R2/Qb; let F be the artificial boundary that contains the tumor;

let 9 be the domain enclosed by F; let {lb 2 (lb/Q. We assume

63

60 in 90,

~ 160 in 90,
5 = Eb in Qb, k =
kb in “0‘
€b(1 + q) in Q.
The problem can be modeled as follows:

Aa+kfl1+wu=0 inQb man
Au0+kfid3=0 hlfiw (332)
u = u0 = us + ui on I‘b, (3.3.3)
(9u aus 311i
— = —— —— 1" . .4
an €b<8n + an) on b’ (33 )
Bus .
E—ik0u3=0 on I‘b,

where (3.3.3) and (3.3.4) are the jump conditions on the surface of the breast.

3.3.1 Born Approximation

Define a test function ii. = eikbx ' d, cf: (cos (9, sin 6), 0 E [0, 271']. Hence it satisfies:
A3+fia=o mob @3m

Multiplying the equation (3.3.1) by a, and integrating over Qb on both sides, we

have

A 2 -
uAudLI: + k f 1 + q uudr = 0.
Ag b§%( )

Integration by parts yields

- 2 . ' "
Auuda: + f k 1 + q uudr + /
‘/Qb Qb b ( ) Pb an an

We have by noting (3.3.5) and the jump conditions (3.3.3) and (3.3.4) that
a” as 65 oi
f9 kgqiri‘tdr = [I‘bmsgg — ebirFun—hls + Ab('u25}é — ebuéfl—Ms,

64

where we take into account that q has compact support in 9. Using the special form

of the incident wave and the test function, we then get

2 «11.12.61 . ~ 2112-11” .9 11220113
1;, , b 1d = .1; -d b 1 _ b l—d
[9 561116 :1: [Pbu bn le 11 e an) s

+/I‘ eik017'd2 + ikbx ' d1(ikbn - d1 — ebikon - d3)ds.
b

(3.3.6)
From Lemma 3.2.2 and Remark 3.2.1, for a small wavenumber, the scattered field
is weak and the inverse scattering problem becomes essentially linear. Dropping the
nonlinear term of (3.3.6), we obtain the linearized integral equation
. ~ . ~ . ~ . ~ 3
/ kgqo($)62k0$ - d2 + ikba: - dlda: =/ (ikbn . deikbrz: - dlus _ eikbx - d1%‘_ms
Q Pb n
+ fr eik0$ ' d2 +2"5b”:"11(z'1.~,,n.. d] — ebikon - (13113, (3.3.7)
b
which is the Born approximation.

Since the scatterer q0(k0, x) has a compact support, we use the notation

where (10(6) is the Fourier-Laplace transform of q0(:1:) with g = (kodl + kbdé), due

to the presence of the evanescent waves. Choose
dj =(C086j,81n6j), ] 21,2,

where Hj are spherical angles. It is obvious that the domain [0,27r] of Oj, j = 1, 2,
corresponds to the ball {5 E R2 : lg] 3 k0 + [kb]}. Thus, the Fourier modes of q0(€)
may be beyond the disk with radius 2k. Please refer to [17] for detailed analysis.
The scattering data with the higher wavenumber must be used in order to recover
more modes of the true scatterer.

The integral equation (337) can be written as the operator form

A060, 0; :r.)q(:1:) = f(lcO, 9). (3.3.8)

65

It is implemented by using the method of least squares with Tikhonov regulariza-

tion [29]
(10080193) = (A =1 A + 60-11,

where a is a small positive number, A>1< is the adjoint operator of A. ‘10 is used as

the starting point of the following recursive linearization algorithm.

3.3.2 Recursive Linearization

As discussed in the previous section, when the wavenumber is small, the Born ap-
proximation allows a reconstruction of those Fourier modes less than or equal to
k0 + |ka for the function q(:1:). We now describe a procedure that recursively de-
termines ko at k0 = kj for j = 1, 2, with the increasing wavenumbers. Suppose
now that the scatterer q]; has been recovered at some wavenumber It, and that the
wavenumber k is slightly larger than 12‘. We wish to determine qk, or equivalently,

to determine the perturbation
6q = qk - (1)}.

For the reconstructed scatterer qk, we solve at the wavenumber k the forward

scattering problem

Ait+kg(l+qi~c(k,a:))fl=0 in ob,
A30 + 33110 = 0 in no,
112'1102118+'ui oan,
05 338 011i

_ = _ _ p
(9n PM 871 + (9n) on b’
011-9 -- -3

For the scatterer qk, we have

Au + kg(1+ qk(k,:1:))u = 0 in Qb,

66

AuO + kguo = 0 in 901

u=u0=u8+ul oan,

au. 3113 3111'

_ = _ __ p
on €b( on + an) O“ b’
, S
95%— — 611063 = 0 on rb. (3.3.10)

Subtracting (3.3.9) from (3.3.10) and omitting the second-order smallness in 6q and

in 6a 2 u — a, we obtain

 

 

A611 + 1%(1 + q]; (k, 2))611 = 4636115 in ob,
A6110 + 336110 = 0 in no,
6u = 6110 = 6us on I‘b,
@ — 6 661‘s on I‘
an ' b on b’
S
as: — ikodus = 0 on Pb. (3.3.11)

Z

For the scatterer qk and the incident wave u , we define the map S (qk, uz) by

S(qk,ui) = uO,

0

where u is the total field data corresponding to the incident wave ui. Let 'y be the

trace operator to the boundary Pb. Define the scattering map

111611-1112) = 130113.111)-

For simplicity, denote [M(qk, at) by M (qk). By the definition of the trace Operator,

we have

M(qk) = 1.01%.

Let D111 (qk) be the Fréchet derivative of Aflqk) and denote the residual operator
by

11(1),?) = 601% — 1301135

67

It follows from Theorem 3.2.13 that
DA/I(ql~é)6q = R(q]~€). (3.3.12)

In order to reduce the computation cost and instability, we consider the Landweber

iteration of (3.3.12), which has the form

6q = 61911155131521, (3.3.131

where [3 is a relaxation parameter and DM*(q]~c) is the adjoint operator of DM (91;)
In order to compute the correction 6g, we need some efficient way to compute

DM*(q}~€)R(ql;), which is given by the following theorem.

Theorem 3.3.1. Given residual H(qh)’ there exists a function gt) such that the

adjoint Fréchet derivative DM*(q]~€) satisfies

 

 

=1- : — 1
[DM ((1,2)Rj(qié)l($) — “(131452? (33-14)
where a. is the solution of (3.3.9).
Proof. Let a be the solution of (3.3.9). Consider the following problem
A611 + 163(1 + qEUc, 22))611 = —k§6qa in 11b, (3.3.15)
A6110 + 336110 = 0 in 90,
(5v 2 6u0 = 6u3 on Pb,
86u _ 86115 n F
3n _ 6b 8n 0 b’
s
83:: — ikoéus = 0 on Pb.
and the adjoint problem
Aw + 1613(1 + (jicUc, 113))1b = 0 in Qb,
21110 + 331110 = 0 in 110,
1’10 _
95n— + 610610 = ( 0 — 210)}:g. on rb. (3.3.16)

68

Since the existence and uniqueness of the weak solution for the ad joint problem may
be established by following the same proof of Lemma 3.2.2, we omit the proof here.
Multiplying the equation (3.3.15) with the complex conjugate of 1b and integrat-

ing over 9b on both sides, we obtain

_ 2 T 2 ~ _.
Adar/1dr +/ k 1+ q~ (Suwda: 2/ —k (Squibdr.
job 11,, b( k) 11,, b

Integration by parts yields

—(95u 31b 2/ ~—
——6 s=—k 5 . d .
/b(wan ua—n)dS b Qb quw 33

It follows from (3.3.12) and the boundary conditions of (3.3.15) and the adjoint

problem that

f (110 — 30)k§eb6uds = kg / 61153112,
Pb 11

DM q~ 6 Rq~ e ds=/ bquibdr.
/Pb (pg (k1,, Q

We know from the adjoint operator D111 * (th that

 

(5 e D111 ~ R q~ ds=/ 6qfl.1,bd:c.
[0be (qk11k1 9

Since it holds for any bq and since q has compact support in {2, we have

 

1 ~ _
DM*(q;.)R(<I,~C) = 2311112.

Taking the complex conjugate of the above equation and letting (b = 1.7) yields the

result. CI

Using this theorem, we can rewrite (3.3.13) as

53 (3.3.17)

c—‘lIIrb

69

So for each incident wave and each wavenumber .190, we have to solve one forward
problem (3.3.9) along with one adjoint problem (3.3.16). Since the adjoint problem
has a similar variational form as the forward problem, essentially, we need to com-

pute two forward problems at each sweep. Once 6q is determined, q]; is updated by

q;J + (5g.

3.4 Implementation

In this section, we discuss the numerical solution of the forward scattering problem
and the computational issues of the recursive linearization algorithm.

The scattering data are obtained by numerical solution of the forward scattering
problem. To implement the algorithm numerically, we employ Nystrbm’s method
in the annulus region Qb and add some suitable boundary conditions on F and
Pb. Readers are referred to [40] for a detailed description Of Nystréim’s method.
See also [26] for the implementation of Nystrém’s method on integral equations
generated by Helmholtz equation. Based on Kirsch and Monk’s idea in [37] and [38],
the exterior problem is solved by integral equation with radiation condition.

Define the space

_ _ . Bus .
W(lR\Qb) 2 {us E H110C(R\Qb)| r 11)me (Hr— — ikous) = 0, r = [23]}.

Define the Operators

1.:

‘1 (1‘) x 11—2011,) _. H1(Qb),

1 _
Ge : H_§(I‘b) _. 111mm,)

1 1
by the following boundary problems. Given ’\I‘ E H _§(I‘ ) and Arb E H _§(I‘b),

70

define C2211 = w where w E H162) is the weak solution of

Aw + hg(1+ qié(k,ac))w = O in Q,
:11)?” + ikow— — ’\I‘ on F. (3.4.1)

Define Cm()\ ,/\ ) as the weak solution of
r r),

Aw + kgw = 0 in (lb,

1

6b —%_n + ikOw— — )‘I‘ on F,

1 8w

a8? + ikow - Arb on Fb' (3.4.2)

Similarly define CeAFb = w as the weak solution Of

Aw + kgw = 0 in R\Qb,

6w ,

5— + ikOw = APb on I‘b,

r _leOO \/1_“ (— — ikow)— — 0, r = [51:]. (3.4.3)

To ensure continuity of solution of the forward problem across F and I‘ b1 it suffices

1 1
to choose )‘I‘ E H_§(I‘) and AFb E H_§(I‘b) such that

(911i , 2' aui . i
GiAI‘ + Cz-(g— + ikOu ) = CmAP + Gm(—a; + ikOu ) on F, (3.4.4)
8_u i .

The function ’\I‘ and Arb are approximated by trigonometric polynomials of or-

der N. Represent F by (C(t) = (rpc0361,rrsin61),0 g 01 S 277, Pb by

N — 1
. . ' 9
a:(t) = (er cos92,r1~b s1n02),0 S 62 S 277. Write )‘F = Z bnem 1 and
n = —N
N — 1 .
Arb = 2 6711217192. Thus for (3.4.4), 2N + 1 finite element problems need to
n = —N

be solved on the left hand side, and 2N + 1 Nystrbm problems need to be solved

71

on the right hand side. Similarly for (3.4.5), 2N + 1 NystrOm problems need to
be solved on the left hand side, and since Pb is a circle, we can compute CeAFb

explicitly as a finite linear combination of Hankel functions:

1 ane 2 l
(ta/‘11. = a Z (11 . (11 HA >031
n = _N (HTL )'(k07‘1"b) + ZHn (horrb)

 

As for the adjoint problem, the continuity conditions are:

Opp + G,((—“u — (1)113) = Gmxp — Ema—11 — mg) on r,

Gmpr + Gm((u — 61133) = 061111) — Ce((_u — mg) on rb,

where the data on the artificial boundary F can be obtained from the data on Pb

by a least square method introduced in [55].

3.5 Numerical Experiments

In the following, to illustrate the performance of the algorithm, two numerical ex-

amples are presented for reconstructing the scatterer of the Helmholtz equation in

 

two dimensions. Assume the dielectric constants Etumor(k01 as) = 54 —i2'6,:)714 and
EnormalUCO’I) = 16.29 — iw (see [27]).
Example 1. Define 0
1112
(M — peep-00025? 11:1 < 0.0025 in 11,
(1050133) = 6normal
0 elsewhere in Q.

See Figure 3.4 and 3.5 for the surface plot of the real and imaginary part of scatterer
function in the domain [1‘] < 0.003. Figure 3.8 and 3.9 are the final reconstructions
using the wavenumber k0 = 7.1, which has relative error 6%. Figure 3.6 and 3.7

shows the result of the Born approximation.

72

Figure 3.4. Real part of smooth scatterer function

 

 

x0.01 - x0.01

Example 2. [31] Assume the diameter of the breast is 10cm and a. 6—mm—diameter
tumor is located in the center.

Ctumor _ 1 $2 + 332 g 0.0032 in 91

q(k0, :r) = 6normal 1

0 elsewhere in 9.
See Figure 3.10 and 3.11 for the surface plot of the real and imaginary part Of
scatterer function in the domain 9 = {as : 35% + 33% 3 0.0062}. Figure 3.14 and 3.15
are the final reconstructions using the wavenumber k0 = 6.1, which has relative error
26.66%. Figure 3.12 and 3.13 shows the result of the Born approximation. It is easily

seen that this scatterer is difficult to reconstruct because of the discontinuity.

73

Figure 3.5. Imaginary part of smooth scatterer function

 

 

X001 '0-4 ‘04

Figure 3.6. Born Approximation of the real part of smooth scatterer function

 

 

x0.01 ’04 ‘0-4

x0.01

74

Figure 3.7. Born Approximation of the imaginary part of smooth scatterer function

 

 

~04 —0.4
x0.01 x0.01

Figure 3.8. Final construction of the real part Of smooth scatterer function

 

 

x0.01 -0.4 —o_4

x0.01

75

Figure 3.9. Final construction of the imaginary part of smooth scatterer function

 

 

"‘14 -0.4
x0.01 110.01

Figure 3.10. Real part of piecewise scatterer function

 

 

x001 '1 ‘0-3

x0.01

76

Figure 3.11. Imaginary part of piecewise scatterer function

3\

 

 

‘1 -0.8
10.01 x001

Figure 3.12. Born Approximation of the real part of piecewise scatterer function

 

 

x001 ’1 '03

x0.01

77

Figure 3.13. Born Approximation of the imaginary part of piecewise scatterer func-
tion

3\

2.5\

 

 

x0.01 '1 ‘03

x001

Figure 3.14. Final construction of the real part of piecewise scatterer function

 

 

‘1 -0.B
x0.01 x001

Figure 3.15. Final construction of the imaginary part Of piecewise scatterer function

 

 

x0.01 _1 '08

x0.01

79

BIBLIOGRAPHY

80

BIBLIOGRAPHY

[1] T. Abbound and J. C. Nédélec, Electromagnetic waves in an inhomogeneous
medium, J. Math. Anal. and App., 164(1992), 40—58.

[2] R. A. Adams, Sobolev Spaces, Pure and Applied Mathematics Series, Academic
Press, New York, San Francisco, London, 1975.

[3] H. Ammari, E. Bonnetier, Y. Capdeboscq, M. Tanter, and M. Fink, Electrical
impedance tomography by elastic deformation, preprints.

[4] H. Ammari, R. Griesmaier and M. Hanke, Identification of small inhomo-
geneities, Mathematics of Computation, 76(2007), 1425-1448.

[5] H. Ammari, E. Iakovleva and S. Moskow, Recovery of small inhomogeneities
from the scattering amplitude at a fixed frequency, SIAM J. Math. Anal,
34(4)(2003), 882-900.

[6] I. Babus’ka and A. Aziz, Survey Lectures on Mathematical Foundations of the
Finite Element Method, The Mathematical Foundations of the Finite Element
Method with Application to Partial Differential Equations, A. Aziz, Academic
Press, New York, 1973, 5—359.

[7] G. Bao, Y. Cao and H. Yang, Numerical solution of diffraction problems by a
least-square finite element method, Math. Methods Appl. Sci., 23(2000), 1073-
1092.

[8] G. Bao and Y. Chen, A nonlinear grating problem in diflractive optics, SIAM
J. Math. Anal., 28(2)(1997), 322—337.

[9] G. Bao, Y. Chen, and F. Ma, Regularity and stability for the scattering map of a
linearized inverse medium problem, J. Math. Anal. Appl., 247 (2000), 255-271.

[10] G. Bao, L. Cowsar and W. Masters, Mathematical Modeling in Optical Sciences,
SIAM Frontiers in Appl. Math, SIAM, Philadelphia, 2001.

[11] G. Bao and DC. Dobson, J .A. Cox, Mathematical studies in rigorous grating
theory, J. Opt. Soc. Amer., 12(A)(1995), 1029-1042.

81

[12] G. Bao, K. Huang, and G. Schmidt, Optimal design of nonlinear diffraction
gratings, J. Comp. Phys., 184(2003), 106-121.

[13] G. Bao, Y. Li and H. Wu, Numerical solution of nonlinear diffraction problems,
Journal Of Computational and Applied Mathematics, 190(1)(2006), 170—189.

[14] G. Bao, Y. Li and Z. Zhou, LP estimate of Maxwell’s equations in a bounded
domain, submitted.

[15] G. Bao and J. Liu, Numerical solution of inverse scattering problems with multi-
erperimental limited aperture data, SIAM J. Sci. Comput., 25(3) (2003), 1102-
1117.

[16] G. Bao and P. Li, Inverse medium scattering for the Helmholtz equation at fixed
frequency, Inverse Problems, 21(2005), 1621-1641.

[17] G. Bao and P. Li, Inverse medium scattering problems in near-field optics,
Journal of Computational Mathematics, 25(3)(2007), 252-265.

[18] G. Bao, A. Minut, and Z. Zhou, An LP estimate for Maxwell ’3 equations with
source term, Comptes Rendus Mathematique, 337(5) (2003), 365-370.

[19] S. S. Chaudhary, R. K. Mishra, A. Swarup and J. M. Thomas, Dielectric proper-
ties of normal and malignant human breast tissues at radiowave and microwave
frequencies, Indian J. Biochem. Biophys. 21(1984), 7679.

[20] Y. Chen, Inverse scattering via Heisenberg’s uncertainty principle, Inverse
Problems, 13(1997), 253-282.

[21] Z. Chen and H. Wu, A perfectly matched layer for the absorption of electro-
magnetic waves, J. Comput. Phys, 144(1994), 185-200.

[22] Z. Chen and H. Wu, An adaptive finite element method with perfectly matched
absorbing layers for the wave scattering by periodic structures, SIAM J. Numer.
Anal, 41(3)(2003), 799-826.

[23] M. H. Choi, T. J. Kao, D. Isaacson, G. J. Saulnier and J. C. Newell, A recon-
struction algorithm for breast cancer imaging with electrical impedance tomogra-

phy in mammography geometry, IEEE Transactions on Biomedical Engineering,
54(4)(2007), 700—710.

[24] R. Ciocan and N. Ida, Transmission line matrix model for detection of local

changes in permeability using a microwave technique, IEEE Transactions on
Magnetics, 40(2)(2004), 651—654.

82

 

[25] F. Collino and P. B. Monk, Optimizing the perfectly matched layer, Comput.
Methods Appl. Mech. Engrg., 164(1998), 157-171.

[26] D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering The-
ory, 2nd ed., Appl. Math. Sci., 93, Springer-Verlag, Berlin, 1998.

[27] M. Converse, E. J. Bond, B. D. VanVeen and S. C. Hagness, A computational
study of ultra-wideband versus narrowband microwave heperthermia for breast

cancer treatment, IEEE Transactions on Microwave Theory and Techniques,
54(5)(2006), 2169-2180.

[28] J. Elschner and G. Schmidt, Mathematical Methods in the Applied Sciences,
Math. Meth. Appl. Sci., 21(1998), 1297-1342.

[29] H. Engl, M. Hanke and A. Neubauer, Regularization of Inverse Problems,
Kluwer, Dordrecht, 1996.

[30] E. C. Fear, S. C. Hagness, P. M. Meaney, M. Okoniewski and M. A. Stuchly,
Enhancing breast tumor detection with near-field imaging, IEEE Microwave
Magazine, 3(1)(2002), 48-56.

[31] E. C. Fear, X. Li, S. C. Hagness and M. A. Stuchly, Confocal microwave imaging
for breast cancer detection: localization of tumors in three dimensions, IEEE
Transactions on Biomedical Engineering, 49(8)(2002), 812-822.

[32] A. Franchois and C. Pichot, Microwave imaging-complex permittivity recon-
struction with a Levenberg-Marquardt mehod, IEEE Transactions on Antennas
and Propagation, 45(2)(1997), 203-215.

[33] P. A. Franken, A. E. Hill, C. W. Peters, and G. Weinreich, Generation of optical
harmonics, Phys. Rev. Lett. 7(1961), 118119.

[34] D. Gilbarg and N. S. Trudinger, Elliptic Partial Differential Equation of Second
Order, Springer-Verlag, 1983.

[35] W. D. Hurt, Multiterm Debye dispersion relations for permittivity of muscle,
IEEE Trans. Biomed. Eng. 32(1985), 6064.

[36] J. Jin, The Finite Element Methods in Electromagnetics, John Wiley & Sons,
2002.

[37] A. Kirsch and P. Monk, Convergence analysis of a coupled finite element and

spectral method in acoustic scattering, IMA Journal of Numerical Analysis,
9( 1990), 425-447.

83

 

[38] A. Kirsch and P. Monk, An analysis of the coupling of finite-element and
Nystro'm methods in acoustic scattering, IMA Journal of Numerical Analysis,
14(1994), 523-544.

[39] R. G. Keys and A. B. Weglein, Generalized linear inversion and the first Born
theory for acoustic media, J. Math. Phys., 24(1983), 1444-1449.

[40] R. Kress, Linear Integral Equations, Spring-Verlag, Berlin Heidelberg New
York, 1989.

[41] X. Li and SC. Hagness, A confocal microwave imaging algorithm for breast
cancer detection, IEEE Microwave Wireless Components Lett., 11(2001), 130-
I32.

[42] S. P. Murarka and M. C. Peckerar, Electronic Materials Science and Processing,
Academic Press, Inc., San Diego, CA (1989).

[43] W. Nakagawa, R. Tyan, and Y. Fainman, Analysis of enhanced second har-
monic generation in periodic nanostructures using modified rigorous coupled

wave analysis in the undepleted pump approximation, J. Opt. Soc. Am. A,
19(9)(2002), 1919-1928.

[44] R. Nilavalan, A. Gbedemah, I. J. Craddock, X. Li and S. C. Hagness, Numer-
ical investigation of breast tumor detection using multi-static radar, Electronic
Letters Online, 39(25)(2003).

[45] J. C. Nédélec, Acoustic and Electromagnetic Equations, Springer-Verlag, 2001.

[46] L. Nirenberg, Functional analysis, Notes by Sibner, Courant Institute of Math-
ematical Sciences, New York University, 1960—1961.

[47] R. Petit, Electromagnetic Theory of Gratings, Topics in Current Physics, Vol.
22, Springer-Verlag, Heidelberg, 1980.

[48] E. Popov and M. Neviere, Surface-enhanced second-harmonic generation in
nonlinear corrugated dielectrics: New theoretical approaches, J. Opt. Soc. Amer.
B 11(1994), 1555-1564.

[49] R. Reinish and M. Neviere, Electromagnetic theory of diffraction in nonlinear
optics and surface-enhanced nonlinear optical effects, Phys. Rev. B 28(1983),
1870-1885.

[50] R. Reinish, M. Neviere, H. Akhouayri, J. Coutaz, D. Mayster, and E. Pic, Grat-
ing enhanced second harmonic generation through electromagnetic resonances,

Opt. Engrg 27(1998), 961-971.

84

[51] L. R. Scott and S. Zhang, Finite element interpolation of nonsmooth functions
satisfying boundary conditions, Math. Comp, 54(1990), 483-493.

[52] Y. R. Shen, Principles of Nonlinear Optics, Wiley, 1984.

[53] P. Stoica, Z. Wang and J. Li, Robust Capon Beamforming, IEEE Signal Pro-
cessing Letters, 10(6)(2003), 172-175.

[54] E. Turkel and A. Yefet, Absorbing PML boundary layers for wave-like equations,
Appl. Numer. Math. 27(1998), 533-557.

[55] C. Yu, Z. Zhou and M. Zhuang, Extension of the Helmholtz Equation Least-
Square Method in Inverse Acoustics, to appear.

[56] Y. Zou and Z. Guo, A review of electrical impedance techniques for breast cancer
detection, Med. Eng. Phys, 25(2003), 79—90.

85

 

11111111111111]1111'