OUTGOING SOLUTIONS AND RADIATION BOUNDARY CONDITIONS FOR THE IDEAL ATMOSPHERIC SCALAR WAVE EQUATION IN HELIOSEISMOLOGY

In this paper, we study the time-harmonic scalar equation describing the propagation of acoustic waves in the Sun’s atmosphere under ideal atmospheric assumptions. We use the Liouville change of unknown to conjugate the original problem to a Schrödinger equation with a Coulomb-type potential. This transformation makes appear a new wavenumber, k, and the link with the Whittaker’s equation. We consider two different problems: in the first one, with the ideal atmospheric assumptions extended to the whole space, we construct explicitly the Schwartz kernel of the resolvent, starting from a solution given by Hostler and Pratt in punctured domains, and use this to construct outgoing solutions and radiation conditions. In the second problem, we construct exact Dirichlet-to-Neumann map using Whittaker functions, and new radiation boundary conditions (RBC), using gauge functions in terms of k. The new approach gives rise to simpler RBC for the same precision compared to existing ones. The robustness of our new RBC is corroborated by numerical experiments. Mathematics Subject Classification. 00A71, 33C15, 35J10, 35L05, 35A08, 35B40, 33C55, 85A20. Received August 21, 2019. Accepted December 6, 2019.


Introduction
In this article, we consider the scalar wave equation which models, in helioseismology, the propagation of acoustic waves in the interior and atmosphere of the Sun, Here is the density, is the sound speed. Absorption is modeled by means of a complex frequency ∈ C defined in terms of an absorption constant ≥ 0 and the real frequency > 0 so that 2 = (1 + ) 2 . The above equation is a Helmholtz equation with variable coefficients and is originally obtained under simplifying assumptions from a vectorial problem, cf. [15]. The value at the surface of the Sun of point-source responses in (1.1) is used to simulate the solar power spectrum ( −ℓ diagram) associated to pressure modes, see [13][14][15].
Our work focuses on the atmosphere of the Sun under ideal atmospheric assumptions. A starting point in helioseismology is to describe the interior of the Sun following the model S of [10], which assumes spherical symmetry [13]. Next, to model the atmosphere of the Sun, there are various simplifying models, the most basic one given by a constant extension from the model S ( [15], p. 2). In a more realistic approach, called ideal atmospheric behavior or the Atmo model (cf. [13], Sect. 2.3), the extension into the atmosphere is achieved by taking the density to decay exponentially at the same rate of the end of model S, while is smoothly extended to a constant. This atmospheric model was considered in [6] by the first author and collaborators for the purpose of constructing radiation boundary conditions when > 0. However, they did not have at their disposal the exact Dirichlet-to-Neumann (D-t-N) map or the asymptotic structure of the outgoing solution. These points will be addressed in our work.
Instead of working directly with (1.1) as was done in [6], we first carry out a change of unknown, called the Liouville transformation, = −1/2 u.  Before [2,7], either that the original form of the equation (1.1) was used, cf. [6,13], or the conjugated form (1.3) was employed, however without atmospheric model (i.e. on a bounded domain of the Sun) or with the assumption that the potential q is of compact support 5 , cf. equations (1.1), (1.2) and (2.1) of [1]. We consider two problems associated with the operator H away from the lower end of the spectrum, by requiring in general that, either ̸ = 0 or (︁ = 0 and A more natural assumption is ≥ 0. In the first problem, we extend the Atmo hypothesis to all of R 3 , and consider Our starting point is the function derived by Hostler and Pratt in [17] using the confluent hypergeometric special Whittaker functions, and the fact that it solves (1.7) with = 0 in a punctured domain R 3 ∖ { } for arbitrary . We use this kernel directly to extend to H − k 2 the standard results (e.g. well-posedness of outgoing solution, 3 In helioseismology, ( ) := 1 ( ) is called the density scale height, cf. equation (12) of [13]. 4 As noted, we are in the case of repulsive potential, since the density decreases strictly everywhere giving a strictly positive inverse density scale height . However, some comparison with attractive potential (i.e. when the coefficient of 1/ is negative) are also made in Remarks 2.3, 3.4 and 3.17. 5 This is a common approach for inverse problems dealing with compactly supported inhomogeneity. In fact, the Liouville transformation (1.2) was used in [28] to study the Calderón's problem (posed on a bounded domain) with an inhomogeneous conductivity, see also [4,26]. mapping property of resolvent, Green representation, radiation condition etc.) of Helmholtz operator −∆ − k 2 as found in, e.g. [11,31].
In the second problem, we return to the exterior problem (1.5), now with = 0, for purpose of constructing radiation boundary conditions (RBC), which by definition are boundary conditions placed on an artificial boundary in order to replace or approximate the radiation condition at infinity which chooses physical solutions 6 , cf. [5,8] and the references listed in [6]. We obtain expansion of general solutions in spherical harmonics and analytic expression of the D-t-N map in terms of Whittaker functions. The later is used as a reference map to evaluate the accuracy of the RBC. The overall analysis also shines light on the role played by k in the oscillation of the solutions in the atmosphere, and leads to simpler high-frequency RBCs (as opposed to those obtained in terms of / ∞ for the same precision in [6]).
While the Hostler-Pratt kernel in its closed form (and not just its expansion in spherical harmonics) is commonly referenced among physics literature (dealing with hydrogen atoms) e.g. equations (B2) and (B3) of [3], to the best of our knowledge, it is much less encountered in mathematics literature (with the exception of a brief mention in [16], p. 239), particularly in inverse scattering (in the spirit of [11]). Although the general spectral property of H is well-known, see discussion in Section 3, our work is the first to verify that the Hostler-Pratt function realizes the Schwartz kernel of the physical resolvent constructed by long-range scattering theory, and to use this explicit kernel to reobtain the aforementioned standard results (mapping properties, radiation condition, etc.). These results can find further-reaching applications beyond helioseismology in direct and inverse scattering (e.g. using layer potential technique). Regarding the second problem, as mentioned, an initial work [6] also constructed RBC, however using the original form of the equation and without the exact D-t-N. Our work is the first to give the exact D-t-N map, and to construct RBC using the energy level (wavenumber) of the operator. We also note that while Whittaker functions and Coulomb wave functions are well recognized in general literature, cf. [12,24,27] and the references therein, the usage of these hypergeometric functions rests rare in the context of helioseismology, appearing only in the recent works [2,7]. The later uses Coulomb wave functions, which differ from the Whittaker ones by a normalizing constant, see Remark 2.5.
The article is organized as follows. The first problem with Hostler-Pratt kernel is discussed in Section 3, the main results are stated in Propositions 3.10 and 3.14 which use as main ingredients the properties of Φ k given in Section 3.1.2. The second problem focusing on radiation boundary conditions is discussed in Section 4, where we include numerical experiments to show the performance of the new radiation conditions constructed with k.

Motivation for Whittaker functions
We decompose the solution of (1.7) in basis of spherical harmonics, We have a similar expansion for right-hand side (rhs) with coefficients ℓ . The coefficient ℓ satisfies the ODE, We further remove the first-order term, and work with unknown ℓ which solves the reduced ODE Then, (2.3) is obtained by using the usual technique in ODE to remove all-at-once the first order term −(2/ + ∞ ) d d . In our approach, we first remove − ∞ d (via the Liouville transformation) and then the first order term coming from the Laplacian. The unknowns in the above ODE are related by, We next introduce the change of variable The function ( ) defined by satisfies the ODE, When the right-hand side is zero, (2.8) is a special case of the Whittaker equation,

Notations
Complex frequencies. With > 0 and ∈ R, and under assumption (1.6), the complex frequency and modified wavenumber k are defined as, In defining expression (2.10), the square root branch is used to agree with the usual convention of outgoing in potential scattering theory. With this choice, for We work with k 2 ∈ C ∖ {0}, and distinguish two cases: when k 2 ∈ C ∖ [0, ∞), and when k = k 0 ∈ (0, ∞). Due to the choice (2.11), this range of k 2 corresponds to k ∈ C + ∪ (0, ∞). In terms of the physical parameters, k ∈ C + occurs when either ̸ = 0, or = 0 with ∞ < ∞ 2 , while k ∈ (0, ∞) occurs when = 0 with ∞ > ∞ 2 . The natural assumption ≥ 0 restricts k 2 to the first and second quadrant of its complex plane, i.e. Arg k 2 ∈ [0, ], which is equivalent to restricting k to the first quadrant of its complex plane, i.e. Arg k ∈ [0, 2 ]. When = 0 (in the absence of attenuation), k is denoted k 0 and can either be strictly positive, or purely imaginary (with strictly positive imaginary part).
When ∞ > ∞ 2 , we have the following one-sided limits, When ∞ < ∞ 2 , the limits from both sides agree, Special functions. To give the definition of the Whittaker functions, we will need the Gamma function Γ (cf.  Radial differential operator. We introduce (2.20)

Properties of Whittaker functions near zero
These limiting forms can be found in equations (13.14.16)-(13.14.19) of [30] , where is the Euler constant. In general, The expansion for the derivative is given in Proposition 30 of [7], For the Bucholtz function, we have, cf. Proposition 29 of [7],
(1) When k = k 0 > 0, they look like incoming/outgoing spherical waves. In particular, using the same convention of outgoing as in [11], (2) When Im k > 0, the functions exhibit oscillation coupled with attenuation in magnitude, and they are distinguished by their boundedness at infinity 9 , Using the analyticity of , ( ) with respect to , and the limits (2.16) and (2.17), when

Numerical illustration
We provide some numerical illustrations of the above asymptotic behavior of the functions The last function illustrates the radiating property of w. In Figure 1, we plot w andw and their respective reference functions exp( k ( )) defined in (2.35) and − 1 2 exp( k ( )) respectively, for with and without attenuation. The numerical computation of the Whittaker function, is discussed in Remark 4.7.
We have the following observations.
-The reference function k ( ) (and − 1 2 exp( k ( ))) gives an excellent representation of the functions w (anď w respectively). Surprisingly, does not need to be very large for this agreement to occur, both for a low and high frequency, with and without attenuation. This means that the reference functions also capture the case when w andw are only rapidly decaying (with no oscillation).
-In addition to its oscillatory behavior, for higher frequency /(2 ) = 20 Hz, the functioñ︀ w is also attenuating when = 0, cf. Figure 1b. This can be explained by the fact that̃︀ w is supposed to be an oscillatory term times one decaying in −1 , i.e. exp( ( )) × O( −1 ).
We further explore the decay rate of their symbol parts by factoring out the oscillatory part, and consider the following ratios, which are plotted in Figure 2: ; ; As expected, r 1 levels out to 1 + o(1), and r 2 to − 1 2 + o(1), at infinity, see Figure 2. As a result of this,̃︀ w tends to 0 modulo O( −1 ) at infinity (2.33).
This equation can also be identified with the Whittaker equation (2.8) by the change of variable = 2 2 , cf. equation (1.5) of [18]. A linearly independent basis near the origin is given by the real-valued ℓ ( , ) (called the regular Coulomb wave function) and ℓ ( , ) (also real-valued), with the first one regular at = 0, while the second one irregular. A pair of linear independent solutions in a neighborhood of infinity is given by the irregular functions ± ℓ ( , ) which are complex conjugates with ± ℓ = ℓ ± ℓ . Like the Whittaker functions, they are all defined from the Kummer functions, cf. footnote 8. The functions ± ℓ are normalized version of the Whittaker function , cf. 33.2.7 of [29] ± ℓ ( , ) = exp )︁ 1/2 . The oscillatory phase for the Coulomb funtions is more "balanced" between the incoming and the outgoing one compared to that of , cf. (2.35) and (2.36). In particular, These functions, instead of the Whittaker functions, were used in [2]. △

Physical solutions in R 3
In this section, we extend the hypothesis of the Atmo model to R 3 , and consider problem (1.7). The current assumption also means constant density scale height (1.4). The main goal is to use the Hostler-Pratt function Φ k derived in [17] to construct outgoing solutions and to establish radiation conditions that define the constructed solutions uniquely.
We start our discussion by citing the results in spectral theory concerning the invertibility of operator H ∞ , simply denoted as H, In the Atmo model, constant ∞ is positive, in which case, H ∞ is called a repulsive Halmiltonian. Like that of the Laplacian, the spectrum (H) of H is purely continuous, i.e.
cf. Section 10.2 of [33] and Chapter 10, p. 191 of [23] (with proof) or [16] (stated without proof). This means that for , called the resolvent at k 2 . For k 2 ∈ (0, ∞), one has the limiting absorption principle for H given by the machinery of long-range scattering theory 10 by Ikebe-Saito, cf. [34], which gives that, as k 2 approaches the spectrum [0, ∞), a limit of (H − k 2 ) −1 exists and is given as a bounded map, We unite the resolvent and its limit on the spectrum as the following mapping for k ∈ C + ∪ (0, ∞), Recall the definition of C + is given in (2.12), see also discussion in Remark 2.2. These results can also be rephrased as follows. For k with Im k > 0, and Furthermore, this solution is defined uniquely by a radiation condition associated with k 0 , as in case for the Helmholtz operator −∆ − k 2 0 . We also have The theory has to be slightly modified to take into account the weak singularity at the origin of the Coulomb potential ∞ | | , cf. [35]. For more details, see also discussion in Sections 3.2 and 3.3 from [7].
With the existence of the resolvent guaranteed, our first task is to describe explicitly the Schwartz kernel of ℛ(k) by using the Hostler-Pratt function Φ k stated below in (3.6). From [17], we will also make use of the fact that this function solves the homogeneous equation on punctured domains, cf. (3.7). Note that [17] only contains the derivation of the function without showing that it is actually a fundamental solution. To arrive at the final results, as in the case for Helmholtz equation, one requires the asymptotic expansion and radiating property of Φ k and ( ) Φ k ( , ), when stays in a bounded set and tends towards infinity. These properties of Φ k , stated in Propositions 3.5-3.7, are the main novelty of the section, since they are more technical to derive than those for the outgoing Helmholtz kernel, due to the lack of invariance under translation. However, once these ingredients are obtained, in a second task, similar techniques as in the Helmholtz case are employed to construct outgoing solutions and its characterizing radiation conditions, Proposition 3.12. The latter gives an alternative proof for uniqueness of outgoing solutions to that given by long-range scattering theory. When there are similarities with Helmholtz case, we only include the main statements 11 . We also make some remarks concerning the attractive potential in Remarks 3.4 and 3.17. The adaptation of the main results for the original problem (1.1) is listed in Remark 3.15.

Physical Green kernel
For k 2 ∈ C ∖ {0} (or equivalently k ∈ C + ∪ (0, ∞)), and for all ∈ R 3 , we look for the fundamental solution Φ k to the conjugated problem, Our starting point is the derivation in [17] which gives the explicit form of Φ k , and that it satisfies the homogeneous equation on punctured domains.

Definition
With defined in (2.14), when k ∈ C + , following equation (8) from [17], we define Φ k ( , ) as, We recall the definition of the auxiliary distances and in (2.19) and the Whittaker functions in (2.24) and (2.25). We will assume, from the work of [17], that , 0), and the definition in (3.6) is well-defined for k 2 on the negative axis, or equivalently k on the positive imaginary axis. We next extend the definition in [17], which is only given for Im k > 0, cf. equation (1) of [17], to k ∈ (0, ∞). When k = k 0 > 0, using (2.16), we define the outgoing kernel Φ + k0 and incoming one Φ − k0 as In particular, the k 0 -outgoing kernel is defined as With (3.9), we have extended the definition of the Hostler-Pratt kernel to all k ∈ C + . Since the natural assumption is ≥ 0, cf. Remark 2.2, we will only work with the function and its outgoing extension. For the rest of the article, for k ∈ C + ∪ (0, ∞), Φ k denotes the physical (also called outgoing) kernel given by (3.6) which also includes (3.9) when k ∈ (0, ∞), since ≥ 0.

Properties
We introduce the reduced Green function = k , The outgoing kernel is then , with defined in (2.14).
This is confirmed in equation (9) from [17]. This current form of the fundamental solution when the source position = 0 is also obtained by Meixner, see further discussion of literature and history in [17]. This can also be obtained directly by separation of variables, cf. Section 4.3.1 of [7] The following proposition gives the limiting value for fixed and → , which determines the local integrability of the fundamental kernel, see Appendix A for the proof of the proposition. As a result of this, (3.14) In addition, for each fixed ∈ R 3 , we have This means that the function k 2 ↦ → Φ k has poles at k 2 = − 2 ∞ 4( +1) 2 only when ∞ < 0, which are the values of k 2 at which the equation (−∆ − k 2 + ∞ /| |) = 0 has nonzero 2 eigenfunction, given in Remark 3.17. When ∞ ≥ 0, Φ k is pole-free for k 2 on C ∖ {0}. Equivalently, in terms of k, we have We next consider the asymptotic expansion when one variable goes to infinity and the other stays in a compact set. In addition to the phase function introduced in (2.35), we introducẽ k, ( ) : With auxiliary variable defined in (2.19), the leading factor in the asymptotic of Φ k is given by k ( ) , which can be approximated by For more discussion on this approximation see Proposition 10, Figures 9 and 10 from [7]. We now describe the full structure of Φ k ( , ) including its symbol part, with explicit description of the leading constant. The proof for the proposition is written up in Appendix A.
for in a compact set, when | | → ∞. Here, the constant C( ,ˆ) is defined as, (3.20) In the above expression, the remainder O(| | −1 ) is uniformly bounded for all in a bounded set and all directions ofˆ. As a result, when k ∈ (0, ∞), for fixed , the function ↦ → Φ( , ) is exponentially decaying, and thus is in 2 (R 3 ).
Numerical illustration of Proposition 3.5. We illustrate in Figure 3 the function Φ k ( , 0 ) for a fixed source in 0 = (0, 0, 1) at = (0, sin( ), cos( )), along two directions of approaching infinity along polar angle = 6 and 3 2 in the ( , )-plane. This also means that = /2. For the sake of conciseness, we only show the plot with attenuation, = 1, and see that the oscillatory behavior is coupled with a decrease of the signal envelope, Figure 3a (in the case without attenuation, the envelope does not decrease and the phase exhibits pure oscillation, see [7]).
To investigate the decay rate of the symbol part, we factor out the leading oscillatory part of the reduced Green function k ( , ) (3.10), and plot the following ratios (as functions of ),︀ The resulting functions are plotted in Figure 3b and show the uniformity of the asymptotic expansion with respect to the direction of approaching infinity. All three remainders have the profile of a constant +o(1) (it is similar without attenuation). However, the convergence to the horizontal asymptotes forr 1 andr 3 is much earlier than what is exhibited in the curve forr 2 . Thus, exp( k ( )) and exp(˜k , ( )) give better representation of k , with smaller error and faster convergence. The next proposition states the radiating property of the fundamental solution Φ + k0 ( , ), as | | → ∞ and in a compact set. Although a bit more involved than that for Proposition 3.3, its proof uses the same techniques and as main ingredient the radiating asymptotics (2.33) of − ,1/2 , cf. Propositions 12 and 13 of [7]. Proposition 3.6 (Radiating asymptotic). When

23)
uniformly for in a compact set and | | → ∞.
When k 2 ∈ C∖[0, ∞), we verify that the Φ k ( , ) satisfies the hypothesis of Schur's test, cf., e.g. Lemma 3.7.4, p. 75 from [22], which will be used to show that Φ k gives rise to an 2 (R 3 ) bounded operator. The proof for the proposition is written up in Appendix A.

Verification of the fundamental solution
To show that Φ k is a fundamental solution, by definition of (3.5), we have to show that with The main ingredient is (3.7) assumed from the derivation of Φ k in [17]. A second ingredient is the second Green's formula: for a bounded region ℛ and , ∈ 2 (ℛ), where ( ) is the normal vector along ℛ and points outward of ℛ. Note that, from Proposition 3.3, we already have, when k 2 ∈ C ∖ [0, ∞), Φ k ( , ) ∈ 2 (R 3 ) for in a compact subset of R 3 and vice versa (due to the symmetry in and ). Proof.
With Ω ( , ) defined in (3.25), and using (3.7) on Ω ( , ) with > 0, we can write with ( ) the normal vector pointing outside of Ω ( , ) . Here, we have defined (using a spherical change of variable centered at , i.e. = + with = | |), Here, S (0,1) is the unit sphere centered at the origin. We next consider the limits of each integral as → 0 + . We use (3.14) to obtain that lim →0 I 1 = 0, and (3.15) to obtain that Putting together these two limits, we have By the same technique, we obtain the Green representation on a bounded domain. The proof uses the same ingredient of that for Proposition 3.8. For details, see Proposition 16 of [7].

Construction of physical solutions
We now have the necessary ingredients to construct physical solution to (3.2), For k 2 ∈ C ∖ [0, ∞), from Proposition 3.7, the integral operator defined below in (3.29), with Φ k as its Schwartz kernel, is 2 (R 3 )-bounded. In another word, (3.29) defines a 2 (R 3 ) function for ∈ 2 (R 3 ). When k 0 ∈ (0, ∞), from Proposition 3.5 and Proposition 3.3, for fixed , ↦ → Φ k0 ( , ) ∈ 2 loc (R 3 ) (and the statement is symmetric in and ), thus the integral in (3.29) is well-defined for ∈ 2 comp (R 3 ). Additionally, (3.29) defines analytic functions due to the analyticity of Φ k , cf. (3.16), thus is 2 loc (R 3 ) in both cases. That (3.29) defines a solution follows from Proposition 3.8. It remains to use standard arguments (e.g. Dominated Convergence Theorem) to justify differentiation past the integral sign to show that the solution defined by (3.29) inherits the radiating property of Φ k , the latter given by Proposition 3.6. These are the essential ingredients to Proposition 3.10. The more pedantic details of the proof can be found in Proposition 17 of [7].
is well-defined and provides a solution to (3.28). When k ∈ C + , the solution is 2 (R 3 ). When k ∈ (0, ∞), the solution is 2 loc (R 3 ) and satisfies the uniform Sommerfeld radiation condition, At this point, using the uniqueness guaranteed by spectral theory, we can conclude that the Schwartz kernel of ℛ(k) is given by Φ k . In particular, for ∈ 2 comp (R 3 ) if k ∈ (0, ∞), or for ∈ 2 (R 3 ) if k ∈ C + , we have Remark 3.11. These results hold true for in a weighted Sobolev space 2 (R 3 ) with > 1/2. △

Radiation conditions and uniqueness
In addition to the uniqueness result by means of spectral theory, the uniqueness of 2 (R 3 ) solution can be established by basic techniques involving coercitivity of the corresponding variational form (cf. e.g. [7], Prop. 1). When k 0 > 0, we can also re-obtain uniqueness, as was done for Helmholtz equation (cf. e.g. [31], Thm. 2.3), by working directly with the kernel Φ k0 (= Φ + k0 ). Due to the positive Coulomb potential, the only change is the modified Green's formula, (1) Uniform Sommerfeld radiation condition: (3.33) (2) 2 -radiation condition: (3) Exterior representation formula: The details of the proof can be found in our report ( [7], Prop. 18).
Definition 3.13. For k 0 > 0, a solution ∈ 2 loc (R 3 ) to (3.28) at k 2 0 in the exterior of bounded region Ω is called outgoing if satisfies one of the conditions listed in Proposition 3.12. △ Proposition 3.14 (Uniqueness). With "outgoing" given by Definition 3.13, and k 0 > 0 and ∈ 2 loc (R 3 ), we have the following statements. Proof. If is an outgoing solution in R 3 , then by the exterior representation formula (3.35), )︁ d ( ), (3.36) in R 3 ∖ B (0, ) for all < | |. Let → 0, we obtain that ( ) → 0. For the second statement, suppose 1 and 2 are two outgoing solutions of the inhomogeneous equation, then = 1 − 2 solves the homogeneous equation and is still outgoing; using the first part, ≡ 0. Due to the exponential decay of , for f ∈ 2 (R 3 ), 1/2 f ∈ 2 (R 3 ) with > 1 2 , and is thus in the domain of ℛ(k) by Remark 3.11. The unique outgoing solution to (1.1) with right-hand-side f ∈ 2 (R 3 ) is given by, Remark 3. 16. In conform with the results given by spectral theory discussed at the beginning of the section, the above results confirm the analyticity of the resolvent C + : k ↦ → ℛ k , which is the equivalent to the analyticity of C + ∋ k ↦ → Φ k cf.

}︁ .
These are called the eigenvalues of the hydrogen atom, and the corresponding eigenfunctions are, cf. Theorem 10.10 of [33], Here L is the generalized Laguerre polynomials, which can be defined in terms of the Whittaker functions, cf.

Radiation Boundary Conditions (RBC)
In this section, we consider the problem (1.5) with a focus on radiation boundary conditions. Recall from (2.1) that the original equation (1.1) on | | ≥ r is equivalent to the conjugate problem (4.1), This equivalence on each spherical mode is between the original modal ODE (2.4) and the conjugate modal ODE (2.2). The latter is equivalent to the reduced ODE (2.3) which is what we use to derive the exact D-t-N and new RBC. From relations (2.5), the following remark relates these three approaches when in used with a Robin-type absorbing boundary condition.
Remark 4.1 (Relations among the modal RBC). We denote by Z ℓ • a radiation impedance coefficient in the Robin-type boundary condition imposed at | | = , with ≥ r , used to truncate the reduced ODE (2.3). This approach is equivalent to using the RBC, and the modal ODE (2.2) of the conjugate problem (4.1), and to working with 4) and the modal ODE (2.4) of the original problem (1.1). This last approach was used in [6] for spherical symmetry and [13] for axis symmetry. See also Remark 4.6 which explains why working with the conjugate and ultimately the reduced ODE are the most natural. △ We first obtain the expansion of a generic solution in spherical harmonics, using the Whittaker functions, which gives readily the exact D-t-N map, stated in Proposition 4.4, and the structure of the outgoing solution. This gives a new perspective in approximating the transparent boundary condition listed in (4.12), in working with the reduced ODE (2.3) and the right gauge function in terms of k. In (4.3), we compare the new radiation impedance coefficients with those obtained in [6,13].

General solutions and exact D-t-N map
In (2.1), the conjugate ODE (2.2) is written as Whittaker equation (2.8). A basis of linearly independent solutions on a neighborhood of infinity of the Whittaker equation is given by the pair of Whittaker functions (2.22). This fact readily gives the expansion of a generic solution of (4.1) in spherical harmonics, by the same reasoning used for Helmholtz, cf. Theorem 1.1 and Appendix 1 from [32].
Remark 4.3. The above result is also used to prove Rellich-type lemma, cf. Theorem 1.2 of [24] or Lemma 20 of [7], which ensures the uniqueness of solution to exterior boundary value problem associated with operator Recall the outer D-t-N operator T + : 1/2 (Γ) → −1/2 (Γ) associated with (1.5) and boundary Γ is defined as, T + = | + Γ where is the (unique) outgoing solution to (1.5) on R ∖ Ω which satisfies boundary condition | − Γ = , cf. e.g. Section 3.1.2 of [19]. We write where ℓ are the coefficients of the expansion of | Γ in the spherical harmonic basis Y ℓ , (2.1); and ℓ are called the modal D-t-N coefficients of T + at level ℓ. By standard computation, the expansion (4.5) gives readily the expression of ℓ . For details of the proof, see Proposition 22 of [7].

List of modal radiation impedance coefficients
In the current discussion, k ∈ C + ∪ (0, ∞). As reference radiation conditions, we have at our disposal Z ℓ DtN which comes from the modal D-t-N coefficients ℓ (4.7), In addition we have the nonlocal modal radiation coefficient Z ℓ nonlocal which was introduced in Section 3.1 of [6] to define (truncated) outgoing radiating wave. Working with the reduced ODE (2.3) and k, we have, cf. Section 6.2 from [7], Lack of the true D-t-N, a numerical approximation was employed in [6] to create a reference solution by placing a homogeneous Dirichlet condition at a very large radius. This technique, only applicable in the case of absorption, creates stability problem for very small absorption.
To approximate the Z ℓ nonlocal , we use the same ideas as in [6], however working with small quantities defined in terms of k and not / ∞ . We also distinguish between the high-frequency family HF, taking as small quantity k −1 , and the small angle of incidence family SAI, taking as small quantity ℓ(ℓ + 1)/( k) 2 . The HF approach yields the three following conditions, See Remark 28 of [7] for another way to come about Z S-HF-0 . The SAI family in k coincides with that in / ∞ , and are denoted here as SAI, The main goal is to see how the new HF conditions S-HF perform compared to the old conditions, labeled as A-HF, which were obtained by taking ( / ∞ ) −1 as small quantity in [6], )︂ · (4.14) In our comparison, we will also include Z A-RBC-1 := 1 + Z S-HF-0 , (4.15) and the "naive Sommerfeld" (in our notation)

Numerical experiments
We investigate the performance of the RBC listed in (4.2), in particular that of the S-HF families, to approximate the exact outgoing solutions. We work with the reduced ODE (2.3) and consider the truncated outer Dirichlet problem on r min ≤ ≤ : Here Z • is one of the RBC introduced above, and the corresponding solution is labeled as Z•,ℓ, or simply Z• . We will compare these with the exact solution, which is the unique solution to (4.17) with Z ℓ DtN , given explicitly by, Remark 4.7. For the resolution of (4.17), we employ a finite difference discretization of order four with a parallel code written in Fortran. The computation of the complex-valued Whittaker function relies on the library arb 12 , [20], for efficient computation with high precision of hypergeometric functions [21]. The use of a dedicated computational library is required, as we have observed instability with the intrinsic Whittaker function in Matlab, in particular near the cut-off frequency (4.21). △ In our experiments, we work with varying ℓ, and , and the following fixed parameters: First, for a chosen ℓ = 3 and /(2 ) = 20 Hz, we plot Z• on [r min , ] against the reference solution ref given in (4.18). For the sake of conciseness, we only plot the case without attenuation ( = 0), in Figure 4 (the case with attenuation is similar, see [7]). It is hard to discern any differences among the solutions, except those corresponding to Z A-HF-0 and Z Naive which give a poor approximation and low accuracy. In order to have a better comparison, we compute the solutions for a range of ℓ from 0 to 180 and frequencies /(2 ) from 1 to 100 Hz (i.e. 18 100 test-cases, doubled due to the presence or absence of attenuation). To evaluate the performance of the RBC, we consider the global 2 error and a local error at a fixed point , In Table 1, we list the mean of the error E over the investigated range of frequencies and ℓ, from which we obtain the following order in terms of accuracy: From these results, we have the following observations.
(1) Each Z • yields lower error in the presence of attenuation than without, and the errors remain small for ℓ ≥ 50 and below the cut-off frequency (4.21), (2) In accordance with the comparison in [6], we observe that the coefficients dependent on ℓ behave better than the independent ones, since they are higher order approximations of the nonlocal one, and that SAI family (small angle approximation) represents (in each group) a better approximation than the HF (high frequency) ones. In particular, among those dependent on ℓ, Z ℓ SAI-1 performs better than Z ℓ S-HF-1b and Z ℓ A-HF-1 , while among those independent of ℓ, Z SAI-0 performs better than Z S-HF-1a and Z S-HF-0 . Although having less accuracy, the coefficients independent of ℓ have the advantage of being simpler to implement in a 3D discretization (since dependence on ℓ implies a tangential differential operator).
(3) The most pertinent conclusion of our comparisons is that the new S-HF family performs better than the existing one A-HF in [6]. In particular, among those dependent on ℓ, Z ℓ S-HF-1b performs better than Z ℓ A-HF-1 , while among those independent of ℓ, Z S-HF-1a and Z S-HF-0 give better performance than Z A-HF-0 . This confirms that k should be the correct wavenumber to work with, instead of / ∞ . (4) Within the coefficients that are independent of ℓ, the old coefficient Z SAI-0 and the two new ones Z S-HF-1a and Z S-HF-0 give better performance than Z A-RBC-1 and Z A-HF-0 from [6]. The fact the Z S-HF-0 performs better than Z A-RBC-1 shows it is important to distinguish between the term 13 1/ coming from the first-order radial derivative of the spherical Laplacian (i.e. 2 −1 in (2.2)) and those in the impedance coefficient Z • . Table 1. Mean of error E (4.20) for ℓ between 0 and 180 and /(2 ) from 1 to 100 Hz, i.e. mean of the error pictured in Figure 6 and a corresponding case for = 5. 3), performs extremely poorly, as shown from our numerical results and those in [6]. This option misrepresents completely the oscillatory behavior of the exact solution, which leads to its incapability to distinguish between the incoming and outgoing one.

Conclusion
We have studied two problems associated to the equation obtained from the conjugation of the ideal atmospheric scalar wave equation by Liouville transformation. With the Atmo assumptions extended to the whole space, we constructed explicitly the Schwartz kernel of the resolvent which allows, among other results, the construction of outgoing solutions and radiation conditions. In the second task, we provide the analytic expression of the exact D-t-N map to be used as a reference in numerical experiments and constructed new RBC (Z S-HF-1a and Z S-HF-0 ), which are simpler for the same precision compared to existing ones. These RBC, together with Z SAI-0 , are particularly interesting for a discretization in 3D without spherical symmetry. Additionally, Z S-HF-0 lends itself readily to a time-domain discretization. As future work, theoretical and numerical experiments shall be extended to consider the complete Sun's model S+Atmo, as well as different models of atmosphere and general geometry.
Appendix A. Proof for asymptotic properties of fundamental kernel )︀ . (A.5) The original phase k (| |) can be replaced by˜k , ( ) without changing the asymptotes of the remainder, cf. (3.18).
Part 2. The result for ( ) Φ k is lengthier however follows the same procedure. It needs in addition to Part 1, the asymptotics of ( ) , ( ) , and | − | −1 and its derivative (cf. [7], Appendix F2), as well as that for and , which are also defined in terms of and ℳ and their derivatives. For this reason, their asymptotic are derived as was done for in Part 1. For more details, see Proposition 11 of [7].
Proof for Proposition 3.7. Due to the symmetry in and of Φ k , it suffices to show the first statement. For convenience, we work with the reduced kernel k .
For ≤ and ≤ , we then have The last equality used ≤ . On the other hand, when ≤ and > , using (2.31) and (2.32), )︀ .
With determined by (A.6), we can bound | k | further by Case 2. We next consider > . Since > , we automatically have > . In addition to the asymptotic of in (A.4), we need those for ℳ given by (2.34), (A.10) We write = 1 − 2 , with To treat separately the two components in the asymptotics of ℳ ′ , we further decompose 1 into 1 = 1 + 1 , with )︀ (︀ 1 + O( −1 ) )︀ (︀ǎ + O( −1 ) )︀ ; A similar decomposition is obtained for 2 = 2 + 2 , now working with the two components in the asymptotic of ℳ. We can define a constant that only depends on , ∞ and k so that )︀ .
It remains to bound the exponentials. The one in 1 is rewritten as, And using > , we have