Global uniqueness in a passive inverse problem of helioseismology

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.

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 ω > 0 proposed in [11]: where ξ ω is the spatial matter displacement vector, c is the sound speed, ρ is the density, γ is the attenuation, f ω is the random source field due to turbulent convection and x ∈ R 3 . In this work we consider this model under an additional spherical symmetry assumption: c = c(|x|), ρ = ρ(|x|), γ = γ(|x|). We also assume that c ∈ L ∞ (R + ), c ≥ c min > 0 a.e., ρ ∈ W 2,∞ (R + ), ρ > 0, γ ∈ L ∞ (R + ), (2) where R + = (0, ∞) and W 2,∞ (R + ) denotes the Sobolev space of functions defined on R + which belong to L ∞ (R + ) 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 [5, section 5.4]), and there are no attenuation: where R a = R + h a , R = 6.957 × 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: value meaning h a 500 km altitude at which the interface is located c 0 6855 m s −1 sound speed in the upper solar atmosphere ρ 0 2.886 × 10 −6 kg m −3 density at the interface H 125 km density scale height in the upper atmosphere In reality, the Sun is surrounded by a corona whose base is located at about h c = 2000 km 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 ω 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 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.
where the indices indicating dependence on ω are suppressed, In this article we consider equation (5) with a general complex-valued poten- for some constants α ∈ R, R a > 0.
If the potential v satisfies (7), then the resolvent (L v − k 2 ) −1 is a meromorphic operator-valued function of k ∈ C + = z ∈ C : z > 0 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 ∈ R + of the distributional kernel G v (x, x ) is called the radiation Green's function for equation (5). In ad- where δ x denotes the Dirac delta function centered at x . We also suppose that 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 denotes the Hilbert space of measurable functions u in R 3 with the finite norm Accordingly, for any k ∈ R + \ Σ P v equation (5) with f ∈ L 2 1+ε (R 3 ), ε ∈ (0, 1 2 ), admits a unique radiation (limiting absorption), solution given by 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 halfline R + with non-coupled channels. More precisely, consider the orthogonal expansions in normalized spherical harmonics Y m : where r > 0, ϑ ∈ S 2 1 and S 2 R = x ∈ R 3 : |x| = R}. Plugging these expansions into formula (5), we get the radial equations where ≥ 0, |m| ≤ . Besides, it follows from formulas (9), (10) that if ψ v is a unique radiation solution of equation (5), then ϕ m can be expressed as: where G v, (r, r ) is the coefficient in the spherical harmonics expansion of G v : One can show that G v, is indeed a Green's function for equation (11), that is, (L v, − k 2 )G v, (·, r ) = δ r , 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 ψ v is also a random function, as well as functions ϕ m v, and f m in the spherical harmonics expansions of formula (10). Following [11], we assume that the power spectrum (power spectral density) of ψ v , defined as P m v, (r) = E|ϕ m v, (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 , we assume that it can be measured at two different observation radii R † o > R o ≥ R . 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 E ϕ m v, (r 1 )ϕ m (r 2 ) to the Green's function G v, (r 1 , r 2 ). We prove the following proposition: Proposition 1. Let v be a complex-valued potential satisfying (7) and let k ∈ R + \ Σ P v be fixed. Assume that the random functions f m satisfy the condition for some Π > 0, R ≥ R a . Then the following formula is valid at fixed r 1 , r 2 > 0: Proposition 1 is proved in section 3.3. This proposition is a variation of a wellknown 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( 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 − v, and excited at the surface r = R with a power proportional to k. (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 − v = 2ωγ, 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].

Uniqueness results
Proposition 1 allows to retrieve G v, (r, r) approximately from the power spectral density of noise P m v, (r) = E|ϕ m v, (r)| 2 at fixed r. Next, we prove that 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 G v known on S 2 r × S 2 r at two different r.
Theorem 1. let v 1 , v 2 be two complex-valued potentials satisfying (7) and let k ∈ R + \ (Σ P v 1 ∪ Σ P v 2 ) be fixed. Assume that that one of the following conditions holds true: is a discrete set without finite accumulation points defined by (34b), (35).

Remark 2.
If v is some potential satisfying (7) then the restriction of G v (x, x ) to M 4 R depends on |x − x | 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 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 )ψ = 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, of the radial equations. This relation will allow to extract the diagonal values R . Using this relation, in section 3.4 we show that the scattering matrix elements s v, can be extracted from 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 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.
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

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 As the potential v is spherically symmetric, this equation separates in spherical coordinates. We seek a solution ψ m which leads to the equation together with the condition that ϕ v, vanishes at the origin. One can show that this determines ϕ v, ∈ H 2 loc (R + ) uniquely up to a multiplicative factor, see, e.g., [3]. We impose the boundary condition which fixes ϕ v, uniquely. Note that ϕ v, (r) does not depend on the values of v in the region [r, +∞), 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 if and only if ϕ v,λ (R) = 0 for all integer λ ≥ 0. In addition, this solution is given by the formula 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 [7,1].
The Coulomb wave functions H ± (η, kr), η = α/(2k), are the unique solutions of equation (17) with v(r) = α/r specified by the following asymptotics as r → +∞: where σ (η) = arg Γ( + 1 + iη) is the Coulomb phase shift and Γ denotes the usual gamma function. Functions H + (η, kr) and H − (η, kr) are complex conjugates of each other and are linearly independent. Using (20a) together with the relation given in [21], one can show that Using these properties, we shall prove the following result.

Lemma 2.
Let v be a complex-valued potential satisfying (7), let k ∈ R + \ Σ P v be fixed and let η = α/(2k). Then the function ϕ v, defined by (17), (18) admits in the region r ≥ R a the representation Proof. As the Coulomb wave functions H ± (η, kr) 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 for some a v, , b v, ∈ C. We shall show that a v, = 0, b v, = 0.
Recall from [22] that for any k ∈ R + outside the singular set Σ P v and for any f ∈ L 2 1+ε (R 3 ) the Schrödinger equation (5) admits the unique solution ψ ∈ 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 3 satisfying the radiation condition (23a). This contradicts the assumption k ∈ Σ P v . Now assume that a v, = 0 and let ψ be defined by (19). Then it follows from (20a), (21) , and satisfies (L v − k 2 )ψ = 0 in R 3 together with the radiation condition (23a). Taking into account definition (8), this also contradicts the assumption k ∈ Σ P v and concludes the proof of Lemma 2.

Remark 3.
The coefficient s v, = s v, (k) in Lemma 2 is called the -th scattering matrix element of the potential v.

Green's functions
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 functon G v for equation (5) known at M 4 R . In addition to the regular solution ϕ v, of equation (5) specified by the boundary condition (18), we consider the outgoing solution ϕ + v, which is specified by the asymptotics The outgoing Green's function G v, (r, r ) for equation (17) is defined as a distributional solution to the equation (L v, − k 2 )G v, (·, r ) = δ r specified by the following boundary conditions at fixed r > 0: for some non-zero constant c = c(r , η, k). If the regular solution ϕ v, (r) and the outgoing solution ϕ + v, (r) are linearly independent, this Green's function exists, is unique and is given by the explicit formula , r, r > 0, where r < = min(r, r ), r > = max(r, r ) and the Wronskian 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 ∈ R + \ Σ P v , then G v, is well-defined and is given in the region r ≥ R a , r ≥ R a by the formula Proof. It follows from Lemma 2 that under the assumption k ∈ R + \ Σ P v the functions ϕ v, (r) and H + (η, kr) are linearly independent in the region r ≥ R a . Using the relation [H − (η, kr), H + (η, kr)] = 2ik, given in [1], and formula (22) 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 ∈ (−1, 1). Besides, the Legendre polynomials form a complete orthogonal system in L 2 (−1, 1) such that where δ m is the Kronecker delta. Now we recall [3] that at fixed k ∈ R + \ Σ P v the Green's function converges for x = ±x and G v, (R, R) has the asymptotics Note that series expansion (30) follows from the spherical harmonics expansion (13), taking into account the well known addition theorem (see, e.g., [3]): Lemma 4. Let k ∈ R + \ Σ P v and fix x , x ∈ S 2 R such that x · x = 0. Then for each ≥ 0 Proof. Using (28) and (30) we get a pointwise convergent series expansion 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.

Extracting the Green's function from cross correlations
In this subsection we prove Proposition 1. First, recall that the Green's function G v, satisfies the reciprocity relation G v, (r 1 , r 2 ) = G v, (r 2 , r 1 ), r 1 , r 2 > 0.
To prove (15), we follow a similar scheme. We start from the equations Multiplying the first equation by G v, (·, r 2 ), subtracting the second equation multiplied by G v, (·, r 1 ), integrating over (0, R), R > r 1 , r 2 , we get In a similar way with the proof of formula (32), one can show that the Wronskian vanishes at zero and that Combining this with formulas (12), (14), (33), (32) to compute E ψ m (r 1 )ψ m (r 2 ) , we get formula (15), which concludes the proof of Proposition 1.

Recovering the scattering matrix elements
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 Recall that the Coulomb function H + (η, kr) does not vanish for r > 0, since H + (η, kr) and its complex conjugate H − (η, kr) form a basis of solutions of equation (17) with v(r) = α/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), Next we shall show that the scattering matrix elements s v, can be extracted from G v only, i.e., without knowing G v . However, the values of G v must be given not only on M 4 R o but also on M 4 H + (η, kr) = |H + (η, kr)| exp(iϑ (η, kr)/2 , where r ≥ R a , η = α/(2k).
This justifies the definition of the singular set Let v 1 , v 2 be two complex-valued potentials satisfying (7), . Then s v 1 , = s v 2 , for all ≥ 0. Besides, the set Σ S α,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 ≥ 0 the equality G v 1 , (r, r) = G v 2 , (r, r) holds true for r = R o , R † o . Together with the discussion before Lemma 6 it implies that s v 1 , = s v 2 , for all ≥ 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 locally uniformly in r ∈ R + . Besides, it follows from formula (33.2.11) of [1] and from discussion below it that at fixed ≥ 0 the set of solutions r ∈ 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.

Recovering the Dirichlet-to-Neumann map
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 3 R , R ≥ R a . Note that Lemma 1 justifies the definition of the singular set

Lemma 7.
Let v be a complex-valued potential satisfying (7) and let k ∈ R + be fixed. Then R ∈ R + \ Σ D v,k if and only if for any g ∈ H 3/2 (S 2 R ) the Dirichlet problem is uniquely solvable for ψ ∈ H 2 (B 3 R ). Besides, the set Σ 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 ∈ H 3/2 (S 2 R ) then R ∈ Σ D v,k . Now suppose that R ∈ Σ D v,k . First we shall show that (37) does not admit a non-zero solution for g = 0. Assuming that ψ ∈ 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 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 ψ m vanish and ψ = 0.
Next we recall that the operator [13]. Together with already established uniqueness for the Dirichlet problem (37) for R ∈ Σ D v,k , this proves the first assertion of Lemma 7.
Next we shall show that Σ D v,k is a discrete set without finite accumulation points. It can be shown [3] that the regular solution ϕ v, (r) defined by (17), (18) satisfies the estimate ϕ v, (r) = r +1 (1 + O( 1 )), → +∞, 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 4. Recalling from [13] that the operator (38) is Fredholm of index zero from
for some constant C v,k,R > 0. In addition, the trace theorem [13] leads to the estimate for some constant C v,k,R > 0, where ∂ψ ∂r = x |x| ∇ψ is the derivative of ψ in the radial direction.
Under the assumption that the Dirichlet problem (37) is uniquely solvable for ψ ∈ H 2 (B 3 R ) for all g ∈ H 3/2 (S 2 R ), we define the Dirichlet-to-Neumann map 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 1 , v 2 be two complex-valued potentials satisfying (7) and let k ∈ R + \ ( Proof. It follows from Lemmas 1 and 2 and from continuity of ϕ v, (r) at r = R that Λ v 1 ,R | H ,R = Λ v 2 ,R | H ,R , where H ,R denotes the space of restrictions to S 2 R of harmonic polynomials of degree , spanned by the spherical harmonics Y m = Y m ( x |x| ), |m| ≤ . More precisely, the following explicit formula for Λ v j ,R | H ,R is valid: where η = α/(2k) and the denominator is non-zero as R ∈ Σ D v j ,k . Now let g ∈ H 3/2 (S 2 R ) and denote by ψ v j ∈ H 2 (B 3 R ) the unique solution of the Dirichlet problem (37) with v = v j . Besides, define g N by so that g N → g in H 3/2 (S 2 R ) and, according to Remark 4, Together with (39), which shows that Λ v 1 ,R g N = Λ v 2 ,R g N , this implies that Λ v 1 ,R g = Λ v 2 ,R g and concludes the proof of Lemma 8.

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 , = s v 2 , for all ≥ 0. Using Lemma 8 we conclude that Λ v 1 ,R = Λ v 2 ,R for any R in the non-empty set [R a , ∞) \ (Σ D v 1 ,k ∪ Σ 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.

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, ρ, γ) are perturbations of some known background quantities q 0 = (c 0 , ρ 0 , γ 0 ) such that where R = 6.957 × 10 5 km is the solar radius, and we assume that both parameter sets q and q 0 satisfy (2), (3). Let Ω ⊂ R + be a finite set of admissible frequencies such that Ω ∩ (Σ P q ∪ Σ P q 0 ) = ∅, where v ω is defined according to (8), and the potentials v ω and v 0 ω 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, , ω) = G v ω , (R + h, R + h). As initial data for inversions from exact data we use the imaginary part of the Green's function G q (h, , ω) measured at two different non-negative altitudes h ∈ {h 1 , h 2 }, at angular degrees ∈ {0, . . . , max }, and at all admissible frequencies ω ∈ Ω. From this data we recover the solar parameters q as follows.
At the next step the scattering matrix elements s q ( , ω) are used to recover the map u q : I → R 3 (where I is the interval defined in (40)) which is defined as follows: u q = (u q,1 , u q,2 , u q,3 ), and such that ω 2 u q,1 + u q,2 − 2iωu q,3 = v ω . This reconstruction is done by applying the iteratively regularized Gauss-Newton method, going back to [4], to the forward map 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, ρ, γ) from u q . Note that definitions (41) lead to the following explicit formulas for determining c and γ: Also note that definitions (41) lead to the following problem for determining ρ: which is solved for the unknown function ρ − 1 2 . This step concludes the reconstruction algorithm from exact data.

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 γ 0 = 102.5 µHz inside of the Sun and decays to zero smoothly in the region [R , R + h a ], where R = 6.957 × 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 details 1 . We also assume that the unknown perturbations to the background values of solar parameters are supported in the interval [0.9R , 0.95R ].

Reconstruction scheme for noisy data
The power spectrum P m v ω , (r) = E|ϕ m v ω , (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 P m v ω , (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 ∆ω = 1/T. Besides, if the time resolution (cadence) is equal to ∆t then the sample power spectrum can be computed for frequencies up to 1/(2∆t).
It is a standard assumption going back to [26] that where χ 2 (2N) denotes the chi-squared distribution with 2N degrees of freedom and r, ω, , 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 G q (h, , ω) = G v ω , (R + h, R + h) using formula (15) with the O( 1 R ) 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., [11, section 7.5].
Then the reconstruction proceeds as described in section 4.1.

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 µHz and the maximal frequency is approximately 8.33 mHz. We simulate P m v ω , (R + h) according to (43) with N = 974 (which is the number of three-day intervals constituting eight years of observations) for angular degrees ∈ {0, . . . , 250}, azimuthal degree m = 0, observation heights h ∈ {105, 144} (km) and frequencies ω ∈ {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.
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.

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.