It is well known that for biological tissues the attenuation of acoustic waves is frequencydependent. One way to model this attenuation is to use fractional time derivatives and consequently the representation of the propagation of ultrasound waves by integro-differential equations. Examples of this modeling are frequency power-law attenuation or fractional Szabo models (see for instance [1][2][3][4][5][6][7]) where the traveling wave may be assumed to satisfy an equation of the form γ -2 ∂ 2 t uu + β∂ k+α t u = F(t,x), for some α ∈ (0, 1), k = 1, 2, and where the fractional derivative term can be written as a convolution in time β(x)∂ k+α t u = t -∞ α (ts, x)∂ k+1 s u(s, x)ds. Assuming, as in thermoacoustics, that the wave field vanishes for negative times, and provided that the kernel is bounded and regular enough, we can perform integration by parts and write the previous integral as a convolution of u with a different kernel, plus time-derivatives of u up to order two. In the case k = 2, the sound speed is perturbed resulting in a different speed c -2 = γ -2 + β α (0), which requires conditions on β and α in order to get an effective wave speed c>0. We point out there is a recent definition for derivatives of fractional order which employs such continuous and bounded kernels [8].
In the present article, we study the inverse problem of finding the initial source f in an attenuating medium, provided boundary data u| [0,T ]×∂ and where the acoustic wave u is assumed to satisfies the system
We suppose a, b, c Then P = ∂ 2 tc 2 outside the domain of interest . The Cauchy problem associated with (1) is
since any solution of (2) extended by zero to (-∞, 0) × R n is a solution of (1). Indeed, given a smooth solution u of (2) we consider H(t)u(x, t) where H(t) is the Heaviside function. Then, we can pull out the Heaviside function from the convolution since it integrates on the interval (0, t), thus we get P (H u) = uδ ′ + 2(∂ t u)δ + auδ + (P u)H with the last term vanishing since P u = 0. For an arbitrary test function φ ∈ C ∞ c (R n+1 ) we have the following,
which is the same as problem (1).
The thermoacoustic tomography problem in a medium with convolution-type attenuation can be modeled by the following initial value problem (IVP):
where we aim to recover the initial source f from boundary measurements u| (0,T )×∂ , assuming the waves propagate freely in the space, that is, we suppose the boundary of does not interact with the outgoing waves. This last assumption has been considered for instance in [9][10][11].
The problem of thermoacoustic tomography has been broadly studied by many authors. Several reconstruction methods have been proposed for homogeneous media [12][13][14][15][16][17][18][19], and also for heterogeneous media [6,9,[20][21][22][23][24][25][26][27][28][29]. See also the reviews [30][31][32] for additional references. The theoretical analysis of the so-called time reversal method has gained considerable attention in the past few years, mainly due to the work of Stefanov and Uhlmann in [23,9]. In its initial formulation, the time reversal technique gives an approximate solution that converges to the exact one as the observation time increases. The problem of recovering the initial source for optimally short measurement time was solved in [23] for variable sound speed employing techniques from microlocal analysis.
Recently, the focus of the mathematical analysis has been placed on extensions in the following two areas. First, there is the problem of accounting for attenuating media. Homan in [10] gave a first extension of Stefanov and Uhlmann's work in this direction by considering the damped wave equation with sufficiently small damping coefficients for the time reversal method to work. In the complete data case, those results were extended to more general damping coefficients in [11]. In [33] the authors addressed the TAT problem with thermoelastic attenuation. Second, recent publications have addressed the TAT problem in enclosed domains to model the interaction of acoustic waves with reflectors and sensors. The advantages of working with this setting is that it naturally allows to consider partial data and the inverse problem is closely related with boundary control theory. See for instance [26][27][28][29].
This article falls in the first group. As far as the authors know, the TAT inverse problem with attenuation of integral type and variable sound speed has not been fully considered in the literature from an analytical point of view. From a heuristic point of view, some advances have been made. For the case of constant wave speed and constant coefficient of attenuation, Modgil et al. [34] designed a method based on relating the unattenuated wave field to the attenuated wave field via an integral operator and its subsequent inversion using a singular value decomposition. Treeby et al. [4,35] proposed a reconstruction based on time reversal and the k-space computational method. Attenuation compensation was achieved by separating the absorbing and dispersion terms in the wave equation, and reversing the sign of the absorbing coefficient during the time reversal. This method was modified in [6] where the coefficient of attenuation was allowed to vary within the region of interest, but the exponent of the power-law attenuation was still assumed to be constant. However, in some practical settings such as in the presence of bone and soft-tissue, the domain exhibits regions of varying power-law exponents. An appropriate method needs to be devised to avoid blurring and distortions in the reconstruction. Our work is a step in that direction, where the coefficients a, b, c and the kernel in (1) are allowed to vary, which effectively accounts for power-law attenuation of spatially varying exponent.
Considering attenuation terms of integral type brings some difficulties to the analysis on the propagation of waves. In particular, the equation is no longer reversible and local in time and consequently it is not possible to use techniques such as Tataru's unique continuation to get uniqueness, at least not in a direct way. Moreover, the microlocal properties of this type of integro-differential operators are not well understood. Nevertheless, it is possible to exploit the fact that an integral term of the sort considered here only presents a compact perturbation of the differential operator. This article can be view as a first attempt to understand the TAT problem in media with memory/attenuation coefficients that vary in space. A subsequent step would be to tackle viscoelastic models, and singular kernels as in the standard definition of fractional derivatives.
The paper is structured as follows. In the next section we set the framework under which our analysis is based, such as the well-posedness of the direct problem, the energy space of initial conditions and the hypothesis on the attenuation parameters, namely the damping coefficient and the attenuation kernel. In Section 3 we prove two uniqueness results. The first one is a sharp result on uniqueness for the thermoacoustic inverse problem assuming the distance function from the boundary allows us to foliate the interior of the domain by strictly convex surfaces. In particular we require ∂ to be strictly convex. The second main theorem of this section, which does not require convexity of the boundary, assumes that the sound speed satisfies a frequently used condition related with the convexity of the euclidean spheres in the metric induced by the sound speed. The stability question is addressed in Section 4 and we devote Section 5 to show the existence of a Neumann series reconstruction formula.
Let U ⊂ R n be an open bounded set with smooth boundary,
We say u is a generalized solution of
where •, • and (•, •) stand for the duality product of H -1 and H 1 0 , and the L 2 inner product respectively, and B(•, •) is the bilinear form given by B(u,ϕ) = (∇u, ∇ϕ) + (ac -2 u t ,ϕ)+ (bc -2 u, ϕ) + (c -2 * u, ϕ).
The well-posedness follows from Theorems 2.1 and 2.2 in [36]. We refer to the appendix for a complete proof. In our case, by finite speed of propagation we can take U to be a large ball containing to ensure we have null Dirichlet conditions.
Given a domain U ⊆ R n and a scalar function u(t, x), we define the local energy of u =[u, u t ] at time t as
In order to give problem (3) a physical sense we need to assume some conditions on the attenuation terms since such system must satisfies that its energy decreases over time. The previous is achieved for instance if a(x) ≥ 0 and the kernel is positive-definite, this is T 0 ( * y)ydt ≥ 0 for all y ∈ C([0, T ]). We then assume the following:
The previous condition guarantees the positive-definiteness of the kernel as shown in [37] and [38]. Moreover, if we define
it turns out thatis also a positive-definite kernel since it satisfies the same condition as .
An example of a kernel satisfying Condition 1 is (t, x) = q(x)e -αt , for some positive function q ∈ C(R n ) and α>0. In the recent article [8], the authors introduce a new definition for fractional derivatives whose kernel is of the form e -αt . As a consequence, the analysis carried out in this paper might be applied to fractional models of wave propagation following this new definition of fractional derivatives.
Under Condition 1 we define the extended energy functional at time τ>0, analogously as in [11], to be
where the last two terms take into account the portion of the energy that is lost due to the attenuation process. If we set U = R n , or by finite propagation speed we take U equal to any sufficiently large ball, in the interval [0, T ] the former energy functional E U is non-increasing since we get
and integrating in time we deduce that the extended functional is conserved. We adopt the same functional framework as in previous articles related to thermoacoustic tomography. The energy space H(U ) of initial conditions is defined to be the completion of
We also let H D (U ) denote the completion of C ∞ 0 (U ) under the norm
Notice that H(U ) = H D (U ) ⊕ L 2 (U ; c -2 dx) with the latter space denoting the L 2 functions under the weight c -2 dx.
Denoting by the region of interest and =[0, T ] × ∂ , we introduce the measurement operator : H( ) ∋ f → u| ∈ H 1 ( ), where u satisfies (3).
The first main result of this section, Theorem 1, is a uniqueness theorem for the full data case under a particular foliation condition. We work in this part with the more general hyperbolic operator (14) associated to a Riemannian metric g. We then, assuming g = c -2 dx 2 , provide a condition for the sound speed that guaranteed the existence of a particular foliation suitable for uniqueness. This is our second main result, Theorem 2.
Foliation conditions were introduced in [40] in the context of an inverse source problem for hyperbolic equations. In the field of travel time tomography in seismology, where one aims to recover the inner structure of the Earth by measuring travel times of seismic waves, at the beginning of the 20th century, Herglotz, and Weichert and Zoeppritz considered a special assumption on the isotropic wave speed that turned out to be a particular instance of Stefanov and Uhlmann's foliation condition (see [39, §6] and reference therein). These type of assumptions seem to be the natural conditions under which one could expect to propagate information from the exterior towards the interior of the domain. In particular, it has been applied before in the thermoacoustic setting to prove uniqueness for the inverse problem of recovering a sound speed when the initial source is known [40, §3].
Theorem 1 is a direct consequence of [41, Theorem 1], a boundary unique continuation result for hyperbolic equations with a memory term, where, for the sake of simplicity, the authors restricted their analysis to the Euclidean case. Nevertheless, in our work we need the full strength of such unique continuation so we have included a brief proof in the case of waves traveling in a general Riemannian setting. In few words, the idea introduced in [41] works as follows. For a solution with null Cauchy data at the boundary, it can be shown by compactness that it vanishes near the boundary and for small times, and then extend its zero set up to the characteristic time but still in a neighborhood of the boundary. The proof is then complete by noticing that assuming the levels sets of the distance to the boundary function are strictly convex, one can repeat the same process in a layer stripping fashion. The utilization of this strategy then makes necessary to consider full boundary data and the convexity assumption on the distance to the boundary function, which indeed gives us the internal foliation of the domain. This is the main reason of our hypotheses in Theorem 1. We point out that a similar method was also applied in [40] to obtain uniqueness for partial data and more general foliations. It was of fundamental importance in such proof the possibility of using a partial boundary unique continuation result independent of the foliation (see Proposition 2.1 in that article).
On the other hand, Theorem 2 addresses the isotropic case where we were able to remove the strong convexity requirement on the boundary by imposing a convexity condition on the sound speed, although losing optimality on the time needed to obtain uniqueness.
Theorem 1. Let be a bounded open subset of R n with ∂ smooth and strictly convex for a Riemannian metric g. Let T>0 be such that x n = dist(x, ∂ ) is a smooth function in with non-zero differential for 0 ≤ x n ≤ T and its level surfaces {x n = s}, for 0 ≤ s ≤ T , are strictly convex for the metric g as well.
Remark 2. Under the above hypothesis, this result presents an improvement in the condition imposed on T for uniqueness in the damped wave equation (T>2T 0 ( ) in [10, Theorem 3.1]).
Let Q = (0, T) × and x n = dist(x, ∂ ) the signed distance function defined in a neighborhood of the boundary and such that and ∂ are characterized respectively by x n > 0 and x n = 0. We define the following weight function
which is invariantly defined for any local coordinates (x 1 , ..., x n-1 ) in ∂ . Here α = α( , g) > 0 is sufficiently small and R, r>0 will be chosen large and close to each other. For ǫ ≥ 0we also consider the sets
By taking r close to R, the set Q(0) is a small neighborhood of {0} × ∂ inside Q.
We recall that in boundary normal coordinates, a Riemannian metric g takes the form
for α, β ≤ n -1. We denote g = ( gαβ (x)). Moreover, the strictly convexity of the level surfaces {x n = s} translates into
with κ s > 0t h e smallest eigenvalue (principal curvature) of the second fundamental form in {x n = s}, where R s = κ -1 s can be think as the largest curvature radius of {x n = s}. The analogous condition for convectors follows from the natural isomorphism ξ i = g ij (x)v j and reads
In what follows we consider the more general integro-differential operator
where u ′ = (u x , u t ), g is a Riemannian metric, and the vector-function A, the scalar-function b and the kernel are continuous functions.
Remark 3. The next two lemmas also hold if the coefficients A and b are analytic functions in t .
, and in particular f = 0 in (0).
Proof. Given a point y = (y ′ , 0) ∈ ∂ , let's consider local coordinates (U, (x 1 , ..., x n-1 )) in the boundary near y ′ . For ǫ ≥ 0we define the sets
In what follows we take r = Rδ, for some δ>0s m a l l enough, therefore
The idea is to obtain a well known local Carleman estimate for ũ and later use it, along with a partition of unity, to get an analogous estimate in Q(0).
Let's denote P = u tt -∂ j (g ij (x)∂ i u), the principal part of P . By analyzing the conjugate operator P τ = e τϕ Pe -τϕ , it is possible to deduce (after long computations) a pointwise estimate for v = e τϕ ũ of the form:
for some constant γ>0 depending on the parameter α which is chosen small enough, and with (Y, Z) a vector-valued function depending on lower order derivatives of v and vanishing in ∂Q y (0)\{ϕ = 0}. In fact, the previous follows by decomposing P τ v as the sum of two operators,
and bounding from below the inequality
Here we apply the convexity condition on the level surfaces {x n = s} in (13). By choosing then R large enough and δ small, we arrive to the estimate
Because v = e τϕ ũ, we can bound from below the left hand side of the previous inequality by similar terms but involving now the function u (and the exponential weight function). Integration over Q y (0) and the Gauss-Ostrogradskiȋ formula give us that
where dS denotes the surface measure on Ŵ y (0) = Q y (0) ∩{ϕ = 0}, and the matrix-function X 1 (x, t), the vector-function X 2 (x, t), and the scalar-function X 3 (x, t) are some continuous functions depending on ϕ and Q y (0). Using the continuity of the coefficients in the lower order terms (l.o.t) of P and noticing that
we can choose τ 0 larger if necessary and absorb the second summand in the right hand side above with the left hand side of (19). Then
The analogous inequality in the larger set Q(0) is obtained by considering a partition of unity and using the compactness of ∂ . More precisely, let now u ∈ C ∞ (Q(0)) and let {U i } i be a finite covering of the boundary such that on each U i we can define boundary local coordinates, and let {χ i } i be a finite smooth partition of unity subordinate to {U i } i . We also consider a collection of points y i ∈ U i . Then, denoting u i = χ 1/2 i u, and the measure dσ = dtdVo l (x) on Q, from the previous estimates we get
where notice [P , χ i ] are differential operators of order 1. We absorb the interior integrals with lower order derivatives of u using the left hand side and get
It follows from a density argument that the previous estimate also holds for functions in H 2 (Q) with null Cauchy data in ∂Q(0) ∩ ∂ . Let u be as in the hypothesis of the lemma. Then, u satisfies an inequality of the form (21), without the interior integral in the right hand side. Noticing that the boundary integral does not depends on τ , we let τ goes to infinity and conclude that u = 0in Q(0). ✷
The aim of the second lemma is to extend the time for which u is zero. Based again on Carleman estimates we will be able to succeed until we hit the characteristic surface associated to the principal part of P , this is the surface {(t, x) : Tt = dist(x, ∂ )}.
In a neighborhood of ∂ we define
and for γ>0we consider the sets
We denote by B(t 0 , x 0 ; r) the ball centered at (t 0 , x 0 ) and radius r for the euclidean metric. Given the following Claim. Suppose that for (t 0 , x 0 ) ∈ Q ǫ 2 , u vanishes below the level surface {ψ(x, t) = ψ(t 0 , x 0 )} near (t 0 , x 0 ), this is in Q ǫ 2 ψ(t 0 ,x 0 ) ∩ B(t 0 , x 0 ; r) for some r>0. Then, u = 0 in a neighborhood of (x 0 , t 0 ). the proof of the lemma is complete by the next argument. Let's assume that suppu ∩ Q ǫ 2 = ∅. We can find 0 <γ * ≤ γ 0 such that
The application of the claim on every contact point (t * , x * ) ∈ suppu ∩{(t, x) ∈ Q ǫ 2 : ψ(t, x) = γ * }, contradicts the choice of γ * . Consequently, we deduce that u = 0on every Q ǫ 2 γ , γ>0, and therefore u = 0in Q ǫ 2 .
It only remains to show the previous claim. Here is where Carleman estimates play a fundamental role, and as before we will consider a particular choice of weight function which needs to fulfill a pseudo-convex condition with respect to P = ∂ 2 t -∂ x j g ij ∂ x i • , in the set {(0, ξ) ∈ T * (t 0 ,x 0 ) }. Moreover, we will take it to be linear and non-increasing in time. Provided the above, it is possible to apply a pseudo-differential Carleman estimate introduced in [42] and conclude that u vanishes near (t 0 , x 0 ).
Let's consider local coordinates in ∂ near some y ∈ ∂ such that in those coordinates y = (x ′ 0 , 0). For some δ>0to be appropriately chosen, we define the following weight function
where here | •|stands for the euclidean norm and ψ as in (22). Denoting the principal symbol of P by p(t, x; θ, ξ) =-θ 2 +|ξ | 2 g , where |ξ | 2 g = g ij (x)ξ i ξ j is the norm on covectors induced by g, the pseudo-convexity condition requires to show that ϕ satisfies (1) Re{p, {p, ϕ}}(t 0 , x 0 ; 0, ξ) > 0 for all ξ = 0 such that p(t 0 , x 0 ; 0, ξ) = 0, (2) 1 iτ {p ϕ , p ϕ }(t 0 , x 0 ; 0, ξ ; τ) > 0 for all ξ = 0, τ>0 such that p ϕ (t 0 , x 0 ; 0, ξ, τ) = 0.
Here p ϕ (t 0 , x 0 ; 0, ξ ; τ) = p(t, x; θ Recall that we are working in boundary normal coordinates hence the metric g takes the form (12). The first condition is trivially fulfilled since the principal symbol p is elliptic in the set {θ = 0}. Let's use the following notation: the variable appearing in the subindex means we are differentiating with respect to such variable, for instance ϕ x ′ = ∂ x ′ ϕ and ϕ tx n = ∂ t ∂ x n ϕ. To verify the second condition we notice first that ϕ x ′ (t 0 , x 0 ) = ψ x ′ (t 0 , x 0 ) = 0, ϕ x n (t 0 , x 0 ) = ψ x n (t 0 , x 0 ) = -α and ϕ t (t 0 , x 0 ) = ψ t (t 0 , x 0 ) =-β where α>β>0. In fact,
Also, denoting δ ij the Kronecker delta,
Secondly, it is easy to check that p ϕ (t 0 , x 0 ; 0, ξ, τ) = 0i s equivalent to ξ n = 0 and |ξ ′ | 2 g = τ 2 (α 2β 2 ). Then, after some tedious computations, in the set of points (t 0 , x 0 ; 0, ξ ; τ) such that p ϕ = 0, we get
Let's recall the positive-definiteness of the second fundamental form in (13), and denote κ = min s∈[0,ǫ 2 ] κ s . By choosing δ>0small enough we obtain that
for the weighted norms
and the pseudo-differential operator E := e ǫ 2τ |D t | 2 . This operator can also be considered as the convolution operator
We would like to apply the above Carleman estimate to u and eventually deduce that u vanishes near (t 0 , x 0 ). With that in mind we need first to localize it near (t 0 , x 0 ). As in [41], in
Therefore, by choosing l 1 < 0s m a l l enough in magnitude, the set {ϕ(t, x) >l 1 } ∩{ψ(t, x) < ψ(t 0 , x 0 )} is contained in a sufficiently small vicinity of (t 0 , x 0 ). We then localize u by multiplying it with a function of the form χ(ϕ(t, x)) with χ ∈ C ∞ (R) a nondecreasing function such that
where l 1 <l 2 < 0are small enough in magnitude, then
In what follows we write χ meaning the composition χ • ϕ. Consequently, v = χu satisfies the inequality (23). We include the integral term in the estimates by noticing that
where P 1 is a differential operator of order 1 with coefficients supported in {(t, x)|ϕ(t, x) <l 2 }. Consequently
+ e -dǫτ e τϕ P(χu) 2 + e -dǫτ e τϕ (χu) 2 (1,τ ) .
(
The idea in what remains of the proof is to estimate Ee τϕ (χu) (2,τ ) by a term of the form e lτ , with l<0, and use [42,Proposition 4.1] to conclude that χu = 0i n {(t, x)|ϕ(t, x) >l}. Such estimate is obtained in exactly the same way as in the proof of Lemma 6 in [41], where everything reduces to estimate the term with the convolution since the other terms in the right hand side of the last inequality are easily bounded. For the arguments needed to conclude the claim we refer the reader to [41]. ✷ Proof of Theorem 1. Let u be a solution of P u = 0 with initial conditions [f, -af ] and such that f = 0. Due to our assumption on the coefficients of P , u solves (∂ 2 t -)u = 0i n (0, T) × (R n \ ) with null initial and Dirichlet boundary data. Then, for any x 0 ∈ R n \ , u vanishes in (0, T) × V for some small neighborhood V of x 0 such that V ∩ =∅. The previous is a consequence of a sharp domain of dependence for the wave operator in the exterior problem (see [12,Proposition 2]). Then u = 0i n (0, T) × (R n \ ) which implies null Neumann data, ∂u ∂ν (0,T )×∂ = 0.
Let's set
Note that ūt (t, x) = u(t, x) and ∂ t = . Moreover
We integrate equation (3) on the interval (0, t) for any t>0. It follows from the previous computations that ū solves a system of the form (15) with vanishing Cauchy data. In addition, notice that ūtt = u t ∈ L 2 (Q), so using equation ( 15) we get c 2 ū ∈ L 2 (Q), which by elliptic regularity implies ū ∈ H 2 (Q). We can now apply Lemma 2 on ū and conclude that u = 0i n a set of the form {(t, x) ∈ Q : x n <ǫ0 <t<Tx n }. This implies we have reduced the problem to the smaller domain [0, Tǫ] ×{x ∈ : x n >ǫ}. If ǫ = T we are done, otherwise we can apply again Lemma 2 in the new domain. Iterating this process we conclude the result. ✷
There is a common condition appearing in the literature of Carleman estimates and inverse problems related to the wave equation with variable sound speed. It assumes the existence of some x 0 ∈ R n for which
In geometric terms, (27) says that the spheres with center at x 0 are strictly convex for the metric c -2 dx 2 [40, §3]. Such collection of spheres can then be used to foliate the domain and, as you will see in the next corollary (see also Fig. 1), it allows us to prove unique continuation and consequently uniqueness for the inverse problem without the assumption of and the level surfaces of the distance function, dist(•, ∂ ), being strictly convex. The price we pay by removing the convexity requirement on is the lost of sharpness in the bound of T that guarantee uniqueness. Let ⊂ R n be an open and bounded subset with ∂ smooth, and T>0. We assume the sound speed c(x) satisfies condition (27) and assume the constant c 0 > 0i s a lower bound for the sound speed. Let's denote
and
Theorem 2. Assume , T and c are as above, and as in the TAT problem, we assume We denote r ={x ∈ R b :|xx 0 | = r} the sphere of center x 0 and radius r, and we set r 0 = max{0, c 0 T -D }. By hypothesis, r with r ∈[r 0 , R ] are strictly convex surfaces for the metric c -2 dx 2 that foliate (see Fig. 1(a)). For a given r ∈[r 0 , R ], let's assume that
Since r is strictly convex we can apply Lemma 2 with replaced by B(x 0 , r), the Euclidean ball of center x 0 and radius r, and deduce that u = 0in
for some ǫ>0. Recalling that c 0 is a lower bound for c, we have that
therefore we can find r 1 ∈ (0, r) such that u vanishes in the smaller set
Moreover, u has null Cauchy data on r 1 for all t ∈ (0, Tc -1 0 (Rr 1 )) (see Fig. 1(b)). If we denote by s the infimum of the radius r ≥ r 0 for which u has vanishing Cauchy data in (0, Tc -1 0 (Rr)) × r , by the first paragraph and the previous argument we know s<R (since T>0). Moreover, if s>r 0 , it must also satisfies the same property, this is, u = ∂ ν u = 0 in s for all t ∈ (0, Tc -1 0 (Rs)). Consequently, we can still apply the arguments in the paragraph above which leads us to conclude s = r 0 .
Let now f ∈ H D,a ( ) be as in the hypothesis, and u solution of (3). Analogously to the proof of Theorem 1, the function ū defined in (25) satisfies a system of the form (15) with null Cauchy data. The result then follows directly from the previous. ✷
The stability with complete data follows directly from the analogous results for the damped and undamped case. Due to the microlocal nature of this property, the minimum time needed to recover f in a stable way is usually larger than the uniqueness time. Indeed, it's necessary to capture information coming from every singularity of the initial source. In a non-trapping domain, such lower bound is related to the value
being 1 2 T 1 when there is no damping coefficient and exactly T 1 for the damped case. Notice that T 1 > 2T 0 and in the case c satisfies (27), T 1 /2 ≤ (Rr )/(αc 0 ) with (see [44, Proposition 7.1])
Theorem 3. Let be strictly convex for the metric g = c -2 dx 2 . Assume that and T are as in Theorem 1 (or as in Theorem 2). In addition, assume T 1 ( ) <T <∞ if a = 0 and 1 2 T 1 ( ) < T<∞ otherwise (resp. 2α -1 c -1 0 D <T <∞ and α -1 c -1 0 D <T <∞). Then there exists C>0 such that
Proof. The idea is to compare the observation operator with its analogous for the undamped and damped case, 0 and a respectively. These last two operators are known to be stable maps (see [9] and [10] respectively) and furthermore, from the results of the previous section, we know is injective. The proof then reduces to show that the respective error operators are compact. We only show this for the case a ≡ 0, the proof when there is a damping coefficient is obtained analogously.
From [9] follows there is a constant C>0 such that
Let's denote R = -0 and u the attenuated wave related with . Then, R maps f ∈ H D ( ) to the boundary data w| (0,T )×∂ of the system
By finite propagation speed we can work in a larger domain ′ such that w = u = 0o n its boundary and outside ′ . Due to the higher regularity theorem in [45, §7.2.3 Theorem 5], since F(t, x) =-[ * u](t, x) satisfies F, F t ∈ L 2 ((0, T); L 2 ( ′ )), we obtain that w ∈ C((0, T); H 2 ( ′ )) and w t ∈ C((0, T); H 1 ( ′ )), and consequently the trace of w in ∂ belongs to H 3/2 ((0, T) × ∂ ), with the latter space compactly embedded in H 1 ((0, T) × ∂ ).
The stability inequality is obtained by recalling the injectivity of from Theorem 1 (respectively Theorem 2) and applying the classical result [46,Proposition V.3.1]. ✷
We aim to construct a Neumann series that allow us to recover f in (3) from boundary measurements as it has been done in [23,9,[26][27][28] for the unattenuated case, and in [10,11] for the damped wave equation. However, due to the convolution term we need to modified the equation satisfied by the time reversed wave. Considering the same equation in the backward direction would imply the knowledge of the future. The strategy then is to solve a time reversal problem in such a way that the initial energy of the error function is bounded by the total energy (kinetic, potential and energy lost by attenuation) of the forward wave, inside the domain and at time T , analogously as the argument presented in [11]. Such total energy in the whole space has the attribute of being conserved in time, fact that allows us to reduce the proof to an estimate involving the norm of the initial source and the energy of the forward wave outside (see Proposition 1). The estimate says that at time T a significant portion of the energy lies outside the domain. It was first used in [9] and subsequently applied in [11].
Let's introduce the following convolution-type operator
which is the adjoint operator of * (•) under the L 2 inner product in (0, T), this is, for any
Indeed, denoting by χ I the indicator function in the interval
Following the same approach than the latest results in reconstruction for TAT in the enclosure case as well as in the attenuated case for the damped wave equation, the idea is to consider the right back projection system that will make the error operator to be a contraction. In the same way as in the proof of uniqueness, instead of working with u we set ū(t, x) = t 0 u(s, x)ds, and (t, x) as in (7).
Then, they satisfy
with p(x) = b(x) -(x, 0) ≥ 0. Notice we do not use (26) to obtain an equation as in (15) and we keep a derivative inside the convolution. If ¯ : L 2 ( ; c -2 dx) → H 1 ((0, T) × ∂ ) denotes the observation operator for this problem, this is ¯ f =ū| (0,T )×∂ , by well-posedness of the direct problem we have the following relation,
For the data h = ¯ f , we consider the solution v of the system
with φ the harmonic extension of h(T , •) in . Notice that problem (33) is well-posed. This is due to the convolution term that involves values of v in the interval (t, T), thus by doing the change of variables t → Tt we get an IBVP of the form (32) which is uniquely solvable. We define the Time Reversal operator by
and denote by K the error operator defined as follows,
with w =ūv, the error function that solves problem (34).
In what follows we suppose the domain is non-trapping (i.e. T 0 ( ) < ∞). The main result of this section is the next Theorem 4. Let be strictly convex for the metric g = c -2 dx 2 . Assume that and T are as in Theorem 1 (or as in Theorem 2). In addition, assume T 1 ( ) <T <∞ if a = 0 and 1 2 T 1 ( ) < T<∞ otherwise (resp. 2α -1 c -1 0 D <T <∞ and α -1 c -1 0 D <T <∞, with α as in (28)). Then A ¯ = Id -K with K L(L 2 ( ;c -2 dx)) < 1, and for any initial condition of (3) of the form f = (f, -af ) with f ∈ H D ( ), the thermoacoustic inverse problem has a reconstruction formula given by
Proof. Notice the error function w =ūv satisfies the equation
with ūT =ū(T , •) and ∂ t ūT = ∂ t ū(T , •). Moreover, we can write
We want to estimate the norm of Kf , hence we need to compute the energy of w. Multiplying (34) by 2c -2 w t and integrating over (0, T) × we obtain
Neglecting the integration in the spatial variable in the last two terms for a moment, we can use the identity (31) which makes them cancel each other out. Furthermore, it follows from the same identity and Condition (6) on the kernels (which guarantees positive-definiteness) that
In consequence we get
The choice of the time reversal system (34) helps to minimize the total energy in the dynamic satisfied by the error function w in a similar way as the functions φ helps to minimize the energy of w at time T . Indeed, by integration by parts we have that
From the above relations ( 35) and ( 36), we deduce
where recall the term in the right hand side is the extended energy functional associated to (3) and defined in (8). By conservation of the extended energy in R n ,
The conclusion of the theorem follows from the next proposition which is known to hold when there is no integral term.
Proposition 1. There is C>0 so that for all f ∈ L 2 ( ; c -2 dx) and ū solutions of (32),
An inequality of this form was first proved in [9, Proposition 5.1] (see (5.15) in the same article) for the case of the unattenuated wave equation, and later extended to the damped case in [11, Proposition 2] requiring a larger lower bound for the measurement time though. Such estimate is obtained by microlocalizing near the singularities and studying how their energy is transmitted across the boundary provided they hit the boundary in a transversal way. By considering strictly convex domains we can be sure that all singularities meet that requirement. When there is no damping coefficient the analysis of the singularities can be decoupled to those following the positive sound speed and negative sound speed. The time needed then for the estimate to hold equals the time needed to get at least one signal from each singularity of the initial condition, this is T> 1 2 T 1 ( ). In contrast, the appearance of a damping term makes no longer possible such microlocal decoupling, and therefore it makes necessary to wait until both signals, issued from every singularity of the initial condition, reach the boundary, or in other words T>T 1 ( ).
Let's prove the above proposition. Denote by Ū(x, t) the solution of the damped wave equation
Denoting f =[0, f ] ∈ H( ), from the paragraph above follows there is C>0so that
Furthermore, defining W = Ū -ū we obtain
and letting ū(t) =[ū(t), ūt (t)], W(t) =[ W(t), Wt (t)], the previous inequality implies
where the error function W satisfies the IVP
We claim the bounded map L 2 ( ; c -2 dx) ∋ f → ū(T ) ∈ H 1 ( c ) ⊗ L 2 ( c ) is injective. In fact, it can be decomposed as the composition of two injective bounded maps, the first one being the observation operator ¯ , which is injective since (32) is equivalent (following the computation in ( 26)) to a system of the form (15) where the method used to prove Theorem 1 (resp. Theorem 2) can be applied, and our choice of T> 1 2 T 1 ≥ T 0 (resp. T>α -1 c -1 0 R ≥ 1 2 T 1 ). The second map is the exterior IBVP map that takes Dirichlet boundary data h ∈ H 1 (0) ([0, T ] × ∂ ) to v(T ) ∈ H 1 ( c ) ⊗ L 2 ( c ), where v solves:
tc 2 ) v(t,x) = 0,( t , x ) ∈ (0,T)× R n \ v| t=0 = 0, ∂ t v| t=0 = 0, v| [0,T ]×∂ = h. (41) To see the injectivity of the latter map, consider h ∈ H 1 (0) ([0, T ] × ∂ ) such that v(T ) =v t (T ) = 0, with v solution of (41). By domain of dependence and reversibility in time of the exterior problem, we have that v vanishes in {(t, x) ∈ (0, ∞) × R n \ ¯ : dist e (x, ∂ ) >t} and also in {(t, x) ∈ (0, ∞) × R n \ ¯ : dist e (x, ∂ ) > |T -t|}. Therefore v = 0i n{(t, x) ∈ (0, 3T/2) × R n \ ¯ : dist e (x, ∂ ) > T /2}. Applying Tataru's unique continuation theorem on any p ∈{dist e (x, ∂ ) >T/2}, we deduce that v = 0i n(R n \ ¯ ) ∩{(t, x) ∈ (0, ∞) × R n :|x -p|+|t -3T/4| < 3T/4}, which implies that h vanishes for t ∈ (T /2, T)(see Fig. 2). We can now apply the same argument replacing T by T/2 and get that h is null in the interval (T /4, T/2). Iterating this process we finally conclude that h = 0for all t ∈ (0, T).
Our second claim is that the map
is compact. It is in fact a composition of the bounded maps L 2 ( ; c -2 dx) ∋ f → ūt ∈ L 2 ((0,T); L 2 (R n )), ūt → W(T ) ∈ H 2 ( c ) ⊗ H 1 ( c ) and the compact embedding
The continuity of the second map for those Sobolev spaces is due to [45, §7.2.3 Theorem 5] since denoting F := * ū t , then F, F t ∈ L 2 ((0, T); L 2 ( c )). It follows from [46,Proposition V.3.1] that for a different constant
The proposition is then proved by recalling the finite speed of propagation and applying Poincare's inequality on a large ball minus . ✷
We conclude the proof of Theorem 4 by joining (37), (38) and Proposition 1, hence for some C>1,