Atmosphericradiation boundary conditions for highfrequency waves in timedistance helioseismology
^{1} MaxPlanckInstitut für Sonnensystemforschung, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
email: fournier@mps.mpg.de
^{2} Institut für Astrophysik, GeorgAugustUniversität Göttingen, FriedrichHundPlatz 1, 37077 Göttingen, Germany
^{3} Center for Space Science, NYUAD Institute, New York University Abu Dhabi, PO Box 129188, Abu Dhabi, UAE
^{4} Magique3D, Inria Bordeaux SudOuest, Université de Pau et des Pays de l’Adour, 64013 Pau, France
^{5} Magique3D, Inria Bordeaux SudOuest, Université de Bordeaux, 33400 Talence, France
Received: 31 May 2017
Accepted: 6 September 2017
The temporal covariance between seismic waves measured at two locations on the solar surface is the fundamental observable in timedistance helioseismology. Above the acoustic cutoff frequency (~5.3 mHz), waves are not trapped in the solar interior and the covariance function can be used to probe the upper atmosphere. We wish to implement appropriate radiative boundary conditions for computing the propagation of highfrequency waves in the solar atmosphere. We consider recently developed and published radiative boundary conditions for atmospheres in which soundspeed is constant and density decreases exponentially with radius. We compute the crosscovariance function using a finite element method in spherical geometry and in the frequency domain. The ratio between first and secondskip amplitudes in the timedistance diagram is used as a diagnostic to compare boundary conditions and to compare with observations. We find that a boundary condition applied 500 km above the photosphere and derived under the approximation of small angles of incidence accurately reproduces the “infinite atmosphere” solution for highfrequency waves. When the radiative boundary condition is applied 2 Mm above the photosphere, we find that the choice of atmospheric model affects the timedistance diagram. In particular, the timedistance diagram exhibits doubleridge structure when using a Vernazza Avrett Loeser atmospheric model.
Key words: Sun: helioseismology / Sun: atmosphere / methods: numerical
© 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 freesurface boundary condition (e.g., ChristensenDalsgaard et al. 1996; Bogdan et al. 1996); but this approach forces all waves to reflect and thus is not appropriate for highfrequency 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 userdefined 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^{2}(r)∇·ξ^{ξ}ξ(r,ω), which in the frequency domain is given by (1)where (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 soundspeed c, and the density ρ, which are read from the standard solar model S (ChristensenDalsgaard 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 solarlike 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 timedistance helioseismology.
In Sect. 2, we describe the equations that we are solving for onedimensional (1D) and twodimensional (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 timedistance diagram at high frequencies.
2. Scalar wave equations
2.1. Twodimensional 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 (3)Each mode ψ_{m} is then the solution of a 2D partial differential equation, (4)where and is the 2D spatial gradient operator and 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. 

Open with DEXTER 
2.2. 1D background models
When the background model depends only on radius, the solution can be written in terms of spherical harmonic functions: (5)For each ℓ and m the function ψ_{ℓm}(r) is defined on the 1D domain 0 ≤ r ≤ r_{t} and solves (6)where (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 (ChristensenDalsgaard 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 VALC are given in Fig. 2. A smooth transition between the model S and the VALC 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 (ChristensenDalsgaard 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 VALC atmospheric model up to a height of 2500 km (Vernazza et al. 1981), which will be used in Sect. 6. 

Open with DEXTER 
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 highfrequency waves. This can be done by changing the unknown, (8)resulting in a new form for Eq. (6), (9)where (10)and ω_{c} is a cutoff frequency defined by (11)where (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 lowfrequency waves (ω<ω_{t}) are reflected back in the Sun, while the highfrequency ones (ω>ω_{t}) tunnel through the peak at 100 km below the surface and continue propagating into the upper atmosphere.
Fig. 3 Acoustic cutoff 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 cutoff frequency is constant and equal to ω_{t} = 5.3 mHz. Near − 200 km, ω_{c} is complex, although is always positive and only enters the wave equation. 

Open with DEXTER 
With Eq. (9)in hand, a radiative boundary condition is given by (13)We note that Q can be negative and thus Q^{1 / 2} is imaginary for frequencies above the acoustic cutoff ω_{t}. Using Eqs. (8)and (13), we can write the radiation boundary condition for the initial unknown ψ_{lm} as 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 present in Q should be replaced by the horizontal Laplacian Δ_{h}(14)However, since this term is located under the square root, it leads to a spatially nonlocal operator which cannot be implemented directly. In order to overcome this difficulty, an asymptotic expansion of the squareroot 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 highfrequency approximation by making a firstorder approximation of the square root supposing that ω_{t}/ω is small, 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: The general form of this BC is similar to (Atmo HF 1) since the radial derivative at the boundary 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 weakform (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 while not expanding around 1 /σ^{2}, as was done for (Atmo HF 1). This leads to the last BC we consider in this study: This condition is again a Sommerfeldlike 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 planeparallel isothermal atmosphere (15)The term in 1 /r_{t} is missing due to the planeparallel approximation and the acoustic cutoff ω_{c}_{0} = 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 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: 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 highfrequency or small angle of incidence approximation.
Summary of the different boundary conditions used in this paper with the associated domain of validity and their applicability for 1D and 2D backgrounds.
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 (16)The term containing the Laplacian is then integrated by parts (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^{1} 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 crosscovariance function. Gizon et al. (2017) showed that if the sources are spatially uncorrelated and equipartitioned then the expectation value of the crosscovariance is related to the imaginary part of the Green’s function (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, (19)where ω_{0}/ 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 10day 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 crosscovariance; one for each skip, (20)We measure the energy of the outgoing (first skip) and reflected (second skip) wave packets as follows: (21)where the intervals select the first and secondskip 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, 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 crosscovariance 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 crosscovariance 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 for better readability of nearsurface behavior. The dashed line indicates the solar surface. 

Open with DEXTER 
Figure 4 shows the spatial evolution of the crosscovariance for wave packets centered at ω_{0}/ 2π = 3 mHz and ω_{0}/ 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 lowfrequency 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 lowfrequency 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 lowfrequency 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 nearsurface 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 highfrequency approximation made with the (Atmo HF 1) BC induces a small error (0.6% difference) in the reflection coefficient.
Ratios between outgoing and reflected energies (ℰ_{2}/ ℰ_{1}) and amplitudes (A_{2}/A_{1}) for a source at the center of the Sun, r_{s} = 0, and a receiver at r = 0.8 R_{⊙}, for selected ranges of frequencies around 3 mHz and 6.5 mHz.
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 highfrequency 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 Timedistance diagrams using a (Dirichlet) and (Atmo SAI 1) BCs by selecting the frequencies around 3 mHz (left) and 6.5 mHz (right). The crosscovariance 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. 

Open with DEXTER 
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 phasespeed filter is applied to the power spectrum in order to select skip distances, and thus, angles of incidence. For the crosscovariance at 3 mHz, the phasespeed filter is centered around 125.2 km s^{1}. For the crosscovariance at 6.5 mHz, the phasespeed filter is centered around 250.4 km s^{1}. A timedistance 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 lowfrequency 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 highfrequency 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 highfrequency 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. 

Open with DEXTER 
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 phasespeeds. We compute the angle of incidence corresponding to a specific phasespeed 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. 

Open with DEXTER 
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 timedistance diagrams at high frequencies using South Pole data and observed a doubleridge structure; they argued that it was due to the reflection of waves in this region. In this section, we study the timedistance diagram in a case where the atmosphere is extended using the VALC atmospheric model (see Fig. 2).
To compare our result with the observations by Jefferies et al. (1997), we filter the crosscovariance 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 timedistance diagrams obtained with the atmospheric extension Atmo and the VALC model. The models are compared to observed crosscovariance 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 doubleridge 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 doubleridge structure is due to wave reflection at either 1 or 2 Mm above the photosphere. Sekii et al. (2004) noted that the doubleridge 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 highfrequency 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 highfrequency waves.
Fig. 8 Timedistance diagram for highfrequency 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 VALC atmospheric model (Vernazza et al. 1981) (center), and from MDI observations (right). 

Open with DEXTER 
7. Conclusion
In this paper we have applied the radiative boundary conditions developed by Barucq et al. (2017) to timedistance helioseismology, which were specifically developed for an atmosphere with a constant soundspeed and an exponentially decaying density. The importance of such boundary conditions in computational helioseismology is in their ability to treat outgoing highfrequency 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 doubleridge feature seen in the highfrequency timedistance diagram. These results represent a step towards the inclusion of highfrequency waves in computational local helioseismology.
References
 Barucq, H., Chabassier, J., Duruflé, M., Gizon, L., & Leguèbe, M. 2017, Math. Modell. Num. Anal., submitted http://hal.inria.fr/hal01581834 [Google Scholar]
 Bayliss, A., Gunzburger, M., & Turkel, E. 1982, SIAM J. Appl. Math., 42, 430 [CrossRef] [Google Scholar]
 Berenger, J.P. 1994, J. Comput. Phys., 114, 185 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Bogdan, T. J., Hindman, B. W., Cally, P. S., & Charbonneau, P. 1996, ApJ, 465, 406 [NASA ADS] [CrossRef] [Google Scholar]
 Cameron, R., Gizon, L., & Duvall, Jr., T. L. 2008, Sol. Phys., 251, 291 [NASA ADS] [CrossRef] [Google Scholar]
 ChristensenDalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Gizon, L., Barucq, H., Duruflé, M., et al. 2017, A&A, 600, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gough, D. 1993, in Astrophysical Fluid DynamicsLes Houches 1987, Vol. 1, 399 [Google Scholar]
 Hanasoge, S. M., Komatitsch, D., & Gizon, L. 2010, A&A, 522, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jefferies, S. M., Osaki, Y., Shibahashi, H., et al. 1997, ApJ, 485, L49 [NASA ADS] [CrossRef] [Google Scholar]
 Komatitsch, D., & Martin, R. 2007, Geophysics, 72, SM155 [Google Scholar]
 Santamaria, I. C., Khomenko, E., & Collados, M. 2015, A&A, 577, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schunker, H., Cameron, R. H., Gizon, L., & Moradi, H. 2011, Sol. Phys., 271, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Sekii, T., Shibahashi, H., & Jefferies, S. M. 2004, in SOHO 14 Helio and Asteroseismology: Towards a Golden Future, ed. D. Danesy, ESA SP, 559, 619 [NASA ADS] [Google Scholar]
 Vernazza, J. E., Avrett, E. H., & Loeser, R. 1981, ApJS, 45, 635 [NASA ADS] [CrossRef] [Google Scholar]
 Zienkiewicz, O., Taylor, R., & Zhu, J. 2005, The Finite Element Method Set, sixth edn. (Oxford: ButterworthHeinemann) [Google Scholar]
All Tables
Summary of the different boundary conditions used in this paper with the associated domain of validity and their applicability for 1D and 2D backgrounds.
Ratios between outgoing and reflected energies (ℰ_{2}/ ℰ_{1}) and amplitudes (A_{2}/A_{1}) for a source at the center of the Sun, r_{s} = 0, and a receiver at r = 0.8 R_{⊙}, for selected ranges of frequencies around 3 mHz and 6.5 mHz.
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. 

Open with DEXTER  
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 (ChristensenDalsgaard 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 VALC atmospheric model up to a height of 2500 km (Vernazza et al. 1981), which will be used in Sect. 6. 

Open with DEXTER  
In the text 
Fig. 3 Acoustic cutoff 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 cutoff frequency is constant and equal to ω_{t} = 5.3 mHz. Near − 200 km, ω_{c} is complex, although is always positive and only enters the wave equation. 

Open with DEXTER  
In the text 
Fig. 4 Temporal crosscovariance 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 crosscovariance 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 for better readability of nearsurface behavior. The dashed line indicates the solar surface. 

Open with DEXTER  
In the text 
Fig. 5 Timedistance diagrams using a (Dirichlet) and (Atmo SAI 1) BCs by selecting the frequencies around 3 mHz (left) and 6.5 mHz (right). The crosscovariance 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. 

Open with DEXTER  
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. 

Open with DEXTER  
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. 

Open with DEXTER  
In the text 
Fig. 8 Timedistance diagram for highfrequency 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 VALC atmospheric model (Vernazza et al. 1981) (center), and from MDI observations (right). 

Open with DEXTER  
In the text 