© ESO, 2017

1. Introduction

The interpretation of local helioseismology observations requires a description of wave propagation in the Sun (the forward problem). Solar seismic waves are described by a wave operator and boundary conditions. The boundary conditions are important, in particular if one wants to study waves with frequencies above the acoustic cutoff frequency.

Common practice in helioseismology is to use a free-surface boundary condition (e.g., Christensen-Dalsgaard et al. 1996; Bogdan et al. 1996); but this approach forces all waves to reflect and thus is not appropriate for high-frequency waves. One solution to model outgoing waves is to introduce perfectly matched layers (PML), first introduced by Berenger (1994). This method consists in adding an extra layer to the computational domain in which the waves are artificially damped by some user-defined factor. Applications have included the modeling of waves in stratified media in geophysics (Komatitsch & Martin 2007), as well as solar magnetohydrodynamic wave simulations (Cameron et al. 2008; Hanasoge et al. 2010; Santamaria et al. 2015). Applications of PMLs in helioseismology have shown that, though the boundary layer absorbs the waves as expected, this is at the expense of two compromises. Firstly, the addition of a layer significantly increases the computational burden to capture the small wavelengths in the upper atmosphere. Second, the PMLs substitute a constant infinite medium which is not the case for the Sun. The purpose of this paper is to propose and implement a radiative boundary condition (RBC) that is not subject to these two limitations.

We consider a scalar wave equation for the divergence of the wave displacement ψ(r,ω) = ρ(r)c²(r)∇·ξ^ξξ(r,ω), which in the frequency domain is given by $\mathrm{-}\frac{{\mathit{\sigma }}^{\mathrm{2}}}{\mathit{\rho }{\mathit{c}}^{\mathrm{2}}}\mathit{\psi }\mathrm{-}\mathrm{\nabla }\mathrm{·}\left(\frac{\mathrm{1}}{\mathit{\rho }}\mathrm{\nabla }\mathit{\psi }\right)\mathrm{=}\mathit{s,}$(1)where ${\mathit{\sigma }}^{\mathrm{2}}\mathrm{=}{\mathit{\omega }}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathrm{i}\mathit{\omega \gamma ,}$(2)which depends on mode frequency ω and attenuation γ(ω). This equation, of Helmholtz type, was studied in detail by Gizon et al. (2017) and is appropriate to study acoustic waves that propagate in the solar convection zone. The medium is specified by the sound-speed c, and the density ρ, which are read from the standard solar model S (Christensen-Dalsgaard et al. 1996). The source function s(r,ω) excites the waves. The development of RBCs for the Helmholtz equation began with Bayliss et al. (1982), but was designed under the assumption that a homogeneous medium is situated outside of the computational domain, which is not the case for the Sun. In order to overcome this difficulty, Gizon et al. (2017) extended and smoothed model S up to 4 Mm above the photosphere in order to reach a constant density and sound speed, where a Sommerfeld boundary condition was applied. While this method has produced a solar-like power spectrum, it requires a significant increase in computational cost in order to model the additional layer in the atmosphere. Recently Barucq et al. (2017) proposed a set of radiative boundary conditions for an atmosphere in which the sound speed is constant and the density decays exponentially. In this paper, we study these proposed boundary conditions in the context of time-distance helioseismology.

In Sect. 2, we describe the equations that we are solving for one-dimensional (1D) and two-dimensional (2D) background models. In Sect. 3, we describe the radiative boundary conditions from Barucq et al. (2017). Section 4 details the numerical method used to solve Eq. (1), and Sect. 5 compares the various boundary condition properties and their absorption properties. In Sect. 6, we extend the atmosphere with a Vernazza Avrett Loeser (VAL) atmospheric model and study its influence on the time-distance diagram at high frequencies.

2. Scalar wave equations

2.1. Two-dimensional background models

Without any assumption, the domain in which Eq. (1)is defined is three dimensional (3D) and encompasses the entire Sun. It is however reasonable to assume that the solar background properties are symmetric about the rotation axis (ẑ). Under this assumption Eq. (1)can then be projected onto azimuthal modes $\mathit{\psi }\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\sum _{\mathit{m}}{\mathit{\psi }}_{\mathit{m}}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}{\mathrm{e}}^{\mathrm{i}\mathit{m\phi }}\mathit{.}$(3)Each mode ψ_m is then the solution of a 2D partial differential equation, $\mathrm{-}\frac{{\mathit{\sigma }}_{\mathit{m}}^{\mathrm{2}}}{\mathit{\rho }{\mathit{c}}^{\mathrm{2}}}{\mathit{\psi }}_{\mathit{m}}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ \mathrm{\nabla }\end{array}\mathrm{·}\left(\frac{\mathrm{1}}{\mathit{\rho }}\begin{array}{c}\mathrm{˜}\\ \mathrm{\nabla }\end{array}{\mathit{\psi }}_{\mathit{m}}\right)\mathrm{=}{\mathit{s}}_{\mathit{m}}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\mathit{,}$(4)where ${\mathit{\sigma }}_{\mathit{m}}^{\mathrm{2}}\mathrm{=}{\mathit{\omega }}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathrm{i}\mathit{\omega \gamma }\mathrm{-}{\mathit{m}}^{\mathrm{2}}{\mathit{c}}^{\mathrm{2}}\mathit{/}\mathrm{(}\mathit{r}\sin \mathit{\theta }{\mathrm{)}}^{\mathrm{2}}$ and $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathrm{\nabla }\end{array}\mathrm{=}{r̂}{\mathit{\partial }}_{\mathit{r}}\mathrm{+}{\theta ̂}\frac{\mathrm{1}}{\mathit{r}}{\mathit{\partial }}_{\mathit{\theta }}$ is the 2D spatial gradient operator and $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathrm{\nabla }\end{array}\mathrm{·}{f}\mathrm{=}{\mathit{\partial }}_{\mathit{r}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}{\mathit{f}}_{\mathit{r}}\mathrm{\right)}\mathit{/}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{\partial }}_{\mathit{\theta }}\mathrm{\left(}\sin \mathit{\theta }{\mathit{f}}_{\mathit{\theta }}\mathrm{\right)}\mathit{/}\mathit{r}\sin \mathit{\theta }$ is the 2D divergence operator. The source term s_m represents the projection of s onto the azimuthal modes. Figure 1 shows the resultant geometrical setup. Here, R_⊙ denotes the solar radius and r_t the limit of the computational domain where an artificial boundary condition needs to be applied.

Fig. 1 Geometrical setup for 2.5D computations. The computational domain is delimited by thick lines and the photosphere by the dashed line.

2.2. 1D background models

When the background model depends only on radius, the solution can be written in terms of spherical harmonic functions: $\mathit{\psi }\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\sum _{\mathit{\ell }\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\sum _{\mathit{m}\mathrm{=}\mathrm{-}\mathit{\ell }}^{\mathit{\ell }}{\mathit{\psi }}_{\mathit{\ell m}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{Y}}_{\mathit{\ell }}^{\mathit{m}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\mathit{.}$(5)For each ℓ and m the function ψ_ℓm(r) is defined on the 1D domain 0 ≤ r ≤ r_t and solves $\mathrm{-}\left(\frac{{\mathit{\sigma }}^{\mathrm{2}}}{{\mathit{c}}^{\mathrm{2}}}\mathrm{-}{\mathit{k}}_{\mathrm{h}}^{\mathrm{2}}\right){\mathit{\psi }}_{\mathit{\ell m}}\mathrm{-}\frac{\mathit{\rho }}{{\mathit{r}}^{\mathrm{2}}}\frac{\mathrm{d}}{\mathrm{d}\mathit{r}}\left(\frac{{\mathit{r}}^{\mathrm{2}}}{\mathit{\rho }}\frac{\mathrm{d}{\mathit{\psi }}_{\mathit{\ell m}}}{\mathrm{d}\mathit{r}}\right)\mathrm{=}\mathit{\rho }{\mathit{s}}_{\mathit{\ell m}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{,}$(6)where ${\mathit{k}}_{\mathrm{h}}\mathrm{=}\frac{\sqrt{\mathit{\ell }\mathrm{(}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}}}{\mathit{r}}\mathit{,}$(7)which is the local horizontal wavenumber, and s_ℓm is the projection of the source s onto spherical harmonics.

2.3. Atmospheric model

Here we discuss the background atmosphere used in this study. The sound speed and density are taken from the standard model S (Christensen-Dalsgaard et al. 1996) which extends up to 500 km above the solar surface. We extend this model by supposing that the density continues to decay exponentially at the same rate as at the end of model S (density scale height equals 125 km). Additionally, the sound speed is smoothed to a constant value of 6855 m/s. We refer to this model as Atmo. This model is similar to the one presented by Schunker et al. (2011) but with a different density scale height (125 km compared to 105 km). While Atmo is reasonable to model the lower chromosphere, it does not reproduce the temperature jump at the transition region around 2 Mm proposed in atmospheric models such as the VAL models (Vernazza et al. 1981). Plots of the sound speed and density for Atmo and VAL-C are given in Fig. 2. A smooth transition between the model S and the VAL-C atmosphere is done between 400 and 600 km.

Fig. 2 Sound speed (red) and density (blue) in the upper layers of the Sun as functions of height measured from the photosphere. Below 500 km (vertical dashed line), the sound speed and density are from model S (Christensen-Dalsgaard et al. 1996). Above 500 km, the model Atmo (solid lines) consists of an extension where density decays exponentially (H = 105 km) and sound speed is constant (c = 6855 m/s). The dashed lines show the VAL-C atmospheric model up to a height of 2500 km (Vernazza et al. 1981), which will be used in Sect. 6.

3. Radiative boundary conditions

3.1. Radiative boundary condition for 1D background

It is possible to transform Eq. (6)in order to obtain a Helmholtz equation that is more convenient to explain the behavior of high-frequency waves. This can be done by changing the unknown, ${\mathit{\phi }}_{\mathit{lm}}\mathrm{=}\mathit{r}{\mathit{\rho }}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{\psi }}_{\mathit{lm}}\mathit{,}$(8)resulting in a new form for Eq. (6), $\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{\phi }}_{\mathit{lm}}}{\mathrm{d}{\mathit{r}}^{\mathrm{2}}}\mathrm{+}\mathit{Q}{\mathit{\phi }}_{\mathit{lm}}\mathrm{=}\mathrm{-}\mathit{r}{\mathit{\rho }}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{s}}_{\mathit{lm}}\mathit{,}$(9)where $\mathit{Q}\mathrm{=}\frac{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{c}}^{\mathrm{2}}}{{\mathit{c}}^{\mathrm{2}}}\mathrm{-}{\mathit{k}}_{\mathit{h}}^{\mathrm{2}}\mathit{,}$(10)and ω_c is a cut-off frequency defined by $\frac{{\mathit{\omega }}_{\mathrm{c}}^{\mathrm{2}}}{{\mathit{c}}^{\mathrm{2}}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}{\mathit{H}}^{\mathrm{2}}}\left(\mathrm{1}\mathrm{-}\mathrm{2}\frac{\mathrm{d}\mathit{H}}{\mathrm{d}\mathit{r}}\mathrm{+}\frac{\mathrm{4}\mathit{H}}{\mathit{r}}\right)\mathit{,}$(11)where $\mathit{H}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{-}{\left(\frac{\mathrm{d}\hspace{0.17em}\ln \mathit{\rho }}{\mathrm{d}\mathit{r}}\right)}^{-1}\mathit{,}$(12)which is the density scale height.

A representation of ω_c is given in Fig. 3. With Atmo, the density scale height is constant in the new region, H = 125 km, resulting in a constant acoustic cutoff ω_t/ 2π = 5.3 mHz. The wave behavior is determined by the sign of the real part of Q. The low-frequency waves (ω<ω_t) are reflected back in the Sun, while the high-frequency ones (ω>ω_t) tunnel through the peak at 100 km below the surface and continue propagating into the upper atmosphere.

Fig. 3 Acoustic cut-off frequency ωc as a function of height from the photosphere. Below 500 km, ωc is obtained using model S and Eq. (11); the density scale height is obtained after smoothing the density. Above 500 km, the acoustic cut-off frequency is constant and equal to ωt = 5.3 mHz. Near − 200 km, ωc is complex, although ${\mathit{\omega }}_{\mathrm{c}}^{\mathrm{2}}$ is always positive and only ${\mathit{\omega }}_{\mathrm{c}}^{\mathrm{2}}$ enters the wave equation.

With Eq. (9)in hand, a radiative boundary condition is given by $\frac{\mathrm{d}{\mathit{\phi }}_{\mathit{lm}}}{\mathrm{d}\mathit{r}}\mathrm{=}\mathrm{i}{\mathit{Q}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{\phi }}_{\mathit{lm}}\mathrm{at}\mathit{r}\mathrm{=}{\mathit{r}}_{\mathrm{t}}\mathit{.}$(13)We note that Q can be negative and thus Q^{1 / 2} is imaginary for frequencies above the acoustic cut-off ω_t. Using Eqs. (8)and (13), we can write the radiation boundary condition for the initial unknown ψ_lm as $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{\psi }}_{\mathit{\ell m}}}{\mathrm{d}\mathit{r}}\mathrm{=}\left[\mathrm{i}{\mathit{Q}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{-}\left(\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{t}}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}\mathit{H}}\right)\right]{\mathit{\psi }}_{\mathit{\ell m}}\mathrm{at}\mathit{r}\mathrm{=}{\mathit{r}}_{\mathrm{t}}\mathit{.}\mathrm{AtmoNonLocal}& & \end{array}$This boundary condition (referred to herein as Atmo Non Local) depends upon both the background model (density, sound speed) and the mode ℓ and its implementation is possible for any numerical method used to solve Eq. (6). In the case of a homogeneous medium, it reduces to the Sommerfeld radiation condition where the derivative of the solution is linked to the solution itself via a wavenumber that depends on the properties of the medium.

3.2. Approximate boundary conditions for the 2D problem

The implementation of the (Atmo Non Local) boundary condition is not possible for 2D computations as the term $\mathit{\ell }\mathrm{(}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}\mathit{/}{\mathit{r}}_{\mathrm{t}}^{\mathrm{2}}$ present in Q should be replaced by the horizontal Laplacian Δ_h${\mathrm{\Delta }}_{\mathrm{h}}{\mathit{\psi }}_{\mathit{m}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}\sin \mathit{\theta }}\frac{\mathit{\partial }}{\mathit{\partial \theta }}\left(\sin \mathit{\theta }\frac{\mathit{\partial }{\mathit{\psi }}_{\mathit{m}}}{\mathit{\partial \theta }}\right)\mathrm{-}\frac{{\mathit{m}}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\theta }}{\mathit{\psi }}_{\mathit{m}}\mathit{.}$(14)However, since this term is located under the square root, it leads to a spatially non-local operator which cannot be implemented directly. In order to overcome this difficulty, an asymptotic expansion of the square-root is required. This can be achieved in one of two ways as detailed in Barucq et al. (2017).

The first approach is to take a high-frequency approximation by making a first-order approximation of the square root supposing that ω_t/ω is small, $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{\psi }}_{\mathit{m}}}{\mathit{\partial r}}\mathrm{=}\left[\mathrm{-}\left(\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{t}}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}\mathit{H}}\right)\mathrm{+}\frac{\mathrm{i}\mathit{\sigma }}{\mathit{c}}\left(\mathrm{1}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\frac{{\mathit{\omega }}_{\mathrm{t}}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}\right)\right]{\mathit{\psi }}_{\mathit{m}}\mathrm{-}\frac{\mathrm{i}\mathit{c}}{\mathrm{2}\mathit{\sigma }}{\mathrm{\Delta }}_{\mathrm{h}}{\mathit{\psi }}_{\mathit{m}}\mathit{.}\mathrm{AtmoHF}\mathrm{1}& & \end{array}$Here, the horizontal Laplacian is not under a square root and thus can be implemented, for example, using the finite element method (FEM; see following section).

The other approach is to expand around the angle of incidence with respect to the normal of the spherical surface. If we consider waves with k_h ≪ ω/c, the expansion of the square root leads to: $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{\psi }}_{\mathit{m}}}{\mathit{\partial r}}& \mathrm{=}& \left[\mathrm{-}\left(\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{t}}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}\mathit{H}}\right)\mathrm{+}\mathrm{i}{\left(\frac{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{t}}^{\mathrm{2}}}{{\mathit{c}}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\right]{\mathit{\psi }}_{\mathit{m}}\\ & & \mathrm{-}\frac{\mathrm{i}}{\mathrm{2}}{\left(\frac{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{t}}^{\mathrm{2}}}{{\mathit{c}}^{\mathrm{2}}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}{\mathrm{\Delta }}_{\mathrm{h}}{\mathit{\psi }}_{\mathit{m}}\mathit{.}\mathrm{AtmoSAI}\mathrm{1}\end{array}$The general form of this BC is similar to (Atmo HF 1) since the radial derivative at the boundary $\frac{\mathit{\partial }{\mathit{\psi }}_{\mathit{m}}}{\mathit{\partial r}}\mathrm{\left(}{\mathit{r}}_{\mathrm{t}}\mathrm{\right)}$ is proportional to the sum of two terms: one linked to the solution itself ψ_m and the other to the horizontal Laplacian Δ_hψ_m. The only difference between (Atmo HF 1) and (Atmo SAI 1) lies in the factors in front of these two terms. While the general form of these two boundary conditions is similar, we expect them to behave differently depending on the validity of the hypothesis made in their respective derivations.

In the FEM, Δ_hψ_m is straightforward to compute since it naturally arises from the weak-form (see following Section). However, other methods may encounter difficulties in computing Δ_hψ_m and thus we seek a further expansion of (Atmo Non Local) to ease computational difficulties. A 0th order expansion of (Atmo Non Local) can be done by simply neglecting the term involving ${\mathit{k}}_{\mathrm{h}}^{\mathrm{2}}$ while not expanding around 1 /σ², as was done for (Atmo HF 1). This leads to the last BC we consider in this study: $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{\psi }}_{\mathit{m}}}{\mathit{\partial r}}\mathrm{=}\left[\mathrm{-}\left(\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{t}}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}\mathit{H}}\right)\mathrm{+}\mathrm{i}{\left(\frac{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{t}}^{\mathrm{2}}}{{\mathit{c}}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\right]{\mathit{\psi }}_{\mathit{m}}\mathit{.}\mathrm{AtmoSAI}\mathrm{0}& & \end{array}$This condition is again a Sommerfeld-like boundary condition with a complex wavenumber that depends on the properties of the medium. However, unlike (Atmo Non Local), it does not take into account the horizontal wavenumber k_h. It is interesting to note that this boundary condition is similar to the one obtained by Gough (1993) by using a plane-parallel isothermal atmosphere $\frac{\mathit{\partial }{\mathit{\psi }}_{\mathit{m}}}{\mathit{\partial r}}\mathrm{\left(}{\mathit{r}}_{\mathrm{t}}\mathrm{\right)}\mathrm{=}\left[\frac{\mathrm{1}}{\mathrm{2}\mathit{H}}\mathrm{+}\mathrm{i}{\left(\frac{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{-}{{\mathit{\omega }}_{\mathrm{c}}}_{\mathrm{0}}^{\mathrm{2}}}{{\mathit{c}}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\right]{\mathit{\psi }}_{\mathit{m}}\mathit{.}$(15)The term in 1 /r_t is missing due to the plane-parallel approximation and the acoustic cut-off ω_c0 = c/ 2H is the isothermal one contrary to the one in (Atmo SAI 0).

A summary of the different boundary conditions is given in Table 1. It also includes the (Dirichlet) boundary condition which supposes that the Lagrangian perturbation of the pressure is zero at the surface and is often used in helioseismology $\begin{array}{ccc}{\mathit{\psi }}_{\mathit{m}}\mathrm{\left(}{\mathit{r}}_{\mathrm{t}}\mathrm{\right)}\mathrm{=}\mathrm{0.}\mathrm{Dirichlet}& & \end{array}$This BC is the limit of our model when the density in the atmosphere decays infinitely quickly (H → ∞). To test the boundary conditions presented thus far, we use an “infinite” atmosphere where a Dirichlet boundary condition is imposed at 35 Mm above the photosphere: $\begin{array}{ccc}{\mathit{\psi }}_{\mathit{m}}\mathrm{\left(}{\mathit{r}}_{\mathrm{t}}\mathrm{\right)}\mathrm{=}\mathrm{0} \mathrm{at}{\mathit{r}}_{\mathit{t}}\mathrm{=}{\mathit{R}}_{\mathrm{\odot }}\mathrm{+}\mathrm{35}\mathrm{Mm}\mathit{.}\mathrm{Reference}& & \end{array}$In this way, all the waves traveling in the extended atmosphere never reach the boundary at 35 Mm due to the damping. This boundary condition is denoted as (Reference) for the remainder of this study.

All the proposed boundary conditions can be implemented for a 1D solver and present the same degree of complexity so one can use (Atmo Non Local) as it represents exactly our atmospheric model. For a 2D problem, this exact boundary condition is no longer available and one should choose between a high-frequency or small angle of incidence approximation.

4. Numerical implementation

4.1. Weak formulation of the equations

In order to solve Eq. (4)for an azimuthally symmetric background or Eq. (6)for a spherically symmetric background, we use the FEM. The implementation of the FEM requires writing the weak form of Eqs. (4)or (6). In order to show how the boundary condition is implemented, we write the weak formulation for the 2D case (the derivation is similar in 1D). Let us multiply Eq. (4)by a test function φ and integrate over the whole domain $\mathrm{-}{\mathrm{\int }}_{\mathrm{\Sigma }}\frac{{\mathit{\sigma }}_{\mathit{m}}^{\mathrm{2}}}{\mathit{\rho }{\mathit{c}}^{\mathrm{2}}}{\mathit{\psi }}_{\mathit{m}}\mathit{\phi }\mathrm{d}\mathit{S}\mathrm{-}{\mathrm{\int }}_{\mathrm{\Sigma }}\begin{array}{c}\mathrm{˜}\\ \mathrm{\nabla }\end{array}\mathrm{·}\left(\frac{\mathrm{1}}{\mathit{\rho }}\begin{array}{c}\mathrm{˜}\\ \mathrm{\nabla }\end{array}{\mathit{\psi }}_{\mathit{m}}\right)\mathit{\phi }\mathrm{d}\mathit{S}\mathrm{=}{\mathrm{\int }}_{\mathrm{\Sigma }}{\mathit{s}}_{\mathit{m}}\mathit{\phi }\mathrm{d}\mathit{S}\mathit{.}$(16)The term containing the Laplacian is then integrated by parts ${\mathrm{\int }}_{\mathrm{\Sigma }}\begin{array}{c}\mathrm{˜}\\ \mathrm{\nabla }\end{array}\mathrm{·}\left(\frac{\mathrm{1}}{\mathit{\rho }}\begin{array}{c}\mathrm{˜}\\ \mathrm{\nabla }\end{array}{\mathit{\psi }}_{\mathit{m}}\right)\mathit{\phi }\mathrm{d}\mathit{S}\mathrm{=}\mathrm{-}{\mathrm{\int }}_{\mathrm{\Sigma }}\left(\frac{\mathrm{1}}{\mathit{\rho }}\begin{array}{c}\mathrm{˜}\\ \mathrm{\nabla }\end{array}{\mathit{\psi }}_{\mathit{m}}\right)\mathrm{·}\begin{array}{c}\mathrm{˜}\\ \mathrm{\nabla }\end{array}\mathit{\phi }\mathrm{d}\mathit{S}\mathrm{+}{\mathrm{\int }}_{\mathit{\partial }\mathrm{\Sigma }}\frac{\mathrm{1}}{\mathit{\rho }}{\mathit{\partial }}_{\mathit{n}}{\mathit{\psi }}_{\mathit{m}}\mathit{\phi }\mathrm{d}\mathit{l,}$(17)where the last term is a line integral over half a circle at the computational boundary (see Fig. 1). The normal derivative ∂_nψ_m is then replaced by the appropriate boundary condition. If the boundary condition contains a horizontal Laplacian as (Atmo HF 1) or (Atmo SAI 1) then the surface term is further integrated by parts.

For 1D and 2D backgrounds, the domain is decomposed into cells (segments in 1D, quadrangles in 2D) and the solution is determined as a polynomial on each cell. The degree of this polynomial determines the order of discretization of the method. The solution ψ_m and the test function φ are expressed as a sum of polynomials on each cell which leads to the matrix form that is solved numerically (see e.g., the textbook by Zienkiewicz et al. 2005). For both types of background, we use the Montjoie code, developed at Inria¹ and adapted to helioseismology by Gizon et al. (2017) to solve the scalar wave equation with solar parameters.

4.2. Computational cost

In order to demonstrate the efficiency gained by these new boundary conditions we compare the computational cost when the boundary of the domain r_t is located at a height of 500 km (end of model S) and at 4 Mm as in Gizon et al. (2017). The mesh size is chosen in order to have at least 10 degrees of freedom per wavelength (1 cell with 10th order polynomials).

For the boundary condition at 500 km, we need N_r = 33 mesh points in the radial direction. For the 2D code, the mesh is built so that the cells have similar dimensions in the horizontal and vertical directions. This results in many cells in the θ direction close to the surface (see Fig. 6 in Gizon et al. 2017) producing a 2D mesh consisting of approximately 10 000 cells (see Table 2).

Extending the computational boundary to 4 Mm results in five additional layers of cells in the radial direction (N_r = 38). These cells are of equal height since the sound speed is constant. This extension adds about 50% more cells to the mesh and results in an increase of the total used memory and computational time of 50% as shown in Table 2. Thus, being able to simulate the propagation of the waves in the extended atmosphere and keeping the computational boundary at 500 km above the surface will result in a significant saving in computational burden.

5. Comparison of the different boundary conditions

To compare the different boundary conditions, we compute the expectation value of the cross-covariance function. Gizon et al. (2017) showed that if the sources are spatially uncorrelated and equipartitioned then the expectation value of the cross-covariance is related to the imaginary part of the Green’s function $\overline{\mathit{C}}\mathrm{\left(}{r}\mathit{,}{{r}}_{\mathit{s}}\mathit{,\omega }\mathrm{\right)}\mathrm{=}\mathit{f}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{Im}\mathit{G}\mathrm{\left(}{r}\mathit{,}{{r}}_{\mathit{s}}\mathit{,\omega }\mathrm{\right)}\mathit{.}$(18)The function f(ω) controls the frequency dependence of the source covariance and will be used to filter waves in a given frequency range. Since waves of frequency below or above the acoustic cutoff frequency behave differently at the surface (where ω_t/ 2π ~ 5.3 mHz), we use the function f(ω) to select waves centered around two frequencies, $\mathit{f}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}\exp \left(\mathrm{-}\frac{\mathrm{(}\left|\mathit{\omega }\right|\mathrm{-}{\mathit{\omega }}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}{\mathit{s}}^{\mathrm{2}}}\right)\mathit{,}$(19)where ω₀/ 2π = 3 mHz or 6.5 mHz and s/ 2π = 0.3 mHz.

The Green’s function G is obtained by solving Eq. (1)with a Dirac delta function as source term located at r_s. For each boundary condition, Green’s functions are computed for 7200 frequencies uniformly distributed from 0 to 8.3 mHz, corresponding to a 10-day observational period at 60 s cadence. The waves are damped as a function of frequency according to the power law model of Gizon et al. (2017).

To quantify the reflexivity of the different boundary conditions, we fit two wavelets to the cross-covariance; one for each skip, ${\mathit{C}}_{\mathrm{fit}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{2}}^{}{\mathit{A}}_{\mathit{i}}\exp \left(\mathrm{-}\frac{\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathit{i}}^{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}}\right)\cos \mathrm{(}{\mathit{\omega }}_{\mathit{i}}\mathit{t}\mathrm{+}{\mathit{\phi }}_{\mathit{i}}\mathrm{)}\mathit{.}$(20)We measure the energy of the outgoing (first skip) and reflected (second skip) wave packets as follows: ${\mathrm{ℰ}}_{\mathit{i}}\mathrm{\left(}{r}\mathit{,}{{r}}_{\mathit{s}}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{{\mathit{w}}_{\mathit{i}}}\overline{\mathit{C}}\mathrm{(}{r}\mathit{,}{{r}}_{\mathit{s}}\mathit{,t}{\mathrm{)}}^{\mathrm{2}}\hspace{0.17em}\mathrm{d}\mathit{t,}$(21)where the intervals ${\mathit{w}}_{\mathit{i}}\mathrm{=}\mathrm{[}{\mathit{t}}_{\mathit{i}}^{\mathrm{0}}\mathrm{-}\mathrm{4}{\mathit{\sigma }}_{\mathit{i}}\mathit{,}{\mathit{t}}_{\mathit{i}}^{\mathrm{0}}\mathrm{+}\mathrm{4}{\mathit{\sigma }}_{\mathit{i}}\mathrm{]}$ select the first- and second-skip wave packets.

5.1. Exact atmospheric radiation boundary conditions

As an initial test, the Dirac source is located at the center of the Sun. In this case, ${\mathit{k}}_{\mathrm{h}}^{\mathrm{2}}\mathrm{=}\mathrm{0}$ and the RBCs (Atmo SAI 1) and (Atmo SAI 0) are equivalent to (Atmo Non Local) and thus all boundary conditions should behave similarly.

Fig. 4 Temporal cross-covariance $\mathit{c}\mathrm{\left(}\mathit{r}{\mathrm{)}}^{\mathrm{2}}\overline{\mathit{C}}\mathrm{\right(}\mathrm{0}\mathit{,r,t}\mathrm{\left)}\exp \mathrm{\right(}\mathit{t}\mathit{/}\mathrm{90}\min \mathrm{)}$ computed with (Atmo Non Local) and (Dirichlet) boundary conditions, originating from the center, with frequencies selected around 3 mHz (left panel) and 6.5 mHz (right panel). As a comparison, we also plot the cross-covariance when using the (Reference) BC (dashed cyan line). The signal is multiplied by exp(t/ 90min) to counteract the effect of the damping, and normalized by the maximum of its amplitude at t = 30 min. The amplitude is plotted as a function of the acoustic depth $\mathit{\tau }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{{}^{\mathrm{\int }}}_{\mathit{r}}^{{\mathit{r}}_{\mathit{a}}}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}\mathit{/}\mathit{c}\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}$ for better readability of near-surface behavior. The dashed line indicates the solar surface.

Figure 4 shows the spatial evolution of the cross-covariance for wave packets centered at ω₀/ 2π = 3 mHz and ω₀/ 2π = 6.5 mHz. The wave packet propagates from the right to the left of the figure before being reflected and traveling in the other direction. It is shown that for any boundary condition, the low-frequency signal is identically reflected near the solar surface. The energy ratios between the outgoing and reflected wave packets are given in Table 3 for each boundary condition. For low-frequency waves approximately 22% of the energy of the outgoing packet is reflected (with energy loss due to damping) corresponding to a reduction in amplitude of about 47%. The lack in difference between each boundary condition comes as no surprise since the low-frequency waves are reflected by the density profile before reaching the computational boundary.

At frequencies above ω_t, Dirichlet reflects the same quantity of energy as the 3 mHz case, though the cause for reflection, at this frequency, is the boundary and not the density profile. For the RBCs and (Reference) BC, only 2% comes back with the extended atmosphere. This weak reflection is not due to the RBCs but to the near-surface stratification profile and the fact that the packets are measured below the surface. The solutions using the (Atmo SAI 1) and (Reference) BCs are identical, suggesting that the BC appropriately models the exponential decay of the density in the extended atmosphere as intended. Finally, we note that the high-frequency approximation made with the (Atmo HF 1) BC induces a small error (0.6% difference) in the reflection coefficient.

5.2. Approximate atmospheric radiation boundary conditions

We now consider the case where k_h ≠ 0 by placing the source on the rotation axis at the photosphere. In this case, the exact form of the RBCs cannot be computed in the 2D case and thus need development, as discussed in Sect. 3. Here, we compare the high-frequency and the small angle developments of (Atmo Non Local) with the 1D setup. Specifically, we are interested in the error induced by considering the different approximations of (Atmo Non Local).

Fig. 5 Time-distance diagrams $\overline{\mathit{C}}\mathrm{\left(}\mathrm{\Delta }\mathit{,t}\mathrm{\right)}$ using a (Dirichlet) and (Atmo SAI 1) BCs by selecting the frequencies around 3 mHz (left) and 6.5 mHz (right). The cross-covariance is filtered using a phase speed filter centered at 125.2 km s-1 for 3 mHz and 250.4 km s-1 for 6.5 mHz and a width of 12.3 km s-1. The dashed black lines indicate the location of the cuts shown at the bottom in order to emphasize the amplitude ratios between the first and second skip.

Due to the fact that the source is so close to the computational boundary, it is difficult to distinguish between outgoing and ingoing waves. To address this, a phase-speed filter is applied to the power spectrum in order to select skip distances, and thus, angles of incidence. For the cross-covariance at 3 mHz, the phase-speed filter is centered around 125.2 km s^-1. For the cross-covariance at 6.5 mHz, the phase-speed filter is centered around 250.4 km s^-1. A time-distance diagram for these two cases is shown in Fig. 5 for (Dirichlet) and (Atmo SAI 1) boundary conditions. Figure 5 again demonstrates that the behavior of the low-frequency waves is similar for both boundary conditions due to the reflection from the density. At high frequencies, the waves with Dirichlet BC show only a small drop in amplitude with successive skips. However, for the (Atmo SAI 1) BC, the waves diminish rapidly with successive skips, suggesting that the high-frequency waves are absorbed by the BC, as desired. We note that the amplitude of the first skip is higher for the Dirichlet BC as part of the wave packet (near the source) is reflected back due to the boundary condition.

Table 4 details the energy and amplitude ratios between the first and second skip waves. At 3 mHz, the energy of the second skip relative to the first is approximately 11%, while the amplitude is 30%. The difference between these skips at 3 mHz is due to the damping. However, at 6.5 mHz the RBCs show a smaller ratio between the first and second skips at 1.1% for the energy and 10.6% for the amplitude. In the case of (Dirichlet), a similar ratio is reflected to that at 3 mHz. If we compare the different approximations of the BC, one can see that the approximation of a small angle of incidence gives a more accurate description of the atmosphere compared to the high-frequency approximation. The amplitude of the second peak is about 10% of the first skip which is in agreement with the observations by Jefferies et al. (1997).

Since the (Atmo HF 1) BC is valid for very high frequencies, the accuracy of this BC as a function of frequency must be examined. The error in the reflection coefficient is shown in Fig. 6. As expected, the accuracy increases exponentially with frequency but the boundary conditions derived under the small angle of incidence assumption are still more accurate.

Fig. 6 Relative difference between the reflection coefficients as a function of frequency. The different boundary conditions are compared with respect to the (Reference) BC. A phase speed filter centered at 250.4 km s-1 is used in all cases.

In order to study the impact of the angle of incidence of the wave, Fig. 7 shows the difference between the (Atmo SAI 1) and (Atmo SAI 0) RBCs with the (Atmo Non Local) as reference, for a range of phase-speeds. We compute the angle of incidence corresponding to a specific phase-speed using ray theory for a frequency of 3 mHz. As expected from the expansions described in Sect. 3.2, the difference decreases exponentially with the angle of incidence. However, the differences are very small which explain the tiny differences seen in the reflection coefficients in Table 4.

Fig. 7 Relative difference between the reflection coefficients as reported in Table 4 for phase speed filters ranging from 62.6 km s-1 to 250.4 km s-1, that is, angles of incidence with the normal to the surface from 1.76° to 7.57°. The (Atmo Non Local) RBC energy ratio is used as reference.

6. Importance of the atmospheric model

So far, we have considered an atmospheric extension where the density is exponentially decaying and the sound speed is constant. However, in reality the temperature profile increases by two orders of magnitude in the transition region between the chromosphere and the corona. This increase in temperature results in an increase of sound speed and acoustic cutoff frequency by a factor of ten. Jefferies et al. (1997) studied time-distance diagrams at high frequencies using South Pole data and observed a double-ridge structure; they argued that it was due to the reflection of waves in this region. In this section, we study the time-distance diagram in a case where the atmosphere is extended using the VAL-C atmospheric model (see Fig. 2).

To compare our result with the observations by Jefferies et al. (1997), we filter the cross-covariance with the same frequency filter and the same spatial filter. The frequency filter is a Gaussian filter centered at 6.75 mHz with a standard deviation of 0.75 mHz. The spatial filter is centered at ℓ = 125 with a standard deviation σ_ℓ = 33. Figure 8 shows the time-distance diagrams obtained with the atmospheric extension Atmo and the VAL-C model. The models are compared to observed cross-covariance function using 72 days of MDI observations. The same frequency and spatial filters were used in all three cases. With an upper atmosphere of constant sound speed, each skip ridge is formed of a single ridge, while a double-ridge structure is seen in the model with sound speed stratification and in the MDI observations (as observed in the South Pole data).

Jefferies et al. (1997) suggests that the double-ridge structure is due to wave reflection at either 1 or 2 Mm above the photosphere. Sekii et al. (2004) noted that the double-ridge structure may be the result of an interference between waves below and above the acoustic cutoff frequency. Several improvements in our simulations would be needed to further study this question, for example by implementing a wave attenuation model that varies not only with frequency but also with height. Nevertheless, our simulations clearly indicate that high-frequency waves can provide information about the upper solar atmosphere. They also suggest that the upper atmosphere needs to be modeled correctly in order to probe the solar interior with high-frequency waves.

Fig. 8 Time-distance diagram for high-frequency waves (isolated by a Gaussian filter centered at 6.75 mHz and of FWHM 0.75 mHz) using a constant sound speed in the atmosphere (left), the VAL-C atmospheric model (Vernazza et al. 1981) (center), and from MDI observations (right).

7. Conclusion

In this paper we have applied the radiative boundary conditions developed by Barucq et al. (2017) to time-distance helioseismology, which were specifically developed for an atmosphere with a constant sound-speed and an exponentially decaying density. The importance of such boundary conditions in computational helioseismology is in their ability to treat outgoing high-frequency waves. We found that the boundary condition derived under the assumption of a small angle of incidence, Atmo SAI 1, reproduces most accurately the atmospheric model Atmo and does not increase the computational burden when applied a few hundred km above the photosphere. By studying the energy of outgoing and reflected waves we have a good comparison with observations.

The radiative boundary conditions proposed here are much more elegant and less computationally expensive than extending the atmosphere by an absorbing layer (PMLs) or than transitioning to a constant medium together with a simple Sommerfeld boundary condition (Gizon et al. 2017).

In addition to the study of the boundary conditions, we have briefly examined the consequence of the sharp increase in sound speed in the transition region. Using a VAL model, we have reproduced the double-ridge feature seen in the high-frequency time-distance diagram. These results represent a step towards the inclusion of high-frequency waves in computational local helioseismology.

References

All Tables

All Figures

	Fig. 1 Geometrical setup for 2.5D computations. The computational domain is delimited by thick lines and the photosphere by the dashed line.
In the text

Fig. 2 Sound speed (red) and density (blue) in the upper layers of the Sun as functions of height measured from the photosphere. Below 500 km (vertical dashed line), the sound speed and density are from model S (Christensen-Dalsgaard et al. 1996). Above 500 km, the model Atmo (solid lines) consists of an extension where density decays exponentially (H = 105 km) and sound speed is constant (c = 6855 m/s). The dashed lines show the VAL-C atmospheric model up to a height of 2500 km (Vernazza et al. 1981), which will be used in Sect. 6.

In the text

Fig. 3 Acoustic cut-off frequency ωc as a function of height from the photosphere. Below 500 km, ωc is obtained using model S and Eq. (11); the density scale height is obtained after smoothing the density. Above 500 km, the acoustic cut-off frequency is constant and equal to ωt = 5.3 mHz. Near − 200 km, ωc is complex, although ${\mathit{\omega }}_{\mathrm{c}}^{\mathrm{2}}$ is always positive and only ${\mathit{\omega }}_{\mathrm{c}}^{\mathrm{2}}$ enters the wave equation.

In the text

Fig. 4 Temporal cross-covariance $\mathit{c}\mathrm{\left(}\mathit{r}{\mathrm{)}}^{\mathrm{2}}\overline{\mathit{C}}\mathrm{\right(}\mathrm{0}\mathit{,r,t}\mathrm{\left)}\exp \mathrm{\right(}\mathit{t}\mathit{/}\mathrm{90}\min \mathrm{)}$ computed with (Atmo Non Local) and (Dirichlet) boundary conditions, originating from the center, with frequencies selected around 3 mHz (left panel) and 6.5 mHz (right panel). As a comparison, we also plot the cross-covariance when using the (Reference) BC (dashed cyan line). The signal is multiplied by exp(t/ 90min) to counteract the effect of the damping, and normalized by the maximum of its amplitude at t = 30 min. The amplitude is plotted as a function of the acoustic depth $\mathit{\tau }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{{}^{\mathrm{\int }}}_{\mathit{r}}^{{\mathit{r}}_{\mathit{a}}}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}\mathit{/}\mathit{c}\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}$ for better readability of near-surface behavior. The dashed line indicates the solar surface.

In the text

Fig. 5 Time-distance diagrams $\overline{\mathit{C}}\mathrm{\left(}\mathrm{\Delta }\mathit{,t}\mathrm{\right)}$ using a (Dirichlet) and (Atmo SAI 1) BCs by selecting the frequencies around 3 mHz (left) and 6.5 mHz (right). The cross-covariance is filtered using a phase speed filter centered at 125.2 km s-1 for 3 mHz and 250.4 km s-1 for 6.5 mHz and a width of 12.3 km s-1. The dashed black lines indicate the location of the cuts shown at the bottom in order to emphasize the amplitude ratios between the first and second skip.

In the text

	Fig. 6 Relative difference between the reflection coefficients as a function of frequency. The different boundary conditions are compared with respect to the (Reference) BC. A phase speed filter centered at 250.4 km s-1 is used in all cases.
In the text

	Fig. 7 Relative difference between the reflection coefficients as reported in Table 4 for phase speed filters ranging from 62.6 km s-1 to 250.4 km s-1, that is, angles of incidence with the normal to the surface from 1.76° to 7.57°. The (Atmo Non Local) RBC energy ratio is used as reference.
In the text

	Fig. 8 Time-distance diagram for high-frequency waves (isolated by a Gaussian filter centered at 6.75 mHz and of FWHM 0.75 mHz) using a constant sound speed in the atmosphere (left), the VAL-C atmospheric model (Vernazza et al. 1981) (center), and from MDI observations (right).
In the text