© The Authors 2024

1 Introduction

An accurate forward-modelling of the wave propagation is an essential step in the interpretation of solar oscillations. The modelling of acoustic waves requires appropriate boundary conditions. The surface boundary conditions are especially important at high frequencies for waves that are not trapped in the interior. These evanescent waves can propagate into the atmosphere, but are also reflected in the interior because the temperature in the chromosphere increases strongly, which leads to a signature in the time-distance diagram (see e.g. Jefferies et al. 1997; Fournier et al. 2017). The importance of the atmosphere is not limited to high-frequency waves because the eigenfrequencies are also significantly affected far below the acoustic cut-off (see e.g. Morel et al. 1994).

The most common boundary condition used in helioseismology is the free-surface boundary condition, such that the Lagrangian perturbation of the pressure vanishes at the surface. This assumption ensures the self-adjointness of the adiabatic equations of stellar oscillations (Lynden-Bell & Ostriker 1967). It is commonly used in helioseismology, but does not allow for the treatment of high-frequency waves. The boundary condition can be improved considering an isothermal atmosphere, which allows waves above the acoustic cut-off frequency to propagate into the atmosphere (see e.g. Unno et al. 1989). This type of atmosphere is often used to compute the eigenvalues of solar oscillations (Christensen-Dalsgaard 2008; Townsend & Teitler 2013) and has been extensively studied (Fournier et al. 2017; Barucq et al. 2018; Barucq et al. 2020) for the simplified scalar problem, $$-\frac{{\sigma }^2}{\rho c^2}\psi -\nabla \cdot \left(\frac{1}{\rho }\nabla \psi \right)=s,$$(1)

where ψ = ρc²∇ ⋅ ξ is proportional to the divergence of the wave displacement ξ. In this equation, ρ is the density, c is the sound speed, and σ² = ω² + 2iωγ is the complex frequency. This equation is obtained from the equations of stellar oscillations (Lynden-Bell & Ostriker 1967) without gravity and background flows. Barucq et al. (2018) derived boundary conditions that are relatively easy to implement and often offer sufficient accuracy for applications in time-distance helioseismology (Fournier et al. 2017).

However, all these boundary conditions are tied to a particular model of the atmosphere, which needs to be simple enough for the boundary conditions to be derived analytically. In this paper, we overcome this restriction by using learned infinite elements (LIE) (Hohage et al. 2021), which emulate a boundary condition for any type of radially symmetric profile in the atmosphere. They allow us to truncate the computational domain at any height (above which we are not interested in the solution, for example the measurement height) even when the background parameters still vary. Thus, any 1D model of the atmosphere can be used in addition to a background solar model without significantly increasing the computational cost.

The traditional approaches to derive the boundary conditions for Eq. (1) are presented in Sect. 2, and the new approach using LIE is explained in Sect. 3. The advantages of this approach are numerically illustrated in Sect. 4, and some possible extensions of this work are sketched in Sect. 5.

2 Scalar wave equation and transparent boundary conditions

2.1 Capability of the scalar equation

We build on the computational framework for local helioseismology from Gizon et al. (2017), in which the gravity term was neglected in the equations of solar oscillations (see also Jensen & Pijpers 2003, for a derivation). Under this assumption, the crosscovariance and the power spectrum (for high-degree modes) are close to those in the observations (see Figs. 9 and 13 in Gizon et al. 2017). This model is highly simplified, but still captures most of the physics of acoustic modes and presents a realistic laboratory for solving the numerical challenges that occur in helioseismology. The resolution is computationally expensive for 2D (axisymmetric) problems and becomes highly challenging for 3D problems because of the strong stratification and thus small wavelengths near the surface. Moreover, this costly resolution has to be repeated many times to solve nonlinear inverse problems in order to image strong perturbations (see e.g. Müller et al. 2024, for an application in this setup). The approach presented in this paper reduces the computational burden by truncating the computational domain as early as possible.

However, it is not only a toy model, and this setup has already proven to be useful for helioseismic applications. We present two examples. The detection of active regions on the far side of the Sun has been an active field of research for decades because of its application for space weather. It is routinely performed using helioseismic holography, which was initially proposed by Lindsey & Braun (1997). The computation of the holograms requires a Green function, which usually only takes the background sound speed into account. The Green function of Eq. (1) has led to a significant improvement of our capability to detect active regions on the far side of the Sun (Yang et al. 2023a,b). The imaging of flows in the solar interior is another main goal of helioseismology because they are important for the solar dynamo. Differential rotation has been successfully inverted using global helioseismology (see e.g. Schou et al. 1998), but the eigenfrequencies have only a limited (second-order) sensitivity to meridional circulation. This large-scale flow can be imaged using time-distance helioseismology (Giles 1999), and the observed travel times from the Michelson Doppler Imager (MDI), the Helioseismic and Magnetic Imager (HMI), and the Global Oscillation Network Group (GONG) were inverted by Gizon et al. (2020) using the forward model given by Eq. (1). The structural errors that result from neglecting gravity only affect the mean travel times between two points and not the difference travel times that are used to image flows.

While this equation is relevant for important applications, we wish to emphasise that it leads to a shift in the computed eigenfrequencies, which are the main input of global helioseismology, and it should not be used to infer structural changes. For example, Gizon et al. (2017) modified the eigenvalue solver ADIPLS (Christensen-Dalsgaard 2011) in order to neglect gravity and obtained a frequency of 3.222 mHz without gravity, compared to 3.216 mHz for the full equations when ℓ = 85 and n = 8. This difference is far larger than the measurement noise (of the order of 0.1 µHz, Korzennik et al. 2013). Whereas the discrepancy reduces as ℓ increases (see e.g. Fig. 11 in Pham et al. 2024, for ℓ = 400), care should be taken when this simplified setup is used for low-degree modes. The adaptation of the method presented in this paper to take into account gravity is the subject of future works and is discussed in Sect. 5.

Fig. 1 Background medium considered in this study. In the solar interior (up to ra), it is characterised by its density, sound speed, attenuation, and sources that can depend on the 3D vector r. In the atmosphere (ra ≤ r ≤ rt), the background medium is radially symmetric and does not contain sources. The boundary condition is usually applied at rt, but can be directly applied at ra using learned infinite elements.

2.2 Main numerical steps for solving the wave equation

We would like to solve Eq. (1) in the geometry sketched in Fig. 1. In the interior (r < r_a), all the coefficients of Eq. (1) can depend on the 3D vector r. On the other hand, the atmosphere (r_a ≤ r ≤ r_t) is assumed to be source-free and spherically symmetric (no background flow). This greatly simplifies the solution of Eq. (1) in this region. We introduce the ansatz $$\psi =\{\begin{array}{ll}{\psi }^{\text{i}}\hfill & \text{ for }r\le r_a,\hfill \\ {\psi }^{\text{e}}\hfill & \text{ for }{r}_a\le r\le r_t\hfill \end{array}$$(2)

and decompose the solution in the exterior into spherical harmonics, $${\psi }^{\text{e}}(r,\theta ,\varphi )={\displaystyle \sum _{\ell ,m}{\psi }_{\ell m}^{\text{i}}}\left(r_a\right){\psi }_{\ell m}^{\text{e}}(r)Y_{\ell }^m(\theta ,\varphi ).$$(3)

Here, ${\psi }_{\ell m}^{\text{i}}\left(r_a\right)$ are the coefficients in the spherical harmonics expansion of ψⁱ at r = r_a, that is, $${{\psi }^{\text{i}}|}_{r=r_a}={\displaystyle \sum _{\ell ,m}{\psi }_{\ell m}^{\text{i}}}\left(r_a\right)Y_{\ell }^m(\theta ,\varphi ).$$(4)

The continuity of ψ decomposed as in Eq. (2) is ensured at r = r_a when we require $${\psi }_{\ell m}^{\text{e}}\left(r_a\right)=1.$$(5)

By inserting ansatz (3) for ψ^e into Eq. (1), because ρ and c are radially symmetric, we obtain the ordinary differential equation (ODE) $$-\left(\frac{{\sigma }^2}{c^2}-k_h^2\right){\psi }_{\ell m}^{\text{e}}-\frac{\rho }{r^2}\frac{ \text{d}}{ \text{d}r}\left(\frac{r^2}{\rho }\frac{ \text{d}{\psi }_{\ell m}^{\text{e}}}{ \text{d}r}\right)=0,r_a<r<r_t.$$(6)

We used here that the source s vanishes in the atmosphere and introduced the horizontal wavenumber $k_h:=\sqrt{\ell (\ell +1)}/r$. It is only possible to derive a similar equation for ψⁱ in the interior when the assumption of spherical symmetry is extended to the entire Sun, which is usually not possible because the Sun is heterogeneous.

Before we provide further details of how Eq. (6) is solved, we discuss the problem of matching the interior solution ψⁱ and exterior solution ψ^e in a meaningful way. To do this, we require that the normal derivatives ∂_n (radial derivative in the case of a spherical star) match at the interface, that is, $$-{\partial }_n{\psi }^{\text{i}}=-{\partial }_n{\psi }^{\text{e}}\text{ at }r=r_a\text{. }$$(7)

In view of Eq. (3), this means that $$-{\left({\partial }_n{\psi }^{\text{i}}\right)|}_{r=r_a}={\displaystyle \sum _{\ell ,m}}{\mathit{\text{𝒹t𝓃}}}_{\ell }{\psi }_{\ell m}^{\text{i}}\left(r_a\right)Y_{\ell }^m=:\mathit{\text{𝒟t𝒩}}{{\psi }^{\text{i}}|}_{r=r_a},$$(8)

where 𝒟t𝒩 is called the Dirichlet-to-Neumann operator, and $${\mathit{\text{𝒹t𝓃}}}_{\ell }:=-\frac{\text{d}{\psi }_{\ell m}^{\text{e}}}{ \text{d}r}\left(r_a\right)$$(9)

are its coefficients on the basis of spherical harmonics. Eq. (6) does not depend on m, which implies that dtn is a function of ℓ only.

To summarize, in order to solve Eq. (1) in the geometry of Fig. 1, we proceeded as follows:

• We solve Eq. (6) with the boundary condition given in Eq. (5) to determine 𝒟t𝒩 , that is, the complex numbers 𝒹t𝓃_ℓ. This task is discussed in Sect. 2.3 in more detail.

• Then, we solve Eq. (1) for ψ = ψⁱ with the boundary condition given by Eq. (8). We discuss this task in Sect. 2.4.

2.3 Atmospheric model

To solve Eq. (6), we first discuss the atmospheric model. A commonly used model is an isothermal atmosphere with a constant sound speed and an exponentially decreasing density. We refer to this model as Atmo (Fournier et al. 2017). However, this model does not take the strong increase in the temperature in the transition region into account, which results in an increase in the sound speed by one order of magnitude. We therefore used the empirical VAL-C model of the atmosphere (Vernazza et al. 1981). For the interior, we used the standard solar model S (Christensen-Dalsgaard et al. 1996). A smooth transition between the two models was used between 400 and 600 km above the surface (see Fig. 2 in Fournier et al. 2017). We chose this model to validate our approach, but the method presented in this paper is general and can be applied to any radially symmetric atmospheric model.

We first discuss the Atmo model in more detail. In this case, the atmosphere is not bounded, that is, r_t = ∞. However, the computational domain needs to be truncated for this problem to be solved numerically. To derive outgoing solutions, it is useful to recast Eq. (6) into a Schrödinger-type equation using the change in the variable ${\varphi }_{lm}^e=r{\rho }^{-1/2}{\psi }_{lm}^e$. Equation (6) becomes $$\frac{{\text{d}}^2{\varphi }_{lm}^e}{ \text{d}{r}^2}+\left(\frac{{\sigma }^2-{\omega }_c^2}{c^2}-k_h^2\right){\varphi }_{lm}^e=0,$$(10)

where ω_c is the cut-off frequency, defined as $$\frac{{\omega }_c^2}{c^2}=\frac{1}{4H^2}\left(1-2\frac{ \text{d}H}{ \text{d}r}+\frac{4H}{r}\right),$$(11)

and H is the density scale height, $$H=-{\left(\frac{\text{d\hspace{0.17em}}\ln \rho }{ \text{d}r}\right)}^{-1}.$$(12)

The cut-off frequency is an important quantity because waves with frequencies above ω_c propagate into the atmosphere, while those below this frequency are mostly reflected into the Sun. The study of waves above the cut-off frequency requires an appropriate treatment of the boundary conditions, but we show below that lower frequencies are also affected by the choice of atmosphere.

In the special case of the Atmo model, an analytic expression of the dtn was obtained by Barucq et al. (2019) in terms of the Whittaker functions W, $${\mathit{\text{𝒹t𝓃}}}_{\ell }^{\mathtt{\text{Atmo }}}=\left(\frac{1}{r_a}+\frac{1}{2H}\right)+2i\text{k}\frac{{W^{\prime }}_{-i/(2kH),\ell +\frac{1}{2}}\left(-2i\text{k}{r}_a\right)}{W_{-i/(2kH),\ell +\frac{1}{2}}\left(-2ik{r}_a\right)},$$(13)

where the wavenumber k is $${\text{k}}^2=\frac{{\sigma }^2}{c^2}-\frac{1}{4H^2}.$$(14)

The Whittaker functions are defined as the outgoing solutions of $$\frac{{\text{d}}^2{W}_{\kappa ,\mu }}{\text{d}{r}^2}+\left(-\frac{1}{4}+\frac{\kappa }{r}+\frac{1/4-{\mu }^2}{r^2}\right)W_{\kappa ,\mu }=0$$(15)

and can be evaluated with libraries such as arb (Johansson 2017).

It is not possible to provide an explicit formula for ${\mathit{\text{𝒹t𝓃}}}_{\ell }^{\mathtt{\text{VAL−C}}}$ corresponding to the VAL-C model. In this case, we have to solve Eq. (6) numerically on the interval [r_a, r_t]. As the data for the VAL-C model stop at about 2.5Mm above the surface, the right end r_t of the interval is finite. We imposed a homogeneous Dirichlet boundary condition at r_t to complement the boundary condition (5) at r_a . Due to the strong increase in the temperature in the transition region, waves do not reach the end of the domain, and the choice of boundary condition applied after the end of the VAL-C model does not cause significant differences in the obtained solutions. Equation (6) is a two-point boundary value problem that can be solved in various ways (see e.g. Keller (1992)). We used the finite-element method (FEM) since it is already employed to discretise the solar interior and allows for a convenient evaluation of the derivative of ψ_ℓm at the r_a required to obtain ${\mathit{\text{𝒹t𝓃}}}_{\ell }^{\mathtt{VAL-C}}$ via Eq. (9).

2.4 Implementing the 𝒟t𝒩 condition

We now discuss the implementation of the boundary condition given in Eq. (8) for Eq. (1). We briefly discuss the special case of a spherically symmetric background before we describe the general setting.

2.4.1 Spherically symmetric background

In a spherically symmetric background, we can also decompose the interior solution ψⁱ into spherical harmonics analogously to Eq. (3) and derive the following ODE for its coefficients: $$-\left(\frac{{\sigma }^2}{c^2}-k_h^2\right){\psi }_{\ell m}^{\text{i}}-\frac{\rho }{r^2}\frac{ \text{d}}{ \text{d}r}\left(\frac{r^2}{\rho }\frac{ \text{d}{\psi }_{\ell m}^{\text{i}}}{ \text{d}r}\right)=\rho s_{lm},0\le r\le r_a,$$(16)

where s_ℓm are the coefficients of the source term s in spherical harmonics. This allows an exact implementation of the boundary condition in Eq. (8) of the form $$-\frac{\text{d}{\psi }_{\ell m}^{\text{i}}}{ \text{d}r}\left(r_a\right)={\mathit{\text{𝒹t𝓃}}}_{\ell }{\psi }_{\ell m}^{\text{i}}\left(r_a\right).$$(17)

2.4.2 General background

When the background is not spherically symmetric, for example because of background flows, an exact implementation of Eq. (8) is no longer feasible. We have to replace 𝒟t𝒩 in Eq. (8) by a simpler expression that approximates its action on the spherical harmonics as well as possible, that is, $\mathit{\text{𝒟t𝒩}}{Y}_{\ell }^m={\mathit{\text{𝒹t𝓃}}}_{\ell }{Y}_{\ell }^m$. To do this, a fairly popular choice is $${\mathit{\text{𝒟t𝒩}}}^0{\psi }^{\text{i}}=\alpha {\psi }^{\text{i}}-\beta {\Delta }_{\Gamma }{\psi }^{\text{i}},$$(18)

for some α,β ∈ ℂ, where −Δ_Γ is the horizontal Laplacian on the boundary Γ := {r, ||r|| = r_a}, fulfilling $$-{\Delta }_{\Gamma }{Y}_{\ell }^m={\lambda }_{\ell }{Y}_{\ell }^m$$(19)

for ${\lambda }_{\ell }=\ell (\ell +1)/r_a^2$. This expression is easy to implement using the FEM. Upon multiplying Eq. (1) by a test function ϕⁱ and integrating by parts, the normal derivative −∂_nψⁱ on Γ appears naturally, $${\displaystyle \underset{\left\{r<r_a\right\}}{\int }\left(\frac{1}{\rho }\nabla {\psi }^{\text{i}}\nabla {\varphi }^{\text{i}}-\frac{{\sigma }^2}{\rho c^2}{\psi }^{\text{i}}{\varphi }^{\text{i}}\right)}\text{d}V-{\displaystyle \underset{\Gamma }{\int }\frac{1}{\rho }}{\partial }_n{\psi }^{\text{i}}{\varphi }^{\text{i}} \text{d}S={\displaystyle \int s}{\varphi }^{\text{i}}dV.$$(20)

It can then be replaced using the approximation −∂_nψⁱ ≈ 𝒟t𝒩⁰ψⁱ. This amounts to replacing the second integral in Eq. (20) by $${\displaystyle \underset{\Gamma }{\int }\frac{1}{\rho }}\left(\alpha \psi \varphi +\beta {\nabla }_{\Gamma }\psi {\nabla }_{\Gamma }\varphi \right)\text{d}S.$$(21)

We note that −Δ_Γ was integrated by parts, leading to the surface gradients Δ_Γ = (I − n ⋅ n^T)∇.

The constants α and β should be determined such that the affine function ${\lambda }_{\ell }\mapsto {\mathit{\text{𝒹t𝓃}}}_{\ell }^0$ in $${\mathit{\text{𝒟t𝒩}}}^0{Y}_{\ell }^m={\mathit{\text{𝒹t𝓃}}}_{\ell }^0{Y}_{\ell }^m,{\mathit{\text{𝒹t𝓃}}}_{\ell }^0=\alpha +\beta {\lambda }_{\ell }$$(22)

fits the exact function 𝒹t𝓃_ℓ in Eq. (9) optimally. The two choices for the Atmo model derived in Barucq et al. (2020) are $${\mathit{\text{𝒹t𝓃}}}^{\text{S}-\text{HF}-1\text{a}}=\left[\left(\frac{1}{r_a}+\frac{1}{2H}\right)-\text{ik}+\frac{\text{i}}{\text{k}H{r}_a}\right],$$(23) $$\begin{array}{cc}{\mathit{\text{𝒹t𝓃}}}^{\text{SAI}-1}\psi =& \left[\left(\frac{1}{r_a}+\frac{1}{2H}\right)-\text{ik}{\left(1-\frac{1}{{\text{k}}^{2}H{r}_a}\right)}^{1/2}\right]\hfill \\ \hfill & +\frac{\text{i}}{2\text{k}}{\left(1-\frac{1}{{\text{k}}^{2}H{r}_a}\right)}^{-1/2}{\lambda }_{\ell }.\hfill \end{array}$$(24)

A representation of the exact 𝒹t𝓃 and these two approximations is given in Fig. 2 for three typical frequencies: below, around, and above the cut-off frequency. The S-HF-1a boundary condition is of order 0 and does not depend on the harmonic degree (β = 0), while SAI-1 is an affine approximation of the 𝒹t𝓃 and scales as ℓ(ℓ + 1). For low harmonic degrees ℓ < both boundary conditions yield good approximations to the 𝒹t𝓃, but for higher values of ℓ, an affine approximation is necessary. While for frequencies below and above the cut-off frequency the affine model is a very good approximation of the exact 𝒹t𝓃, this is not the case around the cut-off frequency for high values of ℓ. However, the affine approximation is very good for ℓ ≤ 300, which is the range often considered in helioseismology to study large-scale flows (see e.g. Gizon et al. 2020, for meridional flow).

Finally, we would like to point out that the derivation of the conditions given in Eqs. (23) and (24) is based on the assumption that the sound speed and density scale height are constant in the atmosphere, which is only valid in the isothermal case. The derivation of analytical expressions for the coefficients α and β in Eq. (18) for other atmospheric models such as VAL-C is generally not possible. We would therefore need to resort to meshing the entire atmosphere for these models, that is, place the boundary Γ at r = r_t instead of r = r_a. This results in a dramatic increase in the computational cost and memory requirements as the wavelength is small in the atmosphere because the sound speed is low. An infinite atmosphere (r_t = ∞) cannot be treated using the FEM, which requires the computational domain to be finite.

Fig. 2 Comparison of 𝒹t𝓃 approximations provided by different transparent boundary conditions at 3.0, 5.2, and 8.0 mHz for the Atmo model at rt = 510 km. The reference is given by Eq. (13), and the boundary conditions S-HF-1a and SAI-1 from Barucq et al. (2020) are recalled in Eqs. (23) and (24). The LIE boundary conditions for N = 0 and N = 1 are obtained by solving the minimisation problem defined in Eq. (29).

3 Learned infinite elements

The method of LIE was introduced by Hohage et al. (2021) to allow an accurate and efficient approximation of the exact 𝒟t𝒩 boundary condition given in Eq. (8). It generalises the approach introduced in Sect. 2.4.2 in two aspects. Firstly, it allows us to obtain optimal parameters α and β in Eq. (22) for any radially symmetric model of the atmosphere, as we describe in Sect. 3.1. Secondly, it permits us to efficiently implement rational approximations of 𝒹t𝓃, that is, we can replace Eq. (22) by $${\mathit{\text{𝒟t𝒩}}}^N{Y}_{\ell }^m={\mathit{\text{𝒹t𝓃}}}_{\ell }^N{Y}_{\ell }^m,$$(25)

where ${\mathit{\text{𝒹t𝓃}}}_{\ell }^N$ is a rational function of λ_ℓ with N simple poles. The superior approximation qualities of rational functions allow us to obtain a very accurate fit of 𝒹t𝓃 even when its behaviour as a function of λ_ℓ is fairly complicated. We explain this advantage in detail in Sect. 3.2.

3.1 Affine approximation

We recall that the objective is to find parameters α and β that minimise the deviation $\left|{\mathit{\text{𝒹t𝓃}}}_{\ell }^0-{\mathit{\text{𝒹t𝓃}}}_{\ell }\right|=\left|\alpha +\beta {\lambda }_{\ell }-{\mathit{\text{𝒹t𝓃}}}_{\ell }\right|$ across a range of 0 ≤ ℓ ≤ L for some L ≥ 1. To do this, it seems natural to set up a least-squares fit. Moreover, it is common practice in helioseismology to apply filters that localize in ℓ over a certain range of primary interest. We can take this into account here by additionally multiplying the deviation with some non-negative weights w_ℓ. Hence, we wish to find α,β ∈ ℂ by minimising $${\Vert \text{Λ}\left(\begin{array}{c}\alpha \\ \beta \end{array}\right)-b\Vert }_2^2,\text{\hspace{1em}for\hspace{0.17em}Λ}=\left(\begin{array}{cc}{w}_0& w_0{\lambda }_0\\ ⋮& ⋮\\ w_L& w_L{\lambda }_L\end{array}\right),b=\left(\begin{array}{c}{w}_0{\mathit{\text{𝒹t𝓃}}}_0\\ ⋮\\ w_L{\mathit{\text{𝒹t𝓃}}}_L\end{array}\right),$$(26)

where || ⋅ ||₂ denotes the Euclidean norm. The minimum norm solution is given by $$\left(\begin{array}{c}\alpha \\ \beta \end{array}\right)={\text{Λ}}^{†}b,$$(27)

where Λ^† := (Λ*Λ)⁻¹Λ* denotes the pseudo-inverse of Λ. The explicit expression of the coefficients α and β as a function of w_ℓ and 𝒹t𝓃_ℓ is given in Appendix B.

3.2 Rational approximation

It is straightforward to generalise the least-squares approach introduced in the previous section to the rational approximation of 𝒹t𝓃 announced in Eq. (25). A rational function of λ_ℓ with N simple poles is given by $${\mathit{\text{𝒹t𝓃}}}_{\ell }^N=A_{00}+\lambda B_{00}-{\displaystyle {\displaystyle \sum }_{j=1}^N}\frac{{\left(A_{0j}+\lambda B_{0j}\right)}^2}{A_{jj}+\lambda }$$(28)

provided A_0j ≠ B_0jA_jj at j = 1,…, N holds. This ansatz is justified by the result established in Preuß (2021, Proposition 3.13), which ensures that complex numbers A_ij and B_ij exist such that the deviation $\left|{\mathit{\text{𝒹t𝓃}}}_{\ell }^N-{\mathit{\text{𝒹t𝓃}}}_{\ell }\right|$ for ℓ ≤ L decreases exponentially fast as N increases. In analogy with Sect. 3.1, we obtain A_ij and B_ij in practice by minimising the least-squares functional $$\frac{1}{2}{\displaystyle \sum _{\ell =0}^L{w}_{\ell }^2}{\left|{\mathit{\text{𝒹t𝓃}}}_{\ell }^N-{\mathit{\text{𝒹t𝓃}}}_{\ell }\right|}^2.$$(29)

We recognise that the problem (29) reduces to the problem (26) for N = 0, and the solution (A₀₀, B₀₀) = (α, β) is given by Eq. (27). For N > 0, on the other hand, problem (29) becomes non-linear and no longer allows an explicit solution. The Lev-enberg–Marquardt algorithm (see e.g. Nocedal & Wright 2006) is a popular option for solving problems like this iteratively. In practice, it is observed to perform well for minimising Eq. (29). For example, Fig. 2 compares the affine (N = 0) and first-order rational approximation (N = 1 ) to the established approximations for the Atmo model. The affine approximation (N = 0) obtained by solving the least-squares problem can at least reproduce the accuracy of the SAI-1 condition. To achieve a better fit of 𝒹t𝓃, the rational approximation (N = 1) can be used. The improvement in accuracy is clearly visible at 3.0 and 5.2 mHz.

The great advantage of LIE is that it is applicable to any spherically symmetric exterior model (by changing the data vector b in (26)), whereas the radiation boundary conditions from Barucq et al. (2018) are tied to the Atmo model. A representation of the 𝒹t𝓃 operator for the VAL-C model applied directly at the solar surface (r_a = R_⊙) is shown in Fig. 3. Its behaviour is more complicated than that for the isothermal atmosphere (Fig. 2) and it cannot be approximated by a low-order boundary condition. The LIE with N = 0 only give a crude approximation of the 𝒹t𝓃, and at least N = 1 is necessary to represent the variations with ℓ.

Fig. 3 Exact 𝒹t𝓃 for the VAL-C atmospheric model and its approximations using LIE for N = 0 and 1. The boundary condition is applied directly at the solar surface ra = R⊙.

3.3 Infinite elements

We call the method based on the rational 𝒟t𝒩^N approximation in Eq. (25) with 𝒹t𝓃^N obtained from solving the minimisation problem in Eq. (29) the learned infinite elements. The term ‘learned’ describes the process of acquiring the coefficients A_ij and B_ij of the rational approximation from the given data describing the atmosphere, that is, the numbers 𝒹t𝓃_ℓ determined by Eq. (6). The term ‘infinite elements’ refers to the algebraic structure that arises when Eq. (25) is implemented in the context of the FEM. The main idea is illustrated in Fig. 4. For N = 0, we added the surface integral in Eq. (21) involving the degrees of freedom (DOFs) on r weighted by the coefficients A₀₀ and B₀₀. For N = 1, we duplicated the DOFs of Γ and added them in the form of an additional layer to the existing discretisation. The corresponding matrix entries were determined from surface integrals like in Eq. (21), weighted with the coefficients A_ij and B_ij. This approach is computationally efficient as it preserves the sparsity structure of the finite elements, and only very few additional such layers (typically N < 5, and often N = 1) are required to achieve high accuracy. We refer to Appendix A for a detailed description of the implementation in the context of FEM and to Sect. 4.5 for a discussion of computational costs.

Fig. 4 Schematic illustration of LIE. The domain is meshed in the solar interior (up to Γ), and the degrees of freedom are defined at the boundary Γ as for classical finite elements. For N > 0, additional degrees of freedom are located in the exterior (see more details in Appendix A).

4 Numerical experiments

In order to test the accuracy of the proposed boundary conditions in terms of physical quantities relevant in helioseismology, we compare power spectrum, cross-covariance, and travel times. We first explain how the Green function is computed before we show its relation to the power spectrum using the setup from Gizon et al. (2017).

4.1 Green’s function

The reference Green function is obtained as described in Sect. 2.4.1, where the source is a delta function δ(r − r₀) located at the surface r₀ = R_⊙. The attenuation model is based on the observed line widths and follows the power law from Gizon et al. (2017). To obtain the exact boundary condition, we used ${\mathit{\text{𝒹t𝓃}}}_{\ell }={\mathit{\text{𝒹t𝓃}}}_{\ell }^{\mathtt{\text{Atmo }}}$ (Eq. (17)) for the Atmo model and ${\mathit{\text{𝒹t𝓃}}}_{\ell }={\mathit{\text{𝒹t𝓃}}}_{\ell }^{\mathtt{\text{VAL−C}}}$ for the VAL-C model. The different transparent boundary conditions were realised by replacing the exact 𝒹t𝓃_ℓ numbers in Eq. (17) by their respective approximations, that is, ${\mathit{\text{𝒹t𝓃}}}_{\ell }^{\text{S}-\text{HF}-1\text{a}},{\mathit{\text{𝒹t𝓃}}}_{\ell }^{\text{SAI}-1}$ or ${\mathit{\text{𝒹t𝓃}}}_{\ell }^N$. For the approximations of the Atmo model, we applied the boundary condition at r_a = R_⊙ + 510 km as the radiation boundary conditions were derived under the hypothesis of constant sound speed and density scale height. On the other hand, LIE can be applied at any height above which the medium is radially symmetric and without sources. Thus, we applied the boundary condition directly at the solar surface (r_a = R_⊙) for the VAL-C atmosphere, which drastically reduced the computational cost (see Sect. 4.5).

4.2 Power spectrum

When the sources in the Sun are equipartitioned, Gizon et al. (2017) showed that the power spectrum is related to the imaginary part of the Green function, $${\mathcal{P}}_{\ell }(\omega )=\frac{\mathrm{\Pi }(\omega )}{2\omega }\text{Im}\left[G_{\ell }\left(r_0,r_0,\omega \right)\right]{\mathcal{F}}_{\ell },$$(30)

where r₀ is the observation height. The function Π(ω) controls the frequency distribution of the sources and was chosen as $$\mathrm{\Pi }(\omega )={\left(1+{\left(\frac{\omega -{\omega }_0}{{\sigma }_0}\right)}^2\right)}^{-1}\text{ with }\frac{{\omega }_0}{2\pi }=3.3\text{\hspace{0.17em}mHz},\frac{{\sigma }_0}{2\pi }=0.6\text{mHz,}$$(31)

so that the distribution of the power spectrum with frequency (∑_ℓ 𝒫_ℓ(ω)) matched the observations reasonably well. Similarly, ℱ_ℓ was chosen so that the distribution of the power spectrum with the harmonic degree (∑_ω 𝒫_ℓ(ω)) agreed with the HMI power spectrum. We fit a polynomial of order 4, $${ℱ}_{\ell }={\displaystyle \sum _{i=0}^4{f}_i}{\left(\ell /L_{\text{max}}\right)}^i,$$(32)

where L_max = 1000 is the maximum harmonic degree used in this study, and the coefficients f_i were ƒ₀ = 0.97, f₁ = −1.65, f₂ = −1.14, ƒ₃ = 3.86, and ƒ₄ = −2.02.

A representation of the power spectrum is given in Fig. 5 with the exact boundary condition for the Atmo model. As already emphasised in Gizon et al. (2017), the hypothesis of source equipartition leads to a power spectrum that compares well with the observations for frequencies below the cut-off. In particular, the high-power peaks are close to the observed eigen-frequencies, with only a small shift due to the surface effect (Rosenthal et al. 1999). For frequencies below the cut-off frequency, the power spectra obtained with the different boundary conditions are similar, even though some differences are already visible by eye for Dirichlet. Above the cut-off frequency, the radiation boundary conditions and LIE perform well, but the Dirichlet boundary condition is not adapted. More surprisingly at first sight, the most challenging task is approximating the boundary condition around the cut-off frequency. For ℓ > 500, differences with the exact 𝒹t𝓃 are visible for low-order radiation boundary conditions and for LIE with N = 0. High-order LIE are necessary to emulate the wave behaviour around the cutoff properly, which can be understood by the complex behaviour of the 𝒹t𝓃 as shown in Fig. 2. These frequencies are interesting as they give information about the lower atmosphere (Vorontsov et al. 1998).

Fig. 5 Comparison of the simulated power spectrum depending on the boundary condition. Left: Atmo model using the exact 𝒹t𝓃. Right: cuts at different frequencies for different transparent boundary conditions.

Fig. 6 Difference between the frequencies computed with the VAL-C atmospheric model (reference) and the Atmo model (green) and Dirichlet boundary condition (blue) for different values of the harmonic degree ℓ.

4.3 importance of the atmospheric model for the power spectrum and the eigenfrequencies

Before we assess the efficiency of the LIE, we recall the importance of the atmosphere for helioseismic (and asteroseismic) studies. The structure of the atmosphere is one of the components that is usually advocated to explain the ‘surface effect’ (Rosenthal et al. 1999). To reduce the differences between the observed and simulated frequencies, a common approach is to patch an atmospheric model (empirical or from an average of 3D numerical simulations) to a standard solar model in the interior (see e..g Morel et al. 1994; Ball et al. 2016). We used a free-surface boundary condition, an isothermal atmosphere, and the empirical VAL-C model to present our results. We computed the power spectrum with the different boundary conditions and measured the frequencies of the resonances by fitting Lorentzian functions to the peaks. The differences with the eigenfrequencies obtained using the VAL-C model (considered as the reference) are shown in Fig. 6. As expected, the frequency shifts are strong and far exceed the measurement error (depending on the mode but of order 0.1 µHz Korzennik et al. 2013). The order of magnitude of the differences in our simplified setup (neglecting gravity) is similar to those that we obtained using GYRE (Townsend & Teitler 2013) between an isothermal atmosphere and a free-surface boundary condition.

4.4 Accuracy of the LiE

In the previous section, we showed that the atmosphere plays a significant role in the observed spectrum. In this section, we show that the VAL-C atmosphere can be emulated with sufficient accuracy using LIE. We refer to Appendix C for a similar accuracy test for the Atmo model. We computed the Green function using LIE of different orders applied directly at the solar surface and compared them with a reference obtained by applying the exact 𝒹t𝓃 from Eq. (6). The error on the eigenfrequencies for the different approximation levels is shown in Fig. 7. The approximation of order 0 is sufficient for frequencies lower than 3 mHz, and a condition of order 1 is necessary for higher frequencies to obtain an error below the measurement error.

For applications in time-distance helioseismology, we estimated the accuracy of the boundary condition in terms of travel times. To do this, we first defined the expectation value of the cross-covariance function. This was obtained by summing the harmonic coefficients of the power spectrum $$C(\theta ,\omega )={\displaystyle \underset{\ell }{{{\displaystyle \sum }}^{\text{}}}}\frac{2\ell +1}{4\pi }{F}_{\ell }(\omega ){\mathcal{P}}_{\ell }(\omega )P_{\ell }(\cos \theta ),$$(33)

where θ is the angular distance between the two observation points that are cross-correlated. P_ℓ are the Legendre polynomials of degree ℓ, and F_ℓ corresponds to a filter to select a certain range of frequencies or type of waves, for example a phase-speed filter. We applied Gaussian frequency filters centred at different frequencies and with a standard deviation σ/2π = 0.6 mHz in order to study the error in different frequency bands. Similar filters have been applied to compute frequency-dependent travel times for meridional flow measurements (Rajaguru & Antia 2020), even though the interpretation remains difficult due to systematic effects (Chen & Zhao 2018). From the filtered cross-covariance, travel times between the reference and the approximate cross-covariances were computed using the linear formula from Gizon & Birch (2002).

Figure 8 shows the travel-time difference between the reference solution and approximate boundary conditions using LIE of different orders. For frequencies below 4 mHz, the LIE with N = 0 are sufficient to reduce the error to about 0.01 s, which is smaller than the observation noise (Liang et al. 2017). Above 4 mHz and in particular around the cut-off frequency, LIE of order 1 should be used to guarantee an error smaller than 0.01 s. Higher-order (N=2) LIE drastically reduce the error (around 10⁻⁶ s). The behaviour of the error is similar for the Atmo model (see Fig. C.1). In this case, however, the LIE with N = 0 are always sufficient to keep the error below 0.01 s.

Fig. 7 Error in the location of the peaks of the modes for different values of the harmonic degree ℓ and for different boundary conditions for the VAL-C atmospheric model.

4.5 Computational cost

The previous section showed that the accuracy provided by LIE placed at the surface r_a = R_⊙ using at most N = 2 is already sufficient to incorporate the response from the complicated VAL-C atmosphere into the simulations. In this section, we discuss the incurred computational cost and provide a comparison with a traditional approach in which the mesh is extended to include the entire atmosphere. To do this, we considered a conforming finite-element discretisation of Eq. (1) in three dimensions (3D) or in an axisymmetric setting (2.5D; see Gizon et al. (2017) for details) up to the surface r_a = R_⊙ using 10 DOFs per wavelength as reference. Table 1 shows the relative increase in the DOFs when LIE of degree N are added to this discretisation. For comparison, the last two columns list the increase when the mesh is extended to the end of model S (+0.5 Mm above the surface) and to the end of the VAL-C model (+2.5 Mm), respectively. The table clearly demonstrates that using LIE is significantly more efficient than meshing the entire atmosphere. The latter would increase the number of DOFs by over 200% in three dimensions, whereas LIE provide sufficient accuracy for an increase smaller than 14%. Moreover, the increase is already significant (around 60%) when the mesh is extended until the end of model S (+0.5 Mm), which is necessary to apply the radiation boundary conditions from Barucq et al. (2018) for the Atmo model. Thus, applying LIE with N = 1 directly at the surface, for instance, is clearly preferable as it only leads to an increase of around 7%.

5 Conclusions and outlook

Helioseismology often relies on the solution of a wave equation in the solar interior with simplified models of the atmosphere (free surface or isothermal). We showed that the choice of the atmospheric model affects the position of the eigenvalues and generates artificial travel-time shifts even for waves at frequencies far below the acoustic cut-off. However, extending the computational domain to take the atmosphere into account leads to a significant increase in the computational cost. For example, the number of degrees of freedom is multiplied by 3 for 3D computations, which renders them unfeasible. We proposed to circumvent this problem by using LIE, which emulate the behaviour of a 1D atmosphere without having to extend the computational domain. This approach is accurate and computationally efficient. It allows us to place the computational boundary directly at the observation height. Placing the boundary condition at the solar surface instead of 500 km above decreases the degrees of freedom by about 50% for 3D computations. This is a significant improvement in this computationally challenging task of simulating wave propagation in a 3D medium with a realistic solar stratification.

The high flexibility of LIE renders them a suitable candidate to be employed in inversions for the solar atmosphere. The exact 𝒹t𝓃 in the minimisation problem (Eq. (29)) could be replaced by observational constraints such as the solar power spectrum. In this case, the synthetic power spectrum generated by solving Eq. (1) would be as close as possible to the observed one. We showed that a general model of the atmosphere such as VAL-C could be emulated with only a few degrees of freedom, and thus, the inversion would have to recover only a few parameters. This makes this method very attractive.

We did not consider background flows in this paper. The method presented here can be applied when flows are confined to the solar interior (e.g. meridional circulation), but they cannot be present in the atmosphere as it breaks the separability of the horizontal and radial variables, which is an essential property in the current design of LIE. This is an important limitation as it prevents us from treating rotation. First steps have been taken in Preuß (2021, chapter 8) for an axisymmetric medium considering only an azimuthal order m = 0, but additional work in this direction is necessary.

We applied the LIE to a simplified scalar equation representative of the propagation of p-modes in the solar interior. The next step is the extension to a more general form of the equation of stellar oscillations as introduced by Lynden-Bell & Ostriker (1967). The well-posedness of this equation was proven by Halla & Hohage (2021), and recent works have studied appropriate discretisations of this problem (Chabassier & Duruflé 2018; Alemán et al. 2022; Halla et al. 2022). These efforts should be combined with a suitable transparent boundary condition. Radiation boundary conditions for an isothermal atmosphere have been proposed by Barucq et al. (2021), who transformed the vector problem into a scalar equation. The general case of a spherically symmetric atmosphere was considered in Halla (2022). Therein, the authors used the Cowling approximation to derive a scalar equation in the atmosphere and coupled it to the vector equation in the interior. This coupling paves the way for using LIE in the equations of Lynden-Bell & Ostriker (1967).

Fig. 8 Cross-covariance for different boundary conditions and resulting error in travel-time measurements. Left: cross-covariance at 6.0 mHz for the different boundary conditions and window function used to select the wave packet to compute travel times (black). Right: travel-time perturbation caused by the approximate boundary conditions for the VAL-C model across the frequency spectrum for two points separated by θ = 30°.

Data availability

The source code to reproduce this research is available at https://zenodo.org/doi/10.5281/zenodo.13271496

Acknowledgements

The authors acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID 432680300 – SFB 1456 (project C04). Additionally, J.P. acknowledges funding by EPSRC grant EP/V050400/1.

Appendix A Implementation in context of FEM

We sketch the implementation of the boundary conditions covered in this article for the FEM. We refer to Hohage et al. (2021); Preuß (2021) for more details . Let us suppose that the discretisation of the first term in Eq. (20), that is, excluding terms on Γ, leads to a linear system of the form $$\left[\begin{array}{c}{L}_{II}{L}_{I\text{Γ}}\\ L_{\text{Γ}I}{L}_{\text{ΓΓ}}^{\text{int }}\end{array}\right]\left[\begin{array}{l}{\underset{\_}{\psi }}_I\hfill \\ {\underset{\_}{\psi }}_{\text{Γ}}\hfill \end{array}\right]=\left[\begin{array}{l}{\underset{\_}{s}}_I\hfill \\ {\underset{\_}{s}}_{\text{Γ}}^{\text{int }}\hfill \end{array}\right].$$(A.1)

Here, we partitioned the DOFs as ${\underset{\_}{\psi }}_I$ belonging to the interior of the domain and a vector ${\underset{\_}{\psi }}_{\Gamma }$ of DOFs on Γ.

A.1 The ideal boundary condition

Let ψ_i for i = 1,…, n_Γ denote the basis functions of the FEM with support on Γ. Further, let M and K denote mass and stiffness matrices of the FEM, that is, $$M_{ij}={\displaystyle {{\displaystyle \int }}_{\text{Γ}}}\frac{1}{\rho }{\psi }_j{\psi }_i\text{d}S,K_{ij}={\displaystyle {{\displaystyle \int }}_{\text{Γ}}}\frac{1}{\rho }{\nabla }_{\text{Γ}}{\psi }_j{\nabla }_{\text{Γ}}{\psi }_i\text{d}S.$$(A.2)

Ideally, one would add a matrix representation M DtN^ext of the second term in Eq. (20) to (A.1), which implements the exact boundary conditions given in (8). This would lead to a linear system of the form $$\left[\begin{array}{cc}{L}_{II}& L_{I\text{Γ}}\\ L_{\text{Γ}I}& L_{\text{ΓΓ}}^{\text{int }}+M{\mathit{\text{𝒟t𝒩}}}^{\text{ext}}\end{array}\right]\left[\begin{array}{l}{\underset{\_}{\psi }}_I\hfill \\ {\underset{\_}{\psi }}_{\text{Γ}}\hfill \end{array}\right]=\left[\begin{array}{c}{\underset{\_}{s}}_I\\ {\underset{\_}{s}}_{\text{Γ}}^{\text{int }}\end{array}\right].$$(A.3)

This approach is, however, not practical as it leads to many issues. Foremost, DtN^ext will be a dense matrix which is computationally very expensive to assemble and incorporate into the otherwise sparse linear system of the FEM.

A.2 Affine approximation

According to Eq. (21), the affine approximation Ɗt𝒩⁰ of Ɗt𝒩 leads to a linear system of the form $$\left[\begin{array}{cc}{L}_{II}& L_{I\text{Γ}}\\ L_{\text{Γ}I}& L_{\text{ΓΓ}}^{\text{int}}+L_{\text{ΓΓ}}\end{array}\right]\left[\begin{array}{l}{\underset{\_}{\psi }}_I\hfill \\ {\underset{\_}{\psi }}_{\text{Γ}}\hfill \end{array}\right]=\left[\begin{array}{l}{\underset{\_}{s}}_I\hfill \\ {\underset{\_}{s}}_{\text{Γ}}^{I\text{nt}}\hfill \end{array}\right]$$(A.4)

for L_ΓΓ = αM + βK. Comparing with (A.3), this corresponds to the approximation L_ΓΓ ≈ M DtN^ext.

A.3 Rational approximation

As described in Sect. 3.3 and Fig. 4, we implement the rational approximation Ɗt𝒩^N of Ɗt𝒩 by introducing additional DOFs ${\underset{\_}{\psi }}_E$ in the exterior. The linear system to be solved is then of the form $$\left[\begin{array}{ccc}{L}_{II}& L_{I\text{Γ}}& 0\\ L_{\text{Γ}I}& L_{\text{ΓΓ}}^{\text{int}}+L_{\text{ΓΓ}}& L_{\text{Γ}E}\\ 0& L_{E\text{Γ}}& L_{EE}\end{array}\right]\left[\begin{array}{c}{\underset{\_}{\psi }}_I\\ {\underset{\_}{\psi }}_{\text{Γ}}\\ {\underset{\_}{\psi }}_E\end{array}\right]=\left[\begin{array}{c}{\underset{\_}{s}}_I\\ {\underset{\_}{s}}_{\text{Γ}}^{\text{int }}\\ 0\end{array}\right]$$(A.5)

where $$\left[\begin{array}{c}{L}_{\text{ΓΓ}}{L}_{\text{Γ}E}\\ L_{E\text{Γ}}{L}_{EE}\end{array}\right]=A\otimes M+B\otimes K.$$(A.6)

Here, ⊗ denotes the Kronecker product of matrices and $$A=\left[\begin{array}{ccccc}{A}_{00}& A_{01}& \cdots & \cdots & A_{0N}\\ A_{01}& \ddots & & 0& \\ ⋮& & \ddots & & \\ ⋮& 0& & \ddots & \\ A_{0N}& & & & A_{NN}\end{array}\right],B=\left[\begin{array}{ccccc}{B}_{00}& B_{01}& \cdots & \cdots & B_{0N}\\ B_{01}& 1& & 0& \\ ⋮& & \ddots & & \\ ⋮& 0& & \ddots & \\ B_{0N}& & & & 1\end{array}\right],$$

where the non-zero matrix entries are obtained from solving the minimisation problem given in Eq. (29). By eliminating the exterior DOFs ${\underset{\_}{\psi }}_E$, it is possible to rewrite the linear system of Eq. (A.5) in the compressed form $$\left[\begin{array}{cc}{L}_{II}& L_{I\text{Γ}}\\ L_{\text{Γ}I}& L_{\text{ΓΓ}}^{\text{int}}+M{\mathit{\text{𝒟t𝒩}}}^N\end{array}\right]\left[\begin{array}{l}{\underset{\_}{\psi }}_I^{\psi }\hfill \\ {\underset{\_}{\psi }}_{\text{Γ}}\hfill \end{array}\right]=\left[\begin{array}{c}{\underset{\_}{s}}_I\\ {\underset{\_}{s}}_{\text{Γ}}^{\text{int}}\end{array}\right]$$(A.7)

for $${\mathit{\text{𝒟t𝒩}}}^N:=M^{-1}\left(L_{\text{ΓΓ}}-L_{\text{Γ}E}{L}_{EE}^{-1}{L}_{E\text{Γ}}\right),{\mathit{\text{𝒟t𝒩}}}^N\approx {\mathit{\text{𝒟t𝒩}}}^{\text{ext}}.$$(A.8)

Comparing with the affine approximation, we see that the term L_ΓΓ approximating M DtN^ext has been enriched by the Schur complement with respect to the exterior DOFs. Note that in practice, it is preferable to solve the sparse system (A.5) rather than (A.7) which contains the dense block M DtN_N on the interface Γ.

Finally, letus explainwhy Eq. (A.8) corresponds to a rational approximation of Ɗt𝒩. In analogy with Eq. (19), let ${\underset{\_}{Y}}_{\ell }$ denote the discrete eigenfunctions of the horizontal Laplacian, which fulfill $K{\underset{\_}{Y}}_{\ell }={\underset{\_}{\lambda }}_{\ell }M{\underset{\_}{Y}}_{\ell }$ for discrete eigenvalues ${\underset{\_}{\lambda }}_{\ell }\approx {\lambda }_{\ell }$. Then it has been shown in (Hohage et al. 2021, Proposition 2.1) that $${\mathit{\text{𝒟t𝒩}}}^N{\underset{\_}{Y}}_{\ell }={\mathit{\text{𝒹t𝓃}}}_{\ell }^N{\underset{\_}{Y}}_{\ell },$$(A.9)

where ${\mathit{\text{𝒹t𝓃}}}_{\ell }^N$ is the rational function from Eq. (28). That is, DtN^N emulates the rational Ɗt𝒩^N operator of Eq. (A.9) at the discrete level.

Appendix B Optimal weights of the LIE for N = 0

The least-square problem to obtain the optimal coefficients α and β for the LIE of order 0 can be solved analytically. The pseudo-inverse Λ^† is defined as $${\text{Λ}}^{†}={\left({\text{Λ}}^T\text{Λ}\right)}^{-1}{\text{Λ}}^{T}b.$$(B.1)

The matrix Λ^TΛ can be computed and inverted analytically and we obtain $$\alpha =\frac{1}{\left|{\Lambda }^T\Lambda \right|}\left[{\displaystyle \sum _{\ell =1}^L{w}_{\ell }^2}{\lambda }_{\ell }^2{\displaystyle \sum _{\ell =1}^L{w}_{\ell }^2}{\mathit{\text{𝒹t𝓃}}}_{\ell }-{\displaystyle \sum _{\ell =1}^L{w}_{\ell }^2}{\lambda }_{\ell }{\displaystyle \sum _{\ell =1}^L{w}_{\ell }^2}{\lambda }_{\ell }{\mathit{\text{𝒹t𝓃}}}_{\ell }\right],$$(B.2) $$\beta =\frac{1}{\left|{\Lambda }^T\Lambda \right|}\left[{\displaystyle \sum _{\ell =1}^L{w}_{\ell }^2}{\lambda }_{\ell }{\displaystyle \sum _{\ell =1}^L{w}_{\ell }^2}{\mathit{\text{𝒹t𝓃}}}_{\ell }-{\displaystyle \sum _{\ell =1}^L{w}_{\ell }^2}{\displaystyle \sum _{\ell =1}^L{w}_{\ell }^2}{\lambda }_{\ell }{\mathit{\text{𝒹t𝓃}}}_{\ell }\right],$$(B.3)

where $$\left|{\Lambda }^T\Lambda \right|=\left({\displaystyle \sum _{\ell =1}^L{w}_{\ell }^2}\right)\left({\displaystyle \sum _{\ell =1}^L{w}_{\ell }^2}{\lambda }_{\ell }^2\right)-{\left({\displaystyle \sum _{\ell =1}^L{w}_{\ell }^2}{\lambda }_{\ell }\right)}^2.$$(B.4)

Appendix C׃ Accuracy of the boundary conditions for the Atmo model

Fig. C.1 Cross-covariance for different boundary conditions and resulting error in travel-time measurements. Left: cross-covariance obtained after applying a frequency filter around 6.0 mHz for different approximate boundary conditions of the Atmo model. Right: associated travel time perturbation using frequency filters centered at different frequencies for a separation distance θ = 30° between the two observation points (the perturbations for the N = 1 and N = 2 conditions are smaller than 10−6 s for frequencies below 4 mHz and higher than 6 mHz and therefore not visible in the plot).

This appendix presents the error in terms of the location of the eigenfrequencies and travel times using approximate boundary conditions for the Atmo model. Figures C.l and C.2 are the equivalent of Figs. 8 and 6 for the Atmo model instead of VAL-C.

Figure C.2 shows the error on the eigenvalues between an approximate boundary condition and the exact 𝒹t𝓃. The frequencies of the resonances are measured by fitting a Lorentzian to the power spectrum for the different boundary conditions. For all the boundary conditions the error is increasing with frequency and harmonic degree. Contrary to the VAL-C model (Fig. 6), the low-order radiation boundary condition SAI-1 or the LIE of order 0 already give errors below the measurement uncertainty.

The error in terms of travel times is shown in Fig. C.1. Travel times in different ranges of frequencies are obtained by applying a Gaussian filter with standard deviation σ/2π = 0.6 mHz. The conclusions are similar to the study of the error in the eigenvalues, that is, the SAI-1 boundary condition or the LIE of order 0 are sufficient to obtain an error smaller than 0.01 s which is below the measurement error for for example meridional circulation measurements (Liang et al. 2017).

Fig. C.2 Error in the location of the peaks of the modes for different values of the harmonic degree ℓ and for different boundary conditions for the Atmo model.

References

All Tables

All Figures

Fig. 1 Background medium considered in this study. In the solar interior (up to ra), it is characterised by its density, sound speed, attenuation, and sources that can depend on the 3D vector r. In the atmosphere (ra ≤ r ≤ rt), the background medium is radially symmetric and does not contain sources. The boundary condition is usually applied at rt, but can be directly applied at ra using learned infinite elements.

In the text

Fig. 2 Comparison of 𝒹t𝓃 approximations provided by different transparent boundary conditions at 3.0, 5.2, and 8.0 mHz for the Atmo model at rt = 510 km. The reference is given by Eq. (13), and the boundary conditions S-HF-1a and SAI-1 from Barucq et al. (2020) are recalled in Eqs. (23) and (24). The LIE boundary conditions for N = 0 and N = 1 are obtained by solving the minimisation problem defined in Eq. (29).

In the text

	Fig. 3 Exact 𝒹t𝓃 for the VAL-C atmospheric model and its approximations using LIE for N = 0 and 1. The boundary condition is applied directly at the solar surface ra = R⊙.
In the text

	Fig. 4 Schematic illustration of LIE. The domain is meshed in the solar interior (up to Γ), and the degrees of freedom are defined at the boundary Γ as for classical finite elements. For N > 0, additional degrees of freedom are located in the exterior (see more details in Appendix A).
In the text

	Fig. 5 Comparison of the simulated power spectrum depending on the boundary condition. Left: Atmo model using the exact 𝒹t𝓃. Right: cuts at different frequencies for different transparent boundary conditions.
In the text

	Fig. 6 Difference between the frequencies computed with the VAL-C atmospheric model (reference) and the Atmo model (green) and Dirichlet boundary condition (blue) for different values of the harmonic degree ℓ.
In the text

	Fig. 7 Error in the location of the peaks of the modes for different values of the harmonic degree ℓ and for different boundary conditions for the VAL-C atmospheric model.
In the text

Fig. 8 Cross-covariance for different boundary conditions and resulting error in travel-time measurements. Left: cross-covariance at 6.0 mHz for the different boundary conditions and window function used to select the wave packet to compute travel times (black). Right: travel-time perturbation caused by the approximate boundary conditions for the VAL-C model across the frequency spectrum for two points separated by θ = 30°.

In the text

Fig. C.1 Cross-covariance for different boundary conditions and resulting error in travel-time measurements. Left: cross-covariance obtained after applying a frequency filter around 6.0 mHz for different approximate boundary conditions of the Atmo model. Right: associated travel time perturbation using frequency filters centered at different frequencies for a separation distance θ = 30° between the two observation points (the perturbations for the N = 1 and N = 2 conditions are smaller than 10−6 s for frequencies below 4 mHz and higher than 6 mHz and therefore not visible in the plot).

In the text

	Fig. C.2 Error in the location of the peaks of the modes for different values of the harmonic degree ℓ and for different boundary conditions for the Atmo model.
In the text