Global uniqueness in a passive inverse problem of helioseismology

TL;DR

This paper proves that solar interior parameters like sound speed, density, and attenuation can be uniquely determined from acoustic observations at specific heights and frequencies, supported by numerical simulations.

Contribution

It establishes the theoretical uniqueness of recovering solar parameters from limited acoustic measurements and confirms it through numerical simulations.

Findings

Observations at two heights and frequencies uniquely determine solar parameters.

Numerical simulations support the theoretical uniqueness result.

Provides a foundation for inverse helioseismology methods.

Abstract

We consider the inverse problem of recovering the spherically symmetric sound speed, density and attenuation in the Sun from the observations of the acoustic field randomly excited by turbulent convection. We show that observations at two heights above the photosphere and at two frequencies above the acoustic cutoff frequency uniquely determine the solar parameters. We also present numerical simulations which confirm this theoretical result.

Equations

$$\begin{gathered}-\nabla\left(\frac{1}{\rho}\nabla\bigl{(}\sqrt{\rho}\psi_{\omega}\bigr{)}\right)-\frac{\sigma^{2}}{c^{2}\sqrt{\rho}}\psi_{\omega}=\frac{f_{\omega}}{\sqrt{\rho}},\\ \psi_{\omega}=\sqrt{\rho}c^{2}\nabla\cdot\xi_{\omega},\quad\sigma^{2}=\omega^{2}+2i\omega\gamma,\end{gathered}$$

$$\text{c∈L∞(R+), c≥cmin>0 a.e.},\text{ρ∈W2,∞(R+), ρ>0},\gamma \in L^{\infty }({\mathbb{R}}_{+}),$$

$$c(r)=c_0,\rho (r)={\rho }_0\exp (-(r-R_a)/H),\gamma (r)=0,r\ge R_a,$$

$$\omega >\frac{c_0}{2H},$$

$$(L_v-k^2)\psi =f,L_v=-\Delta +v,k>0,$$

$$k=\sqrt{\frac{{\omega }^2}{c_0^2}-\frac{1}{4H^2}},v(x)=k^2-\frac{\sigma (|x|{)}^2}{c(|x|{)}^2}+\rho (|x|{)}^{\frac{1}{2}}\Delta (\rho (|x|{)}^{-\frac{1}{2}}),$$

$$\begin{array}{c}v(x)=\tilde{v}(|x|),x\in {\mathbb{R}}^3,\\ \tilde{v}\in L^{\infty }({\mathbb{R}}_{+}),\tilde{v}(r)=\frac{\alpha }{r},r\ge R_a,\\ \text{for some constants α∈R, Ra>0}.\end{array}$$

$$\begin{array}{c}\text{ΣvP⊂R+ is the union of the sets of positive poles k}\\ \text{of functions Gv(x,x′;k) and Gv(x,x′;k)}.\end{array}$$

$$\Vert u{\Vert }_{L_{\delta }^2}={\left[{\int }_{{\mathbb{R}}^d}(1+|x|{)}^{\delta }|u(x){|}^{2}dx\right]}^{1/2},\delta \in \mathbb{R}.$$

$${\psi }_v(x)={\int }_{{\mathbb{R}}^3}{G}_v(x,x^{\prime })f(x^{\prime })d{x}^{\prime }.$$

$${\psi }_v(r\vartheta )=\frac{1}{r}\sum _{\ell \ge 0}\sum _{|m|\le \ell }{\phi }_{v,\ell }^m(r)Y_{\ell }^m(\vartheta ),f(r\vartheta )=\frac{1}{r}\sum _{\ell \ge 0}\sum _{|m|\le \ell }{f}_{\ell }^m(r)Y_{\ell }^m(\vartheta ),$$

$$(L_{v,\ell }-k^2){\phi }_{v,\ell }^m=f_{\ell }^m,L_{v,\ell }=-\frac{d^2}{d{r}^2}+\frac{\ell (\ell +1)}{r^2}+\tilde{v},$$

$${\phi }_{v,\ell }^m(r)={\int }_0^{R_a}{G}_{v,\ell }(r,r^{\prime })f_{\ell }^m(r^{\prime })d{r}^{\prime },$$

$$G_v(r\vartheta ,r^{\prime }\vartheta )=\frac{1}{r{r}^{\prime }}\sum _{\ell \ge 0}\sum _{|m|\le \ell }{G}_{v,\ell }(r,r^{\prime })Y_{\ell }^m(\vartheta )\overline{Y_{\ell }^m}({\vartheta }^{\prime }),r,r^{\prime }>0,\vartheta ,{\vartheta }^{\prime }\in S_1^2.$$

$$\mathbb{E}\bigl{(}\overline{f^{m}_{\ell}(r)}f^{m}_{\ell}(r+\xi)\bigr{)}=\Pi\delta_{0}(\xi)\left(-\Im\widetilde{v}(r)+k\delta_{R}(r)\right),$$

$$\Pi\Im G_{v,\ell}(r_{1},r_{2})=\mathbb{E}\bigl{(}\overline{\varphi^{m}_{v,\ell}(r_{1})}\varphi^{m}_{v,\ell}(r_{2})\bigr{)}+O(\tfrac{1}{R}),\quad R\to+\infty.$$

$$M^{1}_{R}(x_{1},x_{2})=\bigl{\{}(x_{1}\sin\theta+x_{2}\cos\theta,x_{2})\colon\theta\in(0,\pi]\bigr{\}},$$

$$B^{3}_{R}=\bigl{\{}x\in\mathbb{R}^{3}\colon|x|<R\bigr{\}}.$$

$${\psi }_{v,\ell }^m(x)=\frac{1}{|x|}{\phi }_{v,\ell }(|x|)Y_{\ell }^m(\frac{x}{|x|}),$$

$$\bigl{(}L_{v,\ell}-k^{2}\bigr{)}\varphi_{v,\ell}=0,$$

$$\underset{r\to +0}{\lim }{r}^{-\ell -1}{\phi }_{v,l}(r)=1,$$

$$(L_v-k^2)\psi =0\text{in BR3},\psi {|}_{S_R^2}=Y_{\ell }^m,$$

$$\psi (x)=\frac{R}{|x|}\frac{{\phi }_{v,\ell }(|x|)}{{\phi }_{v,\ell }(R)}{Y}_{\ell }^m(\frac{x}{|x|}).$$

$$H_{\ell }^{\pm }(\eta ,kr)=\exp (\pm i{\theta }_{\ell }(\eta ,kr))+O(\frac{1}{r}),$$

$${\theta }_{\ell }(\eta ,kr)=kr-\eta \ln (2kr)-\frac{1}{2}\ell \pi +{\sigma }_{\ell }(\eta ),$$

$$\frac{\partial }{\partial \rho }{H}_{\ell }^{\pm }(\eta ,\rho )=\left(\frac{\ell +1}{\rho }+\frac{\eta }{\ell +1}\right)H_{\ell }^{\pm }(\eta ,\rho )-\sqrt{1+\frac{{\eta }^2}{(\ell +1{)}^2}}{H}_{\ell +1}^{\pm }(\eta ,\rho ),\ell \ge 0$$

$$\frac{\partial }{\partial r}{H}_{\ell }^{\pm }(\eta ,kr)=\pm ik{H}_{\ell }^{\pm }(\eta ,kr)+O(\frac{1}{r}),r\to +\infty .$$

$${\phi }_{v,\ell }(r)=b_{\ell }\left(H_{\ell }^{-}(\eta ,kr)-s_{v,\ell }{H}_{\ell }^{+}(\eta ,kr)\right),r\ge R_a,$$

$${\phi }_{v,\ell }(r)=a_{v,\ell }{H}_{\ell }^{+}(\eta ,kr)+b_{v,\ell }{H}_{\ell }^{-}(\eta ,kr),r\ge R,$$

$${\int }_{|x|\ge 1}(1+|x|{)}^{-1+\epsilon }|\mathcal{D}\psi (x){|}^{2}dx<\infty ,\epsilon \in (0,\frac{1}{2}],$$

$$\mathcal{D}\psi (x)=\nabla \psi (x)+\frac{x}{|x{|}^2}\psi (x)-ik\frac{x}{|x|}\psi (x).$$

Full text

Global uniqueness in a passive inverse problem of helioseismology

A. D. Agaltsov Max-Planck-Institut für Sonnensystemforschung, Justus-von-Liebig-Weg 3, 37077 Göttingen, Germany (agaltsov@mps.mpg.de).   

T. Hohage Institut für Numerische und Angewandte Mathematik, Georg-August-Universität Göttingen, Lotzestr. 16-18, 37083 Göttingen, Germany (hohage@math.uni-goettingen.de) and Max-Planck-Institut für Sonnensystemforschung.   

and R. G. Novikov CMAP, Ecole Polytechnique, CNRS, Université Paris-Saclay, 91128 Palaiseau, France; and IEPT RAS, 117997 Moscow, Russia (novikov@cmap.polytechnique.fr).

We consider the inverse problem of recovering the spherically symmetric sound speed, density and attenuation in the Sun from the observations of the acoustic field randomly excited by turbulent convection. We show that observations at two heights above the photosphere and at two frequencies above the acoustic cutoff frequency uniquely determine the solar parameters. We also present numerical simulations which confirm this theoretical result.

Keywords: inverse scattering problems, uniqueness for inverse problems, helioseismology, passive imaging, long range scattering

AMS subject classification: 35R30, 85A15, 35J10, 65N21

1 Introduction

1.1 Acoustic field in the Sun and its measurement

Turbulent convection in the upper layers of the solar convection zone can reach almost sonic speeds and serves as an efficient driving mechanism for acoustic oscillations [5]. We consider the three-dimensional equation describing these oscillations at fixed frequency $\omega>0$ proposed in [11]:

[TABLE]

where $\xi_{\omega}$ is the spatial matter displacement vector, $c$ is the sound speed, $\rho$ is the density, $\gamma$ is the attenuation, $f_{\omega}$ is the random source field due to turbulent convection and $x\in\mathbb{R}^{3}$. In this work we consider this model under an additional spherical symmetry assumption: $c=c(|x|)$, $\rho=\rho(|x|)$, $\gamma=\gamma(|x|)$. We also assume that

[TABLE]

where $\mathbb{R}_{+}=(0,\infty)$ and $W^{2,\infty}(\mathbb{R}_{+})$ denotes the Sobolev space of functions defined on $\mathbb{R}_{+}$ which belong to $L^{\infty}(\mathbb{R}_{+})$ together with their first two derivatives. We suppose that in the upper atmosphere $|x|\geq R_{a}$ the sound speed is constant, the density is exponentially decreasing (which corresponds to the adiabatic approximation, see [5, section 5.4]), and there are no attenuation:

[TABLE]

where $R_{a}=R_{\odot}+h_{a}$, $R_{\odot}=6.957\times 10^{5}$ km is the solar radius, $h_{a}$ is the altitude at which the (conventional) interface between the lower and upper parts of the atmosphere is located, and $H$ is called density scale height. The first two assumptions of formula (3) follow the model of [10], which extends a standard solar model of [6] to the upper atmosphere. In this article we do not fix exact values of the above parameters, but recall that in [10] they are given by the following table:

[TABLE]

In reality, the Sun is surrounded by a corona whose base is located at about $h_{c}=2000\,$\mathrm{k}\mathrm{m}$$ above the surface and which is highly inhomogeneous. However, it is common to neglect this complication when studying acoustic waves inside of the Sun and in the lower atmosphere, see, e.g., [5, 10]. The adequacy of this simplification have been theoretically justified and numerically confirmed.

The exponential decay of density in the upper atmosphere results in trapping of acoustic waves with frequencies less than the cutoff frequency $\omega_{\text{ctf}}=c_{0}/(2H)$, which is about 5.2 mHz for the Sun, and to quantisation of their admissible frequencies. Observations of these frequencies and corresponding modes at the solar surface provide a common basis for helioseismological studies, see, e.g., [5]. In turn, acoustic waves with frequencies above the cutoff, that is, such that

[TABLE]

propagate into the upper atmosphere. One can expect that the simulated data (oscillation power spectrum) computed from equation (1) is relatively closer to observations for these higher frequency waves. The reason is that convection, granulation and supergranulation significantly contribute to the power of oscillations at lower frequencies (see [12]), but are not captured by the model. Possibility of helioseismic inversions from observations above the cutoff has not been well investigated to date. In this artile we show that observations of acoustic waves, when performed at two frequencies above the cutoff and at two different heights above the solar surface, uniquely determine the sound speed, density and attenuation in the Sun in the spherically symmetric case.

Experimental measurement of solar acoustic waves can be performed through the Doppler shifts in the absorption lines of the solar light, as it is done in the Helioseismic and Magnetic Imager (HMI) onboard the Solar Dynamics Observatory (SDO) satellite, see, e.g., the official SDO website (sdo.gsfc.nasa.gov) for more details and references. HMI observes the full solar disk in the Fe I absorption line at 6173 Å continuously from April 30, 2010. It combines six localised in wavelength photographs (filtergrams) taken in a neighborhood of this spectral line with a cadence of 45 sec to compute the map of Doppler velocities (Dopplergram). Several models show that this Dopplergram can serve as a rough estimate for the line-of-sight matter displacement velocity at about 100 km above the solar surface, which is the formation height for the HMI Dopplergram, see [9].

HMI is the successor of the Michelson Doppler Imager (MDI), which had similar design and purpose, onboard the Solar and Heliospheric Observatory (SOHO) satellite. In contrast to HMI, MDI observed the full solar disk in the Ni I absorption line at 6768 Å with a cadence of 60 sec and only for several months each year during the so-called Dynamic Runs. The formation height for MDI is about 125 km above the solar surface, see [9]. Note that observations from HMI and MDI can be found in the Joint Science Operations Center database at Stanford University (jsoc.stanford.edu).

In the present work we assume that the measurements of the solar acoustic field can be performed at two different heights above the surface. In a rough approximation, HMI and MDI Dopplergrams taken during the Dynamics Run 2010, when both instruments continuously observed the full solar disk, can be used to extract this data. As recently shown in [18], it is also possible to perform multi-height measurements by combining six raw HMI filtergrams in different ways.

The main theoretical results of this work are presented in section 2. Related proofs are given in section 3. Numerical reconstructions confirming our theoretical conclusions are given in section 4.

2 Main results

2.1 Extracting the imaginary part of the Green’s function

Under the assumptions (2), (3), (4) equation (1) at fixed $\omega$ can be rewritten as the Schrödinger equation

[TABLE]

where the indices indicating dependence on $\omega$ are suppressed,

[TABLE]

$v\in L^{\infty}(\mathbb{R}^{3})$, and $v(x)=1/(H|x|)$ for $|x|\geq R_{a}$.

In this article we consider equation (5) with a general complex-valued potential $v$ such that

[TABLE]

If the potential $v$ satisfies (7), then the resolvent $(L_{v}-k^{2})^{-1}$ is a meromorphic operator-valued function of k\in\mathbb{C}_{+}=\bigl{\{}z\in\mathbb{C}\colon\Im z>0\bigr{\}} with the distributional kernel $G_{v}(x,x^{\prime})=G_{v}(x,x^{\prime};k)$ admitting a unique meromorphic continuation across the positive real axis. The restriction to $k\in\mathbb{R}_{+}$ of the distributional kernel $G_{v}(x,x^{\prime})$ is called the radiation Green’s function for equation (5). In addition, $G_{v}(x,x^{\prime})$ is a distributional solution to equation $(L_{v}-k^{2})G_{v}(\cdot,x^{\prime})=\delta_{x^{\prime}}$, where $\delta_{x^{\prime}}$ denotes the Dirac delta function centered at $x^{\prime}$. We also suppose that $k\in\mathbb{R}_{+}\setminus\Sigma^{P}_{v}$, where

[TABLE]

Basic properties of $G_{v}$ can be found in [22, 3]. In particular, at fixed $k$ the function $G_{v}$ is jointly continuous outside the diagonal $\Delta=\{(x,x^{\prime})\in\mathbb{R}^{3}\times\mathbb{R}^{3}\colon x=x^{\prime}\}$. Besides, $(L_{v}-(k+i0)^{2})^{-1}\in\mathcal{L}(L^{2}_{1+\varepsilon}(\mathbb{R}^{3}),L^{2}_{-1-\varepsilon}(\mathbb{R}^{3}))$, $\varepsilon\in(0,\tfrac{1}{2}]$, where $L^{2}_{\delta}(\mathbb{R}^{3})$ denotes the Hilbert space of measurable functions $u$ in $\mathbb{R}^{3}$ with the finite norm

[TABLE]

Accordingly, for any $k\in\mathbb{R}_{+}\setminus\Sigma^{P}_{v}$ equation (5) with $f\in L^{2}_{1+\varepsilon}(\mathbb{R}^{3})$, $\varepsilon\in(0,\tfrac{1}{2})$, admits a unique radiation (limiting absorption), solution given by

[TABLE]

Spherical symmetry of the potential $v$ allows to separate variables in equation (5), reducing it to an equivalent multi-channel Schrödinger equation on the half-line $\mathbb{R}_{+}$ with non-coupled channels. More precisely, consider the orthogonal expansions in normalized spherical harmonics $Y^{m}_{\ell}$:

[TABLE]

where $r>0$, $\vartheta\in S^{2}_{1}$ and S^{2}_{R}=\bigl{\{}x\in\mathbb{R}^{3}\colon|x|=R\}. Plugging these expansions into formula (5), we get the radial equations

[TABLE]

where $\ell\geq 0$, $|m|\leq\ell$. Besides, it follows from formulas (9), (10) that if $\psi_{v}$ is a unique radiation solution of equation (5), then $\varphi^{m}_{\ell}$ can be expressed as:

[TABLE]

where $G_{v,\ell}(r,r^{\prime})$ is the coefficient in the spherical harmonics expansion of $G_{v}$:

[TABLE]

One can show that $G_{v,\ell}$ is indeed a Green’s function for equation (11), that is, $(L_{v,\ell}-k^{2})G_{v,\ell}(\cdot,r^{\prime})=\delta_{r^{\prime}}$, see [3] and section 3.2.

In this article we consider equation (5) with a random source function $f$. In this case the radiation solution $\psi_{v}$ is also a random function, as well as functions $\varphi^{m}_{v,\ell}$ and $f^{m}_{\ell}$ in the spherical harmonics expansions of formula (10). Following [11], we assume that the power spectrum (power spectral density) of $\psi_{v}$, defined as $\mathcal{P}^{m}_{v,\ell}(r)=\mathbb{E}|\varphi^{m}_{v,\ell}(r)|^{2}$, can be measured experimentally. However, in contrast to [11], where the power spectral density is assumed to be known at the solar surface $r=R_{\odot}$, we assume that it can be measured at two different observation radii $R_{o}^{\dagger}>R_{o}\geq R_{\odot}$. These measurements can be roughly achieved by using concurrent MDI and HMI Dopplergrams [9], or multi-height measurements from raw HMI filtergrams [18].

Our first result relates cross correlations \mathbb{E}\bigl{(}\overline{\varphi^{m}_{v,\ell}(r_{1})}\varphi^{m}_{\ell}(r_{2})\bigr{)} to the Green’s function $G_{v,\ell}(r_{1},r_{2})$. We prove the following proposition:

Proposition** 1****.**

Let $v$ be a complex-valued potential satisfying (7) and let $k\in\mathbb{R}_{+}\setminus\Sigma^{P}_{v}$ be fixed. Assume that the random functions $f^{m}_{\ell}$ satisfy the condition

[TABLE]

for some $\Pi>0$, $R\geq R_{a}$. Then the following formula is valid at fixed $r_{1}$, $r_{2}>0$:

[TABLE]

Proposition 1 is proved in section 3.3. This proposition is a variation of a well-known result, see, e.g., [23, 11] and references therein. The main difference is that we consider long range potentials and the radiation Green’s function, whereas in the literature the Green’s function with an artificial boundary condition imposed at $r=R$ and approximating the Sommerfeld radiation condition is used. The approximate radiation boundary condition allows to get rid of the error term $O(\tfrac{1}{R})$ in the formula (15) but complicates the further analysis. In addition, note that in general the Sommerfeld radiation condition does not apply for long range potentials.

Assumption (14) requires that the random sources be uncorrelated in space, excited throughout the volume with a power proportional to $-\Im\widetilde{v}$, and excited at the surface $r=R$ with a power proportional to $k$.

Remark** 1****.**

Recall that equation (5) arises, in particular, by rewriting equation (1) under the assumptions (2), (3), (4). In this case $k$ and $v$ are given by formulas (6), and Proposition 1 has a physical interpretation. Taking into account that $-\Im\widetilde{v}=2\omega\gamma$, condition (14) implies proportionality of the power spectral density of random excitations to the local attenuation (energy dissipation) rate. It has long been known in physics that this condition is related to the possibility to extract the imaginary part of the point-source response function (Green’s function) from the power spectral density of the randomly excited field, which is expressed by relation (15) in our setting. In physical literature similar relations are etablished in fluctuation-dissipation theorems, see, e.g., [16].

2.2 Uniqueness results

Proposition 1 allows to retrieve $\Im G_{v,\ell}(r,r)$ approximately from the power spectral density of noise $\mathcal{P}^{m}_{v,\ell}(r)=\mathbb{E}|\varphi^{m}_{v,\ell}(r)|^{2}$ at fixed $r$. Next, we prove that $\Im G_{v,\ell}(r,r)$ known exactly for all $\ell\geq 0$ and at two different $r$ uniquely determines $v$. Equivalently, taking into account the orthogonal expansion (13), $v$ is uniquely determined by $\Im G_{v}$ known on $S^{2}_{r}\times S^{2}_{r}$ at two different $r$.

Theorem** 1****.**

let $v_{1}$, $v_{2}$ be two complex-valued potentials satisfying (7) and let $k\in\mathbb{R}_{+}\setminus(\Sigma^{P}_{v_{1}}\cup\Sigma^{P}_{v_{2}})$ be fixed. Assume that that one of the following conditions holds true:

1. (A)

$G_{v_{1}}=G_{v_{2}}$* on $M^{4}_{R_{o}}=(S^{2}_{R_{o}}\times S^{2}_{R_{o}})\setminus\Delta$ for some $R_{o}>R_{a}$;* 2. (B)

$\Im G_{v_{1}}=\Im G_{v_{2}}$* on $M^{4}_{R_{o}}\cup M^{4}_{R_{o}^{\dagger}}$ for some $R_{o}^{\dagger}>R_{o}\geq R_{a}$ such that $R_{o}^{\dagger}\not\in\Sigma^{S}_{\alpha,k,R_{o}}$, where $\Sigma^{S}_{\alpha,k,R_{o}}\subset[R_{o},\infty)$ is a discrete set without finite accumulation points defined by (34b), (35).*

Then $v_{1}=v_{2}$ a.e.

Remark** 2****.**

If $v$ is some potential satisfying (7) then the restriction of $G_{v}(x,x^{\prime})$ to $M^{4}_{R}$ depends on $|x-x^{\prime}|$ only because $v$ is spherically symmetric. In particular, Theorem 1 remains valid if the four-dimensional manifolds $M^{4}_{R}$ are replaced by the one-dimensional manifolds

[TABLE]

for some fixed $x_{1}$, $x_{2}\in S^{2}_{R}$ with $x_{1}\cdot x_{2}=0$.

Theorem 1 is proved in section 3 and the proof consists of the following steps presented in sections 3.1, 3.2, 3.4 and 3.5. In section 3.1 we separate variables in the equation $(L_{v}-k^{2})\psi=0$ and establish auxilary results for the regular solutions of the arising radial Schrödinger equations. In section 3.2 we derive an appropriate relation between the Green’s function $G_{v}$ and the Green’s functions $G_{v,\ell}$ of the radial equations. This relation will allow to extract the diagonal values $G_{v,\ell}(R,R)$ from $G_{v}$ on $M^{4}_{R}$. Using this relation, in section 3.4 we show that the scattering matrix elements $s_{v,\ell}$ can be extracted from $G_{v}$ on $M^{4}_{R_{o}}$ or from the imaginary part $\Im G_{v}$ only on $M^{4}_{R_{o}}\cup M^{4}_{R_{o}^{\dagger}}$, under the assumption that $R_{o}^{\dagger}\not\in\Sigma^{S}_{\alpha,k,R_{o}}$. In section 3.5 we prove that the scattering matrix elements $s_{v,\ell}$ determine the Dirichlet-to-Neumann map $\Lambda_{v,R}$ for potential $v$ in some ball $B^{3}_{R}$, $R\geq R_{a}$, where

[TABLE]

In section 3.6 we combine these results together with the uniqueness theorem for the Dirichlet-to-Neumann map from [20] to uniquely determine $v$. This will prove Theorem 1.

Corollary** 1****.**

Let $c$, $\rho$, $\gamma$, and $c^{\prime}$, $\rho^{\prime}$, $\gamma^{\prime}$ be two sets of parameters satisfying (2), (3) and define the corresponding potentials $v=v_{\omega}$, $v^{\prime}=v_{\omega}^{\prime}$ and the wavenumber $k=k_{\omega}$ according to formula (6) at fixed $\omega$. Let $\omega_{1}\neq\omega_{2}$ be two positive frequencies satisfying (4) and such that $k_{\omega}\in\mathbb{R}_{+}\setminus(\Sigma^{P}_{v_{\omega}}\cup\Sigma^{P}_{v_{\omega}^{\prime}})$ for $\omega=\omega_{1}$, $\omega_{2}$. Let $G_{v_{\omega}}$, $G_{v_{\omega}^{\prime}}$ be the radiation Green’s functions at fixed $\omega$ for the potentials $v$, $v^{\prime}$ respectively.

Suppose that $\Im G_{v_{\omega}}=\Im G_{v_{\omega}^{\prime}}$ on $M^{4}_{R_{0}}\cup M^{4}_{R_{o}^{\dagger}}$ for some $R_{o}^{\dagger}>R_{o}\geq R_{a}$ such that $R_{o}^{\dagger}\not\in\Sigma^{S}_{1/H,k_{\omega},R_{o}}$, where $\omega=\omega_{1}$, $\omega_{2}$. Then $c=c^{\prime}$, $\rho=\rho^{\prime}$, $\gamma=\gamma^{\prime}$ a.e.

Proof.

Under the assumptions of Corollary 1 it follows from Theorem 1 that $v_{\omega}=v_{\omega}^{\prime}$ for $\omega=\omega_{1}$, $\omega_{2}$. Using that $\Re v_{\omega}=\Re v_{\omega}^{\prime}$ for $\omega=\omega_{1}$, $\omega_{2}$ one can show that $c=c^{\prime}$, $\rho=\rho^{\prime}$, see, e.g., the proof of [2, Theorem 2.9]. Then, recalling that $\Im v_{\omega}=-2\omega\gamma/c^{2}$, $\Im v_{\omega}^{\prime}=-2\omega\gamma^{\prime}/(c^{\prime})^{2}$, it follows from the equality $\Im v_{\omega}=\Im v_{\omega}^{\prime}$ for $\omega=\omega_{1}$, $\omega_{2}$ together with the equality $c=c^{\prime}$ that $\gamma=\gamma^{\prime}$, concluding the proof of Corollary 1. ∎

In section 4 we shall present numerical simulations which confirm the uniqueness results of Theorem 1 and Corollary 1.

3 Proof of the main results

3.1 Properties of regular radial solutions

In this subsection we shall establish some auxilary results regarding regular solutions of the radial Schrödinger equation which arises by separation of variables in the homogeneous equation $(L_{v}-k^{2})\psi=0$.

As the potential $v$ is spherically symmetric, this equation separates in spherical coordinates. We seek a solution $\psi^{m}_{v,\ell}\in H^{2}_{\text{loc}}(\mathbb{R}^{3})$ of the form

[TABLE]

which leads to the equation

[TABLE]

together with the condition that $\varphi_{v,\ell}$ vanishes at the origin. One can show that this determines $\varphi_{v,\ell}\in H^{2}_{\text{loc}}(\mathbb{R}_{+})$ uniquely up to a multiplicative factor, see, e.g., [3]. We impose the boundary condition

[TABLE]

which fixes $\varphi_{v,\ell}$ uniquely.

Note that $\varphi_{v,\ell}(r)$ does not depend on the values of $\widetilde{v}$ in the region $[r,+\infty)$, see [19, formula (12.4)]. This analysis implies the following lemma.

Lemma** 1****.**

Let $k>0$ and let $v$ be a complex-valued potential satisfying (7). Then the Dirichlet problem

[TABLE]

where $Y^{m}_{\ell}=Y^{m}_{\ell}(x/|x|)$, $x\in S^{2}_{R}$, has a unique solution $\psi\in H^{2}(B^{3}_{R})$ if and only if $\varphi_{v,\lambda}(R)\neq 0$ for all integer $\lambda\geq 0$. In addition, this solution is given by the formula

[TABLE]

Next we shall derive an expression for the regular solution $\varphi_{v,\ell}$ in the domain $r\geq R$ in terms of the Coulomb wave functions and of the so-called scattering matrix element $s_{v,\ell}$. First we recall the definition and some basic properties of the Coulomb wave functions from [7, 1].

The Coulomb wave functions $H^{\pm}_{\ell}(\eta,kr)$, $\eta=\alpha/(2k)$, are the unique solutions of equation (17) with $\widetilde{v}(r)=\alpha/r$ specified by the following asymptotics as $r\to+\infty$:

[TABLE]

where $\sigma_{\ell}(\eta)=\arg\Gamma(\ell+1+i\eta)$ is the Coulomb phase shift and $\Gamma$ denotes the usual gamma function. Functions $H^{+}_{\ell}(\eta,kr)$ and $H^{-}_{\ell}(\eta,kr)$ are complex conjugates of each other and are linearly independent. Using (20a) together with the relation

[TABLE]

given in [21], one can show that

[TABLE]

Using these properties, we shall prove the following result.

Lemma** 2****.**

Let $v$ be a complex-valued potential satisfying (7), let $k\in\mathbb{R}_{+}\setminus\Sigma^{P}_{v}$ be fixed and let $\eta=\alpha/(2k)$. Then the function $\varphi_{v,\ell}$ defined by (17), (18) admits in the region $r\geq R_{a}$ the representation

[TABLE]

with unique $b_{\ell}=b_{\ell}(k)\in\mathbb{C}\setminus\{0\}$ and $s_{v,\ell}=s_{v,\ell}(k)\in\mathbb{C}\setminus\{0\}$.

Proof.

As the Coulomb wave functions $H^{\pm}_{\ell}(\eta,kr)$ are linearly independent, any solution to equation (17) in the region $r\geq R_{a}$ is given by their linear combination with unique coefficients. In particular, $\varphi_{v,\ell}$ can be expressed in the region $r\geq R$ in the form

[TABLE]

for some $a_{v,\ell}$, $b_{v,\ell}\in\mathbb{C}$. We shall show that $a_{v,\ell}\neq 0$, $b_{v,\ell}\neq 0$.

Recall from [22] that for any $k\in\mathbb{R}_{+}$ outside the singular set $\Sigma^{P}_{v}$ and for any $f\in L^{2}_{1+\varepsilon}(\mathbb{R}^{3})$ the Schrödinger equation (5) admits the unique solution $\psi\in H^{2}_{\text{loc}}(\mathbb{R}^{3})\cap L^{2}_{-1-\varepsilon}(\mathbb{R}^{3})$, $\varepsilon\in(0,\tfrac{1}{2}]$, satisfying the radiation condition

[TABLE]

Now assume that $b_{v,\ell}=0$. Then it follows from (20a), (21) that the function defined by (19) is a non-zero solution of class $H^{2}_{\text{loc}}(\mathbb{R}^{3})\cap L^{2}_{-1-\varepsilon}(\mathbb{R}^{3})$, $\varepsilon\in(0,\tfrac{1}{2}]$, to $(L_{v}-k^{2})\psi=0$ in $\mathbb{R}^{3}$ satisfying the radiation condition (23a). This contradicts the assumption $k\not\in\Sigma^{P}_{v}$.

Now assume that $a_{v,\ell}=0$ and let $\psi$ be defined by (19). Then it follows from (20a), (21) that $\overline{\psi}$ is of class $H^{2}_{\text{loc}}(\mathbb{R}^{3})\cap L^{2}_{-1-\varepsilon}(\mathbb{R}^{3})$, $\varepsilon\in(0,\tfrac{1}{2}]$, and satisfies $(L_{\overline{v}}-k^{2})\overline{\psi}=0$ in $\mathbb{R}^{3}$ together with the radiation condition (23a). Taking into account definition (8), this also contradicts the assumption $k\not\in\Sigma^{P}_{v}$ and concludes the proof of Lemma 2. ∎

Remark** 3****.**

The coefficient $s_{v,\ell}=s_{v,\ell}(k)$ in Lemma 2 is called the $\ell$-th scattering matrix element of the potential $v$.

3.2 Green’s functions

In this subsection we shall express the radiation Green’s function $G_{v,\ell}$ for equation (17) in the region $r\geq R_{a}$ in terms of the Coulomb wave functions. Then we shall give a formula for extracting the diagonal values $G_{v,\ell}(R,R)$ from the radiation Green’s functon $G_{v}$ for equation (5) known at $M^{4}_{R}$.

In addition to the regular solution $\varphi_{v,\ell}$ of equation (5) specified by the boundary condition (18), we consider the outgoing solution $\varphi^{+}_{v,\ell}$ which is specified by the asymptotics

[TABLE]

The outgoing Green’s function $G_{v,\ell}(r,r^{\prime})$ for equation (17) is defined as a distributional solution to the equation $(L_{v,\ell}-k^{2})G_{v,\ell}(\cdot,r^{\prime})=\delta_{r^{\prime}}$ specified by the following boundary conditions at fixed $r^{\prime}>0$:

[TABLE]

for some non-zero constant $c=c(r^{\prime},\eta,k)$. If the regular solution $\varphi_{v,\ell}(r)$ and the outgoing solution $\varphi^{+}_{v,\ell}(r)$ are linearly independent, this Green’s function exists, is unique and is given by the explicit formula

[TABLE]

where $r_{<}=\min(r,r^{\prime})$, $r_{>}=\max(r,r^{\prime})$ and the Wronskian

[TABLE]

is independent of $r$. Note that formula (25) is a standard result from the theory of Sturm-Liouville problems, see, e.g., [25, p. 158, formula (5.65)].

Lemma** 3****.**

If $k\in\mathbb{R}_{+}\setminus\Sigma^{P}_{v}$, then $G_{v,\ell}$ is well-defined and is given in the region $r\geq R_{a}$, $r^{\prime}\geq R_{a}$ by the formula

[TABLE]

Proof.

It follows from Lemma 2 that under the assumption $k\in\mathbb{R}_{+}\setminus\Sigma^{P}_{v}$ the functions $\varphi_{v,\ell}(r)$ and $H^{+}_{\ell}(\eta,kr)$ are linearly independent in the region $r\geq R_{a}$. Using the relation $[H^{-}_{\ell}(\eta,kr),H^{+}_{\ell}(\eta,kr)]=2ik$, given in [1], and formula (22) we get $[\varphi_{v,\ell},\varphi^{+}_{v,\ell}]=2ikb_{v,\ell}$. Together with (25), this implies (27). ∎

Next we shall show how the diagonal values $G_{v,\ell}(R,R)$ can be extracted from the Green’s function $G_{v}$ restricted to $M^{4}_{R}$.

First recall that the Legendre polynomials $P_{\ell}$ can be defined using the formal generating identity [24]

[TABLE]

Using the Laplace formula [24, Theorem 8.21.2] and the Dirichlet convergence test one can show that at fixed $t=1$ this series converges pointwise for all $s\in(-1,1)$. Besides, the Legendre polynomials form a complete orthogonal system in $L^{2}(-1,1)$ such that

[TABLE]

where $\delta_{\ell m}$ is the Kronecker delta.

Now we recall [3] that at fixed $k\in\mathbb{R}_{+}\setminus\Sigma^{P}_{v}$ the Green’s function $G_{v}(x,x^{\prime})$ is continuous outside of the diagonal $x=x^{\prime}$ and $G_{v}(x,x^{\prime})=O(|x-x^{\prime}|^{-1})$ as $x\to x^{\prime}$. Besides, the series expansion

[TABLE]

converges for $x\neq\pm x^{\prime}$ and $G_{v,\ell}(R,R)$ has the asymptotics

[TABLE]

Note that series expansion (30) follows from the spherical harmonics expansion (13), taking into account the well known addition theorem (see, e.g., [3]):

[TABLE]

Lemma** 4****.**

Let $k\in\mathbb{R}_{+}\setminus\Sigma^{P}_{v}$ and fix $x^{\prime}$, $x^{\prime\prime}\in S^{2}_{R}$ such that $x^{\prime}\cdot x^{\prime\prime}=0$. Then for each $\ell\geq 0$

[TABLE]

Proof.

Using (28) and (30) we get a pointwise convergent series expansion

[TABLE]

which, in view of (29) and (31), also converges in $L^{2}(-1,1)$. Recalling that Legendre polynomials form a complete orthogonal system in $L^{2}(-1,1)$, we get Lemma 4. ∎

3.3 Extracting the Green’s function from cross correlations

In this subsection we prove Proposition 1. First, recall that the Green’s function $G_{v,\ell}$ satisfies the reciprocity relation

[TABLE]

To prove it, consider the equations

[TABLE]

Multiplying the first equation by $G_{v,\ell}(\cdot,r_{2})$, subtracting the second equation multiplied by $G_{v,\ell}(\cdot,r_{1})$, and integrating over $(0,R)$, $R>r_{1}$, $r_{2}$, we obtain

[TABLE]

where $[-,-]$ denotes the Wronskian defined according to (26) and the notation $x|^{b}_{a}=x(b)-x(a)$ is used. The next step is to show that the term on the left hand side vanishes. Using formulas (18), (25) and the estimate $\tfrac{\partial}{\partial r}\varphi_{v,\ell}(r)=O(r^{\ell})$, $r\to+0$, given in [3, Theorem 3.3], we get $[G_{v,\ell}(\cdot,r_{1}),G_{v,\ell}(\cdot,r_{2})](+0)=0$. Using formulas (20a), (20b), (21), (25), we also get $[G_{v,\ell}(\cdot,r_{1}),G_{v,\ell}(\cdot,r_{2})](R)=O(\tfrac{1}{R})$, $R\to+\infty$. As $R$ tends to $+\infty$, we get formula (32).

To prove (15), we follow a similar scheme. We start from the equations

[TABLE]

Multiplying the first equation by $G_{v,\ell}(\cdot,r_{2})$, subtracting the second equation multiplied by $\overline{G_{v,\ell}(\cdot,r_{1})}$, integrating over $(0,R)$, $R>r_{1}$, $r_{2}$, we get

[TABLE]

In a similar way with the proof of formula (32), one can show that the Wronskian vanishes at zero and that

[TABLE]

Combining this with formulas (12), (14), (33), (32) to compute \mathbb{E}\bigl{(}\overline{\psi^{m}_{\ell}(r_{1})}\psi^{m}_{\ell}(r_{2})\bigr{)}, we get formula (15), which concludes the proof of Proposition 1.

3.4 Recovering the scattering matrix elements

In this subsection we shall show that the scattering matrix elements $s_{v,\ell}$ for the potential $v$ can be extracted from the Green’s function $G_{v}$ on $M^{4}_{R_{o}}$ or from its imaginary part $\Im G_{v}$ only on $M^{4}_{R_{o}}\cup M^{4}_{R_{o}^{\dagger}}$, where $R_{a}\leq R_{o}<R_{o}^{\dagger}$.

Recall that the Coulomb function $H^{+}_{\ell}(\eta,kr)$ does not vanish for $r>0$, since $H^{+}_{\ell}(\eta,kr)$ and its complex conjugate $H^{-}_{\ell}(\eta,kr)$ form a basis of solutions of equation (17) with $\widetilde{v}(r)=\alpha/r$. Together with Lemmas 3 and 4, this leads to the following result.

Lemma** 5****.**

Let $v_{1}$, $v_{2}$ be two complex-valued potentials satisfying (7), let $k\in\mathbb{R}_{+}\setminus(\Sigma^{P}_{v_{1}}\cup\Sigma^{P}_{v_{2}})$ be fixed and let $R_{o}\geq R_{a}$. Suppose that $G_{v_{1}}=G_{v_{2}}$ on $M^{4}_{R_{o}}$. Then $G_{v_{1},\ell}(R_{o},R_{o})=G_{v_{2},\ell}(R_{o},R_{o})$ and $s_{v_{1},\ell}=s_{v_{2},\ell}$ for all $\ell\geq 0$.

Next we shall show that the scattering matrix elements $s_{v,\ell}$ can be extracted from $\Im G_{v}$ only, i.e., without knowing $\Re G_{v}$. However, the values of $\Im G_{v}$ must be given not only on $M^{4}_{R_{o}}$ but also on $M^{4}_{R_{o}^{\dagger}}$ for some $R_{o}^{\dagger}>R_{o}$.

Note that formula (27) implies that

[TABLE]

where $r\geq R_{a}$, $\eta=\alpha/(2k)$. Using (34a) with $r=R_{o}$, $R_{o}^{\dagger}$ one can see that $s_{v,l}$ is uniquely determined from $\Im G_{v,\ell}(R_{o},R_{o})$ and $\Im G_{v,\ell}(R_{o}^{\dagger},R_{o}^{\dagger})$ if and only if $\sin(\vartheta_{\ell}(\eta,kR_{o}^{\dagger})-\vartheta_{\ell}(\eta,kR_{o}))\neq 0$. This justifies the definition of the singular set

[TABLE]

Lemma** 6****.**

Let $v_{1}$, $v_{2}$ be two complex-valued potentials satisfying (7), let $k\in\mathbb{R}_{+}\setminus(\Sigma^{P}_{v_{1}}\cup\Sigma^{P}_{v_{2}})$, and let $R_{o}^{\dagger}>R_{o}\geq R_{a}$ be such that $R_{o}^{\dagger}\not\in\Sigma^{S}_{\alpha,k,R_{o}}$. Suppose that $\Im G_{v_{1}}=\Im G_{v_{2}}$ on $M^{4}_{R_{o}}\cup M^{4}_{R_{o}^{\dagger}}$. Then $s_{v_{1},\ell}=s_{v_{2},\ell}$ for all $\ell\geq 0$. Besides, the set $\Sigma^{S}_{\alpha,k,R_{o}}$ is discrete and does not have finite accumulation points.

Proof.

Under the assumptions of Lemma 6 it follows from Lemma 4 that for all $\ell\geq 0$ the equality $\Im G_{v_{1},\ell}(r,r)=\Im G_{v_{2},\ell}(r,r)$ holds true for $r=R_{o}$, $R_{o}^{\dagger}$. Together with the discussion before Lemma 6 it implies that $s_{v_{1},\ell}=s_{v_{2},\ell}$ for all $\ell\geq 0$. This concludes the proof of the first assertion of Lemma 6.

Next we shall prove the second assertion. Using formulas (33.2.11), (33.5.8), (33.5.9) of [1] one can see that

[TABLE]

locally uniformly in $r\in\mathbb{R}_{+}$. Besides, it follows from formula (33.2.11) of [1] and from discussion below it that at fixed $\ell\geq 0$ the set of solutions $r\in\mathbb{R}_{+}$ of the equation $\sin(\vartheta_{\ell}(\eta,kr)-\vartheta_{\ell}(\eta,kR_{o}))=0$ is discrete and does not have finite accumulation points, as the zero-set of a non-zero analytic function. Together with (36), this concludes the proof of Lemma 6. ∎

3.5 Recovering the Dirichlet-to-Neumann map

In this subsection we shall show that the scattering matrix elements $s_{v,\ell}$ uniquely determine the Dirichlet-to-Neumann map for $v$ is some ball $B^{3}_{R}$, $R\geq R_{a}$.

Note that Lemma 1 justifies the definition of the singular set

[TABLE]

Lemma** 7****.**

Let $v$ be a complex-valued potential satisfying (7) and let $k\in\mathbb{R}_{+}$ be fixed. Then $R\in\mathbb{R}_{+}\setminus\Sigma^{D}_{v,k}$ if and only if for any $g\in H^{3/2}(S^{2}_{R})$ the Dirichlet problem

[TABLE]

is uniquely solvable for $\psi\in H^{2}(B^{3}_{R})$. Besides, the set $\Sigma^{D}_{v,k}$ is discrete and does not have finite accumulation points.

Proof.

It follows from Lemma 1 that if the Dirichlet problem (37) is uniquely solvable for any $g\in H^{3/2}(S^{2}_{R})$ then $R\not\in\Sigma^{D}_{v,k}$.

Now suppose that $R\not\in\Sigma^{D}_{v,k}$. First we shall show that (37) does not admit a non-zero solution for $g=0$. Assuming that $\psi\in H^{2}_{0}(B^{3}_{R})$ is a solution to (37) with $g=0$ one can see (taking into account the spherical symmetry of $v$) that its partial wave components

[TABLE]

belong to $H^{1}_{0}(B^{3}_{R})$ and satisfy (37) weakly. Because of the boundary elliptic regularity [8] they also belong to $H^{2}(B^{3}_{R})$. Then it follows from Lemma 1 that all $\psi^{m}_{\ell}$ vanish and $\psi=0$.

Next we recall that the operator

[TABLE]

is Fredholm of index zero from $H^{2}(B^{3}_{R})$ to $L^{2}(B^{3}_{R})\times H^{3/2}(S^{2}_{R})$, see [13]. Together with already established uniqueness for the Dirichlet problem (37) for $R\not\in\Sigma^{D}_{v,k}$, this proves the first assertion of Lemma 7.

Next we shall show that $\Sigma^{D}_{v,k}$ is a discrete set without finite accumulation points. It can be shown [3] that the regular solution $\varphi_{v,\ell}(r)$ defined by (17), (18) satisfies the estimate

[TABLE]

uniformly in $r\in(0,R]$ at fixed $R>0$. Besides, zeros of each $\varphi_{v,\ell}$ are discrete and do not have finite accumulation points. This concludes the proof of the second assertion of Lemma 7. ∎

Remark** 4****.**

Recalling from [13] that the operator (38) is Fredholm of index zero from $H^{2}(B^{3}_{R})$ to $L^{2}(B^{3}_{R})\times H^{3/2}(S^{3}_{R})$ one can see that if the potential $v\in L^{\infty}(B^{3}_{R})$ is such that the Dirichlet problem (37) with $g\in H^{3/2}(S^{2}_{R})$ is uniquely solvable for $\psi\in H^{2}(B^{3}_{R})$, then

[TABLE]

for some constant $C_{v,k,R}>0$. In addition, the trace theorem [13] leads to the estimate

[TABLE]

for some constant $C^{\prime}_{v,k,R}>0$, where $\tfrac{\partial\psi}{\partial r}=\tfrac{x}{|x|}\nabla\psi$ is the derivative of $\psi$ in the radial direction.

Under the assumption that the Dirichlet problem (37) is uniquely solvable for $\psi\in H^{2}(B^{3}_{R})$ for all $g\in H^{3/2}(S^{2}_{R})$, we define the Dirichlet-to-Neumann map \Lambda_{v,R}\in\mathcal{L}\bigl{(}H^{3/2}(S^{2}_{R}),H^{1/2}(S^{2}_{R})\bigr{)} by $\Lambda_{v,R}\varphi=\tfrac{\partial\psi}{\partial r}|_{S^{2}_{R}}$. Next we shall show that the partial scattering matrix elements $s_{v,\ell}$ known for all $\ell\geq 0$ uniquely determine the Dirichlet-to-Neumann map $\Lambda_{v,R}$.

Lemma** 8****.**

Let $v_{1}$, $v_{2}$ be two complex-valued potentials satisfying (7) and let $k\in\mathbb{R}_{+}\setminus(\Sigma^{P}_{v_{1}}\cup\Sigma^{P}_{v_{2}})$ be fixed. Besides, let $R\in[R_{a},+\infty)\setminus(\Sigma^{D}_{v_{1},k}\cup\Sigma^{D}_{v_{2},k})$. Suppose that $s_{v_{1},\ell}=s_{v_{2},\ell}$ for all $\ell\geq 0$. Then $\Lambda_{v_{1},R}=\Lambda_{v_{2},R}$.

Proof.

It follows from Lemmas 1 and 2 and from continuity of $\varphi_{v,\ell}(r)$ at $r=R$ that $\Lambda_{v_{1},R}|_{\mathcal{H}_{\ell,R}}=\Lambda_{v_{2},R}|_{\mathcal{H}_{\ell,R}}$, where $\mathcal{H}_{\ell,R}$ denotes the space of restrictions to $S^{2}_{R}$ of harmonic polynomials of degree $\ell$, spanned by the spherical harmonics $Y^{m}_{\ell}=Y^{m}_{\ell}(\tfrac{x}{|x|})$, $|m|\leq\ell$. More precisely, the following explicit formula for $\Lambda_{v_{j},R}|_{\mathcal{H}_{\ell,R}}$ is valid:

[TABLE]

where $\eta=\alpha/(2k)$ and the denominator is non-zero as $R\not\in\Sigma^{D}_{v_{j},k}$.

Now let $g\in H^{3/2}(S^{2}_{R})$ and denote by $\psi_{v_{j}}\in H^{2}(B^{3}_{R})$ the unique solution of the Dirichlet problem (37) with $v=v_{j}$. Besides, define $g_{N}$ by

[TABLE]

so that $g_{N}\to g$ in $H^{3/2}(S^{2}_{R})$ and, according to Remark 4, $\Lambda_{v_{j},R}g_{N}\to\Lambda_{v_{j},R}g$ in $H^{1/2}(B^{3}_{R})$ as $N\to\infty$. Together with (39), which shows that $\Lambda_{v_{1},R}g_{N}=\Lambda_{v_{2},R}g_{N}$, this implies that $\Lambda_{v_{1},R}g=\Lambda_{v_{2},R}g$ and concludes the proof of Lemma 8. ∎

3.6 Demonstration of the uniqueness theorem

Now we combine the preliminary results established in sections 3.1, 3.2, 3.4 and 3.5 to prove Theorem 1.

Under the assumptions of Theorem 1 it follows from Lemma 5 in case (A) and from Lemma 6 in case (B) that $s_{v_{1},\ell}=s_{v_{2},\ell}$ for all $\ell\geq 0$. Using Lemma 8 we conclude that $\Lambda_{v_{1},R}=\Lambda_{v_{2},R}$ for any $R$ in the non-empty set $[R_{a},\infty)\setminus(\Sigma^{D}_{v_{1},k}\cup\Sigma^{D}_{v_{2},k})$. It follows from the uniqueness theorem of [20], where the proof does not use that the potential is real-valued, that $v_{1}=v_{2}$ a.e. This proves Theorem 1.

4 Reconstruction

4.1 Reconstruction scheme for exact simulated data

The possibility to use measurements of the solar acoustic field at two heigths above the surface to recover the sound speed, density and attenuation inside of the Sun is confirmed by our numerical simulations. In this subsection we shall briefly describe the reconstruction algorithm that we use.

We assume that the unknown solar parameters $q=(c,\rho,\gamma)$ are perturbations of some known background quantities $q^{0}=(c^{0},\rho^{0},\gamma^{0})$ such that

[TABLE]

where $R_{\odot}=6.957\times 10^{5}$ km is the solar radius, and we assume that both parameter sets $q$ and $q^{0}$ satisfy (2), (3). Let $\Omega\subset\mathbb{R}_{+}$ be a finite set of admissible frequencies such that $\Omega\cap(\Sigma^{P}_{q}\cup\Sigma^{P}_{q^{0}})=\varnothing$, where

[TABLE]

$\Sigma^{P}_{v_{\omega}}$ is defined according to (8), and the potentials $v_{\omega}$ and $v^{0}_{\omega}$ are defined using formula (6) with parameters $q$ and $q_{0}$, respectively. We recall that $c_{0}/(2H)$ is the acoustic cutoff frequency, which separates the regime of oscillations at eigenfrequencies and the scattering regime.

Put $G_{q}(h,\ell,\omega)=G_{v_{\omega},\ell}(R_{\odot}+h,R_{\odot}+h)$. As initial data for inversions from exact data we use the imaginary part of the Green’s function $\Im G_{q}(h,\ell,\omega)$ measured at two different non-negative altitudes $h\in\{h_{1},h_{2}\}$, at angular degrees $\ell\in\{0,\dots,\ell_{\text{max}}\}$, and at all admissible frequencies $\omega\in\Omega$. From this data we recover the solar parameters $q$ as follows.

The first step of the algorithm consists in recovering the scattering matrix elements $s_{q}(\ell,\omega)=s_{v_{\omega},\ell}$, $\ell\in\{0,\dots,\ell_{\text{max}}\}$, $\omega\in\Omega$, which are defined according to Lemma 2. This reconstruction is done by considering equation (34a) with $r=R_{\odot}+h$, $h\in\{h_{1},h_{2}\}$ as a linear system for finding $\Re s_{q}(\ell,\omega)$, $\Im s_{q}(\ell,\omega)$ at each fixed $\ell$, $\omega$.

At the next step the scattering matrix elements $s_{q}(\ell,\omega)$ are used to recover the map $u_{q}\colon I\to\mathbb{R}^{3}$ (where $I$ is the interval defined in (40)) which is defined as follows:

[TABLE]

and such that

[TABLE]

This reconstruction is done by applying the iteratively regularized Gauss-Newton method, going back to [4], to the forward map

[TABLE]

For more details on the iteratively regularized Gauss-Newton method and for sufficient conditions of its convergence see, e.g., [15].

The last step is to determine $q=(c,\rho,\gamma)$ from $u_{q}$. Note that definitions (41) lead to the following explicit formulas for determining $c$ and $\gamma$:

[TABLE]

Also note that definitions (41) lead to the following problem for determining $\rho$:

[TABLE]

which is solved for the unknown function $\rho^{-\frac{1}{2}}$. This step concludes the reconstruction algorithm from exact data.

4.2 Numerical example with exact data

We consider the background sound speed and density from the model of [10], which extends the standard solar model of [6] to the upper atmosphere. We suppose that the background attenuation is equal to $\gamma_{0}=102.5$ $\mu$Hz inside of the Sun and decays to zero smoothly in the region $[R_{\odot},R_{\odot}+h_{a}]$, where $R_{\odot}=6.957\times 10^{5}$ km is the solar radius and $h_{a}=500$ km is the height above the photosphere at which the (conventional) interface between the lower and upper atmosphere is located. Note that this approximate value for the background attenuation can be obtained by analysing the observed full width at half maximum (FWHM) of acoustic modes, see [11, section 7.3] for more details111Our simulations show that another reasonable choice for the attenuation constant in the Sun is $\gamma_{0}=39.8$ $\mu$Hz. For this attenuation coefficient the difference between the observed by HMI and the modeled power spectra (after a linear transformation) is minimal.. We also assume that the unknown perturbations to the background values of solar parameters are supported in the interval $[0.9R_{\odot},0.95R_{\odot}]$.

The initial data for reconstructions is the imaginary part of the radiation Green’s function $\Im G_{q}(h,\ell,\omega)$ at heights $h\in\{105,144\}$ (km), angular degrees $\ell\in\{0,\dots,250\}$, and frequencies $\omega\in\{5.3,5.4\}$ (mHz). Note that observations of the acoustic field at these heights approximately correspond to measuring Doppler velocities of the line center (three-point approximation) and center of gravity (six-point approximation) of the Fe I absorption line at 6173 Å using the data from the Helioseismic and Magnetic Imager (HMI) as described in [18].

Figure 1 shows exact profiles of perturbations $\delta c^{\dagger}$, $\delta\rho^{\dagger}$, $\gamma^{\dagger}$, reconstructed profiles of perturbations $\delta c$, $\delta\rho$, $\delta\gamma$, and relative $L^{2}$ reconstruction errors $e(\delta c,\delta c^{\dagger})$, $e(\delta\rho,\delta\rho^{\dagger})$, $e(\delta\gamma,\delta\gamma^{\dagger})$. Here,

[TABLE]

This reconstruction example confirms the uniqueness results of section 2.2.

4.3 Reconstruction scheme for noisy data

The power spectrum $P^{m}_{v_{\omega},\ell}(r)=\mathbb{E}|\varphi^{m}_{v_{\omega},\ell}(r)|^{2}$ introduced in section 2.1 can not be measured precisely in a real experiment. A standard approach to compute an approximation to the power spectrum is to parse the time series of acoustic oscillations (after an aproppriate preprocessing) into $N$ segments of equal duration $T$, compute for each segment the sample spectrum (periodogram), and then take the arithmetic average of periodograms $\widehat{P}^{m}_{v_{\omega},\ell}(r)$. For the definition and properties of the sample power spectrum see, e.g., [14, section 6.3.1], and for the details on computation of the sample power spectrum from the observed time series of sollar oscillations see [17].

Parameter $T$ is chosen to achieve a desired frequency resolution. We recall that for a time series segment of duration $T$ the frequency resolution of the sample power spectrum is equal to $\Delta\omega=1/T$. Besides, if the time resolution (cadence) is equal to $\Delta t$ then the sample power spectrum can be computed for frequencies up to $1/(2\Delta t)$.

It is a standard assumption going back to [26] that

[TABLE]

where $\chi^{2}(2N)$ denotes the chi-squared distribution with $2N$ degrees of freedom and $r$, $\omega$, $\ell$, $m$ are fixed. We use relation (43) to simulate data for inversions with noisy data. This simulated data is then used to extract the diagonal values of the imaginary part of the radial Green’s function $\Im G_{q}(h,\ell,\omega)=\Im G_{v_{\omega},\ell}(R_{\odot}+h,R_{\odot}+h)$ using formula (15) with the $O(\tfrac{1}{R})$ term dropped and assuming that $\Pi=1$. Note that in reality $\Pi$ is a function of frequency which is not directly accessible to measurements; for a possible model of this function see, e.g., [11, section 7.5].

Then the reconstruction proceeds as described in section 4.1.

4.4 Numerical examples with noisy data

We consider a model situation where the solar oscillations are observed for a total period of eight years with the time resolution of 45 s. In this connection, we recall that the HMI observes the sollar oscillations continuously with this cadence since April 30, 2010. We also assume that the sample power spectra are computed from three-day intervals, so that the frequency resolution is equal to 3.86 $\mu$Hz and the maximal frequency is approximately 8.33 mHz.

We simulate $\widehat{P}_{v_{\omega},\ell}^{m}(R_{\odot}+h)$ according to (43) with $N=974$ (which is the number of three-day intervals constituting eight years of observations) for angular degrees $\ell\in\{0,\dots,250\}$, azimuthal degree $m=0$, observation heights $h\in\{105,144\}$ (km) and frequencies $\omega\in\{5.27,5.29,5.31,5.34,5.36,5.38\}$ (mHz). We use the same background model and the same assumptions on the unknown parameters as in section 4.2.

In contrast to reconstructions with exact data, simultaneous recovery of all the parameters from noisy data fails when these realistic settings are used. The point is that the (numerically computed) singular values of the forward map of formula (42) decrease exponentially fast, which leads to a severe ill-posedness of the inverse problem.

However, if two out of three parameters are known a priori, the third parameter is recovered with reasonable precision. Figure 2 shows parameters $\delta c^{\dagger}$, $\delta\rho^{\dagger}$, $\delta\gamma^{\dagger}$, their reconstructions $\delta c$, $\delta\rho$, $\delta\gamma$ from noisy data and relative $L^{2}$ reconstruction errors for one realization of data. The mean relative $L^{2}$ reconstruction errors for parameters $\delta c^{\dagger}$, $\delta\rho^{\dagger}$, $\delta\gamma^{\dagger}$ are equal to

[TABLE]

Besides, standard deviations of relative $L^{2}$ reconstruction errors are equal to 3%, 12.6%, 3.1%. We emphasize that in each of these examples two out of three parameters are known a priori and fixed, and we reconstruct the remaining parameter.

These simulations show that reconstructions from noisy simulated data and, as a corollary, from experimental data, require a separate and detailed treatment, which is out of scope of the present article.

5 Conclusion

We considered the inverse problem of recovering the radially symmetric sound speed, density and attenuation in the Sun from the measurements of the solar acoustic field at two heights above the photosphere and for a finite number of frequencies above the acoustic cutoff frequency. We showed that this problem reduces to recovering a long range potential (with a Coulomb-type decay at infinity) in a Schrödinger equation from the measurements of the imaginary part of the radiation Green’s function at two distances from zero. We demonstrated that generically this inverse problem for the Schrödinger equation admits a unique solution, and that the original inverse problem for the Sun admits a unique solution when measurements are performed at least two different frequencies above the cutoff frequency. These uniqueness results are confirmed by numerical experiments with simulated data without noise. However, simulations also show that the inverse problem is severly ill-posed, and only separate recovery of one of the solar parameters (i.e. when two other parameters are fixed) using a standard iterative reconstruction method (IRGNM) is reasonably precise for realistic noise levels.