Solving the multi-frequency electromagnetic inverse source problem by the Fourier method

# Solving the multi-frequency electromagnetic inverse source problem by the Fourier method

Guan Wang 111School of Mathematics, Jilin University, Changchun, P. R. China. Email: wangguan@neepu.edu.cn    Fuming Ma 222School of Mathematics, Jilin University, Changchun, P. R. China. Email: mfm@jlu.edu.cn    Yukun Guo 333Corresponding author. Department of Mathematics, Harbin Institute of Technology, Harbin, P. R. China. Email: ykguo@hit.edu.cn    Jingzhi Li 444Department of Mathematics, Southern University of Science and Technology, Shenzhen, P. R. China. Email: li.jz@sustc.edu.cn
###### Abstract

This work is concerned with an inverse problem of identifying the current source distribution of the time-harmonic Maxwell’s equations from multi-frequency measurements. Motivated by the Fourier method for the scalar Helmholtz equation and the polarization vector decomposition, we propose a novel method for determining the source function in the full vector Maxwell’s system. Rigorous mathematical justifications of the method are given and numerical examples are provided to demonstrate the feasibility and effectiveness of the method.

Keywords: inverse source problem, Fourier method, time-harmonic, Maxwell’s equations, multi-frequency

## 1 Introduction

Inverse source problems arise naturally in various application areas such as biomedical tomography, non-invasive detection, target tracking and antenna synthesis. In recent years there have been tremendous advances in the theoretical understanding and numerical treatment of inverse source problems, see, e.g., [2, 5, 7, 22, 13, 6, 23, 19, 14, 15] and the references therein for relevant studies.

The inverse source problems could be posed in the frequency or time domain. The inverse source problem of determining a source in the Helmholtz equation has been studied in [7, 9, 5, 3, 4, 20, 8, 11]. For inverse source problems for time-harmonic Maxwell’s system, we refer to He and Romanov [23], Ammari et al. [13] and Albanese and Monk [22] for the investigation of the localization of brain activities. Rodríguez et al. [1] studied the inverse source problem for the eddy current of Maxwell’s equations. For recent works on dynamic inverse source problem of imaging the trajectory of a moving point source, we refer to [25, 26].

This work is concerned with the inverse source problem of determining the electric current in the time-harmonic Maxwell’s system from the multifrequency near-field measurements. For increasing stability analysis concerning the inverse source problem with multi-frequencies, we refer to [17, 18, 10].

A numerical method based on the Fourier expansion of the source has been proposed to solve the multi-frequency inverse source problem for the Helmholtz equation in [8]. This Fourier method has been extended to the far-field cases of acoustic inverse source problem [24]. The Fourier method is easy to implement with computational efficiency. Our goal in this paper is to extend the Fourier method from the scalar Helmholtz equation to full vector Maxwell’s equations. Due to the complexity of the profound structure of Maxwell’s equations, both the theoretical analysis and numerical implementation in the current study are radically much more challenging than the scalar counterpart presented in [8, 24]. In particular, one would encounter much more complicated technical difficulties in establishing the stability estimates.

The rest of this paper is organized as follows. In the next section, we introduce the model problem and recall some suitable Sobolev spaces. Section 3 is devoted to approximating the source function by the Fourier expansion, a uniqueness result on the approximate source, and then a numerical method for solving the inverse source problem are presented. The result of error estimate of measurements with noise is derived in Section 4. Finally, several numerical examples are provided to show the effectiveness of our method in Section 5.

## 2 Problem setting

In this section, we shall first give a description of the multifrequency inverse source problem under consideration. Then, some notation of relevant Sobolev spaces will be reviewed as the prerequisites for the theoretical analysis in the subsequent sections. Throughout this paper, we adopt the convention of using bold and non-bold fonts for vector and scalar values, respectively.

### 2.1 Model problem

In this paper, we consider the inverse source problem of determining a radiating current density excitation in the time-harmonic Maxwell’s equations in

 ∇×E−ikH=0, (2.1) ∇×H+ikE=J, (2.2)

 lim|x|→∞(∇×E×^x−ikE)=0, (2.3)

which holds uniformly for all directions , with . Here , and denote, respectively, the electric field, the magnetic field and the current density vector. The wavenumber/frequency of the problem is defined by the positive constant , with denoting the angular frequency and and the electric permittivity and the magnetic permeability in vacuum. Eliminating in (2.1)-(2.2) then the time-harmonic Maxwell’s equations can be written as

 ∇×∇×E−k2E=ikJ. (2.4)

We consider the multi-frequency inverse source problem with the following assumptions:
(i) The current density vector is independent of and

 J∈(L2(R3))3,suppJ⊂V0, (2.5)

where ia a cube and .
(ii) The current source is expressed in the following form

 J=pf+p×∇g, (2.6)

where and . Here, the polarization vector is assumed to be known and belongs to the following admissible set

 P:={q∈R3| q×l≠0, ∀l∈Z3∖{0}}.

In what follows, let be the measurement surface. Then the multi-frequency inverse source problem can be stated as follows:

###### Problem 2.1 (Multi-frequency Inverse Source Problem).

Given a fixed polarization vector and a finite set of admissible wavenumbers, reconstruct the source function in the form (2.6) from the measured data

 {E(x;k,p)∣∀x∈ΓR, ∀k∈K}.

where indicates the dependence of the radiated field on the wavenumber and the polarization vector .

In this paper, we derive a novel and effective numerical method for solving the multi-frequency inverse source problem. Our method is based on the Fourier expansion of the source function . We prove that a finite number of Fourier coefficients can be uniquely determined from the measured fields corresponding to a set of wavenumbers which are properly chosen. Meanwhile, we establish the explicit formula for each coefficient.

### 2.2 Tangential vector spaces

In this subsection, we recall some essential ingredients for tangential vector Sobolev spaces and we refer to [21] for more relevant details.

Let be a generic smooth and closed surface. Then the space of tangential vector fields is defined by

 L2t(Γ):={u∈(L2(Γ))3∣u⋅ν=0 on Γ},

where denotes the unit outward normal of .

Let be the unit sphere and be the unit vectors of the spherical coordinates, where , and are the spherical harmonics.

Denote the vector spherical harmonics

 Umn(^x)=1√n(n+1)∇S2Ymn(^x),Vmn(^x)=^x×Umn(^x),^x∈S2,

where the surface gradient is defined by

 ∇S2:=∂∂θeθ+1sinθ∂∂φeφ.

Then the set of all vector spherical harmonics form a complete orthonormal basis of .

Now we introduce the space of three-dimensional vector functions by

 H(curl;BR):={u∈(L2(BR))3∣∇×u∈(L2(BR))3},

and its trace space . We suppose has representation

 ζ(x)=∞∑n=1n∑m=−n(αn,mUmn(^x)+βn,mVmn(^x)),

then its norms in and are defined by

 ∥ζ∥L2t(ΓR):=(∞∑n=1n∑m=−n(|αn,m|2+|βn,m|2))1/2, ∥ζ∥H−1/2(Div,ΓR):=(∞∑n=1n∑m=−n(√n(n+1)|αn,m|2+|βn,m|2√n(n+1)))1/2.

In order to compute the data from , we need the electric-to-magnetic Calderón operator which is defined by

 Gζ:=∞∑n=1n∑m=−n⎛⎝k2Rh(1)n(kR)z(1)n(kR)βn,mUmn+1Rz(1)n(kR)h(1)n(kR)αn,mVmn⎞⎠, (2.7)

where is the spherical Hankel function of the first kind of order , and .

## 3 Fourier approximation and uniqueness

In this section, we will approximate the functions by Fourier expansions and establish explicit formulas for the Fourier coefficients.

We begin this section with some notation and definitions that are used in the paper. Without loss of generality, let

 V0:=(−L2,L2)3,

and introduce the following Fourier basis functions

 ϕl(x):=exp(i2πLl⋅x),l∈Z3,

then the functions and can be written as Fourier expansions of the forms

 f=∑|l|≥0alϕl,g=∑|l|≥1blϕl,

with the Fourier coefficients

 al=1L3∫V0f(x)¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x)dx,|l|≥0, (3.1) bl=1L3∫V0g(x)¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x)dx,|l|≥1, (3.2)

where the overbar denotes the complex conjugate. Thus the source has the following form

 J=pf+p×∇g=a0p+∑|l|≥1(alp+2πiLbl(p×l))ϕl. (3.3)

We introduce the Sobolev space of order

 (Hσp(V0))3:={ Φ∈(Hσ(V0))3∣Φ=pf+p×∇g, f∈Hσ(V0), g∈Hσ+1(V0), p∈S2},

with the norm

 ∥Φ∥p,σ=⎛⎝L3∑|l|≥0(1+|l|2)σ|al|2+4π2L∑|l|≥1(1+|l|2)σ|p×l|2|bl|2⎞⎠1/2, (3.4)

where has the Fourier expansion of the form

 Φ=a0p+∑|l|≥1(alp+2πiLbl(p×l))ϕl.

Now we introduce the definition of admissible polarization vectors in the source function (2.6).

For an admissible polarization vector , the decomposition

 p =p⋅l|l|2l+p⋅l×(p×l)|l×(p×l)|2l×(p×l) =p⋅l|l|2l+1|l|2l×(p×l) (3.5)

is crucial for computing the Fourier coefficients. For simplicity, we denote , and it is clear that form an orthogonal basis of . A direct computation shows that

 ∇×∇×(ulϕl)=−Δ(ulϕl)+∇∇⋅(ulϕl)=4π2L2|l|2ulϕl, (3.6)

where or .

Combining and , we obtain

 J= pf+p×∇g = a0p+∑|l|≥1(al|l|2((p⋅l)l+wl)+2πiLblvl)ϕl. (3.7)

We can approximate by the truncated Fourier expansion

 JN=a0p+∑1≤|l|≤N(al|l|2((p⋅l)l+wl)+2πiLblvl)ϕl. (3.8)

Then we have the following approximation result.

###### Theorem 3.1.

Let be a function in , then the following estimate holds

 ∥J−JN∥p,μ≤Nμ−σ∥J∥p,σ,0≤μ≤σ. (3.9)
###### Proof.

From and , we obtain

 ∥J−JN∥2p,μ = ∑|l|>N(1+|l|2)μ(L3|al|2+4π2L|p×l|2|bl|2) ≤ 1(1+N2)σ−μ∑|l|>N(1+|l|2)σ(L3|al|2+4π2L|p×l|2|bl|2) ≤ 1(1+N2)σ−μ∥J∥2p,σ,

which leads to the estimate . ∎

Motivated by Theorem 3.1, the inverse source problem under concern is to determine an approximate source in the form (3.8) with a fixed admissible polarization vector from the measurements on at a finite set of wavenumbers . To this end, we will establish explicit formulas for the Fourier coefficients and . Before we show that, we first state a uniqueness result.

###### Theorem 3.2.

Let and , then the Fourier coefficients and of in and can be uniquely determined by the measurements .

###### Proof.

Let be a source that produces the data . Using the electric-to-magnetic Calderón operator [21, p.249], we know that the source also produces the data .

For every and , multiplying equation by and and integrating over , then using , we obtain

 ik∫V0J(x)⋅vl¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x)dx = ∫BR(∇×∇×E(x)−k2E(x))⋅vl¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x)dx = ∫BRE(x)⋅(∇×∇×(vl¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x))−k2vl¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x))dx +∫ΓR(^x×∇×E(x)⋅(vl¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x))+^x×E(x)⋅∇×(vl¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x)))ds(x) = ∫ΓR(^x×∇×E(x)⋅(vl¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x))+^x×E(x)⋅∇×(vl¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x)))ds(x),

and

 ik∫V0J(x)⋅wl¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x)dx = ∫ΓR(^x×∇×E(x)⋅(wl¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x))+^x×E(x)⋅∇×(wl¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x)))ds(x).

On the other hand, using , we have

 ∫V0J(x)⋅(wl¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x))dx=L3|wl|2al, ∫V0J(x)⋅(vl¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x))dx=2πiL2|vl|2bl.

Then, we obtain

 al= 1ik|vl|2L3∫ΓR(^x×∇×E(x)⋅(wl¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x)) +^x×E(x)⋅∇×(wl¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x)))ds(x), (3.10) bl= −12πkL2|vl|2∫ΓR(^x×∇×E(x)⋅(vl¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x)) +^x×E(x)⋅∇×(vl¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x)))ds(x), (3.11)

then the Fourier coefficients and of are uniquely determined. This completes the proof. ∎

It is clear that we can reconstruct all Fourier coefficients satisfying from the measurements on the wavenumbers in by and . Now we are in a position to introduce the numerical algorithm to reconstruct .

###### Theorem 3.3.

Let be a small positive constant. Take , , and

 ϕl∗(x):=exp(i2πLl∗⋅x).

Let be defined as

 aN0:= λπsinλπ(1ik∗L3|p×l∗|2∫ΓR(^x×∇×E⋅(w∗¯¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl∗(x)) +^x×E⋅∇×(w∗¯¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl∗(x)))ds(x)−N∑|j|=1sin(j−λ)π(j−λ)πaj), (3.12)

where for . Then the following estimate holds

 |a0−aN0|≤2λ√NL3/2∥J∥p,0. (3.13)
###### Proof.

It is easy to check that for , we have

 ∫V0ϕl(x)¯¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl∗(x)dx=⎧⎪⎨⎪⎩L3sin(j−λ)π(j−λ)π,l=(j,0,0)⊤,0,otherwise, (3.14)

and

 ∇×∇×(w∗ϕl∗)=4π2L2|l∗|2w∗ϕl∗. (3.15)

Similar to the proof of Theorem 3.2, multiplying equation by and integrating over , then using , we obtain

 ∫V0J(x)⋅w∗¯¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl∗(x)dx = 1ik∗∫ΓR(^x×∇×E(x)⋅(w∗¯¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl∗(x))+^x×E(x)⋅∇×(w∗¯¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl∗(x)))ds(x).

On the other hand, using and , we have

 ∫V0J(x)⋅w∗¯¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl∗(x)dx = L3|p×l∗|2sinλπλπa0+∞∑|j|=1L3|p×l∗|2sin(j−λ)π(j−λ)πaj,

Hence we have the following formula for

 a0= λπsinλπ(1ik∗L3|p×l∗|2∫ΓR(^x×∇×E(x)⋅(w∗¯¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl∗(x)) +^x×E(x)⋅∇×(w∗¯¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl∗(x)))ds(x)−∞∑|j|=1sin(j−λ)π(j−λ)πaj), (3.16)

and the approximate formula for

 aN0= λπsinλπ(1ik∗L3|p×l∗|2∫ΓR(^x×∇×E(x)⋅(w∗¯¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl∗(x)) +^x×E(x)⋅∇×(w∗¯¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl∗(x)))ds(x)−N∑|j|=1sin(j−λ)π(j−λ)πaj).

Furthermore,

 |a0−aN0|= = ∞∑|j|=N+1∣∣∣λj−λaj∣∣∣ ≤ ⎛⎝∞∑|j|=N+1∣∣∣λj−λ∣∣∣2⎞⎠1/2⎛⎝∞∑|j|=N+1|aj|2⎞⎠1/2 ≤ 2λ√NL3/2∥J∥p,0.

The proof is completed. ∎

## 4 Stability analysis

Let be the Fourier approximation in case that the measured data is perturbed with the noise level . In this section, we will establish the calculation formula and give an estimate of the error between and .

In general, only data with noise can be measured, and the data cannot be computed from by using the electric-to-magnetic Calderón operator (2.7). In this paper, we follow [8] to overcome this difficulty. For simplicity, we write for with a fixed polarization .

Using spherical coordinates , we can write the measured noisy data on as

 ^x×Eδ(x;k)∣∣ΓR=∞∑n=1n∑m=−n(αδk,n,mUmn+βδk,n,mVmn),

where

 αδk,n,m=∫S2^x×Eδ(x;k)⋅¯¯¯¯¯¯¯¯¯Umnds(x), (4.1) βδk,n,m=∫S2^x×Eδ(x;k)⋅¯¯¯¯¯¯¯¯¯Vmnds(x). (4.2)

Then, for , we define

 ^x×E(x;k)=∞∑n=1n∑m=−n⎛⎝h(1)n(kr)h(1)n(kR)αδk,n,mUmn+Rrz(1)n(kr)z(1)n(kR)βδk,n,mVmn⎞⎠.

Take , and let , we have

 ^x×Eδ∣∣Γρ=∞∑n=1n∑m=−n⎛⎝h(1)n(kρ)h(1)n(kR)αδk,n,mUmn+Rρz(1)n(kρ)z(1)n(kR)βδk,n,mVmn⎞⎠, (4.3) (4.4)

Now, we propose modified formulas for the approximate Fourier coefficients

 aδl=1ik|v|2L3∫Γρ(^x×∇×Eδ⋅(w¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x))+^x×Eδ⋅∇×(w¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x)))ds(x), (4.5) bδl=−12πkL2|v|2∫Γρ(^x×∇×Eδ⋅(v¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x))+^x×Eδ⋅∇×(v¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl(x)))ds(x), (4.6)

for , where , and

 aN,δ0= +^x×Eδ⋅∇×(w∗¯¯¯¯¯¯¯¯¯¯¯¯¯¯ϕl∗(x)))ds(x)−N∑j=1sin(j−λ)π(j−λ)πaδj), (4.7)

where for .

By using and , we define the function by

 JδN(x):=aN,δ0p+∑1≤|l|≤N(aδlp+2πiLbδl(p×l))ϕl(x), (4.8)

which will be an approximation of the source .

Let the measured noisy data satisfy

 ∥^x×E(x;k)−^x×Eδ(x;k)∥L2t(ΓR)≤δ∥^x×E(x;k)∥L2t(ΓR), (4.9)

where . First we give an estimate on .

###### Lemma 4.1.

Let , then there exists a positive constant such that

 ∥^x×E(x;k)∥L2t(ΓR)≤kM∥J∥p,0, (4.10)

where depends only on and .

###### Proof.

Introduce the fundamental solution of the Helmholtz equation

 Φk(x,y)=exp(ik|x−y|)4π|x−y|,x≠y,x,y∈R3.

Then the unique solution to equation with the radiation condition can be written as (see, e.g., [12])

 E(x;k)= ik∫V0(I+∇∇⊤k2)Φk(x,y)J(y)dy = ik∫V0Φk(x,y)((1+ik−1k2|x−y|2)I +(3k2|x−y|2−3ik|x−y|−1)(x−y)(x−y)⊤|x−y|2)J(y)dy, (4.11)

where and denotes the identity matrix.

Let

 τ0:=min{|x−y|∣x∈ΓR, y∈V0}=R−√32L,

then for all and , we obtain

 ∣∣ ∣∣1+ik−1k2|x−y|2∣∣ ∣∣≤1+1k+1k2τ20, ∣∣ ∣∣3k2|x−y|2−3ik