Finding a source inside a sphere

Sixty years ago, Wilson and Bayley [1] began with an assumption 'that the electric field generated by the heart may be regarded as not significantly different from that of a current dipole at the center of a homogeneous spherical conductor', and then went on to examine the effect of moving the dipole away from the center. Point dipoles inside spheres (and ellipsoids) are also used in the context of brain imaging; see Dassios [2] for a recent review and the book by Ammari [3]. Locating point sources and dipoles using surface measurements is an example of an inverse source problem [4].

The basic static problem consists of Laplace's equation in a bounded region V_i with boundary ∂V. There is a source of some kind in V_i and the goal is to identify the source from Cauchy data on ∂V. Uniqueness [5] and stability [6, 7] results are available. For some numerical results, see [8–11] in two dimensions and [12–16] in three dimensions. We will show that some of these problems for balls can be solved exactly by analytical methods. Exact methods for magnetostatic problems have been devised previously [17, 18]; our method is arguably more elementary.

Analogous problems for the Helmholtz equation can be posed [17, 19]. In particular, the problem of low-frequency internal source excitation of a homogeneous sphere has applications in electroencephalography (EEG). The basic principle of EEG lies in the detection and processing of a signal generated by neural activity in order to map certain brain functions. More precisely, the primary source inside the brain is determined from a set of signals measured on the scalp, thus generating electric brain images [20]. The frequency, f, of the measured signal from a human brain is very low and hence the interior excitation of a spherical human head by a low-frequency point source constitutes a suitable EEG model. For example [21], k_ia ≃ 1.3 × 10⁻⁷ for f = 60 Hz and head radius a = 10 cm, where k_i is the interior wavenumber. Further applications arise in magnetic resonance imaging (MRI) [22], brain electrical impedance tomography (EIT) [23], the modeling of antennas implanted inside the head for hyperthermia or biotelemetry [24] and in studies of the operation of wireless devices around sensitive medical equipment (such as cardiac pacemakers) [25].

The direct problem of interior point-source excitation of a layered sphere has been solved in [26] for acoustics and in [27] for electromagnetics. These papers include approximations of the far field in the low-frequency regime and related far-field inverse scattering algorithms. The possibility of using a near-field quantity, namely the scattered field at the source point, in order to obtain inverse scattering algorithms for a small homogeneous soft sphere has been pointed out in [28]; these algorithms are not applicable here as we regard the source point as being inaccessible. For other implementations of near-field inverse problems, see also [29] and [30, p 133].

In this paper, we suppose that there is a point source inside a sphere. There are fields both inside and outside the sphere, with appropriate interface conditions on the sphere. The inverse problem is to determine the location and strength of the source knowing the (total) field on the sphere. For a point dipole, the orientation of the dipole is also to be determined. The internal characteristics (such as wavenumber or conductivity) are also to be found. For static problems, exact and complete results are obtained. Similar results are obtained for acoustic problems, although further approximations are used. The static problem with several point sources is also considered.

Our analytical inversion algorithms make use of the moments obtained by integrating the product of the total field on the spherical interface with spherical harmonic functions. All the information about the primary source and the sphere's interior characteristics is encoded in these moments.

We emphasise that our method is simple, explicit and exact (given exact data). However, as might be expected, it is limited to simple problems (with some extensions mentioned in section 5). We observe that all the inverse problems mentioned above are finite-dimensional, in the sense that the goal is to determine a few numbers (such as the coordinates and strength of the unknown source), not functions. As such, these problems are essentially simpler than many other inverse problems (such as determining the shape of a scatterer), and so it is perhaps not surprising that some of them can be solved analytically. It does not seem to have been noticed that locating a source inside a sphere is one of those problems.

Consider a homogeneous spherical object of radius a, surrounded by an infinite homogeneous medium. Denote the exterior by V_e and the interior by V_i. There is an interior point source of some kind at an unknown location $\mathbf {r}_1\in V_{{\rm i}}$. We are interested in characterizing the source, using information on the spherical interface.

Denote the field outside the sphere by u_e and the (total) field inside by u_i. Then, we can write u_i = u^pr + u^sec, where u^pr is the primary field due to the source (u^pr is singular at $\mathbf {r}_1$) and u^sec is the secondary (regular) field. The field u_e is regular and satisfies an appropriate far-field condition. The fields u_e and u_i are related by continuity (transmission) conditions across the spherical interface.

Introduce spherical polar coordinates (r, θ, ϕ) for the point at $\mathbf {r}$ so that the source is at (r₁, θ₁, ϕ₁) with $r_1=|\mathbf {r}_1|<a$. Then, the transmission conditions are

$$\begin{equation} u_{{\rm e}}=u_{{\rm i}}\quad \mbox{and}\quad \frac{1}{\rho _{{\rm e}}}\frac{\partial u_{{\rm e}}}{\partial r}= \frac{1}{\rho _{{\rm i}}}\frac{\partial u_{{\rm i}}}{\partial r} \quad \mbox{at $r=a$,} \end{equation} \tag{ 1 }$$

where ρ_e and ρ_i are constants. To complete the problem formulation, we must specify the governing differential equations, so we separate into electrostatics and acoustics.

2.1. Static problems

For static problems, both u_e and u^sec are governed by Laplace's equation. The field u_e decays to zero at infinity. In the context of electrostatics, ρ_e and ρ_i are inverse conductivities.

In some applications, the exterior is non-conducting. Then, u_i solves an interior Neumann problem (with (1) replaced by ∂u_i/∂r = 0 on r = a), and u_e solves an exterior Dirichlet problem, with u_e = u_i on r = a.

For the primary field, we could choose a point source

$$\begin{equation} u^{{\rm pr}}(\mathbf {r}; \mathbf {r}_1)=\frac{A}{|\mathbf {r}-\mathbf {r}_1|},\qquad \mathbf {r}\in \mathbb {R}^{3}\backslash \lbrace \mathbf {r}_1\rbrace , \end{equation} \tag{ 2 }$$

where A is a real constant; see, for example, [31, p 49]. Dipole fields will also be considered; see (11).

2.2. Acoustic problems

For time-harmonic acoustic problems, we consider compressible fluids so that ρ_e and ρ_i are densities. The governing equations are Helmholtz equations,

$$\begin{equation*} \big(\nabla ^2+k_{{\rm e}}^2\big)u_{{\rm e}}=0\quad \mbox{in $V_{{\rm e}}$,}\qquad \big(\nabla ^2+k_{{\rm i}}^2\big)u^{\rm sec}=0\quad \mbox{in $V_{{\rm i}}$,} \end{equation*}$$

where k_e and k_i are the respective wavenumbers. The physical fields are given by Re{u_e e^−iωt}, for example; henceforth, we suppress the time dependence. The complex-valued field u_e must satisfy the Sommerfeld radiation condition at infinity. For the primary field, we take

$$\begin{equation} u^{{\rm pr}}(\mathbf {r}; \mathbf {r}_1)= \frac{A}{|\mathbf {r}-\mathbf {r}_1|}\,\exp {(\mathrm{i}k_{{\rm i}}|\mathbf {r}-\mathbf {r}_1|)}, \quad \mathbf {r}\in \mathbb {R}^{3}\backslash \lbrace \mathbf {r}_1\rbrace , \end{equation} \tag{ 3 }$$

generated by a point source with the position vector $\mathbf {r}_1$, where A is a (possibly complex) constant; see, for example, [31, p 144] or [32, p 153].

There is a static source of some kind at $\mathbf {r}_1$ generating the field u^pr. Near the sphere (r₁ < r < a), separation of variables gives the expansion

$$\begin{equation*} u^{{\rm pr}}(\mathbf {r}; \mathbf {r}_1) =\sum _{n=0}^\infty \sum _{m=-n}^n f_n^m(\mathbf {r}_1)(a/r)^{n+1}Y_n^m(\hat{\mathbf {r}}), \end{equation*}$$

where $\hat{\mathbf {r}}=\mathbf {r}/|\mathbf {r}|=(\sin \theta \cos \phi ,\,\sin \theta \sin \phi ,\,\cos \theta )$ and $Y_n^m(\hat{\mathbf {r}})=Y_n^m(\theta ,\phi )$ is a spherical harmonic; we define Y^m_n as in [32, section 3.2].

The quantities f^m_n characterize the source. For a point source, defined by (2),

$$\begin{equation} f_n^m(\mathbf {r}_1)=\frac{4\pi A}{a} \frac{(-1)^m}{2n+1}(r_1/a)^n Y_n^{-m}(\hat{\mathbf {r}}_1); \end{equation} \tag{ 4 }$$

see, for example, [33, equation (3.70)].

The total field inside the sphere is u_i = u^pr + u^sec where

$$\begin{equation*} u^{{\rm sec}}(\mathbf {r})= \sum _{n=0}^\infty \sum _{m=-n}^n \alpha _n f_n^m(\mathbf {r}_1)(r/a)^{n}Y_n^m(\hat{\mathbf {r}}), \quad 0\le r<a \end{equation*}$$

(u^sec must be finite at r = 0) whereas the field outside is given by

$$\begin{equation*} u_{{\rm e}}(\mathbf {r})= \sum _{n=0}^\infty \sum _{m=-n}^n \beta _n f_n^m(\mathbf {r}_1)(a/r)^{n+1}Y_n^m(\hat{\mathbf {r}}), \quad r>a \end{equation*}$$

(u_e must vanish as r → ∞). The transmission conditions at r = a, (1), give

$$\begin{equation*} 1+\alpha _n=\beta _n, \quad -n-1+n\alpha _n=-\varrho (n+1)\beta _n, \end{equation*}$$

where ϱ = ρ_i/ρ_e. These can be solved for α_n and β_n yielding

$$\begin{equation} \alpha _n=\frac{(1-\varrho )(n+1)}{n+\varrho (n+1)}, \qquad \beta _n=\frac{2n+1}{n+\varrho (n+1)}. \end{equation} \tag{ 5 }$$

Note that these quantities do not depend on any characteristics of the source. Also, for ϱ = 1 (no interface at r = a), we can verify that $u_{{\rm e}}(\mathbf {r})=u^{{\rm pr}}(\mathbf {r};\mathbf {r}_1)$ and $u^{\sec }(\mathbf {r})=0$, as expected.

The field on the sphere is

$$\begin{equation*} u_{{\rm surf}}(\theta ,\phi )= \sum _{n=0}^\infty \sum _{m=-n}^n \frac{2n+1}{n+\varrho (n+1)} f_n^m(\mathbf {r}_1)Y_n^m(\theta ,\phi ). \end{equation*}$$

It is this quantity that we will use to find the source.

The spherical harmonics are orthonormal:

$$\begin{equation*} \int _\Omega Y_n^m\overline{Y_\nu ^\mu } \,\mathrm{d}\Omega =\int _0^\pi \int _{-\pi }^\pi Y_n^m(\theta ,\phi ) \overline{Y_\nu ^\mu (\theta ,\phi )} \sin \theta \,\mathrm{d}\phi \,\mathrm{d}\theta =\delta _{n\nu }\delta _{m\mu }, \end{equation*}$$

where Ω is the unit sphere and the overbar denotes complex conjugation. Hence, the moments

$$\begin{equation} M_n^m\equiv \frac{1}{\sqrt{4\pi }} \int _\Omega u_{{\rm surf}} \overline{Y_n^m}\,\mathrm{d}\Omega =\frac{1}{\sqrt{4\pi }}\frac{2n+1}{n+\varrho (n+1)} f_n^m(\mathbf {r}_1) \end{equation} \tag{ 6 }$$

are known, in principle, if u is known on r = a; the double integral over Ω could be approximated using a suitable quadrature rule and corresponding point evaluations of u_surf. The problem now is to determine properties of the source and the interior material (namely, ρ_i = ρ_eϱ) from M^m_n.

3.1. Point source

For a point source, (4) and (6) give

$$\begin{equation} M_n^m=(-1)^m\frac{\tilde{A}\tilde{r}_1^n \sqrt{4\pi }}{n+\varrho (n+1)} Y_n^{-m}(\theta _1,\phi _1),\quad \mbox{with}\quad \tilde{A}=\frac{A}{a}\ \mbox{and}\ \tilde{r}_1=\frac{r_1}{a}. \end{equation} \tag{ 7 }$$

Thus, there are five unknowns, $\tilde{A}$, ϱ, $\tilde{r}_1$, θ₁ and ϕ₁.

As Y⁰₀ = (4π)^−1/2, we obtain

$$\begin{equation*} M_0^0=\tilde{A}/{\varrho }. \end{equation*}$$

This ratio is all that can be recovered if the source is at the sphere's center ($r_1=\tilde{r}_1=0$). So, let us assume now that $\tilde{r}_1\ne 0$.

For n = 1, we can use [32, equation (8.28)]

$$\begin{equation} \fl Y_1^0=\sqrt{\frac{3}{4\pi }} \cos \theta ,\qquad Y_1^1=-\sqrt{\frac{3}{8\pi }} \,\mathrm{e}^{\mathrm{i}\phi }\sin \theta ,\qquad Y_1^{-1}=\sqrt{\frac{3}{8\pi }} \mathrm{e}^{-\mathrm{i}\phi }\sin \theta \end{equation} \tag{ 8 }$$

giving

$$\begin{equation*} M_1^0=\tilde{A}\frac{\tilde{r}_1\sqrt{3}}{1+2\varrho } \cos \theta _1,\qquad M_1^{\pm 1}=\mp \tilde{A}\frac{\tilde{r}_1\sqrt{3/2}}{1+2\varrho } \mathrm{e}^{\mp \mathrm{i}\phi _1}\sin \theta _1. \end{equation*}$$

If M^±1₁ = 0, then θ₁ = 0 or π (the source is on the z-axis and so ϕ₁ is irrelevant); to decide which, note that the sign of M⁰₀M₁⁰ is the sign of cos θ₁. If M^±1₁ ≠ 0, ϕ₁ is determined by noting that the complex number M⁰₀M₁⁻¹ has argument ϕ₁.

If M⁰₁ = 0, then θ₁ = π/2. If M⁰₁ ≠ 0, $\sqrt{2}\,M_1^{-1}/M_1^0=\mathrm{e}^{\mathrm{i}\phi _1}\tan \theta _1$ determines θ₁. Also

$$\begin{equation} (1+2\varrho )^2\big\lbrace \big(M_1^0\big)^2-2M_1^1M_1^{-1}\big\rbrace =3(\tilde{A}\tilde{r}_1)^2. \end{equation} \tag{ 9 }$$

To conclude, we take a measurement with n = 2. Thus,

$$\begin{equation} \varrho (2+3\varrho )M_0^0M_2^m =(\tilde{A}\tilde{r}_1)^2(-1)^m\sqrt{4\pi } Y_2^{-m}(\theta _1,\phi _1). \end{equation} \tag{ 10 }$$

Then, choosing m such that Y^−m₂(θ₁, ϕ₁) ≠ 0 (which means take m = 0 unless P₂(cos θ₁) = 0), eliminate $(\tilde{A}\tilde{r}_1)^2$ between (9) and (10) to give a quadratic equation for ϱ (which is real and positive). Then, $\tilde{A}=A/a=M_0^0\varrho$ and $\tilde{r}_1=r_1/a$ follows from M^m₁ or from (9).

3.2. Point dipole

Let $\hat{{\bf d}}=(\sin \theta _{{\rm d}}\cos \phi _{{\rm d}},\,\sin \theta _{{\rm d}}\sin \phi _{{\rm d}},\,\cos \theta _{{\rm d}})$ be an unknown unit vector, giving the direction of a point dipole at $\mathbf {r}_1$. The field generated is given by

$$\begin{equation} u_{{\mathbf {r}_1}}^{{\rm pr}}(\mathbf {r})=\hat{{\bf d}}\cdot {\rm grad}_1 \left(A/|\mathbf {r}-\mathbf {r}_1|\right), \end{equation} \tag{ 11 }$$

where grad₁ denotes the gradient with respect to $\mathbf {r}_1$ and A is a real constant. As changing the sign of A is equivalent to changing the direction of $\hat{{\bf d}}$, we can assume that A is positive.

From (2) and (4), the source coefficients are given by

$$\begin{equation} f_n^m(\mathbf {r}_1)=\hat{{\bf d}}\cdot {\rm grad}_1\left( \frac{4\pi A}{a} \frac{(-1)^m}{2n+1}(r_1/a)^n Y_n^{-m}(\hat{\mathbf {r}}_1) \right), \end{equation} \tag{ 12 }$$

and then M^m_n is given by (6).

We have f⁰₀ = 0. Then, with n = 1, we have

$$\begin{equation*} \frac{a^2}{\sqrt{4\pi }}f_1^0(\mathbf {r})=\hat{{\bf d}}\cdot {\rm grad}\left( \frac{A}{\sqrt{3}} z\right)=\frac{A}{\sqrt{3}}\cos \theta _{{\rm d}} \end{equation*}$$

and

$$\begin{eqnarray*} &&\fl \frac{a^2}{\sqrt{4\pi }} f_1^{\pm 1} =\hat{{\bf d}}\cdot {\rm grad}\left( \frac{\sqrt{4\pi } A}{3} (-r)Y_1^{\mp 1}\right) \\ &&=\hat{{\bf d}}\cdot {\rm grad}\left( \frac{A}{\sqrt{6}}(\mathrm{i}y\mp x)\right)= \mp \frac{A}{\sqrt{6}}\, \mathrm{e}^{\mp \mathrm{i}\phi _{{\rm d}}}\sin \theta _{{\rm d}}. \end{eqnarray*}$$

If M^±1₁ = 0, then θ_d = 0 or π, which we can determine by examining the sign of M⁰₁, recalling that A > 0. If M^±1₁ ≠ 0, ϕ_d is determined by noting that M⁻¹₁ is a complex number with argument ϕ_d.

If M⁰₁ = 0, then θ_d = π/2. If M⁰₁ ≠ 0, $\sqrt{2}M_1^{-1}/M_1^0 =\mathrm{e}^{\mathrm{i}\phi _{{\rm d}}}\tan \theta _{{\rm d}}$ determines θ_d.

Thus, we have determined the orientation of the dipole, $\hat{{\bf d}}$. Also,

$$\begin{equation} (1+2\varrho )^2a^4\big\lbrace \big(M_1^0\big)^2-2M_1^1M_1^{-1}\big\rbrace =3A^2. \end{equation} \tag{ 13 }$$

To find the dipole's location and strength, we move on to n = 2. Thus,

$$\begin{eqnarray*} &&\fl r^2Y_2^{\pm 2}=\frac{3\sqrt{5}}{\sqrt{4\pi }\sqrt{4!}} (x\pm \mathrm{i}y)^2, \quad r^2Y_2^{\pm 1}=\mp \frac{3\sqrt{5}}{\sqrt{4\pi }\sqrt{3!}} (x\pm \mathrm{i}y)z, \\ &&\fl r^2Y_2^0=\frac{\sqrt{5}}{\sqrt{4\pi }} \left(z^2-\frac{1}{2}(x^2+y^2)\right). \end{eqnarray*}$$

Hence,

$$\begin{eqnarray} &&\fl \frac{a^3}{\sqrt{4\pi }}f_2^0(\mathbf {r}_1) =\hat{{\bf d}}\cdot {\rm grad}_1\left( \frac{A}{2\sqrt{5}} \big(2z_1^2-x_1^2-y_1^2\big)\right) \nonumber \\ &&=\frac{Ar_1}{\sqrt{5}} \left\lbrace 2\cos \theta _1\cos \theta _{{\rm d}}-\sin \theta _1\sin \theta _{{\rm d}}\cos {(\phi _1-\phi _{{\rm d}})}\right\rbrace , \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} &&\fl \frac{a^3}{\sqrt{4\pi }}f_2^{\pm 1}(\mathbf {r}_1) =\mp \hat{{\bf d}}\cdot {\rm grad}_1\left( \frac{3A}{\sqrt{30}} (x_1\mp \mathrm{i}y_1)z_1 \right) \nonumber \\ &&=\mp \frac{3Ar_1}{\sqrt{30}}\lbrace \mathrm{e}^{\mp \mathrm{i}\phi _{{\rm d}}}\cos \theta _1\sin \theta _{{\rm d}}+\mathrm{e}^{\mp \mathrm{i}\phi _1}\sin \theta _1\cos \theta _{{\rm d}}\rbrace , \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} &&\fl \frac{a^3}{\sqrt{4\pi }}f_2^{\pm 2}(\mathbf {r}_1) =\hat{{\bf d}}\cdot {\rm grad}_1\left( \frac{3A}{2\sqrt{30}}(x_1\mp \mathrm{i}y_1)^2 \right) \nonumber \\ &&=\frac{3Ar_1}{\sqrt{30}}\,\mathrm{e}^{\mp \mathrm{i}(\phi _1+\phi _{{\rm d}})} \sin \theta _1\sin \theta _{{\rm d}}. \end{eqnarray} \tag{ 16 }$$

Our first goal is to determine θ₁ and ϕ₁ from M^m₂. Note that if all the moments M^m₂ vanish, then r₁ = 0 (the dipole is at the sphere's center); henceforth, we assume that r₁ > 0.

Let us first dispose of special cases. Suppose that θ_d = 0, in which case

$$\begin{equation*} \frac{a^3}{\sqrt{4\pi }}f_2^0(\mathbf {r}_1) =\frac{2Ar_1}{\sqrt{5}} \cos \theta _1,\qquad \frac{a^3}{\sqrt{4\pi }}f_2^{\pm 1}(\mathbf {r}_1) =\mp \frac{3Ar_1}{\sqrt{30}} \mathrm{e}^{\mp \mathrm{i}\phi _1}\sin \theta _1 \end{equation*}$$

and f^±2₂ = 0. If M^±1₂ = 0, then θ₁ = 0 or π, which the sign of M⁰₂ determines. If M^±1₂ ≠ 0, ϕ₁ is the argument of M⁻¹₂. Then, if M⁰₂ = 0, θ₁ = π/2, otherwise M^±1₂/M₂⁰ determines θ₁. A similar analysis succeeds when θ_d = π.

When θ_d = π/2,

$$\begin{equation*} \fl \frac{a^3}{\sqrt{4\pi }}f_2^{\pm 1}(\mathbf {r}_1) =\mp \frac{3Ar_1}{\sqrt{30}} \mathrm{e}^{\mp \mathrm{i}\phi _{{\rm d}}}\cos \theta _1, \qquad \! \frac{a^3}{\sqrt{4\pi }}f_2^{\pm 2}(\mathbf {r}_1) =\frac{3Ar_1}{\sqrt{30}} \mathrm{e}^{\mp \mathrm{i}(\phi _1+\phi _{{\rm d}})}\sin \theta _1. \end{equation*}$$

So, if M^±2₂ = 0, θ₁ = 0 or π, and the sign of $\mp \mathrm{e}^{\pm \mathrm{i}\phi _{{\rm d}}}M_2^{\pm 1}$ gives the sign of cos θ₁. If M^±2₂ ≠ 0, $\mathrm{e}^{-\mathrm{i}\phi _{{\rm d}}}M_2^{-2}$ has the argument ϕ₁. Then, if M^±1₂ = 0, θ₁ = π/2, otherwise M^±2₂/M₂^±1 determines θ₁.

Now, for θ_d different than 0, π/2 and π, we proceed as follows. If M^±2₂ = 0, then θ₁ = 0 or π, and the sign of $\mp \mathrm{e}^{\pm \mathrm{i}\phi _{{\rm d}}}M_2^{\pm 1}$ gives the sign of cos θ₁. If M^±2₂ ≠ 0, the angle ϕ₁ is determined from (16) by noting that $\mathrm{e}^{-\mathrm{i}\phi _{{\rm d}}}M_2^{-2}$ has the argument ϕ₁. Then, M¹₂/M₂² gives $\cot \theta _1$.

Thus, we have now calculated the angular coordinates of the dipole, θ₁ and ϕ₁, as well as its orientation, θ_d and ϕ_d. It remains to determine A, r₁ and ϱ. Now, from (6) and (12), any non-zero M^m₁ gives A/(1 + 2ϱ) = m₁, say. Similarly, any non-zero M^m₂ gives Ar₁/(2 + 3ϱ) = m₂, say. So, in order to extract all three unknowns, we have to move on to n = 3. (For the point-source problem, discussed in section 3.1, we had non-trivial information from n = 0, so we did not require measurements with n = 3.) Then, any non-zero M^m₃ gives Ar²₁/(3 + 4ϱ) = m₃, say. Eliminating A and r₁ gives

$$\begin{equation*} (1+2\varrho )(3+4\varrho )m_1m_3-(2+3\varrho )^2m_2^2=0, \end{equation*}$$

a quadratic for ϱ. Then, m₁ gives A and m₂ gives r₁. Of course, the details of the calculation depend on the choices of m in M^m_n. If P₃(cos θ₁) ≠ 0, we would select M⁰₃ and use (see [32, section 3.4])

$$\begin{equation*} r^3Y_3^0=\frac{\sqrt{7}}{\sqrt{4\pi }} \left(z^3-\frac{3}{2}z(x^2+y^2)\right), \end{equation*}$$

giving

$$\begin{equation*} \frac{a^4}{\sqrt{4\pi }}f_3^{0}(\mathbf {r}_1) =\hat{{\bf d}}\cdot {\rm grad}_1\left( \frac{A}{\sqrt{7}}\left(z_1^3-\frac{3}{2}z_1(x_1^2+y_1^2)\right) \right) =\frac{3Ar_1^2}{\sqrt{7}}\mathcal {F}, \end{equation*}$$

where $\mathcal {F}= -\sin \theta _{{\rm d}}\cos \theta _1\sin \theta _1\cos (\phi _{{\rm d}}-\phi _1) +\cos \theta _{{\rm d}}(\cos \theta _1-\frac{1}{2}\sin ^2\theta _1)$. Hence,

$$\begin{equation} M_3^0=\frac{\sqrt{4\pi }}{a^4}\frac{3\sqrt{7}}{3+4\varrho }A\,r_1^2\,\mathcal {F}. \end{equation} \tag{ 17 }$$

This completes the determination of all the parameters of the problem.

3.3. N point sources

Suppose there are N point sources, located at $\mathbf {r}_j=(x_j,y_j,z_j)$ with spherical polar coordinates (r_j, θ_j, ϕ_j), j = 1, 2, ..., N. Suppose for simplicity that each source has the same strength A, and that ϱ is known. Thus, there are 3N + 1 unknowns. By linearity and (7), the moments are

$$\begin{equation*} M_n^m=\frac{A (-1)^m \sqrt{4\pi }}{a^{n+1}\lbrace n+\varrho (n+1)\rbrace } \sum _{j=1}^N r_j^nY_n^{-m}(\theta _j,\phi _j). \end{equation*}$$

The moment M⁰₀ determines A.

Now, rⁿY^−m_n(θ, ϕ) is a homogeneous polynomial of degree n in x, y and z [30, section 2.3]. For each n, there are 2n + 1 distinct polynomials. So, if we collect all the moments from n = 1 to n = N_M, we will have N²_M + 2N_M pieces of data from which to recover 3N unknowns. Thus, in principle, we can recover the N locations by solving a system of polynomial equations: exact solutions will be rare.

For two point sources (N = 2), we could use the axisymmetric spherical harmonics (m = 0) in order to determine r₁, r₂, z₁ and z₂; writing down M⁰_n for n = 1, 2, 3 and 4 gives enough equations. For example, ignoring multiplicative constants, rY⁰₁ = z, M⁰₁ = z₁ + z₂, r²Y⁰₂ = 3z² − r², M⁰₂ = 3z₁² + 3z²₂ − r₁² − r²₂, and so on. In this way, we can find a quartic equation for each unknown. We can then use other moments (with m ≠ 0) to recover ϕ₁ and ϕ₂. Evidently, this process becomes more complicated as N is increased.

In this section, we outline how some of the static results can be generalized to acoustic problems. We suppose that there is a point source at $\mathbf {r}_1$ generating the field u^pr defined by (3). Near the sphere (r₁ < r < a), we have the expansion [32, theorem 6.4]

$$\begin{equation*} \fl u^{{\rm pr}}(\mathbf {r};\mathbf {r}_1) =\sum _{n=0}^\infty \sum _{m=-n}^n d_n^m(\mathbf {r}_1)\psi _n^m(k_{{\rm i}},\mathbf {r}), \qquad d_n^m=4\pi \mathrm{i}k_{{\rm i}}A (-1)^m\hat{\psi }_n^{-m}(k_{{\rm i}},\mathbf {r}_1), \end{equation*}$$

where $\psi _n^m(k,\mathbf {r})=h_n(kr)Y_n^m(\hat{\mathbf {r}})$, $\hat{\psi }_n^m(k,\mathbf {r})=j_n(kr)Y_n^m(\hat{\mathbf {r}})$, j_n is a spherical Bessel function and h_n is a spherical Hankel function of the first kind.

The total field inside the sphere is u^pr + u^sec where

$$\begin{equation*} u^{{\rm sec}}(\mathbf {r})= \sum _{n=0}^\infty \sum _{m=-n}^n \alpha _n d_n^m(\mathbf {r}_1)\hat{\psi }_n^m(k_{{\rm i}},\mathbf {r}), \quad 0\le r<a, \end{equation*}$$

whereas the (radiating) field outside is given by

$$\begin{equation*} u_{{\rm e}}(\mathbf {r})= \sum _{n=0}^\infty \sum _{m=-n}^n \beta _n d_n^m(\mathbf {r}_1)\psi _n^m(k_{{\rm e}},\mathbf {r}), \quad r>a. \end{equation*}$$

The transmission conditions at r = a give

$$\begin{equation*} h_n(\kappa _{{\rm i}})+\alpha _nj_n(\kappa _{{\rm i}})=\beta _nh_n(\kappa _{{\rm e}}), \quad h_n^{\prime }(\kappa _{{\rm i}})+\alpha _nj_n^{\prime }(\kappa _{{\rm i}})=\beta _nwh_n^{\prime }(\kappa _{{\rm e}}), \end{equation*}$$

where κ_i = k_ia, κ_e = k_ea and w = (k_eρ_i)/(k_iρ_e). These can be solved for α_n and β_n. (Note again that these quantities do not depend on properties of the source.) For example,

$$\begin{equation} \beta _n=\mathrm{i}\kappa _{{\rm i}}^{-2}\lbrace wj_n(\kappa _{{\rm i}})h_n^{\prime }(\kappa _{{\rm e}}) -j_n^{\prime }(\kappa _{{\rm i}})h_n(\kappa _{{\rm e}})\rbrace ^{-1}. \end{equation} \tag{ 18 }$$

Then, the field on the sphere is

$$\begin{equation*} u_{{\rm surf}}(\theta ,\phi )= \sum _{n=0}^\infty \sum _{m=-n}^n \beta _n d_n^m(\mathbf {r}_1)h_n(\kappa _{{\rm e}})Y_n^m(\theta ,\phi ). \end{equation*}$$

Our inverse problem is to determine the location of the source (r₁, θ₁ and ϕ₁), the strength of the source (A) and the interior properties (κ_i and ρ_i) from u_surf. We assume that κ_e and ρ_e are known.

As the spherical harmonics are orthonormal,

$$\begin{eqnarray*} &&\fl M_n^m \equiv \frac{1}{\sqrt{4\pi } h_n(\kappa _{{\rm e}})} \int _\Omega u_{{\rm surf}} \overline{Y_n^m} \,\mathrm{d}\Omega =\beta _n \frac{d_n^m(\mathbf {r}_1)}{\sqrt{4\pi }} \\ &&=\mathrm{i}k_{{\rm i}}A\beta _n j_n(k_{{\rm i}}r_1)\sqrt{4\pi } (-1)^m Y_n^{-m}(\theta _1,\phi _1). \end{eqnarray*}$$

In principle, M^m_n are known if u is known on r = a. Let us write

$$\begin{equation} B_n=\mathrm{i}k_{{\rm i}}A\beta _n j_n(K)\quad \mbox{with}\quad K=k_{{\rm i}}r_1; \end{equation} \tag{ 19 }$$

these quantities do not depend on θ₁ or ϕ₁.

As Y⁰₀ = (4π)^−1/2, we obtain M⁰₀ = B₀.

For n = 1, we can use (8), giving

$$\begin{equation*} M_1^0=B_1\sqrt{3}\cos \theta _1,\qquad M_1^{\pm 1}=\mp B_1\sqrt{3/2} \,\mathrm{e}^{\mp \mathrm{i}\phi _1}\sin \theta _1. \end{equation*}$$

If all three of these vanish, B₁ = 0, and we move on to n = 2.

If M^±1₁ = 0 but M⁰₁ ≠ 0, then θ₁ = 0 or π, ϕ₁ is irrelevant and B₁ is undetermined.

Suppose now that M^±1₁ ≠ 0. Write B₁ = |B₁|e^iδ. Then, as sin θ₁ > 0, the complex number M⁻¹₁ has the argument δ + ϕ₁ and the complex number −M¹₁ has the argument δ − ϕ₁; thus, we can determine both δ and ϕ₁. For |B₁|, we can use

$$\begin{equation*} \big(M_1^0\big)^2-2M_1^1M_1^{-1}=3B_1^2=3|B_1|^2\mathrm{e}^{2\mathrm{i}\delta }. \end{equation*}$$

Finally, use M⁰₁ to deduce cos θ₁ and hence θ₁.

To make further progress, we make further assumptions. They are either that k_i is known or that k_ia and k_ea are small.

4.1. k_i is known

Suppose that k_i is known. Having already determined ϕ₁ and θ₁, we proceed to determine w, r₁ and A. By considering M⁰_n, we can obtain values for B_n. Then, the recurrence relation for spherical Bessel functions gives

$$\begin{equation*} \frac{1}{K}=\frac{j_{n-1}(K)+j_{n+1}(K)}{(2n+1)j_n(K)} =\frac{B_{n-1}/\beta _{n-1}+B_{n+1}/\beta _{n+1}}{(2n+1)B_n/\beta _n} \end{equation*}$$

for each n ⩾ 1. Equating two of these gives

$$\begin{equation*} \frac{B_{n-1}/\beta _{n-1}+B_{n+1}/\beta _{n+1}}{(2n+1)B_n/\beta _n} =\frac{B_{n}/\beta _{n}+B_{n+2}/\beta _{n+2}}{(2n+3)B_{n+1}/\beta _{n+1}}. \end{equation*}$$

Thus,

$$\begin{equation*} \frac{2n+3}{2n+1} \left( \frac{B_{n-1}}{\beta _{n-1}}+\frac{B_{n+1}}{\beta _{n+1}}\right) \frac{B_{n+1}}{\beta _{n+1}}= \left(\frac{B_{n}}{\beta _{n}}+\frac{B_{n+2}}{\beta _{n+2}}\right) \frac{B_n}{\beta _n}. \end{equation*}$$

For each n ⩾ 1, this is a quadratic equation for w because β⁻¹_n is linear in w; see (18). This equation has the form

$$\begin{equation} \mathcal {A}_nw^2+\mathcal {B}_nw+\mathcal {C}_n=0, \end{equation} \tag{ 20 }$$

where

$$\begin{eqnarray*} &&\fl \mathcal {A}_n= (2n+3)B_{n-1}B_{n+1} j_{n-1}(\kappa _{{\rm i}})h^{\prime }_{n-1}(\kappa _{{\rm e}})j_{n+1}(\kappa _{{\rm i}})h^{\prime }_{n+1}(\kappa _{{\rm e}})\\ &&\mbox{}-(2n+1)B_{n+2}B_{n} j_{n}(\kappa _{{\rm i}})h^{\prime }_{n}(\kappa _{{\rm e}})j_{n+2}(\kappa _{{\rm i}})h^{\prime }_{n+2}(\kappa _{{\rm e}})\\ &&\mbox{} +(2n+3)B^2_{n+1}\big[j_{n+1}(\kappa _{{\rm i}})h^{\prime }_{n+1}(\kappa _{{\rm e}})\big]^2 -(2n+1)B^2_{n}\big[j_{n}(\kappa _{{\rm i}})h^{\prime }_{n}(\kappa _{{\rm e}})\big]^2,\\ &&\fl \mathcal {B}_n= 2(2n+1)B^2_{n}j_{n}(\kappa _{{\rm i}})h^{\prime }_{n}(\kappa _{{\rm e}})j^{\prime }_{n}(\kappa _{{\rm i}})h_{n}(\kappa _{{\rm e}})\\ &&\mbox{} -2(2n+3)B^2_{n+1} j^{\prime }_{n+1}(\kappa _{{\rm i}})h_{n+1}(\kappa _{{\rm e}})j_{n+1}(\kappa _{{\rm i}})h^{\prime }_{n+1}(\kappa _{{\rm e}})\\ &&\mbox{} +(2n+1)B_{n+2}B_{n} \big [j_{n}(\kappa _{{\rm i}})h^{\prime }_{n}(\kappa _{{\rm e}})j^{\prime }_{n+2}(\kappa _{{\rm i}})h_{n+2}(\kappa _{{\rm e}}) \\ &&\mbox{} + j^{\prime }_{n}(\kappa _{{\rm i}})h_{n}(\kappa _{{\rm e}})j_{n+2}(\kappa _{{\rm i}})h^{\prime }_{n+2}(\kappa _{{\rm e}})\big ]\\ &&\mbox{} -(2n+3)B_{n-1}B_{n+1} \big [j_{n-1}(\kappa _{{\rm i}})h^{\prime }_{n-1}(\kappa _{{\rm e}})j^{\prime }_{n+1}(\kappa _{{\rm i}})h_{n+1}(\kappa _{{\rm e}}) \\ &&\mbox{} +j^{\prime }_{n-1}(\kappa _{{\rm i}})h_{n-1}(\kappa _{{\rm e}}) j_{n+1}(\kappa _{{\rm i}})h^{\prime }_{n+1}(\kappa _{{\rm e}})\big ], \\ &&\fl \mathcal {C}_n= (2n+3)B_{n-1}B_{n+1}j^{\prime }_{n-1}(\kappa _{{\rm i}})h_{n-1}(\kappa _{{\rm e}}) j^{\prime }_{n+1}(\kappa _{{\rm i}})h_{n+1}(\kappa _{{\rm e}})\\ &&\mbox{} -(2n+1)B_{n+2}B_{n}j^{\prime }_{n}(\kappa _{{\rm i}})h_{n}(\kappa _{{\rm e}})j^{\prime }_{n+2}(\kappa _{{\rm i}})h_{n+2}(\kappa _{{\rm e}})\\ &&\mbox{} +(2n+3)B^2_{n+1}\big[j^{\prime }_{n+1}(\kappa _{{\rm i}})h_{n+1}(\kappa _{{\rm e}})\big]^2 -(2n+1)B^2_{n}\big[j^{\prime }_{n}(\kappa _{{\rm i}})h_{n}(\kappa _{{\rm e}})\big]^2. \end{eqnarray*}$$

Although the parameter w (which gives the interior density ρ_i) solves the quadratic equation (20) for each integer n, the second solution may depend on n.

The distance r₁ can then be determined using β₀B₁j₀(k_ir₁) = β₁B₀j₁(k_ir₁) (see (19)), while the strength A follows from the expression for B₀.

4.2. k_ia and k_ea are small

Suppose that k_i is unknown. Then analytic progress can be made in the low-frequency zone, using the assumptions κ_i ≪ 1 and κ_e ≪ 1. In that case, the coefficients β_n and B_n have the following leading-order approximations as κ_i → 0 and κ_e → 0:

$$\begin{equation} \beta _n\sim \frac{(2n+1)(k_{{\rm e}}/k_{{\rm i}})^{n+1}}{n+\varrho (n+1)},\qquad B_n\sim \frac{\mathrm{i}A k_{{\rm e}}}{n+\varrho (n+1)}\frac{(k_{{\rm e}}r_1)^n}{c_n}, \end{equation} \tag{ 21 }$$

where c_n = 1 · 3⋅⋅⋅(2n − 1), with c₀ = 1. Then,

$$\begin{equation*} \frac{3B_2B_0}{B_1^2}=\frac{(2\varrho +1)^2}{\varrho (3\varrho +2)} \end{equation*}$$

determines ϱ = ρ_i/ρ_e; as expected, this is equivalent to the static formula, given below (10). Then, B₁/B₀ = ϱk_er₁/(2ϱ + 1) determines r₁ and ϱB₀ = ik_eA determines the strength A.

Evidently, we cannot determine k_i from the leading-order approximation to B_n, (21), because k_i does not appear. Thus, one either should use another kind of measurement, or one could use a higher-order approximation. For example,

$$\begin{equation*} B_0\sim \mathrm{i}A\frac{k_{{\rm e}}}{\varrho }\left( 1 -\frac{\kappa _{{\rm e}}^2}{2}+\frac{\kappa _{{\rm i}}^2}{6} +\frac{\kappa _{{\rm i}}^2}{3\varrho } +\frac{1}{6}(k_{{\rm i}}r_1)^2 \right) \end{equation*}$$

gives an estimate for k_i.

In this paper, we have considered simple inverse problems, locating sources and dipoles using surface data. These are examples of finite-dimensional inverse problems: the goal is to find the numerical values of a few parameters, such as the coordinates and strength of an interior source.

We have shown that the simplest problems (locating one source or one dipole inside a homogeneous sphere with exact data on the sphere) can be solved exactly.

Analytical solutions raise various questions. The first concerns inexact data. As our formulas are exact, the effect of errors in the data can be studied; for a detailed study of related methods, see [34]. Errors can come from the measurements of u_surf and from numerical integration over the unit sphere; the latter can be reduced by using accurate quadrature rules. Note that we use integrals of the surface data over the whole sphere: this is a limitation of our method. There are numerical algorithms that use data gathered over a piece of the sphere [16].

Next, we may consider extensions to other geometries, such as a layered sphere or an ellipsoid (motivated by EEG applications), or to elastic media: these extensions are feasible. We could use different data on the sphere, such as Neumann data (but note that, as we have solved a transmission problem, Neumann and Dirichlet data are related by solving an exterior boundary-value problem).

We have also seen (section 3.3) that increasing the number of sources leads to more complicated results. For example, if we want to locate two static sources, we have to determine nine numbers, namely, the coordinates and strength of each source and the interior density. Formally, we can make progress with this problem but we lose the attractive simplicity of the single source/dipole solution. Therefore, given that these multi-source inverse problems are finite-dimensional, it may then be better to abandon the analytical approach and resort to a numerical method.

For another class of problems, we could replace the point source by a spherical inclusion and then generate a primary field at r = a. The direct problems can be solved exactly, so the inverse problem could be tackled: it will be finite dimensional if the spherical inclusion has constant properties. Problems of this type have been discussed by Bonnetier and Vogelius [35]. We could also consider acoustic problems with a small scatterer embedded in the larger sphere: recall that a small sound-soft object scatters as a source whereas a small sound-hard object scatters as a dipole [32, section 8.2]. We plan to investigate some of these problems in the future.