A&A 484, 851-857 (2008)
DOI: 10.1051/0004-6361:200809728

Analytic approximate seismology of transversely oscillating coronal loops

M. Goossens1 - I. Arregui2 - J. L. Ballester2 - T. J. Wang3


1 - Centre Plasma Astrophysics, Katholieke Universiteit Leuven, Leuven 3001, Belgium
2 - Departament de Física, Universitat de les Illes Balears, 07122 Palma de Mallorca, Spain
3 - Department of Physics, The Catholic University of America and NASA Goddard Space Flight Center, Code 671, Greenbelt, MD 20771, USA

Received 6 March 2008 / Accepted 11 April 2008

Abstract
Aims. We present an analytic approximate seismic inversion scheme for damped transverse coronal loop oscillations based on the thin tube and thin boundary approximation for computing the period and the damping time.
Methods. Asymptotic expressions for the period and damping rate are used to illustrate the process of seismological inversion in a simple and easy to follow manner. The inversion procedure is formulated in terms of two simple functions, which are given by simple closed expressions.
Results. The analytic seismic inversion shows that an infinite amount of 1-dimensional equilibrium models can reproduce the observed periods and damping times. It predicts a specific range of allowable values for the Alfvén travel time and lower bounds for the density contrast and the inhomogeneity length scale. When the results of the present analytic seismic inversion are compared with those of a previous numerical inversion, excellent agreement is found up to the point that the analytic seismic inversion emerges as a tool for validating results of numerical inversions. Actually it helped us to identify and correct inaccuracies in a previous numerical investigation.

Key words: magnetohydrodynamics (MHD) - methods: analytical - Sun: corona - Sun: magnetic fields - Sun: oscillations

  
1 Introduction

References to coronal seismology date back to the 1980s (Roberts et al. 1984) and even the 1970s (Uchida 1970), but coronal seismology remained largely a theoretical concept. This situation changed drastically when observations from space observatories showed that MHD waves are ubiquitous in the solar atmosphere. Opinions might differ, but we are inclined to identify the detection of damped transverse coronal loop oscillations in 1999 by Aschwanden et al. (1999) and Nakariakov et al. (1999) in observations made with the EUV telescope on board of the Transition Region and Coronal Explorer (TRACE) as the real start of coronal seismology. Since then the detection of these oscillations has been confirmed (see e.g. Schrijver et al. 2002). The TRACE oscillations have periods T of the order of $\simeq$2-10 min and comparatively short damping times of the order of $\tau_{\rm d} \simeq 3{-}20$ min. There is general consensus that the TRACE oscillations are fast standing kink mode oscillations. In addition, damped oscillations observed in hot coronal loops by the SUMER instrument on board SOHO have been interpreted as standing slow mode oscillations (Wang et al. 2003,2002) and the measured period has been used to determine the magnetic field strength in the loop (Wang et al. 2007).

Theory shows that when the fast magneto-sonic kink MHD waves have their frequencies in the Alfvén continuum, they couple to local resonant Alfvén waves (Goossens 2008; Wright & Rickard 1995; Goossens et al. 2006; Tirry & Goossens 1996) and get transformed into quasi-modes that are damped by resonant absorption (Tirry & Goossens 1996). This is exactly what happens for a radially stratified cylindrical coronal loop model (Ruderman & Roberts 2002; Goossens et al. 2002). Coupling of the fast magneto-sonic MHD waves to local Alfvén waves is a natural phenomenon in non-uniform coronal loops. This resonant coupling produces a quasi-mode which in a static equilibrium is damped by resonant absorption. The damping is independent of the dissipative coefficients.

This finding does not mean that other damping mechanisms are ruled out. Resonant absorption might not be the only cause of the observed damping but it is definitely operational. Strong support for the robust character of quasi-modes and their resonant damping comes from a recent investigation by Terradas et al. (2008) on MHD waves in multi-stranded coronal loops. An important finding of this investigation is that resonantly damped quasi-modes live on complicated multi-stranded coronal loops. They do not need nice cylindrical magnetic surfaces as might be concluded from studies on simplified 1-dimensional equilibrium models. The message is that 1-dimensional models are a great help to reduce the mathematical complexity but still contain the essential physics of resonantly damped quasi-modes.

Ruderman & Roberts (2002) were the first to suggest that the observed rapid damping of the transverse oscillations of coronal loops could be explained by resonant absorption. In the context of the heating of solar plasmas Hollweg & Yang (1988) have predicted that oscillations in coronal loops are to undergo rapid damping. In the same context Goossens et al. (1992) derived analytical expressions for the frequency and the damping rate of quasi-modes in static and stationary equilibrium models. Ruderman & Roberts (2002) focused on proving the principle of resonant absorption as damping mechanism for the transverse oscillations in coronal loops and considered one specific numerical example. Goossens et al. (2002) looked at the damping times of 11 loop oscillation events and basically confirmed that resonant absorption can explain the observed damping as suggested by Ruderman & Roberts (2002).

The analytic expressions derived by Goossens et al. (1992); Hollweg & Yang (1988); Ruderman & Roberts (2002), and Goossens et al. (2002) for the damping rate are asymptotic in the sense that they are derived in the assumption that the non-uniform layer is thin. This is the so-called thin boundary approximation, in what follows we shall refer to it as the TB-approximation.

The seismological studies on transverse oscillations so far are by Nakariakov (2000), Nakariakov & Ofman (2001), Goossens et al. (2002), Aschwanden et al. (2003), and Arregui et al. (2007). Nakariakov (2000) and Nakariakov & Ofman (2001) used the observed periods and theoretical estimates of the periods based on the long wavelength or thin tube approximation (TT-approximation) for a uniform coronal loop model to derive estimates for the magnetic field strength. The weak link in their analysis is the uncertainties on the density. Goossens et al. (2002) used the observed damping rates and theoretical values of the damping rates based on the TB-approximation to derive estimates for the radial inhomogeneity length-scale. Again, the weak link is the uncertainties on the density. Aschwanden et al. (2003) used the observed damping rates and the damping rates computed by Van Doorsselaere et al. (2004), outside the TB-regime, to determine the density contrast. The first study that used the observational information on both periods and damping times in the context of resonant damping in a consistent manner is by Arregui et al. (2007).

The important finding of Arregui et al. (2007) was that an infinite amount of 1-dimensional cylindrical equilibrium models can reproduce the observed period and damping rate with the internal Alfvén transit time or conversely the internal Alfvén velocity confined to a short range. The study by Arregui et al. (2007) is fully numerical and apparently part of the physics has remained not well understood. The objective of the present paper is to use asymptotic expressions for the period and damping rate to illustrate in a simple and easy to follow manner the process of seismological inversion. The asymptotic expressions are the TT-approximation for the period and the TB-approximation for the damping rate. When both approximations are used simultaneously we shall refer to it as the TTTB-approximation. We are well aware of the fact that in case of strong damping this approximation might give inaccurate results. However, our primary intention is to understand the process of seismic inversion and as we shall see the asymptotic expressions turn out to be accurate far beyond their theoretical range of validity.

2 Asymptotic analytic expressions for period

The analytical expression that we shall use for the period is obtained by (i) adopting the TT-approximation for MHD waves and by (ii) modelling the coronal loop as a uniform cylinder with a straight magnetic field along the z-axis. The TT-approximation means that the period is independent of the radius and that effects due to non-zero radius are absent as far as the period is concerned. The choice of a uniform equilibrium model means that effects due to stratification are absent. The coronal loop is modelled as a cylindrical plasma with constant density $\rho_{\rm i}$ embedded in an external plasma with constant density $\rho_{\rm e}$. The coronal loop is basically a density enhancement with $\rho_{\rm i} > \rho_{\rm e}$. The magnetic field is constant and has the same strength both inside and outside the loop.

Our starting point is the well-known expression for the square of the frequency of the kink mode in a uniform cylinder with a straight magnetic field along the z-axis (see e.g. Edwin & Roberts 1983)

 \begin{displaymath}\omega^2 = \omega_{\rm k}^2 = \frac {\displaystyle \rho_{\rm ... ..._{\rm A,e}^2} {\displaystyle \rho_{\rm i} + \rho_{\rm e}}\cdot \end{displaymath} (1)

The subscripts i and e refer to quantities respectively in the coronal loop and in the external plasma surrounding the loop. $\omega_{\rm A}=k_{z}V_{\rm A}$ is the local Alfvén frequency, with kz the parallel wavenumber, $V_{\rm A}=B/\sqrt(\mu \rho)$ the local Alfvén velocity, and B the magnetic field strength. Hence we can rewrite Eq. (1) as

 \begin{displaymath}\omega^2 = 2 k^{2}_{z} \frac{\displaystyle B^2}{\displaystyle \mu \rho_{\rm e}} \left( 1 + \zeta \right)^{-1}, \end{displaymath} (2)

with $\zeta = \rho_{\rm i}/\rho_{\rm e} > 1$ the density contrast. Now we note that for the observed transverse oscillations, with a wavelength double the length L of the loop, $k_{z} = \pi/L$ and we convert frequencies to periods and rewrite Eq. (2) as (see e.g. Arregui et al. 2008)

 \begin{displaymath}T = \tau_{\rm A,i} \; \sqrt{2} \; \; A(\zeta). \end{displaymath} (3)

Here $\tau_{\rm A,i}=L/V_{\rm A,i}$ is the internal Alfvén travel time and the function $A(\zeta)$ is defined as

 \begin{displaymath}A(\zeta) = \left ( \frac{\displaystyle \zeta +1}{\displaystyle \zeta} \right )^{1/2}\cdot \end{displaymath} (4)

Equation (3) is our first key equation. Let us recall that Eq. (3) has been obtained by use of the TT-approximation, hence effects from non-zero radius and stratification are absent. This equation expresses the period T, which is an observable quantity, in terms of the Alfvén travel time $\tau _{\rm A,i}$ and the density contrast $\zeta $ which are two quantities that we aim to determine from seismic inversion. If we have observed values of the period T and we convince ourselves that Eq. (3) is a good first analytical approximation of the period T then we can invert it for either $\tau _{\rm A,i}$ or $\zeta $. Actually, we shall do both. Let us first solve Eq. (3) for $\tau _{\rm A,i}$. Since we prefer to use dimensionless quantities we introduce y as

 \begin{displaymath}y = \frac{\displaystyle \tau_{\rm A,i}}{\displaystyle T}\cdot \end{displaymath} (5)

From Eq. (3) we obtain

 \begin{displaymath}y = \frac{\displaystyle \tau_{\rm A,i}}{\displaystyle T} = \f... ...splaystyle \zeta}{\displaystyle \zeta + 1} \right )^{1/2}\cdot \end{displaymath} (6)

Since we have not any information on $\zeta $, it might appear that Eq. (6) is not very helpful. However, closer inspection reveals that it contains very valuable information. For a given observed period T, Eq. (6) is a parametric representation of y (or equivalently of $\tau _{\rm A,i}$) in terms of $\zeta $. In order to stress this point we define the function F1 by use of the right hand member of Eq. (6) as

 \begin{displaymath}F_1 \;:\;[1, \; \infty[ \;\rightarrow \mathbf{R},\;\; \zeta \... ...splaystyle \zeta}{\displaystyle \zeta + 1} \right )^{1/2}\cdot \end{displaymath} (7)

For later comparison in Sect. 4 we note that ${\rm d} F_1/{\rm d} \zeta >0$ and ${\rm d}^2 F_1 / {\rm d}\zeta^2 < 0$. Hence $y=F_1(\zeta)$ is a strictly increasing and concave function of $\zeta $. We also note that

\begin{displaymath}F_1(1) = \frac{\displaystyle 1}{\displaystyle 2}, \;\; \lim_{... ...1(\zeta) = \frac{\displaystyle 1}{\displaystyle \sqrt{2}}\cdot \end{displaymath}

This means that
 
$\displaystyle \frac{\displaystyle 1}{\displaystyle 2} \leq \;y < \frac{\displaystyle 1}{\displaystyle \sqrt{2}},$      
$\displaystyle \frac{\displaystyle T}{\displaystyle 2} \leq \;\tau_{\rm A,i} < \frac{\displaystyle T}{\displaystyle \sqrt{2}}\cdot$     (8)

According to inequality (8) the Alfvén travel time is constrained to a narrow range; whatever the density contrast is, the Alfvén travel time is in between $0.5 \;T$ and $0.707 \;T$. If we accept that the density contrast $\zeta $ is not smaller than say 1.5 then $0.5 \;T$ is replaced with $0.548\; T$ narrowing further down the range for $\tau _{\rm A,i}$.

If we are able to determine the length L of the loop, then we can extract $V_{\rm A,i}$ and find

 
    $\displaystyle V_{\rm A,i} = \sqrt{2} \;\; \frac{\displaystyle L}{\displaystyle T} \;\; A(\zeta),$  
    $\displaystyle \sqrt{2} \;\; \frac{\displaystyle L}{\displaystyle T} \; < V_{\rm A,i} \leq \;\; 2 \frac{\displaystyle L}{\displaystyle T}\cdot$ (9)

Hence, here is the narrow range for the local Alfvén velocity that Arregui et al. (2007) found in their numerical inversion. If the density contrast is not smaller than 1.5, $2\; L / T$ is replaced with $1.825 \; L/T$ further reducing the available range for $V_{\rm A,i}$ to less than a 25$\%$ relative margin compared to its possible maximal value.

With the help of the first line of inequality (8) we can refine the definition of F1 and replace (7) with

 \begin{displaymath}F_1 \;:\;[1, \; \infty [ \;\rightarrow \Big[\frac{\displaysty... ...splaystyle \zeta}{\displaystyle \zeta + 1} \right )^{1/2}\cdot \end{displaymath} (10)

Let us now solve Eq. (3) for $\zeta $ and find

 \begin{displaymath}\zeta = \frac{\displaystyle 2 y^2}{\displaystyle 1 - 2 y^2}\cdot \end{displaymath} (11)

Equation (11) is the twin of Eq. (6). For a given observed period T, Eq. (11) is a parametric representation of $\zeta $ in terms of y (or equivalently in terms of  $\tau _{\rm A,i}$). In order to stress this point we define the function G1 by use of the right hand member of Eq. (11) as

 \begin{displaymath}G_1 \;:\;\Big[\frac{\displaystyle 1}{\displaystyle 2}, \;\; \... ...(y) = \frac{\displaystyle 2 y^2}{\displaystyle 1 - 2 y^2}\cdot \end{displaymath} (12)

Since ${\rm d}G_1/{\rm d} y > 0 $and ${\rm d}^2 G_1 / {\rm d} y^2 > 0$, it follows that $\zeta=G_1(y)$ is a strictly increasing and convex function of y. Note now that

\begin{displaymath}G_1(1/2) = 1, \;\;\lim_{y \rightarrow 1/\sqrt{2}} G_1(y)= \infty. \end{displaymath}

With this information on the function G1 we can refine its definition (12). Combined with the definition (10) of F1 we obtain the following prescriptions for the functions F1and G1;
 
$\displaystyle F_1 \;:\;[1, \; \infty[ \;\rightarrow \Big[\frac{\displaystyle 1}... ...{2}} \left( \frac{\displaystyle \zeta}{\displaystyle \zeta + 1} \right )^{1/2},$      
$\displaystyle G_1 \;: \;\Big[\frac{\displaystyle 1}{\displaystyle 2}, \;\; \fra... ...Rightarrow \; G_1(y) = \frac{\displaystyle 2 y^2}{\displaystyle 1 - 2 y^2}\cdot$     (13)

Let us recapitulate what seismic information we have deduced from the observed value of the period. The quantities $\zeta $ and y are in the following intervals

\begin{displaymath}\zeta \in I_{\zeta} = [1, \;\infty[, \;\;\;y \in I_y = \Big[\... ...\;\; \frac{\displaystyle 1}{\displaystyle \sqrt{2}} \Big[\cdot \end{displaymath} (14)

and are related to one another by

 \begin{displaymath}y = F_1(\zeta), \;\;\;\zeta = G_1(y). \end{displaymath} (15)

The functions F1 and G1 are defined by expressions (13). Of course, G1 is the inverse function of F1: G1 = F1-1 and conversely G1-1 = F1. When the period T is known from observations, then there are infinitely many pairs $(\zeta, y)$ that reproduce the observed period. We can let $\zeta $ vary over the interval $[1, \;\infty[$ and compute for each value of $\zeta $ the corresponding value of $y=F_1(\zeta)$, or conversely, let yvary over the interval $[ 1/2, \; 1 /\sqrt{2} [$ and compute for each value of y the corresponding value of $\zeta=G_1(y)$.

3 Asymptotic analytic expressions for damping time

In order for the kink MHD waves to be damped by resonant absorption additional physics has to be introduced in the equilibrium model. The required additional physics is non-uniformity of the local Alfvén velocity. For a constant magnetic field this implies a non-uniform density. Asymptotic expressions for the damping time have been derived by e.g. Hollweg & Yang (1988), Goossens et al. (1992) and Ruderman & Roberts (2002). These asymptotic expressions are derived in the approximation that the non-uniform layer is thin. The true density discontinuity is replaced by a continuous variation in density. Jump conditions are used to connect the solution over the ideal singularity and to avoid solving the non-ideal MHD wave equations. The jump condition for the ideal Alfvén singularity was introduced on an intuitive manner by Hollweg & Yang (1988) and put on a firm mathematical basis by Sakurai et al. (1991), Goossens et al. (1995) and Goossens & Ruderman (1995) for the driven problem, and by Tirry & Goossens (1996) for the eigenvalue problem. The result of this asymptotic analysis is

 \begin{displaymath}\frac{\displaystyle \tau_{\rm d}}{\displaystyle T} = F \;\; \... ... \rho_{\rm e}}{\displaystyle \rho_{\rm i} - \rho_{\rm e}}\cdot \end{displaymath} (16)

Here T and $\tau_{\rm d}$ are the period and the damping time, $\rho_{\rm i}$is the internal density on the axis of the loop and hence in the interval $[0,\; R - \frac{l}{2}]$; $\rho_{\rm e}$ is the constant external density where external refers to $[R + \frac{l}{2}, \; \infty[$, l is the thickness of the non-uniform layer; l/R =0corresponds to a uniform loop and a discontinuous variation in the density at the radius R of the loop. A fully non-uniform loop has l/R = 2 but cannot be studied by an approximate TB expression for the damping time. The numerical factor F depends on the variation in density across the non-uniform layer. For a linear variation as used by Goossens et al. (1992) $F = 4/\pi^2$. Ruderman & Roberts (2002) used a sinusoidal variation in density across the non-uniform layer and found $F = 2/\pi$. In what follows we shall adopt $F = 2/\pi$. A sinusoidal variation in density in the non-uniform transitional layer is probably closer to reality than a linear variation. In addition we shall compare our analytic results with numerical results for fully non-uniform 1-dimensional models of coronal loops obtained by Van Doorsselaere et al. (2004) and by Arregui et al. (2007) by the use of the eigenvalue code LEDA originally designed by Van der Linden (1991). In these studies a sinusoidal variation in density was considered. Note also that in Eq. (16) the period is computed using the TT-approximation.

Rewrite Eq. (16) in terms of $\zeta $ to find

 \begin{displaymath}\frac{\displaystyle \tau_{\rm d}}{\displaystyle T} = \frac {\... ...tyle \zeta - 1} \frac{\displaystyle 1}{\displaystyle l/R}\cdot \end{displaymath} (17)

Equation (17) is the second key equation of the present investigation. Let us now look at this equation from a seismic point of view. If we have observed values of the period T and the damping time  $\tau_{\rm d}$ and we convince ourselves that Eq. (17) is a good first analytical approximation of the damping time divided by period then we can invert Eq. (17) for either $\zeta $ or l/R. Actually we shall do both. Let us first solve Eq. (17) for $\zeta $ and find

 \begin{displaymath}\zeta = \frac{\displaystyle \frac{\displaystyle l}{\displayst... ...rac{\displaystyle \pi \tau_{\rm d}}{\displaystyle T} - 1}\cdot \end{displaymath} (18)

From a seismic point of view, period and damping time can be considered as known from observations and it makes sense to denote $\pi \tau_{\rm d}/ T$ as a constant C. In addition it is convenient to abbreviate l/(2R) as z. Hence

 \begin{displaymath}z = \frac{\displaystyle l}{\displaystyle 2 R}, \;\;C =\frac{\displaystyle \pi \tau_{\rm d}}{\displaystyle T}\cdot \end{displaymath} (19)

Equation (18) can be rewritten as

 \begin{displaymath}\zeta = \frac{\displaystyle C z +1}{\displaystyle C z -1}\cdot \end{displaymath} (20)

Here are several observations to be made. First, since $\zeta >0$ it follows from Eq. (20) that for a given ratio of $\tau_{\rm d}/T$there is a lower bound for the inhomogeneity length scale


 \begin{displaymath}z = \frac{\displaystyle l}{\displaystyle 2 R} > \frac{\displa... ...m d}} = \frac{\displaystyle 1}{\displaystyle C} = z_{\rm min}. \end{displaymath} (21)

Second, realize that Eq. (20) is a parametric representation of $\zeta $ in terms of z= l/(2R). In order to make this point very explicit, we introduce the function G2 defined as

 \begin{displaymath}G_2\;:\;\Big]\frac{\displaystyle 1}{\displaystyle C}, \;\;\; ... ...(z) = \frac{\displaystyle C z + 1}{\displaystyle C z - 1}\cdot \end{displaymath} (22)

It is easy to show that ${\rm d} G_2/{\rm d} z <0$ and ${\rm d}^2 G_2 /{\rm d} z^2 > 0$. Hence $\zeta=G_2(z)$ is a strictly decreasing and convex function of z. It attains its absolute minimum for z = l/(2R) = 1. This minimal value is

 \begin{displaymath}\zeta_{\rm min} = G_2(1) = \frac{\displaystyle C + 1}{\displaystyle C - 1}\cdot \end{displaymath} (23)

Conversely $\zeta $ attains its maximal value of infinity in the limit $z \rightarrow 1/C$. With this information on the function G2 we can refine its definition (22) as follows

 \begin{displaymath}G_2\;:\;\Big] \frac{\displaystyle 1}{\displaystyle C}, \;\;\;... ...(z) = \frac{\displaystyle C z + 1}{\displaystyle C z - 1}\cdot \end{displaymath} (24)

With the minimal value (23) for $\zeta $ we can slightly improve on the lower bound 1/2 for y and T/2 for $\tau _{\rm A,i}$given by inequality (8) and on the upper bound 2 L/T for $V_{\rm A,i}$given by inequality (9). These bounds were found with the sole use of the information on the periods. If we also use the information on the damping rate then we obtain
 
    $\displaystyle y \geq \frac{\displaystyle 1}{\displaystyle 2} \left (\frac{\disp... ...playstyle 2} \left (\frac{\displaystyle C + 1}{\displaystyle C} \right )^{1/2},$  
    $\displaystyle V_{\rm A,i} \leq \frac{\displaystyle 2 L }{\displaystyle T} \left (\frac{\displaystyle C}{\displaystyle C + 1} \right )^{1/2}\cdot$ (25)

With the help of the information on the bounds for y and $\zeta $we can refine the definitions (13) for F1 and G1respectively to their final versions
 
                         $\displaystyle F_1 \;$ : $\displaystyle \;\left[\frac{\displaystyle C + 1}{\displaystyle C - 1}, \;\;\inf... ...;\; \frac{\displaystyle 1}{\displaystyle \sqrt{2}}\Big[ , \right.\right.\right.$  
    $\displaystyle \zeta \; \Rightarrow \; F_1(\zeta) = \frac{\displaystyle 1}{\disp... ...{2}} \left( \frac{\displaystyle \zeta}{\displaystyle \zeta + 1} \right )^{1/2},$  
$\displaystyle G_1 \;$ : $\displaystyle \; \left[\frac{\displaystyle 1}{\displaystyle 2} \left (\frac{\di... ...\displaystyle C + 1}{\displaystyle C - 1}, \;\;\infty\Big[ ,\;\; \right.\right.$  
    $\displaystyle y \; \Rightarrow \; G_1(y) = \frac{\displaystyle 2 y^2}{\displaystyle 1 - 2 y^2}\cdot$ (26)

Let us now solve Eq. (17) for z= l/(2R) and find

 \begin{displaymath}z = \frac{\displaystyle l}{\displaystyle 2 R} = \frac{\displa... ...} \frac{\displaystyle \zeta + 1}{\displaystyle \zeta - 1}\cdot \end{displaymath} (27)

Equation (27) is a parametric representation of z = l/(2R) in terms of $\zeta $. As before we make this point explicit by introducing the function F2 defined as

 \begin{displaymath}F_2 \;:\;\left[\frac{\displaystyle C + 1}{\displaystyle C - 1... ...aystyle \zeta + 1}{\displaystyle \zeta - 1}\cdot\right.\right. \end{displaymath} (28)

It is easy to show that ${\rm d}F_2/{\rm d}\zeta <0$ and ${\rm d}^2 F_2 / {\rm d} \zeta^2 > 0$, hence $z=F_2(\zeta)$ is a strictly decreasing and convex function of $\zeta $. In addition

\begin{displaymath}F_2(\zeta_{\rm min}) = 1, \;\;\lim_{\zeta \rightarrow \infty} F_2(\zeta) = \frac{1}{C}\cdot \end{displaymath}

With this information on the function F2 we can refine its definition (28). Combined with the definition (24) for the function G2 we obtain
 
                            $\displaystyle F_2 \;$ : $\displaystyle \;\Big[\frac{\displaystyle C + 1}{\displaystyle C - 1}, \;\;\; \i... ...le 1}{\displaystyle C} \frac{\displaystyle \zeta + 1}{\displaystyle \zeta - 1},$  
$\displaystyle G_2\;$ : $\displaystyle \;\Big]\frac{\displaystyle 1}{\displaystyle C}, \;\;\; 1\Big] \;\... ...tarrow \; G_2(z) = \frac{\displaystyle C\; z + 1}{\displaystyle C\; z - 1}\cdot$ (29)

Here F2 is the inverse function of G2; F2 = G2-1 and conversely G2 is the inverse function of F2.

Table 1: Left: loop oscillation properties of the analysed events. Right: analytic (A) and numerical (N) inversion results.

  
4 Analytical seismology and comparison with numerical inversion

Let us recapitulate the key results of the previous section. The two quantities that we assume to be known from observations are the period T and the damping time $\tau_{\rm d}$. Analytical theory based on the TTTB-approximation gives us two equations, namely Eqs. (3) and (17) that express the period T and the damping time $\tau_{\rm d}$ in terms of the density contrast $\zeta $, the Alfvén transit time (normalised to the period) $y = \tau_{\rm A,i}/T$ and the inhomogeneity length scale (normalised to the radius of the loop) z = l/(2R). These three quantities $\zeta, y$ and z are the seismic quantities in the sense that they are the quantities that we aim to determine with the use of observed values of the period T and the damping time $\tau_{\rm d}$. Since we have only two equations that relate the three unknown quantities to the two observed quantities there are an infinite number of solutions of the seismic inversion problem, such as first pointed out by Arregui et al. (2007). The seismic variables are constrained to the following intervals

 
$\displaystyle \zeta \in I_{\zeta} = \Big[ \frac{\displaystyle C + 1}{\displaystyle C - 1}, \;\;\infty\Big[$      
$\displaystyle y \in I_y = \Big[\frac{\displaystyle 1}{\displaystyle 2} \left (\... ...le C} \right )^{1/2} , \;\; \frac{\displaystyle 1}{\displaystyle \sqrt{2}}\Big[$      
$\displaystyle z \in I_z = \Big]\frac{\displaystyle 1}{\displaystyle C}, \;\;1\Big]$     (30)

and are related to one another by
 
$\displaystyle y = F_1(\zeta), \;\;\; \zeta \; =\; G_1(y),$      
$\displaystyle z = F_2(\zeta),\;\;\; \zeta \; =\; G_2(z).$     (31)

The functions F1, G1, F2, G2 are defined by (26) and (29).

Of the four functions only two are independent since G1 is the inverse function of F1 and G2is the inverse function of F2. The set (31) gives us the infinitely many solutions of the seismic inversion in parametric form. Each of the three unknowns can be used as parameter and the two remaining unknowns can be expressed in terms of that parameter. For example choose $\zeta $ as parameter. Let $\zeta $take on all values in $I_{\zeta}$ and compute the corresponding values of y and z by the use of $y=F_1(\zeta)$ and $z=F_2(\zeta)$. Or choose y as parameter. Let y take on all values in Iy and then compute the corresponding values of $\zeta $ and z by the use of $\zeta=G_1(y)$ and z = F2(G1(y)). Finally, use z as parameter to define the solutions of the inversion problem. Let z take on all values in Iz and then compute the corresponding values of y and $\zeta $by the use of $\zeta=G_2(z)$ and y = F1(G2(z)).

As an illustrative example we re-analyse loop oscillation event #5 examined by Arregui et al. (2007) in detail using their numerical seismic inversion scheme. The results of the investigation of that loop event are shown on their Fig. 3a. For this loop oscillation event $T = 272 \;\mbox{s}$ and $\tau_{\rm d} = 849 \;\mbox{s}$. The radius R and the length L of this loop are estimated to be $R= 3.65 \;\times 10^6 \; \mbox{m}$ and $L = 1.62 \times 10^8 \;\mbox{m}$respectively. The ratio of the radius to the length of the loop is $R/L \approx 2 \times 10^{-2}.$ This small value is good news for the TT-approximation to the period with effects due to a non-zero radius on period being small.

The constant C= $(\pi\tau_{\rm d})/T$ = 9.81. Hence $\zeta_{\rm {min}} = 1.23, \; y_{\rm {min}} = 0.525, \; y_{\rm {max}} = 0.707$ and $ z_{\rm {min}} =0.102$. The intervals for $\zeta, y, z$ are now

$\displaystyle \zeta \in I_{\zeta} = [1.23, \;\;\infty[,$      
$\displaystyle y \in I_y = [0.525 , \;\;0.707[,$      
$\displaystyle z \in I_z = \;]0.102, \;\;1].$     (32)

 The corresponding interval for $\tau _{\rm A,i}$ is

\begin{displaymath}143 \;\mbox{s} \;\leq \;\tau_{\rm A,i} \;\leq \;192 \;\mbox{s}. \end{displaymath}

If we limit the analysis to $\zeta \geq 1.5$ then the lower bound is $ 150 \;\mbox{s}$ so that

\begin{displaymath}150 \;\mbox{s} \;\leq \;\tau_{\rm A,i} \;\leq \;192 \;\mbox{s}. \end{displaymath}

This can be compared with the interval in Fig. 3a of Arregui et al. (2007) which is (see also their Table 1)

\begin{displaymath}170 \;\mbox{s} \;\leq \;\tau_{\rm A,i} \;\leq \;210 \;\mbox{s}. \end{displaymath}

It is encouraging to note that the overall differences on the Alfvén travel times found here are about $10\%$ although the loop oscillation event is characterised by heavy damping. It is intriguing that the interval as a whole is shifted by about $20\;\mbox{s}$ to shorter Alfvén travel times. This intriguing discrepancy has led us to calculate analytical results for the allowable range of the Alfvén travel time for the remaining 10 loops and to compare them with those obtained by Arregui et al. (2007). It turned out that in all cases there are discrepancies and some of them were far bigger than 10$\%$. This has motivated us to re-examine the numerical results of the investigation of Arregui et al. (2007). It turns out that the values given in Arregui et al. (2007) are inaccurate. We have re-calculated the numerical values for the Alfvén travel time intervals in Arregui et al. (2007) and present them in Table 1. When this issue is taken into account, the corrected numerical interval for the event under consideration is

\begin{displaymath}151 \;\mbox{s} \;\leq \;\tau_{\rm A,i} \;\leq \;187 \;\mbox{s}, \end{displaymath}

which remarkably agrees with the analytic interval. Table 1 shows that analytic and numerical Alfvén travel time intervals agree very well for all events. Analytic lower bounds for $\tau _{\rm A,i}$ are a little below numerical ones. Note that analytic lower bounds for $\zeta $ given in (30) can be slightly lower or larger than the numerical $\zeta=1.5$. When lower, values in parentheses give the analytic $\tau _{\rm A,i}$ for this density contrast. Upper bounds correspond to the limit $\tau_{\rm A,i}(\zeta\rightarrow\infty)=T/\sqrt{2}$, which is nearly approximated at the numerical $\zeta\simeq20$. We have also calculated $z_{\rm min} = 1/C$ for the events studied by Arregui et al. (2007) and the results agree well with the corresponding values of $(l/R)_{\rm min}$ in their Table 1.


  \begin{figure} \par\includegraphics[width=8.5cm,clip]{9728fig1.eps} \end{figure} Figure 1: Analytic inversion (solid lines) and corrected numerical inversion (filled circles) in the ($\zeta $, l/R, $\tau _{\rm A,i}$)-space for loop oscillation event #5 in Table 1.
Open with DEXTER

As an illustration that any of the three seismic quantities can be used as parameter we take z as the parameter. We let z vary and compute for each value of z the corresponding value of $\zeta $ by the use of the function G2defined in (29),

 \begin{displaymath}\zeta = \frac{\displaystyle 9.81 z + 1}{\displaystyle 9.81 z - 1}\cdot \end{displaymath} (33)

Subsequently we compute the corresponding value of y by the use of the function F1 defined in (26),

 \begin{displaymath}y = \frac{\displaystyle 1}{\displaystyle \sqrt{2}} \left( \fr... ...splaystyle \zeta}{\displaystyle \zeta + 1} \right )^{1/2}\cdot \end{displaymath} (34)

Of course, we can only compute $\zeta $ and y for discrete values of z. The results of our computations are summarised on Table 2 and graphically represented in Fig. 1. Recall that the functions F1, G1, F2, G2 (see Eqs. (26) and (29)) are monotonically increasing (F1, G1) or monotonically decreasing (F2, G2) and have concave graphs (F1) and convex graphs ( F2, G1, G2) respectively. This implies f.e. that $y = \tau_{\rm A,i}/T$ is a strictly increasing function of $\zeta $ and a strictly decreasing function of z and conversely that z is a strictly decreasing function of both $\zeta $ and $\tau _{\rm A,i}$. Inspection of the second order derivative of a given quantity with respect to one of the two remaining quantities can tell us that the graph is either concave or convex. For example the graphs of y and of z as function of $\zeta $ are respectively concave and convex. The monotonic variation of $\zeta, y, z$ and the concave or convex appearance of their graphs predicted by our analytic seismic inversion agrees exactly with the behaviour of the numerical inversion. Furthermore, Fig. 1 displays an amazing quantitative agreement between analytic and numerical inversion results.

Table 2: Analytic seismic inversion results for loop #5.

We have not re-analysed in detail loop oscillation event #10 that was examined by Arregui et al. (2007). The inversion for that loop event is shown on their Fig. 3b. A striking property of the solutions is the non-monotonic behaviour of the seismic variables. This is clearly reflected in the pronounced minimum of $\tau _{\rm A,i}$ as function of $\zeta $ and as function of l/R. The decreasing part of $\tau _{\rm A,i}$ as function of $\zeta $ and the increasing part of $\tau _{\rm A,i}$ as function of l/R cannot be recovered by the analytical seismic inversion scheme based on the TTTB-approximation. The analytical TTTB-approximation predicts monotonic variation of the seismic variables and cannot approximate multi-valued solutions with two pairs of $(\zeta, l/R$) corresponding to the same value of $\tau _{\rm A,i}$. The fact that the analytical seismic inversion fails for this loop oscillation event is not disturbing since this event is characterised by extremely heavy damping with $T/\tau_{\rm d} = 0.92$which we anticipated would fall out of the application range of the analytical scheme anyway.

5 Conclusion

In this paper we have presented an analytic approximate seismic inversion scheme based on the TTTB-approximation for computing the period and the damping time. In the TTTB-approximation the period is computed for a uniform loop model in the long wavelength or zero radius approximation. The damping time is computed for relatively weak damping corresponding to thin non-uniform layers. The advantage of this analytical seismic inversion is that it is formulated with the aid of two functions F1 and F2 (and their inverse functions G1 = F1-1 and G2 = F2-1) which are given by simple closed expressions. The practical implementation of the inversion scheme is stunningly simple. The calculations required to obtain solutions can even done with the use of a hand calculator. This analytical scheme seismic inversion clearly shows that the inversion problem has infinitely solutions in the $(\zeta, y, z)$-space as first pointed out by Arregui et al. (2007). It also reveals that the allowable values of y (or Alfvén travel time) are confined to a narrow range. When applied to a loop oscillation event with heavy damping as f.e. loop oscillation event #5 with $T/\tau_{\rm d} = 0.32$ the analytic inversion scheme produces remarkably accurate results. Not only does it recover the overall appearance of the solution curve with the corresponding monotonic behaviour of the seismic variables. In addition, it recovers for a prescribed range of values of $\zeta $ the corresponding values of y (or $\tau _{\rm A,i}$) and z (or l/R).

The disadvantage of the scheme is that (i) it does not take into account the effects of non-zero radius and of radial stratification on the period; (ii) it is strictly speaking only valid for weak damping corresponding to thin non-uniform layers. Corrections due to finite tube radius are of the order of the loop radius to length ratio squared. For the largest observed value of R/L in Table 1 the correction is of the order of 10-3, hence the thin tube approximation does not impose any practical restriction on the applicability of analytical results. As for the thin boundary approximation, for cases with extremely heavy damping, the non-monotonic behaviour displayed by numerical solutions cannot be recovered. This is the price for an analytical scheme. All in all, the accuracy of the results obtained with this analytic inversion is amazing. The final agreement of the analytic seismic inversion with the numerical seismic inversion when the inaccuracies of the numerical inversion are removed is excellent up to the point that the analytic seismic inversion emerges as a tool for validating results of numerical inversions.

Acknowledgements
This research was begun when M.G. was a visitor of the Solar Physics Group at the UIB. It is pleasure for M.G. to acknowledge the warm hospitality of the Solar Physics Group at the UIB and the visiting position from the UIB. M.G. also acknowledges the FWO-Vlaanderen for awarding him a sabbatical leave. IA and JLB acknowledge the funding provided under projects AYA2006-07637 (Spanish MEC) and PRIB-2004-10145 and PCTIB2005GC3-03 (Government of the Balearic Islands). TJW's work was supported by NRL grant N00173-06-1-G033.

References

 
Copyright ESO 2008