A family of inversion formulas in Thermoacoustic Tomography

We present a family of closed form inversion formulas in thermoacoustic tomography in the case of a constant sound speed. The formulas are presented in both time-domain and frequency-domain versions. As special cases, they imply most of the previously known filtered backprojection type formulas.


Introduction and statement of main results
Thermoacoustic tomography (TAT) has recently attracted considerable attention as a promising method of biomedical imaging [28,24,20,21,16]. We give here a brief introduction to the mathematical model of TAT.
An electromagnetic (EM) pulse in visible light or radiofrequency range is sent through a biological object of interest. A fraction of the EM energy is absorbed by the tissues throughout the object. It is known that cancerous tissues absorb much more EM energy than healthy ones (e.g., [28,16]). Thus, knowing the absorption distribution function a(x) would be useful for detecting cancerous tissues. The EM energy absorbed by the tissue causes its thermoelastic expansion, which in turn releases a pressure (ultrasound) wave u(x, t) propagating through the body. This pressure is then measured by transducers located on an observation surface S surrounding the object. The initial pressure f (x) = u(x, 0) is roughly proportional to a(x). One now concentrates on the recovery of f (x) from the measured data g := u| S×[0,∞) . There has been a number of papers devoted to this reconstruction (see reviews in [13,11,17,2]). All known closed form inversion formulas assume that the object is acoustically homogeneous, i.e. the sound speed is constant. We will assume that, after choosing appropriate units, the sound speed is equal to one. The standard mathematical model of TAT is [26,7,22]    u tt (x, t) − △u(x, t) = 0, x ∈ IR n , t ≥ 0, u(x, 0) = f (x), u t (x, 0) = 0, u(y, t) = g(y, t), for y ∈ S, t ≥ 0. (1) Here g(x, t) is the measured data and f is the function to be reconstructed. In other words, one needs to invert the operator T : f → g. Since T is known to be invertible from the left only, different inversion formulas exist. In this text, we develop a family of explicit closed form inversion formulas in the case when S is the unit sphere centered at the origin and f ∈ C ∞ 0 (B), where B is the open unit ball enclosed by S.
We will also need to deal with some other operators closely related to T . Consider the wave equation    u tt (x, t) − △u(x, t) = 0, x ∈ IR n , t ≥ 0, u(x, 0) = 0, u t (x, 0) = f (x), u(y, t) = g(y, t), for y ∈ S, t ≥ 0. (2) One can now define the the operator P that maps function f into g: P(f ) = g.
Then P is related to T as follows: We also introduce the spherical Radon transform R S with centers on S by R S (f )(y, t) = S n−1 f (y + rω)r n−1 dσ(ω), for all (y, t) ∈ S × IR + .
Finally, M S is the spherical means operator with centers on S: In these formulas, dσ(ω) is the standard surface measure on the unit sphere S n−1 ⊂ IR n and ω n is the total measure of the sphere. In various works cited below, different operators from this list were considered. However, due to known explicit connections between T and P, R S , and M S (e.g., [6,15]), inversion formulas for these operators are closely related.
The first such formulas, in odd dimensions, were obtained by Finch, Patch and Rakesh [10] using some trace identities for the wave operator. Xu and Wang [27] derived a different formula for n = 3 by working in the frequency domain. A formula for all dimensions, which coincides with that of [27] when n = 3, was presented by Kunyansky [18]. Its derivation is based upon some symmetry relation for special functions. By applying the same method as in [10], Finch, Haltmeier and Rakesh [9] obtained inversion formulas for even dimensions, which involves the data measured for an infinite time period. The authors of [9] also derived another type of inversion formulas for even n, which uses only the data measured for a finite period of time (which we will refer as "finite-time formulas for even dimensions").
In spite of availability of several types of closed form inversion formulas, some questions have remained unanswered. For instance, it was not clear, what is the relation, if any, between formulas of [9,10] and [18,27], (which are known to be not equivalent outside the range of the operator T ). The same applies to the two types of inversion formulas derived in [9] for even dimensions.
The goal of this work is to obtain a unified family of closed form inversion formulas, which in particular would produce formulas of [10,9,18,27] as particular cases.
We now formulate the main results of the paper. Consider the function is the Hankel function of the first kind. Then (e.g., [1]) Φ(x, y, λ) = G(|x − y|, λ) is the solution for the Helmholtz equation obtained by limiting absorption. For odd n and any s > 0, G(s, .) ∈ C ∞ (IR). For even n, G(s, .) ∈ C ∞ (IR \ {0}) with a logarithmic singularity at zero [1]. Let g(y, t) be the function in (1) that represents the thermoacoustic data and g 0 (y, t) be its even extension with respect to t. Let alsoĝ 0 be the Fourier transform of g 0 with respect to time:ĝ Then, our family of inversion formula reads as follows: Then for any x ∈ B and ξ ∈ IR n , the following equality holds where G is the complex conjugate of G.
One can rewrite this inversion formula without going to the frequency domain. Namely, we define a transform W of a function v ∈ C ∞ [0, ∞) as follows (as long as the expressions involved make sense): Here , if n is even, , if n is odd.
(As noted in Section 7, W is a known transform that intertwines the second derivative and the Bessel operator.) One can now obtain another representation of the inner integral in (5): where W(g)(y, s) := W(g y )(s) with g y (t) := g(y, t). Then Theorem 1.1 is equivalent to Theorem 1.2. Suppose that f ∈ C ∞ 0 (B) and g = T (f ). Then for any x ∈ B and ξ ∈ IR n , the following equality holds One can look at the inner integral in (5) in a different manner yet. Namely, due to the relation between T and R S , it can be shown that where R and I are the real and imaginary parts of G. Then Theorem 1.1 implies Then for any x ∈ B and ξ ∈ IR n , the following equality holds where The key ingredient in the proof of Theorem 1.1 is the following identity, which is equivalent to a range description of operator T (see Remark 4.3): Theorem 1.4. Suppose that Ω is a ball, f ∈ C ∞ 0 (Ω), and u solves (1). Let u 0 be the even extension of u with respect to t. Then for all x ∈ Ω we have ∂Ω IR G(|x − y|, λ)û 0 (y, λ)dλdσ(y) = 0. (10) This statement implies that adding to any inversion formula an expression with an arbitrary operator A also produces an inversion formula (since the added term vanishes on the range of the operator to be inverted). In particular, if A is the operator of multiplication by an arbitrary function ϕ(x), one can add the expression to (8) and (9) to obtain the following inversion formulas: Corrolary 1.5. Suppose that f ∈ C ∞ 0 (B) and g = T (f ). Let ϕ be an arbitrary function defined on B. Then for any x ∈ B and ξ ∈ IR n , the following equalities Now, by suitable choices of ξ and ϕ in (11) and (12), one can recover the inversion formulas known in the literature. If we let ξ = x, and ϕ = −2(n − 2), −2(n − 1), or −2n in (12), we obtain correspondingly the formulas derived in [9,Theorem 4] and [10, Theorem 3]: Here ∂ t = d dt is the derivative with respect to t, and P * is the L 2 −adjoint of P. Choosing ξ = x and ϕ = 0 in (11), one gets a finite-time inversion formula for even dimensions similar to the second one in [9, Theorem 2]: where ω n is the surface area of the unit sphere S n−1 .
In particular, when n = 2, we obtain the second formula in [9, Theorem 1]: Finally, let where J n−2 2 and N n−2 2 are the Bessel and Neumann functions of order n−2 2 . Choosing ξ = 0 and ϕ = 0 in (11), one arrives at: This inversion formula is equivalent to the one obtained in [18]. Indeed, let It can be shown that h(y, s) = 2k n (y, s), and thus Proposition 1.9 implies the following result of [18]: Remark 1.11. Formulas of [10,9,18,27] do not reconstruct f (x) inside S correctly if a part of support of f lies outside S. This feature is absent in other reconstruction methods such as time reversal (the readers are referred to [17,2] for these discussion). I would be interesting to see whether formulas contained in Corollary 1.5 are any different in this regard.
The paper is organized as follows. In Section 2, we derive Theorem 1.1 from Theorem 1.4. In Section 3, we show that Theorem 1.2 equivalent to Theorem 1.1 and Theorem 1.3 is implied by Theorem 1.1. Two proofs of Theorem 1.4 are presented in Section 4. In Section 5, Propositions 1.6, 1.7, 1.9 and 1.10 are derived from Corollary 1.5 using the aforementioned choices of ξ and ϕ. Proofs of some auxiliary results are given in Section 6. Finally, some remarks are provided in Section 7.
2. Derivation of Theorem 1.1 from Theorem 1.4 Let u solves (1) and u 0 be its even extension with respect to t. Then, u 0 solves the wave equation on IR n × IR. Therefore, (19) λ 2û 0 (y, λ) + ∆û 0 (y, λ) = 0. Due to (4) and (19), we have the Green's identity: Since f (x) = u 0 (x, 0), the Fourier inversion gives Due to (20), we get Let S(0, r) be the sphere centered at the origin with radius r. For any function H ∈ C 1 (IR n ), by changing variables, we derive Applying this equality for Due to Theorem 1.4, the right hand side is zero. Thus, Combining this and (21), we get Here G s (s, λ) is the derivative of G(s, λ) with respect to s.
Applying Theorem 1.4 for Ω = B, we get Taking the derivative with respect to x of the above identity along the direction ξ, Subtracting this equality from (23), we conclude that In order to prove this lemma, we need the following explicit formula for G: Proposition 3.2. Let c n be as in (7). For any s > 0, we have This formula must be well known. For completeness, its proof is given in Section 6.
Proof of Lemma 3.1 The relation (26) means that G(s, λ) = W(e iλt )(s). Keeping this fact in mind, one sees that the following proof consists of just changing the order of the operator W with the integral sign: • For even n, due to the decay ofv and Proposition 3.2, the following calculations are valid: • For an odd n, due to Proposition 3.2 and decay ofv, The lemma is proved.
We now prove the following result, which shows that Theorem 1.1 implies Theorem 1.3: 0 (IR n ), g = T (f ) and g 0 is its even extension with respect to t. Then We need the following auxiliary result: Proof: For the expository purpose, we provide here two proofs of this proposition.
1) Letū(x, λ) = This proves the proposition. 2) Consider the transform , we also define B(g)(y, s) = B(g y )(s) where g y (s) = g(y, s). Let g = T (f ), we then have (e.g., [6,8] Substituting the expression of B into the integral and integrating by parts, we obtain Since g 0 is the even extension of g with respect to t, g 0 (y, λ) = 1 2π IR g 0 (y, t)e iλt dt = 1 π Re ∞ 0 g(y, t)e iλt dt.
Due to Proposition 3.4, we get From (27), we arrive at The proof is completed.

Proof of Theorem 1.4
In this section, we present two proofs of Theorem 1.4. One of them relies on some results of [10,9] and works in the time domain while the other one is self-contained and works in the frequency domain. Without loss of generality, we can assume that Ω = B, the open unit ball. We need to prove that if f ∈ C ∞ (B) and g = T (f ) then for all x ∈ B where g 0 is the even extension of g with respect to t. [10,9]. Due to Lemma 3.1, equality (28) is equivalent to S W(g 0 )(y, |x − y|)dσ(y) = 0, for all x ∈ B. Or, since g = g 0 | S×[0,∞) ,

The indirect proof. This proof is somewhat indirect since it uses inversion formulas in
S W(g)(y, |x − y|)dσ(y) = 0.
We introduce the following definition Then the operator P, introduced in Section 1, maps C ∞ 0 (B) into C(S × [0, ∞)). Indeed, one can show it by using the Kirchhoff-Poisson solution formulas for wave equation (e.g., [8]).
We now prove the following auxiliary result: where P * is the L 2 -adjoint of P.
Proof For an even n, a direct calculation (see [13]) gives We observe that Hence, by induction, for all k ≥ 0, Therefore, due to (30), we have For an odd n, a direct calculation (see [10]) gives This proves the lemma.
This lemma tells us that (29) is equivalent to P * (g) = 0. Equivalently, We now derive this equality from inversion formula in [10,9]. Indeed, if n is even, [9,Theorem 4] gives By subtracting these two inversion formulas we obtain the desired equality (31). If n is odd, [10,Theorem 3] gives Again, subtracting these two inversion formulas, we obtain the desired equality (31). The proof of Theorem 1.4 is completed. [12], identity (31) is the complete range description for the operator P (or, equivalently, for T ) when n is odd.
The proof is completed.

Some special cases
In this section, we derive Propositions 1.6, 1.7, 1.9 and 1.10 from Corollary 1.5 by proper choices of ξ and ϕ.  (12), we obtain Here we use the notation ∂ s for d ds .
That is, We now claim an identity, whose proof can be found in Section 6: (38) s∂ s W(g) = W(s∂ s g) − (n − 2)W(g).
Assuming this claim, we see that (37) is equivalent to Due to Lemma 4.2, we obtain • If n = −2(n − 2) then (39) becomes . Propositions 1.6 is proved (in this proof we have used the variable s in place of t).

Proposition 1.7.
Choosing ϕ(x) = 0 in (11), we obtain Hence, Since f is supported inside B, we have R S (f )(y, r) = 0 for r ≥ 2. Therefore, where, as defined in Section 1, From (40), we arrive at Choosing ξ = x, we have That is, We claim that for even n k n (y, s) = (−1) Here the integral is understood in the principal value sense. A proof of this claim will be given in Section 6.4. Assuming it, we obtain   (42), we obtain That is, This proves Proposition 1.9. In order to prove Proposition 1.10, we need only show that h(y, s) = 2k n (y, s) for any s > 0. Indeed, it is a consequence of the following lemma whose proof can be found in Section 6: This confirms (46) for odd n.

6.2.
Proof of identity (36). We prove that for all λ ∈ IR, the function is symmetric with respect to x, z ∈ B. Here we follow the line of reasoning used in [18]. Due to (41), we get where c 0 is a constant depending only on n and λ. We recall the following identities from [18] for |x| > r 0 : wherex = x |x| and Y k l is a spherical harmonic of order k. Consider the spherical harmonic expansion where α = |x| and β = |z|. Due to (49), we obtain Applying (50) and (51), we arrive at Hence, a k,l k,l (α, β) = a k,l k,l (β, α) and K(x, z, λ) = (k,l) a k,l k,l (α, β)Y k l (x)Y k l (ẑ). These two equalities give K(x, z, λ) = K(z, x, λ). The proof is completed.
6.3. Proof of identity (38). We now prove the identity We first recall from ( Here the constant c is the same in these two formulas. Recall that, where e λ (s) = e iλs . Let cos λ (s) = cos(λs) and sin λ (s) = sin(λs), we then have

Remarks
• The operator W is an intertwining operator between the second derivative and Bessel operator (e.g., [23,19,3]), i.e. Thus, W transforms the wave equation into the Darboux equation, which is known to describe spherical means (see [6,14,4,15]). • The symmetry of K(x, z, λ) is similar to that of I(x, z, λ) in [18]. It implies Theorem 1.4 as shown in Section 4.2. It is not clear that the converse is true. However, a close look at the argument in Section 4.2 shows that Theorem 1.4 implies the following weaker symmetry: IR K(x, z, λ)dλ = IR K(z, x, λ)dλ, for all x, z ∈ B.
• In [13], the authors derived an inversion formula using the known Neumann rather than Dirichlet observation data of the solution for wave solution for n = 3 (see [13,Theorem 12]). Applying the results in this paper, one can derive the inversion formulas from Neumann data for any n. Indeed, looking at (22) and (23), we arrive at: f (x) = 2 S IR G(|y − x|, λ) ∂û 0 ∂ν y (y, λ)dλdσ(y).