Multi-wave imaging in attenuating media

We consider a mathematical model of thermoacoustic tomography and other multi-wave imaging techniques with variable sound speed and attenuation. We find that a Neumann series reconstruction algorithm, previously studied under the assumption of zero attenuation, still converges if attenuation is sufficiently small. With complete boundary data, we show the inverse problem has a unique solution, and modified time reversal provides a stable reconstruction. We also consider partial boundary data, and in this case study those singularities that can be stably recovered.


Introduction
A multi-wave imaging technique couples a high-resolution imaging modality with a high-contrast imaging modality through a physical mechanism. In medical imaging, ultrasound waves are often used, as they can provide sub-millimeter resolution in tissue. The typical examples of such techniques are thermoacoustic tomography (TAT) and photoacoustic tomography (PAT). The former uses low-frequency microwaves as the contrast modality; the latter uses near-infrared light. Both are coupled to ultrasound through the photoacoustic effect. Acoustic pressure is then recorded using transducers placed along the boundary of the region of interest.
To reconstruct an image from the acoustic data, one solves two coupled inverse problems. First, one reconstructs the ultrasound source inside the region using the acoustic data on the boundary. Then, the ultrasound source is used to reconstruct the absorption of the contrast modality within the tissue. We do not consider the second problem here. In the case of either TAT or PAT, it is called either quantitative thermoacoustic tomography (QTAT) or quantitative photoacoustic tomography (QPAT).
The problem of determining the source that induces a measured acoustic pressure at the boundary of a region has been well-studied under the assumptions of constant speed and zero attenuation. In this case, the problem is essentially the inversion of a spherical mean operator [1,7]. However, in tissue these assumptions do not hold. Variations in either acoustic speed [10] or attenuation [4][5][6]17] can cause artifacts in reconstruction. This has spurred the development of reconstruction methods using models of ultrasound propagation that include variable acoustic speed and attenuation.
The method of time reversal provides such a reconstruction [8,9]. In this method, one fixes a time T (possibly infinite) and solves a mixed boundary problem for the acoustic wave equation with zero Cauchy data at time T and the measured data along the boundary for t ∈ (0, T ). This exploits the fact that the acoustic wave equation is invariant under time reversal, i.e., t → −t. To ensure the measured data is compatible with the Cauchy data, a cut-off function is often used. This technique has also been extended to the thermoelastic equation [2,11].
It is also possible to choose Cauchy data that is compatible with the measured data, thereby avoiding the use of a cut-off function. In [22], they demonstrate that a modified time reversal of this type provides a stable reconstruction in the acoustic wave model. It has also been extended to the elastic wave equation [25]. One goal of this work is to extend modified time reversal to a model including variable attuenation. An advantage of modified time reversal is that under some circumstances it yields an explicit Neumann series for the ultrasound source. This Neumann series was studied numerically in [18] and was shown to be an improvement over time reversal.
There are several commonly-used models of attenuation in biological tissue [12,13,20]. Generally speaking, the attenuating effect of a medium depends upon the frequency of the wave propagating through it, with higher frequency waves experiencing more attenuation. A common approximation posits a power law between frequency and attenuation [16,21,26]. Here we take the damped wave equation as a simple model of ultrasound attenuation. We follow the undamped model of [3] in assuming that the ultrasound source is pulse of the form f (x)δ ′ (t). This approximation can be made in the case when the contrast media is an electromagnetic wave, as it is in both TAT and PAT. Writing the damped wave operator as a = ∂ 2 t + a∂ t − c 2 ∆, we then consider solutions of in the sense of distributions. Here c is a smooth, positive function representing the sound speed. If Ω is our region of interest, we assume c is equal to one on R n \ Ω. The attenuation coefficient a is assumed to be a smooth, nonnegative function. If we take f to be in H 1 0 (Ω), then solutions of the Cauchy problem are also solutions of (1) when extended by zero to (−∞, 0) × R n . To show this, take u to be a smooth solution of (2) on R n+1 . Then we may consider the distribution H(t)u(t, x), where H(t) is the Heaviside function. We have that a (Hu) = uδ ′ + 2(∂ t u)δ + auδ + ( a u)H.
The last term is zero, as a u = 0. Let φ ∈ C ∞ c (R n+1 ) be a test function. Then which implies that u satisfies (1). This justifies the choice of boundary conditions in (2). Let Γ = ∂Ω be the boundary of region of interest. If we have access to pressure data on the entire boundary, we can model the measurements with the operator Λ a f = u| (0,T )×Γ . For the partial boundary case, we will take a cut-off function χ supported on (0, T )×Γ ′ , with Γ ′ ⊂ Γ, and model the partial boundary data with χΛ a . As (f, −af ) is in the energy space H = H 1 0 (Ω)×L 2 (Ω), it follows from [15] that Λ a f is a bounded operator from H D (Ω) to H 1 ((0, T ) × Γ). The inverse problem is to reconstruct f from the measured data h = Λ a f . In the damped case, the modified time reversal of [22] suggests that we construct a potential pseudoinverse A a h from the backward problem with mixed boundary data, Here φ satisfies ∆φ = 0 on Ω and φ| Γ = h| t=T . By the trace theorem, h| t=T ∈ H 1/2 (Γ), so the problem is well-posed. We then define A a h = v| t=0 .
Our first goal is to determine whether modified time reversal is an effective reconstruction method for the damped wave model. In particular, we are interested in the error operator K a = A a Λ a − I. In the undamped case, under some geometric assumptions, K 0 is a strict contraction [22]. This implies that the Neumann series converges. This series can be used to numerically reconstruct f , and in [18] it was shown to be an improvement over the usual time reversal algorithm. In the second section, we show that for small attenuations the above Neumann series is still convergent. In the third section, we consider the uniqueness and stability of modified time reversal in the presence of arbitrary attenuation, given complete data. In the fourth section, we show that in some cases partial data is sufficient to recover the singularities of f . In the last section we prove an a priori estimate for the damped wave equation with mixed boundary data, following an earlier estimate [15].

Modified time reversal
We begin with some preliminary remarks on the damped wave equation. We will assume that the initial data (f, −af ) lies in the energy space By Poincaré's inequality, the norm of H D (Ω) is equivalent to that of H 1 0 (Ω). Equip H with the inner product which makes H into a Hilbert space.
Let U ⊂ R n . For an arbitrary function w that is in both C 1 (0, T ; L 2 (U )) and C(0, T ; H 1 (U )), define the local energy functional Both solutions of the wave equation and the damped wave equation are such functions, and it is well-known that energy is conserved in the first case and nonincreasing in the second case. For use later on, we state an a priori estimate for solutions of the nonhomogeneous damped wave equation with mixed boundary data. We defer the proof to the last section.

Proposition 1.
Let Ω be a smooth domain, and u be a solution of subject to the compatibility condition h| t=0 = f | Γ . Then there exists a constant C > 0 such that provided the boundary data are sufficiently regular. Now, take u to be the solution of (2) and take v to be the corresponding solution of (3) with h = Λ a f . To study the error operator K a = A a Λ a − I, we study the function w = u − v, which solves the IBVP, Recall φ is the harmonic extension of u| {T }×Γ to Ω, and so the boundary data is compatible. Let W = (w, ∂ t w). Then (4) reduces to the first-order system, where Q a is the operator The following energy estimate for solutions of (4) is used often in what follows.
Proof. Define the operator By [19,Theorem X.48], if Re(Hf, f ) ≥ 0 for all f ∈ D(H), then e −tH is a C 0semigroup of contractions for t ≥ 0. This is true of Q a + ||a|| ∞ as an operator on H. Therefore ||e −tQa || ≤ e t||a||∞ for t ≥ 0. One can also verify directly that ω(t) = e −(T −t)Qa (ω 1 , ω 2 ) is the solution of the problem in the lemma.
Hence for 0 ≤ t ≤ T , yielding the desired energy estimate.
This is a quantitative estimate of the nonincrease of energy in the damped wave equation, in a time-reversed sense. When solved backward, the damped wave equation gains energy at most exponentially. One consequence of this estimate is a bound on ||K a f || HD .
Proof. Again let W = (w, ∂ t w). By Lemma 1, Using the fact that φ is harmonic, we see that Hence By itself, this estimate cannot show K a is a strict contraction on H D (Ω), this being sufficient for the convergence of the Neumann series.

Convergence of the Neumann series
Let T 0 be the supremum of the length of geodesics in (Ω, c −2 dx 2 ). If the metric is trapping, this is infinite. In [22], they show that if T 0 < T < ∞, then ||K 0 || < 1. In that case, the Neumann series converges. This series can be used as an effective numerical reconstruction algorithm [18]. We extend this convergence result by continuity to the damped wave model, given sufficiently small attenuation.
Theorem 1. If T 0 < T < ∞, then there exists a 0 > 0 such that for all a such that ||a|| ∞ < a 0 , the Neumann series converges.
To show this, we use the continuity of the map a → K a .
To prove this, we return to system of coupled PDEs defining K 0 and K a , in order to compare their difference. (4), v solves (3), and u solves (2). Defineū = u 0 − u andw = w 0 − w. Then these functions satisfy the system of PDE Our interest is in E Ω (w, 0), and to determine it we make use of Proposition 1. Note that (φ 0 − φ) is harmonic. It is also equal to zero on Γ, and therefore compatible with the Cauchy data. As before E Ω (ū, T ) is a bound on the energy of the Cauchy data. By Proposition 1, Concerning the first term, by Jensen's inequality, Then by Lemma 2, we have Hence ||a∂ t w|| L 1 (0,T ;L 2 (Ω)) ≤ a 0 (1 + Ca 2 0 ) 1/2 e T a0 ||f || HD To estimate the second term, we use Duhamel's principle to representū as the integralū By nonincrease of energy we have and as before E Ω (u, s) ≤ (1 + C||a|| 2 ∞ )e 2T ||a||∞ ||f || 2 HD . Hence by Jensen's inequality again, Combining the contributions of both terms and Theorem 1, we have which completes the estimate. Note that the constant C does not depend on ||a|| ∞ . Now we return to the proof of Theorem 1.
Proof. By the previous proposition, By assumption ||K 0 || < 1. As a 0 → 0, the second term tends to zero, as C does not depend on a 0 . Hence there exists a 0 sufficiently small such that ||K a || < 1 uniformly.
If ||K a || < 1, then the Neumann series of Theorem 1 also implies stability in the case of complete data.

Complete data
In this section, we use microlocal methods to prove that under some geometric assumptions and with complete boundary data, modified time reversal is stable. This follows directly from the previous section, provided the attenuation is sufficiently small. Here we show stability without this assumption. For convenience, we adopt the notation of [22] for unit speed geodesics. We identify T * Ω \ 0 with T Ω \ 0, via the Euclidean inner product.
In this section and the following one, we assume a priori that supp f is bounded away from Γ. This can always be accomplished in practice by slightly enlarging the region of interest, if necessary.

Uniqueness
First, we show uniqueness. Define d(x, Γ) to be the infimum of the lengths of curves (with respect to c −2 dx 2 ) starting at x and ending at Γ. In the undamped case, [22,Theorem 2] yields uniqueness under the assumption that Ω is strictly convex and that T > T 1 (Ω), where This result exploited the invariance of the wave equation under time reversal. To compensate for the lack of this symmetry in the damped wave model, we show uniqueness under the assumption T > 2T 1 (Ω).
Proof. Let u be the solution of the forward problem (2). We apply Tataru's theorem [23] (actually, the special case [22,Theorem 4]). Let U be a neighborhood of x 0 ∈ R n \ Ω with U ∩ Ω = ∅. By assumption, u = 0 on {0} × (R n \ Ω) and (0, T ) × Γ, and satisfies the damped wave equation on the exterior of (0, T ) × Ω. Therefore u = 0 on (0, T ) × U . By Tataru's theorem, u = 0 on the intersection of the backward light cone at (T, x 0 ) and the forward light cone at (0, x 0 ). Repeating this argument for all x 0 on the exterior of Ω sufficiently close to Γ, we find that u = 0 on

Stability
To show the reconstruction is stable even when attenuation is large, we turn to microlocal analysis. When a = 0, stability was shown in [22,Theorem 3] under the assumption that T was sufficiently large so that for every (x, ξ) ∈ S * Ω, either γ x,ξ or γ x,−ξ hit Γ before time T . Recall from the previous section that T 0 is the supremum of the lengths of geodesics in (Ω, c −2 dx 2 ). Here we present a stability estimate in the special case where c −2 dx 2 is non-trapping, i.e., T 0 < ∞. That this assumption is not necessary for stability follows from the result for partial data in the next section. However, the proof in this case is simpler. Notice that in the non-trapping case, if T > T 0 (Ω), then T > 2T 1 (Ω) as well, and so Theorem 2 implies that Λ a is injective.
In the following theorem, let χ ∈ C ∞ (R × Γ) be a fixed cutoff function equal to one on [0, T 0 ] × Γ and zero in a neighborhood of {T } × Γ. Theorem 3. Assume T 0 < T < ∞. Then the following are true: where R a is a smoothing operator.
Proof. First, we show A a Λ a is Fredholm. In the second section, we considered the error operator K a = f − A a Λ a f . It remains to show that K a is compact.
Let u be the solution of the forward problem (2). By propagation of singularities, WF(u) is the union of null bicharacteristics lying over WF(f ). Let h = Λ a f . As we assume T > T 0 (Ω), WF(u) is empty over a neighborhood of {T } × Ω. Therefore the operator f → (u| t=T − φ, ∂ t u| t=T ) maps H D (Ω) to C ∞ (Ω) × C ∞ (Ω), and hence is compact. Finally, Proposition 1 implies that the operator (u| t=T − φ, ∂ t u| t=T ) → w| t=0 is bounded from H → H D (Ω). Therefore K a is compact, and A a Λ a is Fredholm.
Let v be the solution of the corresponding backward problem (3) with h = χΛ a f , and w = u − v as before. Evidently v| t=T = 0. Therefore w solves (4) with smooth boundary data, compatible to arbitrary order. If we define R a f = w| t=0 , then R a maps H D (Ω) to C ∞ c (Ω), and A a χΛ a = I + R a . To show the first inequality, the previous equation shows ||f || HD ≤ ||A a χΛ a f || HD + ||R a f || L 2 .
By Lemma 2, A a is bounded from H 1 (Γ) → H D (Ω), and so we have To obtain the stability estimate in the theorem, we recall that since T > 2T 1 (Ω), Λ a is injective. So by [24, Proposition V.3.1], there is another constant C > 0 such that ||f || HD ≤ C||Λ a f || H 1 .
Therefore the reconstruction is stable. If in addition Ω is strictly convex (with respect to the Euclidean metric), the results of the next section show that singularities that propagate through the boundary may be stably recovered; the larger T is, the more singularities become visible.
Concerning the second part of the theorem, in principle one could choose a different modification of the Cauchy data in (3) that would be compatible with h to at least first order. In that case, K a would be smoothing of at least first order, and no cut-off function would be necessary.

Partial data
In applications it is often only possible to obtain data on a proper subset of the boundary. Let Γ ′ ⊂ Γ be relatively open. We model partial measurements by the operator χΛ a , where χ ∈ C ∞ (R × Γ) is a cut-off function positive on (0, T ) × Γ ′ and zero elsewhere. In this section, we only show that "visible" singularities (in the sense of [14,22]) may be stably recovered from χΛ a f , by proving an estimate of the form ||f || HD ≤ C(||Λ a f || H 1 + ||f || L 2 ).
If in addition, we know that Λ a is injective, we could prove stability as in the complete data case. In the undamped case, [22,Proposition 2] shows uniqueness provided that d(x, Γ ′ ) < T . However, their result relies on being able to extend solutions of the wave equation to (−T, 0) × Ω by time reversal. It is not clear how to modify their proof to account for attenuation. Therefore, we constrain ourselves to showing stable recovery of visible singularities.
Let K be an open set strictly contained in Ω. By analogy with the complete data case, let T 2 (K, Γ ′ ) be the least time necessary for all (x, ξ) ∈ S * K to have the property that at least one of γ x,ξ and γ x,−ξ reaches Γ ′ before time T 2 (K, Γ ′ ). In other words, we can say that K is visible from Γ ′ . We will assume that T > T 2 (K, Γ ′ ). Our goal in this section is to show the following: Assume Ω is strictly convex (with respect to the Euclidean metric), and that T 2 (K, Γ ′ ) < T < ∞. If in addition, supp f ⊂ K, then the following are true: is an elliptic pseudodifferential operator of order zero.
3. There exists C > 0 independent of a such that Proof. The first part follows directly from the second part of Theorem 3, though one may need to increase T slightly to ensure that χ = 0 near {t = T }.
To show the second part, we construct a geometric optics solution u of (2), and then microlocally construct a parametrix of (3) given h = χΛ a . We defer the proof of this to the following two lemmas.
The principal symbol of a factors as (τ + c(x)|η|)(−τ + c(x)|η|). Our geometric optics solutions will then depend on the existence of two phase functions solving the corresponding eikonal equations. To simplify the proof of the next two lemmas, we will assume that these equations are solvable on the interval [0, T ], i.e., that no caustics form. Afterward, we will remove this assumption. where π : T * x R n → T * x Γ is the tangential projection.
Proof. We will search for a parametrix u for (2) of the form where φ ± are phase functions and A σ are order zero amplitudes with asymptotic expansion where A σ j (t, y, η) is homogeneous of degree −j in η. For u to be a parametrix, we must construct φ σ and A σ so that a u ∈ C ∞ , u| t=0 = f, and ∂ t u| t=0 = −af . To satisfy the first requirement, we calculate a u: For the residual term to be smooth, the integrand must decay rapidly with respect to η. If we impose the eikonal equation and assume for the moment that a solution exists for t ∈ [0, T ], then evidently I 2 = 0.
We recognize the first term of I 1 as the transport operator X σ applied to To control the decay of I 1 + I 1 in η, we write A σ as its asymptotic expansion and consider the coefficients separately. On A σ 0 , we can impose the transport equation The subsequent terms can be taken to satisfy the lower order transport equations These equations reduce to ordinary differential equations along characteristics -the characteristics again being unit speed geodesics -whenever the eikonal equations are solvable, provided we establish Cauchy data at t = 0. To determine the proper Cauchy data to impose, we return to the Cauchy data of (2).
The condition u| t=0 = f implies that when t = 0, and accordingly at t = 0, The condition ∂ t u| t=0 = −af implies that when t = 0, Therefore at t = 0, These two equations yield the following system of boundary conditions at t = 0: Note that the resulting system of equations may be solved recursively, yielding boundary data for each coefficient of the asymptotic expansion of A ± . In particular, A + 0 (0, y, η) = A − 0 (0, y, η) = 1 2 . With this set of boundary data, the recurring system of transport equations stated above may be solved. This completes the construction of u.
Next we perform a similar construction for the backward problem (3), in order to analyze A a χΛ a .
Under these considerations, we search for v solving (3) up to smooth error with the above initial data. We assume v is of the form v(t, y) = (2π) −n e iψ(t,y,η) B(t, y, η)f (η) dη, as suggested by the forward problem.
Propagation of singularities implies that the wavefront set of v should be near the two null bicharacteristics whose projection onto T * (t0,x0) ((0, T ) × Γ) is (τ 0 , ξ 0 ). The assumption of strict convexity implies that there are two distinct geodesics whose tangent at x 0 projects to ξ 0 ; one pointing inward, and one pointing outward.
The other geodesic is the reflection of this one along Γ, according to the laws of geometric optics. Singularities near it propagate further, perhaps reflecting off the boundary at most finitely many more times (as the boundary is strictly convex), until they intersect {t = T }. Since we specified that the Cauchy data is zero there, these reflected broken geodesics carry no singularities.
Since only the null bicharacteristic corresponding to γ y0,η0 carries the singularities of h, we need only construct ψ and B in a small conic neighborhood of this bicharacteristic. The boundary data is of the form v| U = ρχΛ + a f = (2π) −n e iφ + ρχA +f dη.
For v to be a parametrix, we must have This holds in particular if we impose −∂ t v = c(x)|∇ x ψ|. Then ψ and φ + are given by the method of characteristics in a small conic neighborhood of the bicharacteristic. They have the same characteristics and the same boundary data, so in this neighborhood they are equal. Similarly, B solves the same recurring system of transport equations as A + , though with different initial data; on U , B = ρχA + . Since the first transport equation is homogeneous, we have, We repeat the construction for points in T * ((0, T ) × Γ ′ ) with τ 0 < 0, as mentioned before. Passing to a microlocal partition of unity, the above constructs v as a parametrix for (3) with boundary data h = χΛ a f and zero Cauchy data at {t = T }. Restricting v to {t = 0} yields a representation of A a χΛ a as a pseudodifferential operator, whose principal symbol is B 0 . After accounting for both contributions to B 0 (0, y 0 , η 0 ) from singularities propagating along γ y0,η0 and γ y0,−η0 , we see that the principal symbol is Therefore A a χΛ a is of order zero.
Remark. The previous two lemmas assumed that the eikonal equation was solvable on (0, T ) × Ω. In general, solutions with Cauchy data on {t = 0} exist only on a small interval [0, t 1 ]. To continue past {t = t 1 }, we use the previous parametrix to obtain Cauchy data on this hyperplane. Then, we solve the eikonal equations with boundary data φ ± 1 | t=t1 = y · η, and the transport equations with Cauchy data given by the previous parametrix. This will be solvable on some interval (t 1 , t 2 ]. By compactness, there exists a positive lower bound on the length of the interval on which the eikonal equation is solvable. Therefore, after repeating the construction at most finitely many times (when T < ∞), we obtain a parametrix on (0, T ) × Ω. The backward parametrix is then also constructed in a similar way on the same intervals.
We now conclude the proof of Theorem 4. According to Lemma 4, A a χΛ a is a pseudodifferential operator. Its symbol is of order zero by construction. If K is visible from Γ ′ , then the principal symbol comprises two terms; both are nonnegative, and at least one is positive. Hence A a χΛ a is elliptic. If K is not visible, then A a χΛ a is still a pseudodifferential operator, but it is only elliptic at those points of T * Ω \ 0 which are visible from the boundary.
Finally, we note that even if K is not visible from Γ ′ , A a χΛ a remains a pseudodifferential operator, and is still Fredholm by the considerations of Theorem 3.

A priori estimate
We return now to the proof of Proposition 1. We follow [15,Theorem 4.2], using the method of multipliers to develop an a priori estimate. There, they studied operators of the form ∂ 2 t − A, with A uniformly elliptic and in divergence form. For our purposes here we modify their proof to account for the lower order term in the damped wave equation. Recall u is a solution of        a u = F in (0, T ) × Ω, u| t=0 = f, ∂ t u| t=0 = g u| (0,T )×Γ = h subject to the compatibility condition h| t=0 = f | Γ .
Proof. The proof has three steps. First, we establish an a priori estimate of the form above, except that it also depends on the L 1 (0, T ; L 2 (Γ))-norm of the normal derivative of u, along with ||u|| H 1 and ||∂ t u|| L 2 . Next, we remove the first dependency by estimating the normal derivative. At this point the estimate contains a term depending on the L 1 (0, T ; L 2 (Γ))-norm of ∇u, which we separate into its normal and tangential components. The final dependencies are removed via Gronwall's inequality.
Let n be the outward normal vector field of Γ. As Γ is smooth, we can extend n to a smooth vector field in a neighborhood of Γ. Using a cut-off function, we can then extend it to a smooth vector field on a neighborhood of Ω, which we will also write as n. The result does not depend on the particular extension chosen.
The second term can be simplified with the divergence theorem. After integrating from 0 to s and freely applying Young's inequality, we have This agrees with the analogous estimate [15, 4.12]. Here C is a generic constant that is independent of a.
derivative in terms of the boundary data and the L 1 (0, T ; L 2 (Ω))-norms of c −1 (∂ t u) and |∇u|. Notice it is only these last two terms that depend on ||a|| ∞ , and that the dependence is at most linear.
Step 3 : We begin from the inequality derived from applying the estimate for the normal derivative to the estimate from step 1.
Finally, we take the supremum over s ∈ [0, T ] to obtain the estimate of Proposition 1. While we assumed that Ω, c, and a were all smooth, the proof shows that their being C 1 suffices.