*2.3. Two-Step Reconstruction*

One possible approach to solve the inverse problem of FFD-PAT is to directly recover *f* from data in (8) via iterative methods. Typically, each iteration step requires the evaluation of the forward operator **X***R***W***T* and its adjoint (**<sup>X</sup>***R***W***T*)<sup>∗</sup> = **<sup>W</sup>**<sup>∗</sup>*T***X**<sup>∗</sup>*R*. In this paper, we consider a two-step approach where we first invert **X***R* via a direct method and then use an iterative method to invert **W***T*. This avoids repeated and time consuming evaluation of **X***R* and its adjoint. The proposed reconstruction method consists of the following two steps.


$$\mathbf{W}\_{\rm{T}}f(\mathbf{x},y,z) \simeq \mathbf{X}^{\dagger}\underline{g}(\mathbf{x},y,z) := \frac{1}{2\pi^{2}} \int\_{0}^{\pi} \mathrm{P.V.} \int\_{\mathbb{R}} \frac{\left(\partial\_{s}\underline{g}\right)\left(\theta,s,z\right)ds}{\left(\mathbf{x}\cos(\theta) + \mathbf{y}\sin(\theta)\right) - s} d\theta \,\,\boldsymbol{\theta}$$

where P.V. denotes the principal value integral.


$$\text{Recover } f \text{ from data } \quad h = \mathbf{W}\_T f + \text{noise}.\tag{9}$$

To the best of our knowledge, the final time wave inversion problem (9) has not been considered so far and its investigation will be the main theoretical focus of this work.

For solving the wave inversion problem (9), we propose iterative solution methods that are described in Section 3. As main theoretical results, we derive its uniqueness and stability.

Another possible approach for solving the first step would be to work with the exterior Radon transform [43–45]. Instead, we invert the standard Radon transform after setting the missing values of **XW***T f* to zero. We observe numerically that the missing values are indeed close to zero for large enough measurement time *T*. Although we do not have a rigorous mathematical proof for this claim, it is supported by at least two facts. First, in the case of a non-trapping sound speed, the known decay estimate [46] for the solution *p* of (1)–(3) with initial data *f* supported in *B*(0, *R*) states

$$\left|\frac{\partial^{k+|m|}p(\mathbf{x},t)}{\partial^k\_t\bar{\partial}^{|m|}}\right| \le \gamma e^{-\beta t} \left\|f\right\|\_2 \quad \text{for } (\mathbf{x},t) \in B(0,R) \times (T,\infty) \tag{10}$$

for all (*k*, *m*) ∈ N2. Here *β* > 0 is a constant only depending on *c*(**x**) and *T*, and *γ* is a constant depending on *B*(0, *<sup>R</sup>*). Second, in the case of constant sound speed, the X-ray transform reduces (1)–(3) to a 2D wave equation with initial data **X***f* supported in a disc of radius *R*. Thus, in the constant sound speed case, **<sup>X</sup>**(*p*(· , *t*)) describes 2D propagation which rapidly decays inside the disc. Formulating and proving precise statements in such a direction, however, is an interesting unsolved problem.

#### **3. Final Time Wave Inversion Problem for Variable Sound Speed**

In this section we study the final time wave inversion problem (9), where the forward operator **W***T* : *f* → *p*(· , *T*) is defined by (5). According to standard results for the wave equation [47], the forward operator extends to a bounded linear operator **W***T* : *L*<sup>2</sup>(*B*(0, *R*)) → *<sup>L</sup>*<sup>2</sup>(R<sup>3</sup>). Below we establish uniqueness and stability results and derive an iterative reconstruction algorithm. Note that for constant sound speed, recovering the function *f* from the solution at time *T* of (1) with initial data (0, *f*) instead of (*f* , <sup>0</sup>), is equivalent to the inversion from spherical means with fixed radius. Uniqueness and an inversion method for this problem has been obtained in the classical book of Fritz John [48]. However, neither for the case of initial data (*f* , 0) nor in the variable sound speed case we are aware of any similar results.

#### *3.1. Uniqueness and Stability Theorem*

The following theorem is the main theoretical result of this paper and states that the final time wave inversion problem (9) has a unique solution that stably depends on the right-hand side.

**Theorem 1** (Uniqueness and stability of inverting **W***T*)**.** *For a non-trapping speed of sound, the operator* **W***T* : *L*<sup>2</sup>(*B*(0, *R*)) → *L*<sup>2</sup>(R<sup>3</sup>) *is injective and bounded from below.*

The proof of Theorem 1 is presented in the following Section 3.2. See [19,41,49] for related results and methods in the case of standard PAT data. Theorem 1 states that

$$b := \inf \left\{ \frac{\|\mathbf{W} \tau f\|\_{L^2(\mathbb{R}^3)}}{\|f\|\_{L^2(B(0,R))}} \; \middle| \; f \in L^2(B(0,R)) \right\} > 0 \,. \tag{11}$$

In particular, **W***T* : *L*<sup>2</sup>(*B*(0, *R*)) → Ran(**<sup>W</sup>***T*) has a bounded inverse, where Ran(**<sup>W</sup>***T*) denotes the range of **W***T*. The latter implies that gradient based iterative methods for (9) converge linearly, similar to the case of standard PAT data; see [38].

#### *3.2. Proof of Theorem 1*

Let *c* be a non-trapping sound speed. We first prove two crucial Lemmas, from which Theorem 1 follows rather quickly.

#### **Lemma 1.** *Assume that* **W***T f* = 0*, then f* = 0*. That is,* **W***T is injective.*

**Proof.** Let *p* denote the solution of (1)–(3). We will construct a solution *p*¯ of the wave equation which is periodic in time with period 4 *T* such that *p*¯ = *p* on R<sup>3</sup> × [0, *<sup>T</sup>*]. Once this is done, we obtain *f* = *p*¯(· , 0) = *p*¯(· , 4*nT*) for any *n*. Using (10), we arrive at

$$\forall \mathbf{x} \in \mathbb{R}^3 \colon \quad f(\mathbf{x}) = \lim\_{n \to \infty} \not p(\mathbf{x}, \mathbf{4}^n T) = 0 \dots$$

It remains to construct the above-mentioned solution *p*¯ of the wave equation. The idea is to properly reflect the solution *p* in the time variable *t* through the time moments *t* = *T*, 2 *T*, ... , as follows. We first construct *p*¯ on [0, 2 *T*] by the odd reflection of *p* through the moment *t* = *T*: *p*¯(· , *T*) = *p*(· , *T*) for *t* ∈ [0, *T*] and *p*¯(· , *T*) = <sup>−</sup>*p*(·, 2 *T* − *t*) for all *t* ∈ [*<sup>T</sup>*, 2 *<sup>T</sup>*]. Since *p*(· , *T*) = 0 on R3, we obtain that *p*¯ and *p*¯*t* are continuous at *t* = *T*. Therefore, *p* is continuous on [0, 2 *T*] and solves the wave equation on that interval. Next note that *p*¯*t*(·, 2 *T*) = <sup>−</sup>*p*¯*t*(·, 0) = 0 on R3. By the even reflection through *t* = 2*T*: *p*¯(· , *T*) = *p*¯(· , 4 *T* − *t*) for all *t* ∈ [2*T*, 4 *<sup>T</sup>*], we obtain that *p*¯ is a solution of the wave equation in [0, 4 *<sup>T</sup>*]. Finally, we extend the solution by periodicity with period 4 *T*. Noting that *p*¯(· , 0) = *p*¯(· , 4 *T*) and *p*¯*t*(· , 0) = *p*¯*t*(· , 4 *T*) = 0, we obtain that *p*¯ and *p*¯*t* are continuous for all time and *p*¯ satisfies the wave equation in R<sup>3</sup> × R+. This finishes our proof.

**Lemma 2.** *There is a constant C and a pseudo-differential operator* **K***T of order at most* −1 *such that*

$$\forall f \in L^2(B(0, \mathbb{R})): \quad \|f\|\_{L^2(\Omega)} \le 2\left(\|\mathbf{W}\_T f\|\_{L^2(\mathbb{R}^3)} + \|\mathbf{K}\_T f\|\_{L^2(B(0, \mathbb{R}))}\right).$$

**Proof.** Let *p* denote the solution of wave Equations (1)–(3) and recall the parametrix formula *p*(**<sup>y</sup>**, *t*) = 1 (<sup>2</sup>*π*)<sup>3</sup> ∑*<sup>σ</sup>*=<sup>±</sup> R<sup>3</sup> *<sup>a</sup>σ*(**<sup>y</sup>**, *t*, *ξ*)*eiφ*±(**<sup>y</sup>**,*T*,*ξ*) ˆ *f*(*ξ*) *dξ* = ∑*<sup>σ</sup>*=<sup>±</sup> *pσ*(**<sup>y</sup>**, *t*); see [47]. Here, the phase function *φ*± solves the eikonal equation *∂tφ*±(**<sup>y</sup>**, *t*, *ξ*) ± *<sup>c</sup>*(**y**)|∇**y***φ*±(**<sup>y</sup>**, *t*, *ξ*)| = 0 for (**<sup>y</sup>**, *t*) ∈ R<sup>3</sup> × R+ with the initial condition *φ*±(**<sup>x</sup>**, 0, *ξ*) = **x** · *ξ*. The amplitude function is a classical symbol *<sup>a</sup>*±(**<sup>y</sup>**, *t*, *ξ*) = ∑∞*<sup>k</sup>*=<sup>0</sup> *<sup>a</sup>*−*k*,<sup>±</sup>(**<sup>y</sup>**, *t*, *ξ*), where *a*−*k* is homogeneous of order −*k* in *ξ*. Its leading term *<sup>a</sup>*0,± satisfies the transport equation

$$\left(\partial\_{t}\phi\_{\pm}(\mathbf{y},t,\xi)\partial\_{t} - c^{2}(\mathbf{y})\nabla\_{\mathbf{y}}\phi\_{\pm}(\mathbf{y},t,\xi) \cdot \nabla\_{\mathbf{y}} + \mathbb{C}\_{\pm}(\mathbf{y},t,\xi)\right)a\_{0,\pm}(\mathbf{y},t,\xi) = 0,\tag{12}$$

with the initial condition *<sup>a</sup>*±,<sup>0</sup>(**<sup>x</sup>**, 0, *ξ*) = 1/2, see [41]. Here, *<sup>C</sup>*(**<sup>y</sup>**, *ξ*, *t*) only depends on the sound speed *c* and the phase function *φ*<sup>±</sup>. Let us denote by *γ***<sup>x</sup>**,*ξ* the unit speed geodesics originating at **x** along the direction *ξ*. Then, *γ***<sup>x</sup>**,*ξ* is a characteristics curve of the above transport equation; that is, (12) reduces to a homogeneous ODE on each geodesic curve.

We then write

$$\mathbf{W}\_{T}(f)(\mathbf{y}) \quad = \frac{1}{(2\pi)^{3}} \sum\_{\sigma=\pm} \int\_{\mathbb{R}^{3}} a\_{\sigma}(\mathbf{y}, T, \mathfrak{z}^{\mathfrak{z}}) e^{i\boldsymbol{\upphi}\_{\pm}(\mathbf{y}, T, \mathfrak{z}^{\mathfrak{z}})} \widehat{f}(\mathfrak{z}) \, d\mathfrak{z} = \sum\_{\sigma=\pm} \mathbf{W}\_{\sigma}(f)(\mathbf{y}) .$$

Each operator **W**± is a Fourier integral operator (FIO) with the canonical relation given by the pairs (**<sup>y</sup>**±, *λη*±; **x**, *λξ*) for any *λ* ∈ R, *ξ*, *η* unit vectors, **y**± = *γ***<sup>x</sup>**,*ξ* (<sup>±</sup>*T*), and *η*± = *γ*˙ **<sup>x</sup>**,*ξ* (<sup>±</sup>*T*). Let R<sup>3</sup> be equipped with the metrics *c*<sup>−</sup><sup>2</sup>(**x**) *d***x**2. Then, (**<sup>y</sup>**±, *η*±) is obtained by translating (**<sup>x</sup>**, *ξ*) on the geodesic *γ***<sup>x</sup>**,±*ξ* by the distance *T*. From the initial condition of *φ*± and *<sup>a</sup>*0,± we see that, up to lower order terms,

$$p\_{-}(\mathbf{x},0) = p\_{+}(\mathbf{x},0) = \frac{1}{2}f(\mathbf{x})\,. \tag{13}$$

Heuristically, under Equations (1)–(3), each singularity of *f* at (**<sup>x</sup>**, *ξ*) is broken into two equal parts. They propagate along the geodesic *γ***<sup>x</sup>**,*ξ* in the opposite directions ±*ξ* to generate a singularity of **<sup>W</sup>***T*(*f*) at (**<sup>y</sup>**±, *η*±).

From the standard theory of FIOs (see [50]), the adjoint **W**∗± translates (**<sup>y</sup>**±, *η*±) back to (**<sup>x</sup>**, *ξ*) and **<sup>W</sup>**∗±**<sup>W</sup>**± is a pseudo differential operator. On the other hand, **<sup>W</sup>**∗∓**<sup>W</sup>**± is a FIO whose canonical relation consists of the pairs (**<sup>y</sup>**, *η*; **x**, *ξ*) given by **y** = *γ***<sup>x</sup>**,*ξ* (<sup>±</sup>2*<sup>T</sup>*), and *η* = *γ*˙ **<sup>x</sup>**,*ξ* (<sup>±</sup>2*<sup>T</sup>*). That is, **<sup>W</sup>**∗±**<sup>W</sup>**∓ is an infinitely smoothing operator on *B*. Therefore, microlocally, we can write **<sup>W</sup>**<sup>∗</sup>*T***W***<sup>T</sup> f* = **<sup>W</sup>**<sup>∗</sup>+**<sup>W</sup>**+(*f*) + **W**∗ − **<sup>W</sup>**−(*f*). We will show that the principal symbol *<sup>θ</sup>*±(**<sup>x</sup>**, *ξ*) of **<sup>W</sup>**∗±**<sup>W</sup>**± satisfies *<sup>θ</sup>*±(**<sup>x</sup>**, *ξ*) = 1/4. This result can be intuitively understood as follows. Let us consider **<sup>W</sup>**<sup>∗</sup>+**W**<sup>+</sup> and a singularity of *f* at (**<sup>x</sup>**, *ξ*). Under Equations (1)–(3), half of this singularity propagates into the direction *ξ* (corresponding to the function *p*+). At the moment *t* = *T*, it is transformed to a singularity of **W**+(*f*) = *p*+(*T*) at (**<sup>y</sup>**+, *η*+). Under the adjoint Equation (15), half of this singularity propagates back to (**<sup>x</sup>**, *ξ*) at *t* = 0 to generate a singularity of **<sup>W</sup>**<sup>∗</sup>+**<sup>W</sup>**+(*f*). It is natural to believe that this recovered singularity is 1/4 of the original singularity of *f* (due to twice splitting, as described). The proof below verifies this intuition.

Indeed, denote by *q*+ the solution of the time-reversed wave equation, e.g., Equation (15), with the initial condition given by *g*+ = **<sup>W</sup>**+(*f*). Then, by definition (see Theorem 2) **<sup>W</sup>**<sup>∗</sup>*Tg*+ = *q*+(· , <sup>0</sup>)|*<sup>B</sup>*. The solution *q*+ can be decomposed into the sum *q*+ = *q*0 + *q*1. Here, *q*0, *q*1, up to smooth terms, are solutions of the wave equations in R<sup>3</sup> × (0, *T*) and satisfy *q*0(· , 0) = **<sup>W</sup>**<sup>∗</sup>+(*g*+), *q*1(· , 0) = **<sup>W</sup>**∗−(*g*+). We are only concerned with *q*0 since it defines **<sup>W</sup>**<sup>∗</sup>+**W**<sup>+</sup> *f* = *q*0(· , <sup>0</sup>). We can write

$$q\_0(\mathbf{y}, t) = \frac{1}{(2\pi t)^3} \int\_{\mathbb{R}^3} b(\mathbf{y}, t, \boldsymbol{\xi}) e^{i\boldsymbol{\Phi}\_+(\mathbf{y}, t, \boldsymbol{\xi})} \widehat{f}(\boldsymbol{\xi}) \, d\boldsymbol{\xi}.\tag{14}$$

Let *b*0 be the principal part of *b*. Then, the principal symbol *θ*+ of **<sup>W</sup>**<sup>∗</sup>+**W**<sup>+</sup> is given by *<sup>θ</sup>*+(**<sup>x</sup>**, *ξ*) = *b*0(**<sup>x</sup>**, 0, *ξ*). We note that *b*0 satisfies the same equation as *<sup>a</sup>*0,+ (see (12)). Therefore, on each bicharacteristic curve the ratio *b*0/*<sup>a</sup>*0,+ is constant which implies *b*0(**<sup>x</sup>**, 0, *ξ*) =

*<sup>a</sup>*+,<sup>0</sup>(**<sup>x</sup>**, 0, *ξ*)*b*0(**<sup>y</sup>**+, *T*, *η*+)/*<sup>a</sup>*+,<sup>0</sup>(**<sup>y</sup>**+, *T*, *η*+). Similar to the argumen<sup>t</sup> below Equation (13), up to lower order terms, we have

$$q\_0(\mathbf{y}\_+,T) = \frac{1}{2}\mathbf{g}\_+(\mathbf{y}\_+,T) = \frac{1}{(2\pi)^3} \int\_{\mathbb{R}^3} \frac{1}{2} a\_+(\mathbf{y}\_+,T,\xi) e^{i\boldsymbol{\phi}\_+\left(\mathbf{y}\_+,T,\xi\right)} f(\boldsymbol{\xi}) \,d\xi$$

.

This and Equation (14) implies that *b*0(**<sup>y</sup>**+, *T*, *ξ*) = *<sup>a</sup>*+,<sup>0</sup>(**<sup>y</sup>**+, *T*, *ξ*)/2. Therefore, we obtain *b*0(**<sup>x</sup>**, 0, *ξ*) = *<sup>a</sup>*+,<sup>0</sup>(**<sup>x</sup>**, 0, *ξ*)/2 = 1/4. Combining with a similar argumen<sup>t</sup> for **W**∗−**W**−, we obtain that the principal symbol of **<sup>W</sup>**<sup>∗</sup>*T***W***<sup>T</sup>* is *<sup>θ</sup>*(**<sup>x</sup>**, *ξ*) = 1/2. That is, **<sup>W</sup>**<sup>∗</sup>*T***W***<sup>T</sup>* = **I**/2 + **K***T*, where **K***T* is a pseudodifferential operator of order at most −1 and **I** is the identity. We have (**<sup>W</sup>***T f* , **W***T f*)=(**W**<sup>∗</sup>*T***W***<sup>T</sup> f* , *f*)=(*f* , *f*)/2 + (**<sup>K</sup>***T f* , *f*) and therefore conclude *f* 2*L*2 ≤ <sup>2</sup>(**<sup>W</sup>***T f* 2*L*2 + **<sup>K</sup>***T f* 2*L*2 ).

We are now ready to prove Theorem 1.

**Proof of Theorem 1.** Recall from Lemma 2 that *f L*<sup>2</sup>(*B*(0,*R*)) ≤ <sup>2</sup>(**<sup>W</sup>***T f L*<sup>2</sup>(R<sup>3</sup>) + **<sup>K</sup>***T f L*<sup>2</sup>(*B*(0,*R*))) where **K***T* is a pseudo-differential operator of order at most −1. Since **K***T* is compact and **W***T* is injective, applying ([51] Theorem V.3.1), we obtain *f L*<sup>2</sup>(*B*(0,*R*)) ≤ *<sup>C</sup>***<sup>W</sup>***T f L*<sup>2</sup>(R<sup>3</sup>) for some constant *C* ∈ (0, <sup>∞</sup>). This finishes our proof.

#### *3.3. Continuous Adjoint Operator*

Iterative methods for solving (9) require knowledge of the adjoint operator of **W***T*. In this subsection, we compute the continuous adjoint operator **<sup>W</sup>**<sup>∗</sup>*T* : *L*<sup>2</sup>(R<sup>3</sup>) → *L*<sup>2</sup>(*B*(0, *R*)) and prove that it is again given by the solution of a wave equation. More precisely, we have the following result.

**Theorem 2.** *Let g* ∈ *<sup>C</sup>*<sup>∞</sup>0(R<sup>3</sup>)*, consider the time reversed final state problem for the wave equation,*

$$\begin{aligned} q\_{tt}(\mathbf{x},t) - c^2(\mathbf{x})\Delta q(\mathbf{x},t) &= 0, \qquad (\mathbf{x},t) \in \mathbb{R}^3 \times (-\infty,T) \\ q(\mathbf{x},T) &= \mathbf{g}(\mathbf{x}), \quad \mathbf{x} \in \mathbb{R}^3 \\ q\_t(\mathbf{x},T) &= 0 \qquad \mathbf{x} \in \mathbb{R}^3 \end{aligned} \tag{15}$$

*and let χB*(0,*R*) *denote the indicator function of B*(0, *<sup>R</sup>*)*. Then,* **<sup>W</sup>**<sup>∗</sup>*Tg* = *χB*(0,*R*) *q*(· , <sup>0</sup>)*.*

**Proof.** It is clearly sufficient to show **<sup>W</sup>**<sup>∗</sup>*Tg* = *<sup>χ</sup>B*(0,*R*)*ut*(· , <sup>0</sup>), where *u* solves the wave equation *utt* − *c*2Δ*c* = 0 on R<sup>3</sup> × (−∞, *<sup>T</sup>*), with the final state conditions (*u*(· , *<sup>T</sup>*), *ut*(· , *T*)) = (0, *f*). Using the weak formulation (similar to [38]) for the wave equation, shows that *T*0 R<sup>3</sup> *<sup>c</sup>*<sup>−</sup><sup>2</sup>(**x**)*utt*(**<sup>x</sup>**, *<sup>t</sup>*)*v*(**<sup>x</sup>**, *t*) *d***x***dt* + *T*0 R<sup>3</sup> ∇*u*(**<sup>x</sup>**, *t*) · ∇*v* (**<sup>x</sup>**, *t*) *d***x***dt* = 0 for *v* ∈ *<sup>C</sup>*<sup>∞</sup>0 (R<sup>3</sup>). Two times integration by parts, rearranging terms and using the final state conditions for *u* yields R<sup>3</sup> *c*<sup>−</sup><sup>2</sup>(**x**)[ *f*(**x**)*v*(**<sup>x</sup>**, *T*) − *ut*(**<sup>x</sup>**, <sup>0</sup>)*v*(**<sup>x</sup>**, 0) + *<sup>u</sup>*(**<sup>x</sup>**, <sup>0</sup>)*vt*(**<sup>x</sup>**, 0)] *d***x** = *T*0 R<sup>3</sup> *<sup>u</sup>*(**<sup>x</sup>**, *t*) *<sup>c</sup>*<sup>−</sup><sup>2</sup>(**x**)*vtt*(**<sup>x</sup>**, *t*) − <sup>Δ</sup>*v*(**<sup>x</sup>**, *t*) *d***x***dt*. By taking *v* as the solution of (1)–(3) this yields R<sup>3</sup> *c*<sup>−</sup><sup>2</sup>(**x**)*g*(**x**) **<sup>W</sup>***T*(*f*)(**x**) *d***x** = R<sup>3</sup> *<sup>c</sup>*<sup>−</sup><sup>2</sup>(**x**)*ut*(**<sup>x</sup>**, <sup>0</sup>)*f*(**x**) *d***<sup>x</sup>**, which implies **<sup>W</sup>**<sup>∗</sup>*Tg* = *<sup>χ</sup>B*(0,*R*)*ut*(· , 0) = *<sup>χ</sup>B*(0,*R*)*q*(· , 0) and completes the proof.

We can reformulate the adjoint operator as follows.

**Corollary 1.** *For g* ∈ *<sup>C</sup>*<sup>∞</sup>0 (R<sup>3</sup>)*, let q be the solution of*

$$
\rho\_{ll}(\mathbf{x},t) - c^2(\mathbf{x})\Delta q(\mathbf{x},t) = 0 \quad \text{for } (\mathbf{x},t) \in \mathbb{R}^3 \times (0,\infty) \tag{16}
$$

*with initial conditions* (*q*(· , <sup>0</sup>), *qt*(· , <sup>0</sup>)=(*g*, <sup>0</sup>)*. Then,* **<sup>W</sup>**<sup>∗</sup>*Tg* = *χB*(0,*R*) *q*(· , *<sup>T</sup>*)*.*

**Proof.** Clearly *q* solves (16) with initial data (*q*(· , <sup>0</sup>), *qt*(· , <sup>0</sup>)=(*g*, 0) if and only if (*<sup>x</sup>*, *t*) → *q*(*<sup>x</sup>*, *T* − *t*) solves (15). Therefore, the claim follows from Theorem 2.

#### *3.4. Application of the Steepest 6hhod*

We propose solving the finite time wave inversion problem (9) via gradient type methods applied to residual functional Φ(*f*) := 12 **<sup>W</sup>***T f* − *h*2. Using standard PAT data, various iterative methods for PAT accounting for variable sound speed have been investigated in [36–40]. In particular, in [40] the steepest descent method has been demonstrated to be numerically efficient and robust. For the results shown in this paper we will therefore use the steepest-descent method. Our numerical experiments confirm that the steepest-descent method is also efficient for FFD-PAT and the final time wave inversion problem, where it reads as follows.

Because of the injectivity and boundedness of **W***T*, the sequence of iterates generated by Algorithm 1 converges to the unique solution of (9). The stability result (11) even implies that the steepest descent method converges linearly for FFD-PAT. More precisely, the sequence (*fk*)*<sup>k</sup>*∈<sup>N</sup> generated by Algorithm 1 satisfies the estimate *fk*+<sup>1</sup> − *f* 2 ≤ *ck f*0 − *f* 22 for some constant *c* ∈ (0, 1).


## **4. Numerical Simulations**

In this section, we numerical present results using Algorithm 1 for FFD-PAT with variable sound speed. Recall that for constant speed of sound, reconstructions from experimental FFD-PAT data are shown in [30,31]. In the present study, we give numerical results for spatially varying speed of sound. Performing experimental measurements for samples with variable speed of sound and applying our algorithms to such data is subject of future research.
