# Quantitative PAT with simplified $P_{N}$ approximation

Hongkai Zhao  Department of Mathematics, Duke University, NC 27705,
<EMAIL_ADDRESS>.  Yimin Zhong  Department of Mathematics, Duke University,
NC 27705<EMAIL_ADDRESS>.

###### Abstract

The photoacoustic tomography (PAT) is a hybrid modality that combines the
optics and acoustics to obtain high resolution and high contrast imaging of
heterogeneous media. In this work, our objective is to study the inverse
problem in the quantitative step of PAT which aims to reconstruct the optical
coefficients of the governing radiative transport equation from the ultrasound
measurements. In our analysis, we take the simplified $P_{N}$ approximation of
the radiative transport equation as the physical model and then show the
uniqueness and stability for this modified inverse problem. Numerical
simulations based on synthetic data are presented to validate our analysis.

Key words. photoacoustic tomography (PAT), radiative transport equation,
simplified $P_{N}$ method, diffusion approximation, numerical reconstruction

##  1 Introduction

The photoacoustic tomography (PAT) [3, 6, 9, 19, 30, 26] is a emerging hybrid
imaging modality that reconstructs high resolution images of optical
properties of heterogeneous media. The PAT experiment uses a pulse of near-
infra-red (NIR) laser into the medium of interest (e.g. fat, bone, tumor
tissues). These photons propagate inside the medium by following the radiative
transport process. During the propagation, a portion of the photons is
absorbed by the medium and then converted into heat which causes a local
thermoelastic expansion. Such expansion induces a transient pressure change
and leads to the propagation of ultrasound. The ultrasound signals are
measured around the boundary of the medium and we need to infer the optical
properties from the acoustic measurements.

The photon transport process is usually modeled by radiative transport
equation. Let $X=\Omega\times\mathbb{S}^{2}$, where $\Omega$ is the physical
domain and $\mathbb{S}^{2}$ denotes the unit sphere in 3D, the photon density
function $u(\mathbf{x},\mathbf{v})$ satisfies the following

$\displaystyle\mathbf{v}\cdot\nabla
u(\mathbf{x},\mathbf{v})+\sigma_{t}(\mathbf{x})u(\mathbf{x},\mathbf{v})$
$\displaystyle=\sigma_{s}(\mathbf{x})\int_{\mathbb{S}^{2}}p(\mathbf{v}\cdot\mathbf{v}^{\prime})u(\mathbf{x},\mathbf{v}^{\prime})d\mathbf{v}^{\prime}\quad$
in $\displaystyle X\,,$ (1) $\displaystyle u(\mathbf{x},\mathbf{v})$
$\displaystyle=f(\mathbf{x},\mathbf{v})\quad$ on $\displaystyle\Gamma_{-}\,,$

where
$\Gamma_{-}=\\{(\mathbf{x},\mathbf{v})\in\partial\Omega\times\mathbb{S}^{d-1}\mid-{\boldsymbol{\nu}}(\mathbf{x})\cdot\mathbf{v}>0\\}$
is the incoming boundary set. $\sigma_{s},\sigma_{t}$ are the scattering and
total absorption coefficients respectively,
$\sigma_{a}:=\sigma_{t}-\sigma_{s}$ is the intrinsic absorption coefficient.
$f(\mathbf{x},\mathbf{v})$ is the external illumination source. The scattering
phase function $p(\mathbf{v}\cdot\mathbf{v}^{\prime})$ is usually chosen as
the Henyey-Greenstein function

$p(\cos\theta)=\frac{1}{4\pi}\frac{1-g^{2}}{(1+g^{2}-2g\cos\theta)^{3/2}}\,,$
(2)

where $g\in(-1,1)$ is the anisotropy parameter.

The energy absorbed by the medium is
$\sigma_{a}\int_{\mathbb{S}^{2}}u(\mathbf{x},\mathbf{v})d\mathbf{v}$, then the
initial pressure field generated by the _photoacoustic effect_ is:

$H(\mathbf{x}):=\Upsilon(\mathbf{x})\sigma_{a}(\mathbf{x})\int_{\mathbb{S}^{2}}u(\mathbf{x},\mathbf{v})d\mathbf{v},$
(3)

where $\Upsilon(\mathbf{x})$ is the dimensionless Grüneisen coefficient which
measures the efficiency of the photoacoustic effect.

Then the initial pressure field $H(\mathbf{x})$ propagates the ultrasound
wave, which satisfies the following equation [27]:

$\displaystyle\frac{1}{c^{2}(\mathbf{x})}\frac{\partial^{2}p(\mathbf{x},t)}{\partial
t^{2}}-\Delta_{\mathbf{x}}p(\mathbf{x},t)$ $\displaystyle=0,\quad$ in
$\displaystyle\mathbb{R}^{3}\times[0,\infty),$ (4) $\displaystyle
p(\mathbf{x},0)=H(\mathbf{x}),\quad\frac{\partial p}{\partial
t}(\mathbf{x},0)$ $\displaystyle=0,\quad$ in $\displaystyle\mathbb{R}^{3}.$

Here $c(\mathbf{x})$ is the wave speed of the underlying medium. The measured
acoustic signals are $p(\mathbf{x},t)$ on $\partial\Omega\times[0,T]$ for
sufficient large observation time $T$.

The usual reconstruction of PAT is a two-step process. The first step is to
reconstruct the initial pressure field $H(\mathbf{x})$ from the ultrasound
measurements. This problem has been studied extensively by [1, 2, 27, 15, 16]
and the references therein. Here we assume this step has been finished and
recovered the initial pressure field $H(\mathbf{x})$ and we focus on the
second step to reconstruct the optical properties
$(\sigma_{a},\sigma_{s},\Upsilon)$ from the quantity $H(\mathbf{x})$. Under
the diffusion approximation, this quantitative PAT (qPAT) problem has been
well studied [10, 25, 4, 5]. However, with the radiative transport equation
(1), the multi-source inverse problem theory has not been well established
except for _albedo_ type data [20, 4], which requires infinitely many
angularly resolved illumination sources $f(\mathbf{x},\mathbf{v})$. It is
still unclear about the uniqueness and stability of the reconstructions with
finite many source functions or angularly independent sources.

In this paper, we aim to study the qPAT problem with the simplified $P_{N}$
approximation to the equation (1) with angularly independent source functions,
that is, $f(\mathbf{x},\mathbf{v})=f(\mathbf{x})$. The simplified $P_{N}$
approximation is also referred as $SP_{N}$ method, which is utilized to solve
the radiative transport equation by forming a system of elliptic equations
[22, 18]. The $SP_{N}$ approximation with relatively small $N\leq 7$ has been
applied to many optical imaging methods [31, 8, 17] and outperforms the
traditional diffusion approximation ($P_{1}$ method). Theoretically, the
simplified $P_{N}$ approximation is derived from the $P_{N}$ formulation [13]
and the $P_{N}$ approximation converges to the exact solution of RTE as
$N\to\infty$. Under appropriate conditions $SP_{N}$ and $P_{N}$ are
equivalent, see [22], however in general, they are different and not
necessarily converging to the same limit. For the qPAT problem, the case with
$N=3$ has been considered in [11], in our work, we extend the theory to
arbitrary order $N$ under a unified framework.

Under the $SP_{N}$ approximation, the RTE’s solution is expanded with Legendre
polynomials, which derives a weakly coupled diffusion equation system [18].
Formally, the 3D $SP_{N}$ approximation takes the 1D $P_{N}$ equations and
replace the diffusion operators with 3D’s counterpart, which is

$\displaystyle-\left(\frac{n+1}{2n+1}\right)\nabla\frac{1}{\sigma_{n+1}}\nabla\left[\left(\frac{n+2}{2n+3}\right)\phi_{n+2}+\left(\frac{n+1}{2n+3}\right)\phi_{n}\right]$
(5)
$\displaystyle\quad-\left(\frac{n}{2n+1}\right)\nabla\frac{1}{\sigma_{n-1}}\nabla\left[\left(\frac{n}{2n-1}\right)\phi_{n}+\left(\frac{n-1}{2n-1}\right)\phi_{n-2}\right]+\sigma_{n}\phi_{n}=0$

for $n=2,4,\dots,N-1$, where coefficients
$\sigma_{n}=\sigma_{a}+\sigma_{s}(1-g^{n})$, $\sigma_{0}=\sigma_{a}$.
Physically speaking, $\phi_{n}$ represents the $n$-th Legendre moments of the
solution $u(\mathbf{x},\mathbf{v})$ and $\phi_{0}$ will be the angular average
of the solution. The corresponding mixed boundary conditions are derived from
the 1D $P_{N}$ equation’s boundary conditions by replacing
$\frac{\rm{d}}{\rm{d}x}$ with $\mathbf{n}\cdot\nabla$ on the boundary [18].

In the following context, we let

$\displaystyle\varphi_{n}=(2n-1)\phi_{2n-2}+(2n)\phi_{2n},\quad
n=1,2,\dots,(N+1)/2\,,$ (6)

and the column vector
$\Phi=[\varphi_{1},\varphi_{2},\dots,\varphi_{(N+1)/2}]^{\text{T}}$, which
satisfies

$\Phi=\begin{pmatrix}\varphi_{1}\\\ \varphi_{2}\\\ \varphi_{3}\\\ \vdots\\\
\varphi_{(N+1)/2}\end{pmatrix}=M\begin{pmatrix}\phi_{0}\\\ \phi_{2}\\\
\phi_{4}\\\ \vdots\\\ \phi_{N-1}\end{pmatrix}\text{ where
}M:=\begin{pmatrix}1&2&0&\dots&0\\\ 0&3&4&\dots&0\\\ 0&0&5&\dots&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\ 0&0&0&\dots&N\end{pmatrix}\,,$ (7)

since the matrix $M$ is upper triangular, its inverse is also upper
triangular, let the row vector $\mathbf{s}_{k}=[s_{k,1},\dots,s_{k,(N+1)/2}]$
represent the $k$-th row of the inverse matrix $M^{-1}$, according to Lemma
A.1, its entries are

$s_{k,l}=\begin{cases}\frac{1}{2k-1}(-1)^{l-k}\frac{((2l-2))!!}{(2l-1)!!}\frac{(2k-1)!!}{(2k-2)!!},\quad&l\geq
k\,,\\\ 0,&\text{otherwise}\,.\end{cases}$ (8)

Then we derive the $SP_{N}$ diffusion system

$-\nabla\cdot\frac{A_{n}}{(4n-1)\sigma_{2n-1}}\nabla\varphi_{n}-\nabla\cdot\frac{B_{n}}{(4n-5)\sigma_{2n-3}}\nabla\varphi_{n-1}+\sigma_{2n-2}\mathbf{s}_{n}\cdot\Phi=0$
(9)

for $n=1,2,\dots,(N+1)/2$, where the constants $A_{n},B_{n}$ are defined by

$\displaystyle A_{n}$ $\displaystyle=\frac{2n-1}{(4n-3)},\quad
B_{n}=\frac{2n-2}{(4n-3)}.$

We remark that the following matrix differs from $M^{\text{T}}$ by only a
factor of diagonal matrix,

$\begin{pmatrix}A_{1}&0&0&\dots&0\\\ B_{2}&A_{2}&0&\dots&0\\\
0&B_{3}&A_{3}&\dots&0\\\ \vdots&\vdots&\ddots&\ddots&\vdots\\\
0&0&0&B_{(N+1)/2}&A_{(N+1)/2}\end{pmatrix}=LM^{\text{T}}\;\text{ where
}\;L:=\begin{pmatrix}1&0&0&\dots&0\\\ 0&\frac{1}{5}&0&\dots&0\\\
0&0&\frac{1}{9}&\dots&0\\\ \vdots&\vdots&\ddots&\ddots&\vdots\\\
0&0&0&0&\frac{1}{2N-1}\end{pmatrix}\,.$ (10)

The corresponding mixed boundary conditions are

$\sum_{n=1}^{(N+1)/2}\mu_{2m-1,2n-2}\phi_{2n-2}+\frac{1}{(4m-1)\sigma_{2m-1}}\frac{\partial\varphi_{m}}{\partial\mathbf{n}}=k_{2m-1}f,\quad
m=1,2,\dots,(N+1)/2\,.$ (11)

The constants $\mu_{2m-1,2n-2}$ and $k_{2m-1}$ are

$\displaystyle\mu_{2m-1,2n-2}$
$\displaystyle=(4n-3)\int_{0}^{1}P_{2m-1}(x)P_{2n-2}(x)dx$
$\displaystyle=(-1)^{m+n-1}\frac{\Gamma(m+\frac{1}{2})\Gamma(n-\frac{1}{2})}{\pi\Gamma(m)\Gamma(n)(m+n-1)}\frac{(4n-3)}{2n-2m-1}\,,$
$\displaystyle k_{2m-1}$
$\displaystyle=(-1)^{m-1}\frac{(2m-1)!!}{(2m-1)(2m)(2m-2)!!}\,.$

For convenience, we also define the following matrices $\Sigma_{e}$ and $R$
for later uses,

$\displaystyle\Sigma_{e}$
$\displaystyle=\begin{pmatrix}\sigma_{0}&0&0&\dots&0\\\
0&\sigma_{2}&0&\dots&0\\\ 0&0&\sigma_{4}&\dots&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\ 0&0&0&\dots&\sigma_{N-1}\end{pmatrix}$

and

$R=(R_{ij})_{i,j=1,2,\dots,(N+1)/2},\quad\text{with }R_{ij}=\mu_{2i-1,2j-2}.$

In the quantitative photoacoustic tomography, we suppose the datum
$H(\mathbf{x})=\Upsilon(\mathbf{x})\sigma_{a}(\mathbf{x})\phi_{0}(\mathbf{x})$
has been reconstructed from the measured acoustic signals. In the following
sections, we will analyze the uniqueness and stability of reconstruction of
the coefficients $\Upsilon,\sigma_{a},\sigma_{s}$ from the internal data $H$.
We also make the following general assumptions for the rest of paper.

  1. $\mathcal{A}$-i

The physical domain $\Omega$ is simply connected with $C^{2,1}$ boundary.

  2. $\mathcal{A}$-ii

The coefficients $(\sigma_{a},\sigma_{s},\Upsilon)$ are non-negative and
bounded. There exists constants $\underline{c}$ and $\overline{c}$ that

$0\leq\underline{c}\leq\sigma_{a},\sigma_{s},\Upsilon\leq\overline{c}<\infty\,.$
(12)

  3. $\mathcal{A}$-iii

The coefficient $\sigma_{a},\sigma_{s}\in C^{1,1}(\overline{\Omega})$. There
exists a constant $M$ that

$\|\sigma_{a}\|_{C^{1,1}(\Omega)},\|\sigma_{s}\|_{C^{1,1}(\Omega)}\leq
M<\infty\,.$

Moreover, $\sigma_{a},\sigma_{s}$ are both known on $\partial\Omega$.

  4. $\mathcal{A}$-iv

The boundary source function $f\in H^{5/2}(\partial\Omega).$

The rest of the paper is organized as follows. We first present in Section 2
some general properties of the forward problem with $SP_{N}$ approximation.
Then in Section 3 we consider the reconstruction of a single coefficient from
a single data set $H(\mathbf{x})$ and the reconstruction of two coefficients
simultaneously with multiple data in linearized settings. We then demonstrate
some numerical simulations based on synthetic data in Section 4 to validate
some of our theoretical results. Conclusions are found in Section 5.

##  2 General properties

For the forward problem, we establish the wellposedness for the $SP_{N}$
approximation. In order to show there exists a unique weak solution
$\Phi\in[H^{1}(\Omega)]^{(N+1)/2}$ for (9) and (11), we only have to consider
the corresponding variational form for the diffusion system. Let
$\Psi=[\psi_{1},\psi_{2},\dots,\psi_{(N+1)/2}]^{T}\in[H^{1}(\Omega)]^{(N+1)/2}$
be a test function, then the weak form of the $SP_{N}$ system is

$\displaystyle B(\Phi,\Psi)$
$\displaystyle:=\sum_{n=1}^{(N+1)/2}\int_{\Omega}\frac{1}{(4n-1)\sigma_{2n-1}}\nabla\varphi_{n}\cdot\nabla\psi_{n}d\mathbf{x}$
(13) $\displaystyle\quad+\int_{\Omega}\langle
M^{-\text{T}}L^{-1}\Sigma_{e}M^{-1}\Phi,\Psi\rangle
d\mathbf{x}+\int_{\partial\Omega}\langle RM^{-1}\Phi,\Psi\rangle d\mathbf{s}$
$\displaystyle=\mathcal{L}(f,\Psi)\,,$

where $B(\cdot,\cdot)$ is a bilinear form, $\mathcal{L}$ is a linear
functional only involving boundary integrals that

$\mathcal{L}(f,\Psi):=\sum_{n=1}^{(N+1)/2}k_{2n-1}\int_{\partial\Omega}f\psi_{n}d\mathbf{s}\,.$
(14)

We prove the following property of the bilinear form $B(\cdot,\cdot)$.

######  Theorem 2.1.

The bilinear form (13) is bounded and strictly coercive for any $SP_{N}$
approximation.

###### Proof.

The boundedness is obvious since $L$, $M$ are both invertible matrices. We
only need to prove the coerciveness. In the following, we will show that the
matrices $M\Sigma_{e}^{-1}LM^{\text{T}}$ and $M^{T}R$ are positive definite.
For $M\Sigma_{e}^{-1}LM^{\text{T}}$, it is obvious since the diagonal matrix
$\Sigma_{e}^{-1}L$ has all positive entries. For the matrix $M^{\text{T}}R$,
we compute its $(i,k)$-th entry by

$\displaystyle(2i-1)R_{i,k}+(2i-2)R_{i-1,k}$ (15) $\displaystyle=$
$\displaystyle(2i-1)(4k-3)\int_{0}^{1}P_{2i-1}(x)P_{2k-2}(x)dx+(2i-2)(4k-3)\int_{0}^{1}P_{2i-3}(x)P_{2k-2}(x)dx$
$\displaystyle=$
$\displaystyle(4i-3)(4k-3)\int_{0}^{1}xP_{2i-2}(x)P_{2k-2}(x)dx$
$\displaystyle=$ $\displaystyle\int_{0}^{1}q_{i}(x)q_{k}(x)dx\quad\text{with
}q_{i}(x)=(4i-3)\sqrt{x}P_{2i-2}(x)\,,$

where $P_{k}$ is the $k$-th Legendre polynomial and we have used the
recurrence relation

$(2i-1)P_{2i-1}(x)+(2i-2)P_{2i-3}(x)=(4i-3)xP_{2i-2}(x)\,.$ (16)

Hence $M^{\text{T}}R$ is semi-positive definite. On the other hand, if there
is a vector
$\mathbf{z}=[z_{1},\dots,z_{(N+1)/2}]^{\text{T}}\in\mathbb{R}^{(N+1)/2}$ that

$0=\mathbf{z}^{\text{T}}(M^{\text{T}}R)\mathbf{z}=\int_{0}^{1}x\left(\sum_{k=1}^{(N+1)/2}z_{k}(4k-3)P_{2k-2}(x)\right)^{2}dx\,,$
(17)

then for any $x\in[0,1]$, the following polynomial must vanish,

$\sum_{k=1}^{(N+1)/2}z_{k}(4k-3)P_{2k-2}(x)=0\,.$ (18)

Hence the polynomial equals zero for any $x\in\mathbb{R}$ and use the fact
$\\{P_{k}(x)\\}_{k=1}^{(N+1)/2}$ forms an orthogonal basis on $[-1,1]$, then
$\forall k,z_{k}=0$. Therefore $M^{\text{T}}R$ is strictly positive definite,
so is $RM^{-1}=M^{-\text{T}}(M^{\text{T}}R)M^{-1}$. ∎

The wellposedness immediately derives from the Lax-Milgram theorem, there
exists a unique weak solution $\Phi\in[H^{1}(\Omega)]^{(N+1)/2}$ for arbitrary
odd integer $N$. In fact, using the assumptions $\mathcal{A}$-i to
$\mathcal{A}$-iv, the regularity theorem of elliptic systems [23] implies that
the unique solution $\Phi\in[H^{3}(\Omega)]^{(N+1)/2}$, by the Sobolev
embedding, the solution $\Phi\in[C^{1,1/2}(\Omega)]^{(N+1)/2}$.

##  3 Reconstruction under $SP_{N}$ approximation

Generally speaking, if $\sigma_{a}$ is not negligible, the inverse problem is
highly nonlinear and very challenging. Therefore in the following, we only
consider the practical scenario that $\sigma_{a}\ll\sigma_{s}(1-g)$, which
means we can simplify the coefficients
$\sigma_{n}=\sigma_{a}+(1-g^{n})\sigma_{s}\simeq(1-g^{n})\sigma_{s}$ for
$n\geq 1$ and $\sigma_{0}=\sigma_{a}$. This simplification decouples the
coefficients $\sigma_{a}$ and $\sigma_{s}$. In particular, if $g=0$, there is
no need to perform such simplification.

Reconstruction of $\sigma_{a}$ only. Suppose the coefficients
$\Upsilon,\sigma_{s}$ are known on $\Omega$, we consider the reconstruction of
$\sigma_{a}$ from a single measurement datum $H$. Using the assumption that
$\sigma_{n}\simeq(1-g^{n})\sigma_{s}$ for $n\geq 1$, the coefficients
$\sigma_{n}$ are all known for $n\geq 1$. Since $\sigma_{a}=\sigma_{0}$, then
using
$H(\mathbf{x})=\Upsilon(\mathbf{x})\sigma_{a}(\mathbf{x})\phi_{0}(\mathbf{x})$
and $\phi_{0}=\mathbf{s}_{1}\cdot\Phi$, we derive that

$\sigma_{0}(\mathbf{x})=\frac{H(\mathbf{x})}{\Upsilon(\mathbf{x})\mathbf{s}_{1}\cdot\Phi}\,.$
(19)

By isolating the term relevant to $\sigma_{0}$ ($n=0$) in the bilinear form
(13), we can reformulated it as

$\displaystyle B(\Phi,\Psi)$
$\displaystyle=\sum_{n=1}^{(N+1)/2}\int_{\Omega}\frac{1}{(4n-1)\sigma_{2n-1}}\nabla\varphi_{n}\cdot\nabla\psi_{n}d\mathbf{x}+\int_{\Omega}\frac{H(\mathbf{x})}{\Upsilon(\mathbf{x})\mathbf{s}_{1}\cdot\Phi}(\mathbf{s}_{1}\cdot\Phi)(\mathbf{s}_{1}\cdot\Psi)d\mathbf{x}$
(20)
$\displaystyle\quad+\sum_{n=2}^{(N+1)/2}\int_{\Omega}(4n-3)\sigma_{2n-2}(\mathbf{s}_{n}\cdot\Phi)(\mathbf{s}_{n}\cdot\Psi)d\mathbf{x}$
$\displaystyle\quad+\int_{\partial\Omega}\langle RM^{-1}\Phi,\Psi\rangle
d\mathbf{s}\,.$

where the row vector $\mathbf{s}_{k}$ denotes the $k$-th row of $M^{-1}$. We
can establish the following uniqueness and stability result.

######  Theorem 3.1.

Given any $SP_{N}$ approximation, under the assumptions $\mathcal{A}$-i to
$\mathcal{A}$-iv and suppose $(\Upsilon,\sigma_{s})$ are known, $\sigma_{a,1}$
and $\sigma_{a,2}$ are two admissible absorption coefficients, $H_{1},H_{2}$
are the corresponding internal data, respectively. Then $H_{1}=H_{2}$ implies
$\sigma_{a,1}=\sigma_{a,2}$ and the following stability estimate holds

$\displaystyle\|(\sigma_{a,1}-\sigma_{a,2})\frac{H_{1}}{\sigma_{a,1}\Upsilon}\|_{L^{2}(\Omega)}\leq
C\left\|(H_{1}-H_{2})/\Upsilon\right\|_{L^{2}(\Omega)}\,,$ (21)

where $C=C(N,\Omega)$ is a positive constant depending on $N$ and $\Omega$
only.

###### Proof.

Let $\Phi$ and $\widetilde{\Phi}$ be the weak solutions to the $SP_{N}$ system
for the absorption coefficients $\sigma_{a,1}$ and $\sigma_{a,2}$,
respectively. Let
$\delta\Phi:=\Phi-\widetilde{\Phi}=[\delta\varphi_{1},\dots,\delta\varphi_{(N+1)/2}]^{\text{T}}$,
then from the bilinear form (20), we obtain the equation

$\widetilde{B}(\delta\Phi,\Psi)=-\int_{\Omega}\frac{H_{1}-H_{2}}{\Upsilon}(\mathbf{s}_{1}\cdot\Psi)d\mathbf{x}\,,$
(22)

where the above modified bilinear form $\widetilde{B}(\cdot,\cdot)$ is

$\displaystyle\widetilde{B}(\delta\Phi,\Psi)$
$\displaystyle=\sum_{n=1}^{(N+1)/2}\int_{\Omega}\frac{1}{(4n-1)\sigma_{2n-1}}\nabla\delta\varphi_{n}\cdot\nabla\psi_{n}d\mathbf{x}$
(23)
$\displaystyle\quad+\sum_{n=2}^{(N+1)/2}\int_{\Omega}(4n-3)\sigma_{2n-2}(\mathbf{s}_{n}\cdot\delta\Phi)(\mathbf{s}_{n}\cdot\Psi)d\mathbf{x}$
$\displaystyle\quad+\int_{\partial\Omega}\langle RM^{-1}\delta\Phi,\Psi\rangle
d\mathbf{s}\,.$

Since $RM^{-1}$ is strictly positive definite, the coerciveness of
$\widetilde{B}(\cdot,\cdot)$ is immediately deduced from the Poincaré-Sobolev
inequality [32] that $\forall u\in H^{1}(\Omega)$,
$\Omega\subset\mathbb{R}^{3}$,

$\|u\|_{L^{6}(\Omega)}\leq C_{ps}\left(\|\nabla
u\|_{L^{2}(\Omega)}+\|u\|_{L^{1}(\partial\Omega)}\right)\,,$ (24)

where $C_{ps}=C_{ps}(\Omega)$ is a positive constant depending on $\Omega$
only. Therefore there exists another constant $C_{1}(N,\Omega)$ that

$\displaystyle C_{1}(N,\Omega)\|\delta\Phi\|^{2}_{[L^{2}(\Omega)]^{(N+1)/2}}$
$\displaystyle\leq\widetilde{B}(\delta\Phi,\delta\Phi)=-\int_{\Omega}\frac{H_{1}-H_{2}}{\Upsilon}(\mathbf{s}_{1}\cdot\delta\Phi)d\mathbf{x}$
(25)
$\displaystyle\leq\left\|(H_{1}-H_{2})/\Upsilon\right\|_{L^{2}(\Omega)}\|\mathbf{s}_{1}\cdot\delta\Phi\|_{L^{2}(\Omega)}$
$\displaystyle\leq\left\|(H_{1}-H_{2})/\Upsilon\right\|_{L^{2}(\Omega)}\left\|\mathbf{s}_{1}\right\|_{\ell^{2}}\|\delta\Phi\|_{[L^{2}(\Omega)]^{(N+1)/2}}$

by the Hölder inequality that

$\|\mathbf{v}\cdot\Phi\|_{W^{k,p}(\Omega)}\leq\|\mathbf{v}\|_{\ell^{q}}\|\Phi\|_{[W^{k,p}(\Omega)]^{(N+1)/2}},\quad\frac{1}{p}+\frac{1}{q}=1,$

where the norm $\|\cdot\|_{[W^{k,p}(\Omega)]^{(N+1)/2}}$ is defined by

$\|\mathbf{f}\|_{[W^{k,p}(\Omega)]^{(N+1)/2}}^{p}=\sum_{n=1}^{(N+1)/2}\|f_{n}\|_{W^{k,p}(\Omega)}^{p},\quad\mathbf{f}=[f_{1},\dots,f_{(N+1)/2}]^{\text{T}}\,.$

The estimate (25) implies

$\|\delta\Phi\|_{[L^{2}(\Omega)]^{(N+1)/2}}\leq\frac{1}{C_{1}(N,\Omega)}\left\|(H_{1}-H_{2})/\Upsilon\right\|_{L^{2}(\Omega)}\left\|\mathbf{s}_{1}\right\|_{\ell^{2}}\,.$
(26)

Therefore the uniqueness is proved. For the stability estimate, we compute

$\displaystyle\frac{H_{1}-H_{2}}{\Upsilon}$
$\displaystyle=\sigma_{a,1}(\mathbf{s}_{1}\cdot\Phi)-\sigma_{a,2}(\mathbf{s}_{1}\cdot\widetilde{\Phi})$
(27)
$\displaystyle=(\sigma_{a,1}-\sigma_{a,2})(\mathbf{s}_{1}\cdot\Phi)+\sigma_{a,2}(\mathbf{s}_{1}\cdot\delta\Phi)\,.$

Thus using (26), we obtain

$\displaystyle\|(\sigma_{a,1}-\sigma_{a,2})(\mathbf{s}_{1}\cdot\Phi)\|_{L^{2}(\Omega)}$
$\displaystyle\leq\left\|\frac{H_{1}-H_{2}}{\Upsilon}\right\|_{L^{2}(\Omega)}+\|\sigma_{a,2}(\mathbf{s}_{1}\cdot\delta\Phi)\|_{L^{2}(\Omega)}$
(28)
$\displaystyle\leq\left(1+\frac{\overline{c}\|\mathbf{s}_{1}\|_{\ell^{2}}^{2}}{C_{1}(N,\Omega)}\right)\left\|\frac{H_{1}-H_{2}}{\Upsilon}\right\|_{L^{2}(\Omega)}\,.$

Our proof is completed by noticing (19). ∎

The reconstruction algorithm for $\sigma_{a}$ is then naturally divided into
two steps. First, we solve $\Phi$ from the modified bilinear form (20), then
recover $\sigma_{a}$ by the relation (19) whenever
$\mathbf{s}_{1}\cdot\Phi\neq 0$. For the general $SP_{N}$ system, we cannot
guarantee the positivity of $\phi_{0}=\mathbf{s}_{1}\cdot\Phi$ for any
positive source function $f(\mathbf{x})$.

Under appropriate conditions [22], the $SP_{N}$ approximation will be
eventually converging to the radiative transfer model. However, intuitively,
when the order of system $N$ grows, the reconstruction of the coefficients
will be less stable due to the coupling of the Legendre moments in the
solution. In the following, we study the relation of reconstruction’s
stability and the system order $N$. It can be shown that the reconstruction’s
stability estimate’s constant in Theorem 3.1 grows at most proportional to
$N^{11/8}(1+\log N)$.

######  Corollary 3.2.

Under the assumptions $\mathcal{A}$-i to $\mathcal{A}$-iv, suppose
$(\Upsilon,\sigma_{s})$ are known, $\sigma_{a,1}$, $\sigma_{a,2}$ are two
admissible absorption coefficients, $H_{1},H_{2}$ are the corresponding
internal data, respectively. Then $H_{1}=H_{2}$ implies
$\sigma_{a,1}=\sigma_{a,2}$ and the following stability estimate holds

$\displaystyle\|(\sigma_{a,1}-\sigma_{a,2})\frac{H_{1}}{\sigma_{a,1}\Upsilon}\|_{L^{2}(\Omega)}\leq
CN^{11/8}(1+\log N)\left\|(H_{1}-H_{2})/\Upsilon\right\|_{L^{2}(\Omega)}\,,$
(29)

where $C=C(\Omega)$ is a positive constant independent of $N$.

###### Proof.

Recall the estimate (28), we only have to give an estimate for
$C_{1}(N,\Omega)$ and $\|\mathbf{s}_{1}\|_{\ell^{2}}^{2}$ with respect to $N$.
Using the equation (23), we can estimate the lower bound of the coerciveness
for $\widetilde{B}(\cdot,\cdot)$ by neglecting the second term,

$\widetilde{B}(\delta\Phi,\delta\Phi)\geq\frac{1}{(2N+1)\sigma_{s}}\|\nabla\delta\Phi\|_{[L^{2}(\Omega)]^{(N+1)/2}}^{2}+\lambda_{(N+1)/2}(RM^{-1})\|\delta\Phi\|_{[L^{2}(\partial\Omega)]^{(N+1)/2}}^{2}\,,$
(30)

where $\lambda_{n}(K)$ denotes the $n$-th singular value of $K$ ordered from
largest to smallest. Use the inequality introduced in [21], we estimate the
smallest singular value of $RM^{-1}$ that

$\lambda_{(N+1)/2}(RM^{-1})\geq\lambda_{(N+1)/2}(R)\lambda_{(N+1)/2}(M^{-1})=\lambda_{(N+1)/2}(R)/\lambda_{1}(M)\,.$
(31)

Since the largest singular value $\lambda_{1}(M)=\|M\|_{op}$, let
$\mathbf{v}=[v_{1},v_{2},\dots,v_{(N+1)/2}]^{T}\in\mathbb{R}^{(N+1)/2}$ and
take the convention that $v_{(N+3)/2}=0$, then use the Cauchy-Schwartz
inequality, we obtain

$\|M\|_{op}^{2}=\sup_{\|\mathbf{v}\|=1}\|M\mathbf{v}\|^{2}=\sup_{\|\mathbf{v}\|=1}\sum_{k=1}^{(N+1)/2}((2k-1)v_{k}+(2k)v_{k+1})^{2}\leq
2N^{2}\|\mathbf{v}\|^{2}=2N^{2}\,.$ (32)

Therefore we have $\lambda_{1}(M)\leq\sqrt{2}N$. In the next, we only need to
estimate the smallest singular value of $R$. According to the lower bound
estimates introduced in [24, 14], the smallest singular value satisfies

$\lambda_{(N+1)/2}(R)\geq\left(\frac{\frac{N-1}{2}}{\|R\|_{F}^{2}}\right)^{\frac{N-1}{4}}|\det(R)|\,,$
(33)

where $\|\cdot\|_{F}$ denotes the Frobenius norm, then use the results from
Lemma A.2 and Lemma A.3 in Appendix, we have the estimate

$\lambda_{(N+1)/2}(R)\geq\left(\frac{N-1}{N+1}\right)^{\frac{N-1}{4}}|\det(R)|=\mathcal{O}(N^{-3/8})\,.$
(34)

Therefore $C_{1}(N,\Omega)\geq
c\min(\frac{1}{(2N+1)\sigma_{s}},\lambda_{(N+1)/2}(RM^{-1}))=\mathcal{O}(N^{-11/8})$.
To estimate the upper bound of $\|\mathbf{s}_{1}\|^{2}$, we compute the row
vector $\mathbf{s}_{1}=[s_{1,1},\dots,s_{1,(N+1)/2}]$ by performing Gauss
elimination on $M$ and obtain

$|s_{1,k}|=\frac{(2k-2)!!}{(2k-1)!!}=\frac{\sqrt{\pi}}{2}\frac{\Gamma(k)}{\Gamma(k+\frac{1}{2})}\,.$
(35)

Use the Gautschi’s inequality [12] that

$\frac{1}{\sqrt{k+\frac{1}{2}}}<\frac{\Gamma(k)}{\Gamma(k+\frac{1}{2})}<\frac{1}{\sqrt{k-\frac{1}{2}}}\,,$
(36)

we immediately find out

$\|\mathbf{s}_{1}\|_{\ell^{2}}^{2}\leq\frac{\pi}{4}\sum_{i=1}^{(N+1)/2}\frac{1}{i-\frac{1}{2}}=\mathcal{O}(1+\log
N)\,.$ (37)

From the result of Theorem 3.1, the stability estimate now can be formulated
as

$\displaystyle\|(\sigma_{a,1}-\sigma_{a,2})\mathbf{s}_{1}\Phi\|_{L^{2}(\Omega)}\leq\mathcal{O}(N^{11/8}(1+\log
N))\left\|\frac{H_{1}-H_{2}}{\Upsilon}\right\|_{L^{2}(\Omega)}\,.$ (38)

∎

Figure 1: Decay of the smallest singular values with respect to the matrix
sizes. The $x$-axis denotes the matrix size $(N+1)/2$. Left: The red (dashed)
line represents the exact smallest singular value of $R$ and the blue (solid)
line represents the estimated lower bound of the smallest singular value of
$R$ through (33). Right: The red (dashed) line represents the exact smallest
singular value of $RM^{-1}$ and the blue (solid) line represents the estimated
lower bound of the smallest singular value of $RM^{-1}$ through (31).

######  Remark 3.3.

It is possible to improve the above estimate by using a sharper bounded for
the Frobenius norm $\|R\|_{F}$ in Lemma A.3. The simple bounds (31) and (33)
are not sharp for the smallest singular value of $R$, see Fig 1. It seems
possible to achieve better estimate through the calculation of $R^{-1}$’s
Frobenius norm by following the technique in [29]. Use the Wielandt-Hoffman
theorem, we can also deduce that $R$ and the Cauchy-Toeplitz matrix with
$(i,j)$-entry as $\frac{1/\pi}{i-j+\frac{1}{2}}$ share the equally distributed
singular values, i.e., $R$’s singular values have a cluster at $1$, see [28].

######  Remark 3.4.

As $N\to\infty$, the above result shows that the stability estimate’s constant
will also grow to infinity, this seems to give a negative answer to the
uniqueness for qPAT with the radiative transport equation. However, such
estimate is only meant for the worst case, since the boundary source could be
chosen arbitrarily. In practice, if the source function $f$ is sufficiently
smooth, the datum with respect to the $SP_{N}$ model
$H=\Upsilon\sigma_{a}\mathbf{s}_{1}\cdot\Phi$ will converge rapidly and the
high order modes will decay sufficiently fast, which could counter the growth
in the constant. This will be the future work.

Reconstruction of $\Upsilon$ only. Suppose the coefficients
$\sigma_{a},\sigma_{s}$ are known and $\Upsilon$ is unknown, then the $SP_{N}$
equation system is completely known and $\Phi$ could be uniquely solved, so we
can reconstruct $\Upsilon$ explicitly by

$\Upsilon=\frac{H(\mathbf{x})}{\sigma_{0}(\mathbf{x})(\mathbf{s}_{1}\cdot\Phi)}\,.$
(39)

The uniqueness and stability estimate will be straightforward, we conclude in
the following theorem without proof.

######  Theorem 3.5.

Under the assumptions $\mathcal{A}$-i to $\mathcal{A}$-iv, suppose
$(\sigma_{a},\sigma_{s})$ are known, $\Upsilon_{1},\Upsilon_{2}$ are two
Grüneisen coefficients, $H_{1},H_{2}$ are the corresponding internal data,
respectively. Then we have the following stability estimate

$\|(\Upsilon_{1}-\Upsilon_{2})\frac{H_{1}}{\Upsilon_{1}}\|_{L^{2}(\Omega)}\leq
C\|(H_{1}-H_{2})/\sigma_{a}\|_{L^{2}(\Omega)}\,,$ (40)

the constant $C$ does not depend on $N$.

Reconstruction of $\sigma_{s}$ only. Suppose the coefficients
$\Upsilon,\sigma_{a}$ are known and $\sigma_{s}$ is unknown, then

$\mathbf{s}_{1}\cdot\Phi=\frac{H(\mathbf{x})}{\Upsilon(\mathbf{x})\sigma_{0}(\mathbf{x})}$
(41)

is known from the measurement $H$. In addition, we also assume that
$\sigma_{s}$ is known on the boundary $\partial\Omega$. When $N=1$, the
reconstruction process of $\sigma_{s}$ will be solving a linear transport
equation [5] for ${\sigma_{s}}^{-1}$, while larger $N$ will introduce extra
nonlinearity from the coupling of solution components. The linearized case of
$SP_{3}$ approximation has been recently studied in [11], which should be able
to generalize to $SP_{N}$ system with the similar technique. In the following,
we assume $N\geq 3$, for the corresponding nonlinear inverse problem, let
$\mathbf{s}_{k}$ be the $k$-th row of $M^{-1}$ and $\mathbf{u}_{k,n}$ be the
$n$-th row of the rank-one matrix $\mathbf{s}_{k}^{T}\mathbf{s}_{k}$. We
reformulate the equations from the bilinear form (13) as

$-\nabla\cdot\left(\frac{1}{\sigma_{s}}\nabla\varphi_{n}\right)+\sigma_{a}\mathbf{p}_{n}\cdot\Phi+\sigma_{s}\mathbf{q}_{n}\cdot\Phi=0,\quad
n=1,\dots,(N+1)/2\,,$ (42)

where the row vectors $\mathbf{p}_{n}$ and $\mathbf{q}_{n}$ are defined by

$\displaystyle\mathbf{p}_{n}$
$\displaystyle=(4n-1)(1-g^{2n-1})\mathbf{u}_{1,n}\,,$ (43)
$\displaystyle\mathbf{q}_{n}$
$\displaystyle=(4n-1)(1-g^{2n-1})\sum_{k=2}^{(N+1)/2}(4k-3)(1-g^{2k-2})\mathbf{u}_{k,n}\,.$

For convenience, we also denote $\mathbf{P}$ and $\mathbf{Q}$ as the
corresponding matrices with $n$-th rows as $\mathbf{p}_{n}$ and
$\mathbf{q}_{n}$, respectively. In the following, we first prove a simple
lemma to estimate the variation in solutions with respect to changes in the
scattering coefficient.

######  Lemma 3.6.

Under the assumptions $\mathcal{A}$-i to $\mathcal{A}$-iv, suppose
$\sigma_{s,1}$ and $\sigma_{s,2}$ are two admissible scattering coefficients
that $\delta\sigma_{s}=\sigma_{s,1}-\sigma_{s,2}$ satisfies
$\delta\sigma_{s}=0$ on $\partial\Omega$. Let
$\Phi_{i}=[\varphi_{i,1},\dots,\varphi_{i,(N+1)/2}]^{T}\in[H^{1}(\Omega)]^{(N+1)/2}$
the solution associated with scattering coefficient $\sigma_{s,i}$, $i=1,2$,
then

$\|\Phi_{1}-\Phi_{2}\|_{[H^{1}(\Omega)]^{(N+1)/2}}\leq
C\frac{N^{23/8}}{3(1-g)}\|\Phi_{1}\|_{[W^{1,\infty}(\Omega)]^{(N+1)/2}}\|\sigma_{s,1}-\sigma_{s,2}\|_{L^{2}(\Omega)}\,,$
(44)

where $C$ is a positive constant independent of $N$.

###### Proof.

Let $B_{i}(\cdot,\cdot)$ the bilinear form in (13) for coefficient pair
$(\sigma_{a},\sigma_{s,i})$, $i=1,2$, then for any test function
$\Psi\in[H^{1}(\Omega)]^{(N+1)/2}$, we have

$B_{1}(\Phi_{1},\Psi)=B_{2}(\Phi_{2},\Psi)\,.$ (45)

Denote
$\delta\Phi=[\delta\varphi_{1},\dots,\delta\varphi_{(N+1)/2}]^{T}:=\Phi_{1}-\Phi_{2}$,
by Cauchy-Schwartz inequality, we have

$\displaystyle B_{2}(\delta\Phi,\delta\Phi)$
$\displaystyle=\sum_{n=1}^{(N+1)/2}\int_{\Omega}\frac{1}{(4n-1)(1-g^{2n-1})}\left[\frac{\delta\sigma_{s}}{\sigma_{s,1}\sigma_{s,2}}\nabla\varphi_{1,n}\cdot\nabla\delta\varphi_{n}-\delta\sigma_{s}\mathbf{q}_{n}\cdot\Phi_{1}\delta\varphi_{n}\right]d\mathbf{x}$
(46)
$\displaystyle\leq\theta\sum_{n=1}^{(N+1)/2}\left[\|\nabla\delta\varphi_{n}\|_{L^{2}(\Omega)}\left\|\frac{\delta\sigma_{s}\nabla\varphi_{1,n}}{\sigma_{s,1}\sigma_{s,2}}\right\|_{L^{2}(\Omega)}+\frac{3\left\|\delta\varphi_{n}\right\|_{L^{2}(\Omega)}}{4n-1}\left\|\delta\sigma_{s}\mathbf{q}_{n}\cdot\Phi_{1}\right\|_{L^{2}(\Omega)}\right]$
$\displaystyle\leq\theta\kappa\|\delta\sigma_{s}\|_{L^{2}(\Omega)}\sum_{n=1}^{(N+1)/2}\|\delta\varphi_{n}\|_{H^{1}(\Omega)}$
$\displaystyle\leq\theta{\kappa\sqrt{\frac{(N+1)}{2}}}\|\delta\sigma_{s}\|_{L^{2}(\Omega)}\|\delta\Phi\|_{[H^{1}(\Omega)]^{(N+1)/2}}\,,$

where the constants $\theta=(3(1-g))^{-1}$, $\kappa=\sup_{n\geq
1}\left\|\nabla\varphi_{1,n}/(\sigma_{s,1}\sigma_{s,2})\right\|_{L^{\infty}(\Omega)}+\frac{3}{4n-1}\left\|\mathbf{q}_{n}\cdot\Phi_{1}\right\|_{L^{\infty}(\Omega)}$.
Since $\sigma_{s,1},\sigma_{s,2}$ are bounded from below by positive
constants, the first term in $\kappa$ is bounded by
$\|\Phi_{1}\|_{[W^{1,\infty}]^{(N+1)/2}}$. The second term needs to estimate
$\sup_{n\geq 1}\frac{3}{4n-1}\|\mathbf{q}_{n}\|_{\ell^{1}}$. From the
definition of $\mathbf{q}_{n}$ in (43), we can deduce that

$\displaystyle\frac{1}{4n-1}\|\mathbf{q}_{n}\|_{\ell^{1}}$
$\displaystyle\leq\sum_{k=2}^{(N+1)/2}(4k-3)\|\mathbf{u}_{k,n}\|_{\ell^{1}}$
(47) $\displaystyle=\sum_{k=2}^{n}(4k-3)\sum_{j=k}^{(N+1)/2}|s_{k,n}s_{k,j}|$
$\displaystyle<\sum_{k=2}^{n}\sum_{j=k}^{(N+1)/2}\frac{(k+\frac{1}{2})(4k-3)}{(2k-1)^{2}}\frac{1}{\sqrt{j-\frac{1}{2}}\sqrt{n-\frac{1}{2}}}\quad\text{(Gaustchi's
inequality) }$
$\displaystyle\leq\frac{1}{\sqrt{n-\frac{1}{2}}}\sum_{k=2}^{n}\sum_{j=k}^{(N+1)/2}\frac{2}{\sqrt{j-\frac{1}{2}}}\quad(\text{since
}\frac{1}{2}(k+1)(4k-3)\leq(2k-1)^{2})$ $\displaystyle=\mathcal{O}(N).$

This implies
$\kappa\leq\mathsf{c}N\|\Phi_{1}\|_{[W^{1,\infty}(\Omega)]^{(N+1)/2}}$ for
certain constant $\mathsf{c}>0$. On the other hand, from the Corollary 3.2,
there exists a constant $C_{1}$ independent of $N$ that

$B_{2}(\delta\Phi,\delta\Phi)\geq
C_{1}N^{-11/8}\|\delta\Phi\|_{[H^{1}(\Omega)]^{(N+1)/2}}^{2}\,.$ (48)

Combine the estimates (48) and (46), we obtain

$\|\delta\Phi\|_{H^{1}(\Omega)}\leq\frac{\mathsf{c}}{C_{1}}\frac{N^{23/8}}{3(1-g)}\|\Phi_{1}\|_{[W^{1,\infty}(\Omega)]^{(N+1)/2}}\|\delta\sigma_{s}\|_{L^{2}}\,.$
(49)

∎

######  Theorem 3.7.

Under the assumptions $\mathcal{A}$-i to $\mathcal{A}$-iv, suppose
$(\sigma_{a},\Upsilon)$ are known, let $\sigma_{s,1}$ and $\sigma_{s,2}$ be
two admissible scattering coefficients with $\sigma_{s,1}=\sigma_{s,2}$ on the
boundary $\partial\Omega$. Let $\Phi_{1}$ and $\Phi_{2}$ the solutions to the
$SP_{N}$ system with scattering coefficients $\sigma_{s,1}$ and $\sigma_{s,2}$
respectively. $H_{1}$ and $H_{2}$ denote the corresponding internal data for
$\sigma_{s,1}$ and $\sigma_{s,2}$ respectively. Then we have the following
estimate

$\int_{\Omega}\mathcal{V}_{N}(\mathbf{x};\lambda)\left(\frac{\sigma_{s,1}-\sigma_{s,2}}{\sigma_{s,2}}\right)^{2}d\mathbf{x}\leq\frac{\mathsf{C}}{\lambda}N^{2}\left\|\frac{H_{1}-H_{2}}{\Upsilon\sigma_{a}}\right\|^{2}_{H^{2}(\Omega)}\,,$
(50)

where $\lambda\in(0,2\underline{c})$ is an arbitrary constant,
$\mathcal{V}_{N}(\mathbf{x};\lambda)$ is

$\displaystyle\mathcal{V}_{N}(\mathbf{x})$
$\displaystyle:=(\sigma_{s,1}+2\underline{c}-\lambda)(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})^{2}+\kappa_{N}\frac{H_{1}}{\Upsilon}(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})-\frac{1}{\sigma_{s,1}}\nabla\frac{H_{1}}{\Upsilon\sigma_{a}}\cdot\nabla(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})-2\overline{c}^{2}\mathcal{Y}\,,$
(51)

with

$\displaystyle\mathcal{Y}$
$\displaystyle:=C_{3}N^{39/8}\|\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1}\|_{L^{\infty}(\Omega)}\|\Phi_{1}\|_{W^{1,\infty}(\Omega)}\,,$
(52) $\displaystyle\kappa_{N}$
$\displaystyle:=\sum_{n=1}^{(N+1)/2}(4n-1)(1-g^{2n-1})s_{1,n}^{2}\,,$

the constants $\mathsf{C},C_{3}$ are independent of $N$. When
$\mathcal{V}_{N}(\mathbf{x};\lambda)>0$, $\forall\mathbf{x}\in\Omega$, then
$H_{1}=H_{2}$ a.e. on $\Omega$ implies $\sigma_{s,1}=\sigma_{s,2}$ a.e..

###### Proof.

For $n=1,\dots,(N+1)/2$, we have the following $SP_{N}$ systems for $\Phi_{1}$
and $\Phi_{2}$,

$\displaystyle-\nabla\cdot\left(\frac{1}{\sigma_{s,1}}\nabla\varphi_{1,n}\right)+\sigma_{a}\mathbf{p}_{n}\cdot\Phi_{1}+\sigma_{s,1}\mathbf{q}_{n}\cdot\Phi_{1}$
$\displaystyle=0\,,$ (53)
$\displaystyle-\nabla\cdot\left(\frac{1}{\sigma_{s,2}}\nabla\varphi_{2,n}\right)+\sigma_{a}\mathbf{p}_{n}\cdot\Phi_{2}+\sigma_{s,2}\mathbf{q}_{n}\cdot\Phi_{2}$
$\displaystyle=0\,,$

with the mixed boundary condition (11). Since $\Upsilon$ and $\sigma_{a}$ are
known, the measurements are given by

$\mathbf{s}_{1}\cdot\Phi_{i}=\frac{H_{i}}{\Upsilon\sigma_{a}},\quad i=1,2\,.$
(54)

Therefore multiply (53) with $s_{1,n}$ and take summation over $n$, we get

$-\nabla\cdot\left(\frac{1}{\sigma_{s,i}}\nabla\frac{H_{i}}{\Upsilon\sigma_{a}}\right)+\mathbf{s}_{1}\cdot\left(\sigma_{a}\mathbf{P}\Phi_{i}+\sigma_{s,i}\mathbf{Q}\Phi_{i}\right)=0,\quad
i=1,2\,.$ (55)

Take the difference between equations (55) with $\sigma_{s,1}$ and
$\sigma_{s,2}$, respectively. Let
$\delta\sigma_{s}=\sigma_{s,1}-\sigma_{s,2}$, $\delta\Phi=\Phi_{1}-\Phi_{2}$
and $\delta H=H_{1}-H_{2}$, then

$-\nabla\cdot\left(\frac{-\delta\sigma_{s}}{\sigma_{s,1}\sigma_{s,2}}\nabla\frac{H_{1}}{\Upsilon\sigma_{a}}\right)-\nabla\cdot\left(\frac{1}{\sigma_{s,2}}\nabla\frac{\delta
H}{\Upsilon\sigma_{a}}\right)+\mathbf{s}_{1}\cdot\left(\sigma_{a}\mathbf{P}\delta\Phi+\delta\sigma_{s}\mathbf{Q}\Phi_{1}+\sigma_{s,2}\mathbf{Q}\delta\Phi\right)=0\,.$
(56)

Use the following identity,

$\frac{\delta\sigma_{s}}{\sigma_{s,2}}\nabla\cdot\left(\frac{\delta\sigma_{s}}{\sigma_{s,2}}\frac{1}{\sigma_{s,1}}\nabla\frac{H_{1}}{\Upsilon\sigma_{a}}\right)=\frac{1}{2}\left(\frac{\delta\sigma_{s}}{\sigma_{s,2}}\right)^{2}\nabla\cdot\left(\frac{1}{\sigma_{s,1}}\nabla\frac{H_{1}}{\Upsilon\sigma_{a}}\right)+\frac{1}{2}\nabla\cdot\left(\left(\frac{\delta\sigma_{s}}{\sigma_{s,2}}\right)^{2}\frac{1}{\sigma_{s,1}}\nabla\frac{H_{1}}{\Upsilon\sigma_{a}}\right)\,,$
(57)

we multiply (56) with $\delta\sigma_{s}/\sigma_{s,2}$, then

$\displaystyle\frac{1}{2}\left(\frac{\delta\sigma_{s}}{\sigma_{s,2}}\right)^{2}\nabla\cdot\left(\frac{1}{\sigma_{s,1}}\nabla\frac{H_{1}}{\Upsilon\sigma_{a}}\right)+\frac{1}{2}\nabla\cdot\left(\left(\frac{\delta\sigma_{s}}{\sigma_{s,2}}\right)^{2}\frac{1}{\sigma_{s,1}}\nabla\frac{H_{1}}{\Upsilon\sigma_{a}}\right)$
(58)
$\displaystyle-\frac{\delta\sigma_{s}}{\sigma_{s,2}}\nabla\cdot\left(\frac{1}{\sigma_{s,2}}\nabla\frac{\delta
H}{\Upsilon\sigma_{a}}\right)+\frac{\delta\sigma_{s}}{\sigma_{s,2}}\mathbf{s}_{1}\cdot\left(\sigma_{a}\mathbf{P}\delta\Phi+\delta\sigma_{s}\mathbf{Q}\Phi_{1}+\sigma_{s,2}\mathbf{Q}\delta\Phi\right)=0\,.$

The first term can be replaced from (55) that

$\frac{1}{2}\left(\frac{\delta\sigma_{s}}{\sigma_{s,2}}\right)^{2}\nabla\cdot\left(\frac{1}{\sigma_{s,1}}\nabla\frac{H_{1}}{\Upsilon\sigma_{a}}\right)=\frac{1}{2}\left(\frac{\delta\sigma_{s}}{\sigma_{s,2}}\right)^{2}\mathbf{s}_{1}\cdot\left(\sigma_{a}\mathbf{P}\Phi_{1}+\sigma_{s,1}\mathbf{Q}\Phi_{1}\right)\,,$
(59)

then combine (58) and (59), multiply the test function
$\psi=\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1}\in H^{1}(\Omega)$ to (58) and
integrate over $\Omega$. Notice that $\delta\sigma_{s}=0$ on $\partial\Omega$,
we obtain

$\displaystyle\frac{1}{2}\int_{\Omega}\left(\frac{\delta\sigma_{s}}{\sigma_{s,2}}\right)^{2}\left[\mathbf{s}_{1}\cdot\left(\sigma_{a}\mathbf{P}+\sigma_{s,1}\mathbf{Q}+2\sigma_{s,2}\mathbf{Q}\right)\Phi_{1}(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})-\frac{1}{\sigma_{s,1}}\nabla\frac{H_{1}}{\Upsilon\sigma_{a}}\cdot\nabla(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})\right]d\mathbf{x}$
(60)
$\displaystyle-\int_{\Omega}\frac{\delta\sigma_{s}}{\sigma_{s,2}}\left[\nabla\cdot\left(\frac{1}{\sigma_{s,2}}\nabla\frac{\delta
H}{\Upsilon\sigma_{a}}\right)-\mathbf{s}_{1}\cdot(\sigma_{a}\mathbf{P}\delta\Phi+\sigma_{s,2}\mathbf{Q}\delta\Phi)\right](\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})d\mathbf{x}=0\,.$

In the next, observe that

$\displaystyle\mathbf{s}_{1}\cdot\mathbf{P}=\sum_{n=1}^{(N+1)/2}(4n-1)(1-g^{2n-1})s_{1,n}\mathbf{u}_{1,n}=\left(\sum_{n=1}^{(N+1)/2}(4n-1)(1-g^{2n-1})s_{1,n}^{2}\right)\mathbf{s}_{1}\,,$
(61)

therefore in (60), we can replace

$\displaystyle\sigma_{a}\mathbf{s}_{1}\cdot\mathbf{P}\delta\Phi=\left(\sum_{n=1}^{(N+1)/2}(4n-1)(1-g^{2n-1})s_{1,n}^{2}\right)\frac{\delta
H}{\Upsilon}\,,$ (62)
$\displaystyle\sigma_{a}\mathbf{s}_{1}\cdot\mathbf{P}\Phi_{1}=\left(\sum_{n=1}^{(N+1)/2}(4n-1)(1-g^{2n-1})s_{1,n}^{2}\right)\frac{H_{1}}{\Upsilon}\,.$

Define the constant
$\kappa_{N}:=\sum_{n=1}^{(N+1)/2}(4n-1)(1-g^{2n-1})s_{1,n}^{2}$, since
$|s_{1,n}|^{2}=\Theta(n^{-1})$, then $\kappa_{N}=\Theta(N)$, the equation (60)
can be further reduced to

$\displaystyle\frac{1}{2}\int_{\Omega}\left(\frac{\delta\sigma_{s}}{\sigma_{s,2}}\right)^{2}\left[\left(\sigma_{s,1}+2\sigma_{s,2}\right)(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})^{2}+\kappa_{N}\frac{H_{1}}{\Upsilon}(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})-\frac{1}{\sigma_{s,1}}\nabla\frac{H_{1}}{\Upsilon\sigma_{a}}\cdot\nabla(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})\right]d\mathbf{x}$
(63)
$\displaystyle=\int_{\Omega}\frac{\delta\sigma_{s}}{\sigma_{s,2}}\left[\nabla\cdot\left(\frac{1}{\sigma_{s,2}}\nabla\frac{\delta
H}{\Upsilon\sigma_{a}}\right)-\kappa_{N}\frac{\delta
H}{\Upsilon}-\sigma_{s,2}\left(\mathbf{s}_{1}\cdot\mathbf{Q}\delta\Phi\right)\right](\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})d\mathbf{x}\,.$

Due to Lemma A.4, $\|\mathbf{s}_{1}\cdot\mathbf{Q}\|_{\ell^{p}}\leq
C_{2}N^{3/2+1/p}$ for certain constant $C_{2}$ independent of $N$, then
combine with Lemma 3.6, the second term on right-hand-side of (63) is bounded
by

$\displaystyle\big{|}\int_{\Omega}\delta\sigma_{s}\left(\mathbf{s}_{1}\cdot\mathbf{Q}\delta\Phi\right)(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})d\mathbf{x}\big{|}$
$\displaystyle\leq\|\mathbf{s}_{1}\cdot\mathbf{Q}\|_{\ell^{2}}\left\|\delta\sigma_{s}\right\|_{L^{2}(\Omega)}\|\delta\Phi\|_{L^{2}(\Omega)}\|\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1}\|_{L^{\infty}(\Omega)}$
$\displaystyle\leq
C_{2}N^{2}\left\|\delta\sigma_{s}\right\|_{L^{2}(\Omega)}\|\delta\Phi\|_{L^{2}(\Omega)}\|\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1}\|_{L^{\infty}(\Omega)}$
$\displaystyle\leq
C_{3}N^{39/8}\|\Phi_{1}\|_{W^{1,\infty}(\Omega)}\|\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1}\|_{L^{\infty}(\Omega)}\left\|\delta\sigma_{s}\right\|_{L^{2}(\Omega)}^{2}\,,$

where the constant $C_{3}$ is independent of $N$ as well. Let $\mathcal{X}$
and $\mathcal{Y}$ denote the following quantities

$\displaystyle\mathcal{X}$
$\displaystyle:=\kappa_{N}\frac{H_{1}}{\Upsilon}(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})-\frac{1}{\sigma_{s,1}}\nabla\frac{H_{1}}{\Upsilon\sigma_{a}}\cdot\nabla(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})\,,$
(64) $\displaystyle\mathcal{Y}$
$\displaystyle:=C_{3}N^{39/8}\|\Phi_{1}\|_{W^{1,\infty}(\Omega)}\|\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1}\|_{L^{\infty}(\Omega)}\,,$

then we obtain the following inequality,

$\displaystyle\frac{1}{2}\int_{\Omega}\left(\frac{\delta\sigma_{s}}{\sigma_{s,2}}\right)^{2}\left[(\sigma_{s,1}+2\sigma_{s,2})(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})^{2}+\mathcal{X}-2\sigma_{s,2}^{2}\mathcal{Y}\right]d\mathbf{x}$
(65)
$\displaystyle\leq\int_{\Omega}\frac{\delta\sigma_{s}}{\sigma_{s,2}}\left(\nabla\cdot\left(\frac{1}{\sigma_{s,2}}\nabla\frac{\delta
H}{\Upsilon\sigma_{a}}\right)-\kappa_{N}\frac{\delta
H}{\Upsilon}\right)(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})d\mathbf{x}$
$\displaystyle\leq\frac{\lambda}{2}\int_{\Omega}\left[\frac{\delta\sigma_{s}}{\sigma_{s,2}}(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})\right]^{2}d\mathbf{x}+\frac{1}{2\lambda}\int_{\Omega}\left(\nabla\cdot\left(\frac{1}{\sigma_{s,2}}\nabla\frac{\delta
H}{\Upsilon\sigma_{a}}\right)-\kappa_{N}\frac{\delta
H}{\Upsilon}\right)^{2}d\mathbf{x}\,,$

where $\lambda>0$ is an arbitrary number and the last inequality has used the
AM-GM inequality. For the uniqueness, we let $\delta H=0$ in above inequality,
then it becomes

$\int_{\Omega}\left(\frac{\delta\sigma_{s}}{\sigma_{s,2}}\right)^{2}\left[(\sigma_{s,1}+2\sigma_{s,2}-\lambda)(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})^{2}+\mathcal{X}-2\sigma_{s,2}^{2}\mathcal{Y}\right]d\mathbf{x}\leq
0\,.$ (66)

Recall that $\underline{c}\leq\sigma_{s,2}\leq\overline{c}$, therefore if

$(\sigma_{s,1}+2\underline{c}-\lambda)(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})^{2}+\kappa_{N}\frac{H_{1}}{\Upsilon}(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})-\frac{1}{\sigma_{s,1}}\nabla\frac{H_{1}}{\Upsilon\sigma_{a}}\cdot\nabla(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})>2\overline{c}^{2}\mathcal{Y}\,,$
(67)

we could conclude that $\delta\sigma_{s}=0$ a.e.. The stability estimate is
straightforward by noticing $\kappa_{N}=\mathcal{O}(N)$. ∎

######  Remark 3.8.

As $N\to\infty$, the requirement that $\mathcal{V}(x;\lambda)>0$ could be
difficult to fulfill since the growth of $\mathcal{Y}$ is much faster than the
other terms. The estimate could be greatly improved by giving a tighter bound
to
$\int_{\Omega}\delta\sigma_{s}(\mathbf{s}_{1}\cdot\mathbf{Q}\delta\Phi)(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{1})d\mathbf{x}$,
for instance, estimate the Frechét derivative of
$\frac{\delta\Phi}{\delta\sigma_{s}}$.

In the next, we focus on two important cases of simultaneous reconstructions:
$(\Upsilon,\sigma_{a})$ and $(\sigma_{a},\sigma_{s})$ with multiple
illumination sources. In $SP_{1}$ , one can only reconstruct any two
coefficients with the knowledge about the third one [5, 25] and it is
impossible to recover all of them unless additional information is provided.
In $SP_{N}$, it is still unclear whether or not all of the coefficients can be
recovered uniquely.

Reconstruction of $\sigma_{a}$ and $\Upsilon$. In this case, we consider the
simultaneous reconstruction of both $\sigma_{a}$ and $\Upsilon$ with multiple
sources $f_{j}$, $1\leq j\leq J$ ($J\geq 2$). We denote $H_{j}$ the
corresponding measurement for $f_{j}$ from qPAT. The simplest case $SP_{1}$ is
studied in [5, 25] and linearized case of $SP_{3}$ is discussed in [11]. The
key observation is that the ratio of two measurements is independent of the
coefficients, which only implicitly depends on $\sigma_{a}$. We first consider
the linearized setting, suppose the scattering coefficient $\sigma_{s}$ and
the background absorption and Grüneisen coefficients $\sigma_{a},\Upsilon$ are
known and the perturbations are $\delta\sigma_{s},\delta\Upsilon$. For each
$1\leq j\leq J$, suppose $\Phi_{j}$ is the background solution with source
$f_{j}$, the perturbation in the solution is denoted by $\delta\Phi_{j}$, then
by linearizing (42),
$\delta\Phi_{j}=[\delta\varphi_{j,1},\dots,\delta\varphi_{(N+1)/2}]^{\text{T}}$
satisfies the linearized $SP_{N}$ system

$-\nabla\cdot\left(\frac{1}{\sigma_{s}}\nabla\delta\varphi_{j,n}\right)+\sigma_{a}\mathbf{p}_{n}\cdot\delta\Phi_{j}+\sigma_{s}\mathbf{q}_{n}\cdot\delta\Phi_{j}=-\delta\sigma_{a}\mathbf{p}_{n}\cdot\Phi_{j},\quad
n=1,\dots,(N+1)/2\,.$ (68)

For any pair of indices $1\leq i<j\leq J$, we define the following quantity

$\mathcal{H}_{ij}:=(\mathbf{s}_{1}\cdot\Phi_{i})\frac{\delta
H_{j}}{\sigma_{a}\Upsilon}-(\mathbf{s}_{1}\cdot\Phi_{j})\frac{\delta
H_{i}}{\sigma_{a}\Upsilon}=(\mathbf{s}_{1}\cdot\Phi_{i})(\mathbf{s}_{1}\cdot\delta\Phi_{j})-(\mathbf{s}_{1}\cdot\Phi_{j})(\mathbf{s}_{1}\cdot\delta\Phi_{i})\,,$
(69)

which is known and only depends on $\delta\sigma_{a}$ and independent of
$\Upsilon$. Our reconstruction will be a natural two-step process, first solve
$\delta\sigma_{a}$ from the crossing quantity $\mathcal{H}_{ij}$ ($1\leq
i<j\leq J$), then solve $\delta\Phi_{i}$ using the recovered
$\delta\sigma_{a}$, finally if $\sigma_{a}>0$ find $\delta\Upsilon$ through

$\delta\Upsilon=\frac{1}{J}\sum_{j=1}^{J}\frac{\delta
H_{j}-\Upsilon\sigma_{a}\mathbf{s}_{1}\cdot\delta\Phi_{j}-\Upsilon\delta\sigma_{a}\mathbf{s}_{1}\cdot\Phi_{j}}{\sigma_{a}\mathbf{s}_{1}\cdot\Phi_{j}}\,.$
(70)

By taking linearization over the bilinear form (13), for any test function
$\Psi\in[H^{1}(\Omega)]^{(N+1)/2}$,

$B(\delta\Phi_{j},\Psi)=-\int_{\Omega}\left[\delta\sigma_{a}(\mathbf{s}_{1}\cdot\Phi_{j})\mathbf{s}_{1}\right]\cdot\Psi
d\mathbf{x}\,.$ (71)

Take $\Psi=\delta\Phi_{j}$, we can easily conclude that the following
estimates from the coerciveness of $B(\cdot,\cdot)$:

$\displaystyle\|\delta\Phi_{j}\|_{[H^{1}(\Omega)]^{(N+1)/2}}$
$\displaystyle=\mathcal{O}\left(N^{11/8}\|\delta\sigma_{a}(\mathbf{s}_{1}\cdot\Phi_{j})\|_{L^{2}(\Omega)}\|\mathbf{s}_{1}\|_{\ell^{2}}\right)$
(72) $\displaystyle=\mathcal{O}\left(N^{11/8}(1+\log
N)\|\delta\sigma_{a}(\mathbf{s}_{1}\cdot\Phi_{j})\|_{L^{2}(\Omega)}\right)\,,$

On the other hand, multiply (68) with $s_{1,n}$ and sum over $n$,

$\delta\sigma_{a}\kappa_{N}(\mathbf{s}_{1}\cdot\Phi_{j})=-\nabla\cdot\left(\frac{1}{\sigma_{s}}\nabla(\mathbf{s}_{1}\cdot\delta\Phi_{j})\right)+\sigma_{a}\kappa_{N}\mathbf{s}_{1}\cdot\delta\Phi_{j}+\sigma_{s}\mathbf{s}_{1}\cdot\mathbf{Q}\delta\Phi_{j}\,,$
(73)

therefore we have a straightforward estimate

$\|\delta\sigma_{a}(\mathbf{s}_{1}\cdot\Phi_{j})\|_{L^{2}(\Omega)}=\mathcal{O}\left(\frac{\max(\kappa_{N}\|\mathbf{s}_{1}\|_{\ell^{2}},\|\mathbf{s}_{1}\cdot\mathbf{Q}\|_{\ell^{2}})}{\kappa_{N}}\|\delta\Phi_{j}\|_{[H^{2}(\Omega)]^{(N+1)/2}}\right)\,,$
(74)

where the constant
$\kappa_{N}=\sum_{n=1}^{(N+1)/2}(4n-1)(1-g^{2n-1})s_{1,n}^{2}=\Theta(N)$ and
$\|\mathbf{s}_{1}\cdot\mathbf{Q}\|_{\ell^{2}}\leq\mathcal{O}(N^{2})$ from
Lemma A.4. These two estimates imply that

$c^{-1}N^{-11/8}(1+\log
N)^{-1}\|\delta\Phi_{j}\|_{[H_{1}(\Omega)]^{(N+1)/2}}\leq\|\delta\sigma_{a}(\mathbf{s}_{1}\cdot\Phi_{j})\|_{L^{2}(\Omega)}\leq
cN\|\delta\Phi_{j}\|_{[H^{2}(\Omega)]^{(N+1)/2}}.$ (75)

for some constant $c>1$ independent of $N$.

Therefore if there exists two source distinct functions $f_{i},f_{j}$ such
that the linear mapping $\delta\sigma_{a}\mapsto\mathcal{H}_{ij}$ is
invertible and $\mathbf{s}_{1}\cdot\Phi_{j}\neq 0$ over $\Omega$, then one can
recover both $\delta\sigma_{a}$ and $\delta\Phi_{j}$ from $\mathcal{H}_{ij}$
uniquely. In general such problem is ill-posed due to the compactness of the
mapping $\delta\sigma_{a}\mapsto\mathcal{H}_{ij}$, numerical reconstruction of
$\delta\sigma_{a}$ can be done through the following minimization formulation
with regularization,

$\delta\sigma_{a}^{\ast}=\operatorname*{arg\,min}_{\delta\sigma_{a}}\sum_{1\leq
i<j\leq
J}\left[\Big{\|}\mathcal{H}_{ij}-\left[(\mathbf{s}_{1}\cdot\Phi_{i})(\mathbf{s}_{1}\cdot\delta\Phi_{j})-(\mathbf{s}_{1}\cdot\Phi_{j})(\mathbf{s}_{1}\cdot\delta\Phi_{i})\right]\Big{\|}_{L^{2}(\Omega)}^{2}\right]+\alpha\Big{\|}\nabla\delta\sigma_{a}\Big{\|}_{L^{2}(\Omega)}^{2}\,,$

where $\alpha$ is the regularization parameter.

Particularly, when the background absorption coefficient $\sigma_{a}=0$ or
negligible, then we approximately have $\delta
H_{i}=\delta\sigma_{a}\Upsilon\mathbf{s}_{1}\cdot\Phi_{i}$, which does not
contain the perturbation $\delta\Upsilon$, in this case, we can only
reconstruct $\delta\sigma_{a}$, the stability estimate is similar to the
Theorem 3.1.

Without linearization, we take the ratio of two data sets $H_{i}$ and $H_{j}$,
then

$\frac{H_{i}}{H_{j}}=\frac{\mathbf{s}_{1}\cdot\Phi_{i}}{\mathbf{s}_{1}\cdot\Phi_{j}}\,,$
(76)

which only depends on $\sigma_{a}$, therefore our reconstruction strategy is
similar to the linearized case. First, try to solve the minimization problem:

$\sigma_{a}^{\ast}=\operatorname*{arg\,min}_{\sigma_{a}}\sum_{1\leq i<j\leq
J}\|H_{i}(\mathbf{s}_{1}\cdot\Phi_{j})-H_{j}(\mathbf{s}_{1}\cdot\Phi_{i})\|^{2}_{L^{2}(\Omega)}+\alpha\|\nabla\sigma_{a}\|_{L^{2}(\Omega)}^{2}\,.$
(77)

Then compute $\Upsilon^{\ast}=\frac{1}{J}\sum_{1\leq i\leq
J}\frac{H_{i}}{\sigma_{a}^{\ast}\mathbf{s}_{1}\cdot\Phi_{i}}$ with the
reconstructed $\sigma_{a}^{\ast}$.

Reconstruction of $\sigma_{s}$ and $\sigma_{a}$. We consider the simultaneous
reconstruction of both $\sigma_{s}$ and $\sigma_{a}$ from multiple sources
$f_{j}$, $1\leq j\leq J$ provided that $\Upsilon$ is known. Similar to the
previous case, we denote $H_{j}$ the measurement for $f_{j}$ from the qPAT
experiments. Under the linearized setting, let $\sigma_{a}$ and $\sigma_{s}$
the background absorption and scattering coefficients, the corresponding
perturbations are $\delta\sigma_{a}$ and $\delta\sigma_{s}$. For each source
$f_{j}$, let $\Phi_{j}$ the background solution and the perturbation in the
solution is $\delta\Phi_{j}$, the corresponding perturbation in the
measurement is $\delta H_{j}$. Linearize the variational form (20), we obtain
the following equation,

$\displaystyle\widetilde{B}(\delta\Phi_{j},\Psi)=-\int_{\Omega}\frac{\delta
H_{j}}{\Upsilon}(\mathbf{s}_{1}\cdot\Psi)d\mathbf{x}+\sum_{n=1}^{(N+1)/2}\int_{\Omega}\frac{\delta\sigma_{s}}{(4n-1)(1-g^{2n-1})}\left[\frac{1}{\sigma_{s}^{2}}\nabla\varphi_{1,n}\cdot\nabla\psi_{n}-\mathbf{q}_{n}\cdot\Phi_{j}\psi_{n}\right]d\mathbf{x}\,,$

where the bilinear form $\widetilde{B}(\cdot,\cdot)$ is from (23) and
$\mathbf{q}_{n}$ is defined in (43). Hence the perturbation $\delta\Phi_{j}$
only linearly depends on $\delta\sigma_{s}$. On the other hand, since $\delta
H_{j}/\Upsilon=\delta\sigma_{a}\mathbf{s}_{1}\cdot\Phi_{j}+\sigma_{a}\mathbf{s}_{1}\cdot\delta\Phi_{j}$,
the crossing quantity
$\mathcal{H}_{ij}=\frac{1}{\sigma_{a}\Upsilon^{2}}(H_{i}\delta
H_{j}-H_{j}\delta
H_{i})=(\mathbf{s}_{1}\cdot\Phi_{i})(\mathbf{s}_{1}\cdot\delta\Phi_{j})-(\mathbf{s}_{1}\cdot\Phi_{j})(\mathbf{s}_{1}\cdot\delta\Phi_{i})$
only linearly depends on $\delta\sigma_{s}$, therefore we first try to
reconstruct $\delta\sigma_{s}$ from $\mathcal{H}_{ij}$, then find
$\delta\Phi_{j}$ and recover $\delta\sigma_{a}$ using

$\delta\sigma_{a}=\frac{1}{J}\sum_{j=1}^{J}\left[\left(\frac{\delta
H_{j}}{\Upsilon}-\sigma_{a}\mathbf{s}_{1}\cdot\delta\Phi_{j}\right)/(\mathbf{s}_{1}\cdot\Phi_{j})\right]\,.$
(78)

Similar to the previous case, the uniqueness is immediate if the crossing term
$\mathcal{H}_{ij}$ as a linear functional of $\delta\sigma_{s}$ is uniquely
solvable and $\mathbf{s}_{1}\cdot\Phi_{j}\neq 0$ over $\Omega$. However since
$\delta\sigma_{s}\mapsto\delta\Phi_{j}$ is a compact mapping, the inverse
problem is ill-posed. Numerically, we consider the following $L^{2}$
optimization formulation with regularization:

$\delta\sigma_{s}^{\ast}=\operatorname*{arg\,min}_{\delta\sigma_{s}}\sum_{1\leq
i<j\leq
J}\left[\Big{\|}\mathcal{H}_{ij}-\left[(\mathbf{s}_{1}\cdot\Phi_{i})(\mathbf{s}_{1}\cdot\delta\Phi_{j})-(\mathbf{s}_{1}\cdot\Phi_{j})(\mathbf{s}_{1}\cdot\delta\Phi_{i})\right]\Big{\|}_{L^{2}(\Omega)}^{2}\right]+\alpha\Big{\|}\nabla\delta\sigma_{s}\Big{\|}_{L^{2}(\Omega)}^{2}\,.$

Additionally, if we are provided _a priori_ estimate on the perturbation
$\delta\sigma_{a}$ that
$\|\frac{\delta\sigma_{a}H_{i}}{\Upsilon\sigma_{a}^{2}}\|_{H^{2}(\Omega)}\leq\varepsilon\ll
1$ for certain $1\leq i\leq J$, then linearize the equation (55) for
$\Phi_{i}$, we obtain

$-\nabla\cdot\left(\frac{-\delta\sigma_{s}}{\sigma_{s}^{2}}\nabla\mathbf{s}_{1}\cdot\Phi_{i}\right)-\nabla\cdot\left(\frac{1}{\sigma_{s}}\nabla\mathbf{s}_{1}\cdot\delta\Phi_{i}\right)+\kappa_{N}\frac{\delta
H_{i}}{\Upsilon}+\mathbf{s}_{1}\cdot\left(\delta\sigma_{s}\mathbf{Q}\Phi_{i}+\sigma_{s}\mathbf{Q}\delta\Phi_{i}\right)=0\,,$
(79)

where $\mathbf{Q}$ is defined in (43) and
$\kappa_{N}=\sum_{n=1}^{(N+1)/2}(4n-1)(1-g^{2n-1})s_{1,n}^{2}$. Following the
similar approach in Theorem 3.7,

$\displaystyle\frac{1}{2}\int_{\Omega}\left(\frac{\delta\sigma_{s}}{\sigma_{s}}\right)^{2}\left[\kappa_{N}\frac{H_{i}}{\Upsilon}(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{i})+3\sigma_{s}(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{i})^{2}-\frac{1}{\sigma_{s}}\nabla(\mathbf{s}_{1}\cdot\Phi_{i})\cdot\nabla(\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{i})\right]d\mathbf{x}$
(80)
$\displaystyle=\int_{\Omega}\frac{\delta\sigma_{s}}{\sigma_{s}}\left[\nabla\cdot\left(\frac{1}{\sigma_{s}}\nabla\mathbf{s}_{1}\cdot\delta\Phi_{i}\right)-\kappa_{N}\frac{\delta
H_{i}}{\Upsilon}-\sigma_{s}\left(\mathbf{s}_{1}\cdot\mathbf{Q}\delta\Phi_{i}\right)\right](\mathbf{s}_{1}\cdot\mathbf{Q}\Phi_{i})d\mathbf{x}\,,$

Replace $\mathbf{s}_{1}\cdot\delta\Phi_{i}=(\delta
H_{i}-\delta\sigma_{a}H_{i}/\sigma_{a})/(\Upsilon\sigma_{a})$ and we
immediately get the following estimate from the argument of Theorem 3.7 that

$\displaystyle\int_{\Omega}\mathcal{V}_{N}(\mathbf{x})\left(\frac{\delta\sigma_{s}}{\sigma_{s}}\right)^{2}d\mathbf{x}$
$\displaystyle\leq\mathsf{C}\left(N^{2}\left\|\frac{\delta
H_{i}}{\Upsilon\sigma_{a}}\right\|^{2}_{H^{2}(\Omega)}+\left\|\frac{H_{i}\delta\sigma_{a}}{\Upsilon\sigma_{a}^{2}}\right\|^{2}_{H^{2}(\Omega)}\right)$
(81) $\displaystyle\leq\mathsf{C}\left(N^{2}\left\|\frac{\delta
H_{i}}{\Upsilon\sigma_{a}}\right\|^{2}_{H^{2}(\Omega)}+\varepsilon^{2}\right)\,,$

where $\mathcal{V}_{N}(\mathbf{x})$ is the same as in the Theorem 3.7,
$\mathsf{C}$ is a constant independent of $N$.

##  4 Numerical experiments

In this section, we perform our numerical experiments in two phases: (i).
Assume the true model is certain $SP_{N}$ system and then reconstruct the
coefficients with exactly the same model; (ii). Assume the true model is
either radiative transport equation or certain high order $SP_{N}$ system,
then the reconstruction is performed over a low order $SP_{N}$ system.

In all the following numerical experiments, we use the unit square in 2D as
our domain $\Omega$. It is worthwhile to notice that all the previous
arguments are meant for 3D only, the 2D experiments here should be interpreted
as special cases (e.g. infinite tube) of 3D, assuming the solution is
independent of the third dimension. The forward problem is solved though
finite element method on a sufficiently fine mesh and the inverse problem is
solved on a different mesh to avoid inverse crime. The source code for the
numerical experiments is hosted on
GitHub111https://github.com/lowrank/spn_qpat.

###  4.1 Experiment setting

In the following numerical experiments, we will use the boundary source
functions $f_{1}(x,y)=1+x$, $f_{2}(x,y)=1+\sin(4\pi x)$ and the anisotropy
constant $g=0.8$. The coefficients $(\sigma_{a},\sigma_{s},\Upsilon)$ are
selected from the following variable set, see Fig 2.

Figure 2: Coefficient Set. From left to right: absorption coefficient
$\sigma_{a}$ (Shepp-Logan phantom), scattering coefficient $\sigma_{s}$
(Jaszczak phantom), Grüneisen coefficient $\Upsilon$.

###  4.2 Validation of approximation

Before we start to run any of the numerical experiments, we need to verify
that if our $SP_{N}$ model is valid under these settings, which means modeling
error should not be dominating (in practice there are noises in data).
Therefore it is important to compare the quantity
$\phi_{0}^{N}=\mathbf{s}_{1}\cdot\Phi=\frac{H}{\sigma_{a}\Upsilon}$ with the
solution’s angular average
$U(\mathbf{x})=\int_{\mathbb{S}^{d-1}}u(\mathbf{x},\mathbf{v})d\mathbf{v}$ for
the radiative transport equation (1). We summarize the $L^{2}$ relative
errors: $\frac{\|\phi_{0}(\mathbf{x})-U(\mathbf{x})\|}{\|U(\mathbf{x})\|}$
with respect to different $SP_{N}$ models in Tab 1.

Table 1: Relative $L^{2}$ error of the $\phi_{0}^{N}$ from $SP_{N}$ equations with the angular averaged solution $U(\mathbf{x})$ from radiative transport equation for boundary sources $f_{1}$ and $f_{2}$. | $SP_{1}$ | $SP_{3}$ | $SP_{5}$ | $SP_{7}$ | $SP_{9}$ | $SP_{11}$ | $SP_{13}$ | $SP_{15}$ | $SP_{17}$
---|---|---|---|---|---|---|---|---|---
$f_{1}$ | 1.93% | 2.28% | 2.25% | 2.24% | 2.23% | 2.23% | 2.23% | 2.23% | 2.23%
$f_{2}$ | 5.24% | 3.98% | 3.76% | 3.75% | 3.77% | 3.78% | 3.79% | 3.79% | 3.79%

From the table Tab 1, we can see that the modeling error with respect to the
$SP_{N}$ equations indeed stay at a relative low level with the selected
coefficient set. The reason that such modeling error is not converging to zero
might partly attribute to the simplification
$\sigma_{n}\simeq(1-g^{n})\sigma_{s}$ in the computation instead of using
$\sigma_{a}+(1-g^{n})\sigma_{s}$.

It is also informative to look at the convergence rate of $\phi_{0}^{N}$ from
the $SP_{N}$ equations with respect to growing $N$, see Fig 3. If the error
converges sufficiently fast (e.g. exponentially), then we will obtain the
uniqueness of reconstruction for the case that $N=\infty$. However, the theory
about the convergence is still an open problem.

Figure 3: The convergence of $\phi_{0}^{N}$ with respect to the order $N$.
Left: with source function $f_{1}$. Right: with source function $f_{2}$.

###  4.3 Reconstruction of $\sigma_{a}$ only

In this numerical experiment, we are using the algorithm introduced in Section
3 to reconstruct the absorption coefficient $\sigma_{a}$ only. We consider two
scenarios for the reconstruction: (i) The datum $H$ is generated from certain
$SP_{N}$ model. (ii) The datum $H$ is generated from the radiative transport
model. The result is summarized in the following Tab 2 and Tab 3. For all the
reconstructions in the tables below, we have contaminated the datum $H$ with
multiplicative random noises pointwisely by
$H^{\ast}=H(1+\gamma\texttt{random})$ with parameter $\gamma=5\%$ regarded as
the noise level and random is the uniform distributed random variable on
$[-1,1]$.

Table 2: Relative $L^{2}$ error of the reconstructed $\sigma_{a}$ for source function $f_{1}$ with different generating models and reconstruction models. The row label represents the generating model and the column label represents the reconstruction model. | $SP_{1}$ | $SP_{3}$ | $SP_{5}$ | $SP_{7}$ | $SP_{9}$ | $SP_{11}$ | $SP_{13}$ | $SP_{15}$ | $SP_{17}$
---|---|---|---|---|---|---|---|---|---
$SP_{1}$ | 2.89% | 3.13% | 3.15% | 3.14% | 3.14% | 3.14% | 3.14% | 3.14% | 3.14%
$SP_{3}$ | 3.08% | 2.88% | 2.88% | 2.89% | 2.89% | 2.89% | 2.89% | 2.89% | 2.89%
$SP_{5}$ | 3.13% | 2.91% | 2.91% | 2.91% | 2.91% | 2.91% | 2.91% | 2.91% | 2.91%
$SP_{7}$ | 3.09% | 2.89% | 2.88% | 2.88% | 2.88% | 2.88% | 2.88% | 2.88% | 2.88%
$SP_{9}$ | 3.06% | 2.86% | 2.86% | 2.86% | 2.86% | 2.86% | 2.86% | 2.86% | 2.86%
$SP_{11}$ | 3.07% | 2.88% | 2.88% | 2.88% | 2.88% | 2.88% | 2.88% | 2.88% | 2.88%
$SP_{13}$ | 3.11% | 2.91% | 2.91% | 2.91% | 2.91% | 2.91% | 2.91% | 2.91% | 2.91%
$SP_{15}$ | 3.09% | 2.88% | 2.88% | 2.88% | 2.88% | 2.88% | 2.88% | 2.88% | 2.88%
$SP_{17}$ | 3.10% | 2.89% | 2.89% | 2.89% | 2.89% | 2.89% | 2.89% | 2.89% | 2.89%
$RTE$ | 3.58% | 3.14% | 3.16% | 3.18% | 3.19% | 3.19% | 3.19% | 3.19% | 3.20%
Table 3: Same as Tab 2, but for source function $f_{2}$. | $SP_{1}$ | $SP_{3}$ | $SP_{5}$ | $SP_{7}$ | $SP_{9}$ | $SP_{11}$ | $SP_{13}$ | $SP_{15}$ | $SP_{17}$
---|---|---|---|---|---|---|---|---|---
$SP_{1}$ | 2.91% | 3.17% | 3.22% | 3.22% | 3.22% | 3.22% | 3.22% | 3.22% | 3.22%
$SP_{3}$ | 3.10% | 2.88% | 2.89% | 2.90% | 2.90% | 2.90% | 2.90% | 2.90% | 2.90%
$SP_{5}$ | 3.14% | 2.89% | 2.89% | 2.89% | 2.89% | 2.89% | 2.89% | 2.89% | 2.89%
$SP_{7}$ | 3.15% | 2.89% | 2.88% | 2.88% | 2.88% | 2.88% | 2.89% | 2.89% | 2.89%
$SP_{9}$ | 3.16% | 2.91% | 2.89% | 2.89% | 2.89% | 2.89% | 2.89% | 2.89% | 2.89%
$SP_{11}$ | 3.14% | 2.89% | 2.87% | 2.87% | 2.87% | 2.87% | 2.87% | 2.87% | 2.87%
$SP_{13}$ | 3.12% | 2.88% | 2.87% | 2.87% | 2.87% | 2.87% | 2.87% | 2.87% | 2.87%
$SP_{15}$ | 3.17% | 2.91% | 2.90% | 2.89% | 2.89% | 2.89% | 2.89% | 2.89% | 2.89%
$SP_{17}$ | 3.16% | 2.90% | 2.89% | 2.89% | 2.89% | 2.89% | 2.89% | 2.89% | 2.89%
$RTE$ | 8.02% | 7.12% | 6.77% | 6.65% | 6.60% | 6.59% | 6.60% | 6.60% | 6.61%

Observe the diagonals of above tables, one can find out that the
reconstruction error is not growing as $N$ grows, this is because the estimate
in Theorem 3.1 is only meant for the worst boundary source. For the given
source function and coefficient set, the reconstruction based on $SP_{1}$
model (diffusion approximation) appears not as good as the other models when
the data are generated from $SP_{N}$ models. While the performances of most
low order $SP_{N}$ models ($N\leq 7$) are already close to the ones of high
order $SP_{N}$ models ($N\geq 9$).

Particularly, when the datum $H$ is generated from the radiative transport
model, all of the reconstruction errors of $SP_{N}$ models become larger due
to the additional modeling errors, see Section 4.2. We plot some of the
reconstructions in Fig 4. It is not surprising to find that the errors are
relatively larger near boundary since the $SP_{N}$ equation system is still
elliptic over the whole domain, while the behavior of radiative transport
averaged solution is hyperbolic in the vicinity of boundary sets.

Figure 4: The error of reconstruction for the absorption coefficient
$\sigma_{a}$ with datum generated from radiative transport model. Top row,
from left to right: reconstruction error for source $f_{1}$ with $N=1$, $N=7$,
$N=17$. Bottom row is the same but for source $f_{2}$.

###  4.4 Reconstruction of $\sigma_{s}$ only

The reconstruction of $\sigma_{s}$ is a nonlinear problem for $N\geq 3$. In
order to provide a fair comparison of across the $SP_{N}$ models, we will use
the optimization based method to reconstruct the scattering coefficient:

$\sigma_{s}=\operatorname*{arg\,min}_{\sigma_{s}}\frac{1}{2}\int_{\Omega}\left(|H(\mathbf{x})-H^{\ast}|^{2}+|\nabla(H(\mathbf{x})-H^{\ast})|^{2}\right)d\mathbf{x}+\frac{\beta}{2}\int_{\Omega}|\nabla\sigma_{s}|^{2}d\mathbf{x}\,,$
(82)

where $H^{\ast}$ is the measured datum and $\beta$ is the regularization
parameter. The optimization problem is solved by L-BFGS method, the gradient
is computed from the adjoint state technique. The choice of regularization
parameter $\beta$ should depend on the noise level. One should notice that we
use $H^{1}(\Omega)$ norm instead of the traditional $L^{2}(\Omega)$
minimization:

$\sigma_{s}=\operatorname*{arg\,min}_{\sigma_{s}}\frac{1}{2}\int_{\Omega}|H(\mathbf{x})-H^{\ast}|^{2}d\mathbf{x}+\frac{\beta}{2}\int_{\Omega}|\nabla\sigma_{s}|^{2}d\mathbf{x}\,.$
(83)

This is due to the stability estimate in Theorem 3.7, where actually the
$H^{2}(\Omega)$ norm is needed for the stability estimate. However, such
regularity requirement implies that $H(\mathbf{x})$ needs to be globally
$C^{1}$ from Sobolev embedding, which means our finite element space needs to
equip with polynomials of five or higher degrees (e.g. Argyris element). Here
we relax objective functional to $H^{1}(\Omega)$ norm simply to avoid the
extraordinary computational cost. To get a brief impression about the two
optimization schemes, we take the $f_{1}$ source function with $SP_{3}$ model
(both data generation and reconstruction) for an example, the datum $H$ is not
contaminated (noise level $\gamma=0$) and regularization parameter $\beta=0$
as well. The reconstructed scattering coefficients are shown in Fig 5. One can
tell from the images that the coefficient recovered from $L^{2}$ optimization
(83) still contains background artifacts. The reason behind is the relatively
strong smoothing effect of the mapping $\sigma_{s}(\mathbf{x})\mapsto
H(\mathbf{x})$, where the high frequency information in $\sigma_{s}$ could not
be fully recovered if we emphasize equally on $H(\mathbf{x})$’s frequency
information.

Figure 5: Left: The reconstructed scattering coefficeint using $L^{2}$
minimization, the $L^{2}$ relative error is $4.25\%$. Right: The reconstructed
scattering coefficient using $H^{1}$ minimization, the $L^{2}$ relative error
is $1.10\%$.

Therefore the data contamination should be treated carefully for the $H^{1}$
minimization (82). This is because if we still apply the random noises by
$H^{\ast}=H(1+\gamma\texttt{random})$ at each mesh node, then
$\|H-H^{\ast}\|_{H^{1}(\Omega)}=\mathcal{O}(\sqrt{N}\gamma)$, where $N$ is the
total number of nodes assuming the mesh is uniform. Therefore instead of
pointwise multiplicative noise, we aggressively contaminate the datum by
perturbing its Fourier modes

$\mathcal{F}{H^{\ast}}(\boldsymbol{\xi})=\mathcal{F}{H}(\boldsymbol{\xi})(1+\gamma\texttt{random})\,,$
(84)

where we have used $\gamma=2\%$ as the noise level parameter. Note that such
noise will perturb the low frequency modes of $H$ which might cause severe
global artifacts in the reconstruction.

Table 4: Relative $L^{2}$ error of the reconstructed $\sigma_{s}$ for source function $f_{1}$ with different generating models and reconstruction models. The row label represents the generating model and the column label represents the reconstruction model. | $SP_{1}$ | $SP_{3}$ | $SP_{5}$ | $SP_{7}$ | $SP_{9}$ | $SP_{11}$ | $SP_{13}$ | $SP_{15}$ | $SP_{17}$
---|---|---|---|---|---|---|---|---|---
$SP_{1}$ | 12.0% | 19.0% | 19.5% | 19.6% | 19.6% | 19.6% | 19.6% | 19.7% | 19.7%
$SP_{3}$ | 13.9% | 12.1% | 12.5% | 12.8% | 13.0% | 13.2% | 13.3% | 13.4% | 13.5%
$SP_{5}$ | 14.9% | 12.5% | 12.3% | 12.5% | 12.7% | 12.9% | 13.0% | 13.1% | 13.2%
$SP_{7}$ | 15.1% | 13.3% | 12.5% | 12.1% | 11.9% | 11.8% | 11.7% | 11.7% | 11.6%
$SP_{9}$ | 15.1% | 13.3% | 12.5% | 12.3% | 12.2% | 12.3% | 12.3% | 12.4% | 12.4%
$SP_{11}$ | 15.5% | 14.4% | 13.3% | 12.7% | 12.3% | 12.0% | 11.8% | 11.7% | 11.6%
$SP_{13}$ | 15.1% | 13.8% | 12.8% | 12.5% | 12.3% | 12.2% | 12.2% | 12.2% | 12.2%
$SP_{15}$ | 15.7% | 15.1% | 14.0% | 13.2% | 12.7% | 12.3% | 12.1% | 12.0% | 11.8%
$SP_{17}$ | 15.2% | 14.2% | 13.1% | 12.7% | 12.4% | 12.3% | 12.2% | 12.2% | 12.2%
$RTE$ | 25.8% | 25.7% | 18.3% | 21.5% | 36.9% | 37.2% | 39.2% | 39.3% | 39.4%
Table 5: Same as Tab 4 but for source function $f_{2}$. | $SP_{1}$ | $SP_{3}$ | $SP_{5}$ | $SP_{7}$ | $SP_{9}$ | $SP_{11}$ | $SP_{13}$ | $SP_{15}$ | $SP_{17}$
---|---|---|---|---|---|---|---|---|---
$SP_{1}$ | 20.8% | 21.4% | 21.3% | 23.1% | 21.4% | 21.4% | 21.5% | 21.6% | 21.7%
$SP_{3}$ | 23.1% | 18.0% | 18.6% | 19.1% | 19.7% | 19.9% | 20.2% | 20.3% | 20.5%
$SP_{5}$ | 26.9% | 18.3% | 17.7% | 18.0% | 18.3% | 18.5% | 18.6% | 18.6% | 18.7%
$SP_{7}$ | 28.0% | 22.7% | 18.2% | 16.4% | 16.3% | 16.4% | 16.6% | 16.8% | 16.9%
$SP_{9}$ | 29.0% | 19.9% | 18.0% | 17.8% | 17.7% | 17.8% | 17.9% | 18.06% | 18.1%
$SP_{11}$ | 29.4% | 19.5% | 16.9% | 16.4% | 16.3% | 16.3% | 16.4% | 16.4% | 16.5%
$SP_{13}$ | 29.5% | 21.0% | 18.4% | 18.0% | 17.8% | 17.8% | 17.8% | 17.8% | 17.9%
$SP_{15}$ | 30.2% | 20.6% | 17.3% | 16.6% | 16.4% | 16.3% | 16.3% | 16.3% | 16.4%
$SP_{17}$ | 30.5% | 21.5% | 19.4% | 18.1% | 17.9% | 17.8% | 17.7% | 17.8% | 17.8%
$RTE$ | 55.9% | 47.8% | 49.4% | 50.2% | 50.3% | 50.3% | 51.4% | 51.5% | 51.6%

Similar to the previous numerical experiment, we summarize the result in the
Tab 4, where we have fixed the regularization parameter $\beta=10^{-8}$ for
these experiments. From the table, we could clearly see that $SP_{1}$ is not
as good as other models for the reconstruction of $\sigma_{s}$ when the datum
$H$ is coming from higher order models.

The reconstruction errors on the diagonal of the tables look converging as the
order $N$ grows, which indicates that the models converge relatively fast for
the given source functions and the coefficient setting. When the datum $H$ is
generated from the radiative transport equation, the reconstruction error
becomes larger. The reconstructions with respect to source function $f_{2}$
are significantly worse than the ones for $f_{1}$, see Fig 6. However, this
could be explained through the analogue with the $SP_{1}$ model, where the
scattering coefficient’s reconstruction is to solve a transport equation [5]:

$\nabla\phi_{0}\cdot\nabla D(x)+D(x)\Delta\phi_{0}-\sigma_{a}\phi_{0}=0,$ (85)

where $D(x)=1/(3(1-g)\sigma_{s})$, where $\phi_{0}$ and
$\sigma_{s}|_{\partial\Omega}$ are known. Therefore when $\nabla\phi_{0}\neq
0$, the function $D(x)$ could be solved by tracing the characteristics. If
$\nabla\phi_{0}$ vanishes or appears to be small, then the characteristics
could be trapped, where the reconstructions are based on regularization only.

Figure 6: The reconstructions for the scattering coefficient $\sigma_{s}$ with
datum generated from radiative transport model. Top row, from left to right:
reconstruction error for source $f_{1}$ with $N=1$, $N=7$, $N=17$. Bottom row
is the same but for source $f_{2}$.

###  4.5 Reconstruction of $\sigma_{a}$ and $\Upsilon$

In this section, we consider the non-linearized case and follow the
aforementioned two-step reconstruction strategy. Suppose $H_{1},H_{2}$ are the
data sets measured with boundary source functions $f_{1}$ and $f_{2}$,
respectively. Our numerical reconstruction solves the optimization problem:

$\sigma_{a}^{\ast}=\operatorname*{arg\,min}_{\sigma_{a}}\Big{\|}\frac{H_{1}}{H_{2}}-\frac{\mathbf{s}_{1}\cdot\Phi_{1}}{\mathbf{s}_{1}\cdot\Phi_{2}}\Big{\|}_{L^{2}(\Omega)}^{2}+\alpha\|\nabla\sigma_{a}\|_{L^{2}(\Omega)}^{2}\,,$
(86)

where $\Phi_{1},\Phi_{2}$ are the solutions to the $SP_{N}$ equation with the
absorption coefficient $\sigma_{a}$ and source function $f_{1},f_{2}$,
respectively. Intuitively, the ratio ${H_{1}}/{H_{2}}$ should be quite smooth
and weakly depends on $\sigma_{a}$, which means the reconstruction for
$\sigma_{a}$ could be very unstable. In the following, we assume the data sets
$H_{1}$ and $H_{2}$ are generated from the $SP_{N}$ models with multiplicative
noise $H_{i}^{\ast}=H_{i}(1+\gamma\texttt{random})$ for $\gamma=0.1\%$ only,
the regularization parameter is fixed as $\alpha=10^{-8}$. Then we reconstruct
the absorption coefficient using the same model. The reconstructions are shown
in Fig 7, it could be seen that the reconstructions are very unstable even for
small noise, only limited resolution could be obtained.

Figure 7: The reconstructions of absorption coefficients $\sigma_{a}$ with
respect to $SP_{N}$ models. From left to right: $N=1,7,17$. The relative
$L^{2}$ errors are $16.6\%,16.2\%,15.9\%$, respectively.

###  4.6 Reconstruction of $\sigma_{a}$ and $\sigma_{s}$

Similar to the previous case, we consider the reconstruction in a two-step
process as well. Suppose $H_{1},H_{2}$ are the data sets with boundary sources
$f_{1},f_{2}$, respectively. Our algorithm will first construct the scattering
coefficient from the following optimization problem:

$\sigma_{s}^{\ast}=\operatorname*{arg\,min}_{\sigma_{s}}\Big{\|}\frac{H_{1}}{H_{2}}-\frac{\mathbf{s}_{1}\cdot\Phi_{1}}{\mathbf{s}_{1}\cdot\Phi_{2}}\Big{\|}_{H^{1}(\Omega)}^{2}+\alpha\|\nabla\sigma_{s}\|_{L^{2}(\Omega)}^{2}\,,$
(87)

where $\Phi_{1},\Phi_{2}$ are the solutions to the modified $SP_{N}$ equation
(20), where $\sigma_{a}$ has been replaced by
$H_{i}/(\Upsilon\mathbf{s}_{1}\cdot\Phi_{i})$, $i=1,2$. Here we have taken the
$H_{1}$ minimization. In the following numerical experiment, we assume
$H_{1},H_{2}$ are generated from the $SP_{N}$ model, the data sets are
contaminated on the Fourier space through

$\mathcal{F}H^{\ast}_{i}(\boldsymbol{\xi})=\mathcal{F}\mathcal{H}_{i}(\boldsymbol{\xi})(1+\gamma\texttt{random})$

with $\gamma=2\%$. We also fix the regularization parameter
$\alpha=10^{-8}$.The numerical reconstructions are performed over the same
$SP_{N}$ model and the results are illustrated in Fig 8. After the scattering
coefficient has been reconstructed, we will use the recovered scattering
coefficient to find the absorption coefficient following the Experiment 4.4.
The corresponding reconstruction errors of the absorption coefficients are
shown in Fig 9. It can be seen that even the reconstruction of scattering
coefficients contain background artifacts, while the reconstruction errors of
$\sigma_{a}$ are still quite small.

Figure 8: The reconstructions of scattering coefficients $\sigma_{s}$ with
respect to $SP_{N}$ models. From left to right: $N=1,7,17$. The relative
$L^{2}$ errors are $16.1\%,16.7\%,16.8\%$, respectively.

Figure 9: The reconstruction error of absorption coefficients $\sigma_{a}$
with respect to $SP_{N}$ models by using the recovered scattering coefficients
from Fig 8. From left to right: $N=1,7,17$. The relative $L^{2}$ errors are
$1.76\%,5.79\%,5.80\%$, respectively.

##  5 Conclusion

In this work, we studied the quantitative photoacoustic tomography with the
simplified $P_{N}$ approximation model to the radiative transport equation. We
have derived the uniqueness and stability estimates for the reconstruction of
one single coefficient of $(\sigma_{a}$, $\sigma_{s}$, $\Upsilon$) from one
initial pressure datum $H(\mathbf{x})$. For the simultaneous reconstruction of
two coefficients, we have considered the linearized setting and introduced the
optimization based numerical algorithm for the reconstruction. We showed the
numerical simulations based on a synthetic data to validate the mathematical
analysis.

## Acknowledgment

H. Zhao’s research is partially supported by NSF DMS-2048877 and DMS-2012860.

##  Appendix A Appendix

######  Lemma A.1.

Let $s_{k,l}$ be the $(k,l)$-th entry of $M^{-1}$, then

$s_{k,l}=\begin{cases}\frac{1}{2k-1}(-1)^{l-k}\frac{((2l-2))!!}{(2l-1)!!}\frac{(2k-1)!!}{(2k-2)!!},\quad&l\geq
k\,,\\\ 0,&\text{otherwise}\,.\end{cases}$ (88)

###### Proof.

Let $S$ be the matrix with $(k,l)$-th entry as $s_{k,l}$, we then compute the
$(k,j)$-th entry of $SM$ by $a_{k,j}:=s_{k,j}(2j-1)+s_{k,j-1}(2j-2)$. If
$k=j$, $s_{k,j}=\frac{1}{2j-1}$, therefore $a_{k,j}=\frac{1}{2j-1}(2j-1)=1$.
If $j<k$, notice that now $s_{k,j}=0$, we must have $a_{k,j}=0$. If $j>k$, we
compute $a_{k,j}$ directly

$\frac{1}{2k-1}\frac{(2k-1)!!}{(2k-2)!!}\left((-1)^{j-k}\frac{((2j-2))!!}{(2j-1)!!}(2j-1)+(-1)^{j-k-1}\frac{((2j-4))!!}{(2j-3)!!}(2j-2)\right)=0\,.$
(89)

Hence $S=M^{-1}$. ∎

######  Lemma A.2.

$\det(R^{-1})\leq\mathcal{O}(N^{3/8})$, hence
$\det(R)\geq\mathcal{O}(N^{-3/8})$.

###### Proof.

Since

$\displaystyle R_{ij}$
$\displaystyle=(-1)^{i+j-1}\frac{2\Gamma(i+\frac{1}{2})\Gamma(j-\frac{1}{2})}{\pi\Gamma(i)\Gamma(j)}\frac{(4j-3)}{(2i+2j-2)(2j-2i-1)}\,,$
(90)
$\displaystyle=(-1)^{i+j}\frac{\Gamma(i+\frac{1}{2})\Gamma(j-\frac{1}{2})}{\pi\Gamma(i)\Gamma(j)}\frac{2j-\frac{3}{2}}{(i-\frac{1}{4})^{2}-(j-\frac{3}{4})^{2}}\,,$

which gives the factorization

$R=UGV\,,$ (91)

where diagonal matrices

$\displaystyle U$ $\displaystyle={\rm diag}(u_{1},\dots,u_{(N+1)/2})\,,$ (92)
$\displaystyle V$ $\displaystyle={\rm diag}(v_{1},\dots,v_{(N+1)/2})\,,$

with
$u_{i}=(-1)^{i}\sqrt{\frac{2}{\pi}}\frac{\Gamma(i+\frac{1}{2})}{\Gamma(i)},v_{j}=(-1)^{j}\sqrt{\frac{2}{\pi}}\frac{\Gamma(j-\frac{1}{2})}{\Gamma(j)}(j-\frac{3}{4})\,,$
and $G$ is the Cauchy-Toeplitz matrix,

$G=\left(\frac{1}{(i-\frac{1}{4})^{2}-(j-\frac{3}{4})^{2}}\right)_{i,j=1,\dots,(N+1)/2}\,.$
(93)

We denote $x_{i}=(i-\frac{1}{4})^{2}$ and $y_{j}=(j-\frac{3}{4})^{2}$, then
use the Theorem 2.1 in [29],

$G^{-1}=PG^{T}Q\,,$ (94)

where ${P}=\text{diag}(p_{1},\dots,p_{(N+1)/2})$ and
$Q=\text{diag}(q_{1},\dots,q_{(N+1)/2})$ satisfy

$\displaystyle G[p_{1},\dots,p_{(N+1)/2}]^{T}$
$\displaystyle=[1,\dots,1]^{T},$ (95) $\displaystyle
G^{T}[q_{1},\dots,q_{(N+1)/2}]^{T}$ $\displaystyle=[1,\dots,1]^{T}\,.$

Also $p_{i},q_{j}$ are computed explicitly by Cramer’s law,

$\displaystyle p_{i}$
$\displaystyle=\prod_{l=1}^{(N+1)/2}(x_{l}-y_{i})/\prod_{l=1,l\neq
i}^{(N+1)/2}(y_{l}-y_{i})\,,$ (96) $\displaystyle q_{j}$
$\displaystyle=\prod_{l=1}^{(N+1)/2}(x_{j}-y_{l})/\prod_{l=1,l\neq
j}^{(N+1)/2}(x_{j}-x_{l})\,.$

By replacing $x_{i},y_{j}$ with their values, we obtain

$\displaystyle p_{i}$
$\displaystyle=\prod_{l=1}^{(N+1)/2}(l+i-1)(l-i+\frac{1}{2})/\prod_{l=1,l\neq
i}^{(N+1)/2}(l-i)(l+i-3/2)$ (97)
$\displaystyle=\left(\prod_{l=1}^{(N+1)/2}\frac{l+i-1}{l+i-\frac{3}{2}}\right)\left(\frac{1}{2}\prod_{l=1}^{i-1}\left(1-\frac{1}{2l}\right)\prod_{l=1}^{(N+1)/2-i}\left(1+\frac{1}{2l}\right)\right)$
$\displaystyle=\frac{\Gamma(i+\frac{N+1}{2})}{\Gamma(i+\frac{N}{2})}\frac{\Gamma(i-\frac{1}{2})}{\Gamma(i)}\left(2i-\frac{3}{2}\right)\left(\frac{1}{2}\prod_{l=1}^{i-1}\left(1-\frac{1}{2l}\right)\prod_{l=1}^{(N+1)/2-i}\left(1+\frac{1}{2l}\right)\right)\,,$
$\displaystyle q_{j}$
$\displaystyle=\prod_{l=1}^{(N+1)/2}(j+l-1)(j-l+\frac{1}{2})/\prod_{l=1,l\neq
j}^{(N+1)/2}(j-l)(j+l-1/2)$
$\displaystyle=\left(\prod_{l=1}^{(N+1)/2}\frac{j+l-1}{j+l-\frac{1}{2}}\right)\left(\frac{1}{2}\prod_{l=1}^{j-1}\left(1+\frac{1}{2l}\right)\prod_{l=1}^{(N+1)/2-j}\left(1-\frac{1}{2l}\right)\right)$
$\displaystyle=\frac{\Gamma(j+\frac{N+1}{2})\Gamma(j+\frac{1}{2})}{\Gamma(j+\frac{N}{2}+1)\Gamma(j)}\left(2j-\frac{1}{2}\right)\left(\frac{1}{2}\prod_{l=1}^{j-1}\left(1+\frac{1}{2l}\right)\prod_{l=1}^{(N+1)/2-j}\left(1-\frac{1}{2l}\right)\right)\,.$

Then we can easily deduce
$\det(R^{-1})=\sqrt{\det(P)\det(Q)}/(\det(V)\det(U))$, since all matrices
involved are diagonal, let
$S_{i}=\frac{\pi}{2}\prod_{l=1}^{i-1}\left(1+\frac{1}{2l}\right)\prod_{l=1}^{(N+1)/2-i}\left(1-\frac{1}{2l}\right)$,
then from the theory of Gamma functions, we know

$\displaystyle
S_{i}=\frac{\Gamma(i+\frac{1}{2})\Gamma(n-i+\frac{1}{2})}{\Gamma(i)\Gamma(\frac{(N+1)}{2}-i+1)}\,.$
(98)

Hence we can estimate $\det(R^{-1})$’s upper bound by estimating

$\displaystyle\det(R^{-1})$
$\displaystyle=\prod_{i=1}^{(N+1)/2}\frac{\sqrt{p_{i}q_{i}}}{u_{i}v_{i}}$ (99)
$\displaystyle=\prod_{i=1}^{(N+1)/2}S_{i}\sqrt{\frac{\Gamma(i+\frac{N+1}{2})^{2}\Gamma(i)^{2}(i-\frac{1}{4})}{\Gamma(i+\frac{N}{2})\Gamma(i+\frac{N}{2}+1)\Gamma(i+\frac{1}{2})\Gamma(i-\frac{1}{2})(i-\frac{3}{4})}}$
$\displaystyle=\prod_{i=1}^{(N+1)/2}\frac{\Gamma(i-\frac{1}{2})}{\Gamma(i)}\sqrt{\frac{\Gamma(i+\frac{(N+1)}{2})^{2}(i-\frac{1}{2})(i-\frac{1}{4})}{\Gamma(i+\frac{(N+1)}{2}-\frac{1}{2})\Gamma(i+\frac{(N+1)}{2}+\frac{1}{2})(i-\frac{3}{4})}}\,.$

Then by noticing the Chu’s Double Inequality [7],

$\sqrt{x-\frac{1}{4}}<\frac{\Gamma(x+\frac{1}{2})}{\Gamma(x)}<\frac{x}{\sqrt{x+\frac{1}{4}}}\,,$
(100)

the following estimates hold,

$\displaystyle\sqrt{\frac{\Gamma(i+\frac{N+1}{2})^{2}}{\Gamma(i+\frac{N}{2})\Gamma(i+\frac{N}{2}+1)}}$
$\displaystyle\leq\sqrt{\frac{i+\frac{(N+1)}{2}-\frac{1}{2}}{i+\frac{(N+1)}{2}-\frac{1}{4}}}<1\,,$
(101)
$\displaystyle\sqrt{\frac{\Gamma(i-\frac{1}{2})^{2}}{\Gamma(i)^{2}}\frac{(i-\frac{1}{2})(i-\frac{1}{4})}{i-\frac{3}{4}}}$
$\displaystyle\leq\sqrt{\frac{(i-\frac{1}{2})(i-\frac{1}{4})}{(i-\frac{3}{4})^{2}}}\,,$

which implies

$\displaystyle\det(R^{-1})$
$\displaystyle<\prod_{i=1}^{(N+1)/2}\sqrt{\frac{(4i-2)(4i-1)}{(4i-3)^{2}}}$
(102)
$\displaystyle=\exp\left(\frac{1}{2}\sum_{i=1}^{(N+1)/2}\log\left(1+\frac{3}{4i-3}+\frac{2}{(4i-3)^{2}}\right)\right)$
$\displaystyle\leq\exp\left(\frac{1}{2}\sum_{i=1}^{(N+1)/2}\frac{3}{4i-3}+\frac{2}{(4i-3)^{2}}\right)$
$\displaystyle=\mathcal{O}\left(\exp\left(\frac{3}{8}\log\left(\frac{(N+1)}{2}\right)\right)\right)=\mathcal{O}(N^{3/8})\,.$

∎

######  Lemma A.3.

$\|R\|_{F}^{2}\leq\frac{N+1}{2}$.

###### Proof.

We show that for all $m\geq 1$,

$\sum_{j=1}^{m}|R_{j,m}|^{2}+\sum_{j=1}^{m-1}|R_{m,j}|^{2}\leq 1\,.$ (103)

The above estimate is true for $m=1$ since $R_{1,1}=\frac{1}{2}$, we only
focus on the cases that $m\geq 2$. Reformulate $R_{j,m}$ by the Gamma function
as

$|R_{j,m}|=\frac{\Gamma(j+\frac{1}{2})\Gamma(m-\frac{1}{2})}{\pi\Gamma(j)\Gamma(m)(j+m-1)}\frac{(4m-3)}{2m-2j-1}\,.$
(104)

From the Chu’s Double Inequality [7]

$\sqrt{x-\frac{1}{4}}<\frac{\Gamma(x+\frac{1}{2})}{\Gamma(x)}<\frac{x}{\sqrt{x+\frac{1}{4}}}\,,$
(105)

the following estimates hold,

$\frac{\Gamma(j+\frac{1}{2})}{\Gamma(j)}<\frac{j}{\sqrt{j+\frac{1}{4}}}\,,\quad\frac{\Gamma(m-\frac{1}{2})}{\Gamma(m)}<\frac{1}{(m-\frac{3}{4})^{1/2}}\,,$
(106)

we can deduce the estimate

$\displaystyle\sum_{j=1}^{m}|R_{j,m}|^{2}$
$\displaystyle=\frac{1}{\pi^{2}}\sum_{j=1}^{m}\left(\frac{\Gamma(j+\frac{1}{2})\Gamma(m-\frac{1}{2})}{\Gamma(j)\Gamma(m)}\right)^{2}\frac{(4m-3)^{2}}{(j+m-1)^{2}}\frac{1}{(2m-2j-1)^{2}}$
(107)
$\displaystyle\leq\frac{1}{\pi^{2}}\sum_{j=1}^{m}\frac{(4m-3)^{2}(j-\frac{1}{4}+\frac{1}{4(1+4j)})}{(m-\frac{3}{4})(j+m-1)^{2}(2m-2j-1)^{2}}$
$\displaystyle\leq\frac{4}{\pi^{2}}\sum_{j=1}^{m}\frac{(4j-1+\frac{1}{(1+4j)})(4m-3)}{(2j+2m-2)^{2}(2m-2j-1)^{2}}$
$\displaystyle=\frac{4}{\pi^{2}}\sum_{j=1}^{m}\left[\frac{1}{(2m-2j-1)^{2}}-\frac{1}{(2j+2m-2)^{2}}\right]$
$\displaystyle\quad+\frac{4}{\pi^{2}}\sum_{j=1}^{m}\frac{(4m-3)}{(1+4j)(2j+2m-2)^{2}(2m-2j-1)^{2}}$
$\displaystyle\leq\frac{4}{\pi^{2}}\left[\left(1+\frac{\pi^{2}}{8}-\frac{1}{4m}\right)+\frac{1}{5m(2m-3)^{2}}+\frac{1}{9m}\sum_{j=2}^{m}\frac{1}{(2m-2j-1)^{2}}\right]$
$\displaystyle\leq\frac{4}{\pi^{2}}\left[\left(1+\frac{\pi^{2}}{8}\right)+\left(\frac{11}{180}+\frac{\pi^{2}}{72}\right)\frac{1}{m}\right]\,.$

The other part can be estimated in a similar way,

$\displaystyle\sum_{j=1}^{m-1}|R_{m,j}|^{2}$
$\displaystyle=\frac{1}{\pi^{2}}\sum_{j=1}^{m-1}\left(\frac{\Gamma(m+\frac{1}{2})\Gamma(j-\frac{1}{2})}{\Gamma(m)\Gamma(j)}\right)^{2}\frac{(4j-3)^{2}}{(j+m-1)^{2}}\frac{1}{(2j-2m-1)^{2}}$
(108)
$\displaystyle\leq\frac{1}{\pi^{2}}\sum_{j=1}^{m-1}\frac{(4j-3)^{2}(m-\frac{1}{4}+\frac{1}{4(1+4m)})}{(j-\frac{3}{4})(j+m-1)^{2}(2j-2m-1)^{2}}$
$\displaystyle\leq\frac{4}{\pi^{2}}\sum_{j=1}^{m-1}\frac{(4m-1+\frac{1}{1+4m})(4j-3)}{(2j+2m-2)^{2}(2j-2m-1)^{2}}$
$\displaystyle=\frac{4}{\pi^{2}}\sum_{j=1}^{m-1}\left[\frac{1}{(2m+1-2j)^{2}}-\frac{1}{(2j+2m-2)^{2}}\right]$
$\displaystyle\quad+\frac{4}{\pi^{2}}\sum_{j=1}^{m-1}\frac{4j-3}{(1+4m)(2j+2m-2)^{2}(2m+1-2j)^{2}}$
$\displaystyle\leq\frac{4}{\pi^{2}}\left[\left(\frac{\pi^{2}}{8}-1-\frac{1}{4m}\right)+\sum_{j=1}^{m-1}\frac{1}{4m^{2}(2m+1-2j)^{2}}\right]$
$\displaystyle\leq\frac{4}{\pi^{2}}\left[\left(\frac{\pi^{2}}{8}-1-\frac{1}{4m}\right)+\left(\frac{\pi^{2}}{8}-1\right)\frac{1}{4m^{2}}\right]\,.$

Hence we can estimate that

$\sum_{j=1}^{m}|R_{j,m}|^{2}+\sum_{j=1}^{m-1}|R_{m,j}|^{2}\leq
1-\frac{4}{\pi^{2}}\left(\frac{17}{90}-\frac{\pi^{2}}{72}\right)\frac{1}{m}+\frac{1}{\pi^{2}}\left(\frac{\pi^{2}}{8}-1\right)\frac{1}{m^{2}}\,,$
(109)

the above estimate is strictly less than $1$ for $m\geq 2$, then the Frobenius
norm’s square of $R$ is estimated by

$\|R\|_{F}^{2}=\frac{1}{2}+\sum_{m=2}^{(N+1)/2}\left(\sum_{j=1}^{m}|R_{j,m}|^{2}+\sum_{j=1}^{m-1}|R_{m,j}|^{2}\right)\leq\frac{N}{2}\,.$
(110)

∎

######  Lemma A.4.

Let $\mathbf{s}_{k}$ the $k$-th row of the matrix $T^{-1}$ in (7), and matrix
$\mathbf{Q}$ is defined in (43), then

$\displaystyle\|\mathbf{s}_{1}\cdot\mathbf{Q}\|_{\ell^{p}}=\mathcal{O}(N^{3/2+1/p})\,.$
(111)

###### Proof.

Let $s_{k,j}$ denote the $j$-th entry of vector $\mathbf{s}_{k}$. Combine the
Lemma A.1 and Gautschi’s inequality [12], we have

$|s_{k,j}|<\frac{1}{2k-1}\frac{\sqrt{k+\frac{1}{2}}}{\sqrt{j-\frac{1}{2}}}\,.$
(112)

For the matrix $\mathbf{Q}$, we denote its $(n,j)$-th entry by $Q_{nj}$, which
can be estimated by

$\displaystyle|Q_{nj}|$ $\displaystyle\leq(4n-1)\sum_{k\geq
2}^{\min(j,n)}(4k-3)|s_{k,n}s_{k,j}|\leq(4n-1)\sum_{k\geq
2}^{\min(j,n)}\frac{(4k-3)}{(2k-1)^{2}}\frac{k+\frac{1}{2}}{\sqrt{j-\frac{1}{2}}\sqrt{n-\frac{1}{2}}}$
(113)
$\displaystyle\leq\frac{(4n-1)}{\sqrt{j-\frac{1}{2}}\sqrt{n-\frac{1}{2}}}\sum_{k\geq
2}^{\min(j,n)}\frac{(4k-3)(k+\frac{1}{2})}{(2k-1)^{2}}\,.$

Since $(4k-3)(k+\frac{1}{2})\leq 2(2k-1)^{2}$ for all $k$, we will have
$|Q_{n,j}|\leq\frac{(8n-2)}{\sqrt{n-\frac{1}{2}}\sqrt{j-\frac{1}{2}}}(\min(j,n)-1)$,
therefore, by simple calculations,
$\|\mathbf{s}_{1}\cdot\mathbf{Q}\|_{\ell^{\infty}}$ is bounded by

$\displaystyle\|\mathbf{s}_{1}\cdot\mathbf{Q}\|_{\ell^{\infty}}$
$\displaystyle\leq\sup_{j\geq 1}\sum_{n\geq
1}|s_{1,n}||Q_{n,j}|\leq\sup_{j\geq 1}\sum_{n\geq
1}\frac{(8n-2)(\min(j,n)-1)}{(n-\frac{1}{2})\sqrt{j-\frac{1}{2}}}$ (114)
$\displaystyle=\mathcal{O}(N^{3/2})\,.$

Then the $\ell^{p}$ estimate is

$\|\mathbf{s}_{1}\cdot\mathbf{Q}\|_{\ell^{p}}\leq\|\mathbf{s}_{1}\cdot\mathbf{Q}\|_{\ell^{\infty}}\left(\frac{(N+1)}{2}\right)^{1/p}=\mathcal{O}(N^{3/2+1/p})\,.$

∎

## References

  * [1] M. Agranovsky, P. Kuchment, and L. Kunyansky, On reconstruction formulas and algorithms for the thermoacoustic and photoacoustic tomography, Preprint, (2007).
  * [2] H. Ammari, E. Bretin, V. Jugnon, and A. Wahab, Photoacoustic imaging for attenuating acoustic media, in Mathematical modeling in biomedical imaging II, Springer, 2012, pp. 57–84.
  * [3] G. Bal, Hybrid inverse problems and internal functionals, Inverse problems and applications: inside out. II, 60 (2013), pp. 325–368.
  * [4] G. Bal, A. Jollivet, and V. Jugnon, Inverse transport theory of photoacoustics, Inverse Problems, 26 (2010), p. 025011.
  * [5] G. Bal and K. Ren, Multi-source quantitative photoacoustic tomography in a diffusive regime, Inverse Problems, 27 (2011), p. 075003.
  * [6] P. Beard, Biomedical photoacoustic imaging, Interface focus, 1 (2011), pp. 602–631.
  * [7] J. Chu, A modified wallis product and some applications, The American Mathematical Monthly, 69 (1962), pp. 402–404.
  * [8] M. Chu, K. Vishwanath, A. D. Klose, and H. Dehghani, Light transport in biological tissue using three-dimensional frequency-domain simplified spherical harmonics equations, Physics in Medicine & Biology, 54 (2009), p. 2493.
  * [9] B. Cox, J. Laufer, and P. Beard, The challenges for quantitative photoacoustic imaging, in Photons Plus Ultrasound: Imaging and Sensing 2009, vol. 7177, International Society for Optics and Photonics, 2009, p. 717713.
  * [10] T. Ding, K. Ren, and S. Vallélian, A one-step reconstruction algorithm for quantitative photoacoustic imaging, Inverse Problems, 31 (2015), p. 095005.
  * [11] C. Frederick, K. Ren, and S. Vallélian, Image reconstruction in quantitative photoacoustic tomography with the simplified p_2 approximation, SIAM Journal on Imaging Sciences, 11 (2018), pp. 2847–2876.
  * [12] W. Gautschi, Some elementary inequalities relating to the gamma and incomplete gamma function, J. Math. Phys, 38 (1959), pp. 77–81.
  * [13] E. M. Gelbard, Application of spherical harmonics method to reactor problems, Bettis Atomic Power Laboratory, West Mifflin, PA, Technical Report No. WAPD-BT-20, (1960).
  * [14] A. D. Güngör, Erratum to “an upper bound for the condition number of a matrix in spectral norm”[j. comput. appl. math. 143 (2002) 141–144], Journal of Computational and Applied Mathematics, 234 (2010), p. 316.
  * [15] M. Haltmeier, T. Schuster, and O. Scherzer, Filtered backprojection for thermoacoustic computed tomography in spherical geometry, Mathematical methods in the applied sciences, 28 (2005), pp. 1919–1937.
  * [16] Y. Hristova, Time reversal in thermoacoustic tomography—an error estimate, Inverse Problems, 25 (2009), p. 055008.
  * [17] A. D. Klose and E. W. Larsen, Light transport in biological tissue based on the simplified spherical harmonics equations, Journal of Computational Physics, 220 (2006), pp. 441–470.
  * [18] , Simplified spherical harmonics methods for modeling light transport in biological tissue, in Biomedical Topical Meeting, Optical Society of America, 2006, p. MH3.
  * [19] P. Kuchment and L. Kunyansky, Mathematics of thermoacoustic and photoacoustic tomography, preprint, (2007).
  * [20] A. V. Mamonov and K. Ren, Quantitative photoacoustic imaging in radiative transport regime, arXiv preprint arXiv:1207.4664, (2012).
  * [21] A. W. Marshall, I. Olkin, and B. C. Arnold, Inequalities: theory of majorization and its applications, vol. 143, Springer, 1979.
  * [22] R. G. McClarren, Theoretical aspects of the simplified pn equations, Transport Theory and Statistical Physics, 39 (2010), pp. 73–109.
  * [23] W. McLean and W. C. H. McLean, Strongly elliptic systems and boundary integral equations, Cambridge university press, 2000.
  * [24] G. Piazza and T. Politi, An upper bound for the condition number of a matrix in spectral norm, Journal of Computational and Applied Mathematics, 143 (2002), pp. 141–144.
  * [25] K. Ren, H. Gao, and H. Zhao, A hybrid reconstruction method for quantitative pat, SIAM Journal on Imaging Sciences, 6 (2013), pp. 32–55.
  * [26] O. Scherzer, Handbook of mathematical methods in imaging, Springer Science & Business Media, 2010.
  * [27] P. Stefanov and G. Uhlmann, Thermoacoustic tomography with variable sound speed, Inverse Problems, 25 (2009), p. 075011.
  * [28] E. E. Tyrtyshnikov, Cauchy-toeplitz matrices and some applications, Linear Algebra and its Applications, 149 (1991), pp. 1–18.
  * [29] , Singular values of cauchy-toeplitz matrices, Linear algebra and its applications, 161 (1992), pp. 99–116.
  * [30] L. V. Wang, Photoacoustic imaging and spectroscopy, CRC press, 2017\.
  * [31] S. Wright, M. Schweiger, and S. Arridge, Reconstruction in optical tomography using the pn approximations, Measurement Science and Technology, 18 (2006), p. 79.
  * [32] W. P. Ziemer, Weakly differentiable functions: Sobolev spaces and functions of bounded variation, vol. 120, Springer Science & Business Media, 2012.