Instability of the linearized problem in multiwave tomography of recovery both the source and the speed

In this paper we consider the linearized problem of recovering both the sound speed and the thermal absorption arising in thermoacoustic and photoacoustic tomography. We show that the problem is unstable in any scale of Sobolev spaces.


Introduction
In multiwave tomography, one sends one type of wave to the body of a patient, most often electromagnetic (thermoacoustic tomography) or optical radiation (photoacoustic tomography) which interacts with the tissue, and measures the acoustic signal on the boundary generated by this interaction. This combines the high contrast of the incoming waves with the high resolution of the measured ultrasound ones. The mathematical model of the emitted ultrasound wave is the following. Let u solve the problem where the sound speed c = c(x) > 0 and T > 0 are fixed. Assume that f and c − 1 are supported inΩ, where Ω ⊂ R n is some bounded domain with a smooth convex boundary. The measurements are modeled by the operator (2) Λ 1 f := u| [0,T ]×∂Ω .
The first step in multiwave imaging is to recover f given Λ 1 f . The speed c is usually assumed to be known but in practice, it is not. Then the natural question is whether we can recover both c and f . The answer is still unknown. A discussion of this problem, together with a partial local result stating that c can be determined up to a constant scaling can be found in [10]. David Finch observed a link between this problem and the transmission eigenvalues. The problem of recovery of f , given c has received a lot of attention in the past years. We refer to [1,2,3,4,5,6,7,8,9,10,13,14,15,17,18,23,24,28,29] for some works in this direction. If c is constant, then the inversion of Λ 1 is actually an integral geometry problem as well. For variable c, there is always uniqueness if T ≫ 1, and there is stability if c is non-trapping and T ≫ 1. We refer to [23,20] for details.
When an inaccurate speed is used for reconstruction by time reversal, the resulting images are distorted, and full of "artifacts", see for example the images in [12] or [10]. The mathematical structure of the artifacts is easy to understand: the forward map Λ 1 is a Fourier Integral Operator (FIO) with a canonical relation given by the graph of the map determined by the geodesics rays from the interior to the boundary, see [23] for more details. The time reversal is another FIO with a Date: May 5, 2014. First author partly supported by a NSF Grant DMS-0800428. Second author partly supported by NSF and a Walker Family Endowed Professorship. canonical relation given by the graph of the inverse of the former one. When the speed is different, we get that time reversal with a wrong speed is an FIO with a canonical relation not the diagonal (the graph of the identity) but the graph of the composition of the forward one and a backward one with a different speed. In particular, since the canonical relation of Λ 1 and its reversal have two disconnected components, one can see double images of the same singularity, well visible in the examples in [12], for instance. An algorithm to tune in the speed by maximizing the sharpness of the reconstructed f is proposed in [26]. This is related to the FIO description of the artifacts but good mathematical understanding of this algorithm is lacking.
One of the ways to recover the speed is to take additional measurements and to recover c from travel times. This is the method proposed in [27] in thermoacoustic tomography. The travel time problem is stable under some geometric assumptions on c (c −2 dx 2 being a simple metric in the domain is enough) which are satisfied when c is close enough to a constant, in particular, see, e.g., [19] and the references there. Then f can be recovered stably as well. Additional data can be provided either by placing ultrasound sources around the body, or by placing passive absorbing (tissue imitating) objects around the body which become ultrasound sources by the thermo-or the photo-acoustic effect. We refer to the recent paper [12] for references and numerical and experimental implementations of that method. As can be expected, the results are very good.
Simultaneous reconstruction of f and c aside from the above mentioned work [26], have been tried with various success in [31,30,32], for example.
In this paper, we study the linearization δΛ 1 and we show that the latter is unstable. In particular, we prove the following. We also show that a conditional type of stability estimate cannot hold either, see Remark 1. This suggests instability of the non-linear problem as well but does not imply it directly. Stability of the linearization in some Sobolev norms, even if not in the sharp ones, does imply (conditional Hölder) stability of the non-linear problem [22]. The converse however is a much more delicate question. To prove the instability of the linearization, we show that the latter is a smoothing operator for (δf, δc 2 ) belonging to an explicitly defined infinite dimensional linear space, see (39). We refer to section 4 for more details.

Preliminaries
Notice first that c 2 ∆ is formally self-adjoint w.r.t. the measure c −2 dx. Given a domain U , and a function u(t, x), define the energy In particular, we define the space H D (U ) to be the completion of C ∞ 0 (U ) under the Dirichlet norm It is easy to see that H D (U ) ⊂ H 1 (U ), if U is bounded with smooth boundary, therefore, H D (U ) is topologically equivalent to H 1 0 (U ). If U = R n , this is true for n ≥ 3 only. By the finite speed of propagation, the solution with compactly supported Cauchy data always stays in H 1 even when n = 2. The energy norm for the Cauchy data [f 1 , f 2 ], that we denote by · H is then defined by This defines the energy space Here and below, L 2 (U ) = L 2 (U ; c −2 dx). Note also that The wave equation then can be written down as the system where u = [u, u t ] belongs to the energy space H. The operator P then extends naturally to a skewselfadjoint operator on H if c ∈ L ∞ , and c −1 ∈ L ∞ . In this paper, we will deal with either U = R n or U = Ω. In the latter case, the definition of H D (U ) reflects Dirichlet boundary conditions. We generalize next the results in [23] to the inverse problem with general Cauchy data (f 1 , f 2 ) in (1) with g not necessarily zero. What we really need later is Proposition 2 only. Let u solve the problem Assume that f is supported inΩ, where Ω ⊂ R n is some smooth bounded domain. Set The trace Λf is well defined in C (0) [0, T ]; H 1/2 (∂Ω) , where the subscript (0) indicates that we take the subspace of functions h so that h = 0 for t = 0. For a discussion of other mapping properties, we refer to [11]. When f is restricted to functions in H, supported in a fixed compact K ⊂ Ω, then Λ is an FIO with a canonical relation of graph type, and maps f continuously into [23].
Then we define the following pseudo-inverse By [16], is a continuous map. Note that the mapping properties above allow us to apply A to Λf only when f is compactly supported in Ω but the theorem above shows that AΛ extends continuously to the whole H. Let T (Ω) be the length of the longest geodesic inΩ, when (Ω, c −2 dx 2 ) is non-trapping.
, and Λ has an explicit left inverse of the form where u solves (6) with a given f ∈ H. Let v be the solution of (8) with h = Λf . Then v + w solves the same initial boundary value problem in [0, T ] × Ω that u does (with initial conditions at t = T ), therefore u = v + w. Restrict this to t = 0 to get . We will show now that K extends to a compact operator. Since T > T (Ω), all singularities starting fromΩ leaveΩ at t = T . Therefore, u(T, ·) and u t (T, ·), restricted toΩ, are C ∞ . Moreover, considered as linear operators of f , they are operators with smooth Schwartz kernels. Then so is φ, see (9), by elliptic regularity. Therefore, the map defined a priori on C ∞ 0 (Ω) × C ∞ 0 (Ω) functions, extends to a compact one because it is an operator with smooth kernel onΩ. Since the solution operator of (12) from t = T to t = 0 is unitary in H(Ω), we get that the map H(Ω) ∋ f → [w(0, ·), w t (0, ·)] ∈ H D (Ω)) is compact, too, as a composition of a compact and a bounded one.
We know remove the smoothness restriction on f , and let it be any element in H. In what follows, (·, ·) H D (Ω) is the inner product in H D (Ω), see (3), applied to functions that belong to H 1 (Ω) but maybe not to H D (Ω) (because they may not vanish on ∂Ω). Set u T := u(T, ·). By (4) and the fact that . Therefore, the energy of the initial conditions in (12) satisfies the inequality (13) E . Since the Dirichlet boundary condition is energy preserving, we get Therefore, . We show next that actually the inequality above is strict, i.e., (15) Kf H(Ω) < f H(Ω) , f = 0.
Assume the opposite. Then for some f = 0, all inequalities leading to (14) are equalities. In particular, E Ω (w, T ) = E R n (u, T ). Then By the finite domain of dependence then One the other hand, we also have Therefore, In [23], we used the fact that u extends to an even function of t that is still a solution of the wave equation because f 2 = 0 there. Then one gets that (18) actually holds for |t| < 3T /2. Then one concludes by Tataru's unique continuation theorem [25] that u = 0 on [0, T ] × Ω, therefore, f = 0.
We also noted there that the time interval [0, T ] is actually larger (twice as large) than what we need for the Neumann series to converge, see [24,20], where T > T (Ω)/2 only. In the case under consideration, f 2 does not necessarily vanish. We modify the arguments as follows. From John's theorem (equivalent to Tataru  Proof. Note first, that A has the mapping properties above, by the results in [16]. Next, (Id−K) −1 is a bounded map in H by Theorem 2.
Let B 1,2 be the components of B (that sends scalar functions to vector functions), i.e., Bh = (B 1 h, B 2 h). Let Λ 1,2 be the components of Λ (that sends vector functions to scalar functions, i.e., Λf = Λ 1 f 1 + Λ 2 f 2 . We can think of B as a 2x1 matrix, and of Λ as an 1x2 matrix. Then Set ∂ −1 t g = t 0 g(s) ds and let ∆ D be the Dirichlet realization of ∆ in Ω. We have the following.
Relations (21) can be used to show that Λ 2 has a left inverse based on the result in [23] only.
Analyzing the mapping properties of B ′ 2 and Λ 2 as in [23], one gets that the composition B 1 Λ 1 is well defined and is a bounded operator on L 2 (Ω); that therefore has to be identity.
3. Recovery of the speed c when f is known; The linearization δΛ 1 /δc 2 w.r.t. the speed The non-linear problem of recovery of c when f is known, and its linearization were studied in detail by the authors in [21]. We showed there that if δc 2 is a priori supported in a compact subset K, then the linearization δΛ 1 /δc 2 is Fredholm. We also gave conditions for uniqueness for both the linear and the non-linear problem. The geometric requirement is that there exists a foliation ofΩ by strictly convex surfaces. If n = 2, non-trapping implies that condition. In all dimensions, if c is close enough to a constant, that condition holds. In particular, if (22) x · ∂c < c inΩ, then the parts of the spheres |x| = R, R > 0, intersectingΩ form a foliation of surfaces convex w.r.t. the metric c −2 dx 2 , see [20,21]. There is an another condition as well: we require that ∆f (x) = 0, ∀x ∈ K. This condition fits well into our analysis; it is an if and only if condition for the operator Q N below to be elliptic.

Analysis of the linearized operator δΛ 1
We derive a representation of the differenceΛ 1f −Λ 1 f first. If one is interested in the linearization only, the computations are not really simpler -one just has to drop the tilde in most of the formulas.
Let (c, f ), (c,f ) be two pairs, and let u,ũ be the corresponding solutions of (1). Then We have and let P D be the Dirichlet realization of P in Ω. Recall the notation ∂ −1 t g = t 0 g(s) ds. Denote by Ru with u = [u 1 , u 2 ] the restriction of u 1 to [0, T ] × ∂Ω. We have Take the t derivative to get (26) w Differentiate W to get where we used the fact that ∂ tũ = 0 for t = 0, therefore ∂ t ∆ũ = 0 for t = 0. Differentiate one more time to get Relations (26) and (27) show that W = O(t 2 ) at t = 0. Therefore, By Proposition 2, for the first term on the right we get The second term on the r.h.s. of (29) can be written in the form We claim that I 2 = I ′ 2 , where Indeed, I ′ 1 = 0 for t = 0, and direct differentiation shows that because RP −1 D = 0. Therefore, ∂ t I ′ 1 = 0 for t = 0 as well. Differentiate one more time to get I 2 = I ′ 2 . Therefore, Compare with (26) and repeat these arguments to get for any N = 1, 2, . . .
We want to emphasize that in all those formulas, h is considered as an operator of multiplication, i.e., Λ 2 P −k D hP k+1f = Λ 2 (P −k D (hP k+1f )), etc. Combine this with (24), (26), to get the following.
SinceΛ 1f − Λ 1 f = 0 for t = 0, integrating w.r.t. t, and applying Proposition 2, we get Finally, using the arguments above, we write the last term on the right as We therefore proved the following.
Corollary 1. For any N = 1, 2, . . . , the linearized operator δΛ 1 {δf, δc 2 } at (c, f ) has the form Proof of Theorem 1. Note that for (c, f ) ∈ C ∞ , the operator R N is smoothing 2N + 2 degrees, i.e., , and Q N is a ΨDO (in the interior of Ω) with principal symbol ∆f (x)|ξ| −2 . So for δc 2 supported in a compact set in Ω, we can write where R ∞ is smoothing, and l.o.t. stands for "lower order terms", i.e., for a pseudo-differential operator of order −3, which we can always assume to have a proper support. This result is not surprising. All singularities of the kernel of δΛ 1 are of conormal type, at dist(x, y) = t, where dist is the distance in the metric c −2 dx 2 . This can be seen from the representation The operator Q N can be explained by those singularities. Next, R N depends on the smooth part of the kernel, and in particular, the behavior of the kernel inside that geodesic ball.
An important observation is that if the expression in the parentheses in (33) vanishes (and this is an explicit condition on δf and δc 2 ), then the linearization is very smooth. If for a moment we ignore R ∞ , we get a kernel of infinite dimension. Indeed, given δc 2 , compactly supported in Ω, we can always find δf so that δf + c −2 ∆f |D| −2 + l.o.t. δc 2 = 0. If we are given a compactly supported (in Ω) function δf , and if ∆f = 0 on supp δf , we can find δc 2 at least away from a finitely dimensional space, so that the expression above vanishes as well. The effect of R ∞ will be the following. Even though we still do not know whether we can determine both δf and δc 2 , we can claim that even of we could, the problem would be unstable.
More precisely, given N > 1, let V N be the linear space This can easily be generalized to (δh, f ) belonging to negative Sobolev spaces. We then complete the proof of Theorem 1 by a well known argument, see, e.g., [22] or the remark below.