© ESO 2019

1. Introduction

Seismology using periods and damping times of transverse standing magnetohydrodynamic (MHD) waves in coronal loops has become a tool for obtaining estimates of physical quantities that are difficult to measure with other means. The first example of the use of the period of a standing MHD wave is by Nakariakov (2000) and Nakariakov & Ofman (2001). Nakariakov & Ofman (2001) used the classic approximate expression for the period P derived from the phase speed v_ph = v_k, with v_k the kink speed, to obtain estimates of the magnetic field strength. This classic expression is obtained for an equilibrium model with a piece-wise constant density in the direction across a constant vertical magnetic field in the thin tube approximation. If the interest is solely in the period as a possible tool for seismology, then the constant density model suffices for a first comparison between observation and theory.

Observations show that these standing MHD waves are damped with damping times of the order of 3 to 10 periods (see e.g. Aschwanden et al. 1999; Nakariakov et al. 1999, and others). The term damping time is the time that is needed to reduce the amplitude of the wave by a factor e. It was realised that resonant absorption offers an explanation for this observed damping. For resonant absorption see for example Ionson (1978), Hollweg & Yang (1988), Goossens et al. (1992), Ruderman & Roberts (2002). However, for resonant absorption to work the equilibrium has to be non-uniform in the cross-field direction. Equilibrium models with a non-uniform intermediate layer were introduced in the seismological investigations. In the non-uniform layer of thickness l the density varies in a continuous manner from its internal value ρ_i to its external value ρ_e. Goossens et al. (2002) used observed values and an approximate expression for the exponential damping time τ_D to compute the thickness of the non-uniform layer, l/R with R the radius of the loop, for 11 events. Arregui et al. (2007; fully numerical) and Goossens et al. (2008; analytical) used observed values of the period P and the exponential damping time τ_D to obtain estimates for the density contrast ζ = ρ_i/ρ_e, the longitudinal Alfvén travel time τ_Ai and the inhomogeneity length scale l/R. Arregui et al. (2007) and Goossens et al. (2008) noted that their seismological inversion problem had ∞¹ solutions. As a matter of fact both Arregui et al. (2007) and Goossens et al. (2008) were too optimistic. The seismological inversion has ∞² solutions but the unknown steepness of the density profile was removed from the inversion scheme by adopting a sinusoidal variation of the equilibrium density.

Resonant absorption is a damping mechanism that works for both standing and propagating waves. Observations by the COronal Multi-channel Polarimeter (CoMP) show that propagating MHD waves are ubiquitous in the solar atmosphere (Tomczyk et al. 2007). This is in stark contrast to standing MHD waves that are rather rare since they need an energetic event as excitation. This triggered interest in the spatial damping of propagating waves by resonant damping. Terradas et al. (2010) derived an analytical expression for the damping length and confirmed its accuracy by numerical simulations. Verth et al. (2010) used the theory developed by Terradas et al. (2010) to successfully explain the disparity between outward and inward wave power of the waves observed with CoMP (Tomczyk & McIntosh 2009).

The damping time and the damping length discussed so far are exponential as the theory adopts a model with wave variables undergoing an exponential decay in time in case of standing waves or in space in case of propagating waves. This comes naturally since the observed damped signals were modelled as exponentially damped harmonic functions. Numerical experiments on propagating waves suggest that the initial damping is better described by a Gaussian rather than an exponential decay. Gaussian damping is first studied for spatial damping of propagating waves by Pascoe et al. (2012). Analytical theory for spatial damping of propagating waves is developed by Hood et al. (2013). Analytical theory for temporal damping of standing waves is developed by Ruderman & Terradas (2013). Pascoe et al. (2013) used analytical theory to do seismology on propagating MHD waves. Pascoe et al. (2016) present a seismological study that uses values of the period, exponential damping time and Gaussian damping time. They make two statements that concern the essence of the seismological inversion problem. First, that previous seismological studies were concerned with ill-posed problems while they deal with a well-posed problem. Second, that their seismological study leads to a unique solution. Both statements have to be interpreted with sufficient caution. If ill posed means that the problem has infinitely many solutions then the first statement is correct. The inversion problem of standing coronal loop oscillations has as a rule infinitely many solutions. As far as the second statement is concerned we note that a seismological study that leads to a single solution has made additional assumptions on the variation of the equilibrium quantities and/or adopted values for specific characteristic quantities. But this single solution can not be unique.

2. Limitations of coronal seismology and ill posed problems

Let us have a look at the limitations of coronal seismology of standing MKD waves. The period P and the exponential and Gaussian damping times, τ_D and τ_G are determined by the equilibrium model. In the β = 0 approximation they are determined by the variation of the local Alfvén frequency. In case of a straight and constant magnetic field it is the variation of density ρ₀(r) and the strength of the constant magnetic field B that determine the period P and the damping times τ_D, τ_G. Inspired by the success of helioseismology we might hope to determine the variation of ρ₀(r) in the interval [0, R]. The piece-wise constant density model with density ρ_i in the interior of the loop and density ρ_e exterior to the loop was and still is popular in solar physics. Unfortunately it excludes resonant damping. It was replaced by a model in which the density ρ₀(r) varies in a continuous manner in a non-uniform layer R − l/2 ≤ r ≤ R + l/2 with thickness l from its interior value ρ_i to its exterior value ρ_e.

The big difference between helioseismology and seismology of coronal loops is that in helioseismology many ( >10⁶) modes are observed while in seismology of coronal loops there is one or occasionally two or more waves observed. The consequence is that seismology of coronal loops is far from being able to determine the radial structure of equilibrium density ρ₀(r). Only a few characteristic quantities related to the equilibrium density can be determined. It is out of the question to find a unique solution for the density distribution ρ₀(r). The classic seismological problem using observations and theoretical expressions for the periods and damping times of transverse standing MHD waves in coronal loops is better referred to as a reduced seismological problem.

The inversion of the reduced seismological problem has more unknowns than equations. The system is underdetermined and has infinitely many solutions. As we shall point out in what follows, there are four unknowns in what should be called a reduced seismological analysis. The unknowns are the longitudinal Alfvén travel time τ_Ai, the density contrast ζ = ρ_i/ρ_e, the radial inhomogeneity length scale l and the steepness α of the density profile at the point of the resonance. This means that there are three equations for four unknowns and hence there is not a unique solution but ∞¹ solutions. Pascoe et al. (2016) prescribe the value of α. They take α = 1 for a linear variation. In the direct problem it is a good first attempt to adopt a linear variation for density to reduce the mathematical complexity in order to understand the physics. In the inverse problem the determination of the equilibrium density is a prime object of investigation and to assume from the start that the variation of density is linear misses somehow the point of seismology. It also affects the inferred value of l/R. In reduced seismology it means that the number of unknowns is reduced from four to three.

3. The direct problem for period and exponential damping

Let us consider the direct problem before focussing on the inverse problem. The direct problem is concerned with the analytical and/or numerical computation of periods and damping times of MHD waves of well defined theoretical equilibrium models. There is no discussion about the equilibrium distribution of the physical quantities and the geometric structure of the equilibrium. They are prescribed. The equilibrium models that have played an important role in the study of MHD waves are the Cartesian slab model and the straight 1D cylindrical model.

In its initial version the Cartesian model has a piece wise constant density with density ρ_i to the left of x = 0 and density ρ_e to the right of x = 0. It has a constant magnetic field in the z-direction and constant plasma pressure. This equilibrium supports surface Alfvén waves. In the incompressible limit the square of the frequency of this surface Alfvén wave is the weighted mean of the squares of the local Alfvén frequencies

(1)

where

(2)

is the square of the local Alfvén frequency. The indices {i, e} refer to x <  0 and x > 0 respectively. Expression (1) is given by for example Chen & Hasegawa (1974a), Ionson (1978); see also Hasegawa & Uberoi (1982). When the restriction of incompressibility is removed expression (1) is recovered for waves with k_y ≫ k_z.

The initial version of the cylindrical model is a straight cylinder with piece wise constant density ρ_i for 0 ≤ r ≤ R and ρ_e for r > R and constant magnetic field in the z-direction. The non-axisymmetric MHD waves with azimuthal wave number m = 1 are of particular importance since they displace the axis and the magnetic tube as a whole. In the long wavelength or thin tube (TT) limit k_zR ≪ 1 the frequency of the fundamental radial mode of the non-axisymmetric waves is given by the same expression as that for the Cartesian surface Alfvén wave.

The fact that the frequency of the fundamental radial mode of the non-axisymmetric waves with m = 1 is equal to that of the surface Alfvén wave of the Cartesian piece wise constant density model is not surprising. Wentzel (1979) suggested that the non-axisymmetric waves of the cylindrical model were surface waves. Goossens et al. (2012) gave a detailed discussion of the fundamental radial mode of the non-axisymmetric kink (m = 1) waves and showed that they are indeed surface Alfvén waves with properties that are the same as those of the surface Alfvén waves of the Cartesian model.

Expression (1) is given by Eq. (30) of Spruit (1982) in the limit k_zR → 0. Expression (1) can be written as a simple expression for the period P that depends on the density contrast ζ = ρ_i/ρ_e and the Alfvén travel time τ_Ai = L/v_Ai (see Goossens et al. 2008). Use the relations to obtain

(3)

with L the length of the loop and . Recall that the expression (3) for the period is for a model with a piece wise constant density. Expression (3) implies that seismology based on observations of the period can only give information about τ_Ai and ζ.

The interest for the surface Alfvén wave was due to the possible heating because of the Alfvén resonance for thermonuclear fusion, application in geophysics and heating of solar and astrophysical plasmas. Early studies on resonant damping at the Alfvén resonance include Chen & Hasegawa (1974b,c,a), Tataronis & Grossmann (1973), Grossmann & Tataronis (1973), Ionson (1978), Hollweg & Yang (1988). See a review in Hasegawa & Uberoi (1982). However, for resonant absorption to work the equilibrium has to be non-uniform. Since

(4)

the frequency of the surface Alfvén wave is in the Alfvén continuum, the surface Alfvén wave undergoes resonant damping. The non-uniformity was introduced into the equilibrium model by replacing the discontinuity at the surface x = 0 in the Cartesian case or at the boundary surface r = R in the cylindrical case by a continuous variation of equilibrium quantities in a transitional layer. The local Alfvén frequency σ_A varies in a continuous manner from its value σ_Ai to its value σ_Ae and defines the Alfvén continuum. This can be done in several ways in the direct problem. Chen & Hasegawa (1974bc,a), Ionson (1978), Hasegawa & Uberoi (1982) introduced a linear variation for the quantity ϵ defined as

(5)

The transitional layer is [ − a/2,  a/2]. When B is constant Eq. (5) implies that the equilibrium density is a linear function in the transitional layer. Hollweg & Yang (1988) take a constant magnetic field and adopt a linear variation for density in the transitional layer. In the Cartesian case the plasma extends to the left of x = 0 to −∞ and to the right of x = 0 to +∞. So there is no radius involved. Hollweg & Yang (1988) in their numerical example introduced R by use of k_y = m/R with m = 1 for the non-axisymmetric kink wave. Tataronis & Grossmann (1973) and Grossmann & Tataronis (1973) adopted a linear variation for the equilibrium density ρ₀(r) and for the square of the axial magnetic field B_z, 0(r) so that the quantity ϵ(r) defined as

(6)

is a linear function of r. In these studies damping rates were obtained by adopting the so-called thin boundary approximation (see e.g. Hollweg & Yang 1988). The studies by Tataronis & Grossmann (1973) and Grossmann & Tataronis (1973) were for incompressible motions. Goossens et al. (1992) studied the damping of axisymmetric (m = 0) (sausage) waves and non-axisymmetric waves. They adopted the thin tube (TT) approximation and the thin boundary (TB) approximation. The use of the thin boundary (TB) approximation l/R ≪ 1, liml/R → 0 makes it possible to obtain an analytical expression for τ_D. The big advantage of the thin boundary approximation is that numerical solutions of the ideal MHD equations are not required but analytical solutions can be used to the left and to the right of the resonant point. Connection formulae (see e.g. Sakurai et al. 1991; Goossens et al. 1995; Tirry & Goossens 1996) are used to connect the ideal solutions across the resonant point for the determination of τ_D. The disadvantage of the TB approximation is that, strictly speaking, it only works for lim l/R → 0.

An important result of Goossens et al. (1992) is their Eq. (77) which gives an expression of the damping rate for damping due to resonant absorption for non-axisymmetric waves (m ≥ 1). When corrected for a typo it is transformed in Eq. (30) of Goossens et al. (2009)

(7)

In Eq. (7) Δ_A(r_A) is

(8)

In the TB approximation the effect of the non-uniform layer on the damping is reduced to the value of Δ_A. Equation (7) was derived without making any assumptions about the variation of the equilibrium quantities in the non-uniform layer. Both density ρ₀(r) and B_z, 0(r) can vary with r. Consider an equilibrium with given values of ρ_i,  ρ_e,  ω_Ai,  ω_Ae. In the thin tube (TT) and thin boundary (TB) approximations the damping rate due to resonant absorption is determined by the derivative of at r_A with r_A the position of the resonant point. In the TB approximation it can be taken as r_A = R. Expression (7) indicates that observations of the damping rate can only give additional information on the value of (dω_A/dr)_rA.

In what follows we assume that B_z, 0 is constant. This assumption simplifies the mathematical expressions and it is a reasonable approximation of reality. For a constant axial magnetic field Eq. (7) can be simplified to (see Eq. (31) of Goossens et al. 2009)

(9)

Expression (9) is remarkably simple. In the case of a straight and constant magnetic field the effect of the non-uniform layer is reduced to the value of the derivative of the equilibrium density at the resonant position

(10)

The variation of density at the position r_A can be characterised by the dimensionless quantity G defined as

(11)

G relates the local variation of density ρ₀ at r = r_A to the global variation of ρ₀ from the value ρ_i to the value ρ_e over the distance R. The expression for the damping rate Eq. (9) (for m = 1) can then be written as

(12)

Now we can use the density contrast ζ = ρ_i/ρ_e and rewrite Eq. (12) as an expression for the damping time τ_D as

(13)

Equation (13) tells us that for a cylindrical straight equilibrium model with density contrast ζ the ratio of the damping time to the period is completely determined by the dimensionless quantity G. This is a remarkable result. All information on the variation of density is collapsed into this one quantity G. This is a big leap. From the point of coronal seismology it implies that the original seismological problem in the present formulation has infinitely many solutions. There is an unlimited number of density distributions that have the same value of the quantity G.

In the direct problem the equilibrium configuration is freely chosen and the non-uniform variation of density can be confined to a layer of thickness l with steepness α. With that prescription

(14)

A linear variation of density corresponds to α = 1. For a sinusoidal variation of density, introduced by Ruderman & Roberts (2002), α = π/2.

With the use of Eq. (14) expression (11) for G becomes

(15)

Equation (15) implies that there are infinitely many couples (α, l/R) that produce the same value of τ_D/P. For example the value for τ_D/P obtained for a linear variation for a layer with thickness (l/R)_L is also recovered for the sinusoidal profile and a layer with thickness (l/R)_L × π/2. Conversely, the value for τ_D/P obtained for a sinusoidal variation for a layer with thickness (l/R)_S is also recovered for the linear profile and a layer with thickness (l/R)_S × 2/π. In the TTTB approximation different models can produce the same ratio τ_D/P. This does not affect the direct problem since α, and l/R are prescribed. However, for the inverse problem this means that it is not possible to distinguish between different couples of (α, l/R).

With the use of expression (15) we can rewrite Eq. (13) as

(16)

Let us recapitulate what we have so far. Equation (3) is an expression for the period of a piece wise constant density model in the long wave length or thin tube approximation. It depends on the density contrast ζ = ρ_i/ρ_e and the Alfvén travel time τ_Ai = L/v_Ai. It does not take into account a possible non-uniform variation of density. Equation (13) has been derived in the thin tube and the thin boundary (TTTB) approximations. It shows that the damping time for exponential damping due to resonant absorption depends on ζ and G, with G the ratio of α and l/R.

In the direct problem the approximations of thin tube and thin boundary can be dropped and numerical solutions can be obtained for the periods, the damping times and eigenfunctions of the resonantly damped MHD waves. For general equilibrium models it is in general no longer possible to obtain analytical solutions, but in the direct problem that does not matter.

Expression (3) for the period neglects the effect of the variation of the density in the non-uniform layer. This effect is not negligible a priori, in particular when the layer is thick. The effect of the non-uniform layer on the period has been studied by Soler et al. (2013).

A common argument to use the thin boundary approximation for the damping time is that the numerical simulations by Van Doorsselaere et al. (2004) and Arregui et al. (2005) show that the dependence of τ_D on the width of the non-uniform layer is not very strong. Arregui et al. (2008) have shown, in the context of prominence threads oscillations, that different combinations of wavelength, density contrast and width of the non-uniform layer lead to percentage deviations up to 20% with respect to the thin boundary approximation for the damping. However, the numerical simulations by Van Doorsselaere et al. (2004), and Arregui et al. (2005, 2008) were for a sinusoidal variation of density. Soler et al. (2013) have shown that other variations of density result in a different dependence on the width. Soler et al. (2014) and Arregui et al. (2015) discussed the consequences of different density profiles on coronal seismology for a scheme using period and exponential damping times.

4. The direct problem for Gaussian damping

Gaussian damping is first studied for spatial damping of propagating waves by Pascoe et al. (2012) in numerical simulations. Analytical theory for spatial damping of propagating waves is developed by Hood et al. (2013). They studied the resonant damping of propagating kink waves in a semi-infinite magnetic tube, where the waves are excited by a driver at the tube end. They found that the wave amplitude is approximately described by a Gaussian function at sufficiently small distances from the tube end. At larger distances they found that the wave amplitude decays exponentially. Also the distance where the transition from the Gaussian to the exponential damping occurs, increases with the width of the non-uniform layer. Finally, they found that the classical theory of resonant absorption overestimates the damping distance and suggests a slower damping. An expression for the Gaussian damping distance L_G for small k_zz is given in their Eq. (49). We note that Hood et al. (2013) use the TTTB approximation and also a linear variation of density in the non-uniform layer (α = 1).

Analytical theory for the temporal damping of standing waves is developed by Ruderman & Terradas (2013). Similarly to Hood et al. (2013), Ruderman & Terradas (2013) use the TTTB approximation and a linear variation of density in the non-uniform layer (α = 1). An analytical expression for the Gaussian damping time τ_G is not given by Ruderman & Terradas (2013). They solve their governing Eq. (30) and compare numerical values of the damping time with the damping time predicted by the classical theory of resonant absorption, τ_D. They find that τ_D exceeds the corresponding numerical values of the actual damping time by about 9% for ϵ = 0.1, about 14% for ϵ = 0.2, and about 18% for ϵ = 0.3. The quantity ϵ here is different from what is used in Eqs. (5) and (6). Here it is the ratio of the thickness of the non-uniform layer l to the radius R: ϵ = l/R. The classical theory of resonant absorption underestimates the actual damping as the classical damping time is longer than the corresponding Gaussian damping time. The error increases when ϵ increases. However, in the view of Ruderman & Terradas (2013) the error is not very large, and the damping time given by the classical theory of resonant absorption can be taken as a reasonable approximation.

Ruderman & Terradas (2013) note that there is an obvious similarity between their results and those obtained by Hood et al. (2013) when distance is interchanged with time. The oscillation amplitude is approximately described by a Gaussian function for sufficiently short time. For standing waves the oscillation amplitude decays exponentially at later times and the time when the transition from the Gaussian to the exponential amplitude profile occurs increases for larger values of ϵ. Finally, the classical theory of resonant absorption underestimates the actual damping.

Analytical studies of Gaussian damping have used the TTTB approximation and also a linear variation of density in the non-uniform layer (α = 1). There are no analytical expressions for the Gaussian damping time τ_G or interpolation formulae for that quantity derived from results of numerical simulations.

5. Seismology for exponentially damped standing MHD waves

Let us have a look at seismology with classical exponential damping. First, consider observed values of period and the corresponding Eq. (3). There is one equation and two unknowns: the density contrast ζ and the Alfvén travel time τ_Ai. Hence there are ∞¹ solutions. The problem is ill-posed since there are more unknowns (two) than observed quantities (one). When the density contrast ζ is prescribed one solution for τ_Ai is found. This is done by Nakariakov & Ofman (2001). They took ζ = 10 and found one solution for τ_Ai. Nakariakov & Ofman (2001) choose a specific value for L and for ρ_i to determine an estimate for B.

Secondly, consider observed values of τ_D and use Eq. (13). There is one equation. The unknowns are the density contrast ζ and the quantity G, with G = α/(l/R), so that there are three unknowns: the density contrast ζ, the radial inhomogeneity length scale, l/R, and the steepness of the density variation α. Hence the reduced seismological inversion scheme has ∞² solutions. The problem is ill-posed because there are more unknowns (three) than observed quantities (one). When the density contrast ζ and α are prescribed one solution for l/R is found. This is done by Ruderman & Roberts (2002) and Goossens et al. (2002). Both investigations used ζ = 10 and assumed that the density variation in the intermediate non-uniform layer is sinusoidal, so that α = π/2. Ruderman & Roberts (2002) considered one numerical example to proof the principle. Goossens et al. (2002) computed l/R for 11 events. The same results as far as the damping times are concerned can be obtained with the use of a linear variation of density and non-uniform layers that are thinner by a factor 2/π than those with the sinusoidal variation.

Thirdly, consider observed values of P and τ_D and use Eqs. (3) and (13). There are two equations. The unknowns are the Alfvén travel time τ_Ai, the density contrast ζ and the quantity G, with G = α/(l/R), so that there are four unknowns: the Alfvén travel time τ_Ai, the radial inhomogeneity length scale, l/R, the density contrast ζ and the steepness of the density variation α. Hence there are ∞² solutions. The problem is ill-posed since there are more unknowns (four) than observed quantities (two). Arregui et al. (2008; fully numerical) and Goossens et al. (2008; analytical) removed one unknown, namely α, from the analysis by adopting a fixed sinusoidal variation of density so that α = π/2. That left them with two equations for three unknowns and hence there are ∞¹ solutions. We note that Arregui et al. (2007) solved the dissipative MHD equations without making the TT approximation and the TB approximation. Again in the context of the TTTB approximations the results of Arregui et al. (2007) and Goossens et al. (2008) can be recovered by the use of a linear variation of density and non-uniform layers that are thinner by a factor 2/π than those with the sinusoidal variation. Actually, other prescriptions of the variation of density; meaning that other values of α can be used together with the corresponding values of l/R.

6. Seismology for Gaussian damped standing MHD waves

Pascoe et al. (2013) used analytical theory to do seismology of propagating waves. Their Fig. 2 shows that Gaussian damping is only important for low density contrasts and for a short distance. For a density contrast of ζ = 3, Gaussian damping is present for only two (or fewer) wavelengths. The statement that the exponential regime is present for long/distances is too strong. It is already there after two (or fewer) wavelengths. As far as Gaussian damping in time is concerned, the last two paragraphs of Sects. 4 and 5 of Ruderman & Terradas (2013) are very helpful. They point out the similarities between spatial and temporal damping. In particular “the oscillation amplitude is approximately described by the Gaussian function at sufficiently small time”. “The damping time given by the classical theory can be used as a reasonable approximation”.

Pascoe et al. (2016) do seismology of standing MHD waves. They introduce a third quantity, the Gaussian damping time τ_G, and claim that the inversion problem is now well posed. This claim is not correct. First, they forget that they are dealing with a reduced seismological problem and that the full seismological problem a fortiori has not a unique solution. Secondly, not even their reduced seismological problem has a unique solution. There are three observed quantities: the period P and the damping times τ_D, τ_G and four unknowns: the Alfvén travel time τ_Ai, the density contrast ζ, the radial inhomogeneity length scale l/R and the steepness of the density profile at the point of the resonance α. Hence there are ∞¹ solutions. In the same manner as Arregui et al. (2007) and Goossens et al. (2008), Pascoe et al. (2016) prescribe the value of α. They take α = 1 for a linear variation. Since α is given a prescribed value, a single solution for the three remaining unknowns is found. A different choice for α would lead to a different solution for the three remaining unknowns.

For the period P and the exponential damping time τ_D Pascoe et al. (2016) use the classic expressions obtained in the TTTB approximation for a linear profile (α = 1). In order to obtain an expression for τ_G Pascoe et al. (2016) use the relation (see their Eqs. (4) and (5))

(17)

with L_G the distance of Gaussian damping of propagating waves and λ the wavelength. For exponential damping the damping time τ_D for standing waves can be related to the damping length L_D of propagating waves by

(18)

with v_gr = ∂ω/∂k_z the group velocity and v_ph = ω/k_z the phase velocity. When the frequency ω varies linearly with k_z then v_gr/v_ph = 1. It is not clear whether this relation (18) also holds for Gaussian damping. Nowhere could we find any evidence in that respect in the literature. So far we have been unable to derive an analytical expression for τ_G.

7. Seismology from Bayesian inference

So far our discussion was focused on the mathematical procedure to obtain algebraic estimates for unknown parameters from observable wave properties. This section provides arguments from the perspective of Bayesian inference. Bayesian inference is increasingly applied to astrophysical problems (Asensio Ramos & Arregui 2018) and is also implemented in the context of coronal seismology (see review by Arregui 2018).

When confronted with the solution of an inverse problem (see Tarantola 1987, for an in-deep discussion on inverse problems), the first difficulty is to find at least one model (in our case, a discrete set of parameter values) that is consistent with the data (a discrete set of observables). The discussion above has shown that this is not possible when the amount of unknown parameters outnumbers that of observables. The situation will of course depend on what exactly one means by consistent with the data and by the amount of information one is willing to accept (as certain) about the particular value of a given model parameter. Recall that Pascoe et al. (2016), in the same manner as Arregui et al. (2008) and Goossens et al. (2008), prescribe the value of α, by adopting a particular profile for the variation of density in the non-uniform layer. A second difficulty, by all means more complicated that the first, is to quantify the uniqueness of the obtained solution. This is the reason why most applied inverse problems across different research areas use optimisation methods trying to find a so-called best fit model.

The Bayesian approach to the solution of inverse problems gets rid of these issues by assuming from the outset that there is no such a thing as a best fit model. By accepting that because the inversion procedure must use information that is incomplete and uncertain the best we can do is to quantify the degree of plausibility of alternative parameter values/models. This degree of plausibility can only be relative. Hence the solution to the inverse problem cannot be unique and is expressed by probability density functions that tell us how the plausibility is distributed over the considered parameter values/models.

This being said, let us now discuss how the Bayesian framework helps us support our mathematical arguments while at the same time offers some additional insight into the relationship between the algebraic and probabilistic results. We do so by computing posterior probability density distributions, p([θ]|M, D), for the vector of unknown parameters [θ], conditional on the assumed theoretical model M and the observed data D. Once the full posterior is computed, marginal posteriors for each unknown in the parameter vector can be computed by marginalisation.

7.1. Seismology from exponential time damping

The Bayesian solution to the problem of inferring the two parameters that define the cross-field density structure (ζ and l/R) from the damping ratio of resonantly damped oscillations is presented by Arregui & Asensio Ramos (2011; by numerical sampling of the posteriors using Markov chain Monte Carlo methods) and Arregui & Asensio Ramos (2014; by direct numerical integration). In both studies, a value of α = π/2 (sinusoidal density profile) is adopted. Arregui et al. (2015) showed how, in strong damping regimes, inference results computed with different profiles for the density in the non-uniform layer can substantially differ. They also showed how obtaining evidence in favour of a particular density model using Bayesian model comparison can be difficult from the practical point of view. The reason is that sufficiently large Bayes factors in support of a given model can only obtained in very strong damping regimes.

We consider the Bayesian inversion of the direct problem given by Eq. (16) for two density models, α = 1 (linear density profile) and α = π/2 (sinusoidal model), by taking a hypothetical observed damping ratio with its measurement error, D = {τ_D/P} = 2.0 ± 0.2, uniform priors for the two unknowns, [θ] = {ζ, l/R}, over the ranges ζ ∈ [1.1, 10] and l/R ∈ [0.01, 2], and a Gaussian likelihood function (as described in Arregui & Asensio Ramos 2014). The top two panels in Fig. 1 display the joint probability density function for ζ and l/R for the case α = 1 (Fig. 1a) and α = π/2 (Fig. 1b). This joint probability density, p({ζ, l/R}|M, D), expresses how the plausibility of different combinations of ζ and l/R is distributed in the two-dimensional parameter space of unknowns. The dark and light grey shaded regions indicate the 68% and 95% credible regions, respectively. They delimit the boundaries of the areas with combinations of ζ and l/R with probabilities larger that those two particular values. It is clear that, although similar, these two plausibility distributions for different values of α are different, that is p({ζ, l/R}| M_α = 1, D) ≠ p({ζ, l/R}|M_{α = π/2}, D). Even if we consider separately each panel (a or b), corresponding to a given α (1 or π/2), the solution is not unique. Some combinations are simply more plausible than others.

Fig. 1. Results from the comparison of inferences using two values of α in Eq. (16). Panels a and b: two-dimensional joint posteriors in the (ζ, l/R) parameter space for α = 1 and α = π/2, respectively. Panels c and d: one-dimensional cuts of the joint posteriors along fixed values of l/R = 0.6 and ζ = 3, respectively. Same line-styles are used to identify cut directions in panels a and b with the corresponding results in panels c and d. The vertical solid lines in panels c and d show the algebraic inversion results using Eq. (16) and the fixed values of l/R = 0.6 in panelc and ζ = 3 in paneld. Panels e and f: marginal posteriors for ζ and l/R, respectively. In these calculations r = τD/P = 2 and σ = 0.1r.

Open with DEXTER

In algebraic seismology using discrete values for observables and parameters it has been common practice to adopt a fixed value for one of the parameters, thus reducing the number of unknowns. For example, one could fix the value of ζ in Eq. (16) and solve for l/R. We can mimic this procedure in our probabilistic inference example by taking cuts of the joint probability density functions displayed in Figs. 1a and b along fixed values of ζ or l/R. By so doing, we obtain the one-dimensional posterior density functions displayed in the two middle panels of Fig. 1. Figure 1c shows posteriors along cuts at a fixed value of l/R, p(ζ|l/R = 0.6, M_α = 1, D) and p(ζ|l/R = 0.6, M_{α = π/2}, D). Fig. 1d shows posteriors along cuts at fixed value of ζ, p(l/R|ζ = 3, M_α = 1, D) and p(l/R|ζ = 3, M_{α = π/2}, D).

Let us first focus on Fig. 1c. By assuming a given value for l/R, two one-dimensional posteriors for ζ are found, one for each value of α. It is clear that they differ, both in central tendency and dispersion. There is shift in the maximum a posteriori estimates. Even if we were to consider only one of the two posteriors in Fig. 1c, the solution is not unique. Some values of ζ are simply more plausible than others.

We can repeat the same arguments based on the results displayed in Fig. 1d. By assuming a given value for ζ, two one-dimensional posteriors for l/R are found, one for each value of α. The two posteriors also differ. There is shift with the maximum a posteriori estimate of the posterior for α = π/2 being a factor π/2 larger than the one corresponding to the posterior for α = 1 (recall our previous discussion after Eq. (15)). Again, the solution is not unique because it depends on the adopted value of α. Even for a given α, some values of l/R are more plausible than others.

The dispersion of the obtained posteriors depends on the assumed uncertainty on the damping ratio. Decreasing the observational error leads to narrower posteriors. In the limit of very small uncertainty in the damping ratio, σ → 0, the posteriors narrow down around the discrete estimates obtained from algebraic inversion of Eq. (16). These discrete estimates for ζ are ζ(τ_d/P = 2, l/R = 0.6, α = 1) = 2.02 and ζ(τ_d/P = 2, l/R = 0.6, α = π/2) = 3.22. For the transverse inhomogeneity length scale, the discrete estimates are l/R(τ_d/P = 2, ζ = 3, α = 1) = 0.41 and l/R(τ_d/P = 2, ζ = 3, α = π/2) = 0.63. Notice the factor π/2 between the two analytical estimates for l/R. The analytical estimates are shown as vertical solid lines in Figs. 1c and d.

The solutions obtained by adopting a given value of ζ or l/R completely remove the uncertainty on the unknown for which a fixed value was taken¹. When little is known about the fixed parameter, this might lead to the false impression that we know more than what we really know for both of them. One of the advantages of the Bayesian approach is that one must specify explicitly what one is willing to accept about the parameters in the problem. The mechanism to do so is by specifying prior probabilities. When some information on one of the parameters is available, one can incorporate this knowledge by using more informative priors, for example a Gaussian function (Arregui & Asensio Ramos 2011). If no additional information is present (or is willing to be accepted), the proper way to infer ζ (alternatively l/R) is to collapse all the information in Figs. 1a and b by marginalising the joint posterior over l/R (alternatively ζ). The resulting marginal posteriors are shown in Figs. 1e and f. They summarise all we can say about the two parameters of interest. The resulting posteriors for the two values of α differ. The difference is more marked for the posterior for l/R (Fig. 1e) than for ζ (Fig. 1f). Once more, even in we consider a fixed value for α, the solution is not unique and the obtained posteriors indicate the different levels of plausibility of all considered alternative solutions. Also, notice that the maximum a posteriori estimates for the posteriors in Figs. 1e and f, do not match the discrete estimates in Figs. 1c and d. The reason is that the full marginal posteriors for a given parameter take into account all the uncertainty on the other parameter, uncertainty that was removed by considering the fixed parameter as certain. In other words, the discrete algebraic estimate or the corresponding Bayesian cut of the joint posteriors are not good estimators of the marginalised posteriors.

Figures 1e and f show how different the solution to the inverse problem can be when considering two alternative density profiles at the non-uniform layer. They were computed for a fixed value of the error on the observable damping ratio. We next computed the same marginal posteriors for ζ and l/R for the two density models and for another damping ratio, r = τ_D/P = 6 representing week damping. The results for both damping ratios (r = 2 and r = 6) are compared in Fig. 2 (we note that panels Figs. 2a and b are the same as panels Figs. 1e and f). The new computations show that the weaker the damping is the less the adopted density model (α) affects the inference results and the more closer the posteriors for the two values of α are. This is in agreement with the conclusion reached by Arregui et al. (2015) who computed Bayes factors trying to gather evidence in favour/against a particular model. Arregui et al. (2015) found that evidence in favour of a given model is easier to obtain in strong damping regimes. Our results show that this is because the posteriors are more differentiated for stronger damping. In any case, the solution to the inverse problem still depends on the adopted value of α and even for a fixed value of α the solution is never unique.

	Fig. 2. Marginal posteriors for ζ and l/R from the comparison of inferences using two values of α and for two values of damping ratio: panels a and b: r = τD/P = 2; panels c and d: r = τD/P = 6. In both cases σ = 0.1r.
Open with DEXTER

7.2. Inference of α and G

We can now move one step further and consider the inversion of the direct problems given by Eqs. (13) and (16) by considering G or α as parameters. Instead of considering that α takes on discrete values arising from the adoption of particular theoretical density models, for example α = 1 for a linear profile; α = π/2 for a sinusoidal profile; or even for a parabolic profile as in Arregui et al. (2015), we can let α and also G to be continuously varying parameters. Notice that from the Bayesian perspective adopting a particular density model and hence fixing the value of α (or G) translates into considering that the uncertainty of α or G is zero, a strong model assumption. The Bayesian approach enables to relax this assumption by considering α or G as parameters and by explicitly imposing priors on them.

We now have two possible problems with one observable, r = τ_d/P. In the first case, with direct problem Eq. (13) and considering G as a new parameter, there are two explicit unknowns (G and ζ) with the additional l/R contained in G. In the second case, with direct problem Eq. (16) and considering α as a new parameter, there are three unknowns (α, ζ, l/R). Solving these two problems can be seen as performing a model comparison application between alternative density profiles with the difference that they are now characterised by the two continuous parameters α and G. In our example applications, we considered α ∈ [0.1, 25] and G ∈ [0.1, 15]. The computed full posteriors p({α, ζ, l/R}|D) or alternatively p({G, ζ}|D) are then marginalised by integrating with respect to the remaining parameters, ζ ∈ [1.1 − 10] and l/R ∈ [0.01−2] in the first case and with respect to ζ only in the second case.

The results are shown in Fig. 3a for the first problem giving the marginal posterior of α and in Fig. 3b for the second problem giving the marginal posterior of G. The different curves correspond to different values of the observable damping ratio, trying to cover week and strong damping cases. The first result is that regardless of whether the inference of α or G is attempted, both can be inferred and show well-constrained probability density functions. In both cases too, the distributions shift towards smaller values of the inferred parameter for stronger damping regimes.

Fig. 3. Marginal posteriors for panels aα and bG for different values of the damping ratio r = τd/P obtained by inversion of Eqs. (16) and (13) and by considering them as additional parameters together with ζ and l/R. In panel a the vertical dashed lines indicate the values α = 1 (linear density model) and α = π/2 (sinusoidal density model). The mean and error at 68% credible intervals for the posteriors in panel a are: for r = 2; for r = 3; for r = 4; and for r = 5. The mean and error at 68% credible intervals for the posteriors in panel b are: for r = 2; for r = 3; for r = 4; and for r = 5. In all computations the error in damping ratio is fixed to σ = 0.2.

Open with DEXTER

Concerning α, in Fig. 3a we over-plotted vertical lines indicating the values of α for the linear and sinusoidal density models. For all the curves shown in Fig. 3a, the maximum a posteriori estimates of α are well above the one corresponding to the values α = 1 or α = π/2. This implies that the relative plausibility of the linear and sinusoidal models is lower than that of other possible models with larger α.

Regarding G, we mentioned before that Eq. (15) implies that there are infinitely many couples (α, l/R) that produce a given damping ratio. The curves in Fig. 3b simply show how the plausibility of those infinite number of values of G is distributed over the considered range.

So far we have been unable to give an interpretation to the large values of α or G at which their posteriors show the largest probability values.

8. Summary and conclusions

The classic seismological problem using observations and theoretical expressions for the periods and damping times of transverse standing magnetohydrodynamic (MHD) waves in coronal loops deals with performing algebraic inversions using discrete values for observables and parameters. In general, only a small number of characteristic quantities of the equilibrium profiles can be determined, hence these problems are better referred to as reduced seismological problems. In particular, there is no reason to adopt a priori a linear or sinusoidal variation of cross-field density. Even if a particular density model is adopted, the reduced seismological problem does not allow a unique solution. In particular, infinitely many couples of (α, l/R) produce the same value for the damping ratio.

These statements were demonstrated by analysing the direct and inverse problems for transverse kink waves in cylindrical straight equilibrium models damped by resonant absorption. The mathematical formulation was generalised by deriving an expression for the damping ratio as a function of a dimensionless quantity, G, that relates the local variation of density at the resonant position to the global variation of the density over the mean radius of the cylinder. All the information on the variation of density is collapsed into this quantity. From the point of coronal seismology this result implies that the original seismological problem has infinitely many solutions.

There is nothing wrong with infinitely many solutions. We only need a tool to quantify their relative plausibility. We further support our arguments by considering the inversion problem in the Bayesian framework. Here the solution can never be unique because the inversion procedure involves the computation of posterior density functions that quantify the level of plausibility between alternative parameters/models. In our example application, the solution to a problem with a fixed density profile (α) with two parameters (ζ and l/R) and one observable (τ_d/P) leads to a two-dimensional posterior distribution for different combinations of the two parameters. When comparing the inference results for two density models, these joint distributions; their cuts along fixed values of one of the parameters and the marginal posteriors of the parameters differ.

Moving one step further, we performed seismology with exponential damping by considering the quantities α and G as additional parameters. The inference results show that well constrained posterior probability density functions can be obtained for them. Of particular relevance is the inference of G, which contains all the relevant information on the variation of density in the non-uniform layer. We found that the posteriors for both α and G depend on the damping ratio and computed them for different values of this observable.

In addition to offering support to our arguments, our Bayesian computations offer insight into the relationship between algebraic solutions using discrete values of observables and parameters and probabilistic plausibility distributions. Seismological inference is not the only source of information about the transverse density profile. Additional information can be extracted from the EUV transverse intensity profile observed by EUV imagers such as SDO/AIA. For a particular event, Pascoe et al. (2018) found that linear and sinusoidal models are the most consistent with the observed EUV intensity profile. Synthetic and observed EUV profiles could be incorporated into a Bayesian analysis to obtain further information.

Acknowledgments

We acknowledge support by the Spanish Ministry of Economy and Competitiveness (MINECO) through projects AYA2014-55456-P (Bayesian Analysis of the Solar Corona) and AYA2014-60476-P (Solar Magnetometry in the Era of Large Solar Telescopes) and FEDER funds.

References

All Figures

	Fig. 1. Results from the comparison of inferences using two values of α in Eq. (16). Panels a and b: two-dimensional joint posteriors in the (ζ, l/R) parameter space for α = 1 and α = π/2, respectively. Panels c and d: one-dimensional cuts of the joint posteriors along fixed values of l/R = 0.6 and ζ = 3, respectively. Same line-styles are used to identify cut directions in panels a and b with the corresponding results in panels c and d. The vertical solid lines in panels c and d show the algebraic inversion results using Eq. (16) and the fixed values of l/R = 0.6 in panelc and ζ = 3 in paneld. Panels e and f: marginal posteriors for ζ and l/R, respectively. In these calculations r = τD/P = 2 and σ = 0.1r.
Open with DEXTER
In the text

	Fig. 2. Marginal posteriors for ζ and l/R from the comparison of inferences using two values of α and for two values of damping ratio: panels a and b: r = τD/P = 2; panels c and d: r = τD/P = 6. In both cases σ = 0.1r.
Open with DEXTER
In the text

	Fig. 3. Marginal posteriors for panels aα and bG for different values of the damping ratio r = τd/P obtained by inversion of Eqs. (16) and (13) and by considering them as additional parameters together with ζ and l/R. In panel a the vertical dashed lines indicate the values α = 1 (linear density model) and α = π/2 (sinusoidal density model). The mean and error at 68% credible intervals for the posteriors in panel a are: for r = 2; for r = 3; for r = 4; and for r = 5. The mean and error at 68% credible intervals for the posteriors in panel b are: for r = 2; for r = 3; for r = 4; and for r = 5. In all computations the error in damping ratio is fixed to σ = 0.2.
Open with DEXTER
In the text