BAYESIAN ANALYSIS OF MULTIPLE HARMONIC OSCILLATIONS IN THE SOLAR CORONA

Observations of magnetohydrodynamic (MHD) oscillations in solar coronal loops indicate the simultaneous presence of multiple mode harmonics in the same structure (Verwichte et al. 2004; De Moortel & Brady 2007; Van Doorsselaere et al. 2007, 2009; O'Shea et al. 2007). Their potential use for the diagnostic of the plasma conditions in the corona has attracted considerable interest (see Andries et al. 2009 for a review). The idea was first put forward by Andries et al. (2005) and Goossens et al. (2006), who found that the ratio of the fundamental mode period to twice that of its first overtone in the longitudinal direction depends on the density structuring along magnetic field lines. It is equal to unity in longitudinally uniform and unbounded tubes, but is smaller than one when density stratification is present. Using observational measurements of period ratios, Andries et al. (2005) and Van Doorsselaere et al. (2007) obtained estimates for the coronal density scale height.

An alternative hypothesis was formulated by Verth & Erdélyi (2008), whereby the expansion of the magnetic loop produces a deviation from unity in the period ratio, increasing its value. The effect is relevant for obtaining accurate estimates of the density scale height using multiple mode period ratios (Verth & Erdélyi 2008; Verth et al. 2008; Ruderman et al. 2008). Observational evidence for magnetic tube expansion has been reported by, e.g., Klimchuk (2000) in soft X-ray loops observed with Yohkoh and by Watko & Klimchuk (2000) in non-flare and postflare loops observed with Transition Region and Coronal Explorer. Observations of period ratios larger than one have been reported by O'Shea et al. (2007) and De Moortel & Brady (2007).

Since then, period ratio studies have analyzed different models for the density structuring in coronal loops (McEwan et al. 2006, 2008), the influence of the elliptic shape of loops (Morton & Erdélyi 2009), the twist of the magnetic field (Karami & Bahari 2012), the effect of the environment (Orza et al. 2012), and the temporal expansion of loops (Ballai & Orza 2012). Period ratios have also been analyzed for slow MHD modes (Macnamara & Roberts 2010) and for sausage modes (Macnamara & Roberts 2011) and their use has been suggested in the context of prominence seismology (Díaz et al. 2010; Arregui et al. 2012).

The application of Bayesian analysis techniques to coronal seismology is in its infancy. A first attempt of parameter inference using damped loop oscillations was presented by Arregui & Asensio Ramos (2011). In this Letter, we present the first application of Bayesian model comparison techniques to coronal seismology. We first perform Bayesian parameter inference for the coronal density scale height and the magnetic tube expansion, under the two hypotheses of density stratification and magnetic field divergence, using multiple mode oscillations. Then we assess which one of the two hypotheses better explains the observations for given values of the period ratio.

The deviation from unity for the ratio between the fundamental and the first overtone transverse kink oscillation periods has been attributed to two main physical effects. Each one constitutes a hypothesis to explain the data.

In the model by Andries et al. (2005), coronal density stratification produces a decrease of the period ratio. This model projects a vertically stratified isothermal atmosphere onto a semicircular loop. An analytical expression for the dependence of the period ratio on density scale height was obtained by Safari et al. (2007). This expression can be rewritten as

$$\begin{equation} r_1=\frac{P_1}{2P_2}=1-\frac{4}{5}\left(\frac{\eta }{\eta +3\pi ^2}\right), \end{equation} \tag{ 1 }$$

with P₁ and P₂ the periods of the fundamental and first overtone modes, η = L/πH the ratio of the loop height at the apex to the density scale height H, and L the loop length. We have checked that Equation (1) provides us with a good approximation to the numerical results by Andries et al. (2005).

In the model by Verth & Erdélyi (2008), magnetic tube expansion produces an increase of the period ratio given by

$$\begin{equation} r_2=\frac{P_1}{2P_2}=1+\frac{3(\Gamma ^2-1)}{2\pi ^2}, \end{equation} \tag{ 2 }$$

with the expansion defined as Γ = r_a/r_f, where r_a is the radius at the apex and r_f is the radius at the footpoint. Equation (2) was obtained under the assumption that |Γ − 1| ≪ 1. Moreover, it is the result for a particular expanding loop model. Another equilibrium state will produce a quantitatively different result, although Ruderman et al. (2008) anticipate a qualitatively similar result.

Both longitudinal stratification and magnetic tube expansion have forward models that relate one observed quantity, the period ratio P₁/2P₂, to one physical quantity to be inferred, η or Γ.

To perform the inference using existing estimates for the period ratio, we employ Bayes' theorem (Bayes & Price 1763)

$$\begin{equation} p(\mbox{{\boldsymbol $\theta $}} | D, M)=\frac{p(D | \mbox{{\boldsymbol $\theta $}}, M)p(\mbox{{\boldsymbol $\theta $}}|M)}{\int d\mbox{{\boldsymbol $\theta $}}p(D|\mbox{{\boldsymbol $\theta $}},M) p(\mbox{{\boldsymbol $\theta $}}|M)}, \end{equation} \tag{ 3 }$$

which gives the solution to the inverse problem in terms of the posterior probability distribution, p(θ|D, M), that describes how probability is distributed among the possible values of the unknown parameter, θ, given the data D and the assumed model M. The function p(D|θ, M) is the likelihood of obtaining a data realization actually observed as a function of the parameter vector and provides a measure of how well the data are predicted by the model. The prior probability p(θ|M) encodes any prior information we might have on the model parameters without taking into account the observed data. The denominator is the evidence, an integral of the likelihood over the prior distribution. This quantity plays no role in parameter inference, but will become central in the model comparison described in Section 4.

3.1. Longitudinally Stratified Loops

For longitudinally stratified loops, model M₁, the forward problem, is given in Equation (1). To evaluate the likelihood, we assume that the model is true. Then, the period ratio measurement (r) will differ from the prediction (r₁) because of measurement uncertainties (e), so that r = r₁ ± e. The probability of obtaining the measured value is equal to the probability of the error. Assuming Gaussian errors, the likelihood for model M₁ is then expressed in the following manner

$$\begin{equation} p(r | \eta ,M_1)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left[-\frac{(r-r_1)^2}{2\sigma ^2}\right], \end{equation} \tag{ 4 }$$

with σ² the variance associated to the observed period ratio. In the following we assign observed period ratio errors to the standard deviation σ.

The prior indicates our level of knowledge (ignorance) before considering the observed data. We have adopted a uniform prior distribution for the unknown, η, over a given range, so that we can write

$$\begin{equation} p(\eta |M_1)=\frac{1}{\eta ^{\mathrm{max}}-\eta ^{\mathrm{min}}} \mbox{\hspace{5.69046pt}}\mbox{for} \mbox{\hspace{5.69046pt}}\eta ^{\mathrm{min}}\le \eta \le \eta ^{\mathrm{max}}, \end{equation} \tag{ 5 }$$

and zero otherwise. We only consider loops with density decreasing with height. We know from observations that very large coronal loops with the apex height above several atmospheric scale heights exist, so we consider the range η ∈ [0, 8] in the prior above.

Parameter inference is next performed using Bayes' theorem (3). We use period ratio measurements reported in observations by Verwichte et al. (2004). We must note that unambiguous identification of measured periods with particular overtones is a difficult task. Early observations were limited to the measurement of different periodicities at a single point without spatial information (Verwichte et al. 2004). Even in the case in which information at different cuts along the loops is available (De Moortel & Brady 2007), it is unclear how one should assign those periods to theoretical overtones. As our Letter aims at providing a method to perform parameter inference and model comparison, we have considered period ratio values discussed in the review by Andries et al. (2009), noting that the reliability of the results is closely related to that of the adopted assumptions and theoretical interpretations.

Figure 1(a) shows posterior probability distributions for η computed using Equation (3) with likelihood and prior given by Equations (4) and (5) for two period ratio measurements by Verwichte et al. (2004). Well constrained distributions are obtained. For the measured period ratios r ∼ 0.79 and r ∼ 0.91, the inversion leads to density scale heights of H = 21 Mm and H = 56 Mm, respectively, for a loop with a height at the apex of L/π = 70 Mm.

Zoom In Zoom Out Reset image size

3.2. Expanding Magnetic Loops

For expanding loops, model M₂, the forward problem is given in Equation (2). To evaluate the likelihood, we follow the same procedure as above, which leads to

$$\begin{equation} p(r|\Gamma ,M_2)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left[-\frac{(r-r_2)^2}{2\sigma ^2}\right]. \end{equation} \tag{ 6 }$$

We also adopt a uniform prior distribution for the unknown Γ, so that we can write

$$\begin{equation} p(\Gamma |M_2)=\frac{1}{\Gamma ^{\mathrm{max}}-\Gamma ^{\mathrm{min}}} \mbox{\hspace{5.69046pt}}\mbox{for} \mbox{\hspace{5.69046pt}}\Gamma ^{\mathrm{min}}\le \Gamma \le \Gamma ^{\mathrm{max}}, \end{equation} \tag{ 7 }$$

and zero otherwise. We consider Γ ∈ [1, 2.5].

Figure 1(b) shows posterior probability distributions for Γ computed using Equation (3) with likelihood and prior given by Equations (6) and (7) for two period ratio measurements. Again, well constrained distributions are obtained. Andries et al. (2009) discuss period ratio measurements in Table 1 by De Moortel & Brady (2007). Assuming that either the most power is in the fundamental mode or in the first overtone, mean values for the period ratio of r ∼ 1.07 and r ∼ 1.38 are obtained. For those values, the inversion leads to tube expansion factors that are compatible with the estimates by Klimchuk (2000) and Watko & Klimchuk (2000). Note, however, that, according to Figure 1(b), a period ratio of r ∼ 1.38 requires an expansion of the tube by a factor of Γ ∼ 1.85, while observations by Watko & Klimchuk (2000) seem to indicate that in only very few cases does this parameter approach or exceed a value of 2.

For both theoretical models, the Bayesian framework makes use of all the available information in a consistent manner and enables us to consistently propagate errors from observations to inferred parameters.

Parameter inference was performed under the hypothesis that observed period ratios are the effect of either density stratification or magnetic tube expansion. The probability distributions in Figure 1 are conditional on the specific models used to explain the observations. The Bayesian framework enables us to present different models to the same data and assess in a quantitative manner which one is favored. We have compared three models: M₀ for a uniform density and magnetic field strength tube, M₁ for a longitudinally stratified loop, and M₂ for a expanding magnetic loop.

As model M₀ predicts a period ratio of one, regardless of η or Γ, we can write

$$\begin{equation} p(r|M_0)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left[-\frac{(r-1)^2}{2\sigma ^2}\right], \end{equation} \tag{ 8 }$$

for the likelihood in this case. Note that p(r|M₀) = p(r|η, M₀) = p(r|Γ, M₀).

To determine the plausibility of models M₁ and M₂ between them and with respect to model M₀, we evaluate the posterior probabilities to ascertain the relative merits of two models, M_i and M_j. This is done by applying Bayes' theorem (Equation (3)) to the two models and considering posterior ratios of the form (Jeffreys 1961)

$$\begin{equation} \frac{p(M_i|r)}{p(M_j|r)}=\frac{p(r|M_i)}{p(r|M_j)}\frac{p(M_i)}{p(M_j)}. \end{equation} \tag{ 9 }$$

The first ratio on the right-hand side is the Bayes factor. It expresses how well the observed data are predicted by model M_i, compared to model M_j. The second ratio, the prior odds ratio, measures how much our initial beliefs favored M_i over M_j, before considering the data. As we have no particular a priori preference for one model over the other, before considering the data, we take p(M_i) = p(M_j) = 1/2. Our assessment of the plausibility of models is then based on the computation of the Bayes factor of M_i against M_j given by

$$\begin{equation} BF_{ij}=\frac{p(r|M_i)}{p(r|M_j)}. \end{equation} \tag{ 10 }$$

In model selection, we are interested in the most probable model, independently of the parameters, i.e., we should marginalize out all parameters. This is achieved by performing an integral of the likelihood over the full parameter space. The marginal likelihood for a given model M_i is then

$$\begin{equation} \hspace{-2pt}p(r|M_i)=\int ^{\theta ^{\mathrm{max}}}_{\theta ^{\mathrm{min}}}\!p(r,\theta |M_i)d\theta =\int ^{\theta ^{\mathrm{max}}}_{\theta ^{\mathrm{min}}}\!p(r|\theta ,M_i)p(\theta |M_i)d\theta , \end{equation} \tag{ 11 }$$

where θ ∈ [θ^min, θ^max] represents the particular parameter of the model and we have used the product rule to expand the probability of r and θ, given model M_i.

Figure 2(a) displays the marginal likelihoods for the three considered models. For a given observed period ratio, the plausibility of one model over the other is given by the ratio of these two quantities at the measured period ratio r. The uniform model has a marginal likelihood that is maximum at one. The model with density stratification is clearly favored for observations of period ratios below unity. Figure 2(b) quantifies the relative performance of models M₁ and M₀ by computing the Bayes factor. The lower r, the more plausible M₁ is against M₀. Jeffreys' scale (Jeffreys 1961; Kass & Raftery 1995) assigns different levels of evidence to the values of the Bayes factor in natural logarithm units. They are shown using different darkening options in Figures 2(b)–(d). According to Figure 2(b), and given the assumed uncertainty of σ = 0.08, a period ratio measurement should be considered as positive evidence for model M₁ against model M₀ only if it is below 0.87 and as strong evidence only if it is below 0.78. A period ratio below 0.71 indicates very strong evidence for model M₁.

Zoom In Zoom Out Reset image size

A similar comparison was made for models M₂ and M₀ (Figure 2(c)). The model for magnetic tube expansion is clearly favored for observations of period ratios above unity. Model M₀ cannot be ruled out by the simple fact that r > 1. Positive evidence for model M₂ against model M₀ exists for period ratio measurements above 1.16. If r > 1.23, there is strong evidence for model M₂. A period ratio above 1.28 indicates very strong evidence for model M₂.

Finally, we have compared model M₁ against model M₂ in view of observed data. Figure 2(d) shows that the evidence is the same for both models if the period ratio is unity. Around this value, there is a region, 0.96 < r < 1.06, in which no firm conclusion can be established. Positive, strong, and very strong evidence for M₁ occur below 0.96, 0.87, and 0.80, respectively. Positive, strong, and very strong evidence for M₂ occur above 1.06, 1.15, and 1.21, respectively.

In our model comparison, σ = 0.08 has been selected so as to clearly show the different regimes for the evidence. An increase (decrease) of σ produces a decrease (increase) of any evidence. The inference in Figure 1(a) for r = 0.91 with uncertainty of σ = 0.04 falls into the region of positive evidence for hypothesis M₁. The inference in Figure 1(b) for r = 1.07 should have an uncertainty of σ = 0.03 (close to the reported error) to be considered done under positive evidence for hypothesis M₂.

We have presented a simple and straightforward technique to perform parameter inference and model comparison using period ratios of kink oscillations in the Bayesian framework. Parameter inference enables us to obtain estimates for the coronal density scale height and the magnetic loop expansion, using all the available information and with correctly propagated uncertainties. Density scale height estimates are compatible with previous studies. Magnetic tube expansion factors are compatible with the estimates by Klimchuk (2000) and Watko & Klimchuk (2000). Even if Equations (1) and (2) imply well-posed inversion problems from a mathematical point of view, the unknown parameters cannot be obtained uniquely due to the uncertainties associated to observed data. The Bayesian framework enables us to consistently deal with this problem.

Bayesian model comparison enables us to assess the plausibility of the inferences, which are conditional on theoretical models. Deviations of the period ratio below (above) unity do not necessarily imply density stratification (tube expansion) to be preferred in front of the null hypothesis. The degree of evidence for the two theories so far invoked can be quantitatively assessed with the use of Bayes factors, to precisely decide which one of the two hypotheses is more plausible, for a given measured period ratio and the associated uncertainty.

Our analysis provides a simple and easy to use method to perform inference and model comparison in the presence of incomplete and uncertain information. Measured period ratios and their uncertainties determine the strength of the evidence in favor of a particular hypothesis and, therefore, the support of the performed inferences.

The Bayesian formalism for inference and model comparison is the only fully correct way we have to obtain information about physical parameters and the plausibility of hypotheses from observations (see, e.g., Trotta 2008; von Toussaint 2011) under incomplete and uncertain information. Arregui & Asensio Ramos (2011) performed the first Bayesian parameter inference in coronal seismology. This Letter presents the first application of Bayesian techniques to model comparison. When analytical forward problems are available, both parameter inference and model comparison reduce to a simple and straightforward evaluation of the marginal posteriors and the Bayes factors, which can be obtained by solving simple integrals. The methods outlined in this Letter can be directly applied to most of the seismology inversion problems in which other physical effects and parameters are involved. For instance, additional effects that influence period ratios could be compared to the already considered hypotheses.

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 from the Spanish MICINN/MINECO and FEDER funds through project AYA2011-22846. A.A.R. and A.J.D. acknowledge support by the Spanish MINECO through project AYA2010-18029 (Solar Magnetism and Astrophysical Spectropolarimetry). A.A.R. also acknowledges support from the Consolider-Ingenio 2010 CSD2009-00038 project. We are grateful to the referee and to Ramón Oliver for valuable comments.