Photoacoustic tomography (PAT) is an emerging imaging modality that combines two types of physical fields. The domain to be imaged is illuminated with a short laser pulse that gets absorbed by the medium. The absorbed energy triggers an expansive pressure wave whose initial amplitude is proportional to the optical absorption coefficient of the tissues within the domain. The pressure waves propagate to the domain's boundary where they are measured and processed to recover the initial state of the pressure field. The advantage of this multiwave modality is that biological tissues exhibit high contrast in optical absorption whereas the acoustic waves carry high resolution. Thus, high contrast and high resolution can be achieved simultaneously which offers great potential for biomedical imaging [7,12,[46][47][48][49].
Most of the reconstruction methods for PAT assume that the acoustic properties of the medium are known. For a homogeneous wave speed, explicit formulas are available [15,18,[24][25][26]33]. Corrections of analytical formulas, valid for asymptotically small variations of sound speed, were investigated in [13,22]. For heterogeneous media, iterative methods have been designed [1,2,8,11,19,20,[36][37][38]41]. These methods also require the precise knowledge of the acoustic parameters. However, in practice the wave speed is not known exactly. For soft biological tissues, variations of the wave speed can be as great as 10% [22]. If not accounted for in the reconstruction methods, these variations cause blurring and misplacement of features in the reconstructed image.
It has been noted that the photoacoustic measurements carry information not only about the absorbed optical energy, but also about the wave speed of the medium. Based on this observation, the main question is whether the initial state of the pressure field u 0 and the wave speed c can be recovered simultaneously from a single photoacoustic measurement. In general terms, this question of uniqueness is still open. This problem, which is linear in u 0 and nonlinear in c, is very challenging. However, some progress has been made. Under certain practical assumptions, if one of the two components in (u 0 , c) is known, the other can be recovered from boundary measurements. See [40] and references therein. Liu and Uhlmann [30] gave sufficient conditions to recover both the initial pressure profile u 0 and sound speed c. More precisely, given another pair (ũ 0 , c), then c -2 u 0 = c-2 ũ0 if either ϕ = c -2 u 0 -c-2 ũ0 is a harmonic function or ϕ is independent of one variable in R 3 . Oksanen and Uhlmann studied how a modeling error in the wave speed affects the accuracy of the reconstruction of the pressure [35]. Stefanov and Uhlmann concluded that the linearized version of this problem is unstable in any scale of Sobolev spaces [39]. Kirsch and Scherzer [23] proposed an approach to simultaneously identifying the optical absorbing density and speed of sound based on a family of sectional photoacoustic illuminations and corresponding measurements. These specific illuminations must be focused on cross-sectional planes and special detectors should be employed to neglect out-of-plane signals.
Computational studies have also been carried out. Treeby et al proposed a method to select a sound speed that maximizes the sharpness of the reconstructed image [43]. Matthews et al proposed a joint reconstruction method based on an optimization framework and a low-dimensional parametrization of the sound speed [32]. Matthews and Anastasio offered an approach based on combining PAT measurements with ultrasound tomography measurements to estimate the wave speed concurrently with the pressure field [31]. A numerical invest igation was also performed by Huang et al [21]. The common conclusion from these numerical studies is that severe ill-conditioning is observed for the optimization-based methods if no regularization terms are incorporated.
Motivated by the physical scenario encountered in the PAT problem, we propose some assumptions to ensure the unique recovery of the initial pressure state and wave speed. These assumptions are stated in section 2. We show that under those conditions, for each pair (u 0 , c), there is small region of the spaces in which they reside, from which no other pair (ũ 0 , c) can induce the same acoustic measurements. Unfortunately, this region is not a neighborhood of (u 0 , c). However, the region does coincide with the conditions implicitly assumed in practical/ computational scenarios, namely, that the wave speed c can be estimated a priori by a known profile c with much more relative accuracy than the initial state of the pressure field. This result is stated in precise terms at the end of section 2 and the proof is provided in section 3. Some final remarks concerning iterative reconstruction algorithms are offered in section 4.
We consider the initial boundary value problem governed by the wave equation in a domain Ω ⊂ R d , for d = 2, 3, with smooth boundary ∂Ω. The pressure field satisfies,
on a sufficiently large window of time (0, T) where 0 < T < ∞. The impedance γ>0 models the presence of partially absorbing ultrasound sensors on the boundary ∂Ω and ∂ ν denotes the outward normal derivative. Consider also the forward mapping given by
where u solves the system (1a)-(1c). In the PAT scenario, it is common to assume that u 1 = 0. However, the mathematical analysis allows us to consider non-vanishing u 1 .
We also consider a possibly different media characterized by a wave speed c with corresponding wave field ũ that satisfies,
The wave speed c induces the definition of the corresponding forward map,
where ũ solves the system (3a)-(3c). Notice that the boundary conditions (1b) and (3b) are the same, which is implied if the media properties on ∂Ω are identical for both problems. The goal of this paper is to provide reasonable conditions on the wavespeed c, the initial state (u 0 , u 1 ), the domain Ω and time T to guarantee that the data Λ c (u 0 , u 1 ) determines the triplet (c, u 0 , u 1 ) uniquely. This cannot be done in general. Therefore, we restrict our attention to a problem satisfying certain conditions that we list as follows.
Assumption 2.1. Let the wave speeds c and c be smooth in Ω and the impedance γ be positive and smooth in ∂Ω. Moreover, let c c(x) and c c(x) for all x ∈ Ω (5)
for some c > 0. Also assume that
is the Sobolev space defined by
The next assumption is a geometric condition that ensures the observability of waves from the boundary. This is the so-called geometric control or non-trapping condition from Bardos-Lebeau-Rauch [6]. See also [9,10,16,17,27,29,44,45] for references on controllability and observability theory for hyperbolic equations.
Let Ω be a simply-connected domain with smooth boundary ∂Ω. For the Riemannian manifold (Ω, c -2 dx 2 ), let geodesic rays have finiteorder contact with the boundary. Assume there exists T o < ∞ such that any geodesic ray originating from any point in Ω at t = 0, reaches ∂Ω at a non-diffractive point before time
It will become clear that we will rely on the following restrictions on the unknown initial state of the pressure field.
Let the following bounds hold ∇u 0 2
and k ũ0
2
for some constants
Before we state the main result, we wish to comment on the relevance of these three assumptions. In the context of PAT, the assumption 2.1 is quite reasonable because, even if the actual wave speed is unknown, lower and upper bounds are readily available due to the nature of biological tissues. In other words, the types of tissue are known (muscular, granular, stromal, cancerous or fatty tissues, and blood or cerebrospinal fluid), but not their distribution within the domain of interest.
Assumption 2.2 is essential from the physical and mathematical point of view, since it ensures that the energy of the unknown initial pressure reaches the boundary, in finite time, where it can be measured. This assumption allows us to observe, in a stable manner, the initial state of the pressure field from the boundary.
Assumption 2.3 is verifiable for photoacoustic imaging because ũ1 = 0 and the initial pressure profile is given by
where G is the Grüneisen coefficient, σ is the optical absorption coefficient of the medium and I is the intensity of the probing light. The first two factors have well-known upper bounds in biological tissues. The intensity of light is chosen at the boundary and satisfies an elliptic equation in the diffusive regime. Therefore, it is also bounded above. This means that it is possible to find a finite constant K for bound (7). Now, energy must be absorbed by the medium to trigger a measurable acoustic wave. Hence, as long as the Grüneisen and absorption coefficients and the intensity of light are bounded away from zero, then it is possible to find a nonvanishing constant k for bound (8). The discrepancy in the norms employed in bounds ( 7) and ( 8) is a technicality that we were not able to avoid. In brief, assumption 2.3 means that even though the initial pressure state is unknown, its energy has known upper and lower bounds. Under these conditions, we obtain the main result of this paper.
Let assumptions 2.1 and 2.2 hold for the domain Ω, the wave speeds c and c, and the time T o < T < ∞. Let assumption 2.3 hold for the initial states of the pressure fields. There exist positive ǫ = ǫ(Ω, c, c, c, γ, k, K) so that if
where C = C(Ω, c, c, c, γ). In particular, Λ c (u 0 , u 1 )=Λ c(ũ 0 , ũ1 ) implies that u 0 = ũ0 and u 1 = ũ1 .
In the context of PAT, where u 1 = ũ1 = 0, this theorem states that given two different pairs (c, u 0 ) and (c, ũ1 ), if the relative difference between the wave speeds is much smaller than the relative difference between the initial pressure profiles, then these two pairs of data cannot induce the same acoustic measurements at the boundary. Condition (10) is illustrated in figure 1.
Notice that theorem 2.4 is not a statement of uniqueness in the usual sense. The estimate ( 11) is only valid under condition (10). This assumption does not include a neighborhood of the triplet (c, u 0 , u 1 ). It only considers a small conical region as illustrated in figure 1. No other triple (c, ũ0 , ũ1 ) in that region can produce the same boundary measurements as the triple (c, u 0 , u 1 ).
We start by setting up the initial value problem for the contrast w = u -ũ. Notice that w satisfies,
where f = c-2 -c -2 ü, w 0 = u 0 -ũ0 , w 1 = u 1 -ũ1 , and m =Λ c (u 0 , u 1 ) -Λ c(ũ 0 , ũ1 ). In terms of regularity, we have that ü ∈ H 0 ((0, T); H -1 (Ω)) because ũ is a weak solution to the wave equation as implied by the regularity of the initial conditions in assumption 2.3 [14]. Similarly, m ∈ H 1 ((0, T) × ∂Ω) and m| t=0 = 0 because (u 0 -ũ0 ) ∈ H 1 0 (Ω) as required by assumption 2.3.
By virtue of linearity, we can decompose w as follows w = w (1) + w (2) where
w (1) = w 0 and ẇ(1
and
Now we proceed to state some lemmas concerning these initial boundary value problems. In what follows, the constant C > 0 will be a generic constant that change its value from inequality to inequality. However, C does not depend on the profile of c (only on its lower and upper bounds) and does not depend on the solution to any of the boundary value problems.
The first lemma is a well-known result concerning boundary observability for hyperbolic equations. See classical references [6,9,10,29] or more recent related works [3,4,16,17,28,44,45]. Lemma 3.1 (Boundary Observability). Let w (1) solve the initial boundary value problem (13a)-(13c) and assumption 2.2 hold. Then w (1) satisfies the following boundary observability estimate,
for some positive constant C = C(Ω, c).
The next are stability results in less-regular spaces for initial boundary value problems for the wave equation using the transposition method [28]. Lemma 3.2. Let w (2)
for some constant
Illustration of condition (10) for the special case u 1 = ũ1 = 0. For a fixed pair (u 0 ,c), then a different pair (ũ 0 , c) from the pre-image of the shaded region cannot produce the same boundary measurements. The norms are understood as in theorem 2.4. The slope of the tilted line defining the shaded region is ǫ 1/2 > 0.
) and Λ c (u 0 , u 1 ) be given by (2). Then the following stability estimate holds
for some constant C = C(Ω, c, γ).
With the above lemmas, we are ready to prove the main result of this paper.
Proof of theorem 2.4. In order to apply lemmas 3.1 and 3.2, it only remains to estimate the norm of f . Recall that
so that f = gc -2 ü. Now take v ∈ H 0 ((0, T); H 1 0 (Ω)) with v H 0 ((0,T);H 1 0 (Ω))
1 and consider that
for v is arbitrary. Using the standard energy estimate [14,28] for the term ∇ũ and the assumed bound (7), we obtain
where C = C(Ω, c, c). Now, after some algebra, we obtain that
From the assumed bounds ( 5) and ( 6), we obtain a constant C = C(c, c) such that c -2
C and c-2
C. Consequently, there is another constant C = C(c, c) such that,
Plugging ( 16) into (15) and using the assumed bound (10), we obtain
where C = C(Ω, c, c). Now, combining lemmas 3.1 and 3.2, the estimate (17), and the decomposition w = w (1) + w (2) , we obtain
where C = C(Ω, c, c, c, γ), and ∂ ν w = -γ ẇ because w = u -ũ where u and ũ satisfy (1b) and (3b), respectively, on the boundary. Rearranging some terms and using lemma 3.3 and the assumed bounds ( 7) and ( 8), we get
where C = C(Ω, c, c, does not depend on ǫ or c. Therefore, for a choice ǫ<k/ (KC), we obtain the desired result. □
We end by offering some remarks concerning the practical significance of the theoretical results from the previous sections. In particular, we are interested in the successful design of iterative algorithms for the joint reconstruction of the pressure profile and wave speed. The main concern is to ensure that iterates remain within the region of uniqueness defined by (10). We seek to recover (u 0 , c) and assume that u 1 = 0 as in the case in PAT. Let (u (n) , c (n) ) for n = 0, 1, 2, ... be a sequence of iterates converging to (u 0 , c). Many of the iterative algorithms display linear convergence (such as Banach fixed-point, gradient descent, Landweber or conjugate gradient iterations [5,42]), meaning that the error behaves like
for n = 0, 1, 2, ... and for some constants 0 <α u <β u < 1 and 0 <α c <β c < 1, in the appropriate norms. Using the estimates (18a) and (18b), then we arrive at
provided that (u (n) , c (n) ) satisfies (10). Hence, the next iterate (u (n+1) , c (n+1) ) remains in the region of uniqueness defined by (10) if
We can show inductively that all the iterates (u (n) , c (n) ) for n = 1, 2, 3, ... satisfy (10) provided that (20) holds and that the initial guess (u (0) , c (0) ) satisfies (10).
In practical terms this means that for a given linear iterative algorithm, the iterates u (n) for the pressure profile must be relaxed so that they do not converge faster than the iterates c (n) for the wave speed. Or alternatively, the convergence for the wave speed iterates should be accelerated. Acceleration of linear convergence can be achieved by methods of Nesterov [34] or Aitken [42, section 5.10]. Acceleration can also be obtained with additional constraints on the wave speed. For instance, if the wave speed is assumed to belong to finite-dimensional parametric spaces, then the compactness of the projection accelerates methods such as the conjugate gradient or Landweber methods. In the computational setting, both u 0 and c belong to finite-dimensional spaces following the discretization of the domain Ω and governing differential equations. In that case, the parametrization of c should belong to a space with much lower dimension than that of the parametrization of u 0 . This can be achieved with a two-mesh approach, one mesh being much coarser than the other. Another approach was employed in [32] where each pixel value of the wave speed was assumed to belong to one of a few tissue types with each tissue type having a uniform sound speed.