Lipschitz stability at the boundary for time-harmonic diffuse optical tomography

ABSTRACT We study the inverse problem in Optical Tomography of determining the optical properties of a medium , with , under the so-called diffusion approximation. We consider the time-harmonic case where Ω is probed with an input field that is modulated with a fixed harmonic frequency , where c is the speed of light and k is the wave number. We prove a result of Lipschitz stability of the absorption coefficient at the boundary in terms of the measurements in the case when the scattering coefficient is assumed to be known and k belongs to certain intervals depending on some a-priori bounds on , .


Introduction
Although Maxwell's equations provide a complete model for the light propagation in a scattering medium on a micro scale, on the scale suitable for medical diffuse Optical Tomography (OT) an appropriate model is given by the radiative transfer equation (or Boltzmann equation) [1]. If is a domain in R n , with n ≥ 2 with smooth boundary ∂ and radiation is considered in the body , then it is well known that if the input field is modulated with a fixed harmonic frequency ω, the so-called diffusion approximation leads to the complex partial differential equation (see [2]) for the energy current density u − div (K∇u) + (μ a − ik)u = 0, in .
( 1 ) Here k = ω/c is the wave number, c is the speed of light and, in the anisotropic case, the so-called diffusion tensor K, is the complex matrix-valued function where B ij = B ji is a real matrix-valued function, I is the n × n identity matrix and I−B is positive definite [2][3][4] on . The spacially dependent real-valued coefficients μ a and μ s are called the absorption and the scattering coefficients of the medium respectively and represent the optical properties of . It is worth noticing that many tissues including parts of the brain, muscle and breast tissue have fibrous structure on a microscopic scale which results in anisotropic physical properties on a larger scale. Therefore the model considered in this manuscript seems appropriate for the case of medical applications of OT (see [33]). Although it is common practise in OT to use the Robin-to-Robin map to describe the boundary measurements (see [2]), the Dirichlet-to-Neumann (D-N) map will be employed here instead. This is justified by the fact that in OT, prescribing its inverse, the Neumnannto-Dirichlet (N-D) map (on the appropriate spaces), is equivalent to prescribing the Robin-to-Robin boundary map. A rigorous definition of the D-N map for Equation (1) will be given in Section 2.
It is also well known that prescribing the N-D map is insufficient to recover both coefficients μ a and μ s uniquely [5] unless a-priori smoothness assumptions are imposed [6]. In this paper we consider the problem of determining the absorption coefficient μ a in a medium ⊂ R n , n ≥ 3, that is probed with an input field which is modulated with a fixed harmonic frequency ω = k/c, with k = 0 (time-harmonic case) and whose scattering coefficient μ s is assumed to be known. More precisely, we show that μ a , restricted to the boundary ∂ , depends upon the D-N map of (1), K,μ a , in a Lipschitz way when k is chosen in certain intervals that depend on a-priori bounds on μ a , μ s and on the ellipticity constant of I−B (Theorem 2.4). The static case (k = 0), for which (1) is a single real elliptic equation, was studied in [7], where the author proved Lipschitz stability of μ a and Hölder stability of the derivatives of μ a at the boundary in terms of K,μ a . In the present paper we show that in the time-harmonic case, for which (1) is a complex elliptic equation, a Lipschitz stability estimate of μ a at the boundary ∂ in terms of K,μ a still holds true if k is chosen within certain ranges. The case where μ a is assumed to be known and the scattering coefficient μ s is to be determined, can be treated in a similar manner. The choice in this paper of focusing on the determination of μ a rather than the one of μ s is driven by the medical application of OT we have in mind. While μ s varies from tissue to tissue, it is the absorption coefficient μ a that carries the more interesting physiological information as it is related to the global concentrations of certain metabolites in their oxygenated and deoxygenated states.
Our main result (Theorem 2.4) is based on the construction of singular solutions to the complex elliptic equation (1), having an isolated singularity outside . Such solutions were first constructed in [8] for equations of type div(K∇u) = 0, in , when K is a real matrix-valued function belonging to W 1,p ( ), with p > n and they were employed to prove stability results at the boundary in [8], [9], [10] and [11] in the case of Calderón's problem (see [12]) with global, local data and on manifolds. The singular solutions introduced in [8] were extended in [13] to equations of type with real coefficients, where K is merely Hölder continuous. Singular solutions were also studied in [14]. In this paper we extend the singular solutions introduced in [8] to the case of elliptic equations of type (1) with complex coefficients. Such a construction is done by treating (1) as a strongly elliptic system with real coefficients, since K ≥λ −1 I > 0, whereλ is a positive constant depending on the a-priori information on μ s , B and μ a . We wish to stress out, however, that in [8] the author constructed singular solutions to (3) which have an isolated singularity of arbitrary high order, where the current paper extends such construction to singular solutions to the complex Equation (1) having an isolated singularity of Green's type only. This is sufficient to prove the Lipschitz continuity of the boundary values of μ a in terms of the D-N map. The more general construction of the singular solutions with an isolated singularity of arbitrary high order for elliptic complex partial differential equations will be material of future work. This paper is stimulated by the work of Alessandrini and Vessella [15], where the authors proved global Lipschitz stability of the conductivity in a medium in terms of the D-N map for Calderón's problem, in the case when the conductivity is real, isotropic and piecewise constant on a given partition of . This fundamental result was extended to the complex case in [16] and in the context of various inverse problems for example in [17], [18], [19] and [20], [21], [22] in the isotropic and anisotropic settings, respectively. The machinery of the proof introduced in [15] is based on an induction argument that combines quantitative estimates of unique continuation together with a careful asymptotic analysis of Green's functions. The initial step of their induction argument relies on Lipschitz (or Hölder) stability estimates at the boundary of the physical parameter that one wants to estimate in terms of the boundary measurements, which is the subject of the current manuscript. Our paper also provides a first step towards a reconstruction procedure of μ a by boundary measurements based on a Landweber iterative method for nonlinear problems studied in [23], where the authors provided an analysis of the convergence of such algorithm in terms of either a Hölder or Lipschitz global stability estimates (see also [24]). We also refer to [25] and [32] for further reconstruction techniques of the optical properties of a medium.
The paper is organised as follows. Section 2 contains the formulation of the problem (Subsections 2.1 and 2.2) and our main result (Subsection 2.3, Theorem 2.4). Section 3 is devoted to the construction of singular solutions of Equation (1) having a Green's type isolated singularity outside . The proof of our main result (Theorem 2.4) is given in Section 4.

Main assumptions
We rigorously formulate the problem by introducing the following notation, definitions and assumptions. For n ≥ 3, a point x ∈ R n will be denoted by x = (x , x n ), where x ∈ R n−1 and x n ∈ R. Moreover, given a point x ∈ R n , we will denote with B r (x), B r (x ) the open balls in R n , R n−1 , centred at x and x respectively with radius r and by Q r (x) the cylinder We will also denote B r = B r (0), B r = B r (0) and Q r = Q r (0).

Definition 2.1:
Let be a bounded domain in R n , with n ≥ 3. We shall say that the boundary of , ∂ , is of Lipschitz class with constants r 0 , L > 0, if for any P ∈ ∂ there exists a rigid transformation of coordinates under which we have P = 0 and where ϕ is a Lipschitz function on B r 0 satisfying We consider, for a fixed k > 0, where K is the complex matrix-valued function and q is the complex-valued function We recall that I denotes the n × n identity matrix, where the matrix B is given by the OT physical experiment and it is such that B ∈ L ∞ ( , Sym n ), where Sym n denotes the class of n × n real-valued symmetric matrices and such that I−B is a positive definite matrix [2][3][4]. In this paper, we assume that the scattering coefficient μ s is also known in and it is the absorption coefficient μ a that we seek to estimate from boundary measurements. We assume that there are positive constants λ, E and E and p > n such that the known quantities B, μ s and the unknown quantity μ a satisfy the two assumptions below respectively.

Assumption 2.1 (Assumption on µ s and B):
and

Assumption 2.2 (Assumption on µ a ):
We state below some facts needed in the sequel of the paper. Most of them are straightforward consequences of our assumptions.
The inverse of K has real and imaginary parts given by the symmetric, real matrix valued-functions on respectively. As an immediate consequence of Assumptions 2.1 and 2.2 we have for a.e. x ∈ and any ξ ∈ R n . Moreover K −1 R and K −1 I commute, therefore the real and imaginary parts of K are the symmetric, real matrix valued-functions on respectively. Assumptions 2.1 and 2.2 also imply that for a.e. x ∈ , for every ξ ∈ R n and the boundness condition for a.e. x ∈ . Moreover K = {K hk } h,k=1,...,n and q satisfy and respectively, where C 1 is a positive constant depending on λ, E, E, k and n. By denoting q = q R + iq I , the complex equation is equivalent to the system for the vector field u = (u 1 , which can be written in a more compact form as − div(C∇u) + qu = 0, in (27) or, in components, as where {C hk lj } h,k=1,...,n is defined by and {q lj } l,j=1,2 is a 2 × 2 real matrix valued function on defined by (20), together with (22) imply that system (26) is uniformly elliptic and bounded, therefore it satisfies the strong ellipticity condition where C 2 > 0 is a constant depending on λ, E, k and n.
is uniformly positive definite on and it satisfies

Definition 2.2:
We will refer in the sequel to the set of positive numbers r 0 , L, λ, E, E introduced above, along with the space dimension n, p > n, the wave number k and the diameter of , diam( ), as to the a-priori data.

The Dirichlet-to-Neumann map
Let K be the complex matrix valued-function on introduced in (6) and q = μ a − ik, satisfying Assumptions 2.1 and 2.2. B and μ s are assumed to be known in and satisfying Assumption 2.1, so that K is completely determined by μ a , satisfying Assumption 2.2, on . Denoting by ·, · the L 2 (∂ )-pairing between H 1/2 (∂ ) and its dual H −1/2 (∂ ), we will emphasise such dependence of K on μ a by denoting K by

Definition 2.3:
The Dirichlet-to-Neumann (D-N) map corresponding to μ a is the operator defined by and ϕ ∈ H 1 ( ) is any function such that ϕ| ∂ = g in the trace sense.

The main result
or k ≥k 0 := where, λ and E are the positive numbers introduced in Assumptions 2.1 and 2.2, then where C > 0 is a constant depending on n, p, L, r 0 , diam( ), λ, E, E and k.

Singular solutions
We consider where K = {K hk } h,k=1,...,n and q are the complex matrix valued-function and the complex function respectively introduced in Section 1 and satisfying Assumptions 2.1 and 2.2 on B R . and furthermore where w satisfies r<|x|<2r Here α is such that 0 < α < 1 − n/p, and C is a positive constant depending only on α, n, p, R, λ, E, E and k.

Remark 3.2:
Since K −1 (0) is a complex matrix, the expression appearing in the leading term in (42) is defined as the principal branch of (45), where a branch cut along the negative real axis of the complex plane has been defined for z 1/2 , z ∈ C. Expressions like (45) will appear in the sequel of the paper and they will be understood in the same way.
Next we consider two technical lemmas that are needed for the proof of Theorem 3.1. The proofs of these results for the case where L = −div(K∇·), with K a real matrix valued-function, are treated in detail in [8] and their extension to the more general case L = −div(K∇·) + q, with K, q a real matrix valued-function and a real function respectively, was extended in [7], therefore only the key points of their proof will be highlighted in the complex case below.   LG(·, y) = δ(· − y)I, forall y ∈ B R (57) in the sense that for every φ = (φ 1 , Moreover where C is a positive constant depending on n, λ, E, E and k and the vector valued-function u = (u 1 , u 2 ) defined by satisfies where f = (f 1 , f 2 ) and For the existence, uniqueness and asymptotic behaviour of the Green's matrix G on B R as in (57)-(59) we refer to [28]. We also refer to [29], [30] and the more recent result [31] for further reading on the issue of the Green's matrix for elliptic systems of the second order. By an argument based on the monotone convergence theorem, one can show that I 1 and I 2 are both bounded from above by C|x| 2−s , where C is a positive constant depending on A, s, n, p, R, λ, E, E and k.
such that |f N | ≤ |f | on B R , therefore ||f N || L p (˜ ) ≤ ||f || L p (˜ ) , for any˜ ,˜ ⊂⊂ B R \ {0}, for any N ≥ 1. By applying interior L p -Schauder estimates to u N and using the fact where C is a positive constant that depends on˜ . By applying a diagonal process we can find a sub- . This limit satisfies both (55) and (56).
We proceed next with the proof of Theorem 3.1.

Proof of Theorem 3.1.: We start by considering
where L 0 := −div(K(0)∇·) on B R . We want to find w such that satisfying (43), (44), where L is defined by (5). We have Therefore for any r, 0 < r < R/2 we have where β = 1 − n/p and C is a positive constant depending on λ, E, E, R and k only. If we take w ∈ W

Proof of the main result
Since the boundary ∂ is Lipschitz, the normal unit vector field might not be defined on ∂ . We shall therefore introduce a unitary vector field ν locally defined near ∂ such that: (i) ν is C ∞ smooth, (ii) ν is non-tangential to ∂ and it points to the exterior of (see [9, Lemmas 3.1-3.3] for a precised construction of ν). Here we simply recall that any point where τ 0 and C depend on L, r 0 only.

Remark 4.1:
Several constants depending on the a-priori data introduced in Definition 2.2 will appear in the proof of the main result below. In order to simplify our notation, we shall denote by C any of these constants, avoiding in most cases to point out their specific dependence on the a-priori data which may vary from case to case.