ABSTRACT

Direct and inverse scattering problems connected with the wave equation in non-homogeneous bounded domains constitute challenging actual subjects for both mathematicians and engineers. Among them one can mention, for example, inverse source problems in seismology, nondestructive archeological probing, mine prospecting, inverse initial-value problems in acoustic tomography, etc. In spite of its crucial importance, almost all of the available rigorous investigations concern the case of unbounded simple domains such as layered planar or cylindrical or spherical structures. The main reason for the lack of the works related to non-homogeneous bounded structures is the extreme complexity of the explicit expressions of the Green’s functions. The aim of the present work consists in discovering some universal properties of the Green’s functions in question, which reduce enormously the difficulties arising in various applications. The universality mentioned here means that the properties are not depend on the geometrical and physical properties of the configuration. To this end one considers first the case when the domain is partially-homogeneous. Then the results are generalized to the most general case. To show the importance of the universal properties in question, they are applied to an inverse initial-value problem connected with photo-acoustic tomography.

Keywords:

Green’s Functions, Inverse Source Problem, Inverse Initial-Value Problem, Tomography, Photo-Acoustic Tomography, Thermo-Acoustic Tomography, Wave Equation

1. Introduction

Let $V\subset {\Re }^3$ denote a bounded region composed of n sub-domains having different (constant) constitutive parameters, say e_j and m_j ( $j=1,\cdots ,n$ ) as shown in Figure 1 below. The constitutive parameters of the surrounding space, say $V_0$ , will be denoted by ${\epsilon }_0$ and ${\mu }_0$ . Let a function $p\left(x,t\right)$ , defined for $x\in {\Re }^3$ and $t\in {\Re }^{}$ , satisfy the wave equation

$\text{div}\left(\frac{1}{\epsilon }\text{grad}p\right)-\mu \frac{{\partial }^{2}p}{\partial t^2}=-p_0\left(x,t\right),x\in {\Re }^3,t\in \Re$ (1)

in the sense of distribution. Here $p_0\left(x,t\right)$ stands for the density of the source which will be assumed located inside V while $\epsilon =\epsilon \left(x\right)$ and $\mu =\mu \left(x\right)$ denote the above-mentioned constitutive parameters:

$\begin{array}{l}\epsilon =\epsilon \left(x\right)=\{\begin{array}{ll}{\epsilon }_j,\hfill & x\in V_j,j=1,\cdots ,n\hfill \\ {\epsilon }_0,\hfill & x\in V_0\hfill \end{array}\\ \mu =\mu \left(x\right)=\{\begin{array}{ll}{\mu }_{j,}\hfill & x\in V_j,j=1,\cdots ,n\hfill \\ {\mu }_0,\hfill & x\in V_0.\hfill \end{array}\end{array}$

The validity of (1) in the sense of distribution involves in itself both the boundary conditions satisfied on the interfaces between the adjacent sub-regions and the initial conditions, if any, at a certain time, say $t=0$ . Therefore it is not necessary to write down here these boundary and initial conditions explicitly (see for ex. (14b)). To determine $p\left(x,t\right)$ uniquely, one has to add to (1) the so- called radiation conditions also. They will be clarified later on (see (4a, b) below).

As is well known, in a direct propagation problem the so-called outgoing Green’s function plays important role in transporting the data known in the source region towards the observation points. This is a wave propagating towards infinity. But in an inverse source problem, the data observed (measured) on a surface S, surrounding the source region, is transported inversely towards the source region. This transportation is accomplished through a Green’s function which propagates in the inverse direction (ingoing Green’s function).

Therefore, in inverse source problems one has to use both of these Green’s functions together. Now our aim is to reveal some universal properties of these functions. Here the universality means that the properties in question do not depend on the geometrical shapes and constitutive parameters of the sub-domains.

In what follows we will assume that $p\left(x,t\right)$ is Fourier transformable with respect to $t$ . The transform of each quantity are denoted by a hat on that quantity. For example one has

$\stackrel{⌢}{p}\left(x,\omega \right)=\frac{\text{1}}{\text{2π}}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}p\left(x,t\right){\text{e}}^{i\omega t}\text{d}t},\omega \in \Re$ (2a)

and, inversely

$\frac{1}{2}\left[p\left(x,t+0\right)+p\left(x,t-0\right)\right]={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\stackrel{⌢}{p}\left(x,\omega \right){\text{e}}^{-i\omega t}\text{d}\omega },t\in \Re ,$ (2b)

Thus the transform of (1) is as follows:

$\text{div}\left(\frac{1}{\epsilon }\text{grad}\stackrel{⌢}{p}\right)+{\omega }^2\mu \stackrel{⌢}{p}=-{\stackrel{⌢}{p}}_0\left(x,\omega \right)$ (3)

The radiation conditions satisfied by $\stackrel{⌢}{p}\left(x,\omega \right)$ are

$\stackrel{⌢}{p}\left(x,\omega \right)=O\left(1/\left|x\right|\right)\text{as}\left|x\right|\to \infty$ (4a)

and

$\partial \stackrel{⌢}{p}/\partial \left|x\right|-i\omega /c_0\stackrel{⌢}{p}=O\left(1/{\left|x\right|}^2\right)\text{as}\left|x\right|\to \infty$ (4b)

where we put

$c_0=1/\sqrt{{\epsilon }_0{\mu }_0}$ (4c)

The outgoing Green’s function associated with (3)-(4c) is the function $G_1\left(x,\eta ,w\right)$ determined uniquely through the following relations:

$di{v}_x\left(\frac{1}{\epsilon \left(x\right)}gra{d}_x{G}_1\right)+{\omega }^2\mu \left(x\right)G_1=-\delta \left(x-\eta \right)$ (5a)

$G_1\left(x,\eta ,\omega \right)=O\left(1/\left|x\right|\right)\text{as}\left|x\right|\to \infty$ (5b)

$\partial G_1/\partial \left|x\right|-i\omega /c_0{G}_1=O\left(1/{\left|x\right|}^2\right)\text{as}\left|x\right|\to \infty$ (5c)

Here $\eta \in {\Re }^3$ stands for any point while $\delta \left(x-\eta \right)$ is the Dirac distribution concentrated at $x=\eta$ . The sub-indices $x$ in ${\text{div}}_x$ and ${\text{grad}}_x$ , which appear in (5a), means that the derivatives are to be computed with respect to the components of $x$ . Notice that because of the radiation condition (5c), the wave associated with G₁ propagates in the direction $\eta \to x$ . This Green’s function permits one to transport the data known at the source point $\eta$ to the observation point $x$ . Indeed, as we will see in Section 2.1, the solution to the problem (3)-(4c) (i.e. the direct problem) is as follows:

$\stackrel{⌢}{p}\left(x,\omega \right)={\displaystyle \underset{{\Re }^3}{\int }{G}_1\left(\eta ,x,\omega \right){\stackrel{⌢}{p}}_0\left(\eta ,\omega \right)\text{d}{v}_{\eta }}$ (6)

As to the ingoing Green’s function $G_2\left(x,\eta ,\omega \right)$ , it satisfies

${\text{div}}_x\left(\frac{1}{\epsilon \left(x\right)}{\text{grad}}_x{G}_2\right)+{\omega }^2\mu \left(x\right)G_2=-\delta \left(x-\eta \right)$ (7a)

in the sense of distribution under the radiation conditions

$G_2\left(x,\eta ,\omega \right)=O\left(1/\left|x\right|\right)\text{as}\left|x\right|\to \infty$ (7b)

and

$\partial G_2/\partial \left|x\right|+i\omega /c_0{G}_2=O\left(1/{\left|x\right|}^2\right)\text{as}\left|x\right|\to \infty .$ (7c)

Notice that (7a) and (7b) are quite identical to (5a) and (5b) while (7c) differs from (5c) only with the sign of (i). Because of this difference, the wave associated with $G_2\left(x,\eta ,\omega \right)$ propagates in the direction $x\to \eta$ and serves to transport the measured data at $x$ to the source point $\eta$ .

Observe that the Green’s functions defined above differ from the classical ones by the factor $\epsilon \left(x\right)$ (see [1] ). In Section 2 we will prove that these Green’s functions have the universal properties stated in the following theorems:

Theorem-1. For all $x,\eta \in {\Re }^3$ the Green’s function $G_1\left(x,\eta ,\omega \right)$ is symmetric with respect to $x$ and $\eta$ :

$G_1\left(x,\eta ,\omega \right)\equiv G_1\left(\eta ,x,\omega \right).$ (8)

Theorem-2. For $\omega \in \Re$ the Green’s functions $G_1\left(x,\eta ,\omega \right)$ and $G_2\left(x,\eta ,\omega \right)$ are interrelated with the following functional relations:

$G_2\left(x,\eta ,\omega \right)\equiv G_1\left(x,\eta ,-\omega \right),$ (9a)

$G_2\left(x,\eta ,\omega \right)\equiv G_1^{*}\left(x,\eta ,\omega \right),$ (9b)

$G_2\left(x,\eta ,\omega \right)\equiv G_2\left(\eta ,x,\omega \right).$ (9c)

Here $(*)$ means the complex conjugate. The relations in (9a) and (9c) are also valid for complex $\omega$ .

Theorem-3. For all $x,\eta \in {\Re }^3$ and $n\to \infty$ one has

${\displaystyle \underset{-{\Omega }_n}{\overset{{\Omega }_n}{\int }}{G}_1\left(x,\eta ,\omega \right)\omega \text{d}\omega }\to \frac{i\text{π}}{\mu \left(\eta \right)}\delta \left(x-\eta \right)$ (10a)

$\to \frac{i}{2\mu \left(\eta \right)}\frac{\delta \left(\left|x-\eta \right|\right)}{{\left|x-\eta \right|}^2}$ (10b)

$\to -\frac{i}{2\mu \left(\eta \right)}\frac{{\delta }^{\prime }\left(\left|x-\eta \right|\right)}{\left|x-\eta \right|}$ (10c)

Here $\left\{{\Omega }_n\right\}$ stands for any sequence such that ${\Omega }_n\to \infty$ as $n\to \infty$ .

Theorem-4. For all $x,\eta \in {\Re }^3$ and $n\to \infty$ one has

${\displaystyle \underset{-{\Omega }_n}{\overset{{\Omega }_n}{\int }}\left[G_2\left(x,\eta ,\omega \right)-G_1\left(x,\eta ,\omega \right)\right]\omega \text{d}\omega }\to -\frac{2\text{π}i}{\mu \left(\eta \right)}\delta \left(x-\eta \right)$ (11a)

$\to -\frac{i}{\mu \left(\eta \right)}\frac{\delta \left(\left|x-\eta \right|\right)}{{\left|x-\eta \right|}^2}$ (11b)

$\to \frac{i}{\mu \left(\eta \right)}\frac{{\delta }^{\prime }\left(\left|x-\eta \right|\right)}{\left|x-\eta \right|}.$ (11c)

Here $\left\{{\Omega }_n\right\}$ stands for any sequence such that ${\Omega }_n\to \infty$ as $n\to \infty$ .

Theorem-5. For all $x\in V$ , $\eta \in V_0$ , $\alpha \in \left[-\text{π}/\text{2},\text{π}/\text{2}\right]$ and $n\to \infty$ one has

$G_1\left(x,\eta ,{\Omega }_n\sin \alpha \right)\to \frac{{\epsilon }_0}{4\text{π}\left|x-\eta \right|}{\text{e}}^{i{\Omega }_n\sin \alpha \left|x-\eta \right|/c_0}\left[1+O\left(1/{\Omega }_n\right)\right]$ (12a)

$G_2\left(x,\eta ,{\Omega }_n\sin \alpha \right)\to \frac{{\epsilon }_0}{4\text{π}\left|x-\eta \right|}{\text{e}}^{-i{\Omega }_n\sin \alpha \left|x-\eta \right|/c_0}\left[1+O\left(1/{\Omega }_n\right)\right]$ (12b)

Here $\left\{{\Omega }_n\right\}$ stands for any sequence such that ${\Omega }_n\to \infty$ as $n\to \infty$ .

Theorem-6. For all $x,y\in V$ , $\eta \in V_0$ and $n\to \infty$ one has

${\displaystyle \underset{-{\Omega }_n}{\overset{{\Omega }_n}{\int }}{G}_1\left(x,\eta ,\omega \right)G_2\left(y,\eta ,\omega \right)\omega \text{d}\omega }\to -\frac{i}{8\text{π}}\frac{{\epsilon }_0}{{\mu }_0}\frac{{\delta }^{\prime }\left(R_1-R_2\right)}{R_1{R}_2}.$ (13)

Here $\left\{{\Omega }_n\right\}$ stands for any sequence such that ${\Omega }_n\to \infty$ as $n\to \infty$ and

$R_1=\left|x-\eta \right|,R_2=\left|y-\eta \right|.$

2. Proofs of the Theorems

Now we will give the proofs of the above-mentioned theorems in turn. But, first of all we prove the Formula (6) which will be used later on in the proofs of the theorems.

2.1. Proof of the Formula (6)

Multiply first (1) with $G_1$ and (5a) with $\stackrel{⌢}{p}$ , and subtract side by side to get

$\begin{array}{l}{\text{div}}_x\left[\frac{1}{\epsilon \left(x\right)}\left\{G_1{\text{grad}}_x\stackrel{⌢}{p}-\stackrel{⌢}{p}{\text{grad}}_x{G}_1\right\}\right]\\ =-{\stackrel{⌢}{p}}_0\left(x,\omega \right)G_1\left(x,\eta ,\omega \right)+\stackrel{⌢}{p}\left(x,\omega \right)\delta \left(x-\eta \right)\text{.}\end{array}$

If we integrate both sides in a sphere of radius R and transform the first side to the surface integral on $\left|x\right|=R$ , then we write

$\begin{array}{l}{\displaystyle \underset{\left|x\right|=R}{\int }\left[\frac{1}{\epsilon \left(x\right)}\left\{G_1{\text{grad}}_x\stackrel{⌢}{p}-\stackrel{⌢}{p}{\text{grad}}_x{G}_1\right\}\right]\cdot \text{d}{S}_x}\\ =\stackrel{⌢}{p}\left(\eta ,\omega \right)-{\displaystyle \underset{\left|x\right|\le R}{\int }{G}_1\left(x,\eta ,\omega \right){\stackrel{⌢}{p}}_0\left(x,\omega \right)\text{d}{v}_x.}\end{array}$

Now let us make $R\to \infty$ . Because of the radiation conditions satisfied by $G_1$ and $\stackrel{⌢}{p}$ , the surface integral tends to zero and yields (6).

2.2. Proof of Theorem-1

Now consider the function ${\stackrel{⌢}{p}}_1\left(x,\omega \right)$ defined as follows:

${\stackrel{⌢}{p}}_1\left(x,\omega \right)={\displaystyle \underset{{\Re }^3}{\int }{G}_1\left(x,\eta ,\omega \right){\stackrel{⌢}{p}}_0\left(\eta ,\omega \right)\text{d}{V}_{\eta }.}$

One can easily check that ${\stackrel{⌢}{p}}_1$ satisfies (3) as well as the radiation conditions (4a, b). Therefore it is equal to $\stackrel{⌢}{p}\left(x,\omega \right).$ If we compare this new expression of $\stackrel{⌢}{p}\left(x,\omega \right)$ with (6), then we write also

${\displaystyle \underset{{\Re }^3}{\int }\left[G_1\left(x,\eta ,\omega \right)-G_1\left(\eta ,x,\omega \right)\right]{\stackrel{⌢}{p}}_0\left(\eta ,\omega \right)\text{d}{V}_{\eta }\equiv 0.}$

Since ${\stackrel{⌢}{p}}_0\left(\eta ,\omega \right)$ is arbitrary, from the latter one gets (8).

2.3. Proof of Theorem-2

Assume $\omega \in \Re$ and make the change $\omega \to \left(-\omega \right)$ in the differential equation and the radiation conditions satisfied by $G_1$ . One obviously obtains the equations satisfied by $G_2$ . From this one concludes (9a).

If under the same assumption $\omega \in \Re$ one considers the complex conjugates of the equations satisfied by $G_1$ , one sees that the resulting equations are nothing but those satisfied by $G_2$ . This proves (9b).

As to (9c), it is a direct consequence of (9a) and Theorem-1. The relations in (9a) and (9c) are obviously valid for complex $\omega$ also.

2.4. Proof of Theorem-3

Consider the Equation (1) for the case when $p_0\left(x,t\right)=p_0\left(x\right)T_0\left(t\right)$ with $T_0\left(t\right)\equiv 0$ for $t<0$ . In this case one obviously has $p\left(x,t\right)\equiv 0$ for $t<0$ . Thus from (2a) and (2b) written for $t=0$ one gets

$\frac{1}{2}\frac{\partial p}{\partial t}\left(x,+0\right)=-i\underset{n\to \infty }{\lim }{\displaystyle \underset{V}{\int }{\displaystyle \underset{{\Omega }_n}{\overset{{\Omega }_n}{\int }}\omega G_1\left(x,\eta ,\omega \right)p_0\left(\eta \right){\hat{T}}_0\left(\omega \right)\text{d}\omega \text{d}{v}_{\eta }}}$ (14a)

Assume first that $T_0\left(t\right)=\delta \left(t\right)$ which yields ${\hat{T}}_0(\omega )=1/\left(2\text{π}\right)$ . In this case from (1) one gets

$\frac{\partial p}{\partial t}\left(x,+0\right)=\frac{p_0\left(x\right)}{\mu \left(x\right)}$ (14b)

A comparison of (14a) with (14b) yields

$\underset{n\to \infty }{\lim }{\displaystyle \underset{V}{\int }{p}_0\left(\eta \right)Q_n\left(\eta ,x\right)\mathrm{d}{v}_{\eta }=\frac{i\text{π}}{\mu \left(x\right)}{p}_0\left(x\right)}$ (14c)

with

$Q_n\left(\eta ,x\right)={\displaystyle \underset{-{\Omega }_n}{\overset{{\Omega }_n}{\int }}{G}_1\left(\eta ,x,\omega \right)\omega \text{d}\omega ,}n=1,2,\cdots .$

From (14c) one concludes that when $n\to \infty ,$ the sequence $Q_n$ defines a generalized function (= distribution) [2] , namely:

$n\to \infty \Rightarrow Q_n\left(\eta ,x\right)\to \frac{i\text{π}}{\mu \left(\eta \right)}\delta \left(\eta -x\right).$

The latter proves (10a).

To prove (10b), let us define a polar coordinate system $\left(R,\theta ,\phi \right)$ with origin at the point $\eta$ and consider the distribution $\delta \left(R\right)/R^2$ . If $\varphi \left(x\right)\equiv \varphi \left(x_1,x_2,x_3\right)$ denotes any test function, then we write

$\langle \frac{\delta \left(R\right)}{R^2},\varphi \left(x\right)\rangle ={\displaystyle \underset{0}{\overset{2\text{π}}{\int }}{\displaystyle \underset{0}{\overset{\text{π}}{\int }}{\displaystyle \underset{0}{\overset{\infty }{\int }}\delta \left(R\right)\psi }}}\left(R,\theta ,\phi \right)\sin \theta \text{d}R\text{d}\theta d\phi ,$ (15a)

where we put

$\psi \left(R,\theta ,\phi \right)=\varphi \left({\eta }_1+R\sin \theta \cos \varphi ,\cdots \right)$

Now it is worthwhile to remark that the transformation to the polar coordinates causes to a confusion in (15a) because the distribution $\delta \left(R\right)$ is defined on the space of the test functions defined for $R\in \left(-\infty ,\infty \right)$ while the function $\psi \left(R,\theta ,\phi \right)$ is known only for $R\in \left[0,\infty \right)$ . To overcome the difficulty, we extend the functions $\psi \left(R,\theta ,\phi \right)$ defined in $0\le R<\infty$ into $\left(-\infty \right)<R<\infty$ as an even function (which is continuous at $R=0$ !!). If the extended function is denoted by $\tilde{\psi }\left(R\right)$ , then we write

$\begin{array}{c}\psi \left(0\right)=\tilde{\psi }\left(0\right)=\langle \delta \left(R\right),\tilde{\psi }\left(R\right)\rangle \\ ={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\delta \left(R\right)\tilde{\psi }\left(R\right)\text{d}R}\\ =\left({\displaystyle \underset{-\infty }{\overset{0}{\int }}+{\displaystyle \underset{0}{\overset{\infty }{\int }}}}\right)\delta \left(R\right)\tilde{\psi }\left(R\right)\text{d}R\\ =2{\displaystyle \underset{0}{\overset{\infty }{\int }}\delta \left(R\right)\psi \left(R\right)\text{d}R}\end{array}$

which yields

${\displaystyle \underset{0}{\overset{\infty }{\int }}\delta \left(R\right)\psi \left(R\right)\text{d}R}=\frac{1}{2}\psi \left(\text{0}\right)$ (15b)

This formula defines the extension of $\delta \left(x\right)$ into $\delta \left(R\right)$ to be used in $R^3$ and $R^2$ .

From (15a) and (15b) one gets

$\langle \frac{\delta \left(R\right)}{R^2},\varphi \left(x\right)\rangle =2\text{π}\varphi \left(\eta \right)=\langle 2\text{π}\delta \left(x-\eta \right),\varphi \left(x\right)\rangle$

and

$\delta \left(x-\eta \right)=\frac{1}{2\text{π}}\frac{\delta \left(\left|x-\eta \right|\right)}{{\left|x-\eta \right|}^2}$ (15c)

which reduces (10a) to (10b).

Finally, let us consider the distribution ${\delta }^{\prime }\left(R\right)/R$ and write

$\begin{array}{c}\langle \frac{{\delta }^{\prime }\left(R\right)}{R},\varphi \left(x\right)\rangle ={\displaystyle \underset{0}{\overset{2\text{π}}{\int }}{\displaystyle \underset{0}{\overset{\text{π}}{\int }}{\displaystyle \underset{0}{\overset{\infty }{\int }}{\delta }^{\prime }\left(R\right)\chi }}}\left(R,\theta ,\phi \right)\sin \theta \text{d}R\text{d}\theta \text{d}\phi \\ =-{\displaystyle \underset{0}{\overset{2\text{π}}{\int }}{\displaystyle \underset{0}{\overset{\text{π}}{\int }}{\displaystyle \underset{0}{\overset{\infty }{\int }}\delta \left(R\right)\frac{\partial }{\partial R}\chi \left(R,\theta ,\phi \right)\sin \theta \text{d}R\text{d}\theta \text{d}\phi }}}\end{array}$

with

$\chi \left(R,\theta ,\phi \right)=R\varphi \left({\eta }_1+R\sin \theta \cos \phi ,\cdots \right)$ .

By virtue of (15b) we write also

$\begin{array}{c}\langle \frac{{\delta }^{\prime }\left(R\right)}{R},\varphi \left(x\right)\rangle =-2\text{π}\varphi \left(\eta \right)\\ =\langle -2\text{π}\delta \left(x-\eta \right),\varphi \left(x\right)\rangle .\end{array}$

This shows that

$\delta \left(x-\eta \right)=-\frac{1}{2\text{π}}\frac{{\delta }^{\prime }\left(\left|x-\eta \right|\right)}{\left|x-\eta \right|}$ (15d)

which reduces (10a) to (10c).

Remark. All the results obtained above can easily be checked in the simplest case of homogeneous space for which one has

$G_1\left(x,\eta ,\omega \right)=\frac{{\epsilon }_0}{4\text{π}R}{\text{e}}^{i\omega R/c_0},$ (16a)

and

$G_2\left(x,\eta ,\omega \right)=\frac{{\epsilon }_0}{4\text{π}R}{\text{e}}^{-i\omega R/c_0},$ (16b)

where $R=\left|x-\eta \right|$ . To this end one has to observe that in ³

$\frac{{\epsilon }_0}{2iR}{\delta }^{\prime }\left(\frac{R}{c_0}\right)=i\frac{\text{π}}{{\mu }_0}\delta \left(x-\eta \right).$

2.5. Proof of Theorem-4

By considering (9a) we get directly

$\begin{array}{c}{\displaystyle \underset{-{\Omega }_n}{\overset{{\Omega }_n}{\int }}{G}_2\left(x,\eta ,\omega \right)\omega \text{d}\omega }={\displaystyle \underset{-{\Omega }_n}{\overset{{\Omega }_n}{\int }}{G}_1\left(x,\eta ,-\omega \right)\omega \mathrm{d}\omega }\\ =-{\displaystyle \underset{-{\Omega }_n}{\overset{{\Omega }_n}{\int }}{G}_1\left(x,\eta ,\omega \right)\omega \text{d}\omega }\end{array}$

Thus from (10a), (10b) and (10c) one gets directly (11a), (11b) and (11c).

2.6. Proof of Theorem-5

Since the function $G_1\left(x,\eta ,\omega \right)$ is symmetrical with respect to $x$ and $\eta$ , as a function of $\eta \in V_0$ , it satisfies the Helmholtz equation

${\Delta }_{\eta }{G}_1+k_0^2{G}_1=0,\eta \in V_0,x\in V$ (17a)

under the radiation conditions

$G_1=O\left(1/\left|\eta \right|\right),\frac{\partial G_1}{\partial \left|\eta \right|}-i{k}_0{G}_1=O\left(1/{\left|\eta \right|}^2\right)$ (17b)

Here we put $k_0=\omega \sqrt{{\epsilon }_0{\mu }_0}$ . Now let us define a polar co-ordinate system $\left(R,\theta ,\phi \right)$ whose origin is at the point $x$ and $R=\left|\eta -x\right|$ . In this system (17a) and (17b) become

$\frac{{\partial }^2{G}_1}{\partial R^2}+\frac{2}{R}\frac{\partial G_1}{\partial R}+\frac{1}{R^2\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial G_1}{\partial \theta }\right)\text{+}\frac{1}{R^2{\sin }^2\theta }\frac{{\partial }^2{G}_1}{\partial {\phi }^2}+\frac{{\omega }^2}{c_0^2}{G}_1=0$ (18a)

and

$R\to \infty \Rightarrow G_1=O\left(1/R\right),$ (18b)

$R\to \infty \Rightarrow \frac{\partial G_1}{\partial R}-i\frac{\omega }{c_0}{G}_1=O\left(1/R^2\right)$ (18c)

Now for each ${\Omega }_n$ let us make the substitutions

$\begin{array}{l}\omega ={\Omega }_n\sin \alpha ,\alpha \in \left[-\text{π}/2,\text{π}/2\right]\\ \Rightarrow \frac{1}{{\Omega }_n}{k}_0=\frac{\sin \alpha }{c_0}\equiv {\tilde{k}}_0,\end{array}$

and

$\begin{array}{c}{\Omega }_{n}R=r\\ \Rightarrow \frac{\partial G_1}{\partial R}={\Omega }_n\frac{\partial G_1}{\partial r},\frac{{\partial }^2{G}_1}{\partial R^2}={\Omega }_n^2\frac{{\partial }^2{G}_1}{\partial r^2}\\ \Rightarrow \frac{1}{{\Omega }_n}{G}_1\left(x,\eta ,\omega \right)=G_{1n}\left(r,\theta ,\phi ,\alpha \right),\end{array}$

which reduce (18 a-c) to

$\frac{{\partial }^2{G}_{1n}}{\partial r^2}+\frac{2}{r}\frac{\partial G_{1n}}{\partial r}+\frac{1}{r^2\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial G_{1n}}{\partial \theta }\right)+\frac{1}{r^2{\sin }^2\theta }\frac{{\partial }^2{G}_{1n}}{\partial {\phi }^2}+{\tilde{k}}_0^2{G}_{1n}=0$ (19a)

and

$r\to \infty \Rightarrow G_{1n}=O\left(1/r\right),$ (19b)

$\frac{\partial G_{1n}}{\partial r}-i{\tilde{k}}_0{G}_{1n}=O\left(1/r^2\right).$ (19c)

Notice that (19a) is satisfied outside the region V, and $G_{1n}$ may depend on ${\Omega }_n$ (and $\alpha$ ) through the boundary conditions on $\partial V$ also.

It is obvious that (19a) and (19b, c) define the radiating (=outgoing) solution to the Helmholtz equation. When $n\to \infty$ , one has ${\Omega }_n\to \infty \Rightarrow r={\Omega }_{n}R\to \infty$ , which reduces the expression of $G_1{}_n$ to the so-called far field pattern of ${\tilde{G}}_{1n}$ [3] :

$n\to \infty \Rightarrow G_{1n}=\frac{1}{r}{\text{e}}^{i{\tilde{k}}_{0}r}\left\{A_{\infty }+O\left(\frac{1}{r}\right)\right\}.$ (20)

Here $A_{\infty }$ stands for the scattering amplitude. It does not depend on r. By using (20), we write successively

$\begin{array}{l}{\displaystyle \underset{-{\Omega }_n}{\overset{{\Omega }_n}{\int }}{G}_1\left(x,\eta ,\omega \right)\omega \text{d}\omega }={\displaystyle \underset{-\text{π}/2}{\overset{\text{π}/2}{\int }}\left[{\Omega }_n{G}_1{}_n\left(r,\theta ,\phi ,\alpha \right)\right]{\Omega }_n^2\sin \alpha \cos \alpha \text{d}\alpha }\\ \to {\displaystyle \underset{-\text{π}/2}{\overset{\text{π}/2}{\int }}\left[{\Omega }_n{A}_{\infty }\frac{1}{r}{\text{e}}^{i{\tilde{k}}_{0}r)}\right]{\Omega }_n^2\sin \alpha \cos \alpha \text{d}\alpha }\\ =\frac{1}{R}{\displaystyle \underset{-{\Omega }_n}{\overset{{\Omega }_n}{\int }}{A}_{\infty }{\text{e}}^{i\omega R/c_0}\omega \text{d}\omega }\\ =\frac{c_0^2}{R}{\displaystyle \underset{-{\Omega }_n/c_0}{\overset{{\Omega }_n/c_0}{\int }}{A}_{\infty }{\text{e}}^{i\omega R}\omega \text{d}\omega .}(21)\end{array}$

If we compare (21) with (10c) written as follows

$\begin{array}{c}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{G}_1\left(x,\eta ,\omega \right)\omega \text{d}\omega }=-\frac{i}{2{\mu }_0}\frac{{\delta }^{\prime }\left(R\right)}{R}\\ =\frac{1}{4\text{π}{\mu }_0}\frac{1}{R}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\text{e}}^{i\omega R}\omega \text{d}\omega },\end{array}$

then we conclude that

$A_{\infty }\equiv \frac{{\epsilon }_0}{4\text{π}}$

and the explicit expression of (20) is

$G_{\text{1}}{}_n\left(x,\eta ,{\Omega }_n\sin \alpha \right)\to \frac{{\epsilon }_0}{4\text{π}}\frac{1}{r}{\text{e}}^{ir\sin \alpha /c_0}$

as $n\to \infty$ . This yields

$n\to \infty \Rightarrow G_1\left(x,\eta ,{\Omega }_n\sin \alpha \right)\to \frac{{\epsilon }_0}{4\text{π}}\frac{1}{R}{\text{e}}^{i{\Omega }_n\sin \alpha R/c_0}.$ (22)

The latter proves (12a). The proof of (12b) is a direct result of (12a) and (9a).

It is important to remark here that (22) is not the far-field expression of $G_{\text{1}}$ because $R$ is finite. It states merely that when computing the integral appearing on the left-hand side of (21), for $n\to \infty$ $G_1$ can be replaced by (22).

2.7. Proof of Theorem-6

Let us now define two polar coordinates $\left(R_1,{\theta }_1,{\phi }_1\right)$ and $\left(R_2,{\theta }_2,{\phi }_2\right)$ whose origins are located at the points $x$ and $y$ , respectively. If we make the substitutions

$R_1=\left|x-\eta \right|,R_2=\left|y-\eta \right|,$

and

$\begin{array}{l}{r}_1={\Omega }_n{R}_1,r_2={\Omega }_n{R}_2,\omega ={\Omega }_n\sin \alpha \\ \frac{1}{{\Omega }_n}{G}_1\left(x,\eta ,\omega \right)=G_{1n}\left(r_1,{\theta }_1,{\phi }_1,\alpha \right),\end{array}$ (23a)

$\frac{1}{{\Omega }_n}{G}_2\left(y,\eta ,\omega \right)=G_{2n}\left(r_2,{\theta }_2,{\phi }_2,\alpha \right),$ (23b)

then from (22) we write

$n\to \infty \Rightarrow {\Omega }_n{\text{G}}_{\text{1}}{}_n\to \frac{{\epsilon }_0}{4\text{π}}\frac{1}{R_1}{\text{e}}^{i\omega R_1/c_0}$ (24a)

and

$n\to \infty \Rightarrow {\Omega }_n{\text{G}}_{\text{2}}{}_n\to \frac{{\epsilon }_0}{4\text{π}}\frac{1}{R_2}{\text{e}}^{-i\omega R_2/c_0}.$ (24b)

Now let us define the distribution defined through the sequence [2]

$W_n={\displaystyle \underset{-{\Omega }_n}{\overset{{\Omega }_n}{\int }}{G}_1\left(x,\eta ,\omega \right)G_2\left(y,\mu ,\omega \right)\omega \text{d}\omega }.$

By virtue of (24a, b) we have

$\begin{array}{l}{W}_n={\displaystyle \underset{-\text{π}/2}{\overset{\text{π}/2}{\int }}\left[{\Omega }_n{\text{G}}_{\text{1}}{}_{\text{n}}\left({\text{r}}_{\text{1}}\text{,}{\theta }_{\text{1}}\text{,}{\phi }_{\text{1}}\text{,}\alpha \right)\right]}\left[{\Omega }_n{G}_{2n}\left(r_2,{\theta }_2,{\phi }_2,\alpha \right)\right]{\Omega }_n^2\sin \alpha \cos \alpha \text{d}\alpha \\ \to {\left(\frac{{\epsilon }_0}{4\text{π}}\right)}^2\frac{1}{R_1{R}_2}{\displaystyle \underset{-{\Omega }_n}{\overset{{\Omega }_n}{\int }}{\text{e}}^{i\omega \left(R_1-R_2\right)/c_0}\omega \text{d}\omega }\\ \to {\left(\frac{{\epsilon }_0}{4\text{π}}\right)}^2\frac{1}{R_1{R}_2}\left(-2\text{π}i\right){\delta }^{\prime }\left(\frac{R_1-R_2}{c_0}\right)\\ \to -\frac{i}{8\text{π}}\frac{{\epsilon }_0}{{\mu }_0}\frac{{\delta }^{\prime }\left(R_1-R_2\right)}{R_2{R}_2}.(24\text{c})\end{array}$

This proves (13).

Remark. (24c) can easily be checked in the simplest case of homogeneous space for which $G_1$ and $G_2$ are given by (16a, b).

3. Application. Inverse Initial-Value Problem of Acoustic Tomography in a Non-Homogeneous Domain

Reconsider the domain shown in Figure 1 above. Such a domain can model, for example, small-seized animals, women’s breasts or other biological domains where different regions are the bone, skin, muscle, fat, cancerous tissues, epileptic tumor etc. If $V$ is exposed to a light (or thermal) pulse, excited at a certain time, say $t=0$ , then the energy stored in various regions creates a pressure wave $p\left(x,t\right)$ which propagates towards $V_0$ . The initial value of the pressure function, say $p\left(x,+0\right)$ , gives an idea about the configuration. The so-called photo-acoustic and thermo-acoustic tomographies are based on this phenomenon. To this end one first measures the pressure intensity $p\left(x,t\right)$ on a closed surface S lying in $V_0$ in a certain time interval (0, T) and then inserts them into an integral on S to obtain $p\left(x,+0\right)$ in V (see Figure 2). That means that the acoustic tomography consists in an inverse initial-value problem. The pressure function $p\left(x,t\right)$ satisfies the Equation (1) with $p_0\left(x,t\right)=p_0\left(x\right){\delta }^{\prime }\left(t\right)$ , namely [4] :

$\text{div}\left(\frac{1}{\epsilon }\text{grad}p\right)-\mu \frac{{\partial }^{2}p}{\partial t^2}=-\frac{1}{\epsilon }{p}_0\left(x\right){\delta }^{\prime }\left(t\right),$ (25)

where $x\in {\Re }^3$ and $t\in \Re$ . It is worthwhile to remark that in realistic problems both the source density $p_0\left(x\right)$ and the constitutive parameters $\left\{\epsilon \left(x\right),\mu \left(x\right)\right\}$ (i.e. the geometrical shapes of the regions $V_j$ ) are not known beforehand. Remark also that all of the previous theoretical works concern the simplest case where the space is assumed to be homogeneous and the aim consists in the determination of the source density $p_0\left(x\right)$ (inverse source problem) [4] - [12] . The basic principles and different implementations along with a long reference list are given in [4] and [5] .

Let $G_1\left(x,\eta ,\omega \right)$ and $G_2\left(x,\eta ,\omega \right)$ be the Green’s functions associated with the configuration shown in Figure 1. The pressure $p\left(\eta ,t\right)$ measured at a point $\eta \in S\in V_0$ is connected to the source density $p_0\left(x\right)$ through the outgoing Green’s function $G_1\left(x,\eta ,\omega \right)$ as given in Formula (6), where

${\hat{p}}_0\text{(}x,\omega \text{)}=-i\omega p_0\left(x\right)/\left[2\text{π}\epsilon \left(x\right)\right]$ .

Now, by using the ingoing Green’s function $G_2\left(x,\eta ,\omega \right)$ we transport the data known at the point $\eta \in S$ to the point $x\in V$ through the following function $q_n\left(x\right)$

$q_n\left(x\right)={\displaystyle \underset{-{\Omega }_n}{\overset{{\Omega }_n}{\int }}{\displaystyle \underset{S}{\int }\frac{1}{\epsilon \left(\eta \right)}\stackrel{⌢}{p}\left(\eta ,\omega \right)\frac{\partial G_2}{\partial n_{\eta }}}}\left(x,\eta ,\omega \right)\text{d}{S}_{\eta }\text{d}\omega .$ (26)

Here $\left\{{\Omega }_n\right\}$ stands for any sequence such that ${\Omega }_n\to \infty$ as $n\to \infty$ , $n_{\eta }$ is the unit outward normal vector to the surface S at the point $\eta$ and $\partial /\partial n_{\eta }$ means the derivative in this direction (see Figure 2). First we will show that when $n\to \infty$ one has

$-2q_n\left(x\right)\to p\left(x,+0\right),x\in V.$ (27)

To prove (27), let us replace $\hat{p}\left(\eta \text{,}\omega \right)$ in (26) by (6) and write

$q_n\left(x\right)=\frac{1}{2\text{π}}{\displaystyle \underset{V}{\int }{p}_0\left(y\right)P_n\left(y,x\right)\text{d}{v}_y}$

where we put

$P_n\left(y,x\right)=-\frac{i}{\epsilon \left(y\right)}{\displaystyle \underset{-{\Omega }_n}{\overset{{\Omega }_n}{\int }}\left\{{\displaystyle \underset{S}{\int }\frac{1}{\epsilon \left(\eta \right)}{G}_1\left(y,\eta ,\omega \right)\frac{\partial G_2}{\partial n_{\eta }}\left(x,\eta ,\omega \right)\text{d}{S}_{\eta }}\right\}\omega \text{d}\omega }$ (28)

(28) can also be written as follows (in what follows $D_S$ shows the domain bounded by S):

$\begin{array}{c}{P}_n\left(y,x\right)\\ =-\frac{i}{\epsilon \left(y\right)}{\displaystyle \underset{-{\Omega }_n}{\overset{{\Omega }_n}{\int }}\left\{{\displaystyle \underset{S}{\int }{G}_1\left(y,\eta ,\omega \right)\frac{1}{\epsilon \left(\eta \right)}{\text{grad}}_{\eta }{G}_2\left(x,\eta ,\omega \right)\cdot \text{d}{S}_{\eta }}\right\}\omega \text{d}\omega }\\ =-\frac{i}{\epsilon \left(y\right)}{\displaystyle \underset{-{\Omega }_n}{\overset{{\Omega }_n}{\int }}\left\{{\displaystyle \underset{D_S}{\int }{\text{div}}_{\eta }\left[G_1\left(y,\eta ,\omega \right)\frac{1}{\epsilon \left(\eta \right)}{\text{grad}}_{\eta }{G}_2\left(x,\eta ,\omega \right)\right]\text{d}{v}_{\eta }}\right\}\omega \text{d}\omega }\\ =-\frac{i}{\epsilon \left(y\right)}{\displaystyle \underset{-{\Omega }_n}{\overset{{\Omega }_n}{\int }}{\displaystyle \underset{D_S}{\int }\left\{Q_1\left(x,y,\eta ,\omega \right)+Q_2\left(x,y,\eta ,\omega \right)\right\}\text{d}{v}_{\eta }\omega \text{d}\omega }}\end{array}$ (29)

where we put

$Q_1\left(x,y,\eta ,\omega \right)=G_1\left(y,\eta ,\omega \right){\text{div}}_{\eta }\left\{\frac{1}{\epsilon \left(\eta \right)}{\text{grad}}_{\eta }{G}_2\left(x,\eta ,\omega \right)\right\}$

and

$Q_2\left(x,y,\eta ,\omega \right)=\frac{1}{\epsilon \left(\eta \right)}\text{grad}{}_{\eta }G{}_1\left(y,\eta ,\omega \right)\cdot {\text{grad}}_{\eta }{G}_2\left(x,\eta ,\omega \right)$

Now reconsider the equations satisfied by $G_1$ and $G_2$ (see (5a) and (7a)) to write

$\begin{array}{c}{Q}_1\left(x,y,\eta ,\omega \right)=-\left[{\omega }^2\mu \left(\eta \right)G_2\left(\eta ,x,\omega \right)+\delta \left(\eta -x\right)\right]G_1\left(y,\eta ,\omega \right)\\ =G_2\left(\eta ,x,\omega \right){\text{div}}_{\eta }\left[\frac{1}{\epsilon \left(\eta \right)}{\text{grad}}_{\eta }{G}_1\left(\eta ,y,\omega \right)\right]\\ +\delta \left(\eta -y\right)G_2\left(\eta ,x,\omega \right)-\delta \left(\eta -x\right)G_1\left(y,\eta ,\omega \right).\end{array}$

Now let us add the left hand-side again to both sides to obtain a more symmetrical expression

$\begin{array}{c}2Q_1=G_2\left(\eta ,x,\omega \right){\text{div}}_{\eta }\left[\frac{1}{\epsilon \left(\eta \right)}{\text{grad}}_{\eta }{G}_1\left(\eta ,y,\omega \right)\right]\\ +G_1\left(y,\eta ,\omega \right){\text{div}}_{\eta }\left[\frac{1}{\epsilon \left(\eta \right)}{\text{grad}}_{\eta }{G}_2\left(\eta ,x,\omega \right)\right]\\ +\delta \left(\eta -y\right)G_2\left(\eta ,x,\omega \right)-\delta \left(\eta -x\right)G_1\left(y,\eta ,\omega \right).\end{array}$

By using the latter in (29) one obtains

$P_n\left(y,x\right)=P_{n1}\left(y,x\right)+P_{n2}\left(y,x\right)$

where

$P_{n1}\left(y,x\right)=-\frac{i}{2\epsilon \left(y\right)}{\displaystyle \underset{-{\Omega }_n}{\overset{{\Omega }_n}{\int }}\left\{G_2\left(y,x,\omega \right)-G_1\left(y,x,\omega \right)\right\}\omega \text{d}\omega }$

and

$P_{n2}\left(y,x\right)=-\frac{i}{2\epsilon \left(y\right)}{\displaystyle \underset{-{\Omega }_n}{\overset{{\Omega }_n}{\int }}{\displaystyle \underset{D_S}{\int }\left\{Q_1+2Q_2+Q_{{}_3}\right\}\text{d}{v}_{\eta }\omega \text{d}\omega }}$ (30)

with

$Q_3=Q_1\left(x,y,\eta ,-\omega \right)$ .

In what follows we will show that if $D_S$ is convex (see Sec.4 below), then

$\underset{n\to \infty }{\lim }{P}_{n2}\left(y,x\right)=0$ (31)

which means that the contribution of $P_{n2}$ to $q_n\left(x\right)$ is naught when $n\to \infty$ . As to the contribution of P_n1, one has

$q_n\left(x\right)\to \frac{1}{2\text{π}}\underset{n\to \infty }{\lim }{\displaystyle \underset{V}{\int }{p}_0\left(y\right)P_{n1}\left(y,x\right)\text{d}{v}_y}\equiv q\left(x\right)$ .

This equation defines a distribution (generalized function), say $f\left(y,x\right)$ , which permits us to write [2]

$q\left(x\right)=\langle f\left(y,x\right),p_0\left(y\right)\rangle$ (32a)

and

$n\to \infty \Rightarrow \frac{1}{2\text{π}}{P}_{n1}\left(y,x\right)\to f\left(y,x\right).$

On the other hand, from Theorem-4 we already know that

$n\to \infty \Rightarrow \frac{1}{2\text{π}}{P}_n{}_1\left(y,x\right)=-\frac{1}{2\epsilon \left(y\right)\mu \left(y\right)}\delta \left(x-y\right)$ ,

which reduces (32a) to

$q\left(x\right)=-\frac{p_0\left(x\right)}{2\epsilon \left(x\right)\mu \left(x\right)}=-\frac{1}{2}p\left(x,+0\right).$ (32b)

The second equality in (32b) is an obvious issue of the distributional validity of (25).

(32b) is the basic formula of the present work. It can be used in photo-acous- tic and thermo-acoustic tomography problems connected with non-homoge- neous media. To obtain a more explicit expression of $p\left(x,+0\right)$ , let us insert the expression of $q_n\left(x\right)$ given by (26) into (27), and make the change of variable $\omega ={\Omega }_n\sin \alpha$ . If we use also (12b), then we write successively

$\begin{array}{c}{q}_n\left(x\right)={\displaystyle \underset{-\text{π}/2}{\overset{\text{π}/2}{\int }}{\displaystyle \underset{S}{\int }\frac{{\Omega }_n\cos \alpha }{{\epsilon }_0}\hat{p}\left(\eta ,{\Omega }_n\sin \alpha \right)}}\frac{\partial }{\partial n_{\eta }}{G}_2\left(x,\eta ,{\Omega }_n\sin \alpha \right)\text{d}{S}_{\eta }\text{d}\alpha \\ ={\displaystyle \underset{-\text{π}/2}{\overset{\text{π}/2}{\int }}{\displaystyle \underset{S}{\int }\frac{{\Omega }_n\cos \alpha }{{\epsilon }_0}\hat{p}\left(\eta ,{\Omega }_n\sin \alpha \right)}}\frac{\partial }{\partial n_{\eta }}\left[{\text{e}}^{-i{\Omega }_n{R}_1\sin \alpha /c_0}\right]\left\{\frac{{\epsilon }_0}{4\text{π}{R}_1}+O\left(\frac{1}{{\Omega }_n}\right)\right\}\text{d}{S}_{\eta }\text{d}\alpha .\end{array}$