Nonlinear quantitative photoacoustic tomography with two-photon absorption
NNonlinear quantitative photoacoustic tomography withtwo-photon absorption
Kui Ren ∗ Rongting Zhang † October 4, 2018
Abstract
Two-photon photoacoustic tomography (TP-PAT) is a non-invasive optical molec-ular imaging modality that aims at inferring two-photon absorption property of het-erogeneous media from photoacoustic measurements. In this work, we analyze aninverse problem in quantitative TP-PAT where we intend to reconstruct optical coef-ficients in a semilinear elliptic PDE, the mathematical model for the propagation ofnear infra-red photons in tissue-like optical media with two-photon absorption, fromthe internal absorbed energy data. We derive uniqueness and stability results on thereconstructions of single and multiple optical coefficients, and present some numericalreconstruction results based on synthetic data to complement the theoretical analysis.
Key words.
Photoacoustic tomography (PAT), two-photon PAT (TP-PAT), two-photon absorp-tion, hybrid inverse problems, semilinear diffusion equation, numerical reconstruction.
AMS subject classifications 2010.
Two-photon photoacoustic tomography (TP-PAT) [35, 36, 50, 52, 55, 56, 57, 58, 59] is animaging modality that aims at reconstructing optical properties of heterogeneous mediausing the photoacoustic effect resulted from two-photon absorption. Here by two-photonabsorption we mean the phenomenon that an electron transfers to an excited state aftersimultaneously absorbing two photons whose total energy exceed the electronic energy bandgap. The main motivation for developing two-photon PAT is that two-photon optical ab-sorption can often be tuned to be associated with specific molecular signatures, such asin stimulated Raman photoacoustic microscopy, to achieve label-free molecular imaging. ∗ Department of Mathematics and ICES, University of Texas, Austin, TX 78712; [email protected] . † Department of Mathematics, University of Texas, Austin, TX 78712; [email protected] . a r X i v : . [ m a t h . A P ] A p r herefore, TP-PAT can be used to visualize particular cellular functions and molecularprocesses inside biological tissues.The principle of TP-PAT is the same as that of the regular PAT [12, 18, 38, 53], exceptthat the photoacoustic signals in TP-PAT are induced via two-photon absorption in additionto the usual single-photon absorption. In TP-PAT, we send near infra-red (NIR) photonsinto an optically absorbing and scattering medium, for instance a piece of biological tissue,Ω ⊆ R n ( n ≥ u ( x ), solvesthe following semilinear diffusion equation: −∇ · γ ( x ) ∇ u ( x ) + σ ( x ) u ( x ) + µ ( x ) | u | u ( x ) = 0 , in Ω u ( x ) = g ( x ) , on ∂ Ω (1)where the function g ( x ) models the incoming NIR photon source, the function γ ( x ) is the dif-fusion coefficient of the medium, σ ( x ) is the usual single-photon absorption coefficient of themedium, and µ ( x ) is the intrinsic two-photon absorption coefficient. The total two-photonabsorption coefficient is given by the product µ ( x ) | u | where the absolute value operation istaken to ensure that the total two-photon absorption coefficient is non-negative, a propertythat needs to be preserved for the diffusion model (1) to correctly reflect the physics.The medium absorbs a portion of the incoming photons and heats up due to the absorbedenergy. The heating then results in thermal expansion of the medium. The medium coolsdown after the photons exit. This cooling process results in contraction of the medium. Theexpansion-contraction of the medium generates ultrasound waves. The process is called thephotoacoustic effect. The initial pressure field generated by the photoacoustic effect can bewritten as [25] H ( x ) = Γ( x ) (cid:104) σ ( x ) u ( x ) + µ ( x ) | u | u ( x ) (cid:105) , x ∈ Ω . (2)where Γ is the Gr¨uneisen coefficient that describes the efficiency of the photoacoustic effect.This initial pressure field generated by single-photon and two-photon absorption processesevolves, in the form of ultrasound, according to the classical acoustic wave equation [11, 25].The data we measure in TP-PAT are the ultrasound signals on the surface of the medium.From these measured data, we are interested in reconstructing information on the opticalproperties of the medium. The reconstruction is usually done in two steps. In the firststep, we reconstruct the initial pressure field H in (2) from measured data. This step is thesame as that in a regular PAT, and has been studied extensively in the past decade; see,for instance, [4, 14, 16, 24, 28, 31, 32, 33, 40, 48] and references therein. In the second stepof TP-PAT, we attempt to reconstruct information on the optical coefficients, for instance,the two-photon absorption coefficient µ , from the result of the first step inversion, i.e. theinternal datum H in (2). This is called the quantitative step in the regular PAT [3, 9, 11,17, 26, 37, 39, 42, 44, 45, 46, 60].Due to the fact that the intrinsic two-photon absorption coefficient is very small, it isgenerally believed that events of two-photon absorption in biological tissues can only happenwhen the local photon density is sufficiently high (so that the total absorption µ | u | is largeenough). In fact, the main difficulty in the development of TP-PAT is to be able to measurethe ultrasound signal accurate enough such that the photoacoustic signal due to two-photon2bsorption is not completely buried by noise in the data. In recent years, many experimentalresearch have been conducted where it is shown that the effect of two-photon absorption canbe measured accurately; see, for instance, the study on the feasibility of TP-PAT on variousliquid samples in [56, 57, 58] (solutions), [35, 58] (suspensions) and [36] (soft matter).Despite various experimental study of TP-PAT, a thorough mathematical and numericalanalysis of the inverse problems in the second step of TP-PAT is largely missing, not tomention efficient reconstruction algorithms. The objective of this study is therefore to pursuein these directions. In the rest of the paper, we first recall in Section 2 some fundamentalmathematical results on the properties of solutions to the semilinear diffusion equation (1).We then develop in Section 3 the theory of reconstructing the absorption coefficients. InSection 4 we analyze the linearized problem of simultaneously reconstructing the absorptioncoefficients and the diffusion coefficient. Numerical simulations are provided in Section 5to validate the mathematical analysis and demonstrate the quality of the reconstructions.Concluding remarks are offered in Section 6. To prepare for the study of the inverse coefficient problems, we recall in this section somegeneral results on the semilinear diffusion model (1). Thanks to the absolute value operatorin the quadratic term µ | u | u in the equation, we can follow the standard theory of calculusof variation, as well as the theory of generalized solutions to elliptic equations in diver-gence form, to derive desired properties of the solution to the diffusion equation that wewill need in the following sections. The results we collected here are mostly minor modifica-tions/simplifications of classical results in [2, 5, 23, 27]. We refer interested readers to thesereferences, and the references therein, for more technical details on these results.We assume, in the rest of the paper, that the domain Ω is smooth and satisfies the usualexterior cone condition [27]. We assume that all the coefficients involved are bounded in thesense that there exist positive constants θ ∈ R and Θ ∈ R such that0 < θ ≤ Γ( x ) , γ ( x ) , σ ( x ) , µ ( x ) ≤ Θ < ∞ , ∀ x ∈ ¯Ω . (3)Unless stated otherwise, we assume also that( γ, σ, µ ) ∈ [ W , ( ¯Ω)] , and , g ( x ) is the restriction of a C ( ¯Ω) function on ∂ Ω . (4)where W , (Ω) denotes the usual Hilbert space of L (Ω) functions whose first weak derivativeis also in L (Ω). Note that we used W , (Ω) instead of H (Ω) to avoid confusion with the H we used to denote the internal data in (2).Technically speaking, in some of the results we obtained below, we can relax part of theabove assumptions. However, we will not address this issue at the moment. For convenience,we define the function f ( x , z ) and the linear operator L , f ( x , z ) = σ ( x ) z + µ ( x ) | z | z, and L u = −∇ · γ ∇ u. (5)3ith our assumptions above, it is clear that L is uniformly elliptic, and f ( x , z ) is con-tinuously differentiable with respect to z on ¯Ω × R . Moreover, f z ( x , z ) := ∂ z f ( x , z ) = σ ( x ) + 2 µ ( x ) | z | ≥ θ > ∀ z ∈ R .We start by recalling the definition of weak solutions to the semilinear diffusion equa-tion (1). We say that u ∈ W ≡ { w | w ∈ W , (Ω) and w | ∂ Ω = g } is a weak solution to (1)if (cid:90) Ω γ ( x ) ∇ u · ∇ v + σ ( x ) u ( x ) v ( x ) + µ ( x ) | u | u ( x ) v ( x ) d x = 0 , ∀ v ∈ W , (Ω) . We first summarize the results on existence, uniqueness and regularity of the solution to (1)in the following lemma.
Lemma 2.1.
Let ( γ, σ, µ ) satisfy (3) , and assume that g ∈ C ( ∂ Ω) . Then there is a uniqueweak solution u ∈ W , (Ω) such that u ∈ C α (Ω) ∩ C ( ¯Ω) for some < α < . If we assumefurther that ( γ, σ, µ ) and g satisfy (4) , then u ∈ W , (Ω) ∩ C ( ¯Ω) .Proof. This result is scattered in a few places in [2, 5] (for instance [5, Theorem 1.6.6]). Weprovide a sketch of proof here. For any function w ∈ W , we define the following functionalassociated with the diffusion equation (1): I [ w ] = (cid:90) Ω L ( x , w, Dw ) d x = (cid:90) Ω (cid:20) γ |∇ w | + 12 σw + 13 µ | w | w (cid:21) d x . It is straightforward to verify that I [ w ] : W → R is strictly convex (thanks again to theabsolute value in the third term) and differentiable on W with I (cid:48) [ w ] v = (cid:90) Ω (cid:104) γ ( x ) ∇ w · ∇ v + σ ( x ) wv + µ ( x ) | w | wv (cid:105) d x . We also verify that the function L ( x , z, p ) satisfies the following growth conditions: | L ( x , z, p ) | ≤ C (1 + | z | + | p | ) , | D z L ( x , z, p ) | ≤ C (1 + | z | ) , | D p L ( x , z, p ) | ≤ C (1 + | p | ) , for all x ∈ Ω, z ∈ R and p ∈ R n . It then follows from standard results in calculus ofvariations [2, 5, 23] that there exists a unique u ∈ W satisfies I [ u ] = min w ∈W I [ w ] , and u is the unique weak solution of (1). By Sobolev embedding, when n = 2 ,
3, thereexists q > n , such that u ∈ L q (Ω). This then implies that f ( x , u ) ∈ L q/ (Ω) with theassumption (3). Let us rewrite the diffusion equation (1) as − ∇ · ( γ ∇ u ) = f ( x , u ) , in Ω , u = g, on ∂ Ω . (6)Following standard results in [23, 27], we conclude that f ∈ L q/ (Ω) implies u ∈ C α (Ω)for some 0 < α <
1, where α = α ( n, Θ /θ ). Moreover, when g ∈ C ( ∂ Ω), u ∈ C ( ¯Ω). Ifwe assume further that ( γ, σ, µ ) and g satisfy (4), then f ∈ W , thanks to the fact that u ∈ C ( ¯Ω). Equation (6) then implies that u ∈ W , (Ω) ∩ C ( ¯Ω) [23, 27].4e now recall the following comparison principle for the solutions to the semilineardiffusion equation (1). Proposition 2.2. (i) Let u, v ∈ W , (Ω) ∩ C ( ¯Ω) be functions such that L u + f ( x , u ) ≤ and L v + f ( x , v ) ≥ in Ω , and u ≤ v on ∂ Ω . Then u ≤ v in Ω . (ii) If, in addition, Ω satisfies the exterior cone condition or u, v ∈ W , (Ω) , then either u ≡ v or u < v .Proof. For t ∈ [0 , u t = tu + (1 − t ) v and define a ( x ) = (cid:90) f z ( u t , x ) dt . It is thenstraightforward to check that a ( x ) ≥ θ > f z ≥ θ > u ∈ C ( ¯Ω) and v ∈ C ( ¯Ω), we conclude that u t is bounded from above when t ∈ [0 , a ( x ) ≤ Λ < ∞ for some Λ >
0. We also verify that f ( u, x ) − f ( v, x ) = a ( x )( u − v ).Let w = u − v , we have, from the assumptions in the proposition, that L w + a ( x ) w ≤ , in Ω , w ≤ , on ∂ Ω . Since L + a is uniformly elliptic, by the weak maximum principle for weak solutions [27,Theorem 8.1], w ≤ u ≤ v in Ω.If we assume in addition that u, v ∈ W , (Ω), we can use the strong maximum principleto conclude that w ≡ w (0) = 0 for some x ∈ Ω. Therefore, either w ≡
0, in which case u = v , or w <
0, in which case u < v . If u, v ∈ W , (Ω) and Ω satisfies the exterior conecondition, we can use [27, Theorem 8.19] to draw the same conclusion.The above comparison principle leads to the following assertion on the solution to thesemilinear diffusion equation (1). Proposition 2.3.
Let u j be the solution to (1) with boundary condition g j , j = 1 , . As-sume that γ , σ , µ and { g j } j =1 satisfy the assumptions in (3) and (4) . Then the followingstatements hold: (i) if g j ≥ , then u j ≥ ; (ii) sup Ω u j ≤ sup ∂ Ω g j ; (iii) if g > g , then u ( x ) > u ( x ) ∀ x ∈ Ω .Proof. (i) follows from the comparison principle in Proposition 2.2 and the fact that u ≡ g = 0. (ii) By (i), u j ≥
0. Therefore f ( x , u j ) ≥
0. Therefore, we can have −∇ · ( γ ∇ u j ) = − f ( x , u j ) ≤ , in Ω . By the maximum principle, sup Ω u j ≤ sup ∂ Ω g j . (iii) is a direct consequence of part (ii) ofProposition 2.2.In the study of the inverse problems in the next sections, we sometimes need the solutionto the semilinear diffusion equation to be bounded away from 0. We now prove the followingresult. Theorem 2.4.
Let u be the solution to (1) generated with source g ≥ ε > for some ε .Then there exists ε (cid:48) > such that u ≥ ε (cid:48) > . roof. We follow the arguments in [1]. We again rewrite the PDE as −∇ · γ ∇ u = − f ( x , u ) , in Ω , u = g, on ∂ Ω . Then by classical gradient estimates, see for instance [29, Proposition 2.20], we know thatthere exists
K >
0, depending on γ , |∇ γ | and Ω, such that | u ( x ) − u ( x ) | ≤ K | x − x | , ∀ x ∈ Ω , x ∈ ∂ Ω . Using the fact that g ≥ ε , we conclude from this inequality that there exists a d > u ( x ) ≥ ε/ , ∀ x ∈ Ω \ Ω d , where Ω d = { x ∈ Ω : dist( x , ∂ Ω) > d } . Therefore, sup Ω d/ u ≥ ε/ c ( x ) = σ ( x ) + µ ( x ) | u ( x ) | . Due to the fact that u is nonnegative and bounded fromabove, we have that 0 < θ ≤ c ( x ) ≤ Θ(1 + sup ∂ Ω | g | ). We then have that u solves −∇ · γ ∇ u + cu = 0 , in Ω , u = g, on ∂ Ω . By the Harnack inequality (see [27, Corollary 8.21]), we have that there exists a constant C , depending on d , γ , c , Ω, and Ω d/ , such that C inf Ω d/ u ≥ sup Ω d/ u. Therefore, inf Ω d/ u ≥ ε C . The claim then follows from inf Ω u ≥ min { inf Ω d/ u, inf Ω \ Ω d u } ≥ ε { /C, } ≡ ε (cid:48) .We conclude this section by the following result on the differentiability of the datum H with respect to the coefficients in the diffusion equation. This result justifies the linearizationthat we perform in Section 4. Proposition 2.5.
The datum H defined in (2) generated from an illumination g ≥ on ∂ Ω , viewed as the map H [ γ, σ, µ ] : ( γ, σ, µ ) (cid:55)→ Γ( σu + µ | u | u ) W , (Ω) × L ∞ (Ω) × L ∞ (Ω) → W , (Ω) (7) is Fr´echet differentiable when the coefficients satisfies (3) and (4) . The derivative at ( γ, σ, µ ) in the direction ( δγ, δσ, δµ ) ∈ W , (Ω) × L ∞ (Ω) × L ∞ (Ω) is given by H (cid:48) γ [ γ, σ, µ ]( δγ ) H (cid:48) σ [ γ, σ, µ ]( δσ ) H (cid:48) µ [ γ, σ, µ ]( δµ ) = Γ σv + 2 µuv δσu + 2 µ | u | v δσv + 2 µ | u | v + δµ | u | u , (8) where v j ( ≤ j ≤ ) is the solution to the diffusion equation − ∇ · ( γ ∇ v j ) + ( σ + 2 µ | u | ) v j = S j , in Ω , v j = 0 , on ∂ Ω (9) with S = ∇ · δγ ∇ u, S = − δσu, S = − δµ | u | u. roof. We show here only that u is Fr´echet differentiable with respect to γ , σ and µ . Therest of the result follows from the chain rule.Let ( δγ, δσ, δµ ) ∈ W , (Ω) × L ∞ (Ω) × L ∞ (Ω) be such that ( γ (cid:48) , σ (cid:48) , µ (cid:48) ) = ( γ + δγ, σ + δσ, µ + δµ ) satisfies the bounds in (3). Let u (cid:48) be the solution to (1) with coefficients ( γ (cid:48) , σ (cid:48) , µ (cid:48) ), anddefine (cid:101) u = u (cid:48) − u . We then verify that (cid:101) u solves the following linear diffusion equation −∇ · ( γ ∇ (cid:101) u ) + (cid:2) σ + µ ( u + u (cid:48) ) (cid:3)(cid:101) u = ∇ · δγ ∇ u (cid:48) − δσu (cid:48) − δµu (cid:48) , in Ω (cid:101) u = 0 , on ∂ Ωwhere we have used the fact that u ≥ u (cid:48) ≥ g ≥ ∂ Ω). Note also that both u and u (cid:48) are bounded from above by Proposition 2.3. Therefore, σ + µ ( u + u (cid:48) ) is bounded from above. Therefore, we have the following standard estimate [27] (cid:107) (cid:101) u (cid:107) W , (Ω) ≤ C (cid:0) (cid:107) δγ ∇ u (cid:48) (cid:107) L (Ω) + (cid:107) δσu (cid:48) (cid:107) L (Ω) + (cid:107) δµu (cid:48) (cid:107) L (Ω) (cid:1) ≤ C (cid:48) ( (cid:107) δγ (cid:107) L ∞ (Ω) + (cid:107) δσ (cid:107) L ∞ (Ω) + (cid:107) δµ (cid:107) L ∞ (Ω) ) . (10)Let (cid:101)(cid:101) u = u (cid:48) − u − ( v + v + v ) with v , v and v solutions to (9). Then we verify that (cid:101)(cid:101) u satisfies the equation −∇ · ( γ ∇ (cid:101)(cid:101) u ) + (cid:2) σ + 2 µu (cid:3)(cid:101)(cid:101) u = ∇ · δγ ∇ (cid:101) u − δσ (cid:101) u − δµ ( u (cid:48) + u ) (cid:101) u, in Ω (cid:101)(cid:101) u = 0 , on ∂ ΩTherefore, we have the following standard estimate (cid:107) (cid:101)(cid:101) u (cid:107) W , (Ω) ≤ C (cid:0) (cid:107) δγ ∇ (cid:101) u (cid:107) L (Ω) + (cid:107) δσ (cid:101) u (cid:107) L (Ω) + (cid:107) δµ (cid:101) u (cid:107) L (Ω) (cid:1) ≤ C (cid:48) (cid:0) (cid:107) δγ (cid:107) L ∞ (Ω) (cid:107)∇ (cid:101) u (cid:107) L (Ω) + (cid:107) δσ (cid:107) L ∞ (Ω) (cid:107) (cid:101) u (cid:107) L (Ω) + (cid:107) δµ (cid:107) L ∞ (Ω) (cid:107) (cid:101) u (cid:107) L (Ω) (cid:1) . (11)We can thus combine (10) with (11) to obtain the bound (cid:107) (cid:101)(cid:101) u (cid:107) W , (Ω) ≤ C (cid:0) (cid:107) δγ (cid:107) L ∞ (Ω) + (cid:107) δσ (cid:107) L ∞ (Ω) + (cid:107) δµ (cid:107) L ∞ (Ω) (cid:1) . This concludes the proof.We observe from the above proof that differentiability of H with respect to σ and µ can be proven when viewed as a map L ∞ (Ω) × L ∞ (Ω) → L ∞ (Ω), following the maximumprinciples for solutions (cid:101) u and (cid:101)(cid:101) u . The same thing can not be done with respect to γ since wecan not control the term (cid:107)∇ · δγ ∇ u (cid:48) (cid:107) L ∞ (Ω) with (cid:107) δγ (cid:107) L ∞ (Ω) without much more restrictiveassumptions on δγ . We now study inverse problems related to the semilinear diffusion model (1). We firstconsider the case of reconstructing the absorption coefficients, assuming that the Gr¨uneisencoefficient Γ and the diffusion coefficient γ are both known .7 .1 One coefficient with single datum We now show that with one datum set, we can uniquely recover one of the two absorptioncoefficients.
Proposition 3.1.
Let Γ and γ be given. Assume that g ≥ ε > for some ε . Let H and (cid:101) H be the data sets corresponding to the coefficients ( σ, µ ) and ( (cid:101) σ, (cid:101) µ ) respectively. Then H = (cid:101) H implies ( u, σ + µ | u | ) = ( (cid:101) u, (cid:101) σ + (cid:101) µ | (cid:101) u | ) provided that all coefficients satisfy (3) . Moreover, wehave (cid:107) ( σ + µ | u | ) − ( (cid:101) σ + (cid:101) µ (cid:101) | u | ) (cid:107) L ∞ (Ω) ≤ C (cid:107) H − (cid:101) H (cid:107) L ∞ (Ω) , (12) for some constant C .Proof. The proof is straightforward. Let w = u − (cid:101) u . We check that w solves − ∇ · ( γ ∇ w ) = −
1Γ ( H − (cid:101) H ) , in Ω , w = 0 , on ∂ Ω . (13)Therefore H = (cid:101) H implies w = 0 which is simply u = (cid:101) u . This in turn implies that Hu = (cid:101) H (cid:101) u ,that is σ + µ | u | = (cid:101) σ + (cid:101) µ | (cid:101) u | . Note that the condition g ≥ ε > u, (cid:101) u ≥ ε (cid:48) > ε (cid:48) following Theorem 2.4. This makes it safe to take the ratios H/u and (cid:101) H/ (cid:101) u , and toomit the absolute values on u and (cid:101) u .To derive the stability estimate, we first observe that | ( σ + µ | u | ) − ( (cid:101) σ + (cid:101) µ | (cid:101) u | ) | = 1Γ | Hu − (cid:101) H (cid:101) u | = | H ( (cid:101) u − u ) + ( H − (cid:101) H ) u Γ u (cid:101) u | . Using the fact that u and (cid:101) u are both bounded away from zero, and the triangle inequality,we have, for some constants c and c , (cid:107) ( σ + µ | u | ) − ( (cid:101) σ + (cid:101) µ | (cid:101) u | ) (cid:107) L ∞ (Ω) ≤ c (cid:107) (cid:101) u − u (cid:107) L ∞ (Ω) + c (cid:107) H − (cid:101) H (cid:107) L ∞ (Ω) . (14)On the other hand, classical theory of elliptic equations allows us to derive, from (13), thefollowing bound, for some constant c , (cid:107) u − (cid:101) u (cid:107) L ∞ (Ω) ≤ c (cid:107) H − (cid:101) H (cid:107) L ∞ (Ω) . (15)The bound in (12) then follows by combining (14) and (15).The above proof provides an explicit algorithm to reconstruct one of σ and µ from onedatum. Here is the procedure. We first solve − ∇ · ( γ ∇ u ) = − H, in Ω , u = g, on ∂ Ω (16)for u since Γ and γ are known. We then reconstruct σ as σ = H Γ u − µ | u | , (17)8f µ is known, or reconstruct µ as µ = H Γ u | u | − σ | u | , (18)if σ is known.The stability estimate (12) can be made more explicit when one of the coefficients in-volved is known. For instance, if µ is known, then we have | σ − (cid:101) σ | = 1Γ (cid:12)(cid:12) Hu − µ | u | − (cid:0) (cid:101) H (cid:101) u − µ | (cid:101) u | (cid:1)(cid:12)(cid:12) = 1Γ (cid:12)(cid:12) ( H ( (cid:101) u − u ) + ( H − (cid:101) H ) uu (cid:101) u − µ ( | u | − | (cid:101) u | ) (cid:12)(cid:12) . This leads to, using the triangle inequality again, (cid:107) σ − (cid:101) σ (cid:107) L ∞ (Ω) ≤ c (cid:48) (cid:107) (cid:101) u − u (cid:107) L ∞ (Ω) + c (cid:48) (cid:107) H − (cid:101) H (cid:107) L ∞ (Ω) . Combining this bound with (15), we have (cid:107) σ − (cid:101) σ (cid:107) L ∞ (Ω) ≤ C (cid:48) (cid:107) H − (cid:101) H (cid:107) L ∞ (Ω) , (19)for some constant C (cid:48) . In the same manner, we can derive (cid:107) µ − (cid:101) µ (cid:107) L ∞ (Ω) ≤ C (cid:48)(cid:48) (cid:107) H − (cid:101) H (cid:107) L ∞ (Ω) , (20)for the reconstruction of µ if σ is known in advance. We see from the previous result that we can reconstruct σ + µ | u | when we have one datum.If we have data generated from two different sources g and g , then we can reconstruct σ + µ | u | and σ + µ | u | where u and u are the solutions to the diffusion equation (1)corresponding to g and g respectively. If we can choose g and g such that | u | − | u | (cid:54) = 0almost everywhere, we can uniquely reconstruct the pair ( σ, µ ). This is the idea we have inthe following result. Proposition 3.2.
Let Γ and γ be given. Let ( H , H ) and ( (cid:101) H , (cid:101) H ) be the data sets corre-sponding to the coefficients ( σ, µ ) and ( (cid:101) σ, (cid:101) µ ) respectively that are generated with the pair ofsources ( g , g ) . Assume that g i ≥ ε > , i = 1 , , and g − g ≥ ε (cid:48) > for some ε and ε (cid:48) . Then ( H , H ) = ( (cid:101) H , (cid:101) H ) implies ( σ, µ ) = ( (cid:101) σ, (cid:101) µ ) provided that all coefficients involvedsatisfy (3) . Moreover, we have (cid:107) σ − (cid:101) σ (cid:107) L ∞ (Ω) + (cid:107) µ − (cid:101) µ (cid:107) L ∞ (Ω) ≤ (cid:101) C (cid:16) (cid:107) H − (cid:101) H (cid:107) L ∞ (Ω) + (cid:107) H − (cid:101) H (cid:107) L ∞ (Ω) (cid:17) , (21) for some constant (cid:101) C . roof. Let w i = u i − (cid:101) u i , i = 1 ,
2. Then w i solves − ∇ · ( γ ∇ w i ) = −
1Γ ( H i − (cid:101) H i ) , in Ω , w i = 0 , on ∂ Ω . (22)Therefore H i = (cid:101) H i implies u i = (cid:101) u i and σ + µ | u i | = (cid:101) σ + (cid:101) µ | u i | . Collecting the results for both data sets, we have (cid:18) | u | | u | (cid:19) (cid:18) σµ (cid:19) = (cid:18) | u | | u | (cid:19) (cid:18) (cid:101) σ (cid:101) µ (cid:19) . (23)When g and g satisfy the requirements stated in the proposition, we have u − u ≥ ε (cid:48) > ε (cid:48) . Therefore, the matrix (cid:18) | u | | u | (cid:19) is invertible. We can then remove this matrixin (23) to show that ( σ, µ ) = ( (cid:101) σ, (cid:101) µ ).To get the stability estimate in (21), we first verify that( σ − (cid:101) σ ) + ( µ − (cid:101) µ ) | u i | = H i u i − (cid:101) H i (cid:101) u i − (cid:101) µ ( | u i | − | (cid:101) u i | ) , i = 1 , . This leads to, (cid:18) | u | | u | (cid:19) (cid:18) σ − (cid:101) σµ − (cid:101) µ (cid:19) = H u − (cid:101) H (cid:101) u − (cid:101) µ ( | u | − | (cid:101) u | ) H u − (cid:101) H (cid:101) u − (cid:101) µ ( | u | − | (cid:101) u | ) . Therefore, we have (cid:18) σ − (cid:101) σµ − (cid:101) µ (cid:19) = (cid:18) | u | | u | (cid:19) − H ( (cid:101) u − u ) + ( H − (cid:101) H ) u u (cid:101) u − (cid:101) µ ( | u | − | (cid:101) u | ) H ( (cid:101) u − u ) + ( H − (cid:101) H ) u u (cid:101) u − (cid:101) µ ( | u | − | (cid:101) u | ) . It then follows that (cid:107) σ − (cid:101) σ (cid:107) L ∞ (Ω) + (cid:107) µ − (cid:101) µ (cid:107) L ∞ (Ω) ≤ c (cid:16) (cid:107) H − (cid:101) H (cid:107) L ∞ (Ω) + (cid:107) H − (cid:101) H (cid:107) L ∞ (Ω) + (cid:107) u − (cid:101) u (cid:107) L ∞ (Ω) + (cid:107) u − (cid:101) u (cid:107) L ∞ (Ω) (cid:17) . (24)Meanwhile, we have, from (22), (cid:107) u i − (cid:101) u i (cid:107) L ∞ (Ω) ≤ c (cid:48) (cid:107) H i − (cid:101) H i (cid:107) L ∞ (Ω) , i = 1 , . (25)The bound in (21) then follows from (24) and (25).10 Reconstructing absorption and diffusion coefficients
We now study inverse problems where we intend to reconstruct more than the absorptioncoefficients. We start with a non-uniqueness result on the simultaneous reconstructions ofall four coefficients Γ, γ , σ , and µ . (Γ , γ, σ, µ ) Let us assume for the moment that γ / ∈ C (Ω). We introduce the following Liouvilletransform v = √ γu. (26)We then verify that the semilinear diffusion equation (1) is transformed into the followingequation under the Liouville transform:∆ v − (cid:18) ∆ γ / γ / + σγ + µγ / | v | (cid:19) v = 0 , in Ω , v = γ / g, on ∂ Ω (27)and the datum H is transformed into H ( x ) = Γ( x ) (cid:18) σγ / v ( x ) + µγ v ( x ) (cid:19) . (28)Let us now define the following functionals: α = ∆ γ / γ / + σγ , β = µγ / , ζ = Γ σγ / , ζ = Γ µγ . (29)The following result says that once ( α, β, ζ ) or ( α, β, ζ ) is known, introducing new datawould not bring in new information. Theorem 4.1.
Let γ / | ∂ Ω be given and assume that γ / ∈ C (Ω) . Assume that either ( α, β, ζ ) or ( α, β, ζ ) is known, and H is among the data used to determine them. Then forany given new illumination (cid:101) g , the corresponding datum (cid:101) H is uniquely determined by ( (cid:101) g, H ) .Proof. Let us first rewrite the datum as H = ζ v + ζ v . When α and β are known, we knowthe solution v of (27) for any given g . If ζ is also known, we know also ζ v . We thereforecan form the ratio (cid:101) H − ζ (cid:101) vH − ζ v = ζ (cid:101) v ζ v = (cid:101) v v . We then find (cid:101) H as (cid:101) H = (cid:101) v v ( H − ζ v ) + ζ (cid:101) v . If ζ is not known but ζ is known, we can formthe ratio (cid:101) H − ζ (cid:101) v H − ζ v = ζ (cid:101) vζ v = (cid:101) vv . This gives (cid:101) H = (cid:101) vv ( H − ζ v ) + ζ (cid:101) v . The proof is complete.11he above theorem says that we can at most reconstruct the triplet ( α, β, ζ ) or thetriplet ( α, β, ζ ). Neither triplet would allow the unique determination of the four coefficients(Γ , γ, σ, µ ). Once one of the triplets is determined, adding more data is not helpful in termsof uniqueness of reconstructions.Similar non-uniqueness results were proved in the case of the regular PAT [8, 9]. In thatcase, it was also shown that if the Gr¨uneisen coefficient Γ is known, for instance from multi-spectral measurements [10, 39], one can uniquely reconstruct the absorption coefficient andthe diffusion coefficient simultaneously. In the rest of this section, we consider this case,that is, Γ is known, for our TP-PAT model. ( γ, σ, µ ) We study the problem of reconstructing ( γ, σ, µ ), assuming Γ is known, in linearized settingfollowing the general theory of overdetermined elliptic systems developed in [21, 47]. For thesake of the readability of the presentation below, we collect some necessary terminologiesin the theory of overdetermined elliptic systems in Appendix A. We refer interested readersto [7, 34, 54] for overviews of the theory in the context of hybrid inverse problems andreferences therein for more technical details on the theory. Our presentation below followsmainly [7].We linearize the nonlinear inverse problem around background coefficients ( γ, σ, µ ), as-suming that we have access to data collected from J different illumination sources { g j } Jj =1 .We denote by ( δγ, δσ, δµ ) the perturbations to the coefficients. Let u j be the solution to (1)with source g j and the background coefficients. We then denote by δu j the perturbation tosolution u j . Following the calculations in Proposition 2.5, we have, for 1 ≤ j ≤ J , −∇ · ( δγ ∇ u j ) − ∇ · ( γ ∇ δu j ) + δσu j + δµ | u j | u j + ( σ + 2 µ | u j | ) δu j = 0 , in Ω (30) δσu j + δµ | u j | u j + ( σ + 2 µ | u j | ) δu j = δH j / Γ , in Ω (31)To simplify our analysis, we rewrite the above system into, 1 ≤ j ≤ J , −∇ · ( δγ ∇ u j ) − ∇ · ( γ ∇ δu j ) = − δH j / Γ , in Ω (32) u j δσ + | u j | u j δµ + ( σ + 2 µ | u j | ) δu j = + δH j / Γ , in Ω (33)This is a system of 2 J differential equations for J + 3 unknowns { δγ, δσ, δµ, δu , . . . , δu J } .The system is formally overdeterminted when J > { δu j } Jj =1 are given already. They are ho-mogeneous Dirichlet conditions since g does not change when the coefficients change. Theboundary conditions for ( δγ, δσ, δµ ) are what need to be determined. In the case of single-photon PAT, it has been shown that one needs to have γ | ∂ Ω known to have uniqueness inthe reconstruction [8, 44]. This is also expected in our case. We therefore take δγ | ∂ Ω = φ for some known φ . The boundary conditions for σ and µ are given by the data. In fact, onthe boundary, u = g . Therefore, we have, from (33) which holds on ∂ Ω, that g j δσ + | g j | g j δµ = δH j / Γ , on ∂ Ω .
12f we have two perturbed data sets { δH , δH } with g and g sufficiently different, we canthen uniquely reconstruct ( δσ | ∂ Ω , δµ | ∂ Ω ): δσ | ∂ Ω = δH | g | g − δH | g | g Γ g g ( | g | − | g | ) ≡ φ , δµ | ∂ Ω = δH g − δH g Γ g g ( | g | − | g | ) ≡ φ . Therefore, we have the following Dirichlet boundary condition for the unknowns( δγ, δσ, δµ, δu , . . . , δu J ) = ( φ , φ , φ , , · · · , . (34)Let us introduce v = ( δγ, δσ, δµ, δu , . . . , δu J ), S = ( − δH , δH , . . . , − δH J , δH J ) / Γ, and φ = ( φ , φ , φ , , · · · , A ( x , D ) v = S , in Ω B ( x , D ) v = φ, on ∂ Ω (35)where A is a matrix differential operator of size M × N , M = 2 J and N = 3 + J , while B is the identity operator. The symbol of A is given as A ( x , i ξ ) = − i V · ξ − ∆ u γ | ξ | − i ξ · ∇ γ . . . u | u | u σ + 2 µ | u | . . . − i V J · ξ − ∆ u J . . . γ | ξ | − i ξ · ∇ γ u J | u J | u J . . . σ + 2 µ | u J | , (36)with V j = ∇ u j , ≤ j ≤ J and ξ ∈ S n − ( S n − being the unit sphere in R n ).It is straightforward to check that if we take the associated Douglis-Nirenberg numbersas { s i } Ji =1 = (0 , − , . . . , , − , { t j } J +3 j =1 = (1 , , , , . . . , , (37)the principal part of A is simply A itself with the − i ξ · ∇ γ and − ∆ u j (1 ≤ j ≤ J ) termsremoved.In three-dimensional case, we can establish the following result. Theorem 4.2.
Let n = 3 . Assume that the background coefficients γ ∈ C (Ω) , σ ∈ C (Ω) ,and µ ∈ C (Ω) satisfy the bounds in (3) . Then, there exists a set of J ≥ n + 1 illuminations { g j } Jj =1 such that A is elliptic. Moreover, the corresponding elliptic system ( A , B ) , withboundary condition (34) , satisfies the Lopatinskii criterion.Proof. Let us first rewrite the principal symbol A as A ( x , i ξ ) = − i V · ξ γ | ξ | . . . i V · ξ γ | ξ | ( σ + 2 µ | u | ) u u | u | u . . . − i V J · ξ . . . γ | ξ | i V J · ξ γ | ξ | ( σ + 2 µ | u J | ) u J u J | u J | u J . . . .
13t is then straightforward to check that A ( x , i ξ ) is of full-rank as long as the followingsub-matrix is of full-rank: (cid:101) A ( x , i ξ ) = i V · ξ γ | ξ | ( σ + 2 µ | u | ) u u | u | u ... ... ... i V J · ξ γ | ξ | ( σ + 2 µ | u J | ) u J u J | u J | u J . To simplify the calculation, we introduce Σ j = σ + 2 µ | u j | , (cid:98) F j = V j · ξ . We also eliminatethe non-zero common factor i γ | ξ | from the first column and u j from each row. Withoutloss of generality, we check the determinant of first 3 (since J ≥ n + 1 = 4) rows of thesimplified version of the submatrix ˜ A ( x , i ξ ). This determinant is given asdet( (cid:101) A ) = (cid:98) F Σ u ( | u | − | u | ) + (cid:98) F Σ u ( | u | − | u | ) + (cid:98) F Σ u ( | u | − | u | )= Σ Σ Σ u u u (cid:18) (cid:98) F u u ( | u | − | u | )Σ Σ + (cid:98) F u u ( | u | − | u | )Σ Σ + (cid:98) F u u ( | u | − | u | )Σ Σ (cid:19) . With the assumptions on the background coefficients, we can take u j to be the complexgeometric optics solution constructed following Theorem 6.6 (in Appendix B) for ρ j withboundary condition g j . Then we have (cid:98) F k u i u j ( | u i | − | u j | )Σ i Σ j = u i u j ( | u i | − | u j | )Σ i Σ j ∇ u k · ξ = u i u j u k ( | u i | − | u j | )Σ i Σ j ( ρ k + O (1)) · ξ . This gives us,det( (cid:101) A ) ∼ (cid:16) Σ ( | u | − | u | ) ρ + Σ ( | u | − | u | ) ρ + Σ ( | u | − | u | ) ρ (cid:17) · ξ . (38)Let us define α k = Σ k ( | u i | − | u j | ) with ( k, i, j ) ∈ { (1 , , , (2 , , , (3 , , } . Then we have α + α + α = 0. Let ( e , e , e ) an orthonormal basis for R . Then ξ = (cid:80) k =1 c k e k with | c | + | c | + | c | = 1. We take ρ = β (cid:0) τ e + i (cid:101) τ e (cid:1) , ρ = β (cid:0) τ e + i (cid:101) τ e (cid:1) , ρ = β (cid:0) τ e + i (cid:101) τ e (cid:1) , where | τ k | = | (cid:101) τ k | , ∀ ≤ k ≤
3. It is straightforward to verify that ρ k · ρ k = 0, | ρ k | = √ | τ k || β k | for all 1 ≤ k ≤
3. We now deduce from (38) that det( (cid:101) A ) ∼ (cid:3) R + i (cid:3) I where (cid:3) R = α τ β c + α τ β c + α τ β c , (cid:3) I = α (cid:101) τ β c + α (cid:101) τ c + α (cid:101) τ β c . Take β = β = β , τ k = (cid:101) τ k = 1, 1 ≤ k ≤
3. Then det( A ) (cid:54) = 0 unless c = c = c . [ This isbecause if det( A ) = 0, we have (cid:3) R = 0, (cid:3) I = 0, and a + a + a = 0. That is c c c c c c α α α = . { α k } k =1 ]. Let us now take ρ = 2 ρ . Then the submatrix formed by u , u and u will have full rank when c = c = c . Thereforethe submatrix formed by u , u , u and u is of full rank for any ξ .To show that ( A , B ) satisfies the Lopatinskii criterion given in Definition 6.4 for a set ofwell chosen u j , we first observe that since B = I , we have, from the definition in (55), { η j } J +3 j =1 = {− , · · · , − } (39)with the selection of the Douglis-Nirenberg numbers { s i } Ji =1 and { t j } J +3 j =1 in (37), and theprincipal part of B has components B , = 1 and B ,k(cid:96) = 0 otherwise. Therefore, the systemof differential equations in (56) and (57) takes the following form( V j · ζ − i V j · ν ddz ) δγ ( z ) − γ ( −| ζ | + d dz ) δu j ( z ) = 0 , z > u j δσ ( z ) + | u j | u j δµ ( z ) + ( σ + 2 µ | u j | ) δu j ( z ) = 0 , z > δγ = 0 , z = 0 (42)where γ , σ , µ , u j and V j (1 ≤ j ≤ J , are all evaluated at y ∈ ∂ Ω. We first deduce from (41)that, for 1 ≤ j ≤ J , δu j = − u j Σ j δσ − | u j | u j Σ j δµ, z > . Plugging this into (40), we obtain that, for 1 ≤ j ≤ J ,( V j · ζ − i V j · ν ddz ) δγ + γ ( −| ζ | + d dz ) (cid:18) u j Σ j δσ + | u j | u j Σ j δµ (cid:19) = 0 , z > . Without loss of generality, we consider the system formed by u , u , u . Let (cid:101) F j = V j · ν , p j = u j Σ j and q j = u j Σ j . We look for eigenvalues of the system as the root ofdet (cid:98) F − i λ (cid:101) F p γ ( λ − | ζ | ) q γ ( λ − | ζ | ) (cid:98) F − i λ (cid:101) F p γ ( λ − | ζ | ) q γ ( λ − | ζ | ) (cid:98) F − i λ (cid:101) F p γ ( λ − | ζ | ) q γ ( λ − | ζ | ) = 0 . We observe first that the above equation admits two repeated roots λ , = ±| ζ | . Besidesthat, we have another root λ = − i (cid:98) F ( p q − p q ) + (cid:98) F ( p q − p q ) + (cid:98) F ( p q − p q ) (cid:101) F ( p q − p q ) + (cid:101) F ( p q − p q ) + (cid:101) F ( p q − p q ) . Moreover, the eigenvectors corresponding to λ , are of the form π , = xy x and y arbitrary. Therefore, δγ ( z ) = ce i | λ | z . Using the boundary condition δγ (0) = 0and the decay condition δγ ( z ) → z → ∞ , we conclude that δγ ( z ) ≡
0. This in turnimplies, from (41), that δσ ( z ) ≡ δµ ( z ) ≡
0. The proof is complete.For the set of Douglis-Nirenberg numbers { s i } and { t j } in (37), as well as the parameters { η k } given in (39), we defined the function space, parameterized by (cid:96) > n + , W (cid:96) = W (cid:96) − s , (Ω) × . . . × W (cid:96) − s J , (Ω) × W (cid:96) − η − , ( ∂ Ω) × . . . W (cid:96) − η − , ( ∂ Ω) . We have the following uniqueness and stability result.
Theorem 4.3.
Under the same conditions of Theorem 4.2, let { δH j } Jj =1 and { (cid:103) δH j } Jj =1 bethe data sets generated with ( δγ, δσ, δµ ) and ( (cid:102) δγ, (cid:102) δσ, (cid:102) δµ ) respectively. Assume that the dataare such that ( S , φ ) ∈ W (cid:96) and ( (cid:101) S , (cid:101) φ ) ∈ W (cid:96) . Then there exists a set of J ≥ n + 1 boundaryilluminations, { g j } Jj =1 , such that { δH j } Jj =1 = { (cid:103) δH j } Jj =1 (resp. ( S , φ ) = ( (cid:101) S , (cid:101) φ ) ) implies ( δγ, δσ, δµ ) = ( (cid:102) δγ, (cid:102) δσ, (cid:102) δµ ) (resp. v = ˜ v ) if δγ | ∂ Ω = (cid:102) δγ | ∂ Ω . Moreover, the following stabilityestimate holds: J +3 (cid:88) j =1 (cid:107) v j − (cid:101) v j (cid:107) W (cid:96) + tj, (Ω) ≤ C (cid:16) J (cid:88) i =1 (cid:107)S i − (cid:101) S i (cid:107) W (cid:96) − si, (Ω) + (cid:88) k =1 (cid:107) φ k − (cid:101) φ k (cid:107) W (cid:96) − ηk − , ( ∂ Ω) (cid:17) , (43) for all (cid:96) > n + .Proof. We start with the uniqueness result. Let δH j = 0, 1 ≤ j ≤
3, we then have that u j δσ + | u j | u j δµ + ( σ + 2 µ | u j | ) δu j = 0 , ≤ j ≤ . We can eliminate the variables δσ and δµ to have, with E = { (1 , , , (2 , , , (3 , , } , (cid:88) ( i,j,k ) ∈E u i u j ( u j − u i )( σ + 2 µ | u k | ) δu k = 0 . (44)Let G be the Green function corresponding to the operator −∇ · γ ∇ with the homogeneousDirichlet boundary condition. We can then write (44) as, using δγ | ∂ Ω = 0 as well as δu j | ∂ Ω =0, (cid:88) ( i,j,k ) ∈E u i u j ( u j − u i )( σ + 2 µ | u k | ) (cid:90) Ω δγ ( y ) ∇ u k ( y ) · ∇ G ( x ; y ) d y = 0 . Take u k to be the complex geometric optics solution we constructed in Theorem 6.6 withcomplex vector ρ k , using the fact that u k ∼ e ρ k · x (1 + ϕ k ) (and ϕ k decays as | ρ k | − ) and ∇ u k = u k ( ρ k + O (1)), we can rewrite the above equation as, for | ρ k | sufficiently large, (cid:88) ( i,j,k ) ∈E u i u j ( u j − u i )( σ + 2 µ | u k | ) (cid:90) Ω δγ ( y ) u k ( y ) ρ k · ∇ G ( x ; y ) d y = 0 . ρ k such that (cid:60) ρ k < |(cid:60) ρ k | issufficiently large to simplify the above equation further to (cid:90) Ω δγ ( y ) v ( x ; y ) · ∇ G ( x ; y ) d y = 0 . (45)with v the vector given by v = (cid:88) ( i,j,k ) ∈E σ ( x ) (cid:16) u i u j ( u j − u i ) (cid:17) ( x ) u k ( y ) ρ k . (46)We now need the following lemma. Lemma 4.4.
Let v be such that: (i) there exists c > such that | v | ≥ c > for a.e. x ∈ Ω ;and (ii) v ∈ [ W , ∞ (Ω)] n at least. Then (45) implies that δγ ≡ .Proof. Let u be the solution to − ∇ · γ ∇ u − ∇ · ( δγ v ) = 0 , in Ω , u = 0 , on ∂ Ω (47)Then u ( x ) = (cid:90) Ω δγ ( y ) v ( x ; y ) · ∇ G ( x ; y ) d y . Therefore (45) implies that u ≡
0. Therefore − ∇ · δγ v = 0 , in Ω , δγ = 0 , on ∂ Ω . (48)This is a transport equation for δγ that admits the unique solution δγ ≡ v satisfying the assumed requirements; see for instance [8, 13, 15, 20, 30] and referencestherein.It is straightforward to check that we can select { ρ j } j =1 such that the vector field v defined in (46) satisfies the requirement in Lemma 4.4. We then conclude that δγ ≡ φ k = (cid:101) φ k = 0 when 4 ≤ k ≤ J + 3. Note also that the following simplificationcan be made in (43): J (cid:88) i =1 (cid:107)S i − (cid:101) S i (cid:107) W (cid:96) − si, (Ω) ≤ J (cid:88) i =1 (cid:107) δH j − (cid:102) δH j (cid:107) W (cid:96) +2 , (Ω) . The proof is complete. 17
Numerical simulations
We present in this section some preliminary numerical reconstruction results using syntheticinternal data. We restrict ourselves to two-dimensional settings only to simplify the compu-tation. The spatial domain of the reconstruction is the square Ω = (1 , . All the equationsin Ω are discretized with a first-order finite element method on triangular meshes. In allthe simulations in this section, reconstructions are performed on a finite element mesh withabout 6000 triangular elements. The nonlinear system resulted from the discretization ofthe diffusion equation (1) is solved using a quasi-Newton method as implemented in [43].To generate synthetic data for inversion, we solve (1) using the true coefficients. Weperformed reconstructions using both noiseless and noisy synthetic data. For the noisydata, we added random noise to the data by simply multiplying each datum by (1 + √ (cid:15) × − random ) with random a uniformly distributed random variable taking values in [1 , (cid:15) being the noise level (i.e. the size of the variance in percentage).We will focus on the reconstruction of the absorption coefficients σ and µ . We presentreconstruction results from two different numerical methods. Direct Algorithm.
The first method we use is motivated from the method of proofs ofPropositions 3.1 and 3.2. When we have J ≥ { H j } Jj =1 from J illuminations { g j } Jj =1 , we first reconstruct, for each j , u ∗ j as the solutions to −∇ · ( γ ∇ u ∗ j ) = − H ∗ j Γ in Ω , u ∗ j = g j on ∂ Ω . We then reconstruct σ + µ | u ∗ j | = H j Γ u ∗ j . Collecting this quantity from all data, we have, foreach point x ∈ Ω, | u ∗ | ... ...1 | u ∗ J | (cid:18) σµ (cid:19) = H ∗ Γ u ∗ ... H ∗ J Γ u ∗ J . We then reconstruct ( σ, µ ) by solving this small linear system, in least square sense, at eachpoint x ∈ Ω. Therefore, the main computational cost of this algorithm lies in the numericalsolution of the J linear equations for { u ∗ j } . Least-Square Algorithm.
The second reconstruction method that we will use is basedon numerical optimization. This method searches for the unknown coefficient by minimizingthe objective functionalΦ( σ, µ ) ≡ J (cid:88) j =1 (cid:90) Ω (Γ σu j + Γ µ | u j | u j − H ∗ j ) d x + κR ( σ, µ ) , (49)where we use the functional R ( σ, µ ) = (cid:0)(cid:82) Ω |∇ σ | d x + (cid:82) Ω |∇ µ | d x (cid:1) together with the pa-rameter κ to add regularization mechanism in the reconstructions. We use the BFGS quasi-18ewton method that we developed in [43] to solve this minimization problem. It is straight-forward to check, following Proposition 2.5, that the gradient of the functional Φ( σ, µ ) withrespect to σ and µ are given respectively byΦ (cid:48) σ [ σ, µ ]( δσ ) = (cid:90) Ω (cid:40) J (cid:88) j =1 (cid:104) z j Γ u j + v j u j (cid:105) δσ + κ ∇ σ · ∇ δσ (cid:41) d x (50)Φ (cid:48) µ [ σ, µ ]( δµ ) = (cid:90) Ω (cid:40) J (cid:88) j =1 (cid:104) z j Γ | u j | u j + v j | u j | u j (cid:105) δµ + κ ∇ µ · ∇ δµ (cid:41) d x (51)where z j = Γ( σu j + µ | u j | u j ) − H ∗ j and v j solves − ∇ · γ ∇ v j + ( σ + 2 µ | u j | ) v j = − z j Γ( σ + 2 µ | u j | ) , in Ω , v j = 0 , on ∂ Ω (52)Therefore, in each iteration of the optimization algorithm, we need to solve J semilineardiffusion equations for { u j } Jj =1 and then J adjoint linear elliptic equations for { v j } Jj =1 toevaluate the gradients of the objective function with respect to the unknowns.Figure 1: The true coefficients, γ (left), σ (middle), µ (right), used to generate syntheticdata for the reconstructions.Figure 2: The absorption coefficient µ reconstructed using synthetic data containing differentlevels ( (cid:15) = 0 , , , Direct Algorithm is used in thereconstructions.
Experiment I.
We start with a set of numerical experiments on the reconstruction of thetwo-photon absorption coefficient µ assuming that the single-photon absorption coefficient σ is known. We use data collected from four different sources { g j } j =1 , { H j } j =1 . We perform19econstructions using the Direct Algorithm . In Fig. 2 we show the reconstruction resultsfrom noisy synthetic data with noise levels (cid:15) = 0, (cid:15) = 1, (cid:15) = 2, and (cid:15) = 5. The truecoefficients used to generate the data are shown in Fig. 1.To measure the quality of the reconstruction, we use the relative L error. This error isdefined as the ratio between (i) the L norm of the difference between the reconstructed co-efficient and the true coefficient and (ii) the L norm of the true coefficient, expressed in per-centage. The relative L errors in the reconstructions of µ in Fig 2 are 0 . . . .
23% for (cid:15) = 0, (cid:15) = 1, (cid:15) = 2 and (cid:15) = 5 respectively.Figure 3: Same as in Fig. 2 except that the reconstructions are performed with the
Least-Square Algorithm . Experiment II.
One of the main limitations on the
Direct Algorithm is that it requiresthe use of illumination sources that are positive everywhere on the boundary. This is difficultto implement in practical applications. The
Least-Square Algorithm , however, does not havesuch requirement on the optical sources (but it is computationally more expensive). Here werepeat the simulations in Experiment I with the
Least-Square Algorithm . The reconstructionresults are shown in Fig 3. We observe that, with the same (not exactly the same since therealizations of the noise are different) data sets, the reconstructions from the two differentalgorithms are of very similar quality. The relative L errors for the reconstructions in Fig 3are 0 . . . .
36% respectively for the four cases.
Experiment III.
In the third set of numerical experiments, we study the simultaneousreconstructions of the single-photon and two-photon absorption coefficients, σ and µ . Weagain use data collected from four different sources. In Fig. 4, we show the reconstructionsfrom data containing different noise levels using the Direct Algorithm . The relative L error in the reconstructions of ( σ, µ ) are (0 . , . . , . . , . . , . (cid:15) = 0, (cid:15) = 1, (cid:15) = 2 and (cid:15) = 5. Experiment IV.
We now repeat the simulations in Experiment III with the
Least-SquareAlgorithm . The results are shown in Fig. 5. The relative L errors in the reconstructions arenow (0 . , . . , . . , . . , . (cid:15) = 0, (cid:15) = 1, (cid:15) = 2, and (cid:15) = 5. The quality of the reconstructions is slightlylower than, but still comparable to, that in the reconstructions in Experiment III.20igure 4: The absorption coefficient pair σ (top row) and µ (bottom row) reconstructedusing the Direct Algorithm with data at different noise levels ( (cid:15) = 0 , , , Least-Square Algorithm . 21e observe from the above simulation results that, in general, the quality of the recon-structions is very high. When we have the illumination sources that satisfy the positivityrequirement on the whole boundary of the domain, the
Direct Algorithm provides an effi-cient and robust reconstruction method. The
Least-Square Algorithm is less efficient butis as robust in terms of the quality of the reconstructions. The reconstructions with the
Least-Square Algorithm are done for a fixed regularization parameter that we selected by acouple of trial-error test. It is by no means the optimal regularization parameter that canbe selected through more sophisticated algorithms [22]. However, this is an issue that wethink is not important at the current stage of this project. Therefore, we did not pursuefurther in this direction.We also observe that the performances of the reconstructions of σ and µ are different inboth the direct and the least-square methods: the reconstructions of σ is visually better thanthose of µ in general. This phenomenon is not manifested in our current stability results.Further analysis, for instance to refine the stability results in Section 3, are needed to fullyunderstand this difference in numerical stability. We studied in this paper inverse problems in quantitative photoacoustic tomography withtwo-photon absorption. We derived uniqueness and stability results in the reconstruction ofsingle-photon and two-photon absorption coefficients, and proposed explicit reconstructionmethods in this case with well-selected illumination sources. We also studied the inverseproblem of reconstructing the diffusion coefficient in addition to the absorption coefficientsand obtained partial results on the uniqueness and stability of the reconstructions for thelinearized problem. We presented some numerical studies based on the explicit reconstruc-tion procedures as well as numerical optimization techniques to demonstrate the type ofquality that can be achieved in reasonably controlled environments (where noise strength inthe data is moderate).Our focus in this paper is to study the mathematical properties of the inverse problems.There are many issues that have to be address in the future. Mathematically, it would benice to generalize the uniqueness and stability results in Section 4, on multiple coefficientreconstructions in linearized settings, to the fully nonlinear problem. Computationally,detailed numerical analysis, in three-dimensional setting, need to be performed to quantifythe errors in the reconstructions in practically relevant scenarios. It is especially importantto perform reconstructions starting from acoustic data directly, following for instance theone-step reconstruction strategy developed in [19], to see how sensitive the reconstructionof the two-photon absorption coefficient is with respect to noise in the acoustic data. Onthe modeling side, it is very interesting to see if the current study can be generalized toradiative transport type models for photon propagation.22 cknowledgments
We would like to the anonymous referees for their useful comments that help us improve thequality of the paper. This work is partially supported by the National Science Foundationthrough grants DMS-1321018 and DMS-1620473.
Appendix A: Terminologies in overdetermined ellipticsystems
We recall here, very briefly, some terminologies and notations related to overdeterminedlinear elliptic systems, following the presentation in [7, 54]. Let M, (cid:102) M , N be three positiveintegers such that
M > N . we consider the following system of M differential equations for N variables { v , · · · , v N } with (cid:102) M boundary conditions: A ( x , D ) v = S , in Ω (53) B ( x , D ) v = φ, on ∂ Ω (54)Here A ( x , D ) is a matrix differential operator whose ( i, j ) element, denoted by A ij ( x , D )(1 ≤ i ≤ M , 1 ≤ j ≤ N ), is a polynomial in D for any x ∈ Ω. B ( x , D ) is a matrixdifferential operator whose ( k, (cid:96) ) element, denoted by B k(cid:96) ( x , D ) (1 ≤ k ≤ (cid:102) M , 1 ≤ (cid:96) ≤ N ),is a polynomial in D for any x ∈ ∂ Ω.We associate an integer s i (1 ≤ i ≤ M ) to row i of A and an integer t j to column j (1 ≤ j ≤ N ) of A . Definition 6.1.
We call the integers { s i } Mi =1 and { t j } Nj =1 the Douglis-Nirenberg numbersassociated to A if: (a) s i ≤ , ≤ i ≤ M ; (b) when s i + t j ≥ , the order of A ij ( x , D ) isnot greater than s i + t j ; and (c) when s i + t j < , A ij ( x , D ) = 0 . Definition 6.2.
The principal part of A , denoted by A , is defined as the part of A suchthat the degree of A ,ij ( x , D ) is exactly s i + t j . We say that A is elliptic, in the sense of Douglis-Nirenberg, if the matrix A ( x , ξ ) is ofrank N for all ξ ∈ S n − ( S n − being the unit sphere in R n ) and all x ∈ Ω.Let b k(cid:96) be the order of B k(cid:96) and define η k = max ≤ (cid:96) ≤ N ( b k(cid:96) − t (cid:96) ) , ≤ k ≤ (cid:102) M . (55)
Definition 6.3.
The principal part of B , denoted by B , is defined as the part of B suchthat the order of B ,k(cid:96) is exactly η k + t (cid:96) . Let B ( x , D ) be the principal part of B . Fix y ∈ ∂ Ω, and let ν be the inward unit normalvector at y . Let ζ ∈ S n − be a vector such that ζ · ν = 0 and | ζ | (cid:54) = 0. We consider on the23alf line y + z ν , z > A ( y , i ζ + ν ddz ) (cid:101) u ( z ) = 0 , z > , (56) B ( y , i ζ + ν ddz ) (cid:101) u ( z ) = 0 , z = 0 . (57) Definition 6.4.
If for any y ∈ ∂ Ω , the only solution to the system (56) - (57) such that (cid:101) u ( z ) → as z → ∞ is (cid:101) u ≡ , then we say that ( A , B ) satisfies the Lopatinskii criterion. It is well-established that [7, 47, 54] when ( A , B ) satisfies the Lopatinskii criterion, thesystem (53)-(54) can be solved up to possibly a finite dimensional subspace. Moreover, ageneral a priori stability estimate can be established for the system. Define the functionspace W (cid:96) = W (cid:96) − s , (Ω) × . . . × W (cid:96) − s M , (Ω) × W (cid:96) − η − , ( ∂ Ω) × . . . W (cid:96) − σ (cid:102) M − , ( ∂ Ω) , for some (cid:96) > n + . Then it can be shown that, if ( S , φ ) ∈ W (cid:96) , N (cid:88) j =1 (cid:107) v j (cid:107) W (cid:96) + tj, (Ω) ≤ C (cid:16) M (cid:88) i =1 (cid:107)S i (cid:107) W (cid:96) − si, (Ω) + (cid:102) M (cid:88) i =1 (cid:107) φ i (cid:107) W (cid:96) − ηi − , ( ∂ Ω) (cid:17) + (cid:101) C (cid:88) t j > (cid:107) v j (cid:107) L (Ω) , (58)provided that all the quantities involved are regular enough. The last term in the estimatecan be dropped when uniqueness of the solution can be proven. More details of this theorycan be found in [7] and references therein. Appendix B: CGO solutions to equation (1)
This appendix is devoted to the construction of complex geometric optics (CGO) solu-tions [49, 51] to our model equation (1). We restrict the construction to the three-dimensionalsetting ( n = 3). We start by revisiting CGO solutions to the classical diffusion problem thatwas first developed in [49]: − ∇ · ( γ ∇ u ) + σu = 0 , in Ω (59)with the assumption that γ ∈ C (Ω) and σ ∈ C (Ω). Let u ∗ be a solution to this equation,then the Liouville transform defined in (26) shows that ˜ u ∗ = √ γu ∗ solves∆˜ u ∗ − (cid:16) ∆ √ γ √ γ + σγ (cid:17) ˜ u ∗ = 0 , in Ω . (60)The following result is well-known. Theorem 6.5 ([6, 49, 51]) . Let γ ∈ C (Ω) and σ ∈ C (Ω) . For any ρ ∈ C n such that ρ · ρ = 0 and | ρ | sufficiently large, there is a function g such that the solution to (60) , withthe boundary condition ˜ u ∗ | ∂ Ω = g , takes the form ˜ u ∗ = e ρ · x (1 + ϕ ( x )) , (61)24 ith ϕ ( x ) satisfying the estimate | ρ |(cid:107) ϕ (cid:107) W , (Ω) + (cid:107) ϕ (cid:107) W , (Ω) ≤ C (cid:13)(cid:13)(cid:13)(cid:13) ∆ √ γ √ γ + σγ (cid:13)(cid:13)(cid:13)(cid:13) W , (Ω) . (62)The function ˜ u ∗ is called a complex geometric optics solution to (60) and u ∗ = γ − / e ρ · x (1 + ϕ ( x )) (63)is a complex geometric optics solution to (59). With the regularity assumption on γ and (62),it is easy to verify that ∇ u ∗ ∼ u ∗ ( ρ + O (1)) . (64)We now show, using the Newton-Kantorovich method [41], that we can construct a CGOsolution to our semilinear diffusion model that is very close, in W , (Ω), to u ∗ for some ρ . Theorem 6.6.
Let γ ∈ C (Ω) , σ ∈ C (Ω) and µ ∈ C (Ω) . Let ρ ∈ C n be such that ρ · ρ = 0 and | ρ | sufficiently large. Assume further that ρ and Ω satisfy − (cid:101) κ | ρ | ≤ (cid:60) ( ρ · x ) ≤ − κ | ρ | , ∀ x ∈ ¯Ω , (65) for some < κ < (cid:101) κ < ∞ . Then, there exists a function g such that the solution to (1) takesthe form u ( x ) = u ∗ ( x ) + v ( x ) , (66) with v such that (cid:107) v (cid:107) W , (Ω) ≤ ce − κ (cid:48) κ | ρ | , (67) for some constant c and some κ (cid:48) ∈ (1 , .Proof. We first observe that the assumption in (65) allows us to bound the CGO solution u ∗ and its gradient as (cid:107) u ∗ (cid:107) L ∞ (Ω) ≤ ce − κ | ρ | , (68) (cid:107)∇ u ∗ (cid:107) L ∞ (Ω) ≤ ce − κ | ρ | ( | ρ | + 1) . (69)Let P ( x , D ) be the differential operator defined in (1). We verify that, with the assump-tions on the coefficients involved, P (cid:48) u ∗ , the linearization of P at u ∗ , P (cid:48) u ∗ := −∇· γ ∇ + σ +2 µ | u ∗ | ,admits a bounded inverse as a linear map from W , (Ω) to W , (Ω).We observe from the construction that u ∗ is away from 0. Therefore, there exists aconstant r > B r ( u ∗ ) (in the W , (Ω) metric) contains functions that areaway from 0. Let u ∈ B r ( u ∗ ), u ∈ B r ( u ∗ ) and v ∈ W , (Ω) be given, we check that P (cid:48) u v − P (cid:48) u v = 2 µ ( | u | − | u | ) v. (70)This leads to (cid:107)P (cid:48) u v − P (cid:48) u v (cid:107) L (Ω) = (cid:107) µ ( | u | − | u | ) v (cid:107) L (Ω) ≤ c (cid:107) u − u (cid:107) W , (Ω) (cid:107) v (cid:107) W , (Ω) . (71)25e can also bound (cid:107)∇ ( P (cid:48) u v − P (cid:48) u v ) (cid:107) L (Ω) as follows. We first verify that (cid:107)∇ ( P (cid:48) u v − P (cid:48) u v ) (cid:107) L (Ω) = 2 (cid:107)∇ (cid:0) µ ( | u | − | u | ) v (cid:1) (cid:107) L (Ω) ≤ (cid:107) µ ( | u | − | u | ) ∇ v (cid:107) L (Ω) + 2 (cid:107)∇ (cid:0) µ ( | u | − | u | ) (cid:1) v (cid:107) L (Ω) ≤ (cid:107) µ ( | u | − | u | ) (cid:107) L ∞ (Ω) (cid:107)∇ v (cid:107) L (Ω) + 2 (cid:107)∇ (cid:0) µ ( | u | − | u | ) (cid:1) (cid:107) L ∞ (Ω) (cid:107) v (cid:107) L (Ω) ≤ C (cid:107) v (cid:107) W , (Ω) (cid:16) (cid:107) u − u (cid:107) L ∞ (Ω) + (cid:107)∇ (cid:0) µ ( | u | − | u | (cid:1) (cid:107) L ∞ (Ω) (cid:17) . (72)We then perform the expansion, using the fact that | u j | > j = 1 , ∇ (cid:0) µ ( | u | − | u | ) (cid:1) = ( | u | − | u | ) ∇ µ + µ (cid:60) (cid:0) u | u | ∇ (¯ u − ¯ u ) + ( u | u | − u | u | ) ∇ ¯ u (cid:1) . This gives us the bound (cid:107)∇ (cid:0) µ ( | u | − | u | ) (cid:1) (cid:107) L ∞ (Ω) ≤ c (cid:107) u − u (cid:107) L ∞ (Ω) + c (cid:107)∇ ( u − u ) (cid:107) L ∞ (Ω) . (73)We can then combine (72) with (73) and use Sobolev embedding, for instance [27, Corollary7.11], to conclude that (cid:107)∇ ( P (cid:48) u v − P (cid:48) u v ) (cid:107) L (Ω) ≤ C (cid:107) u − u (cid:107) W , (Ω) (cid:107) v (cid:107) W , (Ω) . (74)We then have, from the bounds in (71) and (74), the following bound on the operatornorm of P (cid:48) u − P (cid:48) u by (cid:107)P (cid:48) u − P (cid:48) u (cid:107) L (cid:0) W , (Ω) ,W , (Ω) (cid:1) ≤ c (cid:107) u − u (cid:107) W , (Ω) . (75)Let w be the solution to P (cid:48) u ∗ ( x , D ) w = P ( x , D ) u ∗ (such that w = ( P (cid:48) u ∗ ) − P ( x , D ) u ∗ ), thatis, − ∇ · γ ∇ w + ( σ + 2 µ | u ∗ | ) w = µ | u ∗ | u ∗ , in Ω , w = 0 , on ∂ Ω . (76)It then follows from classical elliptic theory that (cid:107) ( P (cid:48) u ∗ ) − P ( x , D ) u ∗ (cid:107) W , (Ω) ≤ c (cid:107)| u ∗ | u ∗ (cid:107) W , (Ω) ≤ ˜ ce − κ | ρ | ( | ρ | + 1) ≤ ˜˜ ce − κ (cid:48) κ | ρ | , (77)for some κ (cid:48) ∈ (1 , | ρ | is sufficientlylarge, there exists a solution to (1) in the ball of radius r (cid:48) = ˜˜ ce − κ (cid:48) κ | ρ | centered at u ∗ , in W , (Ω). The solution is of the form (66).Let us remark that the assumption (65) on ρ and Ω is feasible. Let ( e , e , e ) be a basisof R . Since Ω is bounded, we can place Ω in a box in the first octant of the ( e , e , e )system. Let us take (cid:60) ( ρ ) = a e + a e + a e for some { a k } k =1 . We check that we canhave ˜ κ |(cid:60) ( ρ ) | ≤ (cid:60) ( ρ ) ≤ κ |(cid:60) ( ρ ) | for some constant ˜ κ and κ . If we further take |(cid:61) ( ρ ) | to beproportional to |(cid:60) ( ρ ) | , we can have (65); see Section 4 for the choices of different ρ in ouranalysis. 26 eferences [1] G. Alessandrini, M. Di Cristo, E. Francini, and S. Vessella , Stability for quantitativephoto acoustic tomography with well-chosen illuminations , Annali di Matematica, 196 (2017),pp. 395–406.[2]
A. Ambrosetti and A. Malchiodi , Nonlinear Analysis and Semilinear Elliptic Problems ,Cambridge University Press, Cambridge, 2007.[3]
H. Ammari, E. Bossy, V. Jugnon, and H. Kang , Mathematical modelling in photo-acoustic imaging of small absorbers , SIAM Rev., 52 (2010), pp. 677–695.[4]
H. Ammari, E. Bretin, V. Jugnon, and A. Wahab , Photo-acoustic imaging for attenu-ating acoustic media , in Mathematical Modeling in Biomedical Imaging II, H. Ammari, ed.,vol. 2035 of Lecture Notes in Mathematics, Springer-Verlag, 2012, pp. 53–80.[5]
M. Badiale and E. Serra , Semilinear Elliptic Equations for Beginners , Springer, London,2011.[6]
G. Bal , Introduction to Inverse Problems . Department of Applied Physics and AppliedMathematics, Columbia University, 2012.[7] ,
Hybrid inverse problems and redundant systems of partial differential equations , in In-verse Problems and Applications, P. Stefanov, A. Vasy, and M. Zworski, eds., vol. 615 ofContemporary Mathematics, American Mathematical Society, 2013, pp. 15–48.[8]
G. Bal and K. Ren , Multi-source quantitative PAT in diffusive regime , Inverse Problems,27 (2011). 075003.[9] ,
Non-uniqueness result for a hybrid inverse problem , in Tomography and Inverse Trans-port Theory, G. Bal, D. Finch, P. Kuchment, J. Schotland, P. Stefanov, and G. Uhlmann, eds.,vol. 559 of Contemporary Mathematics, Amer. Math. Soc., Providence, RI, 2011, pp. 29–38.[10] ,
On multi-spectral quantitative photoacoustic tomography in diffusive regime , InverseProblems, 28 (2012). 025010.[11]
G. Bal and G. Uhlmann , Inverse diffusion theory of photoacoustics , Inverse Problems, 26(2010). 085010.[12]
P. Beard , Biomedical photoacoustic imaging , Interface Focus, 1 (2011), pp. 602–631.[13]
F. Bouchut and G. Crippa , Uniqueness, renormalization and smooth approximations forlinear transport equations , SIAM J. Math. Anal., 38 (2006), pp. 1316–1328.[14]
P. Burgholzer, G. J. Matt, M. Haltmeier, and G. Paltauf , Exact and approximativeimaging methods for photoacoustic tomography using an arbitrary detection surface , Phys. Rev.E, 75 (2007). 046706.[15]
F. Colombini, G. Crippa, and J. Rauch , A note on two-dimensional transport withbounded divergence , Comm. Partial Differential Equations, 31 (2006), pp. 1109–1115.[16]
B. T. Cox, S. R. Arridge, and P. C. Beard , Photoacoustic tomography with a limited-aperture planar sensor and a reverberant cavity , Inverse Problems, 23 (2007), pp. S95–S112. B. T. Cox, S. R. Arridge, K. P. K¨ostli, and P. C. Beard , Two-dimensional quantita-tive photoacoustic image reconstruction of absorption distributions in scattering media by useof a simple iterative method , Applied Optics, 45 (2006), pp. 1866–1875.[18]
B. T. Cox, J. G. Laufer, and P. C. Beard , The challenges for quantitative photoacousticimaging , Proc. of SPIE, 7177 (2009). 717713.[19]
T. Ding, K. Ren, and S. Vallelian , A one-step reconstruction algorithm for quantitativephotoacoustic imaging , Inverse Problems, 31 (2015), p. 095005.[20]
R. J. DiPerna and P.-L. Lions , On the Cauchy problem for Boltzmann equations: globalexistence and weak stability , Ann. Math., 130 (1989), pp. 321–366.[21]
A. Douglis and L. Nirenberg , Interior estimates for elliptic system of partial differentialequations , Comm. Pure Appl. Math., 8 (1955), pp. 503–508.[22]
H. W. Engl, M. Hanke, and A. Neubauer , Regularization of Inverse Problems , KluwerAcademic Publishers, Dordrecht, The Netherlands, 1996.[23]
L. C. Evans , Partial Differential Equations , American Mathematical Society, Providence,RI, 2010.[24]
D. Finch, M. Haltmeier, and Rakesh , Inversion of spherical means and the wave equationin even dimensions , SIAM J. Appl. Math., 68 (2007), pp. 392–412.[25]
A. R. Fisher, A. J. Schissler, and J. C. Schotland , Photoacoustic effect for multiplyscattered light , Phys. Rev. E, 76 (2007). 036604.[26]
H. Gao, S. Osher, and H. Zhao , Quantitative photoacoustic tomography , in MathematicalModeling in Biomedical Imaging II: Optical, Ultrasound, and Opto-Acoustic Tomographies,H. Ammari, ed., vol. 2035 of Lecture Notes in Mathematics, Springer, 2012, pp. 131–158.[27]
D. Gilbarg and N. S. Trudinger , Elliptic Partial Differential Equations of Second Order ,Springer-Verlag, Berlin, 2000.[28]
M. Haltmeier, T. Schuster, and O. Scherzer , Filtered backprojection for thermoacousticcomputed tomography in spherical geometry , Math. Methods Appl. Sci., 28 (2005), pp. 1919–1937.[29]
Q. Han and F. Lin , Elliptic Partial Differential Equations , American Mathematical Society,Providence, 1997.[30]
M. Hauray , On two-dimensional Hamiltonian transport equations with l p loc coefficients , Ann.IHP. Anal. Non Lin., 20 (2003), pp. 625–644.[31] Y. Hristova , Time reversal in thermoacoustic tomography - an error estimate , Inverse Prob-lems, 25 (2009). 055008.[32]
A. Kirsch and O. Scherzer , Simultaneous reconstructions of absorption density and wavespeed with photoacoustic measurements , SIAM J. Appl. Math., 72 (2013), pp. 1508–1523.[33]
P. Kuchment and L. Kunyansky , Mathematics of thermoacoustic tomography , Euro. J.Appl. Math., 19 (2008), pp. 191–224. P. Kuchment and D. Steinhauer , Stabilizing inverse problems by internal data , InverseProblems, 28 (2012).[35]
Y.-H. Lai, S.-Y. Lee, C.-F. Chang, Y.-H. Cheng, and C.-K. Sun , Nonlinear photoa-coustic microscopy via a loss modulation technique: from detection to imaging , Optics Express,22 (2014), pp. 525–536.[36]
G. Langer, K.-D. Bouchal, H. Gr¨un, P. Burgholzer, and T. Berer , Two-photonabsorption-induced photoacoustic imaging of Rhodamine B dyed polyethylene spheres using afemtosecond laser , Optics Express, 21 (2013), pp. 22410–22422.[37]
J. Laufer, B. T. Cox, E. Zhang, and P. Beard , Quantitative determination of chro-mophore concentrations from 2d photoacoustic images using a nonlinear model-based inversionscheme , Applied Optics, 49 (2010), pp. 1219–1233.[38]
C. Li and L. Wang , Photoacoustic tomography and sensing in biomedicine , Phys. Med. Biol.,54 (2009), pp. R59–R97.[39]
A. V. Mamonov and K. Ren , Quantitative photoacoustic imaging in radiative transportregime , Comm. Math. Sci., 12 (2014), pp. 201–234.[40]
L. V. Nguyen , A family of inversion formulas in thermoacoustic tomography , Inverse Probl.Imaging, 3 (2009), pp. 649–675.[41]
J. M. Ortega , The Newton-Kantorovich theorem , Amer. Math. Monthly, 75 (1968), pp. 658–666.[42]
A. Pulkkinen, B. T. Cox, S. R. Arridge, J. P. Kaipio, and T. Tarvainen , ABayesian approach to spectral quantitative photoacoustic tomography , Inverse Problems, 30(2014). 065012.[43]
K. Ren, G. Bal, and A. H. Hielscher , Frequency domain optical tomography based onthe equation of radiative transfer , SIAM J. Sci. Comput., 28 (2006), pp. 1463–1489.[44]
K. Ren, H. Gao, and H. Zhao , A hybrid reconstruction method for quantitative photoa-coustic imaging , SIAM J. Imag. Sci., 6 (2013), pp. 32–55.[45]
K. Ren, R. Zhang, and Y. Zhong , Inverse transport problems in quantitative PAT formolecular imaging , Inverse Problems, 31 (2015). 125012.[46]
T. Saratoon, T. Tarvainen, B. T. Cox, and S. R. Arridge , A gradient-based method forquantitative photoacoustic tomography using the radiative transfer equation , Inverse Problems,29 (2013). 075006.[47]
V. A. Solonnikov , Overdetermined elliptic boundary-value problems , J. Sov. Math., 1 (1973),pp. 477–512.[48]
P. Stefanov and G. Uhlmann , Thermoacoustic tomography with variable sound speed ,Inverse Problems, 25 (2009). 075011.[49]
J. Sylvester and G. Uhlmann , Global uniqueness theorem for an inverse boundary valueproblem , Ann. Math., 125 (1987), pp. 153–169. G. J. Tserevelakis, D. Soliman, M. Omar, and V. Ntziachristos , Hybrid multiphotonand optoacoustic microscope , Opt. Lett., 39 (2014), pp. 1819–1822.[51]
G. Uhlmann , Electrical impedance tomography and Calderon’s problem , Inverse Problems,25 (2009). 123011.[52]
B. E. Urban, J. Yi, V. Yakovlev, and H. F. Zhang , Investigating femtosecond-laser-induced two-photon photoacoustic generation , J. Biomed. Opt., 19 (2014). 085001.[53]
L. V. Wang , Tutorial on photoacoustic microscopy and computed tomography , IEEE J. Sel.Topics Quantum Electron., 14 (2008), pp. 171–179.[54]
T. Widlak and O. Scherzer , Stability in the linearized problem of quantitative elastography ,Inverse Problems, 31 (2015). 035005.[55]
P. W. Winter, A. G. York, D. D. Nogare, M. Ingaramo, R. Christensen, A. Chit-nis, G. H. Patterson, and H. Shroff , Two-photon instant structured illumination mi-croscopy improves the depth penetration of super-resolution imaging in thick scattering sam-ples , Optica, 1 (2014), pp. 181–191.[56]
Y. Yamaoka, M. Nambu, and T. Takamatsu , Frequency-selective multiphoton-excitation-induced photoacoustic microscopy (MEPAM) to visualize the cross sections of dense objects ,in Photons Plus Ultrasound: Imaging and Sensing, A. A. Oraevsky and L. V. Wang, eds.,SPIE, 2010. 75642O.[57] ,
Fine depth resolution of two-photon absorption-induced photoacoustic microscopy usinglow-frequency bandpass filtering , Optics Express, 19 (2011), pp. 13365–13377.[58]
Y. Yamaoka and T. Takamatsu , Enhancement of multiphoton excitation-induced pho-toacoustic signals by using gold nanoparticles surrounded by fluorescent dyes , in Photons PlusUltrasound: Imaging and Sensing, A. A. Oraevsky and L. V. Wang, eds., SPIE, 2009. 71772A.[59]
C. S. Yelleswarapu and S. R. Kothapalli , Nonlinear photoacoustics for measuring thenonlinear optical absorption coefficient , Optics Express, 18 (2010), pp. 9020–9025.[60]
R. J. Zemp , Quantitative photoacoustic tomography with multiple optical sources , AppliedOptics, 49 (2010), pp. 3566–3572., AppliedOptics, 49 (2010), pp. 3566–3572.