Iterative Methods for Photoacoustic Tomography in Attenuating Acoustic Media

TL;DR

This paper develops and analyzes iterative reconstruction methods for photoacoustic tomography in media with frequency-dependent attenuation, providing stable solutions and demonstrating high accuracy through numerical experiments.

Contribution

It introduces a new iterative regularization framework for photoacoustic tomography in attenuating media with explicit convergence analysis and efficient adjoint computations.

Findings

Stable solutions for attenuating media demonstrated

High reconstruction accuracy shown in numerical tests

Effective in limited view scenarios

Abstract

The development of efficient and accurate reconstruction methods is an important aspect of tomographic imaging. In this article, we address this issue for photoacoustic tomography. To this aim, we use models for acoustic wave propagation accounting for frequency dependent attenuation according to a wide class of attenuation laws that may include memory. We formulate the inverse problem of photoacoustic tomography in attenuating medium as an ill-posed operator equation in a Hilbert space framework that is tackled by iterative regularization methods. Our approach comes with a clear convergence analysis. For that purpose we derive explicit expressions for the adjoint problem that can efficiently be implemented. In contrast to time reversal, the employed adjoint wave equation is again damping and, thus has a stable solution. This stability property can be clearly seen in our numerical results. Moreover, the presented numerical results clearly demonstrate the Efficiency and accuracy of the derived iterative reconstruction algorithms in various situations including the limited view case.

Figures

Equations

$${\left({\mathbf{D}}_{\alpha }+\frac{1}{c_0}\frac{\partial }{\partial t}\right)}^2{p}_{\alpha }(x,t)-\Delta p_{\alpha }(x,t)={\delta }^{\prime }(t)h(x)\text{ for }(x,t)\in {\mathbb{R}}^d\times \mathbb{R}.$$

$$\frac{1}{c_0^2}\frac{{\partial }^2}{\partial t^2}{p}_0(x,t)-\Delta p_0(x,t)={\delta }^{\prime }(t)h(x)\text{ for }(x,t)\in {\mathbb{R}}^d\times \mathbb{R}$$

$$\{\begin{array}{rlrl}& {\left({\mathbf{D}}_{\alpha }+\frac{1}{c_0}\frac{\partial }{\partial t}\right)}^2{p}_{\alpha }(\mathbf{x},t)-\Delta p_{\alpha }(\mathbf{x},t)=s(\mathbf{x},t)& & \text{ for }(\mathbf{x},t)\in {\mathbb{R}}^{d+1},\\ & p_{\alpha }(\cdot ,t)=0& & \text{ for }t<0.\end{array}$$

$$G_{\alpha }(\mathbf{x},t)=\frac{K_{\alpha }\left(\mathbf{x},t-\frac{|\mathbf{x}|}{c_0}\right)}{4\pi \left|\mathbf{x}\right|}\text{ with }{K}_{\alpha }(\mathbf{x},t):={\mathcal{F}}_t^{-1}\left(e^{-|\mathbf{x}|\alpha }\right)(t).$$

$$\forall \omega \in \mathbb{R}:{\mathcal{F}}_t(m_{\alpha }(\cdot ,r))(\omega ):=\frac{\omega }{\omega /c_0+\mathsf{i}\alpha (\omega )}{e}^{\mathsf{i}(\omega /c_0+\mathsf{i}\alpha (\omega ))|r|}.$$

$$\forall (x,t)\in {\mathbb{R}}^d\times (0,\infty ):p_{\alpha }(\mathbf{x},t)={\int }_0^t{m}_{\alpha }(t,r)p_0(\mathbf{x},r)\mathrm{d}r,$$

$$\mathrm{Re}[\alpha (\omega )]$$

$$\mathrm{Im}[\alpha (\omega )]$$

$$\mathrm{Re}[\alpha (\omega )]$$

$$\mathrm{Im}[\alpha (\omega )]$$

$$\alpha :\mathbb{R}\to \mathbb{C}:\omega \mapsto a_0(-\mathsf{i}\omega {)}^{\gamma }+b_0(-\mathsf{i}\omega ).$$

$$\alpha (\omega )=\frac{a_0(-\mathsf{i}\omega )}{c_{\infty }\sqrt{1+(-\mathsf{i}{\tau }_1\omega {)}^{\gamma -1}}}+b_0(-\mathsf{i}\omega )\text{ for }\gamma \in (1,2].$$

$$\alpha (\omega )=\frac{(-\mathsf{i}\omega )}{c_{\infty }}\left(\frac{c_{\infty }}{c_0}\sqrt{\frac{1+(c_0/c_{\infty }{)}^2(-\mathsf{i}{\tau }_1\omega )}{1+(-\mathsf{i}{\tau }_1\omega )}}-1\right).$$

$${\mathbf{W}}_{\alpha }:C_0^{\infty }({\Omega }_0)\subseteq L^2({\Omega }_0)\to L^2(\Gamma \times (0,T)):h\mapsto p_{\alpha }{|}_{\Gamma \times (0,T)},$$

$${\mathbf{M}}_{\alpha }g(\cdot ,t):={\int }_0^t{m}_{\alpha }(t,r)g(\cdot ,r)\mathrm{d}r$$

$$\forall (y,t)\in \Gamma \times (0,T):\mathbf{W}h\left(y,t\right)=\frac{1}{2\pi c_0}\frac{\partial }{\partial t}{\int }_0^{c_{0}t}{\int }_{{\mathbb{S}}^{d-1}}\frac{rh\left(y+r\phi \right)}{\sqrt{c_0^2{t}^2-r^2}}\mathrm{d}\phi \mathrm{d}r.$$

$$\forall x\in {\Omega }_0:\left({\mathbf{W}}^{\ast }g\right)\left(x\right)=-\frac{1}{2\pi }{\int }_{\Gamma }{\int }_{|x-y|}^T\frac{{\partial }_{t}g\left(y,t\right)}{\sqrt{c_0^2{t}^2-|x-y{|}^2}}\mathrm{d}t\mathrm{d}S(y).$$

$${\left({\mathbf{D}}_{\alpha }^{\ast }-\frac{1}{c_0}\frac{\partial }{\partial t}\right)}^2{q}_{\alpha }(x,t)-\Delta q_{\alpha }(x,t)=-{\delta }_{\Gamma }(x)g(x,t)\text{ on }{\mathbb{R}}^d\times \mathbb{R},$$

$${\mathbf{W}}_{\alpha }^{\ast }(g)=\frac{\partial q_{\alpha }}{\partial t}(\cdot ,0).$$

$${\left({\mathbf{D}}_{\alpha }+\frac{1}{c_0}\frac{\partial }{\partial t}\right)}^2{q}_{\alpha }^{\ast }(x,t)-\Delta q_{\alpha }^{\ast }(x,t)=-{\delta }_{\Gamma }(x)g(x,T-t)\text{ on }{\mathbb{R}}^d\times \mathbb{R}.$$

$${\mathbf{W}}_{\alpha }^{\ast }(g)=-\frac{\partial q_{\alpha }^{\ast }}{\partial t}(\cdot ,T).$$

$$\forall n\in \mathbb{N}:h_{n+1}^{\delta }=h_n^{\delta }-\lambda {\mathbf{W}}_{\alpha }^{\ast }\left({\mathbf{W}}_{\alpha }{h}_n^{\delta }-g^{\delta }\right).$$

$$\Vert h_{n+1}^{\delta }-g^{\delta }\Vert \le \tau \delta$$

$$h_{n+1}^{\delta }={\mathbf{P}}_C\left(h_n^{\delta }-\lambda {\mathbf{W}}_{\alpha }^{\ast }\left({\mathbf{W}}_{\alpha }{h}_n^{\delta }-g^{\delta }\right)\right),$$

$$h[i]$$

$$x_i$$

$$\mathtt{g}[k,\ell ]$$

$${\phi }_k$$

$$t_{\ell }$$

$$\mathtt{W}:{\mathbb{R}}^{(N_x+1)\times (N_x+1)}\to {\mathbb{R}}^{N_{\phi }\times (N_t+1)},$$

$${\mathtt{M}}_{\alpha }:{\mathbb{R}}^{N_{\phi }\times (N_t+1)}\to {\mathbb{R}}^{N_{\phi }\times (N_t+1)},$$

$$\mathtt{B}:{\mathbb{R}}^{N_{\phi }\times (N_t+1)}\to {\mathbb{R}}^{(N_x+1)\times (N_x+1)},$$

$${\mathtt{P}}_{\mathtt{C}}:{\mathbb{R}}^{(N_x+1)\times (N_x+1)}\to {\mathbb{R}}^{(N_x+1)\times (N_x+1)}.$$

$$\forall n\in \mathbb{N}:{\mathtt{h}}_{n+1}^{\delta }={\mathtt{P}}_{\mathtt{C}}\left({\mathtt{h}}_n^{\delta }-\lambda \mathtt{B}{\mathtt{M}}_{\alpha }^{\ast }\left({\mathtt{M}}_{\alpha }\mathtt{W}{\mathtt{h}}_n^{\delta }-{\mathtt{g}}^{\delta }\right)\right).$$

Full text

Key words: Photoacoustic tomography, image reconstruction, acoustic attenuation, Landweber method, regularization methods.

AMS subject classification: 44A12, 65R10, 92C55.

1 Introduction

Photoacoustic tomography (PAT) is an emerging coupled-physics imaging modality that combines the high spatial resolution of ultrasound imaging with the high contrast of optical imaging (the basic principles are illustrated in Figure 1.1). Potential medical applications include imaging of tumors, visualization of vasculature or scanning of melanoma [6, 46, 57, 75]. In this article we consider PAT using the following general model for acoustic wave propagation in attenuating media,

[TABLE]

Here $h\colon\mathbb{R}^{d}\to\mathbb{R}$ is the photoacoustic (PA) source, $c_{0}>0$ is a constant, and $\mathbf{D}_{\alpha}$ is the time convolution operator associated with the inverse Fourier transform of the complex valued attenuation function $\alpha\colon\mathbb{R}\to\mathbb{C}$. Dissipative pressure wave equation models that can be cast in the form (1.1) can be found in [16, 29, 38, 39, 42, 43, 51, 66, 67]. The particular form of $\alpha$ depends on the used acoustic attenuation model and various different models have been proposed for PAT (see [42] for an overview).

The inverse problem of PAT consists in recovering the source term $h$ from observations of $p_{\alpha}$ on an observation surface $\Gamma\subseteq\mathbb{R}^{d}$ outside its support and for times $t\in(0,{T})$. Taking attenuation into account is essential for high resolution PAT since ignoring attenuation may significantly blur the reconstructed image.

1.1 Our approach

The inverse problem of PAT can be formulated as the problem of estimating $h$ from approximate data $g^{\delta}\simeq\mathbf{W}_{\alpha}h$, where $\mathbf{W}_{\alpha}$ maps the PA source $h$ to the solution of (1.1) restricted to $\Gamma\times(0,T)$. In this paper, we propose the use of regularization methods for stably inverting the operator $\mathbf{W}_{\alpha}$. In particular, we apply the Landweber method, which is a well established regularization method. A main ingredient in the Landweber method is the numerical evaluation of the adjoint $\mathbf{W}_{\alpha}^{*}$. For that purpose, we derive two explicit expressions for the adjoint. The first one takes the form of an explicit formula for the adjoint operator and will be used in our numerical implementation. The second one involves the solution of an adjoint attenuated wave equation. We emphasize that our inversion approach is universal, in the sense, that it can be applied to a wide range of different attenuation models, a general measurement geometry as well as limited data problems.

1.2 Comparison to previous and related work

In the case of vanishing attenuation $\alpha=0$, the attenuated wave equation (1.1) reduces to the standard wave equation

[TABLE]

with sound speed $c_{0}$. Recovering the source term $h$ in (1.2) from boundary data is the standard problem in PAT and various methods for its solution have been derived in the recent years. These approaches can be classified in direct methods, time reversal and iterative approaches. Direct methods are based on explicit solutions for the inverse problem that have been derived in the Fourier domain [2, 27, 48, 78] as well as in the spatial domain [18, 19, 21, 22, 26, 47, 49, 54, 55, 60, 64, 77]. In the time reversal technique, the wave equation (1.2) is solved backwards in time where the measured data are used as boundary values in the time reversed wave equation [11, 19, 32, 56, 70]. Discrete iterative approaches, on the other hand, are usually based on a discretization of the forward problem together with numerical solution methods for solving the resulting system of linear equations [14, 61, 62, 63, 79, 73, 74]. Recently, iterative schemes in a Hilbert space settings have also been introduced and studied; see [5, 7, 25]. In this paper we generalize the iterative Hilbert space approach to attenuating media.

The case of non-vanishing attenuation is much less investigated and existing methods are very different from our approach. One class of reconstruction methods uses the following two-stage procedure: In a first step, by solving an ill-posed integral equation the (idealized) un-attenuated pressure data $p_{0}(z,\,\cdot\,)$ are estimated from the attenuated data $p_{\alpha}(z,\,\cdot\,)$. In the second step, the standard PAT problem is solved. Such a two step method has been proposed and implemented for the power law in [50, 51], and later been used in [4, 42] for various attenuation laws. Compared to two stage approaches, in the single step approach it is easier to include prior information available in the image domain, such as positivity of the PAT source (compare Section 3.1). Furthermore, in the limited data case, where the measurement surface does not fully enclose the PA source, the second step in the two-stage approach is again a non-standard problem for PAT, for which iterative methods can be applied. In such a situation it seems reasonable to directly apply iterative methods to the attenuated data, as considered in the present paper.

A different class of algorithms extends the time reversal technique to the attenuated case (see [1, 3, 10, 31, 35, 40, 59, 72]). Note that the time reversal of the attenuated wave equation yields a noise amplifying equation. Therefore regularization methods have to be incorporated in its numerical solution. Opposed to the time reversal, the adjoint wave equation used in our approach is again damping and no regularization is required for its stable solution. This yields a clear convergence analysis for our method by using standard regularization theory [17, 36, 65]. We are not aware of similar existing results for PAT in attenuating acoustic media. The approaches which are closest to our work seem [33, 34]. In [33] discrete iterative methods are considered, where the problem is first discretized and the adjoint is computed from the discretized problem. Further, both works [33, 34] use attenuation models based on the fractional Laplacian (see [12, 71]) which yields an equation that is non-local in space. It is not obvious how to extend these approaches to model (1.1) which can also include memory.

1.3 Notation

For $k\in\mathbb{N}$, we write $\mathcal{S}(\mathbb{R}^{k+1})$ for the Schwartz space of rapidly decreasing functions $f\colon\mathbb{R}^{k+1}\to\mathbb{C}$, and $\mathcal{S}^{\prime}(\mathbb{R}^{k+1})$ for its dual, the space of tempered distributions. Further we write $\mathcal{F}_{t}$ for the Fourier transform in the temporal variable, defined by $(\mathcal{F}_{t}f)(x,\omega)=\int_{\mathbb{R}}f(t)e^{\mathsf{i}\omega t}f(x,t)\mathrm{d}t$ for $f\in\mathcal{S}(\mathbb{R}^{k+1})$ and extended by duality to tempered distributions. A tempered distribution in $\mathcal{S}^{\prime}(\mathbb{R}^{k+1})$ will be called causal (in the last component) if it vanishes for $t<0$. Finally, for $\alpha\in\mathcal{S}(\mathbb{R})$ we denote by $\mathbf{D}_{\alpha}$ the time convolution operator with kernel $\mathcal{F}^{-1}_{t}(\alpha)$.

1.4 Outline

In Section 2 we formulate the forward operator of the PAT in attenuating acoustic media in a Hilbert space framework. We show that it is continuous between $L^{2}$-spaces and give an explicit expression for its solution. We further derive two expressions for the adjoint operator. In Section 3 we solve the corresponding inverse problem using the Landweber regularization, present convergence results, and give details for its actual implementation. Numerical results are presented in Section 4, and a conclusion is given in Section 5. Finally, in the appendix we present details for the wave equation formulation of the adjoint operator.

2 PAT in attenuating acoustic media

Throughout this paper we assume that $\alpha\colon\mathbb{R}\to\mathbb{C}$ is a weakly causal attenuation function, defined as follows.

Definition 2.1** (Weakly causal attenuation function).**

A function $\alpha\colon\mathbb{R}\to\mathbb{C}$ is called weakly causal attenuation function, if the following assertions hold true:

1. (A1)

$\operatorname*{Re}(\alpha)$* is even and $\operatorname*{Im}(\alpha)$ is odd;* 2. (A2)

$\omega\mapsto\operatorname*{Re}(\alpha(\omega))$* is monotonically increasing for positive $\omega$;* 3. (A3)

$\mathcal{F}^{-1}_{t}(\alpha)(t)$* vanishes for $t<0$.*

Note that (A1) implies that the inverse Fourier transforms of $\alpha$ and $e^{-\alpha(\omega)\,\lvert\mathbf{x}\rvert}$ are real valued. The second condition reflects increasing attenuation with increasing frequency. It is not essential and may be replaced by a similar property. The condition (A3) implies causality of the Greens function $G_{\alpha}$ (i.e. $G_{\alpha}(\,\cdot\,,t)=0$ for $t<0$) and further is equivalent to the Kramers-Krönig relations (see (2.5), (2.6) below). Examples for weakly causal attenuation functions are given in Subsection 2.2.

2.1 Attenuated wave equations

We describe acoustic waves in attenuation media by the integro-differential equation

[TABLE]

Here $s$ is a source term and $\alpha\colon\mathbb{R}\to\mathbb{C}$ a weakly causal attenuation function. For any causal $s\in\mathcal{S}^{\prime}(\mathbb{R}^{d+1})$, the attenuated wave equation (2.1) has a causal solution $p_{\alpha}\in\mathcal{S}^{\prime}(\mathbb{R}^{d+1})$. In particular, this implies the existence of causal Greens function that takes the form (see [42])

[TABLE]

The Greens function represents a spherical wave in attenuating acoustic media that originates at location $\mathbf{x}=0$ and time $t=0$. It satisfies (1.1) with right hand side $s(\mathbf{x},t)=-\delta(\mathbf{x})\delta(t)$. If $K_{\alpha}(\mathbf{x},\,\cdot\,-\tfrac{\left|x\right|}{c_{0}})$ is causal for every $\mathbf{x}$, then the dissipative Green function $G_{\alpha}$, defined in (2.2), has a finite wave front speed $\leq c_{0}$. Attenuation laws with finite wave front speed are called strongly causal in [42].

Throughout we refer to the convolution of the source $s$ with $G_{\alpha}$ as the causal solution of (2.1). For the model of [43] uniqueness is shown in [39]. Note that the proof of [39] can be generalized to any weakly causal attenuation law. This implies uniqueness of a solution of (2.1).

Definition 2.2**.**

Let $\alpha\colon\mathbb{R}\to\mathbb{C}$ be a weakly causal attenuation function. Then, for any $r\in\mathbb{R}$ we define $m_{\alpha}(\,\cdot\,,r)\in\mathcal{S}^{\prime}(\mathbb{R})$ by

[TABLE]

The following result derived in [39, 42] will be frequently used in this paper.

Lemma 2.3** (Relation between attenuated and un-attenuated pressure).**

Let $\alpha$ be a non-vanishing weakly causal attenuation function. Then $m_{\alpha}$ is $C^{\infty}$ on $\mathbb{R}^{2}\setminus\{(0,0)\}$. Moreover,

[TABLE]

where $p_{\alpha}$ and $p_{0}$ denote the causal solutions of (1.1) and (1.2), respectively.

Proof.

See [39, Theorem 1 and Lemma 1]. ∎

2.2 Examples for causal attenuation laws

In this subsection, we give particular examples for causal attenuation laws that we use in this paper: the power law (see [68, 67, 76]), the model of Kowar, Scherzer, Bonnefond (see [42, 43]) and the model of Nachman, Smith and Waag (see [53]).

Remark 2.4** (Kramers-Kronig relations).**

A central property that should be satisfied by (1.1) is causality of the corresponding Greens function $G_{\alpha}$. Causality of $G_{\alpha}$ is equivalent to assumption (A3), the causality of $\mathcal{F}^{-1}_{t}(\alpha)$ (see [42]). Widely use criteria for the causality of $\mathcal{F}^{-1}_{t}(\alpha)$ are the Kramers-Kronig relations (see [44, 45])

[TABLE]

In fact, according to Titchmarsh’s theorem [69, Theorem 95] for square integrable $\alpha$, the causality of $\mathcal{F}^{-1}_{t}(\alpha)$ is equivalent to (2.5) as well as to (2.6). In such a situation, if the imaginary part of the weakly causal attenuation function is known, then its real part is uniquely determined by (2.5).

*Typical acoustic attenuation laws, however, are not square integrable (see the examples below). In this case, the Kramers-Kronig relations cannot be applied directly. Nevertheless, the method of subtractions allows extension to attenuation functions with $\alpha(\omega)=\mathcal{O}(\omega^{n})$ as $\omega\to\infty$ with $n\in\mathbb{N}$ (see [58, Section 1.7]). In such a situation, given the imaginary part $\operatorname{Im}[\alpha]$, the Kramers-Kronig relation (2.5) determines the real part $\operatorname{Re}[\alpha]$ up to $n+1$ additive constants. As a concrete example, consider the case where $\alpha(\omega)=\mathcal{O}(\omega)$. Then the Kramers-Kronig relations yield *

[TABLE]

In particular, the imaginary part of the attenuation function determines its real part provided that $\operatorname{Im}[\alpha(\omega_{0})]$, for some fixed $\omega_{0}$, is given. For the general case $\alpha(\omega)=\mathcal{O}(\omega^{n})$ see [58].

In the following $\alpha_{0}$, $c_{0}$, $c_{\infty}$, $\tau_{1}$ and $\gamma$ denote positive constants.

Example 2.5** (Power law).**

In the power law model, the complex attenuation function takes the form

[TABLE]

Here $(-\mathsf{i}\,\omega)^{\gamma}\coloneqq\left|\omega\right|^{\gamma}\exp(-\mathsf{i}\pi\gamma\operatorname{sign}(\omega)/2)$ and $a_{0},b_{0}$ are arbitrary positive constants. This equation has been considered, for example, in [68, 67, 76]. For tissue, the exponent $\gamma$ in (2.9) is in the range $(1,2]$. Note that for positive $\gamma$ that is not an integer the power law model is weakly causal. In [42] it has been shown that the dissipative waves modeled by (2.9) are strongly causal only if $\gamma\in(0,1)$.

Example 2.6** (Model of Kowar, Scherzer and Bonnefond).**

The model proposed by Kowar, Scherzer and Bonnefond (KSB model) [43] reads

[TABLE]

The KBS model is strongly causal. Because strong causality implies weak causality [43]) the KBS model is also weakly causal. It satisfies the small frequency approximation $\operatorname{Re}[\alpha(\omega)]\asymp a_{0}\sin(\frac{\pi}{2}(\gamma-1))/\left(2\,c_{\infty}\,\tau_{1}\right)\,\left|\tau_{1}\,\omega\right|^{\gamma}$ as $\omega\to 0$. Thus (2.10) behaves as a power law for small frequencies. In fact, the KBS has been proposed as a strongly causal alternative to the power law for the range $\gamma\in(1,2]$, where the power law fails being strongly causal.

Example 2.7** (Model of Nachman, Smith and Waag).**

In the model of Nachman, Smith and Waag (NSW model) with a single relaxation process, the complex attenuation function takes the form (see [53])

[TABLE]

Equation (2.11) and its generalization using $N$ relaxation processes have been derived in [53] based on sound physical principles. The resulting attenuated wave equation is causal and can even be reformulated as differential equation of order $N+2$. In [53, 42] it is shown that the model (2.11) is strongly (and thus weakly) causal provided that $c_{0}<c_{\infty}$. Then the wave front speed is bounded from above by $c_{\infty}$. We note that the attenuation law (i.e. the real part of $\alpha$ of the NSW model (2.11)), satisfies a power law with exponent $\gamma=2$ as $\omega\to 0$.

2.3 The forward operator

In the sequel, we assume the PA source $h$ to be supported in an open set $\Omega_{0}$. We assume that measurements are taken on a piecewise smooth surface $\Gamma\subseteq\partial\Omega$ where $\Omega\subseteq\mathbb{R}^{d}$ is an open set with $\overline{\Omega}_{0}\subseteq\Omega$. Further, let ${T}\geq\operatorname{diam}\left(\Omega\right)/c_{0}$ denote the final measurement time and suppose ${t_{0}}$ is a positive number with ${t_{0}}c_{0}<\operatorname{dist}(\Omega_{0},\Gamma)$.

Definition 2.8** (PAT forward operator).**

We define the PAT forward operator with weakly causal attenuation law $\alpha$ (see Definition 2.1) by

[TABLE]

where $p_{\alpha}$ denotes the causal solution of (1.1). In the case of vanishing attenuation, we write $\mathbf{W}\coloneqq\mathbf{W}_{0}$.

According to Lemma 2.3, the kernel $m_{\alpha}(t,r)$ is smooth for $(t,r)\neq(0,0)$. Moreover,

[TABLE]

defines a bounded linear operator $\mathbf{M}_{\alpha}\colon L^{2}(\Gamma\times({t_{0}},{T}))\to L^{2}(\Gamma\times(0,{T}))$. The representation of the attenuated pressure in terms of the un-attenuated pressure given in Lemma 2.3 therefore shows that $\mathbf{W}_{\alpha}$ is well defined. Moreover, the following result shows that it can be extended to a bounded linear operator on $L^{2}(\Omega_{0})$.

Theorem 2.9** (Mapping properties of the PAT forward operator).**

**

1. (a)

$\mathbf{W}_{\alpha}$* has a unique bounded extension $\mathbf{W}_{\alpha}\colon L^{2}(\Omega_{0})\to L^{2}(\Gamma\times(0,{T}))$;* 2. (b)

$\mathbf{W}_{\alpha}=\mathbf{M}_{\alpha}\circ\mathbf{W}$.

Proof.

(a), (b): It is known that the operator $\mathbf{W}\colon L^{2}(\Omega_{0})\to L^{2}(\Gamma\times({t_{0}},{T}))\colon h\mapsto{p_{0}}|_{\Gamma\times({t_{0}},{T})}$ is well defined, linear, bounded, injective and has closed range (see, for example, [25]). Now suppose $\alpha\neq 0$. Because $m_{\alpha}$ is smooth on $\{(t,r)\neq(0,0)\}$, it follows that $\mathbf{M}_{\alpha}\colon L^{2}(\Gamma\times({t_{0}},{T}))\to L^{2}(\Gamma\times(0,{T}))$ is well defined and bounded. Together with (2.4) this gives (b) and implies the boundedness of $\mathbf{W}_{\alpha}$. In particular, $\mathbf{W}_{\alpha}$ has a unique bounded extension $\mathbf{W}_{\alpha}\colon L^{2}(\Omega_{0})\to L^{2}(\Gamma\times(0,{T}))$. ∎

The well known explicit solution formulas for the standard wave equation give explicit expressions for $\mathbf{W}$. The precise forms of these expression depend on the spatial dimension. For example, in two spatial dimensions we have

[TABLE]

Together with $\mathbf{W}_{\alpha}=\mathbf{M}_{\alpha}\circ\mathbf{W}$ this gives an explicit formula for $\mathbf{W}_{\alpha}$ that can be implemented efficiently.

2.4 The adjoint operator

Because the forward operator $\mathbf{W}_{\alpha}\coloneqq L^{2}(\Omega_{0})\to L^{2}(\Gamma\times(0,\infty))$ is linear and bounded its adjoint $\mathbf{W}_{\alpha}^{*}$ exists and is linear and bounded. In this subsection we give two expressions for the adjoint that can be used for the solution of the inverse problem.

First, we derive an expression for $\mathbf{W}_{\alpha}^{*}$ in the form of an explicit formula. This representation will be used in our numerical reconstruction algorithm.

Theorem 2.10** (Adjoint operator in integral form).**

**

1. (a)

$\mathbf{W}_{\alpha}^{*}\colon L^{2}(\Gamma\times(0,{T}))\to L^{2}(\Omega_{0})$* is well defined and bounded;* 2. (b)

$\mathbf{W}_{\alpha}^{*}=\mathbf{W}^{*}\circ\mathbf{M}_{\alpha}^{*}$; 3. (c)

$\forall g\in L^{2}(\Gamma\times({t_{0}},{T}))\,\forall t\in(0,{T})\colon\mathbf{M}_{\alpha}^{*}g(\,\cdot\,,r)=\int_{r}^{T}m_{\alpha}(t,r)g(\,\cdot\,,t)\mathrm{d}t$.

Proof.

(a): According to Theorem 2.9 (a), $\mathbf{W}_{\alpha}\colon L^{2}(\Omega_{0})\to L^{2}(\Gamma\times(0,{T}))$ is linear and bounded. Therefore, its adjoint is well defined and bounded, too.

(b): Follows from Theorem 2.9 (b).

(c): According to the definition of $\mathbf{M}_{\alpha}\colon L^{2}(\Gamma\times({t_{0}},{T}))\to L^{2}(\Gamma\times(0,{T}))$ we have $\mathbf{M}_{\alpha}g(y,t)\coloneqq\int_{0}^{t}m_{\alpha}(t,r)g(y,r)\mathrm{d}r$. Therefore $\mathbf{M}_{\alpha}^{*}\colon L^{2}(\Gamma\times(0,{T}))\to L^{2}(\Gamma\times({t_{0}},{T}))$ is given by $\mathbf{M}_{\alpha}^{*}g(y,r)=\int_{r}^{T}m_{\alpha}(t,r)g(y,t)\mathrm{d}t$. ∎

Note that the adjoint $\mathbf{W}^{*}$ in the absence of attenuation can be given by an explicit expression and therefore the attenuated adjoint $\mathbf{W}_{\alpha}^{*}=\mathbf{W}^{*}\circ\mathbf{M}_{\alpha}^{*}$ is also given by an explicit formula. The actual expressions for $\mathbf{W}^{*}$ depends on the spatial dimension. For example, in two spatial dimensions, we have

[TABLE]

Our next results show that the adjoint operator can additionally be described by an attenuated wave equation. In absence of attenuation, similar formulations for the adjoint have been derived in [5, 7, 25]. For that purpose, we denote by $\delta_{\Gamma}$ the tempered distribution on $\mathbb{R}^{d}\times\mathbb{R}$ defined by $\left\langle\delta_{\Gamma},\phi\right\rangle=\int_{\mathbb{R}}\int_{\Gamma}\phi(x,t)\,\mathrm{d}S(x)\,dt$ for $\phi\in C_{0}^{\infty}(\mathbb{R}^{d}\times\mathbb{R})$. Furthermore, we denote by $\mathbf{D}_{\alpha}^{\ast}$ the formal $L^{2}$-adjoint of $\mathbf{D}_{\alpha}$ given by the time convolution with the time reversed kernel $\mathcal{F}^{-1}_{t}(\alpha)(-t)$.

Theorem 2.11** (Adjoint operator in wave equation form).**

Let $\alpha$ be any weakly causal attenuation function. For $g\in C_{0}^{\infty}(\Gamma\times(0,{T}))$ let $q_{\alpha}$ be the solution of the adjoint attenuated wave equation

[TABLE]

with $q_{\alpha}(\,\cdot\,,t)=0$ for $t>{T}$. Then,

[TABLE]

Proof.

See Appendix A. ∎

From Theorem 2.11, we immediately obtain the following alternative form.

Corollary 2.12** (Adjoint operator in time reversed wave equation form).**

Suppose the assumptions of Theorem 2.11 are satisfied and let $q_{\alpha}^{*}$ be the causal solution of

[TABLE]

Then we have

[TABLE]

Proof.

Follows from Theorem 2.11 with $q_{\alpha}^{*}(x,t)=q_{\alpha}(x,{T}-t)$. ∎

Note that (2.16), (2.17) and (2.18), (2.19) have a similar form to the time reversal used in [35, 3, 40, 41, 10, 72]. However, unlike the time-reversed wave equation, where the corresponding wave blows up, the adjoint formulation has the same stability properties as the forward equation. Therefore, in contrast to the time reversal procedure, there is no need to include regularization to implement (2.16) or (2.18). Accurate numerical solution of dissipative wave equations for realistic parameters is challenging and numerically quite expensive. Let us consider this issue for the wave equation of Nachmann, Smith and Waag with only one relaxation process. The relaxation time $\tau_{1}$ for fluids is about $100\,\mathrm{ns}$ and the discretization step size should be at least close the relaxation time. This results in a time discretization much finer than usually employed for the simulation of un-attenuated waves, and therefore an increased numerical cost.

3 Solution of the inverse problem

In this section we solve the inverse problems of PAT with attenuation using iterative regularization methods. Our method comes with a clear convergence analysis. We further present details on its actual implementation. Throughout the rest of this paper, we write $\left\|\;\cdot\;\right\|{}$ for the regular $L^{2}$-norms on $L^{2}(\Omega_{0})$ and $L^{2}(\Gamma\times(0,\infty))$ as well as for the operator norm between these two spaces.

3.1 The Landweber method

The Landweber method for the solution of $\mathbf{W}_{\alpha}h\simeq g^{\delta}$ is defined by

[TABLE]

Here $g^{\delta}$ are the noisy data, $0<\lambda\leq\|\mathbf{W}_{\alpha}\|^{-2}$ is the step size, and $h^{\delta}_{0}\coloneqq h^{\delta}$ some initial guess. The superscript $\delta$ indicates the noise level, which means that an estimate $\left\|\mathbf{W}_{\alpha}h^{\star}-g^{\delta}\right\|\leq\delta$ is available, where $h^{\star}$ is the unknown true solution.

The Landweber iteration will be combined with Morozov’s discrepancy principle. According to the discrepancy principle, the iteration is terminated at the index $n=n(\delta,y^{\delta})$, when for the first time

[TABLE]

with some fixed $\tau>1$. From Theorem 2.9 and the general theory of iterative regularization methods, we obtain the following result.

Theorem 3.1** (Convergence of the Landweber iteration).**

Suppose $h\in L^{2}(\Omega_{0})$, $\epsilon>0$ and let $g^{\delta}\in L^{2}((0,{T})\times\Gamma)$ satisfy $\|g^{\delta}-\mathbf{W}_{\alpha}h\|\leq\delta$.

1. (a)

Exact data:* If $\delta=0$, then $(h_{n})_{n\in\mathbb{N}}$ strongly converges to the $h$.* 2. (b)

Noisy data:* Let $(\delta_{m})_{m\in\mathbb{N}}\in(0,\infty)^{\mathbb{N}}$ converge to zero and let $(h_{m})_{m\in\mathbb{N}}\in L^{2}((0,{T})\times\Gamma)^{\mathbb{N}}$ satisfy $\|h_{m}-\mathbf{W}_{\alpha}h\|\leq\delta_{m}$. Then the following hold:*

• $\blacksquare$

The stopping indices $n_{m}:=n_{*}(\delta_{m},h_{m})$ are well defined by (3.2);

• $\blacksquare$

We have $\|h_{n_{m}}^{\delta_{m}}-h\|\to 0$ as $m\to\infty$.

Proof.

According to the Theorem 2.9, the operator $\mathbf{W}_{\alpha}$ is bounded. The claims therefore follow from standard results for iterative regularization methods (see, for example, [17, 36]). ∎

The Landweber method is the most basic iterative regularization method for the solution of inverse problems and behaves very stable due to the smoothing effect of the adjoint. On the other it is quite slow in some applications. Accelerated version include $\nu$-methods [8, 17], the CG Algorithm [37, 28], preconditioned Landweber iterations [15] or Kaczmarz-type iterations [24, 20]. Here we chose the Landweber iteration because the main aim of the present paper is demonstrating the effectiveness of iterative methods for PAT with acoustic attenuation. Furthermore, for the considered application already 10 iterative steps provide very accurate results. Nevertheless, note that generalization to other iterative regularization methods such as the steepest descent or the conjugate gradient method is straight forward.

Another advantage of the Landweber method is that it can easily be combined with a projection step to improve performance. The resulting projected Landweber method reads

[TABLE]

where $\mathbf{P}_{C}$ denotes the projection on a closed convex set $C\subseteq L^{2}(\Omega)$, Actually, in our numerical implementation, we use the projected Landweber method with $C$ being the cone of non-negative functions, which turned out to produce slightly better results than the pure Landweber method with a comparable numerical complexity.

3.2 Implementation of the (projected) Landweber iteration

We outline the implementation for the case that $\Gamma\coloneqq\partial B_{R}(0)$ is a circle of radius $R$ in two spatial dimensions. Extension for more general geometries and higher dimensions are straight forward. Our approach uses the relations $\mathbf{W}_{\alpha}=\mathbf{W}\circ\mathbf{M}_{\alpha}$ and $\mathbf{W}_{\alpha}^{\ast}=\mathbf{M}_{\alpha}^{*}\circ\mathbf{W}^{*}$ (see Theorems 2.9 and Theorems 2.10) that relate the attenuated pressure to the un-attenuated pressure in the direct and adjoint problems, respectively.

For that purpose the PA source $h\colon\mathbb{R}^{2}\to\mathbb{R}$ is represented by a discrete vector ${\tt h}\in\mathbb{R}^{(N_{x}+1)\times(N_{x}+1)}$ obtained by uniform sampling

[TABLE]

Here $(N_{x}+1)^{2}$ is the number of spatial discretization points on an equidistant Cartesian grid. Further, any function $g\colon\partial\Omega\times[0,{T}]\to\mathbb{R}$ is discretely represented by a vector ${\tt g}\in\mathbb{R}^{N_{\varphi}\times(N_{t}+1)}$, with

[TABLE]

Here $N_{\varphi}$ is the number of angular samples (detector locations) and $N_{t}+1$ the number of temporal samples. The sampling conditions obtained in [23] imply that $\Delta x\simeq c_{0}\,\Delta t\simeq R\,\Delta\varphi$, where $\Delta x\coloneqq{2R}/{N_{x}}$ $\Delta t\coloneqq{{T}}/{N_{t}}$ and $\Delta\varphi\coloneqq{2\pi}/{N_{\varphi}}$ yield aliasing free sampling.

The Landweber iteration (3.1) and its projected version (3.3) are implemented by replacing $\mathbf{W}$, $\mathbf{W}^{*}$, $\mathbf{M}_{\alpha}$ and $\mathbf{P}_{C}$ with discrete counterparts

[TABLE]

The resulting discrete (projected) Landweber method is then defined by

[TABLE]

Note that for the sake of computational efficiency the operator $\boldsymbol{\tt{B}}$ will be implemented by a filtered backprojection procedure that which is the exact discrete adjoint of $\boldsymbol{\tt{W}}$. On the other hand, as numerical approximation of $\mathbf{M}_{\alpha}^{*}$ we take the exact adjoint of the discretization of $\boldsymbol{\tt{M}}_{\alpha}$. Finally, the discrete projection will simply be taken as $\boldsymbol{\tt{P}}_{\tt{C}}({\tt h})\coloneqq\operatorname{max}\{0,{\tt h}\}$, which is the projection on convex cone $\tt{C}\coloneqq[0,\infty)^{(N_{x}+1)\times(N_{x}+1)}\subseteq\mathbb{R}^{(N_{x}+1)\times(N_{x}+1)}$. How to implement the other operators will be described in the following subsections.

Remark 3.2** (Numerical complexity).**

Under the reasonable assumption $N_{x}\sim N_{\varphi}\sim N_{t}$, one iterative step (3.13) requires $\mathcal{O}(N_{x}^{3})$ floating point operations (FLOPS) with small leading constants and comparable to the effort of a standard FBP reconstruction algorithm. Since a small number of around 10 turned out to be sufficient, our algorithm is numerically quite efficient.

3.3 Implementation of $\mathbf{W}$ and its adjoint

For numerically approximating the un-attenuated wave operator $\mathbf{W}$, we discretize the explicit formula (2.14). For that purpose, we write (2.14) in the form $\mathbf{W}h=c_{0}^{-1}\partial_{t}\mathbf{A}\mathbf{M}h$, where

[TABLE]

for $(y,t)\in\Gamma\times(0,{T})$. The operator $\mathbf{M}$ is the spherical mean Radon transform and $\mathbf{A}$ the Abel transform (in the second component). We compute discrete spherical means by approximately evaluating (3.14) at the all discretization points $(R\varphi_{k},c_{0}t_{j})$ using the trapezoidal rule for discretizing the integral over $\beta$. The values of $h$ for applying the trapezoidal rule are obtained by using bilinear interpolation of ${\tt h}$. Next for any $k$, the Abel transform is approximately computed by replacing $g(y_{k},\,\cdot\,)$ with the linear interplant through the data pairs $(c_{0}t_{\ell},g(y_{k},c_{0}t_{\ell}))$. Finally, we approximate the time derivative $\partial_{t}$ by finite differences. Inserting these approximations to $\mathbf{W}=c_{0}^{-1}\partial_{t}\mathbf{A}\mathbf{M}$ yields the discretization $\boldsymbol{\tt{W}}$.

The adjoint wave operator $\mathbf{W}^{*}$ is implemented in a similar manner using (2.15) which can be written in the form $\mathbf{W}=-c_{0}^{-1}\mathbf{M}^{*}\mathbf{A}^{*}\partial_{t}$. The operators $\mathbf{A}^{*}$ and $\partial_{t}$ are discretized as above. The adjoint $\mathbf{M}^{*}$ of the spherical mean operator is implemented using a backprojection procedure described in detail in [9, 18].

3.4 Implementation of $\mathbf{M}_{\alpha}$

Recall that for any $(x,t)\in\Gamma\times[0,{T}]$, we have

[TABLE]

The operator $\mathbf{M}_{\alpha}$ is discretized based on these relations by approximately computing $m_{\alpha}(t_{\ell},t_{\ell^{\prime}})$ using (3.17) and then discretizing (3.16). This yields the discrete approximation

[TABLE]

Here $\operatorname{FFT}$ denotes the discrete Fourier transform in the first component and the discrete kernel ${\tt m}_{\alpha}[\ell,\ell^{\prime}]$ is obtained by applying the inverse fast Fourier transform in the first component. Moreover, $\omega[\ell]=\ell\Delta\omega+N_{t}\pi/{T}$ with $\Delta\omega=2\pi/{T}$. Finally, the adjoint $\mathbf{M}_{\alpha}^{\ast}$ is implemented by the adjoint $(\boldsymbol{\tt{M}}_{\alpha}^{\ast}{\tt g})[k,\ell]\coloneqq\Delta t\;\sum_{\ell^{\prime}=0}^{N_{t}}{\tt m}_{\alpha}[\ell^{\prime},\ell]\,{\tt g}[k,\ell^{\prime}]$ of the discrete operator (3.18), (3.19).

4 Numerical results

In this section, we present numerical simulations for PAT with and without attenuation. For all numerical results presented below, the region $\Omega$ is a disc of radius $R$. The final measurement time is taken as ${T}=2R/c_{0}$, where $c_{0}=1540\,\mathrm{m/s}$ is taken as the sound speed in water. For all reconstruction results we take $N_{x}=N_{t}=N_{\varphi}$ in the reconstruction. In order to avoid inverse crime, the data have been computed using a finer temporal discretization.

4.1 Pressure simulation

For the reconstruction results using attenuation data, we will employ the NSW model. It has quadratic frequency dependence for small frequencies. This describes attenuation of water that has an exponent close to $2$ for small frequencies [42, 53]. We use $c_{\infty}=1623\,\mathrm{m/s}$. For simplicity, we restrict ourselves to a single relaxation process in the NSW model. For the relaxation time $\tau_{1}$, we consider two cases, for which we also consider different radii of the measurement circle:

• $\blacksquare$

Case 1: $R=5\,\mathrm{cm}$ and $\tau_{1}=100\,\mathrm{ns}$;

• $\blacksquare$

Case 2: $R=5\,\mathrm{mm}$ and $\tau_{1}=1\,\mathrm{ns}$.

Figure 3.1 shows the corresponding kernel function $m_{\alpha}(t,\,\cdot\,)$ for the above different relaxation times. We see that the support of $m_{\alpha}(t,\,\cdot\,)$ increases with the relaxation time indicating that the larger relaxation time corresponds to stronger attenuation.

In Figure 4.1, we compare the (noisy) un-attenuated pressure data measured $p_{0}(\mathbf{x},\,\cdot\,)$ at location $\mathbf{x}=(R,0)$ with (noisy) attenuated pressure data $p_{\alpha}(\mathbf{x},\,\cdot\,)$ according to the NSW model for the phantom shown in Figure 4.2. We also compare it to the pressure data obtained with the KSB model and the power law with exponent $\gamma=2$. The parameter settings of the KSB and the power law models have been chosen such that the real and imaginary parts of the complex attenuation laws are as close as possible to the one of the NSW law for small frequencies. For the strong attenuation case (Figure 4.1, right), we see that all attenuated pressure data are very similar. Indeed, if we simulate noisy data via the NSW model and then estimate the initial data via the power law, the KSB law or the NSW law, then the results would hardly be distinguishable. However, the left picture in Figure 4.1 shows that this is not true in the small attenuation case where the different attenuation laws yield quite different attenuated pressure signals. Note that for the power law we actually implemented a causal form, where we have truncated $m_{\alpha}(t,r)$ for $r>t$ after evaluating (3.17).

4.2 Reconstruction results for strong attenuation

The numerical simulations that we present in this subsection correspond to strong attenuation with a relaxation time $\tau_{1}=100\,\mathrm{ns}$. The radius of the region of interest is taken as $R=5\,\mathrm{cm}$ and we take $N_{x}=N_{t}=N_{\varphi}=600$.

The numerical phantom and the numerical results using the projected Landweber iteration with and without attenuation are presented in Figure 4.2. We see that the reconstructions using the NSW model (top right) yields a smoother results than in the absence of attenuation (bottom left). In the case with attenuation the thin concentric annuli cannot be resolved, they appear as single thick blurred annulus. Moreover, small or thin structures are blurred and provide less contrast in the case of attenuation. We also applied the projected Landweber iteration using the un-attenuated wave equation to the attenuated data. The reconstruction shown in the bottom right image in Figure 4.2 indicates that thin and long structures are strongly blurred and displaced. Actually, details with small diameter cannot be estimated reliably which clearly indicates that attenuation has to be taken into account. This also reflects that attenuated data are not only smoothed but also displaced. The artifacts in the mixed reconstruction might be reduced by shifting the pressure data appropriately. Indeed, heuristic rules performing such a shift are often applied in practice. However, as the shift depends on the location and the frequency content of the object applying a reasonable shift seems a non-trivial issue that is not required at all in our approach.

We are not aware how to exactly choose the free parameters $c_{0}$ and $c_{\infty}$ in order to accurately model acoustic attenuation in water or soft tissue. To investigate the effect of changing these parameters, we also perform simulations with a significantly increased value of $c_{\infty}=3080\,\mathrm{m/s}$. From the results showing in Figure 4.3, one observes significantly increased attenuation compared to the value $c_{\infty}=1623\,\mathrm{m/s}$ (see top right image in Figure 4.2).

4.3 Reconstruction results for weak attenuation

Now we present simulations of the NSW model for weak attenuation case with relaxation time $\tau_{1}=1\,\mathrm{ns}$. As a consequence, we have to use a finer time discretization for calculating ${\tt m}_{\alpha}$ in (3.18), (3.19). In order to keep the computational expenses reasonable, we decreased the radius to $R=5\,\mathrm{mm}$. From the numerical results presented in Figure 4.4, we see that the attenuated case again yields a smoother reconstructions than in the absence of attenuation. In contrast to the strong attenuation case, the very thin concentric annuli located in the upper half of the image of $f$ can still be resolved; the contrast now even seems better than in the dissipation free case. Also, the very small elliptic structures can be estimated with high quality. The thinner concentric annuli located in the lower half of the image of $f$ cannot be resolved. As in the case of strong attenuation, if the standard wave reconstruction is applied attenuated data, then the reconstruction of thin and long structures is blurred and displaced.

4.4 Reconstruction results for limited view data

Finally, we perform reconstructions using limited view (or limited angle) data, where the detector positions $\mathbf{x}=R(\cos\varphi,\sin\varphi)$ are located on a half circle with $\varphi\in[0,\pi]$. The reconstruction results using the projected Landweber method are shown in Figure 4.5. The top row considers the the strong attenuation case and the bottom row the weak attenuation case. In both cases, the result are compared to the cases without attenuation. All reconstructions show the typical limited view artefacts, even in the absence of attenuation. Further, one notices that the results using un-attenuated data yield a little better contrast for long thin structure and much better contrast for structures with small diameters than the ones with attenuated data. Again, we see if the attenuation is not too strong, then attenuation leads to smoother images and even partly better results than in the absence of attenuation.

5 Conclusion

In this paper, we developed iterative regularization methods for PAT in attenuating media. This comes with a clear convergence theory in the Hilbert space framework that is not shared by any other existing approach. For the sake of clarity, we focused on the Landweber method. Generalizations to other regularization techniques such the CG method or Tikhonov regularization are subject of future research. A main ingredient of these regularization methods is the evaluation adjoint of the forward operator. For that purpose, we developed two formulations for the adjoint: One takes the form of an explicit formula whereas the second one involves the solution of an adjoining wave equation. While the proposed method can equally be applied for general admissible attenuation models, in or numerical numerical results, we focused on the widely accepted attenuation model of Nachman, Smith and Waag [53]. A detailed comparison of reconstructions with different attenuation laws and different reconstruction algorithms is intended for future research. The presented numerical results clearly demonstrate that for moderate attenuation even small structures are estimated well with our method. On the other hand, our results show that not accounting for attenuation yields severe artifacts due to dispersion. This clearly demonstrated the necessity of taking correct attenuation models into account in the inversion process. Moreover, our numerical experiments indicate that for weak attenuation the results are even better than in the absence of attenuation.

Acknowledgement

L.V. Nguyen’s research is partially supported by the NSF grants DMS 1212125 and DMS 1616904. He also thanks the University of Innsbruck for financial support and hospitality during his visit.

Appendix A Adjoint attenuated wave equation

Let $\Omega$ be an open set such that $\partial\Omega$ is a closed smooth surface in $\mathbb{R}^{d}$. For notational convenience, we will denote $\Omega_{c}\coloneqq\mathbb{R}^{d}\setminus\overline{\Omega}$. For $g\in C_{0}^{\infty}(\partial\Omega\times\mathbb{R})$ consider the equation

[TABLE]

Here the notation $u(\,\cdot\,,t)=0$ for $t\ll 0$ means that there exists some $t_{0}\in(-\infty,0]$ such that $u(\,\cdot\,,t)=0$ for all $t<t_{0}$. In this appendix, we show regularity of solutions of (A.1) and demonstrate that its solution defines the adjoint of $\mathbf{W}_{\alpha}$.

A.1 Regularity and classical solution of (A.1)

For any $g\in C_{0}^{\infty}(\partial\Omega\times\mathbb{R})$, the source term $-\delta_{\partial\Omega}(x)\,g(x,t)$ is a tempered distribution that vanishes for sufficiently small $t$. Hence (A.1) has a unique distributional solution

[TABLE]

where $G_{\alpha}\in\mathcal{S}^{\prime}(\mathbb{R}^{d}\times\mathbb{R})$ is the causal Greens function of the attenuated wave equation (2.1). As first result in this section we show that the restrictions of the solution to $\Omega$ and $\Omega_{c}$ are smooth and both can be smoothly extended to $\partial\Omega$.

Theorem A.1** (Regularity of solutions of (A.1)).**

For any $g\in C_{0}^{\infty}(\partial\Omega\times\mathbb{R})$, (A.1) has a unique solution $u\in C^{\infty}((\mathbb{R}^{d}\setminus\partial\Omega)\times\mathbb{R})$. Further, $u$ can be extended continuously to $\mathbb{R}^{d}\times\mathbb{R}$, and ${\nabla u}|_{\Omega}$ and ${\nabla u}|_{\Omega_{c}}$ can be extended continuously to $\overline{\Omega}\times\mathbb{R}$ and $\overline{\Omega_{c}}\times\mathbb{R}$, respectively.

Proof.

Let $u\in\mathcal{S}^{\prime}(\mathbb{R}^{d}\times\mathbb{R})$ denote the unique distributional solution of (2.1) with $s(x,t)=\delta_{\partial\Omega}(x)g(x,t)$, given by A.2. In order to obtain the regularity of $u$ we work in the frequency domain and employ the theory of single and double layer potentials for the Helmholtz equation. For that purpose note that the temporal Fourier transform of $G_{\alpha}$ is given by $\Phi_{\alpha}(x,\omega)={e^{\mathsf{i}k(\omega)\left|\mathbf{x}\right|}}/{\left|\mathbf{x}\right|}$ where $k(\omega)\coloneqq\mathsf{i}\alpha(\omega)+\omega/c_{0}$. Further, write $\hat{u}$ and $\hat{g}$ for the temporal Fourier transform of $u$ and $g$, respectively. Then

[TABLE]

which is recognized as a single layer potential for the Helmholtz equation with density $\hat{g}$. Since $\hat{g}(\,\cdot\,,\omega)\in C^{1}(\partial\Omega)$, the theory of single and double layer potentials (see, for example, [13]) shows the following:

• $\blacksquare$

${[\hat{u}(\,\cdot\,,\omega)]}|_{\partial\Omega}=0$ and ${[\partial_{\nu}\hat{u}(\,\cdot\,,\omega)]}|_{\partial\Omega}=\hat{g}(\,\cdot\,,\omega)$;

• $\blacksquare$

$\nabla\hat{u}(\,\cdot\,,\omega)\in C(\overline{\Omega})$ and $\nabla\hat{u}(\,\cdot\,,\omega)\in C(\overline{\Omega_{c}})$;

• $\blacksquare$

For some function $C(\omega)$ that is at most polynomially growing, we have

[TABLE]

Here and below the bracket $[v]$ denotes the jump of a function $v\in C(\overline{\Omega}\times\mathbb{R})\cap C(\overline{\Omega_{c}}\times\mathbb{R})$ across the surface $\partial\Omega$ (from inside out). Next note that $\omega\mapsto\|\hat{g}(\,\cdot\,,\omega)\|_{C^{1}}$ decays faster than any polynomial. Therefore, (A.4) implies that $u$ is infinitely differentiable on $C(\mathbb{R}^{d})$ and that ${\nabla u}|_{\Omega}$ and ${\nabla u}|_{\Omega_{c}}$ are infinitely differentiable on $C(\Omega)$ and $C(\Omega_{c})$ (with respect to the time variable). ∎

In the following, we call $u\in C^{\infty}((\mathbb{R}^{d}\setminus\partial\Omega)\times\mathbb{R})$ a classical solution of (A.1), if $u$ can be extended continuously to $\mathbb{R}^{d}\times\mathbb{R}$, ${\nabla u}|_{\Omega}$ and ${\nabla u}|_{\Omega_{c}}$ can be continuously extended to $\partial\Omega$, and

[TABLE]

Using Theorem A.1, we one can show the following existence and uniqueness result for solutions of (A.5).

Corollary A.2** (Existence and uniqueness of (A.5)).**

Any classical solution of (A.5) is a distributional solution of (A.1), and vice versa. In particular, for any $g\in C_{0}^{\infty}(\mathbb{R}\times\partial\Omega)$, (A.5) is uniquely solvable.

Proof.

According to Theorem A.1 and its proof, any solution of (A.1) is a solution of (A.5). Conversely, note that $u$ is a distributional solution of (A.1) if and only

[TABLE]

where $\Phi$ is the solution of \bigl{(}\mathbf{D}_{\alpha}^{*}-\tfrac{1}{c_{0}}\tfrac{\partial}{\partial t}\bigr{)}^{2}\Phi-\Delta\Phi=\phi and $\Phi(\,\cdot\,,t)=0$ for $t>T$. Using integration by parts one verifies that any solution of (A.5) satisfies (A.6) and therefore also (A.1). ∎

A.2 Proof of Theorem 2.11

Let $h\in C^{\infty}_{0}(\Omega)$ and let $p_{\alpha}$ be the solution of the attenuated wave equation (1.1). Multiplication of (1.1) with a test function $\psi\in C_{0}^{\infty}(\mathbb{R}^{d}\times\mathbb{R})$ and integrating by parts shows

[TABLE]

Now suppose $g\in C_{0}^{\infty}(\Gamma\times(0,{T}))$ and let $q_{\alpha}\in C^{\infty}(\mathbb{R}^{d}\times\mathbb{R})$ be the solution of the adjoint attenuated wave equation (2.16). Because $p_{\alpha}$ vanishes for $t<0$ and $q_{\alpha}$ vanishes for $T>0$ , identity (A.7) also holds for $\psi=q_{\alpha}$. This gives

[TABLE]

Using two times integration by parts in the first term in (A.8) with respect to $t$, using the definition of $\mathbf{D}_{\alpha}^{*}$, and recalling that $q_{\alpha}$ solves the adjoint attenuated wave equation (2.16) shows

[TABLE]

Consequently, for every $g\in C_{0}^{\infty}(\Gamma\times(0,\infty))$, it holds that

[TABLE]

As the last identity holds on a dense subset of $L^{2}(\Omega)$, this shows the expression (2.16), (2.17) for $\mathbf{W}_{\alpha}^{*}g$ and concludes the proof of Theorem 2.11.