Bayesian Evidence for a Nonlinear Damping Model for Coronal Loop Oscillations

Given a model M with parameter vector θ proposed to explain observed data D, a relational measure of the quality of the model derives from the integral of the joint distribution p( θ , D∣M) over the full parameter space

$$\begin{eqnarray}\begin{array}{rcl}p(D| M) & = & {\displaystyle \int }_{{\boldsymbol{\theta }}}p({\boldsymbol{\theta }},D| M)\,d{\boldsymbol{\theta }}\\ & = & {\displaystyle \int }_{{\boldsymbol{\theta }}}p(D| {\boldsymbol{\theta }},M)\,p({\boldsymbol{\theta }}| M)\,d{\boldsymbol{\theta }}.\end{array}\end{eqnarray} \tag{ 3 }$$

This is the so-called marginal likelihood. In this expression, p( θ ∣M) is a prior distribution over the parameter space and p(D∣ θ , M) is the likelihood of obtaining a specific data realization as a function of the parameter vector. This measure of evidence is relational because a relation is examined to quantify how well the data D are predicted by the model M.

To apply Equation (3) to the damping models M_NL and M_RA a two-dimensional grid over synthetic data space ${ \mathcal D }=(\eta ,{\tau }_{d}/P)$ is constructed. The grid covers the ranges in oscillation amplitude and damping ratio in the observational and theoretical studies by Nechaeva et al. (2019) and Van Doorsselaere et al. (2021). Possible data realizations over the grid are generated using the theoretical predictions given by Equations (1) and (2) for models M_RA and M_NL, respectively. Under the assumption of a Gaussian likelihood function and adopting an error model for the damping ratio alone

$$\begin{eqnarray}&&\begin{array}{l}p({ \mathcal D }| {\boldsymbol{\theta }},{M}_{\mathrm{RA},\mathrm{NL}})\\ \ \ =\frac{1}{\sqrt{2\pi }\sigma }\exp \left\{-\frac{{\left[\frac{{\tau }_{d}}{P}-\frac{{\tau }_{d}}{P}({\boldsymbol{\theta }}){\parallel }_{{M}_{\mathrm{RA},\mathrm{NL}}}\right]}^{2}}{2{\sigma }^{2}}\right\},\end{array}\end{eqnarray} \tag{ 4 }$$

with σ the uncertainty in damping ratio. As for the priors, we choose uniform priors for the unknown parameters over ranges similar to those in Van Doorsselaere et al. (2021), additionally considering the statistical and seismological results by Aschwanden & Peter (2017); Goddard et al. (2017); Pascoe et al. (2018): ${ \mathcal U }(R[\mathrm{Mm}],0.5,7)$, ${ \mathcal U }(\zeta ,1.1,6)$, and ${ \mathcal U }(l/R,0.1,2)$.

Figure 1 shows the resulting marginal likelihoods for the two damping models over the grid of synthetic data in damping ratio and oscillation amplitude. The magnitude of the marginal likelihood at each point over the surface is a measure of how well a particular combination of damping ratio and oscillation amplitude is predicted by each model. This magnitude depends on the functional form of the damping models on their parameters and is not necessarily connected with their linear/nonlinear nature. Certain damping ratio and oscillation amplitude combinations are predicted more often than others. The predictions from each model clearly differ.

Zoom In Zoom Out Reset image size

For nonlinear damping (Figure 1, left panel), the significant magnitudes of marginal likelihood are distributed over the triangular region with right angle in the lower-left corner of the domain. The filled contour plot shows a convex structure with a triangular shape. The area with the highest marginal likelihood corresponds to strong damping regimes with oscillation amplitudes in the range ∼[5, 20] Mm. The obtained marginal likelihood distribution offers a straightforward explanation of the spreading of samples of pairs of oscillation amplitudes and damping ratios in the Monte Carlo analysis by Van Doorsselaere et al. (2021). In addition, it quantifies their relative plausibility. For linear resonant absorption (Figure 1, right panel), the region with the highest marginal likelihood corresponds to low damping ratio values, independently of the oscillation amplitude. The predictive accuracy of the model decreases for weaker damping.

The marginal likelihood values are influenced by the quality of the data, the accuracy of the models and the prior assumptions. The damping times predicted by the nonlinear model M_NL show certain differences with the numerical results of Magyar & Van Doorsselaere (2016). The resonant damping model M_RA describes the long-term asymptotic state of long wavelength oscillations in tubes with thin boundaries. It overestimates the damping for thick layers (Van Doorsselaere et al. 2004; Soler et al. 2014). Linear MHD simulations by Pascoe et al. (2019) show that M_RA significantly overestimates the damping that would occur in the first few cycles of an oscillation when l/R is not small, resulting in a model with an over-preference for low damping ratio values. This means that M_RA has the largest marginal likelihood precisely in the strong damping regime in which it provides the least accurate description of the oscillations, which inevitably overemphasizes the weakness of M_RA and hence the strength of M_NL. Lack of knowledge about the radial density profile (Pascoe et al. 2017b; Arregui & Goossens 2019) is an additional source of inaccuracy. Our prior sensitivity analysis showed that further increasing the upper limit for density contrast to ζ = 9.5 increases the magnitude of $p({ \mathcal D }| {N}_{\mathrm{NL}})$ and the marginal likelihood surface in Figure 1 (left panel) extends toward combinations with slightly larger amplitude and low damping ratio values. Because larger loop radii correspond to a lower nonlinearity parameter η/R, increasing the upper limit for R enables a greater range of τ_d /P to be consistent with a particular observed amplitude.

Overall, there is qualitative agreement between the regions with high marginal likelihood for the nonlinear damping model and the broad location of observed data (compare Figure 1 (left panel) with Figure 6 in Nechaeva et al. 2019). This is suggestive of relational evidence between the observations and the nonlinear damping model.

The marginal likelihood only quantifies the evidence for a model in relation to the data that it predicts. In the general model comparison problem the aim is to assess the relative evidence between alternative models in explaining the same observed data. This is achieved with the calculation of the posterior ratio $p({M}_{\mathrm{NL}}| { \mathcal D })/p({M}_{\mathrm{RA}}| { \mathcal D })$. If the two models are equally probable a priori, p(M_NL) = p(M_RA) and by application of Bayes rule the posterior ratio reduces to the ratio of marginal likelihoods of the two models

$$\begin{eqnarray}&&{B}_{\mathrm{NLRA}}=2\mathrm{log}\frac{p({ \mathcal D }| {M}_{\mathrm{NL}})}{p({ \mathcal D }| {M}_{\mathrm{RA}})}=-{B}_{\mathrm{RANL}},\end{eqnarray} \tag{ 5 }$$

where the logarithmic scale is used for convenience in the evidence assessment.

The Bayes factors B_NLRA and B_RANL defined in Equation (5) are a measure of relative evidence. They quantify the relative plausibility of each of the two models to explain the same data. To assess the levels of evidence the empirical table by Kass & Raftery (1995) is employed. For instance, the evidence in favor of model M_NL in front of model M_RA is inconclusive for values of B_NLRA from 0 to 2; positive for values from 2 to 6; strong for values from 6 to 10; and very strong for values above 10. A similar tabulation applies to B_RANL.

Figure 2 shows the resulting Bayes factor distributions over the ${ \mathcal D }$-space. By construction, regions in this space where B_NLRA and B_RANL reach the different levels of evidence are mutually exclusive and cannot overlap. There is a clear separation between the regions where $p({ \mathcal D }| {M}_{\mathrm{NL}})\gt p({ \mathcal D }| {M}_{\mathrm{RA}})$, thus B_NLRA > 0, from those where $p({ \mathcal D }| {M}_{\mathrm{RA}})\gt p({ \mathcal D }| {M}_{\mathrm{NL}})$, hence B_RANL > 0.

Zoom In Zoom Out Reset image size

Figure 2 (left panel) shows that the evidence supports nonlinear damping in a particular region of data space. The region contains combinations with small oscillation amplitudes in a narrow band below ∼5 Mm and large damping ratio and extends toward combinations with smaller damping ratio and larger oscillation amplitudes in the broader range ∼[2, 23] Mm. Because $p({ \mathcal D }| {M}_{\mathrm{RA}})$ is independent of oscillation amplitude, there is a conformity between the regions with high Bayes factor B_NLRA and marginal likelihood $p({ \mathcal D }| {M}_{\mathrm{NL}})$, although their structures differ. Figure 2 (right panel) shows that the evidence supports resonant absorption in two regions. The largest one extends toward the right-hand side of the domain. Here, the large Bayes factor values in favor of resonant damping are due to the corresponding low values of $p({ \mathcal D }| {M}_{\mathrm{NL}})$ relative to $p({ \mathcal D }| {M}_{\mathrm{RA}})$ (see Figure 1). Because the Bayes factor represents relative strengths, one model being poor at a particular region in data space makes the alternative model seem better in comparison. The evidence in favor of resonant damping is very strong in the region restricted to combinations of very small amplitude oscillations with strong damping, close to the lower-left corner in data space.

Overall, there is qualitative agreement between the regions with high Bayes factor for the nonlinear damping model and the broad location of observed data (compare Figure 2 (left panel) with Figure 6 in Nechaeva et al. 2019). This is suggestive of evidence in favor of the nonlinear damping model relative to linear resonant absorption.

The catalog by Nechaeva et al. (2019) contains information about 223 oscillating loops observed by SDO/Atmospheric Imaging Assembly (AIA) at the 171 Å extreme ultraviolet channel in the period 2010–2018 (see their Table 1). It is an extension of the catalog by Goddard et al. (2016). A subset of 101 cases contains information about both the damping ratio and the oscillation amplitude. The methods described above to compute the marginal likelihood and the Bayes factor are applied to them to assess the strength of the evidence for the nonlinear damping model relative to that for resonant absorption.

Table A1 in the Appendix shows a listing of the cases, with the event and loop identification numbers in Nechaeva et al. (2019); the oscillation period; the damping time; the damping ratio; and the oscillation amplitude. Instead of a grid over synthetic data space, the analysis now considers the specific pairs of measured values for oscillation amplitude and damping ratio for each event, $D={D}_{i}={\{{\eta }_{i},{({\tau }_{d}/P)}_{i}\}}_{i=1}^{101}$. Equations (1) and (2) with the same prior assumptions as before are used to compute the predictions from the models. Predictions are confronted to observations in the computation of the likelihood function in Equation (4) and marginal likelihoods of the models using Equation (3). The ratio of marginal likelihoods in Equation (5) gives the Bayes factor for each case. The results obtained for B_NLRA are given in the rightmost column of Table 1. A positive value implies evidence in favor of the nonlinear damping model. A negative value implies evidence in favor of resonant absorption. The strength of the evidence is assessed with the empirical table by Kass & Raftery (1995).

The marginal likelihood for the nonlinear damping model is larger than the marginal likelihood for resonant damping for the majority of cases, 91 out of 101 with B_NLRA > 0. The opposite happens for the remaining 10 cases with B_NLRA < 0 (B_RANL > 0). The level of evidence is conclusive in 71 cases, with Bayes factors above 2. The level of evidence is inconclusive in 30 cases, with Bayes factors below 2. The number of cases with positive evidence for the nonlinear damping model is 65 (2 < B_NLRA < 6). There is a single case with positive evidence for resonant damping (2 < B_RANL < 6). The number of cases with strong or very strong evidence for resonant damping is 5 (B_RANL > 6). The numerical values of B_NLRA depend on the chosen priors. Increasing the upper limit for ζ slightly favors M_NL while decreasing the upper limit for R tends to favor M_RA.

Figure 3 shows Bayes factor, oscillation amplitude, and damping ratio values for the cases with conclusive evidence supporting either the nonlinear damping model (left panel) or resonant absorption (right panel). The 65 cases with positive evidence for the nonlinear model are grouped in the region in data space where the previous analysis of the structure of the evidence favored nonlinear damping. The six cases with positive or strong evidence for resonant damping are distributed in the complementary region where the low marginal likelihood for nonlinear damping favors resonant absorption.

Zoom In Zoom Out Reset image size

Figure 4 shows the results obtained for all 101 loop oscillations, regardless of the conclusive or inconclusive nature of the evidence. The symbols and their colors indicate the corresponding levels of evidence. The pattern of evidence distribution in the general structure analysis in Section 3.1 is maintained. Among the cases with inconclusive evidence, those with larger marginal likelihood for nonlinear damping (edge-colored red circles) are located at the two sides of the main area with positive evidence for nonlinear damping. Those with larger marginal likelihood for resonant damping (edge-colored blue circles), spread further to the sides, with three of them corresponding to combinations of small amplitudes with strong damping.

Zoom In Zoom Out Reset image size

The amplitudes in the observational data may be consistently underestimated, depending on the orientation of the loop and the distance of the observational slit from the loop apex, because they represent the apparent amplitude projected in the plane of the sky. Since the range of damping ratio values which can be accurately described by M_NL decreases for larger amplitudes (Figure 1, left panel), the effect of underestimating the actual amplitude would generally benefit this model, while having no effect on M_RA.