DETERMINATION OF TRANSVERSE DENSITY STRUCTURING FROM PROPAGATING MAGNETOHYDRODYNAMIC WAVES IN THE SOLAR ATMOSPHERE

Observations have shown that propagating magnetohydrodynamic (MHD) waves are present in the solar atmosphere. Small amplitude propagating transverse MHD waves have received particular attention because of their potential for energy storage, transfer, and deposition in the context of wave-based plasma heating processes. A key problem yet to be solved is the quantification of the role of waves in coronal heating. The solution requires reliable knowledge of the physical properties of magnetic and plasma structures. Because physical properties, such as the magnetic field strength, the plasma density, and their field-aligned and cross-field structuring, cannot be directly measured, seismology of MHD waves offers an alternative method to probe coronal plasmas.

The presence of ubiquitous coronal transverse waves has been reported by Tomczyk et al. (2007) and Tomczyk & McIntosh (2009) in observations taken with the Coronal Multi-Channel Polarimeter. The disturbances have amplitudes of the order of 0.3 km s⁻¹ and propagate at speeds of about 0.6 Mm s⁻¹. They are interpreted as transverse MHD waves of Alfvénic character. The measured discrepancy between inward and outward power associated with the disturbances (Tomczyk & McIntosh 2009) is thought to be an indication of in situ wave damping. Pascoe et al. (2010) has proposed resonant absorption of kink waves to explain this damping. This process, studied in the context of standing kink waves in coronal loops (see, e.g., Goossens et al. 1992, 2002, 2006; Ruderman & Roberts 2002) and in prominence fine structures (Arregui et al. 2008; Soler et al. 2009), occurs because of transverse inhomogeneity of the plasma across the waveguides. The resonant absorption process predicts selective damping as a function of frequency (Terradas et al. 2010). The good agreement between the observationally measured and the theoretically predicted outward/inward power ratio as a function of wave frequency strongly supports the idea that resonant absorption operates on these waves (Verth et al. 2010).

The use of wave damping properties to perform seismology inversions with transverse waves was proposed by Arregui et al. (2007) and Goossens et al. (2008) in the context of standing kink waves in coronal loops, and by Arregui & Ballester (2011) in the context of prominence fine structures (see reviews by Goossens 2008, Arregui et al. 2012, and Arregui 2012). Because the number of unknowns is usually larger than that of observables, density measurements need to be used as additional observational information to fully constrain the parameters of interest (Arregui & Asensio Ramos 2011).

The mechanism of resonant absorption does not make a distinction between standing and propagating MHD kink waves. The dynamics corresponds to a surface Alfvén wave damped by resonant absorption (Goossens et al. 2009, 2012a). For standing kink waves, the observational consequence is the attenuation of the amplitude in time. For propagating kink waves, the observational consequence is the attenuation of wave amplitude in space.

Goossens et al. (2012b) have presented a seismology inversion scheme for propagating MHD waves damped by resonant absorption. In their analysis, a solution curve that relates the Alfvén speed, the density contrast, and the transverse inhomogeneity length scale is obtained. The Alfvén speed is constrained to a narrow range, but the two parameters that define the cross-field density structuring cannot be inferred.

Pascoe et al. (2012) have pointed out that a more general damping profile consisting of a Gaussian and an exponential profile better reproduces the spatial decay of the amplitude for propagating kink waves. A Gaussian seems to better reproduce the amplitude behavior initially, while the exponential profile properly describes the dynamics after a few wavelengths. Their numerical solutions point to a transition between the two damping regimes. Currently, only indirect observational evidence about in situ damping of propagating transverse waves is available, with no measurement of the damping length scales. The confirmation on the existence of two damping regimes would make information on two damping length scales and the height at which the transition between the two regimes occurs accessible.

In this Letter, we demonstrate how this information can be used to fully constrain the cross-field density structuring in coronal waveguides.

The theory for the spatial damping of propagating kink waves due to resonant absorption has been developed by a number of studies (Terradas et al. 2010; Verth et al. 2010; Soler et al. 2011a, 2011b, 2011c; Hood et al. 2013). Numerical simulations have confirmed the obtained damping properties by analyzing spatial and temporal properties of the mode coupling process in coronal loops and arbitrary inhomogeneous coronal structures (Pascoe et al. 2010, 2011, 2012, 2013). The classic theoretical model assumes that transverse kink waves are guided by plasma structures. The waveguides are considered to be density tubes in the zero plasma-β approximation with a uniform magnetic field directed along the axis of a straight cylindrically symmetric structure. They have a uniform internal density, ρ_i, and a uniform external density, ρ_e, with ρ_i > ρ_e, connected in the transverse direction by a non-uniform region of thickness l/R, with R the tube radius. Because of the non-uniformity, resonant absorption produces the decay of the wave amplitude in space as the wave propagates.

For propagating waves with a real frequency ω, an analytical expression for the wavelength in the thin tube approximation is (see Terradas et al. 2010)

$$\begin{equation} \lambda =\sqrt{2} v_{\rm A,i} T \left(\frac{\zeta +1}{\zeta }\right)^{-1/2}. \end{equation} \tag{ 1 }$$

Here v_{A, i} is the internal Alfvén velocity, T is the period, and ζ = ρ_i/ρ_e is the density contrast. This solution is valid provided ωR/v_{A, i} ≪ 1. It expresses a relation between two observable quantities (wavelength and period) and two quantities to be determined (Alfvén velocity and density contrast).

Because of resonant absorption, spatial damping occurs and the transverse velocity amplitude decays with an exponential profile of the form exp (− z/L_d). Under the thin tube and thin boundary (l/R ≪ 1) approximation, an expression for the damping length, L_d, as a function of the relevant physical parameters can be obtained. In units of the wavelength, this expression is (see Terradas et al. 2010)

$$\begin{equation} \frac{L_{\rm d}}{\lambda }= \left(\frac{2}{\pi }\right)^2\left(\frac{R}{l}\right)\left(\frac{\zeta +1}{\zeta -1}\right). \end{equation} \tag{ 2 }$$

The first factor is due to the assumed linear density profile at the non-uniform layer. Note that the right-hand side of this expression is identical to the one for the damping time over the period for standing kink waves. The reason is that resonant absorption does not make any distinction with respect to the standing or propagating character of the wave.

The exponential profile obtained for standing (e.g., Ruderman & Roberts 2002; Goossens et al. 2002) and propagating (e.g., Terradas et al. 2010) kink waves describes the asymptotic state of the damping behavior, i.e., at large times or distances. Pascoe et al. (2012) demonstrated with numerical simulations that the initial damping stage can be described by a Gaussian profile of the form $\exp (-z^2/L^2_{\rm g})$, with L_g the Gaussian damping length scale. Hood et al. (2013) considered the problem analytically and produced an expression for the full nonlinear spatial damping profile, which can be approximated as Gaussian for low heights and exponential at large heights. Instead, Pascoe et al. (2013) proposed a general spatial damping profile composed of a Gaussian damping profile at low heights and the usual exponential profile at large heights. An example of the spatial dependence of the velocity amplitude from numerical simulations and the double profile fitting for such a general damping profile is displayed in Figure 1. The accuracy of this approximate damping profile was demonstrated by the parametric study performed by Pascoe et al. (2013). This study shows that the Gaussian damping length scale can be well described by the expression

$$\begin{equation} \frac{L_{\rm g}}{\lambda }= \left(\frac{2}{\pi }\right) \left(\frac{R}{l}\right)^{1/2}\left(\frac{\zeta +1}{\zeta -1}\right). \end{equation} \tag{ 3 }$$

This equation expresses the Gaussian damping length as a function of the same two parameters that determine the exponential damping length. This means that the observational identification of two damping regimes and the measurement of their associated length scales would provide us with additional information without the inclusion of new model parameters. The height, h, at which the damping regime changes from Gaussian to exponential is given by (see Pascoe et al. 2013)

$$\begin{equation} h=\frac{L^2_{\rm g}}{L_{\rm d}}=\lambda \left(\frac{\zeta +1}{\zeta -1}\right). \end{equation} \tag{ 4 }$$

Zoom In Zoom Out Reset image size

The nonlinear (non-exponential) damping rate is an inherent property of resonantly damped kink oscillations, as derived analytically by Hood et al. (2013) for spatial damping. The previous analyses by, e.g., Ruderman & Roberts (2002) and Goossens et al. (2002) only attempt to describe the damping in the asymptotic state and so they only obtain the later exponential damping stage.

The observational identification of this feature would be strong support for the resonant damping mechanism. However, no measurement of the axial damping spatial scales is available yet. Our analysis shows that, if present, Equations (2)–(4) can be used to fully determine the transverse density structuring in oscillating coronal waveguides.

Out of the three equations, only two are independent. In our inversion scheme, Equations (3) and (4) will be used to infer ζ and l/R from data given by L_g and h. The two algebraic equations can be directly solved for the two unknowns. However, the problem we intend to solve is not exact since real data are noisy and measurements have an associated uncertainty. To accommodate a proper propagation of uncertainty, we make use of Bayesian analysis, a framework that has now started developing in coronal seismology (Arregui & Asensio Ramos 2011; Arregui et al. 2013).

To perform our inversion, we use Bayes' theorem (Bayes & Price 1763)

$$\begin{equation} p(\mbox{{\boldsymbol {\theta }}} | d)=\frac{p(d | \mbox{{\boldsymbol {\theta }}})p(\mbox{{\boldsymbol {\theta }}})}{p(d)}, \end{equation} \tag{ 5 }$$

which gives the solution to the inverse problem in terms of the posterior probability distribution, p(θ|d), that describes how probability is distributed among the possible values of the unknown parameter, θ, given the data d. The posterior is a combination of what is known about the parameters before taking into account the data, the prior p(θ), and the likelihood function, p(d|θ), which tells us how well the model predicts the observed data. The denominator, p(d), is the so-called evidence that normalizes the likelihood and plays no role in parameter inference.

Once the posterior is known, we calculate how a particular parameter is constrained by data computing its marginal posterior. This is done by performing the following integral of the full posterior with respect to the rest of parameters:

$$\begin{equation} p(\theta _i|d) = \int p(\mbox{{\boldsymbol {\theta }}} | d) d\theta _1 \ldots d\theta _{i-1} d\theta _{i+1} \ldots d\theta _{N}. \end{equation} \tag{ 6 }$$

The result provides us with all the information for model parameter θ_i available in the priors and the data. This method also enables us to correctly propagate uncertainties from data to inferred parameters.

We next specify the likelihood function and the priors. In what follows, we assume the observed data are given by d = (L_g, h), where both observed length scales are normalized to the wavelength. The unknowns are gathered in the vector θ = (ζ, l/R). Under the assumption that observations are corrupted with Gaussian noise and they are statistically independent, the likelihood can be expressed as

$$\begin{eqnarray} p(d|{\boldsymbol {\theta }}) &=&(2\pi \sigma _{L_{\rm g}} \sigma _{h})^{-1} \nonumber\\ &&\times\exp \left\lbrace \frac{\left[L_{\rm g}-L^\mathrm{syn}_{\rm g}({\boldsymbol {\theta }})\right]^2}{2 \sigma _{L_{\rm g}}^2} + \frac{\left[h-h^\mathrm{syn}({\boldsymbol {\theta }})\right]^2}{2 \sigma _{h}^2} \right\rbrace,\nonumber\\ \end{eqnarray} \tag{ 7 }$$

with $L_{\rm g}^\mathrm{syn}({\boldsymbol {\theta }})$ and h^syn(θ) given by Equations (3) and (4), respectively. Likewise, $\sigma _{L_{\rm g}}^2$ and $\sigma _{h}^2$ are the variances associated with the Gaussian damping length and the height h, respectively.

The priors indicate our level of knowledge (ignorance) before considering the observed data. We have adopted uniform prior distributions for both unknowns over the given ranges, so that we can write

$$\begin{equation} p(\theta _i)=\frac{1}{\theta ^{\mathrm{max}}_i-\theta ^{\mathrm{min}}_i} \mbox{\hspace{5.69046pt}}\mbox{for} \mbox{\hspace{5.69046pt}}\theta ^{\mathrm{min}}_i\le \theta \le \theta ^{\mathrm{max}}_i, \end{equation} \tag{ 8 }$$

and zero otherwise. For the minimum and maximum values the intervals ζ ∈ (1, 2] and l/R ∈ (0, 2] have been taken, respectively. This choice of priors expresses our belief that the unknown parameters are constrained to those ranges, with all combinations being equally probable.

Computations were repeated for different uniform prior ranges for ζ, with ζ_min = 1.1 and ζ_max = 5, 10, 50, and for l/R, with (l/R)_min = 0.01 and (l/R)_max = 1, 1.5, 2. We also used strongly informative priors (Gaussian distributions in ζ and l/R centered far away the synthetic parameter values). For strong prior information on l/R, convergence to the correct posterior is achieved through repeated application of Bayes' theorem using the posterior of one iteration as a prior for the next one and multiplying again by the likelihood. This means that a few consistent observations of the data lead to a converged posterior independently of the prior. Our sensitivity analyses lead us to conclude that posteriors are dominated by the information contained in the data that overwhelms the prior information.

Posteriors are evaluated for different combinations of parameters using Bayes' theorem (Equation (5)) and inferred distributions are obtained from the computation of the marginal posteriors using Equation (6). To solve the one-dimensional integrals, numerical quadratures were employed using an adaptive Gauss–Kronrod quadrature. All computations were checked by comparison with the results obtained with the MCMC code employed by Arregui & Asensio Ramos (2011).

We first evaluated the performance of our inversion scheme by making the inference under controlled conditions. We generated predictions for the length scales L_g and h for different combinations of the equilibrium parameters, ζ = 1.5, 2, 3, 4 and l/R = 0.05, 0.15, 0.2, 0.4, using Equations (3) and (4). These synthetic data were treated as observed data in the Bayesian inversion. A 10% uncertainty on the data was considered and the posterior distributions for ζ and l/R were computed using the likelihood function (Equation (7)) and the uniform priors. Once the posteriors were known, the median and the variances associated with the 68% confidence level were calculated. Table 1 displays the inversion results for some parameter combinations. In all the cases, the inversion method correctly recovers the values for the physical parameters. The larger the density contrast, the shorter the two length scales L_g and h. This increases the errors in the inferred density contrast, while errors in l/R are not affected that much. For the combinations with the largest ζ = 10 and l/R = 0.5, 1, 1.5, L_g and h are comparable to the wavelength. This will make the observational identification of the two damping regimes very problematic.

Table 1. Inversion of Synthetic Data Using the Analytical Forward Model

Synthetic Parameters		Synthetic Data		Inversion Results
ζ	l/R	Lg/λ	h/λ	ζ	l/R
1.5	0.05	14.2	5.0	1.51$^{+0.08}_{-0.06}$	0.05$^{+0.02}_{-0.01}$
1.5	0.15	8.2	5.0	1.50$^{+0.07}_{-0.06}$	0.16$^{+0.05}_{-0.04}$
1.5	0.2	7.1	5.0	1.51$^{+0.07}_{-0.06}$	0.21$^{+0.06}_{-0.05}$
1.5	0.4	5.0	5.0	1.50$^{+0.07}_{-0.05}$	0.44$^{+0.13}_{-0.11}$
3	0.05	5.7	2.0	3.11$^{+0.59}_{-0.38}$	0.05$^{+0.02}_{-0.01}$
3	0.15	3.3	2.0	3.09$^{+0.61}_{-0.40}$	0.15$^{+0.05}_{-0.04}$
3	0.2	2.9	2.0	3.13$^{+0.58}_{-0.41}$	0.19$^{+0.07}_{-0.05}$
3	0.4	2.0	2.0	3.10$^{+0.60}_{-0.41}$	0.42$^{+0.15}_{-0.12}$
4	0.05	4.8	1.7	4.31$^{+1.52}_{-0.79}$	0.05$^{+0.02}_{-0.01}$
4	0.15	2.7	1.7	4.39$^{+1.47}_{-0.85}$	0.15$^{+0.05}_{-0.04}$
4	0.2	2.4	1.7	4.38$^{+1.69}_{-0.85}$	0.19$^{+0.08}_{-0.06}$
4	0.4	1.7	1.7	4.38$^{+1.55}_{-0.86}$	0.38$^{+0.14}_{-0.11}$
10	0.5	1.1	1.2	11.54$^{+4.58}_{-3.88}$	0.51$^{+0.16}_{-0.11}$
10	1.0	0.8	1.2	11.55$^{+4.69}_{-3.81}$	1.02$^{+0.29}_{-0.22}$
10	1.5	0.6	1.2	12.29$^{+4.32}_{-3.89}$	1.45$^{+0.29}_{-0.28}$

Download table as:  ASCIITypeset image

Then the simulations of transverse kink wave propagation in a magnetic flux tube were performed using a numerical code (see Pascoe et al. 2013 for details). A Lax–Wendroff code is used to solve the linear MHD equations in cylindrical coordinates. The lower boundary is driven harmonically with velocity perturbations corresponding to the loop footprint moving back and forth about its equilibrium. The simulation ends after 10 periods of oscillation and the spatial damping profile is investigated by considering the radial velocity component, v_r, as a function of z at the center of the loop (Figure 1). From the behavior of the amplitude of the excited kink waves at different heights, the damping profile was fitted and values for L_g and h were obtained. Using those fitted values as data, we repeated the inversion procedure. For the sake of comparison, parameter spaces that overlap with those in Table 1 were considered.

Figure 2 displays an example of the marginal posterior distributions and the joint probability distribution for ζ and l/R. For both parameters, well-defined probability distributions are obtained. For each parameter, the median of the marginal posterior and errors given at the 68% credible interval are used to compute the estimates given in Table 2. The table shows the values for the physical parameters used in the simulations, the fitted length scales, and the inferred physical parameters. Numerical and analytical forward models give similar results. This issue is discussed in detail by Pascoe et al. (2013, see their Figures 8–10). Our Bayesian inference technique properly returns the physical parameters of interest. This means that the algebraic expressions for the forward problem are accurate enough since their use correctly returns the values for the physical parameters that were actually used in the numerical setup. As with synthetic data in Table 1, large density contrast values tend to produce larger errors in their determination by inversion. The main problem lies in obtaining the parameters L_g, L_d, and h from the data, and specifically in determining h accurately, which determines the accuracy of the density estimate.

Zoom In Zoom Out Reset image size

Table 2. Inversion of Numerical Data from Simulations

Simulation Parameters		Fitted Data		Inversion Results
ζ	l/R	Lg/λ	h/λ	ζ	l/R
1.5	0.05	11.5	3.8	1.73$^{+0.12}_{-0.09}$	0.05$^{+0.02}_{-0.01}$
1.5	0.15	7.9	4.6	1.56$^{+0.08}_{-0.07}$	0.15$^{+0.05}_{-0.04}$
1.5	0.2	7.0	4.8	1.53$^{+0.08}_{-0.06}$	0.21$^{+0.07}_{-0.05}$
1.5	0.4	5.0	4.9	1.52$^{+0.07}_{-0.06}$	0.39$^{+0.09}_{-0.08}$
3	0.05	5.5	2.1	2.88$^{+0.46}_{-0.33}$	0.06$^{+0.02}_{-0.02}$
3	0.15	3.5	2.2	2.74$^{+0.44}_{-0.32}$	0.16$^{+0.06}_{-0.04}$
3	0.2	3.1	2.2	2.74$^{+0.41}_{-0.30}$	0.21$^{+0.07}_{-0.05}$
3	0.4	2.1	2.0	3.09$^{+0.57}_{-0.40}$	0.38$^{+0.13}_{-0.11}$
4	0.05	4.9	1.7	4.17$^{+1.32}_{-0.74}$	0.05$^{+0.02}_{-0.01}$
4	0.15	3.1	1.9	3.19$^{+0.64}_{-0.42}$	0.16$^{+0.06}_{-0.05}$
4	0.2	2.7	1.9	3.33$^{+0.74}_{-0.43}$	0.21$^{+0.07}_{-0.06}$
4	0.4	2.3	2.2	2.73$^{+0.43}_{-0.29}$	0.38$^{+0.12}_{-0.10}$

Download table as:  ASCIITypeset image

The general spatial damping profile remains an accurate description of the damping behavior (Pascoe et al. 2013) for large density contrasts. The weak link is the errors in the least-squares fit of the damping profile to the data. The errors are larger for higher density contrasts due to the Gaussian stage being very short, and the results are more sensitive to the errors in h. The lesson is that large contrasts are a challenge from an observational point of view, if the theoretical predictions by Hood et al. (2013) and Pascoe et al. (2013) and the inversion technique presented here are to be employed.

The cross-field density structuring of magnetic waveguides cannot be currently estimated using seismology techniques. Yet such an understanding is crucial to assess the possible role of MHD waves in plasma heating processes. Recent theoretical results for the spatial damping of propagating MHD kink waves predict the existence of two different damping regimes for the spatial dependence of the velocity amplitude along the waveguides. In this Letter, we have explained how the observational identification of those regimes and the measurement of the associated damping length scales can be used to fully constrain the cross-field density structuring of magnetic waveguides.

Such measurements are currently unavailable and observational efforts toward the observational confirmation of the theoretical prediction and the measurement of the associated spatial scales should be devised. They would not only provide us with hitherto undetermined information on density, but would strongly endorse resonant absorption as damping/heating mechanism since the presence of two damping regimes is an inherent property of this mechanism.

With the use of Bayesian inference, such measurements would ensure that the inverse problem is consistently solved taking into account all the relevant information and with correct propagation of uncertainty. The range of parameters used in our inversion experiments covers the reasonable range of expected values used in previous studies (Goossens et al. 2002, 2012b; Ruderman & Roberts 2002; Arregui et al. 2007), so the measurement and inversion of observations seem feasible.

The application to real observations would require high-quality data so that a reliable fitting can be performed to obtain the damping length scales and to be able to clearly discriminate between Gaussian and exponential damping regimes. As pointed out by Pascoe et al. (2013), the Gaussian damping stage only applies for one wavelength in the limit of density contrast tending to infinity. Fitting to such a small signal can lead to significant error bars. The method of fitting both damping regimes is applicable for low-density contrast waveguides. Our application is limited to inference for a single oscillating magnetic tube. Observations by Tomczyk et al. (2007) and Tomczyk & McIntosh (2009) show that wave dynamics is ubiquitous over large coronal regions. Our method could be extended in order to apply the inversion technique to such extended regions, thus obtaining overall information about the cross-field density structuring of the corona where these waves propagate.

Recent studies have considered effects on the wave amplitude of propagating transverse kink waves from the consideration of a background field-aligned flow (Soler et al. 2011b) and of longitudinal density stratification (Soler et al. 2011c). These two effects affect the wave amplitude and might partially or totally compensate the damping by resonant absorption, thus producing significant differences in the inferred parameters obtained in this study. Future seismology schemes should include these ingredients in their inversion process.

I.A. and A.A.R. acknowledge support by Ramón y Cajal Fellowships by the Spanish Ministry of Economy and Competitiveness (MINECO). I.A. acknowledges the support by the Spanish MICINN/MINECO and FEDER funds through project AYA2011-22846. A.A.R. acknowledges support by the Spanish MINECO through project AYA2010-18029 (Solar Magnetism and Astrophysical Spectropolarimetry) and from the Consolider-Ingenio 2010 CSD2009-00038 project. The computational work for this Letter was carried out on the joint STFC and SFC (SRIF) funded cluster at the University of St Andrews (Scotland, UK).