{'fraction_non_alphanumeric': 0.05455286081343611, 'fraction_numerical': 0.02030456852791878, 'mean_word_length': 4.094173982442139, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlinear quantitative photoacoustic tomography with two-photon absorption

Two-photon photoacoustic tomography (TP-PAT) is a non-invasive optical molecular imaging modality that aims at inferring two-photon absorption property of heterogeneous media from photoacoustic measurements. In this work, we analyze an inverse problem in quantitative TP-PAT where we intend to reconstruct optical coefficients in a semilinear elliptic PDE, the mathematical model for the propagation of near infra-red photons in tissue-like optical media with two-photon absorption, from the internal absorbed energy data. We derive uniqueness and stability results on the reconstructions of single and multiple optical coefficients, and present some numerical reconstruction results based on synthetic data to complement the theoretical analysis.


Introduction
Two-photon photoacoustic tomography (TP-PAT) [35,36,50,52,55,56,57,58,59] is an imaging modality that aims at reconstructing optical properties of heterogeneous media using the photoacoustic effect resulted from two-photon absorption. Here by two-photon absorption we mean the phenomenon that an electron transfers to an excited state after simultaneously absorbing two photons whose total energy exceed the electronic energy band gap. The main motivation for developing two-photon PAT is that two-photon optical absorption can often be tuned to be associated with specific molecular signatures, such as in stimulated Raman photoacoustic microscopy, to achieve label-free molecular imaging. Therefore, TP-PAT can be used to visualize particular cellular functions and molecular processes inside biological tissues.
The principle of TP-PAT is the same as that of the regular PAT [12,18,38,53], except that the photoacoustic signals in TP-PAT are induced via two-photon absorption in addition to the usual single-photon absorption. In TP-PAT, we send near infra-red (NIR) photons into an optically absorbing and scattering medium, for instance a piece of biological tissue, Ω ⊆ R n (n ≥ 2), where they diffuse. The density of the photons, denoted by u(x), solves the following semilinear diffusion equation: in Ω u(x) = g(x), on ∂Ω (1) where the function g(x) models the incoming NIR photon source, the function γ(x) is the diffusion coefficient of the medium, σ(x) is the usual single-photon absorption coefficient of the medium, and µ(x) is the intrinsic two-photon absorption coefficient. The total two-photon absorption coefficient is given by the product µ(x)|u| where the absolute value operation is taken to ensure that the total two-photon absorption coefficient is non-negative, a property that 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 absorbed energy. The heating then results in thermal expansion of the medium. The medium cools down after the photons exit. This cooling process results in contraction of the medium. The expansion-contraction of the medium generates ultrasound waves. The process is called the photoacoustic effect. The initial pressure field generated by the photoacoustic effect can be written as [25] H(x) = Γ(x) σ(x)u(x) + µ(x)|u|u(x) , x ∈ Ω. (2) where Γ is the Grüneisen coefficient that describes the efficiency of the photoacoustic effect. This initial pressure field generated by single-photon and two-photon absorption processes evolves, 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 optical properties of the medium. The reconstruction is usually done in two steps. In the first step, we reconstruct the initial pressure field H in (2) from measured data. This step is the same 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 step of 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. the internal 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 is generally believed that events of two-photon absorption in biological tissues can only happen when the local photon density is sufficiently high (so that the total absorption µ|u| is large enough). In fact, the main difficulty in the development of TP-PAT is to be able to measure the ultrasound signal accurate enough such that the photoacoustic signal due to two-photon absorption is not completely buried by noise in the data. In recent years, many experimental research have been conducted where it is shown that the effect of two-photon absorption can be measured accurately; see, for instance, the study on the feasibility of TP-PAT on various liquid samples in [56,57,58] (solutions), [35,58] (suspensions) and [36] (soft matter).
Despite various experimental study of TP-PAT, a thorough mathematical and numerical analysis of the inverse problems in the second step of TP-PAT is largely missing, not to mention efficient reconstruction algorithms. The objective of this study is therefore to pursue in these directions. In the rest of the paper, we first recall in Section 2 some fundamental mathematical 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. In Section 4 we analyze the linearized problem of simultaneously reconstructing the absorption coefficients and the diffusion coefficient. Numerical simulations are provided in Section 5 to validate the mathematical analysis and demonstrate the quality of the reconstructions. Concluding remarks are offered in Section 6.

The semilinear diffusion model
To prepare for the study of the inverse coefficient problems, we recall in this section some general results on the semilinear diffusion model (1). Thanks to the absolute value operator in the quadratic term µ|u|u in the equation, we can follow the standard theory of calculus of variation, as well as the theory of generalized solutions to elliptic equations in divergence form, to derive desired properties of the solution to the diffusion equation that we will need in the following sections. The results we collected here are mostly minor modifications/simplifications of classical results in [2,5,23,27]. We refer interested readers to these references, 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 usual exterior cone condition [27]. We assume that all the coefficients involved are bounded in the sense that there exist positive constants θ ∈ R and Θ ∈ R such that Unless stated otherwise, we assume also that where W 1,2 (Ω) denotes the usual Hilbert space of L 2 (Ω) functions whose first weak derivative is also in L 2 (Ω). Note that we used W 1,2 (Ω) instead of H 1 (Ω) 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 the above 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, With our assumptions above, it is clear that L is uniformly elliptic, and f (x, z) is continuously differentiable with respect to z onΩ × R.
Proof. This result is scattered in a few places in [2,5] (for instance [5, Theorem 1.6.6]). We provide a sketch of proof here. For any function w ∈ W, we define the following functional associated with the diffusion equation (1): It is straightforward to verify that I[w] : W → R is strictly convex (thanks again to the absolute value in the third term) and differentiable on W with We also verify that the function L(x, z, p) satisfies the following growth conditions: for all x ∈ Ω, z ∈ R and p ∈ R n . It then follows from standard results in calculus of variations [2,5,23] that there exists a unique u ∈ W satisfies and u is the unique weak solution of (1). By Sobolev embedding, when n = 2, 3, there exists q > n, such that u ∈ L q (Ω). This then implies that f (x, u) ∈ L q/2 (Ω) with the assumption (3). Let us rewrite the diffusion equation (1) as Following standard results in [23,27], we conclude that f ∈ L q/2 (Ω) implies u ∈ C α (Ω) for some 0 < α < 1, where α = α(n, Θ/θ). Moreover, when g ∈ C 0 (∂Ω), u ∈ C 0 (Ω). If we assume further that (γ, σ, µ) and g satisfy (4), then f ∈ W 1,2 thanks to the fact that u ∈ C 0 (Ω). Equation (6) then implies that u ∈ W 3,2 (Ω) ∩ C 0 (Ω) [23,27].
We now recall the following comparison principle for the solutions to the semilinear diffusion equation (1).
If, in addition, Ω satisfies the exterior cone condition or u, v ∈ W 2,2 (Ω), then either u ≡ v or u < v.
. With the assumption that u ∈ C 0 (Ω) and v ∈ C 0 (Ω), we conclude that u t is bounded from above when t ∈ [0, 1]. Therefore, a(x) ≤ Λ < ∞ for some Λ > 0. We also verify that f (u, Let w = u − v, we have, from the assumptions in the proposition, that Since L + a is uniformly elliptic, by the weak maximum principle for weak solutions [27,Theorem 8.1], w ≤ 0 in Ω. This then implies that u ≤ v in Ω. If we assume in addition that u, v ∈ W 2,2 (Ω), we can use the strong maximum principle to conclude that w ≡ 0 if 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 1,2 (Ω) and Ω satisfies the exterior cone condition, 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 the semilinear diffusion equation (1). Proposition 2.3. Let u j be the solution to (1) with boundary condition g j , j = 1, 2. Assume that γ, σ, µ and {g j } 2 j=1 satisfy the assumptions in (3) and (4). Then the following statements hold: Proof. (i) follows from the comparison principle in Proposition 2.2 and the fact that u ≡ 0 is a solution to (1) with homogeneous Dirichlet condition g = 0. (ii) By (i), u j ≥ 0. Therefore f (x, u j ) ≥ 0. Therefore, we can have By the maximum principle, sup Ω u j ≤ sup ∂Ω g j . (iii) is a direct consequence of part (ii) of Proposition 2.2.
In the study of the inverse problems in the next sections, we sometimes need the solution to the semilinear diffusion equation to be bounded away from 0. We now prove the following result.
Proof. We follow the arguments in [1]. We again rewrite the PDE as Then by classical gradient estimates, see for instance [29,Proposition 2.20], we know that there exists K > 0, depending on γ, |∇γ| and Ω, such that Using the fact that g ≥ ε, we conclude from this inequality that there exists a d > 0 such that where Due to the fact that u is nonnegative and bounded from above, we have that 0 < θ ≤ c(x) ≤ Θ(1 + sup ∂Ω |g|). We then have that u solves By the Harnack inequality (see [27,Corollary 8.21]), we have that there exists a constant C, depending on d, γ, c, Ω, and Ω d/2 , such that 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 linearization that we perform in Section 4.
Proposition 2.5. The datum H defined in (2) generated from an illumination g ≥ 0 on ∂Ω, viewed as the map is Fréchet differentiable when the coefficients satisfies (3) and (4). where with Proof. We show here only that u is Fréchet differentiable with respect to γ, σ and µ. The rest of the result follows from the chain rule.
Let (δγ, δσ, δµ) ∈ W 1,2 (Ω)×L ∞ (Ω)×L ∞ (Ω) be such that (γ , σ , µ ) = (γ +δγ, σ +δσ, µ+ δµ) satisfies the bounds in (3). Let u be the solution to (1) with coefficients (γ , σ , µ ), and define u = u − u. We then verify that u solves the following linear diffusion equation where we have used the fact that u ≥ 0 and u ≥ 0 following Proposition 2.3 (since g ≥ 0 on ∂Ω). Note also that both u and u are bounded from above by Proposition 2.3. Therefore, σ +µ(u+u ) is bounded from above. Therefore, we have the following standard estimate [27] and v 3 solutions to (9). Then we verify that u satisfies the equation Therefore, we have the following standard estimate We can thus combine (10) with (11) to obtain the bound 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 maximum principles for solutions u and u. The same thing can not be done with respect to γ since we can not control the term ∇ · δγ∇u L ∞ (Ω) with δγ L ∞ (Ω) without much more restrictive assumptions on δγ.

Reconstructing absorption coefficients
We now study inverse problems related to the semilinear diffusion model (1). We first consider the case of reconstructing the absorption coefficients, assuming that the Grüneisen coefficient Γ and the diffusion coefficient γ are both known.

One coefficient with single datum
We now show that with one datum set, we can uniquely recover one of the two absorption coefficients.
Proof. The proof is straightforward. Let w = u − u. We check that w solves Therefore H = H implies w = 0 which is simply u = u. This in turn implies that that is σ + µ|u| = σ + µ| u|. Note that the condition g ≥ ε > 0 implies that u, u ≥ ε > 0 for some ε following Theorem 2.4. This makes it safe to take the ratios H/u and H/ u, and to omit the absolute values on u and u.
To derive the stability estimate, we first observe that Using the fact that u and u are both bounded away from zero, and the triangle inequality, we have, for some constants c 1 and c 2 , On the other hand, classical theory of elliptic equations allows us to derive, from (13), the following bound, for some constant c 3 , The bound in (12) then follows by combining (14) and (15).
The above proof provides an explicit algorithm to reconstruct one of σ and µ from one datum. Here is the procedure. We first solve for u since Γ and γ are known. We then reconstruct σ as if µ is known, or reconstruct µ as if σ is known.
The stability estimate (12) can be made more explicit when one of the coefficients involved is known. For instance, if µ is known, then we have This leads to, using the triangle inequality again, Combining this bound with (15), we have for some constant C . In the same manner, we can derive for the reconstruction of µ if σ is known in advance.

Two coefficients with two data sets
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 1 and g 2 , then we can reconstruct σ + µ|u 1 | and σ + µ|u 2 | where u 1 and u 2 are the solutions to the diffusion equation (1) corresponding to g 1 and g 2 respectively. If we can choose g 1 and g 2 such that |u 2 | − |u 1 | = 0 almost everywhere, we can uniquely reconstruct the pair (σ, µ). This is the idea we have in the following result.
Proposition 3.2. Let Γ and γ be given. Let (H 1 , H 2 ) and ( H 1 , H 2 ) be the data sets corresponding to the coefficients (σ, µ) and ( σ, µ) respectively that are generated with the pair of sources (g 1 , g 2 ). Assume that g i ≥ ε > 0, i = 1, 2, and g 1 − g 2 ≥ ε > 0 for some ε and ε . Then for some constant C. Therefore Collecting the results for both data sets, we have When g 1 and g 2 satisfy the requirements stated in the proposition, we have u 1 − u 2 ≥ ε > 0 for some ε . Therefore, the matrix 1 |u 1 | 1 |u 2 | is invertible. We can then remove this matrix in (23) to show that (σ, µ) = ( σ, µ).
To get the stability estimate in (21), we first verify that This leads to, Therefore, we have It then follows that Meanwhile, we have, from (22), The bound in (21) then follows from (24) and (25).

Reconstructing absorption and diffusion coefficients
We now study inverse problems where we intend to reconstruct more than the absorption coefficients. We start with a non-uniqueness result on the simultaneous reconstructions of all four coefficients Γ, γ, σ, and µ.
4.1 Non-uniqueness in reconstructing (Γ, γ, σ, µ) Let us assume for the moment that γ 1/2 ∈ C 2 (Ω). We introduce the following Liouville We then verify that the semilinear diffusion equation (1) is transformed into the following equation under the Liouville transform: and the datum H is transformed into Let us now define the following functionals: The following result says that once (α, β, ζ 1 ) or (α, β, ζ 2 ) is known, introducing new data would not bring in new information.
Proof. Let us first rewrite the datum as H = ζ 1 v + ζ 2 v 2 . When α and β are known, we know the solution v of (27) for any given g. If ζ 1 is also known, we know also ζ 1 v. We therefore can form the ratio We then find H as H This The proof is complete.
The above theorem says that we can at most reconstruct the triplet (α, β, ζ 1 ) or the triplet (α, β, ζ 2 ). 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 terms of uniqueness of reconstructions.
Similar non-uniqueness results were proved in the case of the regular PAT [8,9]. In that case, it was also shown that if the Grüneisen coefficient Γ is known, for instance from multispectral measurements [10,39], one can uniquely reconstruct the absorption coefficient and the diffusion coefficient simultaneously. In the rest of this section, we consider this case, that is, Γ is known, for our TP-PAT model.

Linearized reconstruction of (γ, σ, µ)
We study the problem of reconstructing (γ, σ, µ), assuming Γ is known, in linearized setting following the general theory of overdetermined elliptic systems developed in [21,47]. For the sake of the readability of the presentation below, we collect some necessary terminologies in the theory of overdetermined elliptic systems in Appendix A. We refer interested readers to [7,34,54] for overviews of the theory in the context of hybrid inverse problems and references therein for more technical details on the theory. Our presentation below follows mainly [7].
The system is formally overdeterminted when J > 3.
To supplement the above system with appropriate boundary conditions, we first observe that the boundary conditions for the solutions {δu j } J j=1 are given already. They are homogeneous Dirichlet conditions since g does not change when the coefficients change. The boundary conditions for (δγ, δσ, δµ) are what need to be determined. In the case of singlephoton PAT, it has been shown that one needs to have γ |∂Ω known to have uniqueness in the reconstruction [8,44]. This is also expected in our case. We therefore take δγ |∂Ω = φ 1 for some known φ 1 . The boundary conditions for σ and µ are given by the data. In fact, on the boundary, u = g. Therefore, we have, from (33) which holds on ∂Ω, that g j δσ + |g j |g j δµ = δH j /Γ, on ∂Ω.
It is straightforward to check that if we take the associated Douglis-Nirenberg numbers as the principal part of A is simply A itself with the −iξ · ∇γ and −∆u j (1 ≤ j ≤ J) terms removed.
In three-dimensional case, we can establish the following result.
Proof. Let us first rewrite the principal symbol A 0 as It is then straightforward to check that A 0 (x, iξ) is of full-rank as long as the following sub-matrix is of full-rank: To simplify the calculation, we introduce Σ j = σ + 2µ|u j |, F j = V j · ξ. We also eliminate the non-zero common factor i γ|ξ| 2 from the first column and u j from each row. Without loss of generality, we check the determinant of first 3 (since J ≥ n + 1 = 4) rows of the simplified version of the submatrixÃ 0 (x, iξ). This determinant is given as With the assumptions on the background coefficients, we can take u j to be the complex geometric optics solution constructed following Theorem 6.6 (in Appendix B) for ρ j with boundary condition g j . Then we have This gives us, Let us define α k = Σ k (|u i | − |u j |) with (k, i, j) ∈ {(1, 2, 3), (2, 3, 1), (3, 1, 2)}. Then we have α 1 + α 2 + α 3 = 0. Let (e 1 , e 2 , e 3 ) an orthonormal basis for R 3 . Then ξ = 3 k=1 c k e k with |c 1 | 2 + |c 2 | 2 + |c 3 | 2 = 1. We take [ This is because if det(A 0 ) = 0, we have R = 0, I = 0, and a 1 + a 2 + a 3 = 0. That is This contradicts the construction of {α k } 3 k=1 ]. Let us now take Then the submatrix formed by u 1 , u 2 and u 4 will have full rank when c 1 = c 2 = c 3 . Therefore the submatrix formed by u 1 , u 2 , u 3 and u 4 is of full rank for any ξ.
To show that (A, B) satisfies the Lopatinskii criterion given in Definition 6.4 for a set of well chosen u j , we first observe that since B = I, we have, from the definition in (55), with the selection of the Douglis-Nirenberg numbers {s i } 2J i=1 and {t j } J+3 j=1 in (37), and the principal part of B has components B 0,11 = 1 and B 0,k = 0 otherwise. Therefore, the system of differential equations in (56) and (57) takes the following form where γ, σ, µ, u j and V j (1 ≤ j ≤ J, are all evaluated at y ∈ ∂Ω. We first deduce from (41) that, for 1 ≤ j ≤ J, Plugging this into (40), we obtain that, for 1 ≤ j ≤ J, Without loss of generality, we consider the system formed by u 1 , u 2 , u 3 . Let F j = V j · ν, We look for eigenvalues of the system as the root of det    We observe first that the above equation admits two repeated roots λ 2,3 = ±|ζ|. Besides that, we have another root Moreover, the eigenvectors corresponding to λ 2,3 are of the form x y   with x and y arbitrary. Therefore, δγ(z) = ce i|λ 1 |z . Using the boundary condition δγ(0) = 0 and the decay condition δγ(z) → 0 as z → ∞, we conclude that δγ(z) ≡ 0. This in turn implies, from (41), that δσ(z) ≡ 0 and δµ(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 > n + 1 2 , We have the following uniqueness and stability result.
Even though it is not necessary here, we can select ρ k such that ρ k < 0 and | ρ k | is sufficiently large to simplify the above equation further to Ω δγ(y)v(x; y) · ∇G(x; y)dy = 0. (45) with v the vector given by We now need the following lemma. Proof. Let u be the solution to Then Therefore (45) implies that u ≡ 0. Therefore − ∇ · δγv = 0, in Ω, δγ = 0, on ∂Ω.
This is a transport equation for δγ that admits the unique solution δγ ≡ 0 with the vector field v satisfying the assumed requirements; see for instance [8,13,15,20,30] and references therein.
It is straightforward to check that we can select {ρ j } 3 j=1 such that the vector field v defined in (46) satisfies the requirement in Lemma 4.4. We then conclude that δγ ≡ 0.
The conditions assumed on the background coefficients ensure the ellipticity of the system as proven in Theorem 4.2. The stability result then follows from (58); see more discussions in [7] and references therein. Note that the simplification inthe last term in (43) is due to the fact that φ k = φ k = 0 when 4 ≤ k ≤ J + 3. Note also that the following simplification can be made in (43): The proof is complete.

Numerical simulations
We present in this section some preliminary numerical reconstruction results using synthetic internal data. We restrict ourselves to two-dimensional settings only to simplify the computation. The spatial domain of the reconstruction is the square Ω = (1, 1) 2 . All the equations in Ω are discretized with a first-order finite element method on triangular meshes. In all the simulations in this section, reconstructions are performed on a finite element mesh with about 6000 triangular elements. The nonlinear system resulted from the discretization of the 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. We performed reconstructions using both noiseless and noisy synthetic data. For the noisy data, we added random noise to the data by simply multiplying each datum by (1 + √ 3 × 10 −2 random) with random a uniformly distributed random variable taking values in [1,1], 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 present reconstruction results from two different numerical methods.
Direct Algorithm. The first method we use is motivated from the method of proofs of Propositions 3.1 and 3.2. When we have J ≥ 2 data sets {H j } J j=1 from J illuminations {g j } J j=1 , we first reconstruct, for each j, u * j as the solutions to in Ω, u * j = g j on ∂Ω.
We then reconstruct σ + µ|u * j | = H j Γu * j . Collecting this quantity from all data, we have, for each point x ∈ Ω,    Newton method that we developed in [43] to solve this minimization problem. It is straightforward to check, following Proposition 2.5, that the gradient of the functional Φ(σ, µ) with respect to σ and µ are given respectively by where z j = Γ(σu j + µ|u j |u j ) − H * j and v j solves Therefore, in each iteration of the optimization algorithm, we need to solve J semilinear diffusion equations for {u j } J j=1 and then J adjoint linear elliptic equations for {v j } J j=1 to evaluate the gradients of the objective function with respect to the unknowns.  Experiment I. We start with a set of numerical experiments on the reconstruction of the two-photon absorption coefficient µ assuming that the single-photon absorption coefficient σ is known. We use data collected from four different sources {g j } 4 j=1 , {H j } 4 j=1 . We perform reconstructions using the Direct Algorithm. In Fig. 2 we show the reconstruction results from noisy synthetic data with noise levels = 0, = 1, = 2, and = 5. The true coefficients used to generate the data are shown in Fig. 1.
To measure the quality of the reconstruction, we use the relative L 2 error. This error is defined as the ratio between (i) the L 2 norm of the difference between the reconstructed coefficient and the true coefficient and (ii) the L 2 norm of the true coefficient, expressed in percentage. The relative L 2 errors in the reconstructions of µ in Fig 2 are 0.00%, 2.45%,4.98%, and 12.23% for = 0, = 1, = 2 and = 5 respectively. Experiment II. One of the main limitations on the Direct Algorithm is that it requires the use of illumination sources that are positive everywhere on the boundary. This is difficult to implement in practical applications. The Least-Square Algorithm, however, does not have such requirement on the optical sources (but it is computationally more expensive). Here we repeat the simulations in Experiment I with the Least-Square Algorithm. The reconstruction results are shown in Fig 3. We observe that, with the same (not exactly the same since the realizations of the noise are different) data sets, the reconstructions from the two different algorithms are of very similar quality. The relative L 2 errors for the reconstructions in Fig 3 are 0.00%, 2.44%,4.62%, and 9.36% respectively for the four cases.
Experiment IV. We now repeat the simulations in Experiment III with the Least-Square Algorithm. The results are shown in Fig. 5. The relative L 2 errors in the reconstructions are now (0.22%, 2.38%), (1.21%, 6.43%),(2.34%, 10.98%), (5.64%, 22.06%) respectively for data with noise levels = 0, = 1, = 2, and = 5. The quality of the reconstructions is slightly lower than, but still comparable to, that in the reconstructions in Experiment III.   Fig. 4 except that the reconstructions are performed using the Least-Square Algorithm.
We observe from the above simulation results that, in general, the quality of the reconstructions is very high. When we have the illumination sources that satisfy the positivity requirement on the whole boundary of the domain, the Direct Algorithm provides an efficient and robust reconstruction method. The Least-Square Algorithm is less efficient but is 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 a couple of trial-error test. It is by no means the optimal regularization parameter that can be selected through more sophisticated algorithms [22]. However, this is an issue that we think is not important at the current stage of this project. Therefore, we did not pursue further in this direction.
We also observe that the performances of the reconstructions of σ and µ are different in both the direct and the least-square methods: the reconstructions of σ is visually better than those 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 fully understand this difference in numerical stability.

Conclusion and further remarks
We studied in this paper inverse problems in quantitative photoacoustic tomography with two-photon absorption. We derived uniqueness and stability results in the reconstruction of single-photon and two-photon absorption coefficients, and proposed explicit reconstruction methods in this case with well-selected illumination sources. We also studied the inverse problem of reconstructing the diffusion coefficient in addition to the absorption coefficients and obtained partial results on the uniqueness and stability of the reconstructions for the linearized problem. We presented some numerical studies based on the explicit reconstruction procedures as well as numerical optimization techniques to demonstrate the type of quality that can be achieved in reasonably controlled environments (where noise strength in the 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 be nice to generalize the uniqueness and stability results in Section 4, on multiple coefficient reconstructions in linearized settings, to the fully nonlinear problem. Computationally, detailed numerical analysis, in three-dimensional setting, need to be performed to quantify the errors in the reconstructions in practically relevant scenarios. It is especially important to perform reconstructions starting from acoustic data directly, following for instance the one-step reconstruction strategy developed in [19], to see how sensitive the reconstruction of the two-photon absorption coefficient is with respect to noise in the acoustic data. On the modeling side, it is very interesting to see if the current study can be generalized to radiative transport type models for photon propagation.
The functionũ * is called a complex geometric optics solution to (60) and 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)).
It then follows from the Newton-Kantorovich theorem [41] that, when |ρ| is sufficiently large, there exists a solution to (1) in the ball of radius r =ce −κ κ|ρ| centered at u * , in W 3,2 (Ω). The solution is of the form (66).
Let us remark that the assumption (65) on ρ and Ω is feasible. Let (e 1 , e 2 , e 3 ) be a basis of R 3 . Since Ω is bounded, we can place Ω in a box in the first octant of the (e 1 , e 2 , e 3 ) system. Let us take (ρ) = a 1 e 1 + a 2 e 2 + a 3 e 3 for some {a k } 3 k=1 . We check that we can haveκ| (ρ)| ≤ (ρ) ≤ κ| (ρ)| for some constantκ and κ. If we further take | (ρ)| to be proportional to | (ρ)|, we can have (65); see Section 4 for the choices of different ρ in our analysis.