1. Introduction. The inverse source problem for the Maxwell equations is of fundamental interest and considerable practical importance, with applications ranging from geophysics to biomedical imaging [1,12,13,24,38]. The problem is usually stated in the following form: determine the electric current density from boundary measurements of the electric and magnetic fields. It is well known that this problem is underdetermined and does not admit a unique solution, due to the existence of socalled nonradiating sources [20]. However, if the source is spatially localized or some other a priori information is available, it is often possible to characterize the source to some extent [2]. Such a method is applied to the localization of low-frequency electric and magnetic signals originating from current sources in the brain or heart [34].

In this paper, we propose an alternative approach to the electromagnetic inverse source problem. In this approach, which is an extension of the authors' previous work on the acousto-electric inverse source problem for static fields [30], a wavefield is used to control the material properties of a medium of interest, which is then probed by a second wavefield. Also, see related work on hybrid imaging [3,7,8,9,10,11,14,16,17,18,19,23,25,26,27,28,29,30,32,33,35,36,37]. Here the electric current density as well as the conductivity, electric permittivity, and magnetic permeability are spatially modulated by an acoustic wave. In this manner, we find that it is possible to uniquely recover the current density from boundary measurements of the fields with Lipschitz stability.

The remainder of this paper is organized as follows. In section 2 we introduce a model for the acoustic modulation of the current density and the material parameters. In sections 3 and 4 this model is used to formulate the inverse source problem and thereby derive an internal functional from which the source may be recovered. Numerically simulated reconstructions are given in section 5. Finally, our conclusions are presented in section 6.

2. Model. We begin by developing a simple model for acoustic modulation of the electrical current density and material parameters, following the approach of [30]. We begin by considering the time-harmonic Maxwell equations in a bounded domain \Omega \subset R 3 :

(1)

i\omega \varepsi E + \nabla \times H = J + \sigma E in \Omega , - i\omega \mu H + \nabla \times E = 0 in \Omega .

We also impose the impedance boundary condition

(2) H \times \n -\lambda (\n \times E) \times \n = 0 on \partial\Omega , which arises since \Omega is taken to be enclosed by a good conductor. Here the vector functions J, E, and H are the current density, the electric field, and the magnetic field, respectively. The scalar functions \varepsi , \mu , \sigma , and \lambda are the electric permittivity, magnetic permeability, conductivity, and surface impedance, respectively. The vector \ n is the outward unit outward normal to \Omega and \omega is a fixed frequency. Note that in the above, we do not write the equations governing the divergence of E and H which are not needed in what follows.

The inverse source problem is to reconstruct the source J from boundary measurements, assuming that the coefficients \mu , \varepsi , \sigma , \lambda are known. A typical measurement is the tangential electric field on the boundary [1]: (3) g := (\n \times E) \times \n| \partial\Omega .

This problem does not have a unique solution [20]. That is, distinct sources may give rise to the same boundary measurements.

Remark 2.1. An alternative measurement is h := \n \times H| \partial\Omega . Knowledge of g is equivalent to knowledge of h when the impedance boundary condition ( 2) is taken into account, since h = - \lambda g on \partial\Omega .

We now examine the effect of acoustic modulation. Two approaches have been proposed in the literature: modulation with plane waves or spherical waves. The latter has been utilized in electromagnetic tomography [4] and acousto-optic imaging [5,6]. The former has been employed in acousto-optic imaging [10], bioluminescence tomography [7,11,19], optical tomography [28,29], and acousto-electric imaging [30]. In this paper, we make use of plane-wave modulation for ease of mathematical analysis. We note that both types of modulation are equivalent since spherical waves can be synthesized from plane waves by the Weyl formula.

Following [10,30], we consider a system of charge carriers in a fluid, in which a small-amplitude acoustic plane wave propagates. It follows that the current density J \delta is modulated according to

where J is the conductivity in the absence of the acoustic wave, \delta \ll 1 is a small parameter that is proportional to the acoustic pressure, \gamma J is the elasto-electric constant, k is the wave vector of the acoustic wave and \varphi is its phase. Likewise, the conductivity \sigma \delta and permittivity \varepsi \delta are also modulated:

where \sigma and \varepsi are the unmodulated conductivity and permittivity, and the constants \gamma \varepsi , \gamma \sigma are known as the elasto-electric constants. For simplicity we assume that the impedance \lambda is not affected by the acoustic modulation. It follows that the modulated electric and magnetic fields E \delta and H \delta satisfy the modified Maxwell equations (5) i\omega \varepsi \delta E \delta + \nabla \times H \delta = J \delta + \sigma \delta E \delta in \Omega , - i\omega \mu H \delta + \nabla \times E \delta = 0 in \Omega , together with the boundary condition (6) H \delta \times \n -\lambda (\n \times E \delta ) \times \n = 0 on \partial\Omega .

The corresponding boundary measurement becomes (7) g \delta := (\n \times E \delta ) \times \n| \partial\Omega .

3. Internal functional. In this section, we derive the internal functional from boundary measurements of the electric field. We also introduce the necessary function spaces and specify certain technical requirements on the conductivity and permittivity.

3.1. Function spaces. We will use the following standard spaces to discuss the wellposedness of the Maxwell's equations [15]. Let \Omega \subset R 3 be an open bounded set with a C 1,1 boundary, and let

The norm on H(curl, \Omega ) is given by

The two tangential trace maps \Gamma \tau and \Pi \tau have the following definitions:

\Gamma \tau : H(curl, \Omega ) \rightar H - 1/2 (div, \partial\Omega ), u \mapsto \rightar \n \times u| \partial\Omega and \Pi \tau : H(curl, \Omega ) \rightar H - 1/2 (curl, \partial\Omega ), u \mapsto \rightar (\n \times u) \times \n| \partial\Omega , where

Here div \partial\Omega is the surface divergence and curl \partial\Omega is the surface curl. The two spaces H - 1/2 (div, \partial\Omega ) and H - 1/2 (curl, \partial\Omega ) are dual to each other. To handle the impedance boundary condition, we define the tangential trace of a vector field (8) \bfitphi T := \Pi \tau (\bfitphi ) = (\n \times \bfitphi ) \times \n| \partial\Omega , and the space X = \bigl\{ u \in H(curl, \Omega ) : u T \in (L 2 (\partial\Omega )) 3 \bigr\} .

The norm on X is

We denote the dual space of X by X \ast . We denote the (L 2 (\Omega )) 3 -inner product by

and the (L 2 (\partial\Omega )) 3 inner product by \langle u, v\rangle (L 2 (\partial\Omega )) 3 := \int \partial\Omega u \cdot vdx, u, v \in (L 2 (\partial\Omega )) 3 , where v denotes the complex conjugate of v. We denote the dual paring of u \in H - 1/2 (div, \partial\Omega ) and v \in H - 1/2 (curl, \partial\Omega ) by \langle u, v\rangle .

3.2. Assumptions and weak formulation. We will make the following assumptions throughout this paper.

A-1. The domain \Omega is an open bounded connected domain in R 3 with C 1,1 boundary. A-2. The medium is nonmagnetic with \mu = \mu 0 in \Omega , where \mu 0 is the magnetic permeability in vacuum. The coefficients \varepsi and \sigma are real piecewise H 3 (\Omega ) functions. A-3. There exists positive constants K 1 and K 2 , such that (9)

and the conductivity \sigma is nonzero. A-4. The source J is an (L 2 (\Omega )) 3 vector field and is compactly supported in \Omega .

Here the space H 3 (\Omega ) is defined as follows. The domain \Omega can be partitioned into finitely many connected subdomains of Lipschitz boundary, denoted by \Omega p for p = 1, . . . , N , and a function is piecewise H 3 if its restriction on each \Omega p is H 3 (\Omega p ).

Remark 3.1. We conclude from A-2 that \varepsi and \sigma are piecewise C 1 by the Sobolev embedding theorem. We conclude from A-3 that K 1 > \varepsi \delta > K 2 > 0 and \sigma \delta \geq 0, so long as \delta is sufficiently small. The modulated Maxwell equations (5) and the impedance boundary condition (6) can be written in terms of only the electric field: (10) \nabla \times \nabla \times E \delta -\mu (\omega 2 \varepsi \delta + i\omega \sigma \delta )E \delta = i\omega \mu J \delta in \Omega , which is subject to the impedance boundary condition (11) \biggl( 1 \mu \nabla \times E \delta \biggr) \times \n -i\omega \lambda (\n \times E \delta ) \times \n = 0 on \partial\Omega .

We say E \delta \in X is a weak solution of (10) obeying the impedance boundary condition (11) if for all \bfitphi \in X, \biggl( 1 \mu \nabla \times E \delta , \nabla \times \bfitphi \biggr)

It follows from Assumptions A1--A4, that the weak solution E \delta \in X exists and is unique [31].

3.3. Internal functional. We now derive the internal functional for both classical and weak solutions. To proceed, we consider the fields F and G which obey the Maxwell equations without sources: (13) i\omega \varepsi F + \nabla \times G = \sigma F in \Omega ,

- i\omega \mu G + \nabla \times F = 0 in \Omega , along with the impedance boundary condition

where g \in (L 2 (\partial\Omega )) 3 . Equivalently, ( 14)

\times \n -\lambda (\n \times F) \times \n = g on \partial\Omega .

Note that ( 13) are explicitly solvable, since the required coefficients are known. Next, we take the inner product of ( 5) with F, the inner product of ( 14) with E \delta , and then subtract to obtain \nabla \times \nabla \times E \delta \cdot F -\nabla \times \nabla \times

Integrating the above result over \Omega and using the vector identity (\nabla \times A) \cdot B = \nabla \cdot (A \times B) + (\nabla \times B) \cdot A, we find that

We now integrate by parts the divergence terms, which, using the relations 1 \mu \nabla \times E \delta = i\omega H \delta and 1 \mu \nabla \times

Note that the boundary integral only depends on the tangential components of the fields H \delta , E \delta , F, and G, which are known from the boundary measurements (7). Therefore, the left-hand side can be determined from the boundary measurements. For the right-hand side, we consider the asymptotic expansion in the small quantity \delta . The O(1) term is

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 04/11/23 to 132.174.252.171 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy

The O(\delta ) term is of the form ( 16)

Varying k and \varphi in ( 16), and performing the inverse Fourier transform, we obtain the internal functional ( 17)

which is known at every point in \Omega .

Remark 3.2. In place of ( 6), where we assumed that the impedance \lambda of the boundary is not modulated by the acoustic field, we may alternatively consider the case that \lambda is modulated. In the latter case, the internal function (17) remains the same, since the right-hand-side of (15) remains known.

We make the following hypothesis to extract more information from the internal function (17).

Hypothesis 3.3. There exists a finite open cover \{ \Omega \alpha \} \alpha \in \Lambda of \Omega , such that for each \alpha \in \Lambda , there exist three solutions to (14) in \Omega , denoted F 1\alpha , F 2\alpha , and F 3\alpha , that are linearly independent on \Omega \alpha .

The hypothesis means that, in each \Omega \alpha , we can form the nonsingular matrix [F 1\alpha , F 2\alpha , F 3\alpha ], where F j\alpha is the jth column, j = 1, 2, 3. Let Q j\alpha be the internal functional defined as in (17), with F replaced by F j\alpha . Given the row vector

where we view E and J as column vectors, and T denotes the transpose. Therefore, if we define Q \in (L 2 (\Omega )) 3 by specifying its restrictions according to

then Q is well defined since both E and J are global vector fields over \Omega , and we have

Note that we view Q as a vector-valued internal functional.

4. Inverse problem and internal functional. It follows from the above discussion that the inverse problem consists of recovering the source current J from the internal functional Q. In this section we will derive a reconstruction procedure that uniquely recovers J with Lipschitz stability. The analysis depends critically on whether the constant \gamma J vanishes.

4.1. Case I: \bfitgam \bfitJ = 0. In this situation, the equality (18) does not involve J directly.

Proposition 4.1. Suppose the assumptions A1--A4 and the hypothesis (3.3) hold. If \gamma J = 0, then we have the following two subcases:

(I.1) If \Omega \subsete supp (\omega 2 \varepsi \gamma \varepsi + i\omega \sigma \gamma \sigma ), then the source J is uniquely determined with the stability estimate

for some constant C > 0 independent of J, \J . (I.2) If \Omega \not \subsete supp (\omega 2 \varepsi \gamma \varepsi +i\omega \sigma \gamma \sigma ), then the source J cannot be uniquely determined. Moreover, whenever J is uniquely determined, there are explicit reconstruction procedures.

Proof. If \Omega \subsete supp (\omega 2 \varepsi \gamma \varepsi + i\omega \sigma \gamma \sigma ), then (18) implies

This uniquely determines the weak solution E \in X everywhere in \Omega . Consequently, H and J are also uniquely determined via the Maxwell's equations (1). Note that all of these procedures are constructive: given Q, we compute E from the above equality, and then J from (1).

The stability can be derived as follows. If there is another source \J with corresponding electric field \Ẽ , and vector internal functional \Q defined as in (18), then E -\Ẽ is a weak solution of the Maxwell equations. That is,

for all \bfitphi \in X. As the coefficients in this weak formulation are all bounded, there exists a constant C > 0 such that

We deduce that

If \Omega supp (\omega 2 \varepsi \gamma \varepsi + i\omega \sigma \gamma \sigma ), there exists an open set D \subsete \Omega \setminu supp (\omega 2 \varepsi \gamma \varepsi + i\omega \sigma \gamma \sigma ). For any compactly supported smooth function \phi \in C \infty c (D), if (E, H) solves (1), then (E + \nabla \phi , H) solves (1) with J replaced by J + (i\omega \varepsi -\sigma )\nabla \phi . Moreover, since (E + \nabla \phi , H)| \partial\Omega = (E, H), these two pairs both satisfy the boundary condition (2) and produce identical measurement (3). This means that sources of the form J \phi := (i\omega \varepsi - \sigma )\nabla \phi are nonradiating. Thus the source J cannot be uniquely determined from the boundary measurement (3).

4.1.1. Increased regularity. The stability estimate for the subcase \gamma J = 0 and \Omega \subsete supp (\omega 2 \varepsi \gamma \varepsi + i\omega \sigma \gamma \sigma ) is in terms of the X \ast norm, which follows because J was obtained from E using a weak formulation. When the reconstructed E and H are smooth enough, for example, when E \in (H 2 (\Omega )) 3 , we can utilize the strong formulation to control J in (L 2 (\Omega )) 3 in terms of the higher order derivatives of the internal data. Proposition 4.2. Suppose the assumptions A1--A4 and the hypothesis (3.3) hold. Suppose, in addition, that \varepsi , \sigma \in C 1,1 (\Omega ). If \gamma J = 0 and \Omega \subsete supp (\omega 2 \varepsi \gamma \varepsi + i\omega \sigma \gamma \sigma ), then the following stability estimate holds for any two compactly supported sources J, \J \in (H 2 (\Omega ))

Here \Omega 1 is an open set compactly contained in \Omega such that supp J \subset \Omega 1 and supp \J \subset \Omega 1 , and the constant C > 0 is independent of J, \J .

Proof. Define

Then u solves

The following stability estimate is immediate:

It remains to show that the quantity

is finite. To proceed we will employ an interior regularity estimate for elliptic equations. Denote by a := \omega 2 \varepsi > 0, b := i\omega \sigma \geq 0, and f := i\omega (J -\J ). Then take the divergence of (19) to obtain (20) \nabla \cdot u = -\nabla (a + ib) \cdot u + \nabla \cdot f a + ib .

Using the identity \nabla \times \nabla \times u = \nabla (\nabla \cdot u) -\Delta u, we obtain the following elliptic system:

For each component of u, the left-hand side of the above defines a second order elliptic operator with constant coefficients, so we can apply an interior regularity estimate [ Combing this with the fact that u \in (L 2 (\Omega )) 3 , we obtain from [21, Lemma 6.32] that u \in (H 1 loc (\Omega )) 3 . This increased regularity implies that \nabla \biggl[ \nabla (a + ib) \cdot u + \nabla \cdot f a + ib \biggr] +\mu (a + ib)u + \mu f \in (L 2 loc (\Omega )) 3 .

Applying [21, Lemma 6.32] again, we obtain that u \in (H 2 loc (\Omega )) 3 . Next, choose a smooth cutoff function \chi that is compactly supported in \Omega and equal to one on \Omega 1 . We find that \| u\| (H 2 (\Omega 1)) 3 \leq \| \chi u\| (H 2 (R 3 )) 3 < \infty .

Thus we obtain that u \in H 2 (\Omega 1 ).

Case II: \bfitgam \bfitJ \not = 0. Here (18) implies that ( 21)

Inserting the above into the Maxwell equations ( 1), we obtain an equation of the form

where

Note that there are boundary constraints for E, including the impedance boundary condition ( 2) and the measurement (3).

To analyze the stability of the inverse problem, suppose that there is another source \J with corresponding electric field \Ẽ and vector internal functional \Q , defined by (18). Let u := E -\Ẽ \in X be a weak solution of the equation (23) \nabla \times

and obey the impedance boundary condition

(\n \times u) \times \n| \partial\Omega = g -\g.

It remains to establish the solvability of ( 22) with boundary condition (2) or the solvability of ( 23) with boundary condition (24). Now ( 22) and ( 23) are similar in form to (10). The difference is that in (10), the term (\omega 2 \varepsi \delta + i\omega \sigma \delta ) has a strictly positive real part and a nonnegative imaginary part, which ensures the existence and uniqueness of the weak solution by standard methods [31]. These sign conditions no longer hold for the term a+ib in ( 22) and ( 23), due to the presence of the elasto-electric constants \gamma \varepsi , \gamma \sigma , \gamma J . Therefore, we divide the discussion into several sub-cases. By Assumption A3, we see that a is either identically zero or bounded away from zero, b is either nonpositive or nonnegative. This observation accounts for the following classification of subcases.

Theorem 4.3. Suppose the assumptions A1--A4 and the hypothesis (3.3) hold. If \gamma J \not = 0, we have the following subcases:

(II.1) If \gamma \varepsi = \gamma \sigma = \gamma J , then the source J cannot be uniquely determined.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 04/11/23 to 132.174.252.171 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy (II.2) If \gamma \varepsi = \gamma J , \gamma \sigma \not = \gamma J , and \Omega \subsete supp \sigma , then the source J is uniquely determined. If, in addition, | b| is strictly bounded away from zero, then we have the following stability estimates. If \gamma \sigma /\gamma J < 1,

and if \gamma \sigma /\gamma J > 1,

(II.3) If \gamma \varepsi = \gamma J , \gamma \sigma \not = \gamma J , and \Omega \not \subsete supp \sigma , then the source J cannot be uniquely determined. (II.4) If \gamma \varepsi \not = \gamma J , then the source J is uniquely determined. If \gamma \varepsi /\gamma J > 1, we have the following stability estimate:

and if \gamma \varepsi /\gamma J < 1,

Here C > 0 is a constant independent of J, \J . Moreover, whenever J is uniquely determined, there are explicit reconstruction procedures.

The proof is presented in the next few subsections.

Remark 4.4. It is generally expected that \gamma J \not = \gamma \varepsi because the former is solely a density effect, and the latter is due to density variation and Brillouin scattering [10]. Thus case (II.4) is more likely to occur in practice.

Case Subcase Uniqueness \gamma J = 0 (I.1) \Omega \subsete supp (\omega 2 \varepsi \gamma \varepsi + i\omega \sigma \gamma \sigma ) Y (I.2) \Omega \not \subsete supp (\omega 2 \varepsi \gamma \varepsi + i\omega \sigma \gamma \sigma ) N \gamma J \not = 0 (II.1) \gamma \varepsi = \gamma \sigma = \gamma J N (II.2) \gamma \varepsi = \gamma J , \gamma \sigma \not = \gamma J , \Omega \subsete supp \sigma , \sigma > 0 Y (II.3) \gamma \varepsi = \gamma J , \gamma \sigma \not = \gamma J , \Omega \not \subsete supp \sigma N (II.4) \gamma \varepsi \not = \gamma J Y Remark 4.5. The stability in cases II.2 and II.4 indeed deteriorates when \gamma \sigma /\gamma J or \gamma \varepsi /\gamma J approaches 1, in which case a or b in (23) approaches 0, and the stability constant in cases II.2 and II.4 blows up.

4.2.1. Subcase (II.1): \bfitgam \bfitvar = \bfitgam \bfitsig = \bfitgam \bfitJ . This subcase corresponds to a \equiv 0 and b \equiv 0. For any \phi \in C \infty c (\Omega ), if E satisfies the (22) and the boundary condition (2), so does E+\nabla \phi . This means, as a result of (18), that sources of the form J \phi := (i\omega \varepsi - \sigma )\nabla \phi are nonradiating. Thus the original source J cannot be uniquely determined. 4.2.2. Subcase (II.2 and II.3): \bfitgam \bfitvar = \bfitgam \bfitJ and \bfitgam \bfitsig \not = \bfitgam \bfitJ . This subcase corresponds to a \equiv 0 and b \not \equiv 0. Note that either b \geq 0 everywhere or b \leq 0 everywhere due to assumption A3. In the following discussion, we will keep a as a placeholder for ease of exposition.

We now take the inner product of ( 23) with u and use the vector identity (\nabla \times A) \cdot B = \nabla \cdot (A \times B) + (\nabla \times B) \cdot A to integrate by parts

For the boundary integral, we apply the vector triple product identity (B \times C) \cdot A = C \cdot (A \times B) to obtain Here u T = (\n \times u) \times \n| \partial\Omega is the tangential trace of u as defined in (8). From (24), we obtain \n \times ( 1 \mu \nabla \times u) = - i\omega \lambda u T . Thus the boundary integrand becomes

Therefore, separating the real and imaginary parts of ( 26) we obtain

To prove uniqueness, we set g = \g. Then Q = \Q and (28) implies b| u| = 0. If \Omega \subsete supp \sigma = supp b, we conclude that u \equiv 0 in \Omega . If \Omega \not \subsete supp \sigma = supp b, there exists an open set D \subsete \Omega \setminu supp b. For any compactly supported smooth function \phi \in C \infty c (D), the choice u := \nabla \phi is a nontrivial solution to (23), (24), proving the nonuniqueness. Now we prove stability assuming that \sigma is strictly positive, which implies that b is bounded away from zero. When b < 0, recall that a = 0, so there exists a constant C > 0, independent of u, such that

where \eta > 0 is an arbitrary constant. If we choose \eta so that C 2\eta < 1, then the term C 2\eta \| u\| (L 2 (\Omega )) 3 can be absorbed into the left-hand side, resulting in the following estimate (with a different constant C):

This result, combined with (21), yields the stability estimate

When a = 0 and b > 0 is bounded away from zero, we can obtain a better stability estimate. In this case, the left-hand side of ( 28) is the sum of two nonnegative terms, then the estimate (29) \leq C\| Q -\Q \| (L 2 (\Omega )) 3 \| u\| H(curl,\Omega ) .

Canceling out \| u\| 2 H(curl,\Omega ) and applying the relation ( 21) yields the stability estimate

4.2.3. Subcase (II.4): \bfitgam \bfitvar \not = \bfitgam \bfitJ . This subcase corresponds to a \not = 0. Note that due to assumption A3, there exists a constant c > 0 such that either a \geq c > 0 everywhere or a \leq - c < 0 everywhere in \Omega .

\bullet If a \leq - c < 0, the identities ( 26), ( 27), ( 28) still hold, hence there exists a constant C > 0 such that

where the second equality comes from (27). Suppose Q = \Q , then u = 0 in \Omega , proving uniqueness. The above inequality also implies, by the Cauchy--Schwartz inequality, that

Canceling factors of \| u\| L 2 (\Omega ) and applying the relation ( 21) yields the stability estimate \| J -\J \| (L 2 (\Omega )) 3 \leq C\| Q -\Q \| (L 2 (\Omega )) 3 .

\bullet If a \geq c > 0, we consider u equipped with the Dirichlet boundary condition (30) \delta g := u T = (\n \times u) \times \n = \g -g.

Since E, \Ẽ \in X, we conclude that \delta g \in H - 1/2 (curl, \partial\Omega ). Thus, there exists a function \scrG \in H(curl, \Omega ), such that \scrG T = \delta g and \| \scrG \| H(curl,\Omega ) \leq C\| \delta g\| H - 1/2 (curl,\partial\Omega ) .

Set \ũ := u -\scrG , then \ũ T = 0 and \ũ solves (31) \biggl( 1 \mu \nabla \times \ũ, \nabla \times \bfitphi \biggr)

where H(curl, \Omega ) 0 is the subspace of H(curl, \Omega ) with zero tangential trace, and \leq C(\| Q -\Q \| (L 2 (\Omega )) 3 + \| g -\g\| (L 2 (\partial\Omega )) 3 ).

Remark 4.6. Note that when applying [31,Theorem 4.17], the entire boundary \partial\Omega is equipped with the homogeneous Dirichlet condition. Here we have extended [31,Theorem 4.17], which requires b \geq 0, but is obviously correct for b \leq 0 when all of \partial\Omega has homogeneous Dirichlet boundary condition and b is not constantly zero.

Finally, whenever J is uniquely determined, it can be reconstructed as follows. Given Q and g, solve the boundary value problem ( 22) and ( 3) to obtain E. Then use the Maxwell equations (1) to recover J.

In this section, we present numerical experiments to test the reconstruction of J in cases (I.1) and (II.4). The code is implemented in Python using the finite element PDE solver NGSolve. 1 Numerical experiments are performed on the domain consisting of an infinite cylinder of radius r = 1cm, discretized with a uniform triangular mesh of 19276 triangles. The Maxwell equations (1) are solved with a third-order N\' ed\' elec element.

We denote by \varepsi 0 and \mu 0 the electric permittivity and the magnetic permeability in vacuum, respectively. In a medium with electric permittivity \varepsi and magnetic permeability \mu , we define \varepsi r := \varepsi \varepsi 0 , \mu r = \mu \mu 0 .

We refer to \varepsi r and \mu r as the relative electric permittivity and the relative magnetic permeability, respectively. Moreover, let c be the light speed in vacuum and define (32) \\sigma = 1 c\varepsi 0 \sigma , \Ĵ = c\mu 0 J, \\omega = \omega c .

Using the relation c = 1/ \surd \varepsi 0 \mu 0 , we can rewrite the Maxwell equations (1) as \nabla \times 1 \mu r \nabla \times E -(\\omega 2 \varepsi r + i\\omega \\sigma )E = i\\omega \Ĵ , together with the impedance boundary condition (2), with impedance \lambda = 1.

The physical parameters are chosen as follows. According to assumption A2, \mu r = 1. The frequency is selected such that \\omega = \pi [cm - 1 ], which corresponds to a frequency f = \omega 2\pi \approx 15GHz. Density plots of \varepsi r and \\sigma are displayed in Figure 1. For \varepsi r , the background value is taken to be 37.2 for blood (see [22]) and there are three regions with smaller values of \varepsi r which are 7.79 (top) for fat, 20.2 (left) for nerve, and 36.4 (right) for muscle. The source \Ĵ is a real vector, whose components are shown in Figure 2.

Auxiliary solutions are needed in the reconstruction to compute the vector internal data (18) from the scalar internal data (17). Such solutions are obtained by solving the equation for j = 1, 2: \nabla \times \nabla \times F j -(\\omega 2 \varepsi r + i\\omega \\sigma )F j = 0 in \Omega , where g k is defined by (33)

\nabla \times E j \times n -\lambda (n \times E j ) \times n on \partial\Omega , with E 1 = (e - iky , 0) and E 2 = (0, e - ikx ). Here the wave number k = \sqrt{} (\\omega 2 \varepsi r + i\\omega \\sigma ), where \varepsi r , \\sigma are taken from the background values corresponding to blood. The rationale for the choice of g j is that when the medium is homogeneous, then E 1 and E 2 are mutually orthogonal plane waves. Clearly, such an orthogonality relation may not hold in practice due to the distortion caused by the inhomogeneity. In the following two cases, we measure the relative L 2 error of reconstructions by \| J reconstructed -J true \| (L 2 (\partial\Omega )) 3 \| J true \| (L 2 (\partial\Omega )) 3 .

In this experiment, the modulation parameters are chosen as \gamma J = 0, and \gamma \varepsi = 0.25, \gamma \sigma = 0.35. The scalar internal data Q is obtained by solving the forward problem (1). Then 0.1\% multiplicative noise is added to the signal. The vector internal data Q is calculated from the auxiliary solutions. The reconstruction is carried out using the procedure described in Proposition 4.1. That is, we solve for E from (18) and then recover J from the Maxwell equations ( 1 -\bigl( (\omega 2 \varepsi + i\omega \sigma )E, \bfitphi \bigr) (L 2 (\Omega )) 3 -i\omega \langle \lambda E T , \bfitphi T \rangle .

Here J is solved under the Galerkin framework by treating the left-hand side i\omega (J, \bfitphi ) (L 2 (\Omega )) 3 as the bilinear form and the right-hand side as the linear form with known E. The reconstructed source is shown in Figure 3.

In this experiment, the modulation parameters are chosen as \gamma J = 0.65, \gamma \varepsi = 0.35, and \gamma \sigma = 0.35. The scalar internal data Q is obtained by solving the forward problem (1). Then 1\% multiplicative noise is added to the signal. The vector internal data Q is found using the auxiliary solutions. The reconstruction is performed according to Proposition 4.3. That is, the boundary value problem (22), ( 3) is solved to obtain E. The Maxwell equations (1) are then used to find J through the similar approach of case (I.1). The reconstructed source \Ĵ rec is shown in Figure 4. Remark 5.1. The reconstructions for case (II.4) are much better than those for case (I.1). This can be explained by the corresponding stability estimates. Proposition 4.1 for case (I.1) requires the H(curl) norm of Q to be bounded, which is very sensitive to noise. In contrast, the stability estimate in Proposition 4.3 for case (II.4) only requires that the L 2 norm of Q be bounded.

In this paper, we proposed a constructive approach to recover the electrical current density from the internal data induced by acoustically modulated boundary measurement of time-harmonic electromagnetic fields. We analyzed the uniqueness and stability of the inverse problem, and revealed how these results depend on values of the elasto-electric constants. Based on the analysis, a computational method was developed, implemented, and validated using simulated experiments. The performance of the numerical reconstruction was observed to agree with theory.

An important ingredient in deriving the uniqueness and stability results is the existence of three linearly independent fields (see Hypothesis 3.3). However, it remains an open question regarding what conditions are sufficient to ensure the existence of these fields. Characterization of existence is theoretically significant for understanding solutions of the Maxwell equations, and is critical for designing reliable computational methods for the electromagnetic inverse source problem. This topic will be explored in future work.