BAYESIAN MAGNETOHYDRODYNAMIC SEISMOLOGY OF CORONAL LOOPS

Our model describing resonantly damped transverse loop oscillations contains three free parameters, ζ, l/R, and τ_Ai, that have to be estimated. Let these three parameters be gathered in the vector θ = [ζ, l/R, τ_Ai]^t and let the set d = [P_obs, τ_{d, obs}]^t contain discretized data on period and damping times. When data as a result of measurements determine the values of d, the state of knowledge on θ is represented by p(θ|d), which is given by Bayes' rule (Bayes & Price 1763)

$$\begin{equation} p(\mbox{{\boldmath $\theta $}} | d)=\frac{p(d | \mbox{{\boldmath $\theta $}})p(\mbox{{\boldmath $\theta $}})}{p(d)}. \end{equation} \tag{ 3 }$$

In this expression, p(θ|d) is the so-called posterior probability distribution. The function p(d|θ) is the conditional probability distribution of the data given the parameters. It contains the theoretical relations between parameters and data (including the noise properties) and is a measure of how well the data are predicted by the model. Before d has been observed, it represents the probability distribution associated with possible data realizations for a fixed parameter vector. A posteriori, after observation, p(d|θ) has a very different interpretation. It is the likelihood of obtaining a realization actually observed as a function of the parameter vector and is hence called the likelihood function. Under the assumption that observations are corrupted with Gaussian noise and that they are statistically independent, the likelihood can be expressed as

$$\begin{eqnarray} p(d|\mbox{\boldmath $\theta $}) &=& \left(2\pi \sigma _P \sigma _\tau \right)^{-1}\nonumber\\ &&\times\exp \left\lbrace \frac{\left[P-P^\mathrm{syn}(\mbox{\boldmath $\theta $})\right]^2}{2 \sigma _P^2} + \frac{\left[\tau _{\rm d}-\tau _{\rm d}^\mathrm{syn}(\mbox{\boldmath $\theta $})\right]^2}{2 \sigma _\tau ^2} \right\rbrace, \end{eqnarray} \tag{ 4 }$$

with $P^\mathrm{syn}(\mbox{\boldmath $\theta $})$ and $\tau _{\rm d}^\mathrm{syn}(\mbox{\boldmath $\theta $})$ given by the forward analytical problem of Equation (1). Likewise, σ²_P and σ²_τ are the variances associated with the period and damping times, respectively.

The quantity p(θ) is the prior probability of the model vector. It should contain any prior information we might have on the model parameters, without taking into account the observed data. Three different versions are used in this paper. The normalization constant on the right-hand side of Equation (3) is the marginal likelihood, called "Bayesian evidence" in cosmology applications. The evidence can be written as

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

which is an integral of the likelihood over the prior distribution that normalizes the likelihood and turns it into a probability. This quantity is central for model comparison purposes, since it evaluates the model's performance in the light of data. In Bayesian inference, to compute the posterior distribution one needs to run the model over the entire parameter space of interest and sum that all up.

In the Bayesian formulation the inference of a given parameter set θ is based on the knowledge of p(θ|d), which contains all the information available on θ given the data d and therefore is the solution to the inverse problem. Once the posterior distribution p(θ|d) is known, the position of the maximum value gives the most probable combination of parameters θ that fit the data. Bayes' rule constitutes a simple mathematical formulation of the process of learning from experience. What can be inferred about the model parameters "a posteriori" is a combination of what is known "a priori," independently of the data, and of the information contained in the data. It is an appealing formulation that shows how previous knowledge can be updated when new information becomes available. An additional advantage of the Bayesian approach is that error propagation is consistently taken into account by marginalization. When one needs to know how the data d constrains one of the parameters, say θ_i, it is enough to compute the following integral:

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

The result is known as the marginal posterior, which encodes all information for model parameter θ_i available in the priors and the data. Furthermore, it correctly propagates uncertainties in the rest of the parameters to the one of interest.

The prior choice, p(θ), is a fundamental ingredient of Bayesian statistics, in order to obtain optimal results. The results obtained from the statistical inference should be independent on prior information, provided the prior has a support that is non-zero in regions of the parameter space where the likelihood is large. In this case, repeated application of Bayes' rule will lead to a posterior probability distribution that converges to a common objective inference on the hypothesis (Trotta 2008). Two scientists in the same state of knowledge should assign the same prior, hence the posterior should be identical if they observe the same data. In our particular case, dealing with the inversion of physical parameters in oscillating coronal loops, expressions (2) define the intervals over which the unknown physical parameters are allowed to vary, according to the analytic inversion scheme by Goossens et al. (2008).

When all information we have on unknown parameters is restricted to ranges of variation, a reasonable choice is to assign the same probability to all values contained within that range. This defines a uniform type prior, and we can write

$$\begin{equation} p(\theta _i)=H\big(\theta _i, \theta ^{\mathrm{min}}_i, \theta ^{\mathrm{max}}_i\big), \end{equation} \tag{ 7 }$$

where H(x, a, b) is the top-hat function

$$\begin{equation} H(x, a, b)=\left\lbrace \begin{array}{@{}ll@{}}\frac{1}{b-a}\hspace{22.76228pt} \hbox{$a\le x \le b$},\\ \\ 0 \hspace{31.2982pt} \hbox{otherwise}.\\ \end{array}\right. \end{equation} \tag{ 8 }$$

Another option is to assign a decreasing probability distribution for increasing values of the parameter, over a given range. This can be accomplished by considering a Jeffreys type prior, and we can write

$$\begin{equation} p(\theta _i) = \left[\theta _i \log {\left(\frac{ \theta ^{\mathrm{max}}_i}{\theta ^{\mathrm{min}}_i}\right)}\right]^{-1}. \end{equation} \tag{ 9 }$$

If some additional information about the unknown parameter is available from observations, a Gaussian distribution centered on the measured value can be used, so that

$$\begin{equation} p(\theta _i) = \left(2\pi \sigma _{\theta _{i}}^2\right)^{-1/2} \exp \left[ \frac{-\left(\theta _i-\mu _{\theta _i}\right)^2}{2 \sigma _{\theta _i}^2} \right]. \end{equation} \tag{ 10 }$$

In this expression, the mean $\mu _{\theta _i}$ and the standard deviation $\sigma _{\theta _i}$ would be directly obtained from observations.

In our analysis, we use a uniform prior for the transverse inhomogeneity length scale and the internal Alfvén travel time. The limits for the transverse inhomogeneity are imposed by the modeling, since this parameter must be in the range l/R ∈ [0, 2]. For the internal Alfvén travel time, we have considered uniform prior distributions with minimum and maximum values that enclose the intervals given by I_y in expressions (2), taking into account the observed period. Coronal loops are overdense structures with an internal density that is at most one order of magnitude larger than the external coronal density. Densities two orders of magnitude larger than the coronal density, typical in prominence plasmas, are less likely in coronal loops. Our analysis makes use of the three prior distributions given by Equations (7), (9), and (10), for the density contrast. The virtues and disadvantages of each of them are discussed in Section 4.1. Figure 1 displays an example of the three different types of prior in density contrast. In each case, the integral of the prior probability over the entire range of values must amount to one.

Zoom In Zoom Out Reset image size

The objective of our parameter inference under the Bayesian framework is to sample the full posterior distribution so that, for instance, we can locate the most plausible model that maximizes Equation (3) or calculate the marginalization integrals. To this end, one has to evaluate p(θ|d) for different combinations of parameters. If 10 values for a parameter are to be evaluated and N_par is the number of parameters this means that ${\sim}10^{N_{\rm par}}$ evaluations are needed. This exponential increase of the number of evaluations with the number of parameters is known as the curse of dimensionality. A convenient way to handle this kind of problem is to perform a Markov Chain Monte Carlo (MCMC) simulation, for the sampling of the posterior probability distribution function. The MCMC algorithm we apply allows us to construct a sequence of points in parameter space, called a "chain," whose density is proportional to the posterior distribution function. It therefore provides a method for sampling the posterior distribution up to a multiplicative constant. We use the same code employed by Asensio Ramos et al. (2007) in their inversion of Stokes profiles. The details of the technique and the algorithm are given in that paper, so only the most relevant information is detailed here.

The obtained Markov chain is defined as a sequence of random variables (θ₀, θ₁, ... , θ_n), such that the probability of the θ_i element in the chain only depends on the value of the previous element, θ_{i − 1}. Starting from a given element in the chain, the next point is chosen in such a way that the distribution of the chain asymptotically tends to be equal to the posterior distribution. Markov chains can be shown to converge to a stationary state where successive elements of the chain are samples from the target distribution, in our case the posterior p(θ|d).

Our method uses the Metropolis algorithm (Metropolis et al. 1953; Neal 1993). Starting from a given vector of parameters θ₀, the posterior probability given the data is calculated, p(θ₀|d). This includes the calculation of the priors and the likelihood, which involves the evaluation of the forward problem. Next, a new vector of parameters, θ_i, is obtained sampling from a proposal density distribution, q(θ_i|θ_{i − 1}). Again the posterior is evaluated, p(θ_i|d). The ratio

$$\begin{equation} r=\frac{p(\mbox{{\boldmath $\theta _i$}}|d)q(\mbox{{\boldmath $\theta _i$}}|\mbox{{\boldmath $\theta _{i-1}$}})}{p(\mbox{{\boldmath $\theta _{i-1}$}}|d)q(\mbox{{\boldmath $\theta _{i-1}$}}|\mbox{{\boldmath $\theta _i$}})} \end{equation} \tag{ 11 }$$

is evaluated and θ_i is admitted with probability β = min[1, r]. The process further progresses by obtaining a new vector of parameters from the proposal density distribution.

The key ingredient in this process is the proposal density q(θ_i|θ_{i − 1}). This distribution is used to sample points from the posterior starting from a given point. Ideally, it should be chosen as close as possible to the unknown posterior distribution. Also, it should be simple to sample from it. Since the first condition is impossible to fulfill unless the posterior is known, we choose a sufficiently general distribution function which is easy to sample from. The obvious selection is a multivariate Gaussian distribution. In our case, we allow for a non-diagonal covariance matrix in the Gaussian distribution, which improves the sampling efficiency of the MCMC code (see Appendix A in Asensio Ramos et al. 2007 for more details). We have chosen a uniform distribution for the initial N_unif steps of the chain. Once some information about the posterior is known the algorithm changes to a Gaussian proposal density centered on the current value of the parameter.

In order to assess the performance of our MCMC code, when applied to the problem at hand, the inversion of physical parameters is first performed with synthetic data for period and damping times. These synthetic data are obtained from the solution to the forward problem (expressions (1)) with known values of the loop parameters, so that the numerical inversions are performed under controlled conditions. In addition, this approach enables us to test the performance of the three different prior probability distributions for the density contrast discussed in Section 3.2.

We have considered a model coronal loop with the following physical parameters: ζ = 5, l/R = 0.25, and τ_Ai = 150 km s⁻¹. Forward modeling, according to the theory of resonantly damped kink waves under the thin tube and thin boundary approximations, predicts a period P = 232.4 s and a damping ratio τ_d/P = 3.8. The period and the damping ratio are then used as data for the inversion and the ability of the code to recover the model coronal loop properties is analyzed, for different prior distributions and varying ranges for their definition. In all our simulations, variances of σ_P = 0.1P and $\sigma _{\tau _{\rm d}}=0.1\tau _{\rm d}$ are considered.

Figures 2 and 3 show the inversion results obtained with these data, for fixed prior distributions for the transverse inhomogeneity length scale and internal Alfvén travel time. For the density contrast, the three different prior types defined in Section 3.2 are used in the range ζ ∈ [1.5, 20]. In all three cases, the resulting Markov chains closely follow the same regions outlined by the analytic inversion curve. The distribution and density of the elements in the chain—which directly determine the properties of the posterior—differ in different regions of the parameter space. For a uniform prior distribution in density contrast, the elements in the chain are roughly uniformly distributed along the inversion curve (Figure 2(a)). The relevant quantitative information from the inversion comes from the analysis of the posterior probability distribution functions. Since the elements of the Markov chains are samples from the full posterior it is easy to divide the range for a given parameter into a series of bins and count the number of samples falling within each bin. By computing the marginal posterior for each parameter, we see that the density contrast cannot be constrained (see Figure 3(a)). However, the projection of the chain onto the (l/R, τ_Ai)-plane indicates that a well-defined posterior probability distribution can be obtained for these two parameters. Figures 3(b) and (c) show that, indeed, a very good estimation of these two parameters is possible. When a Jeffreys type prior for the density contrast is used, the elements in the chain are not uniformly distributed. Their density appears to be larger at positions in the parameter space that are close to the model coronal loop properties (Figure 2(b)). The marginal posterior for the density contrast, however, does not show a well-defined probability distribution (Figure 3(a)). The remaining two parameters, though, are well constrained (Figures 3(b) and (c)). Finally, we introduce a hypothetical density contrast measurement for our model coronal loop and use a Gaussian prior distribution for the density contrast centered around that measurement. Figure 2(c) shows that in that case, the inversion produces a Markov chain whose elements are closely packed around the correct loop properties. The marginal posteriors in this case (Figure 3) show well-defined Gaussian-like distributions for the three parameters of interest.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

These results indicate that, in general, the Bayesian inversion scheme will enable us to constrain two out of the three parameters of interest, the transverse inhomogeneity length scale and the internal Alfvén travel time. The density contrast, in general, cannot be constrained. If additional information about the internal and external densities is available from observations, this information can then be used to fully constrain the three unknowns. From the obtained posterior distributions in Figure 3, the median and the computation of the 68% confidence level can be used to obtain estimates with error bars for the inferred parameters.

We found that the obtained results are largely insensitive to the variation of the range assumed for the internal Alfvén travel time, as long as this range is sufficiently wide so it encloses the tails of the resulting posterior distribution. As for the density contrast, additional inversions were performed for different values of the upper limit in the prior distribution for this parameter, considering the three prior types. The results indicate that some prior types perform better than others. Figure 4 displays the joint probability distribution for the internal Alfvén travel time and the transverse inhomogeneity in the form of contours for the 68% and 95% confidence levels, for some illustrative cases. In all plots the known coronal loop parameters are indicated by a symbol. Numerical values for all performed inversions are given in Table 1.

Zoom In Zoom Out Reset image size

When a uniform prior for density contrast is considered, increasing the upper limit of the considered range (Figures 4(a)–(c)) produces a displacement of the joint marginal posterior, in such a way that the inversion underestimates the transverse inhomogeneity length scale and overestimates the Alfvén travel time. Although the real values are for most of the cases enclosed within the error bars, we can see that, for example, for the largest considered range (Figure 4(c)), the real value is already outside the 68% confidence level given by the inversion. Another effect that can be appreciated looking at Figures 4(a)–(c) is the shrinking of the joint posterior for l/R that results in smaller error bars for this parameter (see also Table 1). This result cannot be taken as an improved inversion, as it is solely due to the fact that by extending the range in ζ, the uniform prior gives more weight to elements in the chain with large density contrast, hence lower transverse inhomogeneity length scales and larger Alfvén travel times are obtained. As a matter of fact, the integral of the prior distribution in Equation (7) in the range ζ ∈ [10, 100] is about 10 times the same integral in the range ζ ∈ [1, 10]. The optimal results for the inversion using a uniform prior in density contrast are obtained when the upper limit in ζ is around two or three times the real density contrast value. The problem is that in a real application we can hardly know the real density contrast value.

A better performance and independence of the results from the prior information is obtained when considering a Jeffreys type prior for density contrast. Since this distribution assigns decreasing probability for increasing contrast, the integrals of the prior in Equation (9) in the ranges ζ ∈ [1, 10] and ζ ∈ [10, 100] are of the same order. Figures 4(d)–(f) show the results for this prior. It can be appreciated that the obtained marginal posteriors are almost independent of the assumed range on the prior information, once the upper limit for ζ is sufficiently large, and that optimal results that recover the known coronal loop parameters are obtained (Table 1). For this reason, even if the density contrast itself cannot be constrained, the use of a Jeffreys prior in density contrast is found to be appropriate to perform the inversion of the remaining two parameters in real coronal loops, in the case where information about their density is lacking or uncertain.

When density measurements are available, the Gaussian prior defined in Equation (10) can be used. We have tested the performance of the inversion using this type of prior. As Figures 3, 4(g)–(i), and Table 1 show, in this case a full determination of the three unknowns is obtained that remarkably well recovers the known coronal loop parameters. In addition, the inversion results are independent of the underlying assumptions of the a priori ranges in density contrast. Caution is called to the fact that, in this case, the reliability of the inferred parameters is entirely dependent on the reliability of the density contrast measurement.

Additional inversions under controlled conditions were performed for different combinations for the model coronal loop properties with, e.g., ζ = 2.5, 10, 15 and l/R = 0.5, 0.85. The obtained results do not change in a qualitative way those described above.

Once the performance of our inversion method is evaluated, using synthetic data, we have applied the Bayesian inversion scheme to the 11 loop oscillation events analyzed in, e.g., Ofman & Aschwanden (2002), Goossens et al. (2002, 2008), Aschwanden et al. (2003), and Arregui et al. (2007). Table 2 displays the main oscillation properties for these events and the results from the application of analytic and Bayesian inversion techniques. Estimated values correspond to the median of the obtained distribution and errors are given at 68% confidence level.

Observed period and damping ratios are used as data in our inversion code. Measurement errors in observed periods and damping times for oscillating coronal loops are lacking in the literature. Aschwanden et al. (2002) presented the most complete catalog of events with measurements of periods and damping times for 26 oscillating coronal loops. Uncertainties up to 40%–60% can be found in those events for the measured periods and damping times. The analysis by Van Doorsselaere et al. (2007) drastically reduces the uncertainties to about 1% in period and 3% in damping time. The meaning of such small errors on periods and damping times, with uncertainties that are shorter than the exposure time, is not clear. In our analysis, we consider Gaussian errors of 10% on both quantities. Uniform prior distributions are considered for the transverse inhomogeneity length scale, in the range l/R ∈ [0, 2], and for the internal Alfvén travel time, in a range determined by the oscillation period. We found that the ranges τ_Ai ∈ [1, 400] and τ_Ai ∈ [1, 800] easily accommodate the resulting posterior distributions for the shorter and longer period oscillation events in Table 2, respectively. As for the prior information in density contrast, based on our assessment in Section 4.1, we have performed the inversions using both a Jeffreys prior and a Gaussian prior, centered on the measurements reported by Aschwanden et al. (2003), for the very same 11 loop oscillation events. For the Jeffreys prior the minimum, ζ_min, is computed using the interval for ζ in expressions (2). The maximum value is set to ζ_max = 50, based on our assessment with synthetic results (Table 1). For the Gaussian prior, besides the estimated values we also use the computed error bars (see column "Density ratio" in Table 3 by Aschwanden et al. 2003) to implement the standard deviation in Equation (10) for each case. Uncertainties in density contrast listed by Aschwanden et al. (2003) can be as high as 50%, because of the indirect methods followed to obtain those estimates. As the inversions using Gaussian prior information in density contrast are independent of the assumed ranges, we set the same range as before, ζ ∈ [ζ_min, 50]. The use of a uniform prior in density contrast is discarded, because of its relatively poor performance in front of the Jeffreys prior, in the absence of observational information on this parameter.

Table 2 first lists the analytically obtained intervals for internal Alfvén travel time computed using the corresponding expression for y = τ_Ai/P in expressions (2). This is all the information we can get from that inversion scheme, unless some additional assumption is made. Estimates for the Alfvén travel time using both Bayesian inversions, with Jeffreys and Gaussian priors in density contrast, fit very well within the analytically obtained intervals. These estimates are directly linked to the oscillation period, but in contrast to Nakariakov & Ofman (2001), the inversion now makes use of the information contained in both the period and the damping time. In addition, the Bayesian inversion enables us to compute error bars for this parameter, which are consistently calculated from the uncertainties on data. Error bars in the determination of the Alfvén travel time are found to be rather symmetric.

We next consider the determination of the transverse inhomogeneity length scale. Regardless of the prior type used in the inversion, our estimates for l/R are closely related to the observed damping rate, P/τ_d. Note that the larger this quantity, the stronger the damping by resonant absorption is and the larger the obtained inhomogeneity length scales are. Error bars in the determination of l/R are found to be rather asymmetric, with larger upper errors. The reason is that the posterior for l/R has a more elongated tail toward the right of the maximum, because of the contribution of samples with low values of the density contrast (see the projection onto the [ζ, l/R]-plane in Figure 2(b)). The transverse inhomogeneity length scale is a quantity that is meant to capture the radial variation of the local Alfvén frequency, which is unknown. Our modeling assumes a sinusoidal density variation. This being an assumption, it is not surprising that estimates for l/R are subject to uncertainties. By comparing the estimates of l/R obtained with Jeffreys and Gaussian priors, we find that the numerical estimate considerably differs in some cases, e.g., for loops 4, 5, 9, 10, and 11. When the error bars are taken into account they are rather compatible. More importantly, the concordance between damping ratio and inferred transverse inhomogeneity length scale holds in both cases, in such a way that an ordering of inhomogeneities can be established by just following an ordering of damping ratios.

For the analyzed loop oscillation events, Aschwanden et al. (2003) additionally provide estimates for the transverse inhomogeneity length scale, since the radius (R) and the skin depth (l) are measured (see their Table 2). These measurements lead to values of l/R in between 0.75 and 0.96, which in view of the values displayed on our Table 2 would point to rather more transversely inhomogeneous loops, than those obtained from our inversions. Note that Aschwanden et al. (2003) have typical uncertainties of 20% (up to 50% in some cases) on their estimates of l.

Density contrasts in the last column basically recover the input mean and standard deviation taken from the data measured by Aschwanden et al. (2003) and incorporated in the Gaussian prior. Their value is in the fact that they enable us to collapse the inversion chain, as explained in Section 4.1, and more accurately determine the remaining two parameters, thus producing the smaller error bars in l/R, in comparison with the Jeffreys prior results. As mentioned above, the reliability of the Gaussian inversion, especially in the determination of the transverse inhomogeneity length scale, is entirely dependent on the reliability of the density contrast measurement. For this reason, it is reassuring that estimates of l/R using both prior types—the less informative Jeffreys prior and the more informative Gaussian prior—result in similar values for most of the cases. On the other hand, because density measurements in the solar corona are challenging, we might consider that the inversion results using the less informative Jeffreys prior are more reliable that the ones that make use of uncertain information on density contrast. If that is the case, the differences in the determination of l/R when comparing both Bayesian inversions could be ascribed to inaccurate density contrast measurement.