ENERGY LOSS OF SOLAR p MODES DUE TO THE EXCITATION OF MAGNETIC SAUSAGE TUBE WAVES: IMPORTANCE OF COUPLING THE UPPER ATMOSPHERE

The turbulent motions in the upper convection zone of the Sun are thought to be the source of acoustic energy, stochastically exciting p modes. The convection zone is also the region where many thin magnetic fibril fields reside, continually shaken by these modes.

Solar p modes have been studied extensively using high-quality data gathered by ground- and spaced-based telescopes, namely at the six BiSON sites (Chaplin 1997), the GONG network (Hill et al. 1996), the GOLF experiment, and MDI instrument on board Solar and Heliospheric Observatory (Lazrek et al. 1997; Scherrer et al. 1995), and, more recently, Helioseismic and Magnetic Imager on board Solar Dynamics Observatory (Schou et al. 2012). We now have a clear and consistent picture of the observed linewidths (damping rates) associated with f and p modes (Chaplin et al. 1997; Lazrek et al. 1997; Komm et al. 2000a, 2000b, 2001). The most notable feature of the observed damping rates is the "plateau" between 2.5 and 3.1 mHz followed by a gradual rise. Damping rates indicate a global observational signature of the energy loss experienced by the p modes while a local signature of this energy loss is characterized by an absorption coefficient which is a measure of the difference between the incoming and outgoing p-mode power around a magnetic region. Magnetic features observed within the solar photosphere such as sunspots, plages, and pores are efficient absorbers of p-mode energy. The absorption coefficient was first measured by a Fourier–Hankel decomposition method (Braun et al. 1988; Braun 1995) and in recent years by helioseismic holography (Braun & Birch 2008). The frequency dependence of the absorption coefficient has a distinctive bell curve structure which is consistent over a variety of magnetic regions and p-mode orders, with a peak between 3 and 4 mHz. The saturation of the absorption coefficient seen above 4 mHz is likely a result of the acoustic emission of p-mode power observed around magnetic features at higher frequencies (see, e.g., Jain & Haber 2002).

There have been several theoretical studies to explain the observed frequency dependence of damping rates and absorption coefficients. Bogdan et al. (1996) developed a model whereby p-mode energy is lost to magnetic flux tubes via the excitation of MHD tube waves by p-mode buffeting. They found that tube mode excitation on fibril magnetic fields may be a dominant source of low-frequency p-mode damping rates. Further work on this mechanism of p-mode damping by Hindman & Jain (2008) found that the damping rates are very sensitive to the upper boundary condition. They also demonstrated that p modes efficiently generate tube waves and that an energy flux in excess of 10⁵ erg cm⁻² s⁻¹ can be driven upward through photospheric levels. Jain et al. (2009) presented a method for calculating the absorption coefficient for this mechanism for a single magnetic flux tube and for a collection of identical tubes to simulate a plage region observed by Braun & Birch (2008). This work was extended to incorporate a fixed and random distribution of tubes with varying plasma properties, thus enabling us to measure and estimate the statistical properties of flux tubes within an observed plage region (see Gascoyne et al. 2011; Jain et al. 2011a, 2011b for details).

All the abovementioned studies modeled the magnetic fibrils and their surroundings with a truncated polytrope. Regardless of whether a stress-free (Bogdan et al. 1996) or a maximal-flux (Hindman & Jain 2008) boundary condition was assumed at the truncation height (see also Crouch & Cally 1999 for a radiation boundary condition), the models did not satisfactorily represent wave reflection from the upper photosphere/lower chromosphere. Also, these studies revealed that the magnitude of damping rate, Γ, and absorption coefficient, α, were sensitive to the boundary conditions used. Here, we calculate both damping rates and absorption coefficients for the equilibrium model where the magnetic fibril is extended into an isothermal region above a polytrope. We demonstrate that the coupling of sausage waves and external p modes to the overlying atmosphere is necessary to avoid the sensitivity of boundary conditions applied at the model photosphere (z = −z₀).

The paper is organized as follows. In Section 2, we describe the equilibrium atmosphere of our two-region model and the static configuration of the magnetic flux tubes. Section 3 details the wavefield in the non-magnetic surroundings for both the polytropic and isothermal atmospheres separately and the coupling conditions across the interface. The wavefields describe the propagating (in the polytrope) and evanescent (in the isothermal region) form of the f- and p-mode oscillations. In Section 4, we derive the driven sausage waves of the thin magnetic flux tube and calculate the energy flux of sausage waves at the bottom and top of the polytropic and isothermal atmosphere respectively. In Section 5, we present the damping rates and absorption coefficients of p modes, and discuss our results.

In this section, we describe the equilibrium model of a slender magnetic fibril as an axisymmetric magnetic flux tube embedded in a plane-parallel, gravitationally stratified two-region atmosphere. The lower and upper layers of the model represent the subsurface and the upper atmosphere, respectively. It is therefore appropriate to consider polytropic and isothermal atmospheres for the interior and the upper photosphere, respectively. A schematic representation of the excitation of tube waves along a thin magnetic fibril due to p-mode buffeting in our two-region model is shown in Figure 1.

Zoom In Zoom Out Reset image size

The surrounding non-magnetic medium is modeled as a polytropic atmosphere below the isothermal atmosphere. The vertical variation of the pressure p_e and density ρ_e for the polytropic atmosphere is

$$\begin{equation} p_{e}=p_{0}\left(-\frac{z}{z_{0}}\right)^{a+1},\quad \rho _{e}=\rho _{0}\left(-\frac{z}{z_{0}}\right)^{a}, \end{equation} \tag{ 1 }$$

where ρ₀ is the density at the interface z = −z₀ and p₀ = gz₀ρ₀/(a + 1). The isothermal atmosphere, which represents the upper photosphere/lower chromosphere of the Sun, has its lower boundary at the truncation depth z = −z₀ and extends infinitely upward as z → ∞. The vertical variation of the pressure p_e and density ρ_e in the isothermal atmosphere is

$$\begin{equation} p_{e}=p_{0}e^{-(z+z_{0})/H},\quad \rho _{e}=\rho _{0}e^{-(z+z_{0})/H}, \end{equation} \tag{ 2 }$$

where the scale height H = z₀/(a + 1). The sound speed for both atmospheres is c² = γp_e/ρ_e.

It is clear from Equations (1) and (2) that the equilibrium quantities are continuous through the interface. In the polytropic atmosphere, the pressure, density, and sound speed vary as a power law where a is the polytropic index and is related to the ratio of specific heats γ via γ = (a + 1)/a. In the isothermal atmosphere, the pressure and density decrease exponentially with height at the same rate, thus the sound speed (and the temperature) is constant.

The characteristic physical scales for our model photosphere (z = −z₀) are assumed to coincide with the photospheric reference model of Maltby et al. (1986), and so we adopt ρ₀ = 2.78 × 10⁻⁷ g cm⁻³, p₀ = 1.21 × 10⁵ g cm⁻¹ s⁻², and g = 2.775 × 10⁴ cm s⁻². The choice of polytropic index a = 1.5 yields the truncation depth z₀ = 392 km and a photospheric sound speed of 8.52 km s⁻¹.

Now we define the profile of our magnetic fibril field, first introduced by Bogdan et al. (1996) (see also Hasan 1997; Hindman & Jain 2008; Jain et al. 2009; Gascoyne et al. 2011). We assume the magnetic fibril is a thin, untwisted, axisymmetric, and vertically aligned magnetic flux tube which lacks internal lateral structure. For such a thin magnetic flux tube to be in thermal equilibrium and hydrostatic balance requires that temperature and the total pressure is uniform across the tube and is equal to the external pressure and temperature. Thus, the magnetic pressure has the same scale height as the gas pressure which results in a constant plasma-β (ratio of gas to magnetic pressure) with height within the tube. Since the gas pressure decreases rapidly with height in the solar atmosphere, so does the magnetic pressure; this then results in significant flaring of field lines which inevitably will break the thin flux tube approximation. However, if the skin depth of the wave solutions are sufficiently small, the waves will not sense the uppermost portion of the atmosphere where the approximation becomes invalid. Therefore, we can safely ignore this problem if we restrict our attention to frequencies below the cut off frequency for the isothermal atmosphere. The lateral variation of the magnetic field is ignored (see Bogdan et al. 1996 for details). Thus, the tube's internal pressure P(z), density ρ(z), magnetic field B(z), and cross-sectional area A(z) can be described as follows:

$$\begin{equation} p(z)=\frac{\beta }{\beta +1}p_{e}(z), \end{equation} \tag{ 3 }$$

$$\begin{equation} \rho (z)=\frac{\beta }{\beta +1}\rho _{e}(z), \end{equation} \tag{ 4 }$$

$$\begin{equation} \frac{B^2(z)}{8\pi }=\frac{1}{\beta +1}p_{e}(z), \end{equation} \tag{ 5 }$$

$$\begin{equation} A(z)=\frac{\Theta }{B(z)}=\left(\frac{\beta +1}{8\pi p_{e}(z)}\right)^{1/2}\Theta, \end{equation} \tag{ 6 }$$

where Θ is the constant magnetic flux. All parameters associated with the external non-magnetic atmosphere which have an internal counterpart associated with the magnetic flux tube will be distinguished with the subscript "e."

We assume that the fluid displacement takes the form

$$\begin{equation} \boldsymbol {\xi }_{e}=\tilde{\xi }_{e}(z)e^{i(kx-\omega t)}, \end{equation} \tag{ 7 }$$

where $\tilde{\xi }_{e}(z)$ is the vertical eigenfunction of the displacement vector. For the polytrope, the temperature varies linearly with depth (∝ − gz/a), thus the vertical eigenfunction is represented by the Whittaker function $W_{\kappa,\mu }$ (see Equation (A3)). However, for the isothermal region the temperature is constant and so we also Fourier analyze the governing equation (Equation (A6)) in height. Thus, the vertical wavenumber in the isothermal region is defined as

$$\begin{equation} k^{\prime }_{z}=\sqrt{\left(\omega ^2-\omega ^2_{{\rm ac}}\right)\frac{1}{c^2}-\left(\omega ^2-\omega ^2_{{\rm BV}}\right)\frac{k_{x}^2}{\omega ^2}}, \end{equation} \tag{ 8 }$$

where ω_ac = γg/(2c) and $\omega _{{\rm BV}}=\sqrt{\gamma -1}g/c$ are the acoustic cut-off and Brunt Väisälä frequencies, respectively. We match the wave perturbations in the two atmospheres by ensuring continuity of the vertical displacement ξ_{z, e} and Lagrangian pressure perturbation δp_e across the interface z = −z₀. The upper and lower boundary conditions are that the wave energy density vanishes as z → ±∞. The details of this procedure is described in Appendix A.

In order to satisfy the matching conditions, we solve Equation (A8) for κ, at each dimensionless frequency ν. We obtain a discrete set of eigenvalues which correspond to each individual wave mode (Figure 2) similarly to the truncated (Bogdan et al. 1996) and complete (Fan et al. 1995) polytrope cases. The fundamental mode (n = 0) is the same for all three cases which is expected since this mode is a surface gravity wave, i.e., incompressible ($\mathbf {\nabla }\,{\cdot}\, \boldsymbol {\xi _{e}}=0$) in both the polytropic and isothermal atmospheres. The isothermal atmosphere sets a global frequency cut-off when $k^{\prime }_{z}=0$ (see Equation (8)). For a polytropic index of a = 1.5 and an interface located at z = −z₀ where z₀ = 392 km, the cut-off frequency is ≈4.3 mHz (Figure 2), below this frequency the waves are evanescent whereas above the cut-off the waves are propagating in the isothermal region. In this paper, we are concerned with trapped p modes, thus we calculate eigenvalues and eigenfunctions for the two-region model for frequencies ω/(2π) ⩽ 4.3 mHz. The solid curve in Figure 3 illustrates the vertical displacement for the p₂-mode at 3 mHz, clearly showing the propagating portion of the wave in the polytropic region and the evanescent tail in the isothermal region.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The incident f and p modes, described by Equations (A10) and (A11), drive waves along the thin magnetic tube, thus removing energy from the p-mode cavity. Here we shall concentrate solely on the excitation of sausage waves which can be described as an axisymmetric-pressure-driven wave with motions parallel to the magnetic flux tube. Under the thin flux tube approximation we derive a set of MHD equations for an inviscid and infinitely conducting medium. On linearizing these equations and assuming a time dependence consistent with the external perturbations, we can describe the interaction of the field-free medium with the thin flux tube for the sausage mode in both the polytropic and isothermal atmosphere. Thus, using the formulation of Hindman & Jain (2008), Jain et al. (2009), and Gascoyne et al. (2011) for the polytropic atmosphere and Hasan (1997) for the isothermal atmosphere, we obtain the governing equations for the propagation of sausage waves within our prescribed two-region atmosphere,

$$\begin{equation*} \left[c^2_{{\rm T}}\frac{\partial ^2}{\partial z^2}-\frac{\gamma g}{2}\frac{c^2_{{\rm T}}}{c^2}\frac{\partial }{\partial z}+\omega ^2+N_{g}\right]\xi _{z} \end{equation*}$$

$$\begin{equation} =\frac{c^2_{{\rm T}}}{\rho V^2_{{\rm A}}}\left(\frac{\partial }{\partial z}+\frac{g}{c^2}\right)p_{e}^{\prime }, \end{equation} \tag{ 9 }$$

where $V^2_{{\rm A}}=B^2/(4\pi \rho)$ and $c^2_{{\rm T}}=c^2/(1+c^2/V^2_{{\rm A}})$ are the Alfvén speed and tube speed, respectively. The governing equations for the propagation of sausage waves in the two regions are essentially the same besides the term N_g ($=-\omega _{{\rm BV}}^2\left(c^2+\frac{\gamma }{2}V^2_{{\rm A}}\right)/(c^2+V^2_{{\rm A}})$), which permits the propagation of gravity waves in the isothermal layer; this quantity is zero for the polytrope which is in adiabatic equilibrium.

The solutions to Equation (9) are derived in detail in Appendix B and take the form

$$\begin{equation} \boldsymbol {\xi }_{\Vert }=\xi _{\Vert }(z)e^{-i\omega t}\hat{z}, \end{equation} \tag{ 10 }$$

where ξ_||(z) is written explicitly in Equations (B2) and (B3) for the isothermal and polytropic atmospheres, respectively. The solution in the polytropic atmosphere is a combination of Hankel functions and Greens functions, whereas in the isothermal region the solution can be written in terms of Fourier components, thus the vertical wavenumber for the sausage wave in the isothermal atmosphere can be written as

$$\begin{equation} k_{\Vert }=\sqrt{\frac{\omega ^2-\omega _{s}^2}{c^2_{{\rm T}}}},\quad \omega _{s}^2=\omega _{{\rm BV}}^2+\frac{c^2_{{\rm T}}}{H^2}\left(\frac{3}{4}-\frac{1}{\gamma }\right)^2. \end{equation} \tag{ 11 }$$

The matching conditions for the sausage tube wave, applied at z = −z₀, are analogous to the field-free medium in that the continuity of the vertical displacement ξ_|| and Lagrangian pressure perturbation δp_|| across the interface between the two regions are ensured (discussed in detail in Appendix B.1). The upper and lower boundary conditions we require are that the wave energy density is finite as z → ±∞. The details are given in Appendix B.

The aim is to understand how MHD tube waves can transport energy from p modes in the convection zone down into the deep interior and up into the upper photosphere/lower chromosphere. In Appendix C, we derive expressions for the upward, $E^{(u)}_{\Vert }$, and downward, $E^{(d)}_{\Vert }$, energy flux of sausage waves driven by f and p modes in a two-region, polytropic, and isothermal atmosphere. We quantify the energy loss of f and p modes as a damping rate, Γ, and absorption coefficient, α, which are explicitly calculated in Appendix C.1 from the wave energy flux of sausage waves and the energy contained in f and p modes for the given model.

In Figures 4 and 5, we present damping rates and absorption coefficients for three values of the plasma-β = 0.1, 1, 10. The β dependence of both the damping and of absorption vary non-linearly, with a peak at around β = 1 and decreasing in magnitude with increasing and decreasing β from this peak. The upward and downward wave energy flux are calculated as z → ∞ and z → −∞, respectively. No wave energy escapes the upper domain into the overlying atmosphere since the sausage waves are evanescent in the frequency regime to which we are restricted (see Figure 2 below the dot–dashed line). Therefore, all wave energy that escapes the acoustic cavity in the model is through the lower domain. This property is similar to that of a model with a stress-free condition imposed at the photosphere, except that additional driving can occur above the photosphere and the phase shift of the reflected wave may differ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The low-frequency portion of the damping rates and absorption coefficients are in very good agreement with previous studies where a stress-free boundary condition was applied at the model photosphere (z = −z₀) (see Hindman & Jain 2008; Jain et al. 2009). In this frequency regime, the upper turning point of both sausage waves and p modes are well below the interface, thus the isothermal region behaves very similar to a stress-free boundary at z = −z₀ (see Figure 6). This also implies that for low frequencies most of the p-mode energy lost to the tube escapes down the tube and very little escapes into the isothermal atmosphere, which is apparent when the energy fluxes $E_{\Vert }^{(u)}$ and $E_{\Vert }^{(d)}$ are investigated in detail (see Appendix C). In Figure 7, we plot the upward sausage wave energy flux, $E_{\Vert }^{(u)}$ (Equation (C4)), as a function of height and the corresponding value of $E_{\Vert }^{(d)}$ as a reference for how large the upward energy flux is. The solid curve in Figure 7 shows that the upward wave energy flux decays with height, but the rate of decay decreases with frequency. For frequencies above 3 mHz, a significant amount of energy flux remains in the tube, even a few scale heights above the interface, this energy flux eventually decays to zero and does not leave the domain, i.e., does not contribute to the damping and absorption of p modes. The sausage waves are evanescent in the isothermal region, thus their wave energy flux is not being carried by free oscillations. The flux is the result of driven oscillations from below and the flux decreases with height because the p-mode driver (also evanescent) acts as a sink in the isothermal region removing energy form the tube. However, caution must be taken for frequencies too close to the cut-off frequency, since very close to this frequency the p-mode eigenfunction and, hence the wave driver, has a skin depth that is insufficiently short to ensure that the excitation remains finite with height. This manifests as divergent behavior in $E_{\Vert }^{(u)}$ in the limit z → ∞. Explicitly, the term $\mathcal {S}$ in Equation (C4) diverges with height because of its dependence on ∼exp (1/2[iλ_{z, e} + iλ_|| + (a + 1)/2]). This divergent behavior is clearly unphysical and can be avoided if we ensure that $|k_{z}^{\prime }|+|k_{\Vert }|>({1}/{4H})$, or for the parameters used here, if we restrict our attention to frequencies below 4.23 mHz.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The inclusion of an upper isothermal atmosphere allows us to take into account the coupling of all waves between the two regions which is an arduous procedure but required in order to model the reflective properties of the upper atmosphere with even minimal fidelity. In Figure 6, we compare damping rates for our two-region (infinite domain) model with truncated polytrope models with stress-free and maximal-flux boundary conditions applied at the interface (z = −z₀) (Bogdan et al. 1996; Hindman & Jain 2008; Jain et al. 2009). As noted earlier, the low-frequency behavior of the two-region model is similar to that of the stress-free case; the isothermal atmosphere reflects all waves back down the tube in a similar manner as in the stress-free case. For higher frequencies, the damping of p modes decrease compared to the stress-free case; this is due to the mechanism discussed earlier, that the p modes act as a sink in the isothermal atmosphere extracting energy out of sausage waves. In the stress-free model, the upward propagating wave is forced to reflect at the interface and thus propagates down the tube and this energy loss is measured at the lower boundary as a combination of downward and upward reflected propagating waves. In the present two-region case, this upward propagating wave in the polytrope is no longer reflecting in the same manner at the interface and instead becomes evanescent very close to the interface and a portion of its wave energy is transported into the isothermal region. This energy flux then decreases to zero with height and thus we measure a zero contribution to the escaping wave energy flux, reducing the damping rate at higher frequencies where this effect becomes increasingly dominant. The maximal-flux condition (introduced in Hindman & Jain 2008) is vital in providing an upper limit on the amount of energy lost from f and p modes to magnetic tube waves. However, it does not provide insight into the nature of the interaction of tube waves with a realistic overlying atmosphere (Figure 6, dashed lines).

5.1. Conclusion

This work is supported by STFC (UK). B.W.H. also acknowledges NASA grants NNX09AB04G, NNX14AC05G, and NNZ14AG05.

In this section, we shall derive the wavefields for both the polytropic and isothermal atmospheres. At the interface z = −z₀ the vertical component of the displacement perturbation ξ_{z, e}, and the Lagrangian pressure perturbation $\delta p_{e}=-\gamma p_{e}\mathbf {\nabla }\cdot \boldsymbol {\xi _{e}}$ are both (assumed to be) continuous.

A.1. Polytropic Wavefield

In the field-free polytropic atmosphere we follow the procedure adopted by Bogdan et al. (1996), Hindman & Jain (2008), and Gascoyne et al. (2011). Thus, the wave perturbations $p_{e}^{\prime }$, $\rho _{e}^{\prime }$, $\boldsymbol {\xi _{e}}$, in pressure, density, and fluid displacement, respectively, can be written in terms of a displacement potential Φ,

$$\begin{equation} \boldsymbol {\xi _{e}}=\mathbf {\nabla }\Phi,p_{e}^{\prime }=-\rho _{e}\frac{\partial ^2\Phi }{\partial t^2},\rho _{e}^{\prime }=-\frac{\rho _{e}}{c^2}\frac{\partial ^2\Phi }{\partial t^2}, \end{equation} \tag{ A1 }$$

which together form the wave equation (see Lamb 1945),

$$\begin{equation} \frac{\partial ^2\Phi }{\partial t^2}=c^2\mathbf {\nabla }^2\Phi -g\frac{\partial \Phi }{\partial z}. \end{equation} \tag{ A2 }$$

The incident acoustic wavefield expressed by Equation (A2) supports plane wave solutions of the form

$$\begin{equation} \Phi (\boldsymbol {x},t)=z_{0}\mathcal {A}e^{i(k_{x}x-\omega t)}Q(z), \end{equation} \tag{ A3 }$$

where $\mathcal {A}$ is the complex wave amplitude, ω the temporal frequency, k_x the horizontal wavenumber, and Q(z) the vertical eigenfunction. Thus, (A2) yields

$$\begin{equation} Q(w)=w^{-(\mu +1/2)}W_{\kappa,\mu }(w), \end{equation} \tag{ A4 }$$

where $W_{\kappa,\mu }$ is the Whittaker function of the second kind. Also,

$$\begin{eqnarray} &&w=-2k_{x}z,\quad \nu ^2=\frac{a\omega ^2z_{0}}{g},\quad \kappa =\frac{\nu ^2}{2k_{x}z_{0}},\nonumber\\ && \mu =\frac{(a-1)}{2},\quad \lambda =2k_{x}z_{0}. \end{eqnarray} \tag{ A5 }$$

A.2. Isothermal Wavefield

In the field-free isothermal atmosphere, the linearized hydrostatic equations can be manipulated to form the governing wave equation for the vertical displacement perturbation as

$$\begin{equation} \frac{\partial ^2\boldsymbol {\xi _{e}}}{\partial t^2}=c^2\mathbf {\nabla }\left(\mathbf {\nabla }\cdot \boldsymbol {\xi _{e}}\right)+(\gamma -1)\boldsymbol {g}\left(\mathbf {\nabla }\cdot \boldsymbol {\xi _{e}}\right)+\mathbf {\nabla }\left(\boldsymbol {g}\cdot \boldsymbol {\xi _{e}}\right). \end{equation} \tag{ A6 }$$

Since this partial differential equation has constant coefficients, we Fourier analyze in time and space, i.e., $\boldsymbol {\xi _{e}}=\boldsymbol {A}\exp \left[i\left(\boldsymbol {k}\cdot \boldsymbol {x}-\omega t\right)\right]$. If we refer to Figure 1, the isothermal atmosphere is situated on top of the polytropic atmosphere with its base at z = −z₀, thus substituting z = z + z₀ and denoting $k_{z}=-i({1}/{2H})+k_{z}^{\prime }$ (see Equation (8)) the displacement $\boldsymbol {\xi _{e}}$ can be written as

$$\begin{equation} \boldsymbol {\xi _{e}}=\boldsymbol {A}e^{(z+z_{0})/(2H)}e^{i\left(k_{z}^{\prime }(z+z_{0})+k_{x}x-\omega t\right)}. \end{equation} \tag{ A7 }$$

A.3. Matching Conditions across the Interface

The continuity of Lagrangian pressure perturbation at the interface yields the following relation

$$\begin{eqnarray} && [\nu ^2a-2\kappa ^2(i\lambda _{z,e}+a+1)]W_{\kappa,\mu }(\lambda)\nonumber\\ &&\qquad+\,\kappa (2\kappa -a)(i\lambda _{z,e}-a+1)W_{\kappa,\mu +1}(\lambda)=0, \end{eqnarray} \tag{ A8 }$$

where

$$\begin{equation} \lambda _{z,e}=2k_{z}^{\prime }z_{0}=\sqrt{\frac{\lambda }{\kappa }(4\kappa ^2+a-\nu ^2)-(a+1)^2}. \end{equation} \tag{ A9 }$$

Note that in order to solve Equation (A8), we require that λ_{z, e} is purely imaginary, thus we only consider trapped p-mode solutions. We couple the propagating waves in the polytrope on to evanescent waves in the isothermal atmosphere, thus we only consider the frequencies below the acoustic cut-off frequency which is given by $\lambda _{z,e}^2<0$.

After solving Equation (A8) for the eigenfrequencies, using the non-dimensional parameters (A5) the vertical component of the eigenfunctions can be calculated as

$$\begin{equation} \xi ^{{\rm (I)}}_{z,e}=-z_{0}\frac{\mathcal {A}}{z_{0}}\lambda \left.\frac{dQ}{dw}\right|_{w=\lambda }e^{\frac{1}{2}(i\lambda _{z,e}+a+1)\left(1-\frac{w}{\lambda }\right)},\quad w\le \lambda, \end{equation} \tag{ A10 }$$

$$\begin{equation} \xi ^{{\rm (P)}}_{z,e}=-z_{0}\frac{\mathcal {A}}{z_{0}}\lambda \frac{dQ}{dw},\quad w\ge \lambda, \end{equation} \tag{ A11 }$$

where we have used the notation $\xi ^{(I)}_{z,e}$ and $\xi ^{(P)}_{z,e}$ for the isothermal and polytropic vertical displacements, respectively. Above we have omitted the Fourier components e^{i(kx − ωt)}, which we shall continue to do so for brevity but still take into account the t and x derivatives when encountered.

A.4. Energy Contained in f- and p-mode Oscillations

The wave energy density of f and p modes for the prescribed model is given by (see Bray & Loughhead 1974)

$$\begin{equation} \mathcal {W}=\frac{1}{4}\rho _{e}\left|\frac{\partial \boldsymbol {\xi }_{e}}{\partial t}\right|^2+\frac{|p^{\prime }_{e}|^2}{4\rho _{e}c^2}-\frac{1}{4}\frac{g^2}{c^2}\rho _{e}\left(1+\frac{c^2}{g\rho _{e}}\frac{d\rho _{e}}{dz}\right)|\xi _{z,e}|^2. \end{equation} \tag{ A12 }$$

For a neutrally stable atmosphere, the last term, which corresponds to the gravitational energy, vanishes. Taking the time average, integrating over depth and multiplying by $4\pi R^2_{\odot }$ gives the energy contained in f- and p-mode oscillations which we shall denote E_n. The polytropic segment of this integral is given by Hindman & Jain (2008),

$$\begin{equation} E^{(P)}_{n}=4\pi R^2_{\odot }\frac{g\rho _{0}z_{0}^2}{4a}\lambda ^2\frac{|\mathcal {A}|^2}{z_{0}^2}\mathcal {N}_{n}^{(1)}, \end{equation} \tag{ A13 }$$

where

$$\begin{equation} \mathcal {N}_{n}^{(1)}=\frac{\kappa ^2}{\lambda ^a}\int ^{\infty }_{\lambda }dw\frac{w^a}{\kappa }\left[\left(\frac{dQ(w)}{dw}\right)^2+\left(\frac{1}{4}+\frac{\kappa }{w}\right)Q^2(w)\right]. \end{equation} \tag{ A14 }$$

The isothermal contribution to the wave energy density is calculated as (see also Hasan 1997)

$$\begin{equation} \mathcal {W}^{(I)}=\frac{p_{e}}{4}\left\lbrace \frac{\gamma }{c^2}(\omega ^2+\omega ^2_{{\rm BV}})|\xi _{z,e}|^2+\frac{1}{\gamma }\left(1+\frac{c^2k_{x}^2}{\omega ^2}\right)\left|\frac{p^{\prime }_{e}}{p_{e}}\right|^2\right\rbrace. \end{equation} \tag{ A15 }$$

Integrating Equation (A15) over the domain (which is from the interface z = −z₀ to z = ∞ or w = λ to w = −∞) and multiplying by our model surface area, we obtain

$$\begin{equation} E^{(I)}_{n}=4\pi R^2_{\odot }\frac{g\rho _{0}z^2_{0}}{4a}\lambda ^2\frac{|\mathcal {A}|^2}{z^2_{0}}\mathcal {N}_{n}^{(2)}, \end{equation} \tag{ A16 }$$

where

$$\begin{eqnarray} \mathcal {N}_{n}^{(2)}&=&\left\lbrace \nu ^2+a+\left(1+\frac{\nu ^2}{4\kappa ^2}\right)\left[\frac{2\kappa ^2(i\lambda _{z,e}-a+1)}{\nu ^2-4\kappa ^2}\right]^2\right\rbrace\nonumber\\ &&\times \left(\frac{1}{-i\lambda _{z,e}}\right)\left(\left.\frac{dQ}{dw}\right|_{w=\lambda }\right)^2. \end{eqnarray} \tag{ A17 }$$

The total energy contained in the f- and p-mode oscillations for the two-region model is $E_{n}=E^{(P)}_{n}+E^{(I)}_{n}$.

A.5. Lateral Energy Flux of f- and p-mode Oscillations

The external (to the flux tube) wavefield in both the polytropic and isothermal regions can be converted from Cartesian (x, y, z) to cylindrical polar coordinates (r,ϕ,z) with its origin centered on the flux tube (i.e., r = 0) using the Jacobi–Anger expansion. Thus, in both regions the wavefields functional dependence on r, ϕ, and t is

$$\begin{equation} f(r,\phi,t)=\frac{1}{2}e^{i(m\phi -\omega t)}H^{(2)}_{m}(kr). \end{equation} \tag{ A18 }$$

The energy flux of f and p modes for the prescribed model is given by (see Bray & Loughhead 1974)

$$\begin{equation} \boldsymbol {F}=p^{\prime }_{e}\boldsymbol {v}_{e}. \end{equation} \tag{ A19 }$$

Taking the time average and integrating over a cylinder of radius r gives the lateral energy flux of f- and p-mode oscillations which we shall denote F_n. The polytropic component of the lateral energy flux, i.e., from z = −∞ to z = −z₀, can be written as (see Jain et al. 2009)

$$\begin{equation} F^{(P)}_{n}=-\frac{|\mathcal {A}|^2}{z_{0}^2}\rho _{0}z_{0}^2\left(\frac{gz_{0}}{a}\right)^{3/2}\mathcal {H}^{(1)}_{n}, \end{equation} \tag{ A20 }$$

where

$$\begin{equation} \mathcal {H}^{(1)}_{n}=\frac{\nu ^3}{2}\left(\frac{\kappa }{\nu ^2}\right)^{a+1}\int ^\infty _{\frac{\nu ^2}{\kappa }}dw w^aQ^2(w). \end{equation} \tag{ A21 }$$

We now consider the calculation of the isothermal component of the lateral energy flux, i.e., from z = −z₀ to z = ∞,

$$\begin{equation} \langle F^{(I)}_{r}\rangle =\int ^{2\pi }_{0}rd\phi \int ^{\infty }_{-z_{0}}\frac{1}{4}(p^{\prime }_{e}v^{\prime *}_{r,e}+p^{\prime *}_{e}v^{\prime }_{r,e})dz. \end{equation} \tag{ A22 }$$

First, we calculate $v^{\prime }_{r,e}=-i\omega \xi _{r,e}$ as

$$\begin{equation} v^{\prime }_{r,e}=-\frac{2i}{\omega }\frac{\partial }{\partial r}\left(\frac{p^{\prime }_{e}}{\rho _{e}}\right). \end{equation} \tag{ A23 }$$

In cylindrical polar coordinates $p^{\prime }_{e}$ can be written using (A18) in the form

$$\begin{equation*} p_{e}^{\prime }=-\gamma p_{e}\lambda \left.\frac{dQ}{dw}\right|_{w=\lambda }\frac{\mathcal {A}}{2z_{0}}\left[\frac{2\kappa ^2(i\lambda _{z,e}-a+1)}{\nu ^2-4\kappa ^2}\right] \end{equation*}$$

$$\begin{equation} \times e^{\frac{1}{2}\left(i\lambda _{z,e}+a+1\right)\left(1-\frac{w}{\lambda }\right)}H^{(2)}_{m}(kr)e^{i(m\phi -\omega t)}. \end{equation} \tag{ A24 }$$

Substituting Equations (A23) and (A24) into (A22) and using the Wronskian for the Hankel functions, we obtain

$$\begin{equation} F^{(I)}_{n}=-\frac{|\mathcal {A}|^2}{z_{0}^2}\rho _{0}z_{0}^2\left(\frac{gz_{0}}{a}\right)^{3/2}\mathcal {H}^{(2)}_{n}, \end{equation} \tag{ A25 }$$

where

$$\begin{equation} \mathcal {H}^{(2)}_{n}=\nu ^3\left(\frac{1}{-i\lambda _{z,e}}\right)\left(\left.\frac{dQ}{dw}\right|_{w=\lambda }\right)^2\left[\frac{2\kappa (i\lambda _{z,e}-a+1)}{\nu ^2-4\kappa ^2}\right]^2. \end{equation} \tag{ A26 }$$

The time-averaged lateral energy flux of solar f- and p-mode oscillations for the two-region model is $F_{n}=F^{(P)}_{n}+F^{(I)}_{n}$.

Under the thin flux tube approximation, we derive a set of MHD equations for an inviscid and infinitely conducting medium. On linearizing these equations and assuming a time dependence consistent with the external perturbations, we can describe the interaction of the field-free medium with the thin flux tube for the sausage mode in both the polytropic and isothermal atmosphere. Thus, using the formulation of Hindman & Jain (2008), Jain et al. (2009), and Gascoyne et al. (2011) for the polytropic atmosphere and Hasan (1997) for the isothermal atmosphere, we obtain the governing equations for the propagation of sausage waves for our prescribed two-region atmosphere

$$\begin{equation} \left[c^2_{{\rm T}}\frac{\partial ^2}{\partial z^2}-\frac{\gamma g}{2}\frac{c^2_{{\rm T}}}{c^2}\frac{\partial }{\partial z}+\omega ^2+N_{g}\right]\xi _{z}=\frac{c^2_{{\rm T}}}{\rho V^2_{{\rm A}}}\left(\frac{\partial }{\partial z}+\frac{g}{c^2}\right)p_{e}^{\prime }, \end{equation} \tag{ B1 }$$

where $V^2_{{\rm A}}=B^2/(4\pi \rho)$ and $c^2_{{\rm T}}=c^2/(1+c^2/V^2_{{\rm A}})$ are the squares of Alfvén speed and tube speed, respectively. The governing equations for the propagation of sausage waves in the two regions are essentially the same apart from the term N_g ($=-\omega _{{\rm BV}}^2\left(c^2+\frac{\gamma }{2}V^2_{{\rm A}}\right)/(c^2+V^2_{{\rm A}})$), which permits the propagation of gravity waves in the isothermal layer; this quantity is zero for the polytrope which is in adiabatic equilibrium and therefore gravity waves cannot propagate. For the two atmospheres we obtain the following solutions to Equation (9):

$$\begin{equation} \xi ^{(I)}_{\Vert }=-\mathcal {A}[B_{1}e^{b_{1}(1-s)}+B_{2}e^{b_{2}(1-s)}+\mathcal {K}e^{b_{3}(1-s)}] \end{equation} \tag{ B2 }$$

$$\begin{equation} \xi ^{(P)}_{\Vert }=-\frac{i\pi }{2}\mathcal {A}\lbrace \psi _{\Vert }(s)[\Omega +\mathcal {J}^*(s)]+\psi _{\Vert }^*(s)[\mathcal {I}-\mathcal {J}(s)]\rbrace \end{equation} \tag{ B3 }$$

where * denotes the complex conjugate and

$$\begin{eqnarray} &&b_{1}=\frac{1}{4}(2i\lambda _{\Vert }+a+1),\quad b_{2}=-\frac{1}{4}(2i\lambda _{\Vert }-a-1),\nonumber\\ &&b_{3}=\frac{1}{2}(i\lambda _{z,e}+a+1),\nonumber\\ &&\lambda _{\Vert }=2k_{\Vert }z_{0}=\sqrt{2(\nu ^2-a)(2+\gamma \beta)-\frac{(a-3)^2}{4}} \end{eqnarray} \tag{ B4 }$$

$$\begin{equation} \mathcal {K}=\lambda \left.\frac{dQ}{dw}\right|_{w=\lambda }\frac{2\gamma (1+\beta)(\nu ^2-a)}{\left\lbrace \left(i\lambda _{z,e}+\frac{a+1}{2}\right)^2+\lambda ^2_{z}\right\rbrace }, \end{equation} \tag{ B5 }$$

$$\begin{eqnarray} \psi _{\Vert }(s)&=&s^{-\mu /2}H_{\mu }^{(1)}(2\nu \sqrt{\epsilon s}), \end{eqnarray} \tag{ B6 }$$

$$\begin{eqnarray} \mathcal {J}(s)&=&\int ^{s}_{1}r^{\mu +1}\psi _{\Vert }(r)f(r)dr, \end{eqnarray} \tag{ B7 }$$

$$\begin{eqnarray} \mathcal {I}&=&\int ^{\infty }_{1}r^{\mu +1}\psi _{\Vert }(r)f(r)dr, \end{eqnarray} \tag{ B8 }$$

$$\begin{eqnarray} f(s)&=&-\frac{(1+a)(\beta +1)}{2a}\frac{\nu ^2}{s}\frac{d Q(s)}{d s}. \end{eqnarray} \tag{ B9 }$$

We have also introduced the dimensionless parameters s = −z/z₀ and = (2a + β(1 + a))/(2a). Assuming a finite energy density at the upper boundary (s → ∞), we set B₂ = 0.

B.1. Matching Conditions for the Sausage Tube Waves at the Interface

At the interface s = 1 we require the two solutions $\xi ^{(I)}_{\Vert }$ and $\xi ^{(P)}_{\Vert }$ (Equations (B2) and (B3) respectively) to be equal, thus we obtain a relation for B₁,

$$\begin{equation} B_{1}=\frac{i\pi }{2}\left(H_{\mu }^{(1)}(2\nu \sqrt{\epsilon })\Omega +H_{\mu }^{(2)}(2\nu \sqrt{\epsilon })\mathcal {I}\right)-\mathcal {K}. \end{equation} \tag{ B10 }$$

The only other parameter needing to be specified is the constant Ω, which is calculated by ensuring that the Lagrangian pressure perturbation inside the tube is continuous across the interface. We obtain the following result:

$$\begin{equation} \Omega =\frac{i}{\pi \mathcal {R}_{1}}[(\gamma \nu ^2Q_{0}-\Lambda)(\beta +1)+2(b_{1}-b_{3})\mathcal {K}]-\frac{\mathcal {R}_{2}}{\mathcal {R}_{1}}\mathcal {I} \end{equation} \tag{ B11 }$$

where

$$\begin{equation} \mathcal {R}_{1}=\nu \sqrt{\epsilon }H_{\mu +1}^{(1)}(2\nu \sqrt{\epsilon })-b_{1}H_{\mu }^{(1)}(2\nu \sqrt{\epsilon }), \end{equation} \tag{ B12 }$$

$$\begin{equation} \mathcal {R}_{2}=\nu \sqrt{\epsilon }H_{\mu +1}^{(2)}(2\nu \sqrt{\epsilon })-b_{1}H_{\mu }^{(2)}(2\nu \sqrt{\epsilon }). \end{equation} \tag{ B13 }$$

In Figure 3, we plot the displacement for the incident p₂-mode with its corresponding excited sausage wave. The inset plot for s > 1 shows the familiar polytropic solutions for the p₂-mode and sausage wave. The s < 1 isothermal region shows the evanescent tail which grows with height due to the decrease in density. Figure 3 clearly shows the smooth matching of the wave solutions across the interface s = 1.

The wave energy flux at any point along the tube is given by the following expression (Bray & Loughhead 1974):

$$\begin{equation} \boldsymbol {F}=p^{\prime }\boldsymbol {v}+\frac{1}{4\pi }[(\boldsymbol {B}\cdot \boldsymbol {B}^{\prime })\boldsymbol {v}-(\boldsymbol {B}^{\prime }\cdot \boldsymbol {v})\boldsymbol {B}], \end{equation} \tag{ C1 }$$

where p', $\boldsymbol {v}$, and $\boldsymbol {B}^{\prime }$ are the pressure perturbation, velocity vector, and magnetic field perturbation vector, respectively. Reducing Equation (C1) to obtain the parallel component (with respect to the tube axis),

$$\begin{equation} F_{\Vert }=p_{e}^{\prime }\frac{\partial \xi _{\Vert }}{\partial t}-\frac{BB_{\Vert }^{\prime }}{4\pi }\frac{\partial \xi _{\Vert }}{\partial t}. \end{equation} \tag{ C2 }$$

The perturbed magnetic field $B_{\Vert }^{\prime }$ can be obtained from the derivation of the governing Equation (9) and takes the same form in both layers,

$$\begin{eqnarray} \frac{BB_{\Vert }^{\prime }}{4\pi }&=&\frac{2}{(2+\gamma \beta)}p^{\prime }_{e}+\frac{\gamma \beta }{(2+\gamma \beta)}\frac{B^2}{4\pi }\frac{\partial \xi _{\Vert }}{\partial z}\nonumber\\ &&-g\left(\rho -\frac{\gamma \beta }{(2+\gamma \beta)}\rho _{e}\right)\xi _{\Vert }. \end{eqnarray} \tag{ C3 }$$

Substituting Equations (A10), (A11), and (C3) into Equation (C2) and taking the time average, after substantial manipulation, we obtain the energy flux of sausage tube waves escaping from the magnetic flux tube as

$$\begin{equation} E_{\Vert }^{(u)}=\frac{\gamma \beta }{(2+\gamma \beta)}\frac{2g\rho _{0}z^2_{0}\omega A_{0}}{4(a+1)(\beta +1)}\frac{|\mathcal {A}|^2}{z_{0}^2}[{\rm Re}(\lambda _{\Vert })|B_{1}|^2+\mathcal {S}], \end{equation} \tag{ C4 }$$

$$\begin{equation} E_{\Vert }^{(d)}=-\frac{\gamma \beta }{(2+\gamma \beta)}\frac{\pi g\rho _{0}z_{0}^2\omega A_{0}}{4(a+1)(\beta +1)}\frac{|\mathcal {A}|^2}{z_{0}^2}|\Omega +\mathcal {I}^*|^2, \end{equation} \tag{ C5 }$$

where

$$\begin{eqnarray} \mathcal {S}&=&e^{\frac{1}{2}i\lambda _{z,e}(1-s)}[2\mathcal {K}{\rm Im}(b_{1}B_{1}e^{b_{1}(1-s)})\nonumber\\ &&-\lbrace \Lambda (\beta +1)+2b_{3}\mathcal {K}\rbrace {\rm Im}(B_{1}e^{b_{1}(1-s)})]|_{s\rightarrow -\infty }, \end{eqnarray} \tag{ C6 }$$

$$\begin{equation} \Lambda =-\gamma \lambda Q^{\prime }_{0}\left[\frac{2\kappa ^2(i\lambda _{z,e}-a+1)}{\nu ^2-4\kappa ^2}\right]. \end{equation} \tag{ C7 }$$

Note above that "Re" and "Im" denote the real and imaginary parts, respectively, of the corresponding complex quantities inside the parenthesis. $E_{\Vert }^{(u)}$ and $E_{\Vert }^{(d)}$ denote the wave energy flux calculated at z → ∞ (s → −∞) and z → −∞ (s → ∞), respectively. We have also multiplied the energy fluxes by the cross-sectional area of the tube

$$\begin{equation} A(s)=\left\lbrace \begin{array}{@{}ll@{}} A_{0}e^{(\mu +1)(1-s)} & \quad s\le 1\, {\rm isothermal}\\ A_{0}s^{-(\mu +1)} & \quad s\ge 1\ {\rm polytrope}.\\ \end{array} \right. \end{equation} \tag{ C8 }$$

We find that $E_{\Vert }^{(u)}$ is zero for the frequency range under consideration.

C.1. p-mode Damping Rates and Absorption Coefficients

We define the damping rate Γ for a single tube as (see also Hindman & Jain 2008)

$$\begin{equation} \Gamma =-\frac{1}{2\pi }\frac{E_{\Vert }}{E_{n}}=-\frac{1}{2\pi }\frac{E^{(u)}_{\Vert }+E^{(d)}_{\Vert }}{E^{(P)}_{n}+E^{(I)}_{n}}. \end{equation} \tag{ C9 }$$

Substituting Equations (A13), (A16), (C4), and (C5) into (C9), we obtain

$$\begin{equation} \Gamma =\frac{\beta }{4(\beta +1)\epsilon }\frac{\omega A_{0}}{4\pi R^2_{\odot }\lambda ^2}\frac{\left|\Omega +\mathcal {I}^*\right|^2}{\mathcal {N}_{n}}, \end{equation} \tag{ C10 }$$

where $\mathcal {N}_{n}=\mathcal {N}^{(1)}_{n}+\mathcal {N}^{(2)}_{n}$ (see Equations (A14) and (A17)). Multiplying by the number of tubes on the solar surface $N=4\pi R^2_{\odot }f/A_{0}$ where f is the magnetic filling factor (fraction of the solar surface the tubes occupy), we compute the total damping rate for f and p modes.

We follow the same methodology as in Jain et al. (2009) and Gascoyne et al. (2011) and define the absorption coefficient α for a single tube as

$$\begin{equation} \alpha =\frac{E_{\Vert }}{F_{n}}=\frac{E^{(u)}_{\Vert }+E^{(d)}_{\Vert }}{F^{(P)}_{n}+F^{(I)}_{n}}. \end{equation} \tag{ C11 }$$

Substituting Equations (A20), (A25), (C4), and (C5) into (C11) after simplification, we obtain

$$\begin{equation} \alpha =\frac{\pi \beta }{4(\beta +1)\epsilon }\frac{A_{0}}{z_{0}^2}\frac{\left|\Omega +\mathcal {I}^*\right|^2}{\mathcal {H}_{n}}, \end{equation} \tag{ C12 }$$

where $\mathcal {H}_{n}=\mathcal {H}^{(1)}_{n}+\mathcal {H}^{(2)}_{n}$ (see Equations (A21) and (A26)).