Thermoacoustic Tomography with Circular Integrating Detectors and Variable Wave Speed

TL;DR

This paper investigates the mathematical properties of thermoacoustic tomography with circular detectors under variable wave speed, analyzing how singularities in data relate and determining visibility conditions, supported by numerical experiments.

Contribution

It demonstrates that the measurement operator is a Fourier Integral Operator and characterizes the visibility of singularities in variable wave speed scenarios.

Findings

Measurement operator is a Fourier Integral Operator.

Visibility of singularities depends on the canonical relation.

Numerical results validate theoretical analysis.

Abstract

We explore Thermoacoustic Tomography with circular integrating detectors assuming variable, smooth wave speed. We show that the measurement operator in this case is a Fourier Integral Operator and examine how the singularities in initial data and measured data are related through the canonical relation of this operator. We prove which of those singularities in the initial data are visible from a fixed open subset of the set on which measurements are taken. In addition, numerical results are shown for both full and partial data.

Equations

$$H_D(\Omega )=\left\{f\in {\mathcal{E}}^{\prime }(\Omega )\mid {\int }_{\Omega }|\nabla f{|}^{2}dx<\infty \right\},$$

$$\{\begin{array}{ll}{\partial }_t^{2}u(t,x)=c^2(x)\Delta u(t,x),& (t,x)\in \mathbb{R}\times {\mathbb{R}}^2\\ u(0,x)=f(x),& f(x)\in H_D\left(\Omega \right)\\ {\partial }_{t}u(t,x){\mid }_{t=0}=0,& x\in {\mathbb{R}}^2\end{array}$$

$$u(t,x)=\frac{1}{(2\pi {)}^2}\sum _{\sigma =\pm }\int e^{i{\varphi }_{\sigma }(t,x,\xi )}{a}_{\sigma }(t,x,\xi )\hat{f}(\xi )d\xi .$$

$$\left(({\partial }_t{\varphi }_{\sigma }){\partial }_t-c^2(x)({\nabla }_x{\varphi }_{\sigma })\cdot {\nabla }_x+\frac{1}{2}({\partial }_t^2-c^2(x){\Delta }_x){\varphi }_{\sigma }\right)a_{\sigma }^{(0)}=0$$

$$Mf(t,\theta )=\frac{1}{2\pi }{\int }_0^{2\pi }u(t,C(\theta ,\alpha ))d\alpha ,$$

$$Mf(t,\theta )=\frac{1}{2\pi }{\int }_{{\mathbb{R}}^2}u(t,x)\delta ({\left|x-R\mathbf{\theta }\right|}^2-r^2)dx.$$

$$M_{+}f=\frac{1}{(2\pi {)}^3}{\iint }_{{\mathbb{R}}^2\times {\mathbb{R}}^2}{e}^{i{\varphi }_{+}(t,x,\xi )}{a}_{+}(t,x,\xi )\hat{f}(\xi )d\xi \delta (|x-R\mathbf{\theta }{|}^2-r^2)dx$$

$$M_{-}f=\frac{1}{(2\pi {)}^3}{\iint }_{{\mathbb{R}}^2\times {\mathbb{R}}^2}{e}^{i{\varphi }_{-}(t,x,\xi )}{a}_{-}(t,x,\xi )\hat{f}(\xi )d\xi \delta (|x-R\mathbf{\theta }{|}^2-r^2)dx.$$

$$M_{+}f(t,\theta )=\frac{1}{(2\pi {)}^3}{\iint }_{{\mathbb{R}}^2\times {\mathbb{R}}^2}{e}^{i\varphi (t,x,\xi )}a(t,x,\xi )\hat{f}(\xi )d\xi \delta ({\left|x-R\mathbf{\theta }\right|}^2-r^2)dx.$$

$$M_{+}f(t,\theta )=\frac{1}{(2\pi {)}^4}{\iiint }_{{\mathbb{R}}^2\times \mathbb{R}\times {\mathbb{R}}^2}{e}^{i\varphi (t,x,\xi )+i\lambda ({|x-R\mathbf{\theta }|}^2-r^2)}a(t,x,\xi )\hat{f}(\xi )d\xi d\lambda dx.$$

$$M_{+}f(t,\theta )=\frac{1}{(2\pi {)}^4}{\text{iiiint}}_{{\mathbb{R}}^7}{e}^{i\varphi (t,x,\xi )+i\lambda ({|x-R\mathbf{\theta }|}^2-r^2)-iy\cdot \xi }a(t,x,\xi )f(y)dyd\xi d\lambda dx,$$

$$\Phi (t,\theta ,y;\lambda ,x,\xi )=\varphi (t,x,\xi )+\lambda ({\left|x-R\mathbf{\theta }\right|}^2-r^2)-y\cdot \xi .$$

$$M_{+}f=\frac{1}{(2\pi {)}^4}{\text{iiiint}}_{{\mathbb{R}}^7}{e}^{\tilde{\Phi }(t,\theta ,y;\lambda ,\tilde{x},\xi )}\tilde{a}(t,\tilde{x},\xi ,\lambda )f(y)dyd\xi d\lambda d\tilde{x}$$

$$\tilde{a}(t,\tilde{x},\xi ,\lambda )=a\left(t,\frac{\tilde{x}}{|(\xi ,\lambda )|},\xi \right)|(\xi ,\lambda ){|}^{-2}.$$

$${\tilde{\Phi }}_{\lambda }$$

$${\tilde{\Phi }}_{\tilde{x}}$$

$${\tilde{\Phi }}_{\xi }$$

$${\left|x-R\mathbf{\theta }\right|}^2-r^2$$

$${\varphi }_x+2\lambda (x-R\mathbf{\theta })$$

$${\varphi }_{\xi }-y$$

$${\Sigma }_{+}\ni (t,\mathbf{\theta },y;\lambda ,x,\xi )\mapsto (t,\mathbf{\theta },y;{\Phi }_t,{\Phi }_{\theta },{\Phi }_y)$$

$${\Sigma }_{+}\ni (t,\mathbf{\theta },y;\lambda ,x,\xi )\mapsto (t,\mathbf{\theta },y;-c(y)|\xi |,-\frac{R}{r}c(y)|\xi |(x-(x\cdot \mathbf{\theta })\mathbf{\theta }),-\xi )$$

$${\Sigma }_{\pm }\ni (t,\mathbf{\theta },y;\lambda ,x,\xi )\mapsto (t,\mathbf{\theta },y;\mp c(y)|\xi |,\mp \frac{R}{r}c(y)|\xi |(x-(x\cdot \mathbf{\theta })\mathbf{\theta }),-\xi )$$

$$C_{\pm }:=\left\{(t,\mathbf{\theta },{\Phi }_t,{\Phi }_{\theta };y,\xi )\mid (t,\mathbf{\theta },y;\lambda ,x,\xi )\in {\Sigma }_{\pm }\right\}$$

$$=\left\{\left(t_{\pm ,i}(y,\xi ),{\mathbf{\theta }}_{\pm }(y,\xi ),\mp c(y)|\xi |,\mp \frac{R}{r}c(y)|\xi |({\gamma }_{y,\hat{\xi }}(t_{\pm ,i}(y,\xi ))-(x_{\pm ,i}\cdot {\mathbf{\theta }}_{\pm }){\mathbf{\theta }}_{\pm }(y,\xi ));y,\xi \right)\right\}$$

$$C_{\pm }:(y,\xi )\mapsto \{\begin{array}{c}\left(t_{\pm ,1},\frac{1}{R}{\gamma }_{y,\hat{\xi }}(t_{\pm ,1}\pm r),\mp c(y)|\xi |,\pm \frac{R}{r}c(y)|\xi |\left(\left(x_{\pm ,1}\cdot {\theta }_{\pm }\right){\gamma }_{y,\hat{\xi }}(t_{\pm ,1}\pm r)-{\gamma }_{y,\xi }(t_{\pm ,1})\right)\right)\\ \\ \left(t_{\pm ,2},\frac{1}{R}{\gamma }_{y,\hat{\xi }}(t_{\pm ,2}\mp r),\mp c(y)|\xi |,\mp \frac{R}{r}c(y)|\xi |\left(\left(x_{\pm ,2}\cdot {\theta }_{\pm }\right){\gamma }_{y,\hat{\xi }}(t_{\pm ,2}\mp r)-{\gamma }_{y,\xi }(t_{\pm ,2})\right)\right)\end{array}$$

$$C=C_{+}\cup C_{-}$$

$$C_{\pm }=\left\{\left(t_{\pm ,i},{\mathbf{\theta }}_{\pm },\mp c(y)|\xi |,\pm (-1{)}^i\frac{R}{r}c(y)|\xi |\left({\gamma }_{y,\hat{\xi }}(t_{\pm ,i})-(x_{\pm ,i}(y,\xi )\cdot {\mathbf{\theta }}_{\pm }){\mathbf{\theta }}_{\pm }\right);y,\xi \right)\right\},$$

$$C_{\pm }=\left\{((t_{\pm },{\mathbf{\theta }}_{\pm }(y,\xi ),\mp c(y)|\xi |,\mp \frac{1}{r}c(y)|\xi |({\gamma }_{y,\hat{\xi }}(t_{\pm })-(x_{\pm }\cdot {\mathbf{\theta }}_{\pm }){\mathbf{\theta }}_{\pm }(y,\xi ));y,\xi )\mid (y,\xi )\in T^{\ast }(\Omega )\setminus 0\right\}.$$

$$C=C_{+}\cup C_{-}$$

$$C_{\pm }=\left\{((t_{\pm },{\mathbf{\theta }}_{\pm }(y,\xi ),\mp c(y)|\xi |,\mp \frac{1}{r}c(y)|\xi |({\gamma }_{y,\hat{\xi }}(t_{\pm })-(x_{\pm }\cdot {\mathbf{\theta }}_{\pm }){\mathbf{\theta }}_{\pm }(y,\xi ));y,\xi )\mid (y,\xi )\in T^{\ast }(\Omega )\setminus 0\right\},$$

$$\omega (y_1,{\xi }_1)-\omega (y_2,{\xi }_2)=x_1-(x_1\cdot \mathbf{\theta })\mathbf{\theta }-(x_2-(x_2\cdot \mathbf{\theta })\mathbf{\theta })$$

$$=(x_1-x_2)-[(x_1-x_2)\cdot \mathbf{\theta }]\mathbf{\theta }$$

Full text

Author partly supported by NSF Grant DMS-1600327

1. Introduction

Thermoacoustic Tomography is a medical imaging method in which a short pulse of electromagnetic radiation is used to excite cells in some object we wish to image, typically the organs of a patient. Upon absorbing the EM radiation, the cells in the patient in turn vibrate, creating ultrasonic waves that then propagate out of the patient and are measured by any number of methods. Using this measured data, we then try to reconstruct, in some sense, an image of the inside of the patient. This is a hybrid imaging method which uses high contrast, low resolution EM radiation to excite the cells, and low contrast, high resolution ultrasound waves as measurement [7]. The hope is to be able to get an image with good contrast and resolution by combining these two types of waves. The case of point-wise measurements with constant and variable wave speed in the region of interest has been studied extensively [6, 9, 12]. Other methods of measurement of the ultrasonic waves include measurements with linear integrating detectors [5], planar integrating detectors [2, 10] and circular integrating detectors or cylindrical stacks of circular integrating detectors [13, 14]. Circular integrating detectors have a few advantages over linear integrating detectors and planar integrating detectors, including compactness of the experimental setup [13]. The case of planar integrating detectors was studied in [10], and that work focused on the problem with a smooth, variable wave speed. The case of circular (and cylindrical) integrating detectors with constant wave speed has been studied in [13, 14]. In those works, explicit formulae are given for reconstruction of an initial pressure density using full measurements, i.e. measurements for every circular integrating detector of a fixed radius with center on the unit circle, for all time. That reconstruction is stable in the case that the object being imaged is contained in the interior of the circular integrating detectors, but is unstable for the case that the object lies entirely outside of the circular integrating detector. The present work focuses on the case of circular integrating detectors in the plane with a 2 dimensional region of interest. Further, we do not make a constant wave speed assumption, we only assume that the wave speed $c(x)>0$ is smooth in all of $\mathbb{R}^{2}$ and is $1$ outside of a compact subset of $\mathbb{R}^{2}$. We show that the measurement operator in this case is a Fourier Integral Operator and compute its canonical relation, which happens to be a local diffeomorphism, thus allowing us to determine how singularities in initial data propagate to the measurement data. We also show that this operator is injective and prove stability of the measurement operator, and in addition we prove what singularities in the initial data are visible from a fixed open subset of the set of points on the circle where the measurements are taken in a given time interval. Lastly, we provide numerical results obtained through simulation in Matlab using both full and partial data that support our findings.

2. Setup

We begin by defining the space of distributions that our initial pressure distribution must be in. Let

[TABLE]

where $\Omega\subset\mathbb{R}^{2}$ is open and $\mathcal{E}^{\prime}(\Omega)$ is the space of distributions compactly supported in $\Omega$. This is the natural space in which to take $f$ when the energy of the system is taken into consideration. Let $\lVert f\rVert_{H_{D}(\Omega)}^{2}=\int_{\Omega}|\nabla f|^{2}\,dx$. The space $H_{D}(\Omega)$ is the completion of $C_{0}^{\infty}(\Omega)$ under the given norm. We further suppose that $\bar{\Omega}\subset B_{1}(0)$ where $B_{1}(0)$ is the unit ball in $\mathbb{R}^{2}$ centered at the origin. We view $f$ as an initial pressure distribution of some object to be imaged represented by $\Omega$. Then, after exposing $\Omega$ to microwave radiation, the ultrasonic waves created solve the acoustic wave equation given by

[TABLE]

where $0<c(x)\in C^{\infty}(\mathbb{R}^{2})$ is the smooth wave speed, assumed to be known. Outside of $\Omega$, $c\equiv 1$. The problem of interest is to detect these waves, solutions $u(t,x)$ to the above wave equation, with detectors located on the boundary of $B_{R}(0)$, and then using these measurements, reconstruct the initial pressure distribution $f$. As mentioned in the introduction, extensive research has been done in the constant speed case ($c(x)\equiv 1$ for all $x\in\mathbb{R}^{2}$) and variable speed case with point detectors in which we assume access to $u\mid_{U\times\Gamma}$ where $\Gamma\subset S^{1}$ is open and $U\subset\mathbb{R}$ is some time interval. Research has also been done for linear and planar integrating detectors in both the constant and variable speed case, and also circular and cylindrical integrating detectors in the constant wave speed case. When imaging with these integrating detectors, instead of assuming direct access to $u$ on some open subset of the boundary of $\Omega$, the measured data is an average of $u$ over a circular detector of radius $r$ centered on the boundary of the ball of radius $R$ (which we will choose later), and data is assumed to be collected on an open subset of of this boundary, not necessarily the entirety of the boundary. The present work considers the problem with variable speed in $\Omega$, and circular integrating detectors. We will have two cases to consider, which we will call the large radius detector case and the small radius detector case, which are depicted in Figure 1. The large radius detector case is the experimental setup in which $\Omega$ is on the “inside” of the circular integrating detectors, and the small radius detector case is the setup in which $\Omega$ is on the “outside” of the circular integrating detector.

3. Construction

We are interested in seeing what singularities we can recover from rotating a circular integrating detector around some object that has been heated via microwaves. To start, we recall that solving the wave equation (as above) up to a smooth error, for $x\in\mathbb{R}^{2}$, can be accomplished by use of the geometric optics construction (see Section 3 of [11]). After making this construction, we have that $u$ up to a smooth error is given locally by

[TABLE]

Here, $\phi_{\sigma}$ solves the eikonal equation: $((\phi_{\sigma})_{t})^{2}=c^{2}(x)|\nabla_{x}\phi_{\sigma}|^{2}$ with initial condition $\phi_{\sigma}(0,x,\xi)=x\cdot\xi$. In addition, $\phi_{\sigma}$ is positive homogeneous of order $1$ in $\xi$, i.e. $\phi_{\sigma}(t,x,\lambda\xi)=\lambda\phi_{\sigma}(t,x,\xi)$ for $\lambda>0$. The eikonal equation is only solvable locally in time, which results in our solution $u(t,x)$ being only a local solution in time. This is not actually an issue however, as we simply repeat this procedure with new initial conditions. The resulting Fourier Integral Operator is then actually a composition of Fourier Integral Operators, which is a technical issue we refer the reader to [9] for more details. Because of this, we may assume that the eikonal equation is solvable until geodesics intersect circular integrating detectors. By this construction, we may obtain $a_{\sigma}$ as a classical symbol of order [math] with expansion $a_{\sigma}(t,x,\xi)\sim\sum_{j\geq 0}a_{\sigma}^{(j)}(t,x,\xi)$ where $a_{\sigma}^{(j)}(t,x,\xi)$ is positively homogeneous of order $-j$ for $|\xi|$ large. In particular, $a_{\sigma}^{(0)}(t,x,\xi)$ solves the transport equation

[TABLE]

with initial conditions $a_{\sigma}^{(0)}(0,x,\xi)=1/2$. The last term on the left hand side of this equation acts on $a_{\sigma}^{(0)}$ by multiplication.

Now in the situation of Thermoacoustic Tomography using circular integrating detectors around the object we wish to image, the measurements at the detector are given by the circular Radon transform:

[TABLE]

where $C(\theta,\alpha)=R(\cos\theta,\sin\theta)+r(\cos\alpha,\sin\alpha)$ is a parametrization of the circular detector, $R$ is the distance from the origin to the center of the circle, $r$ is the radius of the circular detector, $\theta$ is the angle made between the positive horizontal axis and the ray from the origin to the center of the circle, and $u(t,x)$ is the solution to the IVP (1). The radius of the circular integrating detector, $r$, must be chosen so that the detector does not intersect $\Omega$. To accomplish this, we must have either $r$ small enough so that $R-r\geq 1$, guaranteeing the detector does not intersect $\Omega\subset B_{1}(0)$, or, we could fix $R=1$ and choose $r\geq 2$, in which case $\Omega$ is contained in the interior of the disc defined by the detector (see Figure 1). For convenience, we define $\boldsymbol{\theta}=(\cos\theta,\sin\theta)$. We can rewrite $Mf$ by using the $\delta$ distribution:

[TABLE]

We now plug in our solution for $u$ obtained via the geometric optics construction and denote by $M_{+}$ and $M_{-}$ the operators taking $\sigma=+$ and $\sigma=-$ respectively after substituting the geometric optics solution in $Mf(t,\theta)$. Then $Mf=M_{+}f+M_{-}f$, where

[TABLE]

We drop subscripts in the integral for now and consider only $M_{+}$,

[TABLE]

We make use of the fact that $\frac{1}{2\pi}\int e^{i\lambda(\left|x-R\boldsymbol{\theta}\right|^{2}-r^{2})}d\lambda=\delta(\left|x-R\boldsymbol{\theta}\right|^{2}-r^{2})$ to say

[TABLE]

Lastly, we unpack the Fourier transform of $f$ to get

[TABLE]

where we have identified $\mathbb{R}^{2}\times\mathbb{R}^{2}\times\mathbb{R}\times\mathbb{R}^{2}$ with $\mathbb{R}^{7}$. This is an indication that the measurement operator $M_{+}$ is a Fourier Integral Operator with phase function

[TABLE]

It can be shown that this phase function is non-degenerate in the sense of [11].

One issue with this phase function is that $\Phi$ is not homogeneous of degree one in the fiber variables $(\lambda,x,\xi)$, but this can be fixed by making a change of variable. Let $\tilde{x}:=x\left|(\xi,\lambda)\right|$ and define $\tilde{\Phi}(t,\theta,y;\lambda,\tilde{x},\xi):=\Phi(t,\theta,y;\lambda,\frac{\tilde{x}}{\left|(\xi,\lambda)\right|},\xi)$. This makes $\tilde{\Phi}$ homogeneous of degree one in the variables $(\lambda,\tilde{x},\xi)$, and so we can proceed. After making this change of variable, we now write $M_{+}$ as

[TABLE]

where

[TABLE]

Also, because $a$ is an amplitude of order [math], by [3] we know that $\tilde{a}(t,\tilde{x},\lambda,\xi)$ is an amplitude of order $-1/2$. Note that this change of variable does not affect the characteristic manifold for $\Phi$, for

[TABLE]

so that $\tilde{\Phi}_{\tilde{x}}=0$ if and only if $\Phi_{x}=0$, from which it is clear that $\tilde{\Phi}_{\lambda,\tilde{x},\xi}=0$ if and only if $\Phi_{\lambda,x,\xi}=0$.

The characteristic manifold of this FIO is defined as the set $\Sigma_{+}=\left\{(t,\theta,y;\lambda,x,\xi)\mid\Phi_{\lambda,x,\xi}=0\right\}$. Taking the derivative of $\Phi$ with respect to $(\lambda,x,\xi)$, we see that this gives the system of equations

[TABLE]

We’ve made no assumption on which experimental setup we’ve chosen to examine so far. There are two different cases, as mentioned above: (1)$R-r\geq 1$, which we will call the small radius case, and (2)$R=1$ with $r\geq 2$, which we will call the large radius case. The analysis of these two cases are largely the same, but with a few key differences. We examine both cases.

3.1. Case 1: Small Radius

From the system of equations obtained by looking at the characteristic manifold of the FIO, we see that $\phi_{\xi}=y$, and by the geometric optics construction, $x$ lies on the geodesic $\gamma_{y,\hat{\xi}}(t)$ issued from $(y,\hat{\xi})$ where $\hat{\xi}=\xi/(c|\xi|)$ is the unit covector in the metric identified with a unit vector, and $(\gamma_{y,\hat{\xi}}(t),c|\xi|\dot{\gamma}_{y,\hat{\xi}}(t))=(x,\phi_{x})$. Now from the first equation, we know that $x=(x_{1},x_{2})$ must lie on the circular integrating detector of radius $r$ with center $R\boldsymbol{\theta}$. So $x$ is the intersection of the geodesic $\gamma_{y,\hat{\xi}}$ with the circle defined by $|x-R\boldsymbol{\theta}|^{2}=r^{2}$. There are two of these points of intersection in general (we will show that the geodesic does not intersect the circular integrating detector tangentially if the singularity is to be detected), which we label $x_{+,1}(y,\xi)$ and $x_{+,2}(y,\xi)$. Also denote the times at which these intersections occur $t_{+,1}(y,\xi)>0$ and $t_{+,2}(y,\xi)>0$ respectively.

Finally, we have from the second equation $\phi_{x}=2\lambda(R\boldsymbol{\theta}-x)$. This tells us that $\dot{\gamma}_{y,\xi}(t)$ is parallel to $x-R\boldsymbol{\theta}$, provided that $\lambda\neq 0$ (in which case $\gamma_{y,\hat{\xi}}(t)$ intersects the circular integrating detector tangentially). Supposing for a moment that $\lambda=0$, then taking magnitudes on both sides of the equation $\phi_{x}(t,x,\xi)=2\lambda(R\boldsymbol{\theta}-x)$, we obtain $|\phi_{x}(t,x,\xi)|=0\implies c(y)|\xi||\dot{\gamma}_{y,\hat{\xi}}(t)|=0$. This in turn means $|\dot{\gamma}_{y,\hat{\xi}}(t)|=0$ as $c(y)\neq 0$ for any $y\in\mathbb{R}^{2}$ and $\xi$ is a non zero vector. But we know near the integrating detectors, that $|\dot{\gamma}_{y,\hat{\xi}}(t)|=1$, a contradiction, so $\lambda\neq 0$ on the characteristic manifold. And so the geodesic intersects the circular detector perpendicularly. We know that $\phi_{x}=c|\xi|\dot{\gamma}_{y,\hat{\xi}}(t)$, And $c\equiv 1$ outside of the region of interest, so that $|\lambda_{+,i}|=c(y)/(2r)|\xi||(\dot{\gamma}_{y,\hat{\xi}}(t_{+,i}))|=c(y)/(2r)|\xi|$. Note that $\lambda_{+,1}>0$ and $\lambda_{+,2}<0$, so that $\lambda_{+,1}=c(y)/(2r)|\xi|$ and $\lambda_{+,2}=-\lambda_{+,1}$. We’ll simply denote $\lambda_{+}:=\lambda_{+,1}$. Because $\gamma_{y,\hat{\xi}}(t)$ intersects the circular integrating detector perpendicularly, it must go through the center of the circular integrating detector, as outside of $\Omega$, we know that $c\equiv 1$ implies that $\gamma_{y,\hat{\xi}}$ is a straight line near the integrating detectors, and so we see that $\boldsymbol{\theta}=\gamma_{y,\hat{\xi}}(t_{+,1}(y,\xi)+r)/R=\gamma_{y,\hat{\xi}}(t_{+,2}(y,\xi)-r)/R$, where $r$ is the fixed radius of the circular integrating detector. This gives then the entire characteristic manifold parametrized by $(y,\xi)$, giving a smooth manifold of dimension 4 consisting of 2 connected parts. Define $\Sigma_{+,1}=\left\{(t_{+,1}(y,\xi),\boldsymbol{\theta}(y,\xi),y;\lambda_{+}(y,\xi),x_{+,1}(y,\xi),\xi)\mid(y,\xi)\in T^{*}(\Omega)\setminus\left\{0\right\}\right\}$ and $\Sigma_{+,2}=\left\{(t_{+,2}(y,\xi),\boldsymbol{\theta}(y,\xi),y;-\lambda_{+}(y,\xi),x_{+,2}(y,\xi),\xi)\mid(y,\xi)\in T^{*}(\Omega)\setminus\left\{0\right\}\right\}$. Then $\Sigma_{+}=\Sigma_{+,1}\cup\Sigma_{+,2}$ as a disjoint union.

We now write down the canonical relation given by

[TABLE]

This mapping can be calculated as

[TABLE]

The analysis for $\sigma=-$ is the same giving us $t_{-,i}(y,\xi)<0,x_{-,i}(y,\xi)$ for $i=1,2$ and $\boldsymbol{\theta}_{-}$. We also see that $\lambda_{-,1}<0$ and $\lambda_{-,2}>0$. Denote $\lambda_{-}:=\lambda_{-,1}$. We see that the maps

[TABLE]

are smooth and of full rank, and so the canonical relations associated to the operators $M_{+}$ and $M_{-}$ are given by

[TABLE]

[TABLE]

such that $(y,\xi)\in T^{*}(\Omega)\setminus{0}$, where $i=1,2$. Here, $\boldsymbol{\theta}_{+}=\gamma_{y,\hat{\xi}}(t_{+,1}(y,\xi)+r)/R$ and $\boldsymbol{\theta}_{-}=\gamma_{y,\hat{\xi}}(t_{-,1}(y,\xi)-r)/R$. Writing this as a mapping, we have

[TABLE]

Where $x_{\pm,i}=x_{\pm,i}(y,\xi)$, $\theta_{\pm}=\theta_{\pm}(y,\xi)$, and $t_{\pm,1}=t_{\pm,1}(y,\xi)$ and similarly, $t_{\pm,2}=t_{\pm,2}(y,\xi)$. The canonical relations for the operators $M_{+}$ and $M_{-}$ are each one to two and of the above form. The above has shown the following:

Theorem 1**.**

For $R-r\geq 1$, the operator $Mf=M_{+}f+M_{-}f$ defined above is a Fourier Integral Operator of order $-1/2$ with canonical relation given by

[TABLE]

where

[TABLE]

($t_{\pm,i}=t_{\pm,i}(y,\xi)$ and $\boldsymbol{\theta}_{\pm}=\boldsymbol{\theta}_{\pm}(y,\xi)$) with $(y,\xi)\in T^{*}(\Omega)\setminus{0}$ and $i=1,2$.

Note that this canonical relation is locally one to four.

3.2. Case 2: Large Radius

In this case, the analysis is almost entirely the same, except there is only one point of intersection of the geodesic $\gamma_{y,\hat{\xi}}(t)$ with the circular integrating detector defined by $|x-\boldsymbol{\theta}|=r$ (see Figure 1). We then have $\Sigma_{\pm}=\{(t_{\pm}(y,\xi),\boldsymbol{\theta}_{\pm}(y,\xi),y;\pm\lambda(y,\xi),x_{\pm}(y,\xi),\xi)\mid(y,\xi)\in T^{*}(\Omega)\setminus\{0\}\}$, and the canonical relations are given by

[TABLE]

Here $t_{\pm}=t_{\pm}(y,\xi)$. The canonical relations are each locally one to one in this case, and we have an analogous result as in the first case:

Theorem 2**.**

For $R=1$ and $r\geq 2$, the operator $Mf=M_{+}f+M_{-}f$ defined above is a Fourier Integral Operator of order $-1/2$ with canonical relation given by

[TABLE]

where

[TABLE]

with $t_{\pm}=t_{\pm}(y,\xi)$.

This canonical relation is locally one to two, as each individual canonical map is locally one to one.

Note that the maps $C_{+}$ and $C_{-}$ are not globally one to 1, although each are locally one to one, for suppose (looking only at $C_{+}$ for a moment) $(t(y_{1},\xi_{1}),\boldsymbol{\theta}(y_{1},\xi_{1}),\tau(y_{1},\xi_{1}),\omega(y_{1},\xi_{1}))=(t(y_{2},\xi_{2}),\boldsymbol{\theta}(y_{2},\xi_{2}),\tau(y_{2},\xi_{2}),\omega(y_{2},\xi_{2}))$ where $\tau(y,\xi)=-c(y)|\xi|$ and $\omega(y,\xi)=-\frac{1}{r}c(y)|\xi|(x-(x\cdot\boldsymbol{\theta})\boldsymbol{\theta}(y,\xi))$ and $x=\gamma_{y,\hat{\xi}}(t(y,\xi))$. We’ll call $x_{i}=x(y_{i},\xi_{i})$ for $i=1,2$. Then clearly we have $t(y_{1},\xi_{1})=t(y_{2},\xi_{2})$ and $\boldsymbol{\theta}(y_{1},\xi_{1})=\boldsymbol{\theta}(y_{2},\xi_{2})$, which we’ll just label $t$ and $\boldsymbol{\theta}$ respectively. We also clearly have $c(y_{1})|\xi_{1}|=c(y_{2})|\xi_{2}|$. Suppose for a moment that $(x_{1}-x_{2})\cdot\boldsymbol{\theta}^{\bot}=0$ where $\boldsymbol{\theta}^{\bot}=(-\sin(\theta),\cos(\theta))$ is the unit vector perpendicular to $\boldsymbol{\theta}$. Then note that

[TABLE]

Now note that $(\omega(y_{1},\xi_{1})-\omega(y_{2},\xi_{2}))\cdot\boldsymbol{\theta}=0$ and $(\omega(y_{1},\xi_{1})-\omega(y_{2},\xi_{2}))\cdot\boldsymbol{\theta}^{\bot}=0$. Because $\boldsymbol{\theta}$ and $\boldsymbol{\theta}^{\bot}$ are linearly independent, this shows that $\omega(y_{1},\xi_{1})-\omega(y_{2},\xi_{2})=0$ or in other words that $\omega(y_{1},\xi_{1})=\omega(y_{2},\xi_{2})$. This shows that, provided $t(y_{1},\xi_{1})=t(y_{2},\xi_{2})$, $\boldsymbol{\theta}(y_{1},\xi_{1})=\boldsymbol{\theta}(y_{2},\xi_{2})$, $c(y_{1})|\xi_{1}|=c(y_{2})|\xi_{2}|$ and $(x_{1}-x_{2})\cdot\boldsymbol{\theta}^{\bot}=0$, that $(y_{1},\xi_{1})$ and $(y_{2},\xi_{2})$ get mapped to the same point under the canonical relation. However, for $c$ close enough to $1$, this will not happen locally.

4. Injectivity

4.1. Case 1: Small Radius

Let $u(t,x)$ be the solution to (1) and $\Gamma\subset S^{1}$ be open. In local coordinates, suppose $\Gamma$ is the open interval given by $\Gamma=(\theta_{1},\theta_{2})$, with $0\leq\theta_{1}<\theta_{2}<2\pi$. Then, for $\boldsymbol{\theta}=(\cos(\theta),\sin(\theta))\in\Gamma$, we have

[TABLE]

Where $R_{0}$ is the distance from the orgin to the center of the circular integrating detector. For fixed $r>0$, we may view $R$ as variable with $R\geq R_{0}\geq 1+r$ (i.e. we may translate the circular integrating detectors away from the region of interest). We denote this by letting $Mf(t,\theta)$ vary with $R$ and denote the operator then as $Mf(t,\theta;R)$ and note $Mf(t,\theta)=Mf(t,\theta;R_{0})$. Let $P(t,\theta,R)=Mf(t,\theta;R)$. Then we see that

[TABLE]

This can be seen as follows: first note to save space, $u$ and all of its partial derivatives are understood to be evaluated at $(t,R\cos(\theta)+r\cos(\alpha),R\sin(\theta)+r\sin(\alpha))$. We have

[TABLE]

and so we have $P_{\theta\theta}$ can be written:

[TABLE]

Noting that

[TABLE]

we see that $P_{\theta\theta}(t,\theta,R)$ is given by

[TABLE]

Remembering that $\Delta u$ above is evaluated at $(t,R\cos(\theta)+r\cos(\alpha),R\sin(\theta)+r\sin(\alpha)),$ where $c\equiv 1$, we have $P_{\theta\theta}(t,\theta,R)=$

[TABLE]

so that

[TABLE]

[TABLE]

and

[TABLE]

So we only need to show that $R(RP_{R})_{R}=$

[TABLE]

This follows from direct calculation and the chain rule.

Rewriting (2) in a more standard form, we have:

[TABLE]

(3) simply says that $P(t,\theta,R)$ is a solution to the constant speed wave equation in polar coordinates with initial condition $P(0,\theta,R)=0$ for $R\geq R_{0}>1+r$ and $P(t,\theta,R_{0})=Mf(t,\theta)$ for $t\geq 0$ and $\theta\in\Gamma$. It is well known that the wave equation on a Riemannian manifold has a unique solution, and so the solution to (3) is unique. To show then that $Mf(t,\theta)\mid_{\theta\in\Gamma}$ is uniquely determined by $f$, we need to show that $Mf(t,\theta)\mid_{\theta\in\Gamma}=0\implies f(x)=0$ by the linearity of $M$. So then, we assume that $Mf(t,\theta)\equiv 0$ for $t\in\mathbb{R}^{+}$ and $\theta\in\Gamma$. By the uniqueness of solutions to the wave equation, this tells us that $P(t,\theta,R)\equiv 0$ for any $(t,\theta,R)\in\mathbb{R}^{+}\times\Gamma\times[R_{0},\infty)$. Let $T\in\mathbb{R}$. We may extend $u(t,x)$ in an even way for $|t|<T$ so that $u$ is still a solution to the wave equation, and so we may extend $P(t,\theta,R)$ in an even way such that $P(t,\theta,R)=0$ for $|t|<T$. By finite speed of propagation, we know that $\textnormal{supp }u(T,\cdot)\subset B_{1+T}(0)$. Let $\theta_{0}\in\Gamma$. Note that the set $\mathcal{A}=\{(\theta,R)\mid\theta\in\Gamma,\,R>R_{0}\}$ is open and connected in $\mathbb{R}^{2}$, because $\Gamma\subset S^{1}$ is an open interval. We know that the circular Radon transform of $u(T,\cdot)$ is 0 for any $\theta\in\Gamma$ and for any $R>R_{0}$. Further, because the interior of these circular integrating detectors is just $D(\theta,R)=B_{r}(R\boldsymbol{\theta})$, we can take $R>R_{0}$ large enough so that $\textnormal{supp }u(T,\cdot)\cap D(\theta_{0},R)=\emptyset$, because $\textnormal{supp }u(T,\cdot)$ is contained in a bounded set. It follows then by Theorem 1.2 in [8] that $\textnormal{supp }u(T,\cdot)$ is disjoint from $\cup_{(\theta,R)\in\mathcal{A}}D(\theta,R)$. So in particular, there is a neighborhood $V$ of $x_{0}=R_{0}\boldsymbol{\theta}_{0}$ such that $u(T,x)=0$ on $V$. $T$ was chosen arbitrarily, so this result holds for all $|t|<T$, and so by Tataru’s unique continuation, $u(t,x)\equiv 0$ in the domain of influence $|t|+\textnormal{dist}(x,x_{0})<T$. So, taking $T$ large enough so that $\textnormal{dist}(x,x_{0})<T$ for all $x\in\Omega$, we see that $u(0,x)=f(x)\equiv 0$, and so $Mf(t,\theta)\mid_{[0,T]\times\Gamma}$ is uniquely determined by $f$.

4.2. Case 2: Large Radius

Again, we consider $u(t,x)$ a solution to (1). We consider only the full data case $\Gamma=S^{1}$. Then in this case, in which $R=1$, we have that the measurement operator $Mf(t,\theta)$ is given by

[TABLE]

where $r_{0}\geq 2$ is the fixed radius of the circular integrating detector. We may however view $r$ as variable, noting the operator with variable detector radius by

[TABLE]

Where $r\geq r_{0}\geq 2$. Let $P(t,\theta,r)=Mf(t,\theta;r)$, where $\theta\in\Gamma$ and $r\geq r_{0}>2$. Note that $P(t,\theta,r_{0})=Mf(t,\theta)$ for $\theta\in\Gamma$. It follows then that $P(t,\theta,r)$ solves the following PDE:

[TABLE]

We see this as follows: (Again note that $u$ and all of its partial derivatives are understood to be evaluated at $(t,\cos(\theta)+r\cos(\alpha),\sin(\theta)+r\sin(\alpha))$.) We have

[TABLE]

We integrate by parts to get

[TABLE]

Then, we use the fact that

[TABLE]

to obtain

[TABLE]

We recall that $u$ and its derivatives are evaluated at $(t,\cos(\theta)+r\cos(\alpha),\sin(\theta)+r\sin(\alpha))$ where $c\equiv 1$, so that $\Delta u=u_{tt}$ there and we see that

[TABLE]

The last line can be seen by direct calculation of $P_{rr}(t,\theta,r)$. Rearranging we see $rP_{tt}=rP_{rr}+P_{r}$ and (4) then follows. This partial differential equation is the wave equation with axial symmetry, and so solutions to this equation are again unique.

To show that $Mf(t,\theta)$ determines $f$ uniquely for $\theta\in\Gamma$, by the linearity of $M$, we need only show that $Mf(t,\theta)=0$ for all $(t,\theta)=[0,T)\times\Gamma\implies f=0$. We note that $Mf(t,\theta)=0$ for all $(t,\theta)\in[0,T)\times\Gamma\implies P(t,\theta,r)\equiv 0$, because $P(t,\theta,r)\equiv 0$ clearly solves the PDE (4) in this case, and this solution is unique. For any $T\in\mathbb{R}$, we know that $\textnormal{supp }u(T,\cdot)\subset B_{1+T}(0)$, and so for any finite $T$, $u(T,\cdot)$ has bounded support. Let $\theta_{0}\in\Gamma$, $T\in\mathbb{R}^{+}$. Let $r_{1}=\sup\{r\geq r_{0}\mid C(\theta_{0},r)\cap\textnormal{supp }u(T,\cdot)\neq\emptyset\}$, where $C(\theta_{0},r)$ is the circle centered at $\boldsymbol{\theta}_{0}$ with radius $r$. We know, because $\textnormal{supp }u(T,\cdot)\subset B_{1+T}(0)$, that $\{r\geq r_{0}\mid C(\theta_{0},r)\cap\textnormal{supp }u(T,\cdot)\neq\emptyset\}$ is bounded above by $2+T$. Assume that this set is nonempty so that $r_{1}\geq r_{0}$ is finite. Then we have $C(\theta_{0},r_{1})\cap\textnormal{supp }u(T,\cdot)\neq\emptyset$ by the compactness of $\textnormal{supp }u(T,\cdot)$. Let $x\in C(\theta_{0},r_{1})\cap\textnormal{supp }u(T,\cdot)$, and let $(x,\xi)\in N^{*}(C(\theta_{0},r_{1}))$. By construction, $\textnormal{supp }u(T,\cdot)$ is on one side of $C(\theta_{0},r_{1})$ at $x$, so that by Theorem 8.5.6 in [4], we have that $(x,\xi)\in\textnormal{WF}_{A}(u(T,\cdot))$. Note that as in the terminology of [8], it is impossible for $x$ to be a $C(\theta_{0},r_{1})$ self mirror point, for tracing $(x,\xi)$ back along the geodesic defined by $(x,\xi)$, we see that the geodesic never intersects the interior of $B_{1}(0)$, which is impossible. So then, there are two cases we must consider. First we consider the case where the $C(\theta_{0},r_{1})$ mirror point of $x$, which we will call $\tilde{x}$, is not in the intersection $C(\theta_{0},r_{1})\cap\textnormal{supp }u(T,\cdot)$. Then, by the compactness of $\textnormal{supp }u(T,\cdot)$, we have $u(T,\cdot)=0$ in a neighborhood of $\tilde{x}$. We also have $P(T,\theta,r)=0$ in a neighborhood of $(\theta_{0},r_{1})\in\Gamma\times[r_{0},\infty)$. It then follows by Proposition 2.4 of [8] that $(x,\xi)\not\in\textnormal{WF}_{A}(u(T,\cdot))$, a contradiction to $(x,\xi)\in\textnormal{WF}_{A}(u(T,\cdot))$. It follows that $\{r\geq r_{0}\mid C(\theta_{0},r)\cap\textnormal{supp }u(T,\cdot)\neq\emptyset\}$ is empty and so $\textnormal{supp }u(T,\cdot)\subset\bar{B}_{r_{0}}(\boldsymbol{\theta}_{0})$.

The second case we consider is that $\tilde{x}$ is in the intersection $C(\theta_{0},r_{1})\cap\textnormal{supp }u(T,\cdot)$. We’ll show that there then exists $\theta_{1}\in S^{1}$ and $r_{2}>r_{1}$, such that $C(\theta_{1},r_{2})\cap\textnormal{supp }u(T,\cdot)=\{\tilde{x}\}$. Assume for now that this is the case, and let $(\tilde{x},\tilde{\xi})\in N^{*}(C(\theta_{1},r_{2}))$. Then again, we have by construction that $\textnormal{supp }u(T,\cdot)$ is on one side of $C(\theta_{1},r_{2})$ at $\tilde{x}$ and so $(\tilde{x},\tilde{\xi})\in\textnormal{WF}_{A}(u(T,\cdot))$. It also follows as before that $u(T,\cdot)$ is zero in a neighborhood of the $C(\theta_{1},r_{2})$ mirror point of $\tilde{x}$. We then see again from Proposition 2.4 of [8] that $(\tilde{x},\tilde{\xi})\not\in\textnormal{WF}_{A}(u(T,\cdot))$, a contradiction. It follows again that $\textnormal{supp }u(T,\cdot)\subset\bar{B}_{r_{0}}(\boldsymbol{\theta}_{0})$. That $f$ is zero then follows from Tataru’s unique continuation as in the small radius case.

Now we show the existence of the circle $C(\theta_{1},r_{2})$ with property that $C(\theta_{1},r_{2})\cap\textnormal{supp }u(T,\cdot)=\{\tilde{x}\}$ mentioned above. We let $C(\theta_{0},r_{1})$, and $x,\tilde{x}$ be as above. We define the sets $C^{-}(\theta_{0},r_{1})=\{y\in C(\theta_{0},r_{1})\mid(y-\boldsymbol{\theta_{0}})\cdot\boldsymbol{\theta_{0}}\leq 0\}$, and $C^{+}(\theta_{0},r_{1})=\{y\in C(\theta_{0},r_{1})\mid(y-\boldsymbol{\theta_{0}})\cdot\boldsymbol{\theta_{0}}\geq 0\}$. It is clear that $B_{1}(0)$ is contained in the interior of the region bounded by $C^{-}(\theta_{0},r_{1})\cup L(\theta_{0},r_{1})$, where $L(\theta_{0},r_{1})$ is the diameter of $C(\theta_{0},r_{1})$ defined by the vector $\boldsymbol{\theta}_{0}^{\bot}$. We may assume without loss of generality that $x\in C^{-}(\theta_{0},r_{1})$, for if not, we may simply swap the roles of $x$ and $\tilde{x}$ in what follows. Now $x\in C^{-}(\theta_{0},r_{1})\implies\tilde{x}\in C^{+}(\theta_{0},r_{1})$, because $x$ cannot be a $C(\theta_{0},r_{1})$ mirror point, as we have shown. The line $\ell(t)$ defined by

[TABLE]

intersects $S^{1}$ at 2 points: $\boldsymbol{\theta}_{0}$ and

[TABLE]

This implies there are two distinct circles with centers on $S^{1}$ such that $(\tilde{x},(\tilde{x}-\boldsymbol{\theta}_{0}))$ is in the conormal bundle to these circles, namely $C(\theta_{0},r_{1})$ and $C(\theta_{1},r_{2})$, where $r_{2}=|\tilde{x}-\boldsymbol{\theta}_{1}|$. Note that

[TABLE]

where the last inequality follows because

[TABLE]

and we see then that $r_{2}=r_{1}+\frac{2(\tilde{x}-\boldsymbol{\theta}_{0})\cdot\boldsymbol{\theta}_{0}}{r_{1}}$. Now, we need only show that $C(\theta_{1},r_{2})\cap\textnormal{supp }u(T,\cdot)=\{\tilde{x}\}$. Clearly by construction, $\tilde{x}$ is in this intersection. Also note that by the choice of $r_{1}$, that $C(\theta_{0},r)\cap\textnormal{supp }u(T,\cdot)=\emptyset$ for all $r>r_{1}$. So let $y\in C(\theta_{1},r_{2})$. Then, we have using the triangle inequality

[TABLE]

so that $|y-\boldsymbol{\theta}_{0}|\geq r_{1}$, and equality holds only when $y-\boldsymbol{\theta}_{0}=\alpha(\boldsymbol{\theta}_{1}-\boldsymbol{\theta}_{0})$ for some $\alpha>0$, but this is only the case when $y=\tilde{x}$. So, we’ve shown that $y\in C(\theta_{1},r_{2})\implies|y-\boldsymbol{\theta}_{0}|>r_{1}$ if $y\neq\tilde{x}$. In other words, $y\neq\tilde{x}\implies y\in C(\theta_{0},r)$ for some $r>r_{1}$, so that so that $y\not\in\textnormal{supp }u(T,\cdot)$, and in particular, $y\not\in C(\theta_{1},r_{2})\cap\textnormal{supp }u(T,\cdot)$ if $y\neq\tilde{x}$. This completes the proof in the large radius case for full data.

5. Stability

It is natural to take $f\in H_{D}(\Omega)\subset H_{0}^{1}(\Omega)$ when taking conservation of energy into consideration. By the above theorems, in both the small and large radius cases, we have that $Mf(t,\theta)$ is an elliptic FIO of order $-1/2$. We will model finite time measurements in $\left[0,T\right]$ by premultiplying $Mf(t,\theta)$ by $\chi\in C_{0}^{\infty}(\mathbb{R})$ with $\textnormal{supp }\chi\subset[-1,T_{1}]$ for $T<T_{1}<\infty$, and $\chi\equiv 1$ on $[0,T]$. Now because $M$ is an FIO of order $-1/2$ associated with the graph of the canonical relation $C$, we have $M^{*}$ is also an FIO of order $-1/2$, and it is associated with the canonical relation $C^{-1}$. $M^{*}\chi M$ then is an elliptic pseudodifferential operator of order $-1$. This implies that a parametrix $B$, which necessarily is an elliptic $\Psi$DO of order 1, exists such that

[TABLE]

where $R$ is a regularizing operator. We may assume that $B$ is a properly supported $\Psi$DO, which means $B,B^{t}:\mathcal{E}^{\prime}(\Omega)\rightarrow\mathcal{E}^{\prime}(\Omega)$. From this we see that

[TABLE]

$B$ is a continuous linear operator, so we have

[TABLE]

for some $C>0$, independent of $f$. And lastly, $M^{*}$ is a continuous linear operator, so we have from [4], Cor. 25.3.2.,

[TABLE]

Note that because we’ve multiplied $Mf(t,\theta)$ by $\chi$ and $Mf(t,\theta)$ has $\theta$ support in $S^{1}$, which is a compact manifold, that $\chi Mf(t,\theta)$ has compact support in $(0,T)\times S^{1}$, and so the norm above in $H^{3/2}([0,T]\times S^{1})$ is finite. By virtue of the injectivity of $\chi Mf(t,\theta)$, we may then write (at the cost of possibly increasing $C^{\prime}$)

[TABLE]

This gives stability of the measurement operator $M$.

Theorem 3**.**

Let $f\in H_{D}(\Omega)$, and $Mf(t,\theta)$ be defined as the either of the measurement operators above. If $\chi\in C_{0}^{\infty}(\mathbb{R})$ with $\chi\equiv 1$ on $[0,T]$, then we have the following stability estimate:

[TABLE]

5.1. Visible Singularities

A singularity $(y,\xi)\in\textnormal{WF}(f)$ is called visible on an open subset $U\times\Gamma$ of $\mathbb{R}\times S^{1}$ for $M$ if it creates a singularity in the measurement data $Mf\mid_{U\times\Gamma}$. Now because $Mf$ is an elliptic FIO associated with a local canonical diffeomorphism $C$ (See Theorems 1,2), we know by [4]

[TABLE]

Let $(t_{0},\boldsymbol{\theta}_{0})\in\mathbb{R}\times S^{1}$ and let $U\times\Gamma$ be an open neighborhood of $(t_{0},\boldsymbol{\theta}_{0})$. By the above arguments (Theorems 1 and 2), we know that singularities $(y,\xi)\in\textnormal{WF}(f)$ split and travel along geodesics $\left(\gamma_{y,\hat{\xi}}(t),c(y)|\xi|\dot{\gamma}_{y,\hat{\xi}}(t)\right)$, and that this will create a singularity at $(t_{0},\boldsymbol{\theta}_{0})$ if and only if the geodesic intersects the circular integrating detector with center $\boldsymbol{\theta}_{0}$ perpendicularly at time $t_{0}$, with no singularity to mask it intersecting the circular integrating detector at mirror points on the circle (i.e. an antipodal point in the small radius case, and a $C(\theta,r)$ mirror point in the large radius case). Therefore, to determine those singularities of $f$ that are visible from $U\times\Gamma$, we simply trace all geodesics that go through $R\boldsymbol{\theta}$ back to $\Omega$ and see if they have nonempty intersection with $\textnormal{WF}(f)$, see Figure 2.

For each $(t,\boldsymbol{\theta})\in U\times\Gamma$, $\lambda\in\mathbb{R}\setminus 0$, let

[TABLE]

and

[TABLE]

These are the sets of all points on geodesics intersecting the half circle $C^{+}(\theta,r)$ (respectively, $C^{-}(\theta,r)$) perpendicularly at time $t$, with tangent vector of magnitude $\lambda$. For

[TABLE]

define $\sim\left(\gamma_{(x,\frac{x-R\boldsymbol{\theta}}{r})}(t),\lambda\hat{\dot{\gamma}}_{(x,\frac{x-R\boldsymbol{\theta}}{r})}(t)\right)=\left(\gamma_{(\tilde{x},\frac{\tilde{x}-R\boldsymbol{\theta}}{r})}(t),\lambda\hat{\dot{\gamma}}_{(\tilde{x},\frac{\tilde{x}-R\boldsymbol{\theta}}{r})}(t)\right)\in A_{t,\theta,\lambda}^{\mp}$ where $\tilde{x}$ is the appropriate mirror point on the circular integrating detector, depending on the experimental setup. Now, $(y,\xi)\in\textnormal{WF}(f)$ is visible from $(t,\theta)\in U\times\Gamma$ if and only if $(y,\xi)\in A^{\pm}_{(t,\theta,\lambda)}$ for some $(t,\theta,\lambda)\in U\times\Gamma\times(\mathbb{R}\setminus 0)$ and $\sim(y,\xi)\not\in\textnormal{WF}(f)$. Let

[TABLE]

It then follows from the above arguments that the set of visible singularities is given by

[TABLE]

We have shown the following:

Theorem 4**.**

Let $U\times\Gamma\subset\mathbb{R}\times S^{1}$ be an open subset, and for each $(t,\theta,\lambda)\in U\times\Gamma\times(\mathbb{R}\setminus 0)$ let $A^{\pm}_{t,\theta,\lambda}$ and $B^{\pm}_{t,\theta,\lambda}$ be defined as above. Then in both the small radius detector case and the large radius detector case, the singularities of $f$ that are visible from $U\times\Gamma$ in the restricted data $Mf\mid_{U\times\Gamma}$ are given by

[TABLE]

From this we see that if $\Gamma=S^{1}$ and $(0,T]\subset U$, where $T=\sup_{x\in\Omega}\inf_{\theta\in[0,2\pi)}\text{dist}(x,C(\theta,r))$ (where the distance is the geodesic distance), then all singularities of $f$ are visible assuming $T<\infty$. We make a note that just because a singularity is not “visible” from some open set $U\times\Gamma$ does not mean that we cannot infer the existence of that singularity from the restricted data $Mf$.

6. Numerical Results

To simulate the collection of forward data, we numerically solve the wave equation with variable wave speed using the implementation of Perfectly Matched Layers (PML) found in [1] for a number of different smooth initial conditions. This ensures that measured data will only come from signals inside the region of interest, and not from reflections at the boundary of the window of computation. Then, we collect simulated measurement data on the unit circle $Mf(t,\theta)$ for $0\leq\theta<2\pi$ and $0\leq t<5$, for a specific initial condition. In general, the amount of time that we collect data should depend on the wave speed inside the medium we are imaging, and for the wave speed we have chosen of $1+0.3\sin(8x)\cos(5y)\eta(x,y)$ with $\eta(x,y)\in C_{0}^{\infty}(B_{1}(0))$, t=5s suffices as an appropriate time range. We’ve show the graph of the wave speed in Figure 3. We then use an iterative solver to reconstruct the smooth initial condition using the simulated data over the given time interval. The reconstruction shown in Figures 4 was made using the model $R=1$ and $r=2$ (the large radius detector model), with data taken on the full unit circle. An almost identical reconstruction is obtained if we use the small radius integrating detector model with full data.

We also run numerical simulations with data taken on an open subset of the unit circle. Here we take data for $\theta\in(-\pi/2,0)$, with the same wave speed interior to the object. We then multiply by a smooth cutoff function so as to not introduce new singularities into the reconstruction. Results for the partial data case are shown in Figure 5.