Global uniqueness in a passive inverse problem of helioseismology

Solar oscillations are excited near the solar surface by turbulent convection. We consider the approximate equations describing these oscillations in three-dimensions at fixed angular frequency ω > 0 proposed in [1]:

$$\begin{equation}\begin{array}{c}-\nabla \left(\frac{1}{\rho }\nabla \left(\sqrt{\rho }{\psi }_{\omega }\right)\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 \quad {\sigma }^{2}={\omega }^{2}+2\text{i}\omega \gamma ,\end{array}\end{equation} \tag{ 1 }$$

where ξ_ω is the spatial displacement vector, c is the sound speed, ρ is the density, γ is the attenuation, f_ω is the random source field due to turbulent convection⁵, ψ_ω describes the pressure oscillations, and ∇ is the gradient with respect to position vector $x\in {\mathbb{R}}^{3}$. In this work we consider this model under an additional spherical symmetry assumption: c = c(|x|), ρ = ρ(|x|), γ = γ(|x|). We also assume that

$$\begin{equation}c\in {L}^{\infty }\left({\mathbb{R}}_{+}\right),\quad c{\geqslant}{c}_{\mathrm{min}}{ >}0,\quad \quad \rho \in {W}^{2,\infty }\left({\mathbb{R}}_{+}\right),\quad \rho { >}0,\quad \quad \gamma \in {L}^{\infty }\left({\mathbb{R}}_{+}\right),\end{equation} \tag{ 2 }$$

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

$$\begin{equation}c\left(r\right)={c}_{0},\quad \quad \rho \left(r\right)={\rho }_{0}\;\mathrm{exp}\left(-\left(r-{R}_{a}\right)/H\right),\quad \quad \gamma \left(r\right)=0,\quad r{\geqslant}{R}_{a},\end{equation} \tag{ 3 }$$

where R_a = R_⊙ + h_a, R_⊙ = 6.957 × 10⁵ km is the solar radius (the radius of the photosphere), 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 [3], which extends a standard solar model of [4] to the upper atmosphere. In this article we do not fix exact values of the above parameters, but recall that in [3] they are given by the following table:

In reality, the Sun is surrounded by a hot corona whose base is located at about h_c = 2000 km above the surface and which is highly inhomogeneous, for more details see [5]. Here we neglect this complication, and assume that the sound speed remains constant and that the density decays exponentially above the interface.

The exponential decay of density in the atmosphere leads to the trapping of acoustic waves with frequencies less than the cutoff frequency ω_ctf = c₀/(2H), which is about ω_ctf/2π ≈ 5.2 mHz for the Sun, and to the quantization of their admissible frequencies. These mode frequencies measured at the solar surface constitute the basic input data for helioseismological studies, see, e.g., [2].

On the other hand, acoustic waves with frequencies above the cutoff, that is, such that

$$\begin{equation}\omega { >}\frac{{c}_{0}}{2H},\end{equation} \tag{ 4 }$$

propagate away into the upper atmosphere. Several motivations for using measurements above the acoustic cutoff frequency for helioseismic inversions have been given in the literature. In particular, it has been shown that this high-frequency data can be used to study the sources of wave excitation [6], as well as the properties of the medium above the solar surface [7], see also [3,8].

The motions near the solar surface due to the solar acoustic waves can be measured through the Doppler shifts of an absorption line in the solar spectrum. Different instruments use different absorption lines. For example, the Michelson Doppler imager (MDI) instrument onboard the solar and heliospheric observatory (SOHO) satellite measures the line-of-sight Doppler velocity 200 km above the phostosphere, while the helioseismic and magnetic imager (HMI) instrument onboard the solar dynamics observatory (SDO) satellite observes solar oscillations at a height of 125 km.

In the present work we assume that the cross-covariance measurements of the solar acoustic field can be performed simultaneously at two different heights above the surface. This was the case during one year in 2010, when SOHO/MDI and SDO/HMI were both in operation. Also, the ground-based network GONG can be used in combination with the satellites. As recently shown in [9], it is also possible to perform multi-height measurements by combining six raw HMI filtergrams in different ways. In the present work we do not discuss details of measurements of the solar acoustic field, for more information on this subject see [1,10,11].

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.

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

$$\begin{equation}\left({L}_{v}-{k}^{2}\right)\psi =f,\quad {L}_{v}=-{\Delta}+v,\;k{ >}0,\end{equation} \tag{ 5 }$$

where the indices indicating dependence on ω are suppressed,

$$\begin{equation}k=\sqrt{\frac{{\omega }^{2}}{{c}_{0}^{2}}-\frac{1}{4{H}^{2}}},\quad v\left(x\right)={k}^{2}-\frac{\sigma {\left(\vert x\vert \right)}^{2}}{c{\left(\vert x\vert \right)}^{2}}+\rho {\left(\vert x\vert \right)}^{\frac{1}{2}}{\Delta}\left(\rho {\left(\vert x\vert \right)}^{-\frac{1}{2}}\right),\end{equation} \tag{ 6 }$$

$v\in {L}^{\infty }\left({\mathbb{R}}^{3}\right)$, and v(x) = 1/(H|x|) for |x| ⩾ R_a. Note that the term k² is included in the definition of v to guarantee that v decays at infinity.

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

$$\begin{equation}\begin{array}{c}v\left(x\right)=\tilde {v}\left(\vert x\vert \right),\quad x\in {\mathbb{R}}^{3},\\ \tilde {v}\in {L}^{\infty }\left({\mathbb{R}}_{+}\right),\quad \quad \tilde {v}\left(r\right)=\frac{\alpha }{r},\quad r{\geqslant}{R}_{a},\\ \text{for}\;\text{some}\;\text{constants}\;\alpha \in \mathbb{R},\quad {R}_{a}{ >}0.\end{array}\end{equation} \tag{ 7 }$$

If the potential v satisfies (7), then the resolvent ${\left({L}_{v}-{k}^{2}\right)}^{-1}$ is a meromorphic operator-valued function of $k\in {\mathbb{C}}_{+}=\left\{z\in \mathbb{C}\,:\Im z{ >}0\right\}$ with the distributional kernel G_v(x, x') = G_v(x, x'; 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') is called the radiation Green's function for equation (5). In addition, G_v(x, x') is a distributional solution to equation (L_v − k²)G_v(⋅, x') = δ_x', where δ_x' denotes the Dirac delta function centered at x'. We also suppose that $k\in {\mathbb{R}}_{+}{\backslash}{{\Sigma}}_{v}^{P}$, where

$$\begin{align}& {{\Sigma}}_{v}^{P}\subset {\mathbb{R}}_{+}\;\text{is}\;\text{the}\;\text{union}\;\text{of}\;\text{the}\;\text{sets}\;\text{of}\;\text{positive}\;\text{poles}\;k\\ & \quad \quad \quad \text{of}\;\text{functions}\;{G}_{v}\left(x,{x}^{\prime };k\right)\;\text{and}\;{G}_{\overline{v}}\left(x,{x}^{\prime };k\right).\end{align} \tag{ 8 }$$

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

$$\begin{equation*}{{\Vert}u{\Vert}}_{{L}_{\delta }^{2}}={\left[{\int }_{{\mathbb{R}}^{d}}{\left(1+\vert x\vert \right)}^{\delta }\vert u\left(x\right){\vert }^{2}\;\;\mathrm{d}x\right]}^{1/2},\quad \delta \in \mathbb{R}.\end{equation*}$$

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

$$\begin{equation}{\psi }_{v}\left(x\right)={\int }_{{\mathbb{R}}^{3}}{G}_{v}\left(x,{x}^{\prime }\right)f\left({x}^{\prime }\right)\;\;\mathrm{d}{x}^{\prime }.\end{equation} \tag{ 9 }$$

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}_{\ell }^{m}$:

$$\begin{equation}{\psi }_{v}\left(r\vartheta \right)=\frac{1}{r}\sum _{\ell {\geqslant}0}\sum _{\vert m\vert {\leqslant}\ell }{\varphi }_{v,\ell }^{m}\left(r\right){Y}_{\ell }^{m}\left(\vartheta \right),\quad \quad f\left(r\vartheta \right)=\frac{1}{r}\sum _{\ell {\geqslant}0}\sum _{\vert m\vert {\leqslant}\ell }{f}_{\ell }^{m}\left(r\right){Y}_{\ell }^{m}\left(\vartheta \right),\end{equation} \tag{ 10 }$$

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

$$\begin{equation}\left({L}_{v,\ell }-{k}^{2}\right){\varphi }_{v,\ell }^{m}={f}_{\ell }^{m},\quad {L}_{v,\ell }=-\frac{{\text{d}}^{2}}{\text{d}{r}^{2}}+\frac{\ell \left(\ell +1\right)}{{r}^{2}}+\tilde {v},\end{equation} \tag{ 11 }$$

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

$$\begin{equation}{\varphi }_{v,\ell }^{m}\left(r\right)={\int }_{0}^{{R}_{a}}{G}_{v,\ell }\left(r,{r}^{\prime }\right){f}_{\ell }^{m}\left({r}^{\prime }\right)\;\text{d}{r}^{\prime },\end{equation} \tag{ 12 }$$

where G_v,ℓ(r, r') is the coefficient in the spherical harmonics expansion of G_v:

$$\begin{equation}{G}_{v}\left(r\vartheta ,{r}^{\prime }{\vartheta }^{\prime }\right)=\frac{1}{r{r}^{\prime }}\sum _{\ell {\geqslant}0}\sum _{\vert m\vert {\leqslant}\ell }{G}_{v,\ell }\left(r,{r}^{\prime }\right){Y}_{\ell }^{m}\left(\vartheta \right)\overline{{Y}_{\ell }^{m}}\left({\vartheta }^{\prime }\right), {r}^{\prime },{\eta }^{\prime }{ >}0,\;\vartheta ,{\vartheta }^{\prime }\in {S}_{1}^{2}.\end{equation} \tag{ 13 }$$

One can show that G_v,ℓ is indeed a radiation Green's function for equation (11), see [12] and section 3.2.

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

Our first result relates cross correlations $\mathbb{E}\left(\overline{{\varphi }_{v,\ell }^{m}\left({r}_{1}\right)}{\varphi }_{v,\ell }^{m}\left({r}_{2}\right)\right)$ to the Green's function G_v,ℓ(r₁, r₂). We prove the following proposition:

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

$$\begin{equation}\mathbb{E}\left(\overline{{f}_{\ell }^{m}\left(r\right)}{f}_{\ell }^{m}\left(r+\xi \right)\right)={\Pi}{\delta }_{0}\left(\xi \right)\left(-\Im \tilde {v}\left(r\right)+k{\delta }_{R}\left(r\right)\right),\end{equation} \tag{ 14 }$$

for some Π > 0, R ⩾ R_a. Then the following formula is valid at fixed r₁, r₂ > 0:

$$\begin{equation}{\Pi}\;\Im {G}_{v,\ell }\left({r}_{1},{r}_{2}\right)=\mathbb{E}\left(\overline{{\varphi }_{v,\ell }^{m}\left({r}_{1}\right)}{\varphi }_{v,\ell }^{m}\left({r}_{2}\right)\right)+O\left(\frac{1}{R}\right),\quad R\to +\infty .\end{equation} \tag{ 15 }$$

Proposition 1 is proved in section 3.3. This proposition is a variation of a well-known result, see, e.g., [1, 16] 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\left(\frac{1}{R}\right)$ 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 \tilde {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)–(4). In this case k and v are given by formulas (6), and proposition 1 has a physical interpretation. Taking into account that $-\Im \tilde {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 established in fluctuation-dissipation theorems, see, e.g., [17].

Proposition 1 allows to retrieve Im G_v,ℓ(r, r) approximately from the power spectral density of noise ${\mathcal{P}}_{v,\ell }^{m}\left(r\right)=\mathbb{E}\vert {\varphi }_{v,\ell }^{m}\left(r\right){\vert }^{2}$ at fixed r. Next, we prove that Im G_v,ℓ(r, r) known exactly for all ℓ ⩾ 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}_{r}^{2}{\times}{S}_{r}^{2}$ at two different r.

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

• (a)  
  ${G}_{{v}_{1}}={G}_{{v}_{2}}$ on ${M}_{{R}_{o}}^{4}=\left({S}_{{R}_{o}}^{2}{\times}{S}_{{R}_{o}}^{2}\right){\backslash}{\Delta}$ for some R_o > R_a;
• (b)  
  $\Im {G}_{{v}_{1}}=\Im {G}_{{v}_{2}}$ on ${M}_{{R}_{o}}^{4}\cup {M}_{{R}_{o}^{{\dagger}}}^{4}$ for some ${R}_{o}^{{\dagger}}{ >}{R}_{o}{\geqslant}{R}_{a}$ such that ${R}_{o}^{{\dagger}}\notin {{\Sigma}}_{\alpha ,k,{R}_{o}}^{S}$, where ${{\Sigma}}_{\alpha ,k,{R}_{o}}^{S}\subset \left[{R}_{o},\infty \right)$ is a discrete set without finite accumulation points defined by (34b), (35).

Then v₁ = v₂ in ${L}^{\infty }\left({\mathbb{R}}^{3}\right)$.

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

$$\begin{equation*}{M}_{R}^{1}\left({x}_{1},{x}_{2}\right)=\left\{\left({x}_{1}\mathrm{sin}\;\theta +{x}_{2}\mathrm{cos}\;\theta ,{x}_{2}\right)\,:\,\theta \in \left(0,\pi \right]\right\},\end{equation*}$$

for some fixed x₁, ${x}_{2}\in {S}_{R}^{2}$ with x₁ ⋅ x₂ = 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²)ψ = 0 and establish auxiliary 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,ℓ of the radial equations. This relation will allow to extract the diagonal values G_v,ℓ(R, R) from G_v on ${M}_{R}^{4}$. Using this relation, in section 3.4 we show that the scattering matrix elements s_v,ℓ can be extracted from G_v on ${M}_{{R}_{o}}^{4}$ or from the imaginary part Im G_v only on ${M}_{{R}_{o}}^{4}\cup {M}_{{R}_{o}^{{\dagger}}}^{4}$, under the assumption that ${R}_{o}^{{\dagger}}\notin {{\Sigma}}_{\alpha ,k,{R}_{o}}^{S}$. In section 3.5 we prove that the scattering matrix elements s_v,ℓ determine the Dirichlet-to-Neumann map Λ_v,R for potential v in some ball ${B}_{R}^{3}$, R ⩾ R_a, where

$$\begin{equation*}{B}_{R}^{3}=\left\{x\in {\mathbb{R}}^{3}\,:\,\vert x\vert {< }R\right\}.\end{equation*}$$

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

Corollary 1. Let c, ρ, γ, and c', ρ', γ' be two sets of parameters satisfying (2), (3) and define the corresponding potentials v = v_ω, v' =v'_ω and the wavenumber k = k_ω according to formula (6) at fixed ω. Let ω₁ ≠ ω₂ be two positive frequencies satisfying (4) and such that ${k}_{\omega }\in {\mathbb{R}}_{+}{\backslash}\left({{\Sigma}}_{{v}_{\omega }}^{P}\cup {{\Sigma}}_{{v}_{\omega }^{\prime }}^{P}\right)$ for ω = ω₁, ω₂. Let ${G}_{{v}_{\omega }}$, ${G}_{{v}_{\omega }^{\prime }}$ be the radiation Green's functions at fixed ω for the potentials v, v' respectively.

Suppose that $\Im {G}_{{v}_{\omega }}=\Im {G}_{{v}_{\omega }^{\prime }}$ on ${M}_{{R}_{0}}^{4}\cup {M}_{{R}_{o}^{{\dagger}}}^{4}$ for some ${R}_{o}^{{\dagger}}{ >}{R}_{o}{\geqslant}{R}_{a}$ such that ${R}_{o}^{{\dagger}}\notin {{\Sigma}}_{1/H,{k}_{\omega },{R}_{o}}^{S}$, where ω = ω₁, ω₂. Then c = c', ρ = ρ', γ = γ' in ${L}^{\infty }\left({\mathbb{R}}_{+}\right)$.

Proof. Under the assumptions of corollary 1 it follows from theorem 1 that v_ω = v'_ω for ω = ω₁, ω₂. Using that $\mathfrak{R}{v}_{\omega }=\mathfrak{R}{v}_{\omega }^{\prime }$ for ω = ω₁, ω₂ one can show that c = c', ρ = ρ', see, e.g., the proof of [19, theorem 2.9]. Then, recalling that Im v_ω = −2ωγ/c², Im v'_ω = −2ωγ'/(c')², it follows from the equality Im v_ω = Im v'_ω for ω = ω₁, ω₂ together with the equality c = c' that γ = γ', concluding the proof of corollary 1. □

Remark 3. In corollary 1 it is assumed for simplicity of proofs and clarity of presentation that the measurements are performed in the homogeneous part of the solar atmosphere. This assumption can be relaxed (the proofs must be modified in a straightforward way) by requiring that the sound speed, density, and attenuation are known functions in the inhomogeneous part of the solar atmosphere, where the measurements are to be performed.

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

In this subsection we shall establish some auxiliary results regarding regular solutions of the radial Schrödinger equation which arises by separation of variables in the homogeneous equation (L_v − k²)ψ = 0.

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

$$\begin{equation}{\psi }_{v,\ell }^{m}\left(x\right)=\frac{1}{\vert x\vert }{\varphi }_{v,\ell }\left(\vert x\vert \right){Y}_{\ell }^{m}\left(\frac{x}{\vert x\vert }\right),\end{equation} \tag{ 16 }$$

which leads to the equation

$$\begin{equation}\left({L}_{v,\ell }-{k}^{2}\right){\varphi }_{v,\ell }=0,\end{equation} \tag{ 17 }$$

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

$$\begin{equation}{\mathrm{lim}}_{r\to +0}{r}^{-\ell -1}{\varphi }_{v,l}\left(r\right)=1,\end{equation} \tag{ 18 }$$

which fixes φ_v,ℓ uniquely.

Note that φ_v,ℓ(r) does not depend on the values of $\tilde {v}$ in the region [r, +∞), see [20, 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

$$\begin{equation*}\left({L}_{v}-{k}^{2}\right)\psi =0\quad \text{in}\text{\,\,}{B}_{R}^{3},\quad {\psi \mid }_{{S}_{R}^{2}}={Y}_{\ell }^{m},\end{equation*}$$

where ${Y}_{\ell }^{m}={Y}_{\ell }^{m}\left(x/\vert x\vert \right)$, $x\in {S}_{R}^{2}$, has a unique solution $\psi \in {H}^{2}\left({B}_{R}^{3}\right)$ if and only if φ_v,λ(R) ≠ 0 for all integer λ ⩾ 0. In addition, this solution is given by the formula

$$\begin{equation}\psi \left(x\right)=\frac{R}{\vert x\vert }\frac{{\varphi }_{v,\ell }\left(\vert x\vert \right)}{{\varphi }_{v,\ell }\left(R\right)}{Y}_{\ell }^{m}\left(\frac{x}{\vert x\vert }\right).\end{equation} \tag{ 19 }$$

Next we shall derive an expression for the regular solution φ_v,ℓ in the domain r ⩾ R in terms of the Coulomb wave functions and of the so-called scattering matrix element s_v,ℓ. First we recall the definition and some basic properties of the Coulomb wave functions from [21,  22].

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

$$\begin{equation}{H}_{\ell }^{{\pm}}\left(\eta ,kr\right)=\mathrm{exp}\left({\pm}\text{i}{\theta }_{\ell }\left(\eta ,kr\right)\right)+O\left(\frac{1}{r}\right),\end{equation} \tag{ 20a }$$

$$\begin{equation}{\theta }_{\ell }\left(\eta ,kr\right)=kr-\eta \;\mathrm{ln}\left(2kr\right)-\frac{1}{2}\ell \pi +{\sigma }_{\ell }\left(\eta \right),\end{equation} \tag{ 20b }$$

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

$$\begin{equation*}\frac{\partial }{\partial \rho }{H}_{\ell }^{{\pm}}\left(\eta ,\rho \right)=\left(\frac{\ell +1}{\rho }+\frac{\eta }{\ell +1}\right){H}_{\ell }^{{\pm}}\left(\eta ,\rho \right)-\sqrt{1+\frac{{\eta }^{2}}{{\left(\ell +1\right)}^{2}}}{H}_{\ell +1}^{{\pm}}\left(\eta ,\rho \right),\quad \ell {\geqslant}0\end{equation*}$$

given in [23], one can show that

$$\begin{equation}\frac{\partial }{\partial r}{H}_{\ell }^{{\pm}}\left(\eta ,kr\right)={\pm}\text{i}k{H}_{\ell }^{{\pm}}\left(\eta ,kr\right)+O\left(\frac{1}{r}\right),\quad r\to +\infty .\end{equation} \tag{ 21 }$$

Note that this derivation of formula (21) uses the identity

$$\begin{equation*}\left(\frac{\eta }{\ell +1}-\text{i}\right){\mathrm{e}}^{\text{i}\frac{\pi }{2}+\text{i}{\sigma }_{\ell }\left(\eta \right)-\text{i}{\sigma }_{\ell +1}\left(\eta \right)}-\sqrt{1+\frac{{\eta }^{2}}{{\left(\ell +1\right)}^{2}}}=0,\end{equation*}$$

which is a corollary of the recurrence relation for the gamma function:

$$\begin{equation*}{\Gamma}\left(z+1\right)=z{\Gamma}\left(z\right),\quad z\in \mathbb{C}{\backslash}\left\{0,-1,-2,\dots \;\right\}.\end{equation*}$$

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}}_{+}{\backslash}{{\Sigma}}_{v}^{P}$ be fixed and let η = α/(2k). Then the function φ_v,ℓ defined by (17), (18) admits in the region r ⩾ R_a the representation

$$\begin{equation}{\varphi }_{v,\ell }\left(r\right)={b}_{\ell }\left({H}_{\ell }^{-}\left(\eta ,kr\right)-{s}_{v,\ell }{H}_{\ell }^{+}\left(\eta ,kr\right)\right),\quad r{\geqslant}{R}_{a},\end{equation} \tag{ 22 }$$

with unique ${b}_{\ell }={b}_{\ell }\left(k\right)\in \mathbb{C}{\backslash}\left\{0\right\}$ and ${s}_{v,\ell }={s}_{v,\ell }\left(k\right)\in \mathbb{C}{\backslash}\left\{0\right\}$.

Proof. As the Coulomb wave functions ${H}_{\ell }^{{\pm}}\left(\eta ,kr\right)$ are linearly independent, any solution to equation (17) in the region r ⩾ R_a is given by their linear combination with unique coefficients. In particular, φ_v,ℓ can be expressed in the region r ⩾ R in the form

$$\begin{equation*}{\varphi }_{v,\ell }\left(r\right)={a}_{v,\ell }{H}_{\ell }^{+}\left(\eta ,kr\right)+{b}_{v,\ell }{H}_{\ell }^{-}\left(\eta ,kr\right),\quad r{\geqslant}R,\end{equation*}$$

for some a_v,ℓ, ${b}_{v,\ell }\in \mathbb{C}$. We shall show that a_v,ℓ ≠ 0, b_v,ℓ ≠ 0.

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

$$\begin{equation}{\int }_{\vert x\vert {\geqslant}1}{\left(1+\vert x\vert \right)}^{-1+\varepsilon }\vert \mathcal{D}\psi \left(x\right){\vert }^{2}\;\;\mathrm{d}x{< }\infty ,\quad \varepsilon \in \left(0,\frac{1}{2}\right],\end{equation} \tag{ 23a }$$

$$\begin{equation}\mathcal{D}\psi \left(x\right)=\nabla \psi \left(x\right)+\frac{x}{\vert x{\vert }^{2}}\psi \left(x\right)-\text{i}k\frac{x}{\vert x\vert }\psi \left(x\right).\end{equation} \tag{ 23b }$$

Now assume that b_v,ℓ = 0. Then it follows from (20a), (21) that the function defined by (19) is a non-zero solution of class ${H}_{\text{loc}}^{2}\left({\mathbb{R}}^{3}\right)\cap {L}_{-1-\varepsilon }^{2}\left({\mathbb{R}}^{3}\right)$, $\varepsilon \in \left(0,\frac{1}{2}\right]$, to (L_v − k²)ψ = 0 in ${\mathbb{R}}^{3}$ satisfying the radiation condition (23a). This contradicts the assumption $k\notin {{\Sigma}}_{v}^{P}$.

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

Remark 4. The coefficient s_v,ℓ = s_v,ℓ(k) in lemma 2 is called the ℓth scattering matrix element of the potential v.

In this subsection we shall express the radiation Green's function G_v,ℓ for equation (17) in the region r ⩾ R_a in terms of the Coulomb wave functions. Then we shall give a formula for extracting the diagonal values G_v,ℓ(R, R) from the radiation Green's function G_v for equation (5) known at ${M}_{R}^{4}$.

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

$$\begin{equation}{\varphi }_{v,\ell }^{+}\left(r\right)={H}_{\ell }^{+}\left(\eta ,kr\right),\quad r{\geqslant}{R}_{a}.\end{equation} \tag{ 24 }$$

The outgoing Green's function G_v,ℓ(r, r') for equation (17) is defined as a distributional solution to the equation (L_v,ℓ − k²)G_v,ℓ(⋅, r') = δ_r' specified by the following boundary conditions at fixed r' > 0:

$$\begin{equation*}{G}_{v,\ell }\left(r,{r}^{\prime }\right)=O\left({r}^{\ell +1}\right),\quad r\to +0,\quad \quad {G}_{v,\ell }\left(r,{r}^{\prime }\right)=c{H}_{\ell }^{+}\left(\eta ,kr\right),\quad r\to +\infty ,\end{equation*}$$

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

$$\begin{equation}{G}_{v,\ell }\left(r,{r}^{\prime }\right)=-\frac{{\varphi }_{v,\ell }\left({r}_{{< }}\right){\varphi }_{v,\ell }^{+}\left({r}_{{ >}}\right)}{\left[{\varphi }_{v,\ell },{\varphi }_{v,\ell }^{+}\right]},\quad r,{r}^{\prime }{ >}0,\end{equation} \tag{ 25 }$$

where r_< = min(r, r'), r_> = max(r, r') and $\left[{\varphi }_{v,\ell },{\varphi }_{v,\ell }^{+}\right]$ is the Wronskian:

$$\begin{equation}\left[{\varphi }_{v,\ell },{\varphi }_{v,\ell }^{+}\right]={\varphi }_{v,\ell }\left(r\right)\frac{\partial {\varphi }_{v,\ell }^{+}\left(r\right)}{\partial r}-\frac{\partial {\varphi }_{v,\ell }\left(r\right)}{\partial r}{\varphi }_{v,\ell }^{+}\left(r\right)\end{equation} \tag{ 26 }$$

By a standard theory of ordinary differential equations, this Wronskian is independent of r since the homogeneous second-order equation satisfied by φ_v,ℓ and ${\varphi }_{v,\ell }^{+}$ does not contain first-order terms, see [24, p 83, lemma 3.11]. Also note that formula (25) is a standard result from the theory of Sturm–Liouville problems, see, e.g., [24, p 158, formula (5.65)].

Lemma 3. If $k\in {\mathbb{R}}_{+}{\backslash}{{\Sigma}}_{v}^{P}$, then G_v,ℓ is well-defined and is given in the region r ⩾ R_a, r' ⩾ R_a by the formula

$$\begin{equation}{G}_{v,\ell }\left(r,{r}^{\prime }\right)=\frac{\text{i}}{2k}\left({H}_{\ell }^{-}\left(\eta ,k{r}_{{< }}\right)-{s}_{v,\ell }{H}_{\ell }^{+}\left(\eta ,k{r}_{{< }}\right)\right){H}_{\ell }^{+}\left(\eta ,k{r}_{{ >}}\right),\quad r,{r}^{\prime }{\geqslant}{R}_{a}.\end{equation} \tag{ 27 }$$

Proof. It follows from lemma 2 that under the assumption $k\in {\mathbb{R}}_{+}{\backslash}{{\Sigma}}_{v}^{P}$ the functions φ_v,ℓ(r) and ${H}_{\ell }^{+}\left(\eta ,kr\right)$ are linearly independent in the region r ⩾ R_a. Using the relation $\left[{H}_{\ell }^{-}\left(\eta ,kr\right),{H}_{\ell }^{+}\left(\eta ,kr\right)\right]=2\text{i}k$, given in [21], and formula (22) we get $\left[{\varphi }_{v,\ell },{\varphi }_{v,\ell }^{+}\right]=2\text{i}k{b}_{v,\ell }$. Together with (25), this implies (27). □

Next we shall show how the diagonal values G_v,ℓ(R, R) can be extracted from the Green's function G_v restricted to ${M}_{R}^{4}$.

First recall that the Legendre polynomials P_ℓ can be defined using the formal generating identity [25]

$$\begin{equation}\frac{1}{\sqrt{1-2st+{t}^{2}}}=\sum _{\ell =0}^{\infty }{P}_{\ell }\left(s\right){t}^{\ell }.\end{equation} \tag{ 28 }$$

Using the Laplace formula [25, theorem 8.21.2] and the Dirichlet convergence test one can show that at fixed t = 1 this series converges pointwise for all s ∈ (−1, 1). Besides, the Legendre polynomials form a complete orthogonal system in L²(−1, 1) such that

$$\begin{equation}{\int }_{-1}^{1}{P}_{\ell }\left(s\right){P}_{m}\left(s\right)\;\;\mathrm{d}s=\frac{2}{2\ell +1}{\delta }_{\ell m},\end{equation} \tag{ 29 }$$

where δ_ℓm is the Kronecker delta.

Now we recall [12] that at fixed $k\in {\mathbb{R}}_{+}{\backslash}{{\Sigma}}_{v}^{P}$ the Green's function G_v(x, x') is continuous outside of the diagonal x = x' and G_v(x, x') = O(|x − x'|⁻¹) as x → x'. Besides, the series expansion

$$\begin{equation}{G}_{v}\left(x,{x}^{\prime }\right)=\frac{1}{4\pi {R}^{2}}\sum _{l=0}^{\infty }\left(2\ell +1\right){G}_{v,\ell }\left(R,R\right){P}_{\ell }\left(x\cdot {x}^{\prime }/{R}^{2}\right),\quad \left(x,{x}^{\prime }\right)\in {M}_{R}^{4},\end{equation} \tag{ 30 }$$

converges for x ≠ ±x' and G_v,ℓ(R, R) has the asymptotics

$$\begin{equation}{G}_{v,\ell }\left(R,R\right)=\frac{R}{2\ell +1}\left(1+O\left(\frac{1}{\ell }\right)\right),\quad l\to +\infty .\end{equation} \tag{ 31 }$$

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

$$\begin{equation*}{P}_{\ell }\left(\vartheta \cdot {\vartheta }^{\prime }\right)=\frac{4\pi }{2\ell +1}\sum _{m=-\ell }^{\ell }{Y}_{\ell }^{m}\left(\vartheta \right)\overline{{Y}_{\ell }^{m}\left({\vartheta }^{\prime }\right)},\quad \vartheta ,{\vartheta }^{\prime }\in {S}_{1}^{2}.\end{equation*}$$

Lemma 4. Let $k\in {\mathbb{R}}_{+}{\backslash}{{\Sigma}}_{v}^{P}$ and fix x', ${x}^{{\prime\prime}}\in {S}_{R}^{2}$ such that x' ⋅ x^'' = 0. Then for each ℓ ⩾ 0

$$\begin{equation*}\begin{array}{l}\hfill {G}_{v,\ell }\left(R,R\right)=\frac{R}{2\ell +1}-2\pi {R}^{2}{\int }_{0}^{\pi }{g}_{v,R}\left(\mathrm{cos}\;\theta \right){P}_{\ell }\left(\mathrm{cos}\;\theta \right)\mathrm{sin}\;\theta \text{d}\theta ,\\ \hfill {g}_{v,R}\left(\mathrm{cos}\;\theta \right)={G}_{v}\left({x}^{\prime }\mathrm{cos}\;\theta +{x}^{{\prime\prime}}\mathrm{sin}\;\theta ,{x}^{\prime }\right)-\frac{1}{4\pi R}\frac{1}{\sqrt{2-2\;\mathrm{cos}\;\theta }},\;\;\theta \in \left(0,\pi \right).\end{array}\end{equation*}$$

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

$$\begin{equation*}{g}_{v,R}\left(s\right)=\frac{1}{4\pi {R}^{2}}\sum _{l=0}^{\infty }\left(2\ell +1\right)\left({G}_{v,\ell }\left(R,R\right)-\frac{R}{2\ell +1}\right){P}_{\ell }\left(s\right),\quad s\in \left(-1,1\right),\end{equation*}$$

which, in view of (29) and (31), also converges in L²(−1, 1). Recalling that Legendre polynomials form a complete orthogonal system in L²(−1, 1), we get lemma 4. □

In this subsection we prove proposition 1. First, recall that the Green's function G_v,ℓ satisfies the reciprocity relation

$$\begin{equation}{G}_{v,\ell }\left({r}_{1},{r}_{2}\right)={G}_{v,\ell }\left({r}_{2},{r}_{1}\right),\quad {r}_{1},{r}_{2}{ >}0.\end{equation} \tag{ 32 }$$

To prove it, consider the equations

$$\begin{equation*}\left({L}_{v,\ell }-{k}^{2}\right){G}_{v,\ell }\left(\cdot ,{r}_{1}\right)={\delta }_{{r}_{1}},\quad \quad \left({L}_{v,\ell }-{k}^{2}\right){G}_{v,\ell }\left(\cdot ,{r}_{2}\right)={\delta }_{{r}_{2}}.\end{equation*}$$

Multiplying the first equation by G_v,ℓ(⋅, r₂), subtracting the second equation multiplied by G_v,ℓ(⋅, r₁), and integrating over (0, R), R > r₁, r₂, we obtain

$$\begin{equation*}{\left.\left[{G}_{v,\ell }\left(\cdot ,{r}_{1}\right),{G}_{v,\ell }\left(\cdot ,{r}_{2}\right)\right]\right\vert }_{+0}^{R}={G}_{v,\ell }\left({r}_{1},{r}_{2}\right)-{G}_{v,\ell }\left({r}_{2},{r}_{1}\right),\end{equation*}$$

where [−, −] denotes the Wronskian defined according to (26) and the notation ${x\vert }_{a}^{b}=x\left(b\right)-x\left(a\right)$ is used. The next step is to show that the term on the left-hand side vanishes. Using formulas (18), (25) and the estimate $\frac{\partial }{\partial r}{\varphi }_{v,\ell }\left(r\right)=O\left({r}^{\ell }\right)$, r → +0, given in [12, theorem 3.3], we get [G_v,ℓ(⋅, r₁), G_v,ℓ(⋅, r₂)](+0) = 0. Using formulas (20a), (20b), (21), (25), we also get $\left[{G}_{v,\ell }\left(\cdot ,{r}_{1}\right),{G}_{v,\ell }\left(\cdot ,{r}_{2}\right)\right]\left(R\right)=O\left(\frac{1}{R}\right)$, R → +∞. As R tends to +∞, we get formula (32).

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

$$\begin{equation*}\left(\overline{{L}_{v,\ell }}-{k}^{2}\right)\overline{{G}_{v,\ell }\left(\cdot ,{r}_{1}\right)}={\delta }_{{r}_{1}},\quad \quad \left({L}_{v,\ell }-{k}^{2}\right){G}_{v,\ell }\left(\cdot ,{r}_{2}\right)={\delta }_{{r}_{2}}.\end{equation*}$$

Multiplying the first equation by G_v,ℓ(⋅, r₂), subtracting the second equation multiplied by $\overline{{G}_{v,\ell }\left(\cdot ,{r}_{1}\right)}$, integrating over (0, R), R > r₁, r₂, we get

$$\begin{align}& {\left.\left[\overline{{G}_{v,\ell }\left(\cdot ,{r}_{1}\right)},{G}_{v,\ell }\left(\cdot ,{r}_{2}\right)\right]\right\vert }_{+0}^{R}-2\text{i}{\int }_{0}^{R}\Im \tilde {v}\left(r\right)\overline{{G}_{v,\ell }\left(r,{r}_{1}\right)}{G}_{v,\ell }\left(r,{r}_{2}\right)\text{d}r\\ & \quad ={G}_{v,\ell }\left({r}_{1},{r}_{2}\right)-\overline{{G}_{v,\ell }\left({r}_{2},{r}_{1}\right)}.\end{align} \tag{ 33 }$$

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

$$\begin{equation*}\left[\overline{{G}_{v,\ell }\left(\cdot ,{r}_{1}\right)},{G}_{v,\ell }\left(\cdot ,{r}_{2}\right)\right]\left(R\right)=2\text{i}k\overline{{G}_{v,\ell }\left(R,{r}_{1}\right)}{G}_{v,\ell }\left(R,{r}_{2}\right)+O\left(\frac{1}{R}\right),\quad R\to +\infty .\end{equation*}$$

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

In this subsection we shall show that the scattering matrix elements s_v,ℓ for the potential v can be extracted from the Green's function G_v on ${M}_{{R}_{o}}^{4}$ or from its imaginary part Im G_v only on ${M}_{{R}_{o}}^{4}\cup {M}_{{R}_{o}^{{\dagger}}}^{4}$, where ${R}_{a}{\leqslant}{R}_{o}{< }{R}_{o}^{{\dagger}}$.

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

Lemma 5. Let v₁, v₂ be two complex-valued potentials satisfying (7), let $k\in {\mathbb{R}}_{+}{\backslash}\left({{\Sigma}}_{{v}_{1}}^{P}\cup {{\Sigma}}_{{v}_{2}}^{P}\right)$ be fixed and let R_o ⩾ R_a. Suppose that ${G}_{{v}_{1}}={G}_{{v}_{2}}$ on ${M}_{{R}_{o}}^{4}$. Then ${G}_{{v}_{1},\ell }\left({R}_{o},{R}_{o}\right)={G}_{{v}_{2},\ell }\left({R}_{o},{R}_{o}\right)$ and ${s}_{{v}_{1},\ell }={s}_{{v}_{2},\ell }$ for all ℓ ⩾ 0.

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

Note that formula (27) implies that

$$\begin{equation}\mathrm{cos}\;{\vartheta }_{\ell }\left(\eta ,kr\right)\mathfrak{R}{s}_{v,\ell }-\mathrm{sin}\;{\vartheta }_{\ell }\left(\eta ,kr\right)\Im {s}_{v,\ell }=\frac{{H}_{\ell }^{-}\left(\eta ,kr\right)-2k\Im {G}_{v,\ell }\left(r,r\right)}{\vert {H}_{\ell }^{+}{\left(\eta ,kr\right)\vert }^{2}},\end{equation} \tag{ 34a }$$

$$\begin{equation}{H}_{\ell }^{+}\left(\eta ,kr\right)=\vert {H}_{\ell }^{+}\left(\eta ,kr\right)\vert \mathrm{exp}\left(\text{i}{\vartheta }_{\ell }\left(\eta ,kr\right)/2\right),\end{equation} \tag{ 34b }$$

where r ⩾ R_a, η = α/(2k). Using (34a) with r = R_o, ${R}_{o}^{{\dagger}}$ one can see that s_v,l is uniquely determined from Im G_v,ℓ(R_o, R_o) and $\Im {G}_{v,\ell }\left({R}_{o}^{{\dagger}},{R}_{o}^{{\dagger}}\right)$ if and only if $\mathrm{sin}\left({\vartheta }_{\ell }\left(\eta ,k{R}_{o}^{{\dagger}}\right)-{\vartheta }_{\ell }\left(\eta ,k{R}_{o}\right)\right)\ne 0$. This justifies the definition of the singular set

$$\begin{equation}{{\Sigma}}_{\alpha ,k,{R}_{o}}^{S}=\left\{r\in \left[{R}_{o},+\infty \right)\,:\,\mathrm{sin}\left({\vartheta }_{\ell }\left(\eta ,kr\right)-{\vartheta }_{\ell }\left(\eta ,k{R}_{o}\right)\right)=0\quad \text{for}\text{\,\,}\text{some}\;\ell {\geqslant}0\right\}.\end{equation} \tag{ 35 }$$

Lemma 6. Let v₁, v₂ be two complex-valued potentials satisfying (7), let $k\in {\mathbb{R}}_{+}{\backslash}\left({{\Sigma}}_{{v}_{1}}^{P}\cup {{\Sigma}}_{{v}_{2}}^{P}\right)$, and let ${R}_{o}^{{\dagger}}{ >}{R}_{o}{\geqslant}{R}_{a}$ be such that ${R}_{o}^{{\dagger}}\notin {{\Sigma}}_{\alpha ,k,{R}_{o}}^{S}$. Suppose that $\Im {G}_{{v}_{1}}=\Im {G}_{{v}_{2}}$ on ${M}_{{R}_{o}}^{4}\cup {M}_{{R}_{o}^{{\dagger}}}^{4}$. Then ${s}_{{v}_{1},\ell }={s}_{{v}_{2},\ell }$ for all ℓ ⩾ 0. Besides, the set ${{\Sigma}}_{\alpha ,k,{R}_{o}}^{S}$ is discrete and does not have finite accumulation points.

Proof. Under the assumptions of lemma 6 it follows from lemma 4 that for all ℓ ⩾ 0 the equality $\Im {G}_{{v}_{1},\ell }\left(r,r\right)=\Im {G}_{{v}_{2},\ell }\left(r,r\right)$ 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 ℓ ⩾ 0. This concludes the proof of the first assertion of lemma 6.

Next we shall prove the second assertion. The following formulas are valid as ℓ → +∞ at fixed η and ρ ≠ 0:

$$\begin{equation*}\begin{array}{c}\hfill \Im {H}_{\ell }^{+}\left(\eta ,\rho \right)\sim {C}_{\ell }\left(\eta \right){\rho }^{\ell +1},\quad \quad \mathfrak{R}{H}_{\ell }^{+}\left(\eta ,\rho \right)\sim \frac{{\rho }^{-\ell }}{\left(2\ell +1\right){C}_{\ell }\left(\eta \right)},\hfill \\ \hfill {C}_{\ell }\left(\eta \right)\sim {\mathrm{e}}^{-\pi \eta /2}\frac{{\mathrm{e}}^{\ell }}{\sqrt{2}{\left(2\ell \right)}^{\ell +1}};\hfill \end{array}\end{equation*}$$

see formulas (33.2.11), (33.5.8), (33.5.9) of [21]. It follows that

$$\begin{equation}\mathrm{sin}\left({\vartheta }_{\ell }\left(\eta ,kr\right)-{\vartheta }_{\ell }\left(\eta ,k{R}_{o}\right)\right)\sim {\mathrm{e}}^{-1-\pi \eta }{\left(\frac{ekr}{2\ell }\right)}^{2\ell +1},\quad \ell \to +\infty ,\end{equation} \tag{ 36 }$$

locally uniformly in $r\in {\mathbb{R}}_{+}$. Besides, it follows from formula (33.2.11) of [21] and from discussion below it that at fixed ℓ ⩾ 0 the set of solutions $r\in {\mathbb{R}}_{+}$ of the equation sin(ϑ_ℓ(η, kr) − ϑ_ℓ(η, 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. □

In this subsection we shall show that the scattering matrix elements s_v,ℓ uniquely determine the Dirichlet-to-Neumann map for v is some ball ${B}_{R}^{3}$, R ⩾ R_a.

Note that lemma 1 justifies the definition of the singular set

$$\begin{equation*}{{\Sigma}}_{v,k}^{D}=\left\{R{ >}0\,:\,{\varphi }_{v,\ell }\left(R\right)=0\quad \text{for}\text{\,\,}\text{some}\;\ell {\geqslant}0\right\}.\end{equation*}$$

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

$$\begin{equation}\left({L}_{v}-{k}^{2}\right)\psi =0\quad \text{in}\text{\,\,}{B}_{R}^{3},\;{\psi \vert }_{{S}_{R}^{2}}=g,\end{equation} \tag{ 37 }$$

is uniquely solvable for $\psi \in {H}^{2}\left({B}_{R}^{3}\right)$. Besides, the set ${{\Sigma}}_{v,k}^{D}$ 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}\left({S}_{R}^{2}\right)$ then $R\notin {{\Sigma}}_{v,k}^{D}$.

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

$$\begin{equation*}{\psi }_{\ell }^{m}\left(x\right)={Y}_{\ell }^{m}\left(x/\vert x\vert \right){\int }_{{S}_{1}^{2}}\psi \left(\vert x\vert \omega \right){\overline{Y}}_{\ell }^{m}\left(\omega \right)\;\text{d}{S}_{\omega }\end{equation*}$$

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

Next we recall that the operator

$$\begin{equation}\psi {\mapsto}\left(\left({L}_{v}-{k}^{2}\right)\psi ,{\psi \vert }_{{S}_{R}^{2}}\right)\end{equation} \tag{ 38 }$$

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

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

$$\begin{equation}{\varphi }_{v,\ell }\left(r\right)={r}^{\ell +1}\left(1+O\left(\frac{1}{\ell }\right)\right),\quad \ell \to +\infty ,\end{equation} \tag{ 39 }$$

uniformly in r ∈ (0, R] at fixed R > 0. Besides, zeros of each φ_v,ℓ are discrete and do not have finite accumulation points. This concludes the proof of the second assertion of lemma 7. □

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

$$\begin{equation*}{{\Vert}\psi {\Vert}}_{{H}^{2}\left({B}_{R}^{3}\right)}{\leqslant}{C}_{v,k,R}{{\Vert}g{\Vert}}_{{H}^{3/2}\left({S}_{R}^{2}\right)}\end{equation*}$$

for some constant C_v,k,R > 0. In addition, the trace theorem [27] leads to the estimate

$$\begin{equation*}{{\Vert}\frac{\partial \psi }{\partial r}{\Vert}}_{{H}^{1/2}\left({S}_{R}^{2}\right)}{\leqslant}{C}_{v,k,R}^{\prime }{{\Vert}g{\Vert}}_{{H}^{3/2}\left({S}_{R}^{2}\right)}\end{equation*}$$

for some constant ${C}_{v,k,R}^{\prime }$C'_v,k,R > 0, where ${C}_{v,k,R}^{\prime }$ is the derivative of ψ in the radial direction.

Under the assumption that the Dirichlet problem (37) is uniquely solvable for $\psi \in {H}^{2}\left({B}_{R}^{3}\right)$ for all $g\in {H}^{3/2}\left({S}_{R}^{2}\right)$, we define the Dirichlet-to-Neumann map ${{\Lambda}}_{v,R}\in \mathcal{L}\left({H}^{3/2}\left({S}_{R}^{2}\right),{H}^{1/2}\left({S}_{R}^{2}\right)\right)$ by ${{\Lambda}}_{v,R}\varphi ={\frac{\partial \psi }{\partial r}\vert }_{{S}_{R}^{2}}$. Next we shall show that the partial scattering matrix elements s_v,ℓ known for all ℓ ⩾ 0 uniquely determine the Dirichlet-to-Neumann map Λ_v,R.

Lemma 8. Let v₁, v₂ be two complex-valued potentials satisfying (7) and let $k\in {\mathbb{R}}_{+}{\backslash}\left({{\Sigma}}_{{v}_{1}}^{P}\cup {{\Sigma}}_{{v}_{2}}^{P}\right)$ be fixed. Besides, let $R\in \left[{R}_{a},+\infty \right){\backslash}\left({{\Sigma}}_{{v}_{1},k}^{D}\cup {{\Sigma}}_{{v}_{2},k}^{D}\right)$. Suppose that ${s}_{{v}_{1},\ell }={s}_{{v}_{2},\ell }$ for all ℓ ⩾ 0. Then ${{\Lambda}}_{{v}_{1},R}={{\Lambda}}_{{v}_{2},R}$.

Proof. It follows from lemmas 1 and 2 and from continuity of φ_v,ℓ(r) at r = R that ${{\Lambda}}_{{v}_{1},R\vert {\mathcal{H}}_{\ell ,R}}={{\Lambda}}_{{v}_{2},R\vert {\mathcal{H}}_{\ell ,R}}$,where ${\mathcal{H}}_{\ell ,R}$ denotes the space of restrictions to ${S}_{R}^{2}$ of harmonic polynomials of degree ℓ, spanned by the spherical harmonics ${Y}_{\ell }^{m}={Y}_{\ell }^{m}\left(\frac{x}{\vert x\vert }\right)$, |m| ⩽ ℓ. More precisely, the following explicit formula for ${{\Lambda}}_{{v}_{j},R\vert {\mathcal{H}}_{\ell ,R}}$ is valid:

$$\begin{equation}{{\Lambda}}_{{v}_{j},R}{Y}_{\ell }^{m}=R\frac{{\left[\frac{\partial }{\partial r}\left(\frac{1}{r}{H}_{\ell }^{-}\left(\eta ,kr\right)\right)-{s}_{{v}_{j},\ell }\frac{\partial }{\partial r}\left(\frac{1}{r}{H}_{\ell }^{+}\left(\eta ,kr\right)\right)\right]}_{r=R}}{{H}_{\ell }^{-}\left(\eta ,kR\right)-{s}_{{v}_{j},\ell }{H}_{\ell }^{+}\left(\eta ,kR\right)}{Y}_{\ell }^{m},\end{equation} \tag{ 40 }$$

where η = α/(2k) and the denominator is non-zero as $R\notin {{\Sigma}}_{{v}_{j},k}^{D}$.

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

$$\begin{equation*}{g}_{N}=\sum _{\ell =0}^{N}\sum _{\vert m\vert {\leqslant}\ell }{g}_{\ell }^{m}{Y}_{\ell }^{m},\quad {g}_{\ell }^{m}={\int }_{{S}_{1}^{2}}g\left(R\omega \right){\overline{Y}}_{\ell }^{m}\left(\omega \right)\text{d}{S}_{\omega }.\end{equation*}$$

so that g_N → g in ${H}^{3/2}\left({S}_{R}^{2}\right)$ and, according to remark 5, ${{\Lambda}}_{{v}_{j},R}{g}_{N}\to {{\Lambda}}_{{v}_{j},R}g$ in ${H}^{1/2}\left({B}_{R}^{3}\right)$ as N → ∞. Together with (40), 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. □

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 6in case (b) that ${s}_{{v}_{1},\ell }={s}_{{v}_{2},\ell }$ for all ℓ ⩾ 0. Using lemma 8 we conclude that ${{\Lambda}}_{{v}_{1},R}={{\Lambda}}_{{v}_{2},R}$ for any R in the non-empty set $\left[{R}_{a},\infty \right){\backslash}\left({{\Sigma}}_{{v}_{1},k}^{D}\cup {{\Sigma}}_{{v}_{2},k}^{D}\right)$. It follows from the uniqueness theorem of [18], where the proof does not use that the potential is real-valued, that v₁ = v₂ a.e. This proves theorem 1.

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, ρ, γ) are perturbations of some known background quantities q⁰ = (c⁰, ρ⁰, γ⁰) such that

$$\begin{equation}\mathrm{s}\mathrm{u}\mathrm{p}\mathrm{p}\left(q-{q}^{0}\right)\subseteq I,\quad I=\left[{A}_{1},{A}_{2}\right]\subseteq \left(0,{R}_{\odot }\right],\end{equation} \tag{ 41 }$$

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

$$\begin{equation*}{{\Sigma}}_{q}^{P}=\left(0,\frac{{c}_{0}}{2H}\right]\cup \left\{\omega { >}\frac{{c}_{0}}{2H}\,:\,{k}_{\omega }\in {{\Sigma}}_{{v}_{\omega }}^{P}\right\},\end{equation*}$$

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

Put ${G}_{q}\left(h,\ell ,\omega \right)={G}_{{v}_{\omega },\ell }\left({R}_{\odot }+h,{R}_{\odot }+h\right)$. In view of lemma 4, as initial data for inversions from exact data we use the imaginary part of the Green's function Im G_q(h, ℓ, ω) measured at two different non-negative altitudes h ∈ {h₁, h₂}, at angular degrees ℓ ∈ {ℓ_min, ..., ℓ_max}, and at all admissible frequencies ω ∈ Ω. 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}\left(\ell ,\omega \right)={s}_{{v}_{\omega },\ell }$, ℓ ∈ {ℓ_min, ..., ℓ_max}, ω ∈ Ω, which are defined according to lemma 2. This reconstruction is done by considering equation (34a) with r = R_⊙ + h, h ∈ {h₁, h₂} as a linear system for finding $\mathfrak{R}{s}_{q}\left(\ell ,\omega \right)$, Im s_q(ℓ, ω) at each fixed ℓ, ω.

At the next step the scattering matrix elements s_q(ℓ, ω) are used to recover the map ${u}_{q}\,:\,I\to {\mathbb{R}}^{3}$(where I is the interval defined in (41)) which is defined as follows:

$$\begin{equation}\begin{array}{c}{u}_{q}=\left({u}_{q,1},{u}_{q,2},{u}_{q,3}\right),\\ {u}_{q,1}=\frac{1}{{c}_{0}^{2}}-\frac{1}{{c}^{2}},\quad {u}_{q,2}={\rho }^{\frac{1}{2}}\left(\frac{{\text{d}}^{2}}{\text{d}{r}^{2}}+\frac{2}{r}\frac{\text{d}}{\text{d}r}\right){\rho }^{-\frac{1}{2}}-\frac{1}{4{H}^{2}},\quad {u}_{q,3}=\frac{\gamma }{{c}^{2}},\end{array}\end{equation} \tag{ 42 }$$

and such that

$$\begin{equation*}{\omega }^{2}{u}_{q,1}+{u}_{q,2}-2\text{i}\omega {u}_{q,3}={v}_{\omega }.\end{equation*}$$

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

$$\begin{equation}F:{L}^{2}\left(I,{\mathbb{R}}^{3}\right)\to {\mathbb{C}}^{{\ell }_{\mathrm{max}}-{\ell }_{\mathrm{min}}+1}{\times}{\mathbb{C}}^{\vert {\Omega}\vert },\quad F\left({u}_{q}\right)={s}_{q},\end{equation} \tag{ 43 }$$

where |Ω| denotes the number of elements in the set Ω. For more details on the iteratively regularized Gauss–Newton method and for sufficient conditions of its convergence see, e.g., [29]. Note that we do not iterate directly over c, ρ, γ as it would require a more complicated forward map defined on a Sobolev-type space, see formula (42).

Remark 6. For constructing v_ω from the Dirichlet-to-Neumann map one can also use methods developed, for example, in [30]. However, in the present work we restrict ourselves to the iteratively regularized Gauss–Newton method.

The last step is to determine q = (c, ρ, γ) from u_q. Note that definitions (42) lead to the following explicit formulas for determining c and γ:

$$\begin{equation*}c={\left({c}_{0}^{-2}-{u}_{q,1}\right)}^{-\frac{1}{2}},\quad \quad \gamma ={c}^{2}{u}_{q,3}.\end{equation*}$$

Also note that definitions (42) lead to the following problem for determining ρ:

$$\begin{equation*}-\left(\frac{{\text{d}}^{2}}{\text{d}{r}^{2}}+\frac{2}{r}\frac{\text{d}}{\text{d}r}\right){\rho }^{-\frac{1}{2}}+\left({u}_{q,2}+\frac{1}{4{H}^{2}}\right){\rho }^{-\frac{1}{2}}=0,\quad \rho \left({A}_{1,2}\right)={\rho }^{0}\left({A}_{1,2}\right),\end{equation*}$$

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

We consider the background sound speed and density from the model of [3], which extends the standard solar model of [4] to the upper atmosphere. We suppose that the background attenuation is equal to γ₀/2π = 102.5 μHz inside of the Sun and decays to zero smoothly in the region [R_⊙, R_⊙ + h_a], where R_⊙ = 6.957 × 10⁵ 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 [1, section 7.3] for more details⁸. We assume that the unknown perturbations to the background values of solar parameters are supported in the interval [0.9R_⊙, 0.95R_⊙]. Besides, the background sound speed used in our reconstructions is by a factor of 1.1 larger in the region [0, 0.5R_⊙] than the sound speed used to generate data.

The initial data for reconstructions is the imaginary part of the radiation Green's function Im G_q(h, ℓ, ω) at heights⁹ h ∈ {105, 144} km, angular degrees ℓ ∈ {50, ..., 215}, and frequencies ω/2π ∈ {5.3, 5.4} mHz. Note that observations of the acoustic field at these heights approximately correspond to measuring Doppler velocities at line center (three-point approximation) and using the center of gravity (six-point approximation) of the HMI absorption line [9].

We consider perturbations δc^†, δρ^†, δγ^† of the sound speed, density and attenuation of maximal relative magnitudes of 5%, 9%, and 100%, respectively. The perturbations are represented by piecewise constant functions with 200 equidistant nodes. Figure 1 shows exact profiles of perturbations δc^†, δρ^†, δγ^†, reconstructed profiles of perturbations δc, δρ, δγ, and relative L² reconstruction errors e(δc, δc^†), e(δρ, δρ^†), e(δγ, δγ^†). Here,

$$\begin{equation*}e\left(f,{f}^{{\dagger}}\right)=\frac{{{\Vert}f-{f}^{{\dagger}}{\Vert}}_{{L}^{2}\left(I\right)}}{{{\Vert}{f}^{{\dagger}}{\Vert}}_{{L}^{2}\left(I\right)}}.\end{equation*}$$

This reconstruction example confirms the uniqueness results of section 2.2.

Zoom In Zoom Out Reset image size

Recall that the reconstructions are obtained using the data for angular degrees ℓ ∈ {50, ..., 215}. Note that:

• (a)  
  data for lower angular degrees are sensitive to perturbations in the region [0, 0.5R_⊙] where our background model is wrong.
• (b)  
  data for larger angular degrees are insensitive to perturbations in the region of interest [0.9R_⊙, 0.95R_⊙].

The reason is that the penetration depth of acoustic modes decreases with ℓ, as indicated by formula (39) and confirmed by plotting the regular solutions φ_v,ℓ for different ℓ, see figure 2.

Zoom In Zoom Out Reset image size

The power spectrum ${P}_{{v}_{\omega },\ell }^{m}\left(r\right)=\mathbb{E}\vert {\varphi }_{{v}_{\omega },\ell }^{m}{\left(r\right)\vert }^{2}$ introduced in section 2.1 can be approximated from the data. A standard approach is to parse the time series of acoustic oscillations into N segments of equal duration T, compute for each segment the sample power spectrum, and then take the arithmetic average ${\hat{P}}_{{v}_{\omega },\ell }^{m}\left(r\right)$.

The duration 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 1/T. Besides, for a time resolution (cadence) equal to ΔT, the sample power spectrum can be computed for frequencies up to the Nyquist frequency 1/(2ΔT).

It is a standard assumption going back to [31] (see also [32]) that

$$\begin{equation}\frac{2N{\hat{P}}_{{v}_{\omega },\ell }^{m}\left(r\right)}{{P}_{{v}_{\omega },\ell }^{m}\left(r\right)}\;\text{is}\text{\,\,}{\chi }^{2}\left(2N\right)\text{\,\,}\text{distributed},\end{equation} \tag{ 44 }$$

where χ²(2N) denotes the chi-squared distribution with 2N degrees of freedom and r, ω, ℓ, m are fixed. We use relation (44) to simulate noisy data. These simulated data are then used to extract the diagonal values of the imaginary part of the radial Green's function $\Im {G}_{q}\left(h,\ell ,\omega \right)=\Im {G}_{{v}_{\omega },\ell }\left({R}_{\odot }+h,{R}_{\odot }+h\right)$ using formula (15) with the $O\left(\frac{1}{R}\right)$ term dropped and assuming that Π = 1. Note that in reality Π is a function of frequency which is not directly accessible to measurements; for a possible model of this function see, e.g., [1, section 7.5].

Then the reconstruction proceeds as described in section 4.1.

We consider a model situation where the solar oscillations are observed for a total period of eight years with the time resolution of ΔT = 45 s. We recall that the HMI instrument observes the solar oscillations continuously with this cadence since April 30, 2010. We assume that the sample power spectra are computed from three-day intervals (that is, T = 3 × 24 × 3600s), so that the frequency resolution is equal to 1/T ≈ 3.86 μHz and the Nyquist frequency is 1/(2ΔT) ≈ 11.11 mHz.

We simulate ${\hat{P}}_{{v}_{\omega },\ell }^{m}\left({R}_{\odot }+h\right)$ according to (44) with N = 974 (which is the number of three-day intervals constituting eight years of observations) for angular degrees ℓ ∈ {50, ..., 250}, azimuthal degree m = 0, observation heights h ∈ {105, 144} km and six equidistant frequencies ω/2π from the range [5.27, 5.38] mHz. We use the same assumptions on the unknown parameters and the same background models for data generation and reconstructions and 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 (43) 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 3 shows parameters δc^†, δρ^†, δγ^†, their reconstructions δc, δρ, δγ from noisy data and relative L² reconstruction errors for one realization of data. The mean relative L² reconstruction errors for parameters δc^†, δρ^†, δγ^† are equal to

$$\begin{equation*}e\left(\delta c,\delta {c}^{{\dagger}}\right)=8.58\%,\quad \quad e\left(\delta \rho ,\delta {\rho }^{{\dagger}}\right)=18.32\%,\quad \quad e\left(\delta \gamma ,\delta {\gamma }^{{\dagger}}\right)=8.44\%.\end{equation*}$$

Besides, standard deviations of relative L² reconstruction errors are equal to 1.85%, 10.62%, 2.1%. The statistics are computed from 60 simulations for each parameter. 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.

Zoom In Zoom Out Reset image size

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

Recall that we use the data for six equidistant frequencies ω/2π in the range [5.27, 5.38] mHz. We do not use a larger range of frequencies because it would require to take into account the frequency dependence of attenuation which is not assumed by our model¹⁰. Our simulations also show that for a power law frequency dependence of attenuation

$$\begin{equation}\begin{array}{c}\gamma \left(\omega \right)={\gamma }_{0}{\left(\frac{\omega }{{\omega }_{0}}\right)}^{\beta }\\ {\gamma }_{0}/2\pi =4.29\;\text{mHz},\;{\omega }_{0}/2\pi =3\;\text{mHz},\;\beta =5.77,\end{array}\end{equation} \tag{ 45 }$$

as in [1], higher frequency waves have smaller penetration and are less sensitive to the medium properties in the region [0.9R_⊙, 0.95R_⊙], see figure 4. Note that the actual frequency dependence of attenuation is a subject of active research in helioseismology.

Zoom In Zoom Out Reset image size