© ESO, 2010

1. Introduction

The choice of appropriate boundary conditions in simulations remains a significant challenge in computational physics. Amidst the vast number of options that exist, a widely invoked construct is that of the absorbing boundary, one that is expected to relieve the computational domain of outgoing waves or other structures without affecting the solution in the region of interest (for a review, see e.g., Colonius 2004). A number of different methods exist to accomplish this, such as an attenuating “sponge” (e.g., Lui 2003; Colonius 2004), characteristics-based boundary conditions (e.g. Thompson 1990), an artificially imposed supersonic outward directed advective flow (e.g., Lui 2003), perfectly matched layers (PMLs; Bérenger 1994), etc. Unfortunately, a number of these methods is afflicted with disadvantages, a primary cause of concern being that of absorption efficiency. For example, the characteristics-based boundary conditions do not perform very well when the incident waves at the boundary are significantly inclined – and a fair degree of reflection is observed. The sponges are better at absorbing the waves and are relatively easy to implement; however, the reflectivity can still be substantial.

The criterion of low reflectivity is perhaps best satisfied by the PML formulation, developed first by Bérenger (1994) as an absorbing layer for the Maxwell equations. The central idea of this technique lies in performing an analytic continuation of the wave vector into the complex plane. Wave vectors perpendicular to the boundary are forced to assume a decaying form in the complex plane, thereby dramatically reducing the amplitudes of the waves in the boundary PML region. In the absence of discretization errors, the PML as set out by Bérenger (1994) is highly absorbent. After numerical discretization however, there are weak associated reflections. An important issue in the classical split formulation of Bérenger (1994) is that the absorption efficiency decreases rapidly at grazing incidence (e.g., Collino & Monk 1998; Winton & Rappaport 2000).

There have been numerous advances in constructing stable analytic continuations of wave vectors in the boundary region. In particular, the stable convolutional PMLs (C-PML; e.g., Roden & Gedney 2000; Festa & Vilotte 2005; Komatitsch & Martin 2007; Drossaert & Giannopoulos 2007), contain a Butterworth filter inside the PML, thereby dramatically improving the absorption efficiency at grazing incidence. This formalism was adopted to the study of waves in anisotropic geophysical media by Komatitsch & Martin (2007); the numerical stability of said technique was studied by Komatitsch & Martin (2007); )Meza-Fajardo & Papageorgiou (2008). Extensions to the poro- and viscoelastic cases have been introduced by Martin et al. (2008) and Martin & Komatitsch (2009), respectively. Our goal in this contribution is to develop a C-PML for the 3D linearized ideal MHD equations in stratified media.

Astrophysical media such as stellar interiors and atmospheres may be strongly stratified and highly magnetized. Simulating wave propagation in such environments is of interest because in understanding the effects of stratification and magnetic fields on the oscillations, we may be able to better interpret observations and constrain certain properties of the object in question. In particular, we focus on the Sun, a star that has been well studied and for which high-quality observations exist. The interior of the Sun is opaque to electromagnetic waves – consequently, only photons that arrive from very close to the surface (photosphere) of the Sun are visible. Helioseismology is the inference of the internal structure and dynamics of the Sun by observing the surface manifestation of its acoustic pulsations (e.g., Christensen-Dalsgaard 2002; Gizon & Birch 2005; Gizon et al. 2010). Armed with accurate interaction theories of waves and measurements of the acoustic field at the surface, one can attempt to determine the sub-surface structure and dynamics of various solar features such as sunspots (see e.g., Gizon et al. 2009), large-scale meridional flow, interior convective length scales, etc. In order to construct these interaction theories, it is important to simulate and study small amplitude (linear) wave propagation in a solar-like medium – filled with scatterers like sunspots or mean flows. The time scales of wave propagation are typically much smaller than the rate at which the scatterer itself evolves; thus one may invoke the assumption of time stationarity of the background medium. Numerical simulations of waves in solar-like media have been performed by numerous groups (e.g., Werne et al. 2004; Hanasoge et al. 2006; Shelyag et al. 2006; Khomenko & Collados 2006; Cameron et al. 2007; Parchevsky & Kosovichev 2007; Hanasoge et al. 2008; Hanasoge 2008; Cameron et al. 2008). Two groups in particular (Khomenko & Collados 2006; Parchevsky & Kosovichev 2007) currently utilize the classical split PML formulation in order to solve the wave equations in stratified media. There exist drawbacks in these formulations. For example, Khomenko & Collados (2006) have noted that there are long term instabilities associated with their method. The technique discussed in Parchevsky & Kosovichev (2007) involves the addition of a small arbitrary constant (see Eq. (21) and related discussion of their article) that could possibly be acting as a sponge; it is unclear if their method is perfectly matched after the introduction of this constant. The modal power spectrum that Parchevsky & Kosovichev (2007) display in Fig. 9 of their article shows substantial reflected power from their lower boundary, atypical of perfectly matched formulations.

The plan of this article is as follows: in Sect. 2, we recall the linearized ideal MHD equations and discuss the C-PML method as was previously developed for the seismic wave equations. The C-PML is then extended to the stratified Euler/MHD equations in Sect. 3 and results from numerical tests are discussed. We summarize and conclude in Sect. 4.

2. The linearized ideal MHD equations

As discussed in the introduction, we consider a system with a steady background, the wave equations are written as small perturbations around this state. The assumption of small amplitude wave perturbations is consistent with observations of solar and stellar oscillations (e.g., Christensen–Dalsgaard 2003; Bogdan 2000). In the equations that follow, all quantities with a subscript “0” denote the steady background properties while the unsubscripted terms represent the fluctuations around the corresponding property. All bold faced terms represent vectors. $\begin{array}{ccc}{\mathit{\partial }}_{\mathit{t}}\mathit{\rho }& \mathrm{=}& \\ {\mathit{\partial }}_{\mathit{t}}\mathrm{\left(}{\mathit{\rho }}_{\mathrm{0}}{v}\mathrm{\right)}& \mathrm{=}& \mathbf{\nabla }\mathrm{·}\left[\frac{{B}{{B}}_{\mathrm{0}}}{\mathrm{4}\mathit{\pi }}\mathrm{+}\frac{{{B}}_{\mathrm{0}}{B}}{\mathrm{4}\mathit{\pi }}\mathrm{-}\left(\mathit{p}\mathrm{+}\frac{{{B}}_{\mathrm{0}}\mathrm{·}{B}}{\mathrm{4}\mathit{\pi }}\right)\underline{{I}}\right]\mathrm{-}\mathit{\rho }{\mathit{g}}_{\mathrm{0}}{{e}}_{\mathit{z}}\\ & \mathrm{+}& \\ {\mathit{\partial }}_{\mathit{t}}\mathit{p}& \mathrm{=}& \\ {\mathit{\partial }}_{\mathit{t}}{B}& \mathrm{=}& \\ \mathbf{\nabla }\mathrm{·}{B}& \mathrm{=}& \end{array}$In order, these are the equations of continuity (Eq. (1)), momentum (Eq. 2)), adiabaticity (Eq. (3)), induction (Eq. (4)), and divergence-free field (Eq. (5)). The terms ρ and p denote density and pressure, B = (B_x,B_y,B_z) is the magnetic field, g₀ is the gravity, c₀ the sound speed, and v = (v_x,v_y,v_z) is the velocity. A Cartesian coordinate system (x,y,z) with unit vectors (e_x,e_y,e_z) is adopted. The identity tensor is . The momentum equations are forced by an arbitrary spatio-temporally smooth source function S. All background terms p₀,ρ₀,c,B₀ are assumed to be heterogeneous, i.e., functions of (x,y,z) and g₀ = g₀(z). Note that we have assumed the background medium to be non-moving, i.e., that v₀ = 0. The extension to the v₀ ≠ 0 situation requires no additional treatment at the boundary as long as the background velocity is zero in this region.

2.1. PMLs and C-PMLs

Consider a wave propagating in the z direction, towards the upper boundary. For now, ignoring the fact that the medium is stratified, we have v_z ~ Ae^{i(k_zz − ωt)}, where A is the amplitude of the wave, k_z is the wavenumber along the vertical direction, ω the frequency, and t, time. The classical PML of Bérenger (1994) prescribes the following transformation for k_z: ${\mathit{k}}_{\mathit{z}}\mathrm{\to }{\mathit{k}}_{\mathit{z}}\left[\mathrm{1}\mathrm{-}\frac{\mathit{d}}{\mathit{i\omega }}\right]\mathit{,}$(6)where d = d(z)f(x,y,z₀) ≥ 0 is a damping parameter and z = z₀ is the vertical coordinate corresponding to the entrance of the layer. We include the term f(x,y,z₀) to account for the possibility of strong wave speed variations in the horizontal (non-PML) directions at z = z₀ (i.e., at the entrance of the PML). Note that the damping term applies only to k_z; i.e. the wavenumbers k_x and k_y remain unchanged. It has been shown that, cast in the form of Eq. (6), the PML is unstable in a number of situations (for f(x,y,z₀) = 1), including in the presence of mean flows (e.g., Hu 2001; Appelö et al. 2006), and in anisotropic media (e.g., Bécache et al. 2003; Komatitsch & Martin 2007, also personal communication, Khomenko 2009). A version of the PML with a more general transformation was introduced by Roden & Gedney (2000) for Maxwell equations and adapted to the seismic wave equation by Festa & Vilotte (2005), Komatitsch & Martin (2007), and Drossaert & Giannopoulos (2007). They applied the following relation: ${\mathit{k}}_{\mathit{z}}\mathrm{\to }{\mathit{k}}_{\mathit{z}}\left[\mathit{\kappa }\mathrm{+}\frac{\mathit{d}}{\mathit{\alpha }\mathrm{-}\mathit{i\omega }}\right]\mathit{,}$(7)where α = α(z) > 0 and κ = κ(z) ≥ 1 are additional parameters. In particular, the inclusion of α, i.e. of a filter, helped improve the issue of damping waves with grazing incidence at the boundaries. We follow the formalism set out in Komatitsch & Martin (2007) and derive the equations for the convolutional PML. The transformation of Eq. (7) corresponds to a grid stretching of the relevant coordinate: $\begin{array}{}\mathrm{˜}\\ \mathit{z}\end{array}\mathrm{=}\left(\mathit{\kappa }\mathrm{+}\frac{\mathit{d}}{\mathit{\alpha }\mathrm{-}\mathit{i\omega }}\right)\mathit{z}\mathit{.}$(8)Derivatives along the vertical direction in Eqs. (1) through (5) are now computed in terms of the stretched coordinate. In other words, the vertical derivative of a generic variable ψ(x,y,z,t) is transformed according to: ${\mathit{\partial }}_{\mathit{z}}\mathit{\psi }\mathrm{\to }{\mathit{\partial }}_{\begin{array}{}\mathrm{˜}\\ \mathit{z}\end{array}}\mathit{\psi }\mathit{.}$(9)The term may be rewritten as: $\begin{array}{ccc}{\mathit{\partial }}_{\begin{array}{}\mathrm{˜}\\ \mathit{z}\end{array}}\mathit{\psi }& \mathrm{=}& \end{array}$(10)The first term inside the parenthesis on the right-hand side of Eq. (10) is easily obtained by dividing the derivative by κ. The second term is more complicated because it involves products in the temporal Fourier space of two terms with frequency content; this results in a convolution in time. In order to compute these convolutions, Roden & Gedney (2000); ); )Festa & Vilotte (2005); )Komatitsch & Martin (2007) use a recursion formula to update the convolution at each time step. Instead of applying a recursive relation, here we use a differential equation to recover (also done in the auxiliary differential equation formulation of Gedney & Zhao 2010; Martin et al. 2010), where $\mathit{s̅}\mathrm{\left(}\mathit{x,y,z,\omega }\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathit{d}\mathit{/}{\mathit{\kappa }}^{\mathrm{2}}}{\mathrm{(}\mathit{d}\mathit{/}\mathit{\kappa }\mathrm{+}\mathit{\alpha }\mathrm{)}\mathrm{-}\mathit{i\omega }}\mathrm{·}$(11)Let the time history of χ be such that χ(x,y,z,t ≤ 0) = 0. In addition, we also require ∂_zψ = 0 for all t ≤ 0. Then, the following differential equation for χ${\mathit{\partial }}_{\mathit{t}}\mathit{\chi }\mathrm{=}\mathrm{-}\frac{\mathit{d}}{{\mathit{\kappa }}^{\mathrm{2}}}{\mathit{\partial }}_{\mathit{z}}\mathit{\psi }\mathrm{-}\left(\frac{\mathit{d}}{\mathit{\kappa }}\mathrm{+}\mathit{\alpha }\right)\mathit{\chi ,}$(12)ensures that χ has the desired frequency response, i.e., $\mathit{\chi }\mathrm{\left(}\mathit{x,y,z,\omega }\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathit{d}\mathit{/}{\mathit{\kappa }}^{\mathrm{2}}}{\mathrm{(}\mathrm{d}\mathit{/}\mathit{\kappa }\mathrm{+}\mathit{\alpha }\mathrm{)}\mathrm{-}\mathit{i\omega }}{\mathit{\partial }}_{\mathit{z}}\mathit{\psi }\mathit{.}$(13)The differential equations for the C-PML may now be written as: $\begin{array}{ccc}{\mathit{\partial }}_{\mathit{t}}\mathit{\rho }& \mathrm{=}& \\ {\mathit{\partial }}_{\mathit{t}}\left({\mathit{\rho }}_{\mathrm{0}}{v}\right)& \mathrm{=}& {\mathbf{\nabla }}_{\mathit{h}}\mathrm{·}\left[\frac{{B}{{B}}_{\mathrm{0}}}{\mathrm{4}\mathit{\pi }}\mathrm{+}\frac{{{B}}_{\mathrm{0}}{B}}{\mathrm{4}\mathit{\pi }}\mathrm{-}\left(\mathit{p}\mathrm{+}\frac{{{B}}_{\mathrm{0}}\mathrm{·}{B}}{\mathrm{4}\mathit{\pi }}\right){I}\right]\\ & \mathrm{+}& {\mathit{\partial }}_{\begin{array}{}\mathrm{˜}\\ \mathit{z}\end{array}}\left[\frac{{\mathit{B}}_{\mathit{z}}{{B}}_{\mathrm{0}}}{\mathrm{4}\mathit{\pi }}\mathrm{+}\frac{{\mathit{B}}_{\mathrm{0}\mathit{z}}{B}}{\mathrm{4}\mathit{\pi }}\mathrm{-}\left(\mathit{p}\mathrm{+}\frac{{{B}}_{\mathrm{0}}\mathrm{·}{B}}{\mathrm{4}\mathit{\pi }}\right){{e}}_{\mathit{z}}\right]\\ & \mathrm{-}& \mathit{\rho }\begin{array}{}\mathrm{˜}\\ {\mathit{g}}_{\mathrm{0}}\end{array}{{e}}_{\mathit{z}}\mathrm{-}\mathit{\sigma }{\mathit{\rho }}_{\mathrm{0}}{v}\mathrm{+}{S}\mathit{,}\\ {\mathit{\partial }}_{\mathit{t}}\mathit{p}& \mathrm{=}& \\ {\mathit{\partial }}_{\mathit{t}}{B}& \mathrm{=}& \end{array}$where represents the derivatives in the non C-PML directions, and are the modified background gravity, pressure and density gradients in the absorption region respectively. A sponge-like decay term  − σρ₀v, where σ = σ(x,y,z) is introduced in the momentum equation in order to stabilize the system. It is important to note that in the absorbing layer, the solution is non-physical and therefore discarded. Thus any modifications applied to the equations, in relation to issues such as altering the stratification or not enforcing ∇·B = 0, is justifiable as long as the solution in the interior domain is unaffected and the formulation is stable. The Eqs. (14) through (17) in conjunction with Eq. (12) provide the following prescription: $\begin{array}{ccc}{\mathit{\partial }}_{\mathit{t}}\mathit{\rho }& \mathrm{=}& \\ {\mathit{\partial }}_{\mathit{t}}\mathrm{\Psi }& \mathrm{=}& \\ {\mathit{\partial }}_{\begin{array}{}\mathrm{˜}\\ \mathit{z}\end{array}}{\mathit{\rho }}_{\mathrm{0}}& \mathrm{=}& \\ {\mathit{\partial }}_{\mathit{t}}\left({\mathit{\rho }}_{\mathrm{0}}{v}\right)& \mathrm{=}& {\mathbf{\nabla }}_{\mathit{h}}\mathrm{·}\left[\frac{{B}{{B}}_{\mathrm{0}}}{\mathrm{4}\mathit{\pi }}\mathrm{+}\frac{{{B}}_{\mathrm{0}}{B}}{\mathrm{4}\mathit{\pi }}\mathrm{-}\left(\mathit{p}\mathrm{+}\frac{{{B}}_{\mathrm{0}}\mathrm{·}{B}}{\mathrm{4}\mathit{\pi }}\right)\underline{{I}}\right]\\ & \mathrm{+}& \frac{\mathrm{1}}{\mathit{\kappa }}{\mathit{\partial }}_{\mathit{z}}\left[\frac{{\mathit{B}}_{\mathit{z}}{{B}}_{\mathrm{0}}}{\mathrm{4}\mathit{\pi }}\mathrm{+}\frac{{\mathit{B}}_{\mathrm{0}\mathit{z}}{B}}{\mathrm{4}\mathit{\pi }}\mathrm{-}\left(\mathit{p}\mathrm{+}\frac{{{B}}_{\mathrm{0}}\mathrm{·}{B}}{\mathrm{4}\mathit{\pi }}\right){{e}}_{\mathit{z}}\right]\\ & \mathrm{+}& {\xi }\mathrm{-}\mathit{\rho }\begin{array}{}\mathrm{˜}\\ {\mathit{g}}_{\mathrm{0}}\end{array}{{e}}_{\mathit{z}}\mathrm{-}\mathit{\sigma }{\mathit{\rho }}_{\mathrm{0}}{v}\mathrm{+}{S}\mathit{,}\\ {\mathit{\partial }}_{\mathit{t}}{\xi }& \mathrm{=}& \mathrm{-}\frac{\mathit{d}}{{\mathit{\kappa }}^{\mathrm{2}}}{\mathit{\partial }}_{\mathit{z}}\left[\frac{{\mathit{B}}_{\mathit{z}}{{B}}_{\mathrm{0}}}{\mathrm{4}\mathit{\pi }}\mathrm{+}\frac{{\mathit{B}}_{\mathrm{0}\mathit{z}}{B}}{\mathrm{4}\mathit{\pi }}\mathrm{-}\left(\mathit{p}\mathrm{+}\frac{{{B}}_{\mathrm{0}}\mathrm{·}{B}}{\mathrm{4}\mathit{\pi }}\right){{e}}_{\mathit{z}}\right]\\ & \mathrm{-}& \\ \begin{array}{}\mathrm{˜}\\ {\mathit{g}}_{\mathrm{0}}\end{array}& \mathrm{=}& \\ {\mathit{\partial }}_{\mathit{t}}\mathit{p}& \mathrm{=}& \mathrm{-}{\mathit{c}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}{\mathbf{\nabla }}_{\mathit{h}}\mathbf{·}{v}\mathrm{-}{\mathit{c}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}\left[\frac{\mathrm{1}}{\mathit{\kappa }}{\mathit{\partial }}_{\mathit{z}}{\mathit{v}}_{\mathit{z}}\mathrm{+}\mathrm{\Psi }\right]\\ & \mathrm{-}& {v}\mathbf{·}{\mathbf{\nabla }}_{\mathit{h}}{\mathit{p}}_{\mathrm{0}}\mathrm{-}{\mathit{v}}_{\mathit{z}}{\mathit{\partial }}_{\begin{array}{}\mathrm{˜}\\ \mathit{z}\end{array}}{\mathit{p}}_{\mathrm{0}}\mathit{,}\\ {\mathit{\partial }}_{\begin{array}{}\mathrm{˜}\\ \mathit{z}\end{array}}{\mathit{p}}_{\mathrm{0}}& \mathrm{=}& \\ {\mathit{\partial }}_{\mathit{t}}{B}& \mathrm{=}& {\mathbf{\nabla }}_{\mathit{h}}\mathrm{·}\mathrm{(}{{B}}_{\mathrm{0}}{v}\mathrm{-}{v}{{B}}_{\mathrm{0}}\mathrm{)}\\ & \mathrm{+}& \frac{\mathrm{1}}{\mathit{\kappa }}{\mathit{\partial }}_{\mathit{z}}\left({\mathit{B}}_{\mathrm{0}\mathit{z}}{v}\mathrm{-}{\mathit{v}}_{\mathit{z}}{{B}}_{\mathrm{0}}\right)\mathrm{+}{\eta }\mathit{,}\\ {\mathit{\partial }}_{\mathit{t}}{\eta }& \mathrm{=}& \end{array}$Note that ξ,η are vector memory variables required for the calculation of the convolutional part of Eq. (10) and that η = (η_x,η_y,0). These are additional variables are needed in order to solve the auxiliary differential Eqs. (19), (22), and (27).

Why is the stratification altered in the boundary region? Consider for example, the derivative term . The differential Eq. (12) when applied to the steady ∂_zp₀ term leads to: $\begin{array}{ccc}{\mathit{\partial }}_{\begin{array}{}\mathrm{˜}\\ \mathit{z}\end{array}}{\mathit{p}}_{\mathrm{0}}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{\kappa }}{\mathit{\partial }}_{\mathit{z}}{\mathit{p}}_{\mathrm{0}}\mathrm{+}\mathit{\chi }\\ {\mathit{\partial }}_{\mathit{t}}\mathit{\chi }& \mathrm{=}& \end{array}$Since we are dealing with a steady background medium, we drop the time derivative in Eq. (29) and obtain ${\mathit{\partial }}_{\begin{array}{}\mathrm{˜}\\ \mathit{z}\end{array}}{\mathit{p}}_{\mathrm{0}}\mathrm{=}\frac{\mathit{\alpha }\mathit{/}\mathit{\kappa }}{\mathit{d}\mathit{/}\mathit{\kappa }\mathrm{+}\mathit{\alpha }}{\mathit{\partial }}_{\mathit{z}}{\mathit{p}}_{\mathrm{0}}\mathit{.}$(30)Note that we can pursue a similar formalism for the gradient term . To maintain hydrostatic stability in the absorption layer, we alter gravity by prefixing it by the same term, . And hence the relations in Eqs. (20), (23), and (25).

In the presence of magnetic field, three wave branches appear: magneto-sonic slow, fast, and Alfvén (the distinctions between these modes apply only if the magnetic and hydrodynamic pressures are not in equipartition). The physics of the propagation of these waves is intricate and inherently highly anisotropic (e.g., Goedbloed & Poedts 2004). Alfvén waves are transverse incompressible oscillations whereas the slow and fast modes are compressive magnetically guided/influenced oscillations. The slow mode dispersion relation exhibits certain pathologies wherein for certain wavenumbers, the phase speed vanishes while the group speed becomes very small (e.g., see pp. 195–214 of Goedbloed & Poedts 2004). This strong anisotropy destabilizes the calculation in the entry regions of the boundary layer. Thus to offset the possibility of this instability, we add the sponge term  − σρ₀v to the momentum equations in regions of the magnetic field.

The choice of the parameters κ,d,α is fairly central to creating an efficient absorption boundary formulation. Following Komatitsch & Martin (2007), we set α = πf₀, where f₀ is a characteristic frequency of the waves. The α decays linearly to zero over the absorption region, being maximum at the start of the layer and falling to zero at the boundary. The damping function d = f(x,y,z₀)   (z / L)^N, is zero at the boundary layer entry point and attains its maximum value at the boundary. We find that for the Euler equations, the following choice for f results in efficient absorption (e.g., Collino & Tsogka 2001): $\mathit{f}\mathrm{\left(}\mathit{x,y,}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathit{N}\mathrm{+}\mathrm{1}}{\mathrm{2}\mathit{L}}\mathit{c̅}{\log }_{\mathrm{10}}{\mathit{R}}_{\mathit{c}}\mathit{,}$(31)where L is the length of the layer, R_c is the amount of reflection tolerated, and is the characteristic sound speed in the C-PML layer. Note that is not a function of z. For the ideal MHD equations, we use: $\mathit{f}\mathrm{\left(}\mathit{x,y,}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathit{N}\mathrm{+}\mathrm{1}}{\mathrm{2}\mathit{L}}{\mathit{c}}_{\mathit{w}}{\log }_{\mathrm{10}}{\mathit{R}}_{\mathit{c}}\mathit{,}$(32)where, is the characteristic propagation speed of the fastest waves in the absorption layer (set for a specific vertical location z = z₀) and R_c,N,L retain the same meaning as for the Euler equations. Wave propagation speeds can vary substantially from strongly magnetized regions to non-magnetic regions; for this reason, we allow c_w to vary horizontally, i.e., c_w = c_w(x,y,z₀). The sound speed is denoted by c and the Alfvén speed by . The κ term acts to damp evanescent waves (Roden & Gedney 2000; Bérenger 2002); although there are evanescent waves in our simulations, we find that instabilities are set into motion when κ is allowed to vary from 1 at the absorption layer entry to 8 at the boundary. For values of κ less than this, the improvement in absorption efficiency is not very noticeable. We therefore set κ = 1. The damping sponge is σ =  −  [ (N + 1)c_A / L ] (z / L)^Nlog ₁₀R_c. The perfectly matched status of the MHD formulation comes into question with the introduction of the sponge term. Note that because σ ∝ c_A the technique is perfectly matched for the stratified Euler equations, i.e. for the c_A = 0 case.

3. Numerical results

3.1. Methods

The simulation code employed in these calculations has been discussed in various previous publications (Hanasoge 2007; Hanasoge et al. 2008; Hanasoge 2008). The validation and verification of the numerical aspects of the code may also be found in these articles. In summary, to capture variations in the vertical direction, we apply sixth order, low dissipation and dispersion compact finite differences (Lele 1992). A non-uniformly spaced grid is used because the eigenfunctions change much more gradually deeper within the interior; the vertical grid spacing varies from several hundred kilometers deeper within to tens of kilometers in the near-surface layers. Because the dissipation of the scheme is low, spectral blocking occurs along the vertical directions (energy buildup at the highest wavenumbers). This is a direct consequence of the strong vertical stratification (terms like c,ρ₀,p₀ vary strongly with z) and the non-uniform grid spacing; in order to avoid instabilities, we regularly de-alias the variables along the z direction according to the technique described in Hanasoge & Duvall (2007). Zero Dirichlet upper and lower boundary conditions are enforced at the outer edge of the C-PMLs for all the variables. A Fourier spectral method is applied in the horizontal directions in order to compute derivatives. Time stepping is achieved through the application of a second-order optimized Runge-Kutta scheme (Hu et al. 1996). The horizontal boundaries are chosen to be periodic while the vertical boundaries are required to be absorbent. We study this case in detail below; there is no loss of generality however, since the method can be extended in order to render the horizontal boundaries absorbent as well.

An inspection of Eqs. (18) to (27) reveals that 6 memory variables are needed: for v_z, the three (vector) Lorentz force components, and two in the induction equation. The memory and computational overheads are relatively small: (1) the boundary layer works well with 8–10 grid points; thus the six memory variables need only be stored over a small number of points; and (2) the additional computations are limited to a small number of multiplications and additions, as required for the evolution of the convolution Eqs. (19), (22), and (27).

3.2. Waves in a non-magnetic stratified fluid

We first demonstrate the ability of the C-PMLs at absorbing outgoing waves. The background medium is a stratified polytrope with index m = 2.15; the vertical extent of the computational domain is such that approximately 2.6 scale heights in density and 3.72 scale heights in pressure fill the domain (i.e. the pressure at the bottom is e^3.72 = 41.5 times the value at the top). The stratification properties of the polytrope are described in Appendix A.1. Waves are locally excited in a small region located at a distance of 18 Mm above the bottom boundary (the full vertical length is 68 Mm). For the C-PML parameters, we choose N = 2,R_c = 0.1%,f₀ = 0.005   Hz in Eqs. (29) and (31) for these calculations. Both the upper and lower C-PMLs are 10 grid points thick. In Fig. 1, we show snapshots of the wavefronts at four instants of time. The scaling in all the plots is identical. It is seen that the waves are almost completely absorbed by the upper boundary. In order to study the absorption more quantitatively, we plot the following energy invariant (summed over the entire grid) as a function of time in Fig. 2, (e.g., Bogdan et al. 1996): $\mathit{e}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left|}\mathrm{\right|}{v}\mathrm{|}{\mathrm{|}}^{\mathrm{2}}\mathrm{+}\frac{{\mathit{p}}^{\mathrm{2}}}{\mathrm{2}\mathit{\gamma }{\mathit{p}}_{\mathrm{0}}}\mathrm{·}$(33)

In order to test the stability of the C-PMLs, we integrate the hydrodynamic wave equations for around 300 wave periods (peak frequency  ~6 mHz; 12 h time integration). In this case, the setup is more realistic; the computational domain now straddles the solar photosphere, a region where the density and pressure drop exponentially with height. The waves propagate in a domain that contains 21 scale heights in density, representing a contrast of 1.3 billion between the bottom and top. For details of the stratification, see Appendix A.2.

We solve these equations in a computational domain of size 200 × 200 × 35 Mm³, (horizontal sides followed by vertical length) (1      Mm = 10⁶m), discretized using 256 × 256 × 300 grid points. Vertically, the box extends from approximately 34 Mm below the solar photosphere to 1 Mm above. Standing waves or normal modes are naturally created by the stratification of the medium; the modes are trapped between the surface and a specific depth below that depends on the frequency and wavenumber of the mode. A consequence of this trapping is that the mode has an evanescent tail that extends from the surface into the overlying atmosphere. We therefore place the upper computational boundary high enough above the surface so that the normal mode evanescent tails are largely unaffected by it.

An additional aspect is that of the wave excitation mechanism: in the Sun, waves are created by the action of the strong near-photospheric turbulence (e.g., Stein & Nordlund 2000). To mimic this process, we add a horizontally homogeneous vertically localized deterministic forcing function (S of the momentum equations).This source term is only vertically directed, i.e., S = S(x,y,t)δ(z − z_e)e_z, where z_e =  − 50 km. In order to generate S, we use a Gaussian random number generator to populate its Fourier domain and the resultant function is multiplied by the average frequency power spectrum of the Sun, modeled as a Gaussian with peak frequency at 3 mHz and a full width of 1 mHz. The inverse transform of this function is S(x,y,t) (e.g., Hanasoge et al. 2006; Hanasoge 2008).

The calculation is stable over the entire time integration. The wavenumber-frequency power spectrum of the vertical velocity component v_z at a height of z = 200 km above the photosphere is plotted in Fig. 3. The horizontal axis represents the non-dimensional horizontal wavenumber k_h   R_⊙, where R_⊙ is the solar radius and k_h is the horizontal wavenumber, while the vertical axis is frequency in mHz.

In order to study the impact of these absorbing layers, we perform simulations with C-PMLs placed adjacent to the upper and lower boundaries (Fig. 3) and damping sponges (Fig. 4). A demonstration of the efficacy of the C-PML is the absence of artifacts related to the lower boundary in the modal power spectrum in Fig. 3. Weak reflections from the lower boundary cause the dispersion relation to change, manifested in the power spectrum as a flattening of the curvature of the ridges, as seen in Fig. 4.

3.3. Stratified MHD fluid

We first perform a test akin to that of Fig. 2. A constant inclined magnetic field is embedded in the polytrope of Appendix A.1; the strength of the field is such that the Alfvén speed is approximately four times the sound speed at the upper boundary. Waves are excited in the vicinity of the upper boundary by horizontally shaking the field lines. The fast and Alfvén waves reach the upper boundary first but the slow modes follow behind, which makes the demonstration of the absorptive properties of the boundary formulation less evident than in Fig. 2. Furthermore, constructing energy invariants for MHD waves in stratified media is not a trivial affair due to the fairly complicated Lorentz force terms. A number of authors have studied this problem (e.g., Bray & Loughhead 1974; Parker 1979; Leroy 1985) and have arrived at differing versions of an invariant (personal communication, P. S. Cally 2009); we use the following approximate form (kinetic + thermal + magnetic energies; e.g., Bray & Loughhead 1974): $\mathit{e}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left|}\mathrm{\right|}{v}\mathrm{|}{\mathrm{|}}^{\mathrm{2}}\mathrm{+}\frac{{\mathit{p}}^{\mathrm{2}}}{\mathrm{2}\mathit{\gamma }{\mathit{p}}_{\mathrm{0}}}\mathrm{+}\frac{\mathrm{\left|}\mathrm{\right|}{B}\mathrm{|}{\mathrm{|}}^{\mathrm{2}}}{\mathrm{8}\mathit{\pi }}\mathrm{·}$(34)We plot the time history of the energy summed over the entire computational domain in Fig. 5.

Our final test consists of studying the long-term stability of the method numerically. A magnetic flux tube (e.g., Moradi et al. 2009) is embedded in the stratified polytrope (Appendix A.2); waves are excited in the non-magnetic regions and allowed to propagate through the flux tube, whose geometry is shown in Fig. 6. The system is stable up to 12 h of time integration. The modal power spectrum from this data looks almost identical to that in Fig. 3 (with the exception of small frequency shifts) and is hence not shown. The maintenance of ∇·B = 0 numerically is an essential task of the scheme (e.g., Tóth 2000). Discretization errors are one of the primary causes for non-zero ∇·B = 0; additionally, in this case, the boundary formulation in non ∇·B = 0 conserving. Consequently, it is important to test if the error grows with time and to quantify its level. In order to measure the normalized errors in maintaining ∇·B = 0, we compute the following quantity: $\begin{array}{ccc}{\mathit{n}}_{\mathit{e}}& \mathrm{=}& \end{array}$(35)where z₁,z₂ are the limits of the physically-relevant domain (i.e., the non-absorbing regions) and the following definition for the L₂ norm is applied: ${\mathit{L}}_{\mathrm{2}}\mathrm{\left[}\mathit{f}\mathrm{\right(}\mathit{x}\mathrm{\left)}\mathrm{\right]}\mathrm{=}\sqrt{\sum _{\mathit{i}}\mathit{f}\mathrm{(}{\mathit{x}}_{\mathit{i}}{\mathrm{)}}^{\mathrm{2}}}\mathit{,}$(36)where the sum is over all the grid points x_i in the computational domain. We display the time evolution of n_e on the left panel of Fig. 7. One can see that n_e is nearly constant, perhaps showing a slight decrease with time. On the right panel we show the time history of the numerator of n_e (Eq. (35)). The gradual increase in the magnetic energy is due to constant input wave excitation via the S term in the vertical momentum equation. The growth rates of the magnetic energy demonstrate the numerical stability of the system over the period of time integration.

4. Conclusions and future work

We have developed a Convolutional Perfectly Matched Layer (C-PML) for the stratified linearized Euler equations and a highly efficient absorption layer for the ideal MHD equations. This boundary formulation is quite useful for the calculation of wave propagation in astrophysical plasmas, where stratification and magnetic fields abound. The absorption layer for MHD waves is a slightly altered version of the convolutional formulation of (e.g., Roden & Gedney 2000), requiring an extra sponge-like term in order to stabilize the system. While the boundary methods as discussed here are perfectly matched for the stratified Euler equations, the inclusion of the sponge term casts some doubt as to whether the same could be said of the MHD equations. Simulations with the absorption layers for the stratified MHD and Euler equations were performed and found to be stable over a time length of 300 wave periods. With increasingly sophisticated boundary formulations, such as that of Hu et al. (2008), who developed the PML for turbulent flows, we may in the near future be able to extend these techniques to fully non-linear MHD turbulence/dynamo calculations.

Acknowledgments

S.M.H. is supported by the German Aerospace Center (DLR) through the grant “German Data Center for SDO”. The study on sunspot modeling is supported by the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement #210949, “Seismic Imaging of the Solar Interior”, to PI L. Gizon.

References

Appendix A: Stratification

A.1. Only polytrope

For the absorption tests in Sects. 3.2 and 3.3, we use the following polytropic stratification prescription: $\begin{array}{ccc}{\mathit{p}}_{\mathrm{0}}\mathrm{\left(}\mathit{z}\mathrm{\right)}& \mathrm{=}& {\mathit{p}}_{\mathrm{poly}}{\left(\mathrm{1}\mathrm{-}\frac{\mathit{z}}{{\mathit{z}}_{\mathit{f}}}\right)}^{\mathit{m}\mathrm{+}\mathrm{1}}\mathit{,}\\ {\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{z}\mathrm{\right)}& \mathrm{=}& {\mathit{\rho }}_{\mathrm{poly}}{\left(\mathrm{1}\mathrm{-}\frac{\mathit{z}}{{\mathit{z}}_{\mathit{f}}}\right)}^{\mathit{m}\mathrm{+}\mathrm{1}}\mathit{,}\\ {\mathit{g}}_{\mathrm{0}}& \mathrm{=}& \\ {\mathrm{\Gamma }}_{\mathrm{1}}& \mathrm{=}& \mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathit{m}}\mathit{,}\\ {\mathit{c}}_{\mathrm{0}}\mathrm{\left(}\mathit{z}\mathrm{\right)}& \mathrm{=}& \sqrt{\frac{{\mathrm{\Gamma }}_{\mathrm{1}}{\mathit{p}}_{\mathrm{poly}}}{{\mathit{\rho }}_{\mathrm{poly}}}\left(\mathrm{1}\mathrm{-}\frac{\mathit{z}}{{\mathit{z}}_{\mathit{f}}}\right)}\mathrm{·}\end{array}$We set p_poly = 1.178 × 10⁵   dynes   cm^-2,ρ_poly = 3.093 × 10^-7   g   cm^-3,z_f =  − 0.450      Mm, and m = 2.150. Note also that z is the height, i.e. the atmospheric layers are described by z > 0 and vice versa; z = 0 is the fiducial surface. The polytrope is truncated at z =  − 31.77 Mm (upper boundary) and z =  − 104.40 Mm (lower boundary) respectively.

A.2. Polytrope + isothermal layer

Waves propagating toward the surface in the Sun are reflected by a sharp fluctuation in the background density gradient (e.g., Christensen–Dalsgaard 2003). This reflection zone is located at a height of z_r ~  −0.050 Mm below the surface. We mimic this by patching the polytropic stratification of Sect. A.1 with an overlying isothermal layer. This patching results in a fluctuation in the density gradient, leading to an acoustic cut-off frequency of approximately 5.4 mHz, similar to the solar value. For z < z_M, we use the prescription of Sect. A.1. For z ≥ z_r, we use the following equations: $\begin{array}{ccc}{\mathit{p}}_{\mathrm{0}}\mathrm{\left(}\mathit{z}\mathrm{\right)}& \mathrm{=}& {\mathit{p}}_{\mathrm{iso}}\exp \left[\frac{{\mathit{z}}_{\mathrm{r}}\mathrm{-}\mathit{z}}{\mathit{H}}\right]\mathit{,}\\ {\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{z}\mathrm{\right)}& \mathrm{=}& {\mathit{\rho }}_{\mathrm{iso}}\exp \left[\frac{{\mathit{z}}_{\mathrm{r}}\mathrm{-}\mathit{z}}{\mathit{H}}\right]\mathit{,}\\ {\mathit{\rho }}_{\mathrm{iso}}& \mathrm{=}& {\mathit{\rho }}_{\mathrm{poly}}{\left(\mathrm{1}\mathrm{-}\frac{{\mathit{z}}_{\mathrm{r}}}{{\mathit{z}}_{\mathit{f}}}\right)}^{\mathit{m}}\mathit{,}\\ {\mathit{p}}_{\mathrm{iso}}& \mathrm{=}& {\mathit{p}}_{\mathrm{poly}}{\left(\mathrm{1}\mathrm{-}\frac{{\mathit{z}}_{\mathrm{r}}}{{\mathit{z}}_{\mathit{f}}}\right)}^{\mathit{m}\mathrm{+}\mathrm{1}}\mathit{,}\\ \mathit{H}& \mathrm{=}& \frac{{\mathit{p}}_{\mathrm{iso}}}{{\mathit{g}}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{iso}}}\mathrm{·}\end{array}$(A.3)The relations for ρ_iso and p_iso arise from the requirement of continuity of the pressure and density at the matching point between the polytropic and isothermal layers. The relation for H is a consequence of enforcing hydrostatic balance. This model is truncated at z =  −34 Mm (lower boundary) and z =  + 1 Mm (upper boundary).

All Figures