Applications of CGO Solutions on Coupled-Physics Inverse Problems

This paper surveys inverse problems arising in several coupled-physics imaging modalities for both medical and geophysical purposes. These include Photo-acoustic Tomography (PAT), Thermo-acoustic Tomography (TAT), Electro-Seismic Conversion, Transient Elastrography (TE) and Acousto-Electric Tomography (AET). These inverse problems typically consists of multiple inverse steps, each of which corresponds to one of the wave propagations involved. The review focus on those steps known as the inverse problems with internal data, in which the complex geometrical optics (CGO) solutions to the underlying equations turn out to be useful in showing the uniqueness and stability in determining the desired information.


Introduction
Coupled-physics inverse problems arise in various hybrid medical imaging and seismic imaging modalities. Usually two or more different types of wave propagations are involved, subsequently triggered through natural energy conversion. Such physical coupling mechanism overcomes limitations of classical single-measurement based tomography techniques and delivers potentially life-saving diagnostic information with both better contrast and higher resolution.
To be more specific, many traditional single-propagation-based imaging methods suffer from either low contrast or low resolution. An example of low-contrast imaging method is Ultrasound Imaging (UI). UI exhibits high resolution due to its hyperbolic nature, or more plainly saying, its richer time-dependent measurements. Yet the reconstructed sound speed of the material does not distinguish healthy tissues from cancerous tissues very well since both tissues have similar acoustic properties. On the other side of the spectrum, a class of methods, such as Optical Tomography (OT), Electrical Impedance Tomography (EIT) and so on, aim at reconstructing optical/electrical properties of the material. Such properties are more sensitive to the intrinsic physiological properties (oxyand deoxy-hemoglobin, water, lipid, and scatter power), hence provide better contrast in imaging soft tissues. However, due to the diffusive nature of such propagations, when the measurements are made outside the object (non-invasively), sharp features of the material have been "smoothed out", resulting in a low resolution. In mathematical terms, this low resolution phenomena manifests the ill-posedness of the inverse problem of reconstructing diffusive (optical/electrical) coefficients from boundary measurements. The fix offered by multi-wave coupled-physics modalities is roughly speaking to carry the internal information, correlated to the optical/electrical properties, stably to the boundary using another wave propagation, e.g., the sound wave.
To be more illustrative about the idea, we compare EIT and the most popular coupledphysics methods known as Photo-acoustic Tomography (PAT) and Thermo-Acoustic Tomography (TAT). In EIT, imaging is based on recovery of the value of conductivity function γ(x) everywhere in the bounded region Ω modeling the human organ. The measurement is the voltage-to-current (or current-to-voltage) map taken on the boundary ∂Ω. The mathematical inverse problem is the classical Calderón problem [29]: to reconstruct γ from the Dirichlet-to-Neumann (DtN) map Λ γ of the elliptic conductivity equation ∇ · γ∇u = 0 (physically Ohm's law). A lot of work has been done in solving this nonlinear inverse problem (see [87] for a thorough review of the problem). It has been shown that when γ is scalar and satisfies certain regularity conditions, one can expect to reconstruct γ uniquely from Λ γ . However, it is also known (e.g., in [2]) that such problem is ill-posed, which accounts for the low resolution behavior mentioned above. It is well understood that this is due to the smoothing effect of the operator γ → Λ γ .
On the other hand, PAT and TAT are based on the photo-acoustic effect [38,44,92]. When an object (usually animal tissue) is exposed to a short pulse of electromagnetic radiation, a fraction of the radiation is absorbed by the medium, resulting in a thermal expansion. This expansion then emits acoustic waves, which propagate to the boundary of the domain. This physical coupling between the absorbed radiation and the emitted acoustic waves is called the photo-acoustic effect. What distinguishes PAT and TAT is the frequency of the radiation used to illuminate: in PAT, high frequency radiation such as near-infrared with sub-µm wavelength is used; while in TAT, low frequency microwave with wavelengths comparable to 1m is used [61]. The inverse imaging process consists of two steps: first, to reconstruct the absorbed radiation inside the medium from acoustic signals measured on the boundary; second, to reconstruct the optical property (exhibiting better contrast) of the tissue from this internal information obtained in the first step. The first step is shown to be a stable inverse source problem for the wave equation (see section 2 for more details and references). The second step is an inverse problem with internal data, which is richer compared to the boundary data used in EIT.
From above examples, we observe that coupled-physics imaging methods involve solving multiple steps of inverse problems, each of which has to be well-posed (stable) to give an overall high resolution in imaging. It is a major feature of these methods that usually the last step is to solve an inverse problem with internal data obtained from previous steps. Such inverse problems are what we focus on in this review. The coupled-physics methods we are going to consider include PAT, TAT, Electro-Seismic Conversion (ESC), Transient Elastrography (TE) and Acousto-Electric Tomography (AET). The internal data, which we denote by H through the paper, obtained in these modalities are usually polynomials of the solutions u to the underlying equations of radiation, or their derivatives. In another word, the information of the coefficients of interest are hidden in the form of the underlying equations and the internal data of solutions to them. Interpreted this way, it is not surprising that a particular type of solutions play a major role in solving this type of inverse problems. Here we explore one of such solutions known as complex geometrical optics (CGO) solutions.
CGO solutions were first introduced in [81] to the conductivity operator ∇ · γ∇, to solve the nonlinear Calderón problem for EIT. Since then, such solutions were successfully constructed to several other equations such as the elasticity equation, Maxwell's equations and so on to solve various inverse problems. (see [31,37,67,68,69,70,88,89]). The strategy of construction usually starts with the reduction to a Schrödinger equation (∆ + q)u = 0 in three or higher dimensions. In case of systems, Maxwell's equations can be reduced to a matrix Schrödinger equation of the same type in a less straightforward fashion (see [70]) while elasticity equations can be reduced to a Schrödinger equation with external Yang-Mills potentials (see [51]). Given a complex vector ζ ∈ C n such that ζ · ζ = 0, a CGO solution to (∆ + q)u = 0 is of the form where ψ ζ satisfies certain decaying property as |ζ| → ∞. In [81], this is done by solving an equation about ψ ζ with leading operator ∆ + 2iζ · ∇ whose inverse is the integral operator with Faddeev kernel, that is, CGO solutions with nonlinear complex phase are also available (see [42,52]) using Carleman estimate. The construction can be manipulated to give solutions vanishing on part of the boundary, which is useful for partial data problems, that is, the inverse problem with measurements taken only on part of the boundary [27,52]. For a thorough review on CGO solutions, we refer the reader to the review [87]. Recently, breakthrough has been made on constructing CGO solutions to equations with less regular parameters, using an averaging technique, see [32,46,47]. In this survey, our emphasis is on applications of CGO solutions in solving inverse problems with internal data arising in various imaging modalities introduced above. It is not an easy task to categorize the modalities based on the usage of CGOs. Instead, we contribute roughly each section to address the application in one modality. The following Table 1 summarizes the modalities, along with the brief information of the underlying equations, formats of the available data, types of CGO solutions used and the obtained results. Roughly speaking, when solving the listed inverse problems, we often find ourselves in two scenarios after certain reduction, where CGO solutions turn out to be useful, • The first scenario seen in QPAT, ESC and TE (Section 2-4) is when the parameters of interest are associated to the gradient of u, where u is the solution to the underlying equation. Simple algebra reduces the problem to solving a transport equation for the unknown parameters. The solvability of the transport equation relies on the density of vector fields β, which is written in terms of ∇u. Another example is in AET (Section 5) where the internal data itself is a functional of ∇u. The overall strategy here is to use CGO solutions so ∇u ∼ iζe iζ·x for |ζ| sufficiently large. With well-chosen ζ's, one obtains sufficiently many linearly independent vectors β so the parameters can be recovered by solving the transport equations. • In the second case, for example for the system model of QTAT (Section 6), (also see [14,16] for anisotropic conductivities in AET and UMEIT, ), the linearized inverse problem is considered. As a result, it is reduced to solving a boundary value problem for a system of sometimes overdetermined (pseudo-) differential equations. By the Douglis-Nirenberg theory, the ellipticity of the boundary value problem provides the stability estimate for the linearized inverse problem. In particular, for the principal symbol of the (pseudo-) differential operator to be Uniqueness and stability in determining (γ, σ) (see [22]). data: u| Γ → σu| Ω (partial boundary illuminations).

Section 4: Transient Elastography
Elasticity system: to the Schrödinger equation with external Yang-Mills potentials.
to the conductivity equation.
non-degenerate, one needs to show again there are sufficient linear independent vector fields. This can be achieved by plugging well-chosen CGO solutions. These are the two commonly seen approaches to apply CGO solutions to inverse problems with internal data. There are other scenarios where CGO solutions can be used on a case-by-case basis suggested by the special structure of underlying equations. For example, CGO is used to give the Fourier transform of the internal data in the first step of AET, and to form a contraction of the unknown parameters in the scalar case of QTAT.
We would like to point out that this paper by no means has covered every aspect or method used in tackling inverse problems with similar types of internal information of the solutions. For example, we do not attempt to include the review on reconstruction methods. However, we would like to mention a local reconstruction scheme introduced in [11,13,23] where a local linear independence condition needs to be satisfied in order to guarantee the reconstruction. In [23], harmonic polynomials are used to show the condition locally. This condition can be extended for global reconstruction by Runge approximation/unique continuation in some situations, e.g., in [12,23] by applying CGO solutions. Such method is versatile in reconstructing anisotropic tensor-valued parameters. Other features of the inverse problems with internal data can also be found in another extensive review [5].

Quantitative Photo-acoustic Tomography
In the present section and the two sections following this, we present methods sharing a general strategy of the CGO application, which was first developed by Bal and Uhlmann in [22]. The idea is to reduce the equations modeling these problems to the Schrödinger's equation or Maxwell's equations, and then inserting sufficiently many internal functions H to obtain a transport equation in one of the unknown parameters. The uniqueness and stability of the recovery of this unknown finally rely on the uniqueness and stability of the solution to the transport equation. Then CGO solutions are used to show the solvability of the transport equation. This idea is initiated by [22] and then further expanded and utilized in the analysis of many coupled physics inverse problems. Here we review some results in photoacoustic tomography [22,34]. The results for electro-seismic conversion [35] is presented in Section 3, and that of transient elastography [60] is presented in Section 4.
In both PAT and TAT, the first step of the recovery procedure is to reconstruct the absorbed radiation from the boundary measurements of the acoustic waves. This step is typically modelled as an inverse source problem for the acoustic wave equation. This problem has been extensively studied in the mathematical literature, see [1,4,43,48,50,55,58,59,71,78,79,80]. We will assume in this section that the absorbed radiation H(x) has been recovered and concentrate on the second step.
The second step of PAT and TAT is modeled by different equations due to the difference of the radiation used. In PAT, the high frequency radiation (near-infrared laser pulses) is modeled by the diffusion equation, while in TAT the low frequency microwave is modeled by Maxwell's equations. The second step in PAT and TAT are usually referred to as Quantitative Photo-acoustic Tomography (QPAT) and Quantitative Thermo-acoustic Tomography (QTAT), respectively. We consider the QPAT model in this section. A detailed discussion of the QTAT model can be found in Section 6.
In QPAT, radiation propagation is modeled by the following diffusion equation Here γ(x) is the diffusion coefficient, σ(x) is the absorption coefficient, f is the illumination on the boundary. The measurement is the absorbed radiation which we assume to be known after solving the first step. The objective of QPAT is to reconstruct (σ, Γ) from the knowledge of H(x) obtained for a given number of illuminations f . For references on QPAT, see e.g., [19,20,22,38,39,74,93].
2.1. Full Data. In this part we review some full data results on the QPAT model due to Bal and Uhlmann [22]. The QPAT model is described by the diffusion equation with internal data. It was observed in [22] that the inverse problem of the diffusion equation with internal data can be reduced to the one of the Schrödinger's equation by the following transform. Define and the internal measurement becomes H(x) = σu = µv. The goal here is to recover q and µ from H(x). Then by the definition of q, one can solve Therefore, it remains to consider the inverse problem for the Schrödinger's equation with internal data. For this purpose, the authors of [22] construct a class of CGO solutions with higher regularity.

Smoother CGO solutions.
We begin by reviewing the main ingredients in the construction of L 2 -CGO solutions initiated by Sylvester and Uhlmann [81], based on which smoother CGO solutions will be constructed for later applications.
Let ζ ∈ C n be a complex vector with ζ · ζ = 0. Define the space L 2 δ (δ ∈ R) to be the completion of C ∞ c (R n ) with respect to the norm Notice that the function u ζ := e iζ·x (1 + ψ ζ (x)) is a solution of the Schrödinger equation (∆ + q)u = 0 if and only if ψ ζ solves It remains to find solutions ψ ζ of (2). To this end, denote Faddeev's Green kernel by where F is the Fourier transform. It is shown [81] that for |ζ| large one has Therefore if |ζ| is sufficiently large, the equation (2) has a unique solution ψ ζ by the fixed point theorem. This is the following result.
). Let q ∈ L ∞ (Ω) and −1 < δ < 0. For any ζ ∈ C n with ζ · ζ = 0, there exists a unique solution to the Schrödinger equation (∆ + q)u = 0 of the form with ψ ζ ∈ L 2 δ . Moreover, ψ ζ satisfies the estimate Smoother CGO solutions in higher order Sobolev spaces can be constructed with minor modifications as follows [21,22]. Introduce the weighted Sobolev space H s δ (s ≥ 0) as the completion of C ∞ c (R n ) with respect to the norm Here (I − ∆) s 2 is a pseudodifferential operator whose symbol is (1 + |ξ| 2 ) s 2 . Noticing that the two constant coefficient operators (∆ + 2iζ · ∇) and (I − ∆) s 2 commute, one obtains Finally, by a Neumann series argument, a solution ψ ζ to (2) is obtained and satisfies when s = n 2 + k + for some k positive integer and > 0. Restricting to the bounded domain Ω where q is compactly supported, and applying Sobolev embedding theorem yield For the consideration of the inverse problem, define the set of admissible parameters (Ω) ≤ M < ∞, and µ is bounded away from 0 .

Uniqueness and stability.
Suppose ∂Ω is of class C k+1 , g j ∈ C k,α (∂Ω; C), j = 1, 2, with α > 1 2 , and (µ, q) ∈ P. Then the following problem admits a unique solution u j ∈ C k+1 (Ω). From this we verify that u 1 ∆u 2 − u 2 ∆u 1 = 0. Taking into consideration u j = Hj µ , one has The unique solvability of this transport equation in µ depends heavily on the behavior of the computable vector field √ 2 |ζ| and α 1 · α 2 = 0, one can compute after some basic algebra that using the estimate (2). This expression indicates that, as |ζ| becomes larger, the direction of β (the imaginary part of β) becomes more consistent with the vector field α2 |α2| . Moreover, it is non-vanishing, implying that each point in Ω is connected to a point on ∂Ω by an integral curve of β. Then the transport equation (10) is uniquely solvable and one can take the imaginary part of it to solve uniquely for µ.
It is easy to see that u ζ = u ζ . Then the analysis above shows that by taking two realvalued boundary illuminations such as { u ζ | ∂Ω , u ζ | ∂Ω } one is able to uniquely determine µ. However, for this to provide a reconstruction scheme, one needs to know the boundary values of the CGO solutions, which are not available. For this to be partly resolved, it is shown in [22] that in fact this set of boundary illuminations can be made larger, that is, close enough to { u ζ | ∂Ω , u ζ | ∂Ω }, using the elliptic regularity theory.
Converting the above unique determination result back to the diffusion equation case, Bal and Uhlmann obtained in [22] the following theorem. Let k ≥ 1 and set be the internal data for the coefficients (γ, σ) and (γ,σ), respectively, and with boundary conditions g := {g 1 , g 2 }. Then there is an open set of illuminations g ∈ (C 1,α (∂Ω)) 2 for some α > 1 2 such that if H =H, then (γ, σ) = (γ,σ). By taking a closer look at the behavior of the integral curve of the vector field β, a Lipschitz type stability is also derived in [22] under certain assumptions on the geometry of Ω from two real-valued measurements. Since the proof is not a direct application of CGO solutions, we refer interested readers to [22] for more details.
With more measurements available, the above idea is further developed in [22] to derive a stability result. More specifically, by imposing on the boundary 2n real-valued (n is the spatial dimension) illuminations that are taken to be the real and imaginary parts of n well-chosen CGO solutions of the diffusion equation, we obtain n internal complex measurements H c := (H c 1 , . . . , H c n ). These internal functions can be used in the same manner as (11) to generate a collection of vector fields (β 1 , . . . , β n ) whose imaginary parts approximately form a basis of R n . As a result, we obtain a system of transport equations, leading to (12) ∇µ Then the uniqueness of µ follows from the uniqueness of the solution of the system (12), and the stability of this reconstruction of µ follows from the continuous dependence of the solution µ on the coefficient Γ(x) as well as (13). Converting back to the diffusion equation yields Theorem 2.3 ( [22]). Under the above notations. Assume that (γ, σ) and (γ,σ) are in M with γ| ∂Ω =γ| ∂Ω . Let H = (H 1 , . . . , H 2n ) andH = (H 1 , . . . ,H 2n ) be the real-valued internal data for the coefficients (γ, σ) and (γ,σ) respectively. Then there is an open set of illuminations g ∈ (C k,α (∂Ω)) 2n for some α > 1 2 and a constant C > 0 such that 2.2. Partial data. Under certain circumstances, imposing radiation on the whole ∂Ω may be either too costly or impossible, thus it is necessary to consider partial data problems. In QPAT, the partial data problem is considered in [34]. The key idea is to apply a special type of CGO solutions that vanish on part of the boundary. These solutions were first constructed by Kenig, Sjöstrand and Uhlmann in [52], which we briefly explain in the following.
To construct such partial data CGO solutions, the linear phase function ρ · x (ρ = ζ) is replaced by a so-called limiting Carleman weight to allow more freedom. A limiting Carleman weight ϕ is a real-valued C ∞ function on Ω such that ∇ϕ is non-vanishing and the following relation holds: ϕ ∇ϕ, ∇ϕ + ϕ ξ, ξ = 0 whenever |ξ| 2 = |∇ϕ| 2 and ∇ϕ · ξ = 0, where ϕ is the Hessian matrix of ϕ. This relation can be viewed as a generalization of the relation ζ · ζ = 0, thus the linear phase ρ · x is a limiting Carleman weight. Another example, which will be used later, is the function log |x − x 0 | where x 0 ∈ R n is a fixed point outside of the convex hull of Ω. It also follows from the definition that if ϕ is a limiting Carleman weight, so is −ϕ.
Taking ϕ(x) = log |x − x 0 | with x 0 outside of the convex hull of Ω, one divides ∂Ω into two parts: the front side and the back side where ν is the unit outer normal vector field on ∂Ω. Let Γ be a neighborhood of ∂ + Ω in ∂Ω, and Γ − = ∂Ω\Γ.

Proposition 2.4 ([52]
). Let q ∈ L ∞ (Ω). There exists a real function ψ ∈ C ∞ (Ω) with |∇ϕ| = |∇ψ| and ∇ϕ · ∇ψ = 0 such that (∆ + q)u = 0 admits solutions of the form Here One may also increase the regularity to obtain smoother CGO solutions [34]. Then the ideas in the previous full data section can be adapted to obtain uniqueness and stability for partial data problems. Let (µ, q) ∈ P. For any two solutions of (∆ + q)u j = 0 (j = 1, 2), (10) is still valid. Instead of (5), we plug in the CGO solutions from Proposition 2.4 Then the vector field β satisfies Notice that the vector x0−x |x0−x| 2 points to x 0 , hence the integral curves of β hit ∂Ω near the front side ∂ + Ω. This along with the fact that u j | Γ− = 0 allow to make measurements only near ∂ + Ω to have a unique solution µ for the transport equation. Notice that to assemble the two CGO solutions, one needs 4 real-valued boundary illuminations. Therefore, one has Similar to the sequence of results in [22], with more measurements we have are the corresponding internal data.

Electro-seismic Conversion
Seismo-Electric (SE) and Electro-Seismic (ES) conversions are phenomenon occurred in fluid-saturated porous media. These conversions couple electromagnetic waves and elastic waves through the electro-kinetic effect. SE conversion employs seismic source to generate electromagnetic waves, while ES conversion emits electromagnetic waves to excite elastic waves. These conversions have been applied in oil prospection as well as other geophysics studies. Detailed description of the physical mechanism underlying these conversions can be found, for instance, in [35]. Theoretical and experimental results on SE conversion have been obtained in, e.g. [28,64,65,83,85]. However we concentrate on ES conversion in this section The governing equations of ES conversion are derived by Pride [72] based on Biot's theory on the elastic wave propagation in porous media [25,26]. Pride's equations are analyzed in [75,90,91] and tested in [49,73,82,84]. The coupled-physics inverse problem to be considered below is the inverse problem of the linearized electro-seismic conversion. It consists of two steps: the first step is to inverse Biot's equations [72] to recover the internal potential from boundary measurements, and the second step is the inversion of Maxwell's equation with internal measurement. The first step has been investigated in [72] and an approximation method is proposed there. Here we review the results in [35] where it is assumed that the first step has been successfully implemented and investigates only the second step.
Let Ω ⊂ R 3 be an open, bounded and connected domain with C 2 boundary ∂Ω. The propagation of the electric fields in ES conversion is modeled by the following Maxwell's equations when there is no source current. (14) ∇ Here ω > 0 denotes a fixed seismic wave frequency, the constant µ 0 > 0 denotes the magnetic permeability, σ = σ(x) the conductivity, ε = ε(x) the relative permittivity, E the electric field and H the magnetic field. In this section it is assumed that ε(x) = ε 0 and σ(x) = 0 for all x outside some sufficiently large ball containing Ω, where ε 0 > 0 is the dielectric constant. The measurement, which is obtained from the first step, is the internal potential given by H := LE in Ω, where L = L(x) is the coupling coefficient. The controllable boundary illumination is the tangential boundary electric field where ν denotes the unit outer normal on ∂Ω. The question is similar to the one in QPAT: by choosing boundary illuminations G, can we uniquely determine the pair (L, σ) from the internal data H? For small ι > 0, it will be shown that this is possible for the set of coefficients where the wave number κ > 0 and the refractive index n(x) are given by The authors of [35] proved the uniqueness and stability to the reconstruction of the pair (L, σ) under appropriate assumptions. The key ingredient is the construction of special CGO solutions. With these solutions it is possible to derive again a transport equation in L, from which the uniqueness and stability results follow. The solutions used in [35] take the form where ζ ∈ C 3 \R 3 , η ∈ C 3 , are constant vectors satisfying (16) ζ · ζ = κ 2 , ζ · η = 0.
Thus β has approximately constant directions for small h, and their integral curves connect every internal point to ∂Ω. This ensures that the above transport equation admits a unique solution. By a boundary perturbation argument based on the regularity theory of Maxwell's equations, it is obtained in [35] that Based on the above procedure, stability results can be established as well, see [35] for more details.
The results reviewed in this section are obtained under the assumption that the magnetic permeability µ 0 is constant. Recently in [33] this assumption has been relaxed to functions by converting the system (14) to a matrix Schrödinger's equation and utilizing the CGO solutions constructed in [70], interested readers are referred to [33] for more details.

Transient Elastography
In this section, we consider a hybrid inverse problem involving elastic measurements called Transient Elastography (TE) which enables detection and characterization of tissue abnormalities.
TE is a non-invasive tool for measuring liver stiffness. The device creases highresolution shear stiffness images of human tissue for diagnostic purposes. Shear stiffness is targeted because shear wave speed is larger in abnormal tissue than in normal tissue. In the experiment tissue initially is excited with pulse at the boundary. This pulse creates the shear wave passing through the liver tissue. Then the tissue displacement is measured using the ultrasound. The displacement is related to the tissue stiffness because the soft tissue has larger deformation than stiff tissue. When we have tissue displacement, we want to reconstruct shear modulus µ and the parameter λ. See [63] and references there for more details.
This modality, TE, also includes two steps. The first step is to solve an inverse scattering problem for a time-dependent wave equation of ultrasound, which is a high resolution step as PAT and TAT. The second step is to recover the Lamé parameters from the knowledge of the tissue displacement obtained from the first step. In this section, we will again focus on the second step which is to quantitatively reconstruct coefficients that display high contrast from the internal data of tissue displacement.
The mathematical problem is formulated as follows. Let Ω ⊂ R n , n = 2, 3, be an open bounded domain with smooth boundary. Let u = (u 1 , . . . , u n ) T be the displacement satisfying the linear isotropic elasticity system (19) ∇ · (λ(∇ · u)I + 2S(∇u)µ) + k 2 u = 0 in Ω, u = g on ∂Ω, where S(A) = (A + A T )/2 denotes the symmetric part of the matrix A. Here (λ, µ) are Lamé parameters and k ∈ R is the frequency. For the forward problem of elasticity system, we refer the readers to [36]. The time-harmonic scalar model for TE was showed by Bal and Uhlmann [23] where the reconstruction of coefficients in scalar second-order elliptic equations was also studied. The time-harmonic Lamé system from internal measurements was also considered in [8]. The reconstruction algorithms for fully anisotropic elasticity tensors from knowledge of a finite number of internal data were derived in [18]. In this section, we concern the linear isotropic elasticity setting. Assume that k is not an eigenvalue of the elasticity system. We impose J ∈ Z + boundary displacements g (j) (1 ≤ j ≤ J), then the set of internal functions obtained by the first step in TE is given by where u (j) denotes the solution to (19) with boundary condition u (j) | ∂Ω = g (j) . Notice that the internal functions for QPAT, Electro-seismic and TE are all linear functionals of solutions to underlying PDE's. In order to recover the Lamé parameters from the internal functions H, the strategy is similar: we reduce the inverse problem to solving a transport equation for µ or λ.
To be more specific, since Ω is bounded, we pick a ball B R for R > 0, such that Ω ⊂ B R and extend λ and µ to R n by preserving their smoothness and also supp(λ), supp(µ) ⊂ B R . We will use the reduced system derived by Ikehata [51]. This reduction was also mentioned in [86]. Then the elasticity equations (19) can be transformed to the following equations. Suppose (w, f ) T satisfy Here V 0 contains the third derivative of µ and The solution to the elasticity system (19) is then given by Here ∇ 2 g denotes the Hessian matrix (∂ 2 g/∂x i ∂x j ) ij . The CGO solutions to (20) are presented in the following lemma.
and B(x) are (n + 1) × (n + 1) matrices. Then for τ sufficiently large, there exists a solution to (22) of the form and det C 0 is never zero, p(z) is an arbitrary polynomial in complex variables z, and The above lemma provides the CGO solutions (w ζ , f ζ ) T to (20). Substituting it into (21), we obtain the CGO solutions u ζ to (19).
Meanwhile, let us consider the two dimensional case and the reconstruction of µ here (the three dimensional case and the reconstruction of λ are similar), in which case we can deduce the following equation from (19), where u is a solution to (19) and with a = ∂ 2 u 1 + ∂ 1 u 2 , b = ∂ 1 u 1 − ∂ 2 u 2 , and u * = (u 1 + u 2 ). The component of u and u can be computed from the known internal data u.
Obverse that the vector u has three different entries. The rest of the proof is to choose enough proper CGO solutions u ζ on some subdomain of Ω such that the first term u · F in (23) vanishes. Then we obtain a transport equation only about µ of the form (24) ∇µ where Γ and Ψ are vector-valued functions that can be computed from the known internal data. The uniqueness and stability of µ follow immediately from the similar argument as before (also see [22]). More specifically, in [60], we choose three suitable CGO solutions u (j) for j = 0, 1, 2 to (19). Now for = , , * , set where χ j (x) is some nonzero function (see [60]). Since λ, µ are real-valued functions, we have Based on the choice of CGO solutions in [60], it is shown that u (1) , u (2) , u (3) are three linearly independent vectors at every x in some subdomain Ω 0 of Ω. Then, for each l = 1, 2 there exist three functions Θ l 1 , Θ l 2 , and Θ l 3 such that Then for l = 1, 2, multiplying (26) by Θ l j and summing over j and equation (25), we have which can be further rewritten in the form Here again the choice of CGO solutions guarantees that β 1 (x) and β 2 (x) are linearly independent for every x ∈ Ω 0 , which implies (24). Therefore, µ can be uniquely and stably reconstructed in Ω 0 . Let us denote It is shown in [60] that the local reconstruction of λ in two dimensional case requires additional four CGO solutions, hence one needs seven CGO solutions to obtain the reconstruction of λ and µ together. The main result is stated in the following theorem. Suppose that the Lamé parameters (λ, µ) and (λ,μ) ∈ P satisfy µ| ∂Ω =μ| ∂Ω . Let u (j) andũ (j) be the solutions of the elasticity system with boundary data g (j) for parameters (λ, µ) and (λ,μ), respectively, and H = (u (j) ) 1≤j≤J ,H = (ũ (j) ) 1≤j≤J be the corresponding internal data for some integer J ≥ 3n + 1 .
Then there is an open set of the boundary data (g (j) ) 1≤j≤J such that H =H implies (λ, µ) = (λ,μ) in Ω. Moreover, we have the stability estimate

Acousto-electric Tomography
In Acousto-Electric Tomography (AET), also known as the ultrasound modulated electrical impedance tomography (UMEIT), acoustic waves sent to a conductive object cause a slight change in the conductivity of the object. This interaction between the electric and acoustic waves is called the acousto-electric effect. Even though the change in conductivity due to acousto-electric effect is small, it can be observed by making electrical boundary measurements [94]. This observable difference in the boundary measurement readings is the data for the inverse problem of AET, and the goal is to recover the background internal conductivity from these measurements [3,6,9,16,30,45,54,56,57,94], Similar to the other coupled physics methods a two-step approach is used in AET. First step is the reconstruction of an internal data from the actual boundary measurements. More precisely, before the ultrasound modulation, voltage-to-current (or alternatively current-to-voltage) measurements are done on the boundary as in the case of EIT. Then the same measurements are done while the object is being scanned by ultrasound waves. The difference in these two current responses is due to the change in the conductivity and it is used to obtain internal information about the medium. The second step is then to obtain the background conductivity from this internal data obtained in the previous step. Therefore the model here is different from that of EIT, and it is shown that the inverse problem of AET has better stability compared to EIT.
The change in the conductivity map due to the ultrasound waves is modelled by where γ(x) is the background conductivity map, and m(x) is the modulation which is assumed to be a known smooth map that depends on the acoustic signal strength [3,6,57].
Let Ω be an open bounded connected domain in R n with smooth boundary ∂Ω and γ be the unknown isotropic background conductivity map which is assumed to be known on the boundary ∂Ω. The conductivity equation is given by (27) ∇ · (γ∇u) = 0 in Ω, u| ∂Ω = f, for a given boundary potential f ∈ H 1/2 (∂Ω). The Dirichlet-to-Neumann (DtN) map Λ γ : f → γ∂ ν u| ∂Ω , also known as the voltage-to-current map, can then be defined by the following quadratic form where u g is an extension of g to Ω. Then the difference of DtN maps due to ultrasound modulation m, evaluated at boundary potential f , is The goal of AET is to recover γ(x) from M f . First, we consider the CGO solutions of the conductivity equation ∇ · γ∇u = 0, where ζ ∈ C n satisfies ζ · ζ = 0 and ψ ζ ∈ H s satisfies (7). The existence and construction of these solutions are given in Section 2. The conductivity equation can be reduced to the Schrödinger's equation (1) using the same transform described in Section 2.1 with q = − ∆ √ γ √ γ . Then it can be seen that the gradient of these special solutions are of the form Then, we set m(x) = e (−iζ−ik)·x , where k ∈ R n is fixed. Here the map M f is extended to complex valued functions. Let u, u m ∈ H 1 (Ω) be the solutions to (27) for the conductivities γ and γ m , respectively. We obtain that for g = u ζ | ∂Ω , where |r(g)| is bounded from above and the bound is independent of ζ. Then it is possible to recover an internal data of the form √ γ∇u in Ω. This roughly follows from the fact that √ γ∇u ζ is almost flat up to a known function thanks to the behavior of CGO functions that ψ ζ H s decays as |ζ| increases, as mentioned above. where u is the solution to the conductivity equation with u| ∂Ω = f . Let ζ j = −iτ(e j + ie n ) for 1 j < n, and ζ n = −iτ(e n + ie 1 ) where {e j } n j=1 denotes the standard basis of R n . Then √ γ(x)∇u(x) can be recovered from where A is a known invertible matrix and F denotes the Fourier transform.
The calculation of s(k, ζ) requires the knowledge of g = u ζ | ∂Ω . We refer to [53] and [66] where these traces of CGO solutions can be recovered from the DtN map Λ γ . Next, assuming the knowledge of this internal data √ γ(x)∇u(x) in Ω, AET is reduced to an inverse problem of reconstructing conductivity γ from this internal data, which can be solved stably.
Theorem 5.2. [54] Let Ω be a bounded domain in R n with smooth boundary. Given two conductivities γ andγ, and two boundary potentials f andf ∈ H 1/2 (∂Ω) satisfying for some positive constant M , denote the corresponding internal data H = √ γ∇u and H = √γ ∇ũ, respectively. Then we have the following stability estimates: (2) For some l 1, suppose that √ γ, √γ are bounded by M in C l,α (Ω). Similarly, suppose f,f are bounded by M in C l,α (∂Ω). If |H(x)| > δ > 0, then there exists a constant C = C(M, Ω) such that It is also possible to extend the above theorem to the case with a more general internal data of the form H = γ s ∇u for s ∈ R. In [54] this is achieved by using similar techniques as in [22].
Another internal data used in AET is power densities of the form γ|∇u| 2 and the second step is now to recover the conductivity from these power densities. CGO solutions are also useful in this case. In [9] a similar inverse problem is studied where the power densities are considered as the internal data. The authors give stability estimates and reconstruction of γ using a data set of m illuminations (where m is equal to the spatial dimension n for even n and n + 1 otherwise). The corresponding set of internal data is assumed to satisfy (31) det(H 1 , . . . , H m ) c 0 > 0 in Ω. The set of solutions that satisfy the above condition is a key tool in their approach and it is constructed by utilizing CGO solutions in the Lemma 3.3 of [9]. The existence of the set of illuminations whose internal data satisfies (32) is obtained by employing the fact that the gradient of CGO solutions can be chosen to be flat up to known functions and terms that can be made negligible. In [16] the authors generalize their results to multi-dimensional settings, using an internal data of the form γ α ∇u j ·∇u k , where α ∈ R, (n − 2)α + 1 = 0 and by constructing CGO solutions to satisfy an analog of (32). In a series of work [10,11,12,14,15,17], these results are extended to the case of anisotropic conductivities and also for the Maxwell's equations.

Quantitative Thermo-Acoustic Tomography
In this section, we consider the second step of the inversion for thermo-acoustic tomography (TAT). Recall that TAT is another coupled-physics imaging method implementing the photo-acoustic effect of electromagnetic radiation. But different from PAT, lower frequency electromagnetic "illuminations" are used, in order to achieve deeper penetration in the tissue. As a result, the underlying equation is either modeled by an approximate scalar Helmholtz type equation or the full Maxwell's system. We explore the applications of CGO type solutions in solving inverse problems for both models, that is, to recover the high contrast conductivity (otherwise known as the absorption coefficient) from the absorbed energy by the tissue. A difference of these inverse problems from the previous ones is that the internal measurements are no longer linear functionals of the solutions but the norm squares of them. Consequently, we can see that the CGO solutions are used in very different fashions in answering uniqueness and stability questions for the two models. In particular, the general uniqueness question is still open for the system case.
6.1. QTAT for scalar equations. In this section we consider the Helmholtz model of radiation propagation, governed by the equation where q(x) = k 2 +ikσ(x) with k denoting the wavenumber, σ(x) the conductivity we wish to reconstruct, and g(x) the boundary illumination. The amount of absorbed radiation by the underlying medium is given by A slightly more detailed derivation of this scalar approximation of electromagnetic propagation can be found in [21].
In [21], the CGO solutions to (33) are used to show that the internal data H leads to a functional of σ that admits a unique fixed point. To be more specific, we take σ ∈ H s (Ω) for s > n/2, so q ∈ H s (Ω) and we have from Section 2 the smoother CGO solutions where ψ ζ solves and satisfies the following estimate, by (7), Moreover, by refining the Neumann series argument step in the construction of these solutions, it is shown in [21] that ψ ζ is actually a contraction of σ, that is, for |ζ| large enough, whereψ ζ is the solution to (35) with σ replaced byσ. Substituting the CGO solutions into the internal data, we obtain Using estimates (36) and (37), and the fact that H s (Ω) is an algebra for s > n/2, we have shown that for |ζ| sufficiently large, there exists a constant c < 1 such that that is, H[σ] is a contraction in H s (Ω). To summarize, we have Theorem 6.1.
[21] Let ζ ∈ C n be such that |ζ| is sufficiently large and ζ · ζ = 0. Let σ be a function in Then there exists a boundary illumination g = u| ∂Ω where u is the CGO solution (34) such that the corresponding internal measurement H(x) uniquely determines σ. Moreover, we have the reconstruction algorithm where H[σ] is defined by (39). LetH be the counterpart replacing σ byσ ∈ M. Then there exists a constant C independent of σ andσ such that The remainder part of the section is original. By an estimate presented in [76], we are able to reduce the regularity assumption on the coefficient σ to L ∞ in the uniqueness and stability result above. Lemma 6.2. [76, Lemma 1] Let G ζ be the Faddeev's Green kernel defined in (3). There exists constant c > 0 depending on σ such that for any f ∈ L 2 (Ω) and for |ζ| > 1, we have that where l < 1 for n = 2 and l < 1 2 for n = 3. Theorem 6.3. Let σ,σ be a bounded measurable function satisfying for some constant M > 0. Assume that ζ ∈ C n satisfies that ζ · ζ = 0 and |ζ| is sufficiently large. Then there exists a boundary illumination g such that the corresponding internal measurement H(x) uniquely determines σ. Moreover, we have the reconstruction algorithm where H is defined by (39). LetH be the counterpart replacing σ byσ, also measurable and satisfying (41). Then there exists a constant C independent of σ andσ such that holds true.
Proof. In the following, C denotes the generic positive constant. Note that, by using the estimate of the operator (∆ + 2iζ · ∇) −1 as in Proposition 2.1, we have that ψ ζ is L 2 -bounded for sufficiently large |ζ|. Applying Lemma 6.2 to ψ ζ , with the aid of the L 2 -boundedness of ψ ζ , we obtain that where C only depends on M and Ω for |ζ| large. Similarly, we can obtain the boundedness of ψζ,ψ ζ andψζ.
Let v = ψ ζ −ψ ζ . It follows that Applying Lemma 6.2 to v, we obtain that The above inequality implies that for some constant C only depends on M and Ω when |ζ| is large enough. The same estimate applies to ψζ −ψζ. It follows that By (43), (44) and (45), we conclude that for some constant C only depends on k, M and Ω. This yields that H is a contraction when one chooses a sufficiently large |ζ|. Recall that The stability estimate (42) follows from (47).
We formulate the stability estimate in L ∞ , which is an algebra as H s . Let us remark however that stability estimate in L 2 as could also be obtained with few adjustments of the proof. Actually, with the L 2 version of (44), (45) and (47), the above estimate follows.
6.2. QTAT for Maxwell system. The above scalar model is an approximation of the full Maxwell system model of electromagnetic radiation. Consider time harmonic electromagnetic waves satisfying with boundary illumination given in terms of the tangential electric field where ν is the unit outer normal to ∂Ω. The amount of absorbed radiation by the underlying tissue is given by The quantitative step of TAT concerns the reconstruction of (n, σ) from knowledge of {H j (x) = H gj (x)} 1≤j≤J obtained from the first step by probing the medium with J illuminations {g j } 1≤j≤J . It was shown in [21] that with the refractive index n(x) being a constant, the conductivity σ(x) can be uniquely and stably reconstructed from a single (well-chosen) internal measurement provided that σ is sufficiently small (compared to k).
Here we review the results in [24] for the general case. The fixed point type analysis above is no longer available due to lack of a contraction estimate as (37) for the CGO solutions to Maxwell's equations. Alternatively, in [24], it is shown that the linearization of the propagation equation and the internal measurements {H j } as an operator of the electric fields and the parameters, is an elliptic matrix-valued differential operator. The ellipticity is shown by plugging in proper CGO solutions to Maxwell's equations. Therefore, with sufficiently many measurements, (n(x), σ(x)) can be uniquely and stably reconstructed with no loss of derivatives.
It is not hard to derive that the linearization of (50) is Here δ q = k 2 δ n + ikδ σ denotes the perturbation of q 0 = k 2 n 0 + ikσ 0 , that is, q = q 0 + δ q for > 0 small, andẼ j = E j + δ Ej + o( ) where E j is the solution corresponding to q 0 . Also, taking the Fréchet derivative of ∆H j (equivalent to taking ∆ of the Fréchet derivative dH j ) yeilds where l.o.t. represents the lower order terms. Therefore, we obtain a system of linear equations: (51), the conjugate of (51) and (52) and A(x, D) is a second order 7J ×(6J + 2) matrix differential operator with the principal part in the Douglis-Nirenberg sense given by (54) A where I k is the k × k identity matrix. Here a j (x, D) and b j (x, D) are second order operators whose symbols are a j (x, ξ) = −|E j | 2 |ξ| 2 + 2k 2 σ 2 0 |q 0 | 2 |E j · ξ| 2 , b j (x, ξ) = 2k 4 σ 0 n 0 |q 0 | 2 |E j · ξ| 2 .
We remark that at this point the reason becomes self-explained for taking the laplacian of the internal function H j in our nonlinear system. Then for A 0 (x, ξ) to have full rank 6J +2 (Here J ≥ 2 so the system is not underdetermined) for every x ∈ Ω and ξ ∈ R 3 \{0}, one has to show that the rank of is 2. This is equivalent to show the following relation for every x ∈ Ω: To this end, we implement CGO solutions to the background Maxwell's equations (49) with q replaced by the background q 0 . These solutions were also used in ESC as mentioned in Proposition 3.1, originally constructed in [37]. Basically, Faddeev kernel allows the construction in higher order Sobolev spaces, which in turn provides an L ∞ bound at our disposal.
It is shown in [7] that the above linear problem is injective (i) when the coefficients v = ({E j , E * j } J j=1 , σ, n) are in a sufficiently small vicinity of an analytic coefficient (with the vicinity depending on that analytic coefficients); and (ii) when the domain Ω is sufficiently small.
The stability estimate presented in the above theorem then extends to the following nonlinear inverse problem (60) F (v) := A tF (v) = A t H in Ω, v| ∂Ω = v δ and ∂ ν v| ∂Ω = j δ .
Then defining v = v 0 + w and linearizing the above inverse problem about v 0 , we observe that the linear equation for w is precisely of the form (58), (59). This and the theory presented in [7] allow us to obtain the following result.
Theorem 6.6. Let us assume that the linear problem defined in (58) and (59) is injective.
Let v andṽ be solutions of (60) with respective source terms H andH and respective boundary conditions v δ andṽ δ as well as j δ andj δ . Then (H , v δ , j δ ) = (H ,ṽ δ ,j δ ) implies that v =ṽ, in other words the nonlinear hybrid inverse problem is injective. Moreover, we have the stability estimate . This estimate holds for C = C s when s > 7 2 . We refer the reader to [7] for similar type of analysis applied to other coupled-physics imaging inverse problems.