**5. Conclusions**

In this paper, we investigated FFD-PAT, where projection data of acoustic pressure are measured. For the first time, the variable speed of sound has been taken into account for this kind of setup. We developed a two-step reconstruction procedure that recovers the 3D pressure *p*(· , *T*) in a first step, which is then used as input for a finite time wave inversion problem in a second step. As the main theoretical contribution of this paper we prove uniqueness and stability results for the finite time wave inversion problem. For the actual solution, we propose the steepest descent iteration which we found to be numerically efficient and stable for FFD-PAT. Moreover, our numerical results demonstrate that ignoring sound speed variations significantly degrades the reconstruction quality. We point out that the novelties of the paper are the development of a reconstruction method together with a mathematical uniqueness and stability analysis of the arising final time wave inversion problem for FFD-PAT with variable sound speed. The experimental feasibility of FFD-PAT has been demonstrated previously in [30,31]. In future work we will perform experiments using FFD-PAT with spatially variable sound

speed and experimentally verify our two-step algorithm. Moreover, in upcoming work we will study limited view problems for FFD-PAT which naturally arise in applications, for instance in the case of breast imaging.

**Author Contributions:** M.H. and G.Z. developed the reconstruction algorithms, the numerical implementation and performed the numerical experiments; L.N. established the uniqueness and stability results; M.H., G.Z., L.N. and R.N. wrote the paper; M.H. and R.N. supervised the project.

**Funding:** G.Z. and M.H. acknowledge support of the Austrian Science Fund (FWF), project P 30747-N32. The research of L.N. is supported by the National Science Foundation (NSF) Grants DMS 1212125 and DMS 1616904. The work of R.N. has been supported by the FWF, project P 28032.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **Appendix A. k-Space Method**

We briefly describe the k-space method for the 3D wave Equations (1)–(3) as we use it for the numerical computation of **W***T* and **W**∗ *T*. The k-space method is an attractive alternative to standard methods using finite differences, finite elements or pseudospectral methods, since it does not suffer from numerical dispersion [52,53]. It utilizes the decomposition *p* (**<sup>x</sup>**, *t*) = *<sup>w</sup>*(**<sup>x</sup>**, *t*) − *<sup>v</sup>*(**<sup>x</sup>**, *t*), where *v*, *w* are defined by *<sup>w</sup>*(**<sup>x</sup>**, *t*) := (*c*2 0/*c*<sup>2</sup>(**x**)) *p*(**<sup>x</sup>**, *t*) and *<sup>v</sup>*(**<sup>x</sup>**, *<sup>t</sup>*)=(*c*<sup>2</sup> 0/*c*<sup>2</sup>(**x**) − <sup>1</sup>)*p*(**<sup>x</sup>**, *t*) where *c*0 := max{*c*(**x**) | **x** ∈ R<sup>3</sup>} denotes maximal speed of sound. It can be checked that with this definition of *v* and *w* the wave equation with variable speed of sound splits into the system *wtt*(**<sup>x</sup>**, *t*) − *c*2 0<sup>Δ</sup> *<sup>w</sup>*(**<sup>x</sup>**, *t*) = −*c*<sup>2</sup> <sup>0</sup>Δ*v*(**<sup>x</sup>**, *t*) and *<sup>v</sup>*(**<sup>x</sup>**, *t*) = *c*2 0<sup>−</sup>*c*<sup>2</sup>(**x**) *c*2 0 *<sup>w</sup>*(**<sup>x</sup>**, *t*). In the k-space method we use the time stepping formula

$$w(\mathbf{x}, t + h\_l) = 2w(\mathbf{x}, t) - w(\mathbf{x}, t - h\_l) - 4\mathbf{F}\_{\tilde{\xi}}^{-1} \left\{ \sin \left( \frac{c\_0 |\tilde{\xi}| h\_l}{2} \right)^2 \mathbf{F}\_{\mathbf{X}} \left\{ w(\mathbf{x}, t) - v(\mathbf{x}, t) \right\} \right\},\tag{A1}$$

where **F***x* and **F**−<sup>1</sup> *ξ* denote the Fourier transform and its inverse with respect to space and frequency variables **x** and *ξ* and *ht* is the time step size. This equivalent formulation motivates the following algorithm for numerically solving the wave equation.

#### **Algorithm A1:** k-space method for numerically solving (1)–(3).

*(S1) Define initial conditions:*

$$\begin{aligned} w(\mathbf{x}, -\mathbf{h}\_{\mathbf{f}}) &= w(\mathbf{x}, 0) = c\_0^2 / c^2(\mathbf{x}) f(\mathbf{x}), \\ v(\mathbf{x}, 0) &= (c\_0^2 / c^2(\mathbf{x}) - 1) f(\mathbf{x}) \end{aligned}$$

*(S2) Set t* = 0 *(S3) Compute <sup>w</sup>*(**<sup>x</sup>**, *t* + *ht*) *according to Equation* (A1) *(S4) Compute <sup>v</sup>*(**<sup>x</sup>**, *t* + *ht*) := (*c*<sup>2</sup>(**x**)/*c*<sup>2</sup> 0 − <sup>1</sup>)*w*(**<sup>x</sup>**, *ht*) *(S5) Compute p*(**<sup>x</sup>**, *t* + *ht*) := *<sup>w</sup>*(**<sup>x</sup>**, *t* + *ht*) − *<sup>w</sup>*(**<sup>x</sup>**, *t* + *ht*)*(S6) Substitute t by t* + *ht and go back to step (3).*
