A&A 430, 1109-1118 (2005)
DOI: 10.1051/0004-6361:20041832

Coronal loop oscillations

Calculation of resonantly damped MHD quasi-mode kink oscillations of longitudinally stratified loops

J. Andries¹ - M. Goossens¹ - J. V. Hollweg ² - I. Arregui¹ - T. Van Doorsselaere¹

1 - Centre for Plasma Astrophysics, K.U. Leuven, Celestijnenlaan 200B, 3001 Leuven, Belgium
2 - Space Science Center, and Institute for the Study of Earth, Oceans and Space University of New Hampshire, Durham, UK

Received 11 August 2004 / Accepted 27 September 2004

Abstract
The observed coronal loop oscillations and their damping are often theoretically described by the use of a very simple coronal loop model, viz. a straight, longitudinally invariant, axi-symmetric, and pressureless flux tube with a different density inside and outside of the loop. In this paper we generalize the model by including longitudinal density stratification and we examine how the longitudinal density stratification alters the linear eigenmodes of the system, their oscillation frequencies, and the damping rates by resonant absorption.

Key words: Sun: corona - Sun: magnetic fields - Sun: oscillations

1 Motivation

In 1999, coronal loop oscillations were observed for the first time by the TRACE spacecraft (Aschwanden et al. 1999). Since then, several oscillating loops have been reported and thoroughly studied (Aschwanden et al. 2002; Schrijver et al. 2002; Nakariakov et al. 1999). For a review on coronal oscillations see e.g. Aschwanden (2004).

The observed coronal loop oscillations have been modelled as fast kink oscillations of a straight longitudinally invariant and axi-symmetric flux tube by e.g. Nakariakov et al. (1999), Ruderman & Roberts (2002) and Goossens et al. (2002) and as phase-mixed torsional Alfvén waves (Ofman & Aschwanden 2002). The rapid damping of the oscillations has been the subject of speculation. Nakariakov et al. (1999) concluded that Reynolds numbers smaller by 8 to 9 orders of magnitude than the classical value of 10¹⁴ are needed to explain the rapid damping. A similar conclusion was drawn by Ofman & Aschwanden (2002) who compared several damping mechanisms and, based on the observed periods and damping times, found phase mixing of torsional Alfén waves to be most likely.

De Pontieu et al. (2001) computed damping rates of Alfvén waves due to footpoint leakage. Under the assumption that their analysis (for Alfvén waves, basically the analysis of Davila (1991)) is also valid for fast waves, De Pontieu et al. (2001) point out that the observed rapid damping can be explained by footpoint leakage within the uncertainties involved in the measurements.

Resonant absorption with coupling to local Alfvén oscillations offers an attractive explanation for the fast damping without the need to invoke substantial changes in the magnetic Reynolds number. The damping of global oscillations by resonant absorption was originally studied in the context of the coronal heating problem (see e.g. Ionson 1978), but in the present context it is rather the damping of the global oscillation that we are interested in. About a decade before coronal loop oscillations were discovered, Hollweg & Yang (1988) predicted that oscillations in coronal loops will be damped in a few oscillation periods due to resonant absorption. This prediction was based on an analysis of quasi-perpendicularly propagating oscillations in a 1D planar slab with a thin non-uniform transitional layer. The analytical result for the damping rate in the original Cartesian geometry was reformulated for a 1D cylindrical loop by relating the Cartesian components of the wave vector to corresponding cylindrical wave numbers. Goossens et al. (1992) derived an analytical expression for the damping rate due to resonant absorption of oscillations in 1D cylindrical flux tubes. However, they did not relate their expression to that obtained by Hollweg & Yang. When the damped coronal loop oscillations were discovered, resonant absorption was suggested as a possible damping mechanism by Ruderman & Roberts (2002). A strong case for resonant absorption was made by Goossens et al. (2002). They used the observed periods and damping rates combined with analytic results for thin loops with thin nonuniform layers to deduce the width of the nonuniform layer for 11 loops. They concluded that resonant absorption had been ruled out by Ofman & Aschwanden (2002), because the latter did not allow for a variation from loop to loop of the radial length scales. However, most of the values for the width of the nonuniform layers are too large for the thin boundary assumption to be an accurate approximation. Goossens et al. (2002) interpreted this as a motivation for an eigenvalue analysis for 1D nonuniform equilibrium states where the non-uniformity is not restricted to a thin layer. Eigenmodes of such highly inhomogeneous loop models have been calculated recently by Van Doorsselaere et al. (2004a). The results of those calculations have been compared with measurements of the inhomogeneity length scales and of the density contrast by Aschwanden et al. (2003). The inhomogeneity length scales are found to be of the order of the width of the loop, and the observed density contrast of the oscillating loops seems to be much lower than that of loops in general. This is in agreement with the computations that show that loops with inhomogeneity length scales of the order of the width of the loop but with a larger density contrast are damped in less than a period and can therefore not be seen to oscillate.

The effect of the loop curvature on the frequencies and damping times of the quasi-mode oscillations has recently been investigated by Van Doorsselaere et al. (2004b)

The main aim of the present paper is to investigate the effect of longitudinal stratification in the equilibrium models on the properties of the coronal loop oscillations. Observations indicate that loops are stretching out into the atmosphere up to heights comparable with the density scale height. Thus, there is a considerable difference in the density at the footpoints and at the loop tops. Together with an almost constant magnetic field this results in a substantial change of the Alfvén velocity along the loop. This in turn, might be anticipated to have a strong effect on the oscillation of the loop.

Our aim here is to use analytic theory to the largest extent possible to understand the effects of longitudinal density variation on the loop oscillations. Part of the method that is used is very similar to that described by Díaz et al. (2002,2001). The present method, however, is more general. Firstly, while Díaz et al. used a discontinuous longitudinal density profile and were thus able to determine the eigenmodes of the local Alfvén operator analytically, we allow a general longitudinal density profile and therefore in general need to obtain the eigenmodes of the Alfvén operator numerically. Secondly, we allow the inclusion of a thin radial transition layer, which causes resonant coupling with local Alfvén waves and hence the damping of the waves by resonant absorption.

In order to make analytical progress we use equilibrium models in which the variation in the radial direction is confined to a thin transitional layer. We are well aware that this approximation is not a very good representation of reality as we have argued in previous papers (Aschwanden et al. 2003; Goossens et al. 2002). However, when variation of the density over the whole loop in both radial and longitudinal directions is included in the equilibrium models, we can only tackle the problem of the damped coronal loop oscillations by numerical means. Hence the emphasis is on analytical insight and we are therefore restricted to models with thin boundaries.

2 Model and linearized equations

The classical model, in which a coronal loop is viewed as a straight cylindrical 1-dimensional flux tube where the equilibrium values only vary with the radial coordinate r, is in this paper extended to a 2-dimensional model where in addition to the radial dependence, longitudinal variation of the density can occur. Thus, in a system of cylindrical coordinates $(r, \varphi, z)$with the z-axis coinciding with the axis of the cylinder (loop), the equilibrium magnetic field $\vec{B} = (0, B_{\varphi}(r), B_z(r))$ and pressure p(r) are functions of the radial distance only, just as in the 1-dimensional case. But for the density, longitudinal variation is taken into account: $\rho(r,z)$. The equilibrium condition is the radial force balance equation

$$\frac{{\rm d}}{{\rm d}r}\left(p+\frac{B^{2}}{2\mu}\right) =-\frac{B_{\varphi}^{2}}{\mu r}\cdot$$ (1)

Just as in the 1-dimensional model the equilibrium condition does not impose any restriction on the density profile. Hence, inclusion of a 2-dimensional density profile in the equilibrium is fairly simple as any density profile may be included. For applications in the lower corona the $\beta=0$approximation can be used, which removes the gas pressure and hence the slow waves form the analysis. In addition we assume a straight magnetic field aligned with the loop $\vec {B}= B(r)\vec{1}_z$ so that Eq. (1) implies a constant magnetic field. The coronal loop is thus a density enhancement in an almost homogeneous field. In addition to the density enhancement in the interior of the loop, the longitudinal density stratification is taken into account in our model.

The observed oscillations are modelled as linear oscillations because the observed velocities are small compared to the local Alfvén speed. This is equivalent with the observation that the displacement is small compared to the length of the loop. As for the oscillation of a string, the fact that the displacement is much larger than the thickness of the loop/string need not be in disagreement with the assumption of linearity.

Since the equilibrium quantities are independent of $\varphi$, the perturbed quantities can be set proportional to $\exp [\imath (m\varphi) ]$ because no coupling between different Fourier modes can occur. Here m (an integer) is the azimuthal wave number. Similarly, we assume a time dependence $\exp[- \imath \omega t]$. For the quasi-modes that are considered the frequency $\omega$ is complex, where $\omega_{\rm r}=\Re (\omega)$ is the oscillation frequency and $\omega_{\rm i}=\Im (\omega)$ is the growth rate. The damping is thus expressed by a negative imaginary part of the frequency. The appropriate theory to justify the proposed time dependence is by means of a Laplace transform rather than a Fourier transform.

These assumptions reduce the governing linearized MHD equations to the following set of partial differential equations in the radial component of the Lagrangian displacement $\xi _r$ and the Eulerian perturbation of the total pressure $p_{\rm T}$:

 

$$L_{\rm A}\frac{1}{r}\frac{\partial r \xi_r}{\partial r} \left(\frac{m^2}{r^2}-\frac{\mu}{B^2}L_{\rm A}\right)p_{\rm T}$$ =

$$\frac{\partial p_{\rm T}}{\partial r} L_{\rm A}\xi_r,$$ = (2)

where we have defined the Alfvén operator:

$$L_{\rm A}=\rho\omega^2+\frac{B^2}{\mu}\frac{\partial^2}{\part... ...ft(\omega^2+v_{\rm A}^2\frac{\partial^2}{\partial z^2}\right),$$

with $v_{\rm A}^2=B^2/(\rho\mu)$ the square of the Alfvén speed. For the azimuthal component of the displacement $\xi _\varphi$ we have the additional equation:

$$L_{\rm A}\xi_\varphi=\imath\frac{m}{r}p_{\rm T},$$

while the longitudinal component of the Lagrangian displacement $\xi_z$ vanishes due to the removal of the slow waves by the $\beta=0$ assumption.

  
3 Derivation of the dispersion relation

  
3.1 Eigenmodes of the Alfvén operator

At any radial location and for a fixed value of the frequency we can search for the set of eigenmodes $\psi^{(k)}(z)$ of the local Alfvén operator:

$$L_{\rm A}\psi^{(k)}(z)=\lambda_k\psi^{(k)}(z) .$$

Notice that these eigenmodes of the Alfvén operator do not directly represent Alfvén waves. But for certain frequencies $\omega_{{\rm A},n}$ the Alfvén operator may have a vanishing eigenvalue. The local Alfvén waves of a certain magnetic surface are described by the set of frequencies $\omega_{{\rm A},n}$ for which the Alfvén operator has a vanishing eigenvalue $\lambda_{k_n}=0$ (there are non-trivial elements in the kernel) and by the corresponding set of eigenmodes (the non-trivial elements of the kernel) $\phi^{(n)}(z)=\psi^{(k_n)}(z)$ since:

$$\mbox{\boldmath {$\xi$ }}(r,z)=\delta(r-r_0)\phi^{(n)}(z)\mathbf{1}_\varphi$$

is indeed a solution because $L_{\rm A}(r_0,\omega_{{\rm A},n})\xi_\varphi=0$.

More generally, however, the eigenmodes $\psi^{(k)}(z)$ for a fixed value of the frequency form a complete set of orthogonal functions in the z direction in the sense:

$$\left\langle \psi^{(k)}\mid\psi^{(l)}\right\rangle\equiv\frac... ..._0^L\overline{\psi^{(k)}(z)}\psi^{(l)}(z){\rm d}z=\delta_{k,l}$$ (3)

where the bar denotes the complex conjugation. Practically, all functions in the z direction can be expressed as a sine series and in particular:

$$\psi^{(k)}(z)=\sum_{n=1}^{+\infty}\psi^{(k)}_n\sin\left(\frac{n\pi}{L}z\right).$$

Hence, we speak of the eigenmodes of the Alfvén operator in terms of the infinite sine series coefficient vectors $\psi^{(k)}_n$. In that formalism the scalar product defined in Eq. (3) reduces to the classical complex scalar product for vectors (the factor 2 in definition (3) is included so that the sine functions have unit length). Likewise the Alfvén operator can be expressed as a matrix operator. In order to determine the expression of the Alfvén matrix operator, the stratified density profile is expressed as a sine series plus a constant term taking into account the density at the footpoints:

$$\rho(r,z)=\rho_0(r)\left [1+\sum_{n=1}^{+\infty}\alpha_n(r)\sin\left(\frac{n\pi}{L}z\right)\right ].$$

Notice that in principle the $\alpha$ parameters can depend on the radial coordinate as the stratification is not necessarily the same inside and outside the loop.

The sine series expansion of the multiplication of two sine functions is:

$$\sin\left(\frac{n\pi}{L}z\right)\sin\left(\frac{k\pi}{L}z\rig... ... \sum_{l=1}^{+\infty}S_{n,k,l}\sin\left(\frac{l\pi}{L}z\right)$$

where S_n,k,l=0 when n+k+l even, but for n+k+l odd:

$$\begin{eqnarray*}S_{n,k,l}&=&\frac{2}{L}\int_0^L \sin\left(\frac{n\pi}{L}z\right... ...}\left [\frac{1}{k^2-(l-n)^2}-\frac{1}{k^2-(l+n)^2}\right ]\cdot \end{eqnarray*}$$

Note that the ordering of the indices does not matter.

The density matrix operator thus becomes:

$$\rho=\rho_0 \left(I+\sum_{n=1}^{+\infty}\alpha_n S_n\right)\equiv\rho_0\widehat{S}$$

where I is the unity matrix, and the matrix $S_n=\left(S_{n,k,l}\right)_{k,l}$. The Alfvén matrix operator therefore becomes:

$$L_{\rm A}=\rho_0\omega^2\widehat{S}-\frac{B^2}{\mu}\partial_z... ...rho_0\left(\omega^2\widehat{S}-v_{\rm A0}^2\partial_z^2\right)$$ (4)

with $v_{\rm A0}$ the Alfvén speed at the footpoints of the field lines and $\partial_z^2$ the diagonal matrix:

$$(\partial_z^2)_{k,l}=\left(\frac{k\pi}{L}\right)^2\delta_{k,l}.$$

The eigenmodes of the Alfvén operator can thus be found as the eigenmodes of an infinite matrix problem. The structure of the S matrices implies that symmetric stratification ( $\alpha_n=0$ for n even), which is approximately the case for loops, only causes coupling within the symmetric modes and within the antisymmetric modes but not between symmetric and anti-symmetric modes, which could also be deduced directly from the symmetry of the problem.

In the absence of longitudinal stratification, the matrix $\widehat{S}$ is simply the unity matrix and the Alfvén operator reduces to a diagonal matrix. The eigenmodes are thus single sines with a longitudinal wavenumber $k_z=\frac{k\pi}{L}$.

  
3.2 Fast eigenmodes of radially uniform loops

Since the eigenmodes of the Alfvén operator form a complete set we can express the variables $\xi_r(r,z)$ and $p_{\rm T}(r,z)$ at any radial location as a superposition of those eigenmodes:

$$\begin{eqnarray*}\xi_r(r,z)=\sum_{k=1}^{+\infty}X^{(k)}(r)\psi^{(k)}(z)\\ p_{\rm T}(r,z)=\sum_{k=1}^{+\infty}P^{(k)}(r)\psi^{(k)}(z) . \end{eqnarray*}$$

When the density is radially independent, the Alfvén operator, as well as its eigenmodes $\psi^{(k)}(z)$ and eigenvalues $\lambda_k$, are radially independent. Using the orthogonality of the eigenmodes, it is straightforward to obtain that the set of partial differential Eqs. (2) is translated to a set of two ordinary differential equations in the radial direction for X^(k)(r) and P^(k)(r) for each k separately:

$$\begin{eqnarray*}\lambda_k\frac{1}{r}\frac{\partial r X^{(k)}}{\partial r}&=&\le... ...{(k)}\\ \frac{\partial P^{(k)}}{\partial r}&=&\lambda_kX^{(k)}. \end{eqnarray*}$$

From these two equations one second order ordinary differential equation in P^(k) can be deduced (Bessel's equation):

$$\frac{{\rm d}^2P^{(k)}}{{\rm d}r^2}+\frac{1}{r}\frac{{\rm d}P... ...}}{{\rm d}r}- \left(\frac{m^2}{r^2}+\kappa_k^2\right)P^{(k)}=0$$ (5)

with:

$$\kappa_k^2=-\frac{\lambda_k}{B^2/\mu}\cdot$$

If the coronal loop is modelled as an internal region and an external region that are radially uniform and separated by a boundary surface r=R, then solutions in both the internal and external (radially) uniform regions can be obtained easily. Inside the loop we have to impose regularity at the axis and we thus find the solution:

$$\begin{eqnarray*}p_{\rm T}^{({\rm in})}(r,z) &=& \sum_{k=1}^{+\infty}A^{({\rm in... ...^{({\rm in},k)}I'_m(\kappa_{{\rm in},k} r)\psi^{({\rm in},k)}(z) \end{eqnarray*}$$

where I_m is the modified Bessel function of the first kind of order m and the prime denotes the derivative with respect to its argument. The $A^{({\rm in},k)}$ are arbitrary coefficients of each of the Alfvén eigenmodes in the solution.

For the external (radially) homogeneous region we impose that the solution should vanish at infinity so that:

$$\begin{eqnarray*}p_{\rm T}^{({\rm ex})}(r,z) &=& \sum_{k=1}^{+\infty}A^{({\rm ex... ...^{({\rm ex},k)}K'_m(\kappa_{{\rm ex},k} r)\psi^{({\rm ex},k)}(z) \end{eqnarray*}$$

where K_m is the modified Bessel function of the second kind of order m and the prime denotes the derivative with respect to its argument. The root $\kappa_k$ has to be taken so that it has a positive real part in order to ensure that the solution vanishes at infinity. If $\kappa_k^2$ is negative the root has to be taken so that the external wave is outgoing. In that case the solution is better expressed in terms of Hankel functions. We will not go into details about this since in our applications $\kappa_k^2>0$ in the external medium. The $A^{({\rm ex},k)}$ are again arbitrary coefficients of each of the Alfvén eigenmodes in the solution. The merit of the splitting into eigenmodes of the Alfvén operator is thus that the solution, in which the radial and longitudinal variables are not separable, is now expressed as a sum of separable terms.

The dispersion relation is obtained by requiring that at the boundary surface r=R both $\xi _r$and $p_{\rm T}$ are continuous (in the radial direction):

$$\begin{eqnarray*}p_{\rm T}^{({\rm ex})}(R,z) &=& p_{\rm T}^{({\rm in})}(R,z)\\ \xi_r^{({\rm ex})}(R,z) &=& \xi_r^{({\rm in})}(R,z). \end{eqnarray*}$$

The variables $\xi _r$ and $p_{\rm T}$ can also be expressed as a sine series expansion:

$$\begin{eqnarray*}p_{\rm T}(r,z) &=& \sum_{n=1}^{+\infty}P_n(r)\sin\left(\frac{n\... ... &=& \sum_{n=1}^{+\infty}X_n(r)\sin\left(\frac{n\pi}{L}z\right). \end{eqnarray*}$$

Practically, the matching condition is then expressed as an infinite set of matching conditions, one for each coefficient in the sine series of $p_{\rm T}(R,z)$ and one for each coefficient in the sine series of $\xi_r(R,z)$. Thus, an infinite set of equations for the variables $A^{({\rm in},k)}$ and $A^{({\rm ex},k)}$ is obtained:

$$\left( \begin{array}{ccccc} \Pi^{({\rm ex},1)}_1 & -\Pi^{({\r... ...rm ex},2)} \\ A^{({\rm in},2)} \\ \vdots \end{array}\right)=0$$ (6)

where $\Pi^{(k)}_j$ is due to the contribution of the jth sine series coefficient of the kth eigenmode in the pressure perturbation at r=R, and $\Xi^{(k)}_j$ is due to the contribution of the jth sine series coefficient of the kth eigenmode in the Lagrangian displacement at r=R. Thus:

$$\begin{eqnarray*}\Pi^{({\rm in},k)}_j &=& I_m(\kappa_{{\rm in},k} R)\psi^{({\rm ... ...a_{{\rm ex},k}}K'_m(\kappa_{{\rm ex},k} R)\psi^{({\rm ex},k)}_j. \end{eqnarray*}$$

The dispersion relation is then given by requiring that the system (6) has non-trivial solutions i.e. its determinant is zero.

In the absence of longitudinal stratification, $\psi^{(k)}_j=\delta_j^k$, so that the matrix of the dispersion relation contains zeroes everywhere except for in $2\times 2$ blocks along the diagonal. The determinant of such a block diagonal system can easily be reduced so that the dispersion relation simplifies to:

$$\prod_j \left\vert \begin{array}{cc} \Pi^{({\rm ex},j)}_j & -... ...rm ex},j)}_j & -\Xi^{({\rm in},j)}_j \end{array}\right\vert =0$$

where each factor yields the dispersion relation for a wave with a specific longitudinal wave number k_z which simplifies to the familiar impedance-matching expression:

$$\frac{\Xi^{({\rm ex},k)}_k}{\Pi^{({\rm ex},k)}_k}=\frac{\Xi^{({\rm in},k)}_k}{\Pi^{({\rm in},k)}_k}\cdot$$

  
3.3 The thin boundary approximation for the fast quasi-modes

When a transition layer of thickness d, in which the density changes continuously in the radial direction, is included in between the radially uniform internal and the external region, the set of Eqs. (2) should be integrated inside the non-uniform layer. As the eigenmodes and eigenvalues of the Alfvén operator are radially dependent in this region it does not help to express the solutions as a superposition of eigenmodes of the Alfvén operator. Instead, in this region we could immediately express the solutions in a sine series. In this way the set of partial differential Eqs. (2) is translated into an infinite set of coupled ordinary differential equations in the radial direction in the sine series coefficient vectors P_n and X_n. When cut down to finite dimensions these equations can be easily integrated numerically. However, just like in the 1D problem, such an integration method fails around certain radial positions in the boundary layer where the equations are singular. This occurs at positions where the Alfvén operator has non-trivial elements in its kernel i.e. it has a vanishing eigenvalue. When expressed as an infinite set of coupled ordinary differential equations the Alfvén operator is a matrix operator, which has a vanishing determinant at resonant positions. The fast wave resonantly interacts with the Alfvén wave corresponding to the non-trivial element of the local Alfvén operator. At those resonant positions the solutions become unbounded and thus cannot be tracked numerically. In a thin layer where the solutions become large, dissipation becomes important and inclusion of dissipative terms in this layer regularizes the solutions. However, the characteristic jump conditions for the radial displacement $\xi _r$ and the Eulerian perturbation of the total pressure $p_{\rm T}$ over the resonant (dissipative) layer can in fact also be retrieved in ideal MHD by analytical continuation of the Greens function. For the 2D problem, expressions for the local ideal solutions around the resonances were obtained by Thompson & Wright (1993). As in the 1-dimensional case a jump in $\xi _r$ occurs in the 2D problem, where, not surprisingly, the jump is z-dependent and exactly proportional to the non-trivial element in the kernel of the local Alfvén operator $\phi(z)$:

$$[\![\xi_r(z)]\!]=-\imath \pi \mbox{sign}(\omega_{\rm r})\frac... ...{\vert \langle\phi\vert L_{{\rm A}1}\phi\rangle\vert}\phi(z) ,$$

where we have denoted the radial derivative of the Alfvén operator as $L_{{\rm A}1}$:

$$L_{{\rm A}1}=\frac{\partial L_{\rm A}}{\partial r}=\frac{\par... ...k}\right]+\rho\omega^2\frac{\partial\alpha}{\partialr}S_{k,k}.$$

Notice that $\langle\phi\vert p_{\rm T}\rangle=P^{({\rm res})}\langle\phi\vert \phi\rangle$ is proportional to $P^{({\rm res})}$, the coefficient of $\phi(z)$ when $p_{\rm T}(z)$ is expressed as a superposition of the eigenmodes of the local Alfvén operator. Just as in 1D it can be shown that $P^{({\rm res})}$ is approximately constant near the resonant position.

These jumps around the resonances are the main contributions of the continuous boundary layer provided that it is thin, i.e. its width d is much smaller than the radial length scale of the solutions. The thin boundary approximation then consists of including the jump contribution of the boundary layer into the dispersion relation. Dispersion relation (6) then has to be modified by replacing the internal contributions by the internal contribution plus the jump contribution:

$$\Xi^{({\rm in},k)}_j \longrightarrow \Xi^{({\rm in},k)}_j+[\![\xi_r(z)]\!]^{(k)}_j$$

with:

$$[\![\xi_r(z)]\!]^{(k)}_j = -\imath \pi \mbox{sign}(\omega_{\r... ...^{(k)}}{\vert \langle\phi L_{{\rm A}1}\phi \rangle\vert}\phi_j$$

where we have used and defined:

$$\begin{eqnarray*}\left\langle\phi\vert p^{({\rm in})}_{\rm T}\right\rangle &=& \... ...sum_k \left(\overline{\phi}\Pi^{({\rm in})}\right)^{(k)}A^{(k)}. \end{eqnarray*}$$

Notice that in absence of longitudinal stratification, $[\![\xi_r(z)]\!]^{(k)}_j=0$ except for k=j. As the eigenmodes of the Alfvén operator are single sine contributions, the factor $\left(\overline{\phi}\Pi^{({\rm in})}\right)^{(k)}$ vanishes except for local Alfvén waves and global fast waves with the same longitudinal wave number. Thus, as before, the dispersion relation reduces to the familiar impedance-matching condition of the 1D case with the jump contribution in the right hand side:

 

$$\frac{\Xi^{({\rm ex},k)}_k}{\Pi^{({\rm ex},k)}_k}-\frac{\Xi^{({\r... ...t/\Pi^{({\rm in},k)}_k}{\vert \langle\phi L_{{\rm A}1}\phi \rangle\vert}\phi_k.$$ (7)

  
4 Linear expansion in ${\rm\alpha _n}$

Although the application to coronal loops involves strong longitudinal stratification and therefore requires values of $\alpha_n$ of order unity, we study in this section the linear expansion of the dispersion relation in $\alpha_n$ and thus impose $\alpha_n\ll 1$ since this limit can be tracked analytically. In a linear approximation in $\alpha_n$, inclusion of several $\alpha_n$ just leads to a sum of the contributions. The subscript on $\alpha_n$ and S_n is dropped for notational convenience.

We start by determining the first order effect on the eigenmodes of the local Alfvén operator. In terms of the sine series coefficient vectors $\psi^{(k,0)}_j$ the zeroth order eigenmodes $\psi^{(k,0)}$ and eigenvalues are clearly:

$$\begin{eqnarray*}\psi^{(k,0)}_j &=& \delta_{k,j} \\ \lambda_{k0} &=& \rho_0\left[\omega^2-v_{{\rm A}0}^2\left(\frac{k\pi}{L}\right)^2\right]\cdot \end{eqnarray*}$$

In order to find the first order effect on the eigenvalues, i.e. the roots $\lambda_k$ of the determinant of the matrix $L_{\rm A}-\lambda_kI$ with $L_{\rm A}$ the matrix operator as expressed in Eq. (4), we express the eigenvalue as $\lambda_k=\lambda_{k,0}+\alpha\lambda_{k,1}$. All off-diagonal terms of $(L_{\rm A}-\lambda_k I)$ are now proportional to $\alpha$ and the kth diagonal term is also proportional to $\alpha$. Therefore, the linear part (in $\alpha$) of the determinant is given by:

$$(\rho_0\omega^2\alpha S_{k,k}-\alpha\lambda_{k,1})\prod_{j\neq k}(\lambda_{j,0}-\lambda_{k,0})$$

which makes it possible to calculate $\lambda_{k1}=\rho_0\omega^2S_{k,k}$. Analogously, we express $\psi^{(k)}=\psi^{(k,0)}+\alpha\psi^{(k,1)}$, and find that the first order (in $\alpha$) part of $(L_{\rm A}-\lambda_k)\psi^{(k)}$ is:

$$\alpha(\lambda_{j,0}-\lambda_{k,0})\psi^{(k,1)}_j+\rho_0\omega^2\alpha (S_{k,j}-S_{k,k}\delta_{kj}).$$

Notice that for j=k both terms vanish. Thus we are only capable of determining $\psi^{(k,1)}_j$for $j\neq k$. This is, however, not a problem as a first order shift parallel to the zeroth order result is indeed meaningless. We obtain:

$$\psi^{(k,1)}_j=-\frac{\rho_0\omega^2 S_{k,j}}{\lambda_{j,0}-\lambda_{k,0}} \ \ \ \ \ \ \ \ \forall j\neq k.$$

We thus conclude that to first order in $\alpha$

$$\begin{eqnarray*}\lambda_k &=& \rho_0\left[\omega^2(1+\alpha S_{k,k})-v_{{\rm A}... ...}}{\lambda_{j,0}-\lambda_{k,0}} \ \ \ \ \ \ \ \ \forall j\neq k. \end{eqnarray*}$$

The eigenmodes are not normalized, and they need not be for dispersion relation (7) to remain valid.

We should, however, be aware that $\alpha$ may be different inside and outside. Therefore we replace $\alpha$ with $\varepsilon\sigma(r)$, where  $\varepsilon$ is the new small parameter independent of radial position, while $\sigma(r)$ is of order unity and $\sigma_{({\rm in})}/\sigma_{({\rm ex})}$ indicates the relative strength of the longitudinal stratification inside compared to outside the loop. We are looking for solutions to the dispersion relation in the neighborhood of the zeroth order solution of the kth quasi-mode so that we propose $\omega=\omega_0+\varepsilon\omega_1$. Note that the zeroth order terms (in $\varepsilon$) in the dispersion matrix only appear in $2\times 2$ blocks along the diagonal. It is not hard to verify that first order (in $\varepsilon$) contributions can then only arise from the product of the determinants of the $2\times 2$ block matrices. Since $\omega_0$ is the zeroth order frequency of the kth quasi-mode, the determinant of the kth $2\times 2$ block is proportional to $\varepsilon$ so that the linear part (in $\varepsilon$) of the dispersion relation is just:

$$\left\vert \begin{array}{cc} \Pi^{({\rm ex},k)}_k & -\Pi^{({\... ...{\rm in},j)}_j-[\![\xi_r]\!]^{(j)}_j \end{array}\right\vert_0.$$

Thus, the familiar impedance-matching expressions (7) for the dispersion relation remain true to first order in $\varepsilon$. In Appendix (A) we derive the first order (in $\varepsilon$) expansions of those impedances and of the jump contribution. This expansion leads to the following frequency shift:

 

$$\omega_1 -\frac{1}{2}\omega_0 S_{k,k} \frac{\mu_{({\rm in})}\sigma_{({\rm ... ...ht)_{({\rm res})}} {\mu_{({\rm in})}+\mu_{({\rm ex})}+\imath~\mu_{({\rm res})}}$$ =

$$-\frac{1}{2}\omega_0 S_{k,k} \overline{\sigma}$$ = (8)

with $\mu_{({\rm in})}$, $\mu_{({\rm ex})}$, $\mu_{({\rm res})}$ and $\nu_{({\rm res})}$ as defined in Appendix A. Apart from the term with $\nu_{({\rm res})}$ (which is due to the radial derivative of the stratification parameter at the resonance) the last factor in expression (8) is a weighted mean of the values of the stratification parameter at the internal and the external region and at the resonance with complex weighting coefficients $\mu_{({\rm in})}$, $\mu_{({\rm ex})}$ and $\mu_{({\rm res})}$. Thus when the stratification is constant we simply have:

$$\omega_1=-\frac{1}{2}\omega_0\sigma S_{k,k}$$

and thus to first order the effect on the oscillation frequency and on the damping is the same so that the observational parameter $\tau_{\rm damping}/{\rm Period}$ does not change to first order. This result can however be dramatically changed when the stratification parameter is not constant. As the weighting coefficients are complex, the mean value of the stratification parameter is complex as well and thus the effect of stratification may act differently on the oscillation frequency and on the damping rhythm. It is not hard to show that for small damping $\mu$ is approximately real and $\mu_{({\rm res})}\ll\mu_{({\rm in},{\rm ex})}$. Let us denote (the subscript ${\rm disc}$ refers to the fact that it represents the values in a discontinuous model):

$$\begin{eqnarray*}\mu_{({\rm disc})}&=&\mu_{({\rm in})}+\mu_{({\rm ex})}\\ \sigm... ...rm ex})}\sigma_{({\rm ex})}} {\mu_{({\rm in})}+\mu_{({\rm ex})}} \end{eqnarray*}$$

and find that to first order in $\mu_{({\rm res})}/\mu_{({\rm disc})}$:

$$\overline{\sigma}\approx\sigma_{({\rm disc})}+\imath~\frac{\m... ...u_{({\rm disc})}}(\sigma_{({\rm res})}-\sigma_{({\rm disc})}).$$

The dominant part in the shift of the oscillation frequency becomes clearly:

$$\omega_{{\rm r}1}\approx-\frac{1}{2}S_{k,k} \omega_{{\rm r}0}\sigma_{({\rm disc})}.$$

But for the dominant part of the imaginary part of the frequency we obtain:

$$\omega_{{\rm i}1}\approx-\frac{1}{2}S_{k,k} \left(\omega_{{\r... ...m disc})}}(\sigma_{({\rm res})}-\sigma_{({\rm disc})})\right).$$

Because $\omega_{{\rm i}0}\ll\omega_{{\rm r}0}$ the second term may be of the same order as the first term even though $\mu_{({\rm res})}\ll\mu_{({\rm disc})}$. When looking at kink modes that are trapped in the loop structure, $\mu_{({\rm disc})}$ and $\sigma_{({\rm disc})}$ are dominated by the internal contribution so that when the stratification in the boundary layer is more or less the same as in the internal region the second term vanishes and the influence of longitudinal stratification is again the same on the oscillation frequency and on the damping. However, when the stratification in the boundary layer is more like that in the external region the second term comes to play an important role and may dramatically change the observational parameter $\tau_{\rm damping}/{\rm Period}$.

Figure 1: Oscillation frequency $\omega _{\rm r}$ (panel a)), damping rate $\omega _{\rm i}$ (panel b)), and normalized damping rate $-\omega _{\rm i}/\omega _{\rm r}$ (panel c)) as a function of the stratification parameter $\alpha$. The straight line is the prediction from the linear analysis, the curved line is the non-linear numerical result. When $\alpha$ is varied, $\rho _0$ is kept constant.

Open with DEXTER

5 Numerical results for loop parameters

In order to implement the procedure outlined in Sects. 3 and 4 into a numerical code we need to cut down the infinite matrix problem to finite dimensions. It is not hard to see that the S_n matrix is dominant along its $\pm n$th off diagonal, thereby coupling the jth sine contribution mainly with the $j\pm n$th sine contributions. It can thus be estimated that the results obtained by reducing the system to finite dimensions (2l) are only accurate for the computation of eigenmodes (both of the Alfvén operator and of the full fast wave dispersion relation) that are dominant in the jth sine contribution provided that $l\gg j+n$, where n is the order of the highest sine contribution in the density.

As we have the application to coronal loop oscillations in mind we will study the effect of the stratification on the m=1 kink mode and look for solutions that have a dominant fundamental sine contribution (j=1). The stratification parameters that are appropriate for coronal loops can be estimated by assuming a semicircular coronal loop and an exponentially stratified atmosphere. For a coronal loop of typical length $200~{\rm Mm}$ (thus extending to a height $200/\pi {\rm Mm}\approx 64~{\rm Mm}$) and a density scale height of typically $50~{\rm Mm}$ we obtain:

$$\alpha_1=-0.78 \ \ \ \alpha_3=-0.076 \ \ \alpha_5=-0.017 \ \ \ldots$$

The effect of the longitudinal stratification will therefore be dominated by the fundamental component of the stratification. In our calculations we therefore only impose fundamental stratification and thus we set $\alpha_i=0$ for $i\neq 1$. Notice that $\alpha_1$ is negative because the density is higher at the tops than at the footpoints. Our computations show that inclusion of more than 7 sine contributions hardly changes the results which is in accordance to the estimate mentioned above.

In the following presentation of the results we use non-dimensional parameters. Lengths (length of the loop L and thickness of the boundary layer l) are scaled with respect to the inner loop radius. $\rho _0$ is scaled with respect to its values in the surrounding plasma. The magnetic field is constant and thus set to 1. Times (and thus frequencies) are therefore expressed in terms of Alfvén crossing times (using the loop radius and the external Alfvén speed). Figure 1 shows the frequencies and damping rates as a function of $\alpha$ with parameters $\rho_0^{(\rm in)}=2$, L=100 and l=0.01. In all three panels it can be seen clearly that the analytical linear prediction is indeed correct as the curves are tangent for $\alpha=0$. However, very strong deviations between the two curves are present for higher values of $\alpha$. For the loop value of $\alpha$ (-0.78) the frequency and damping rate is about 70% larger than for the unstratified case and about 30% larger than the linear prediction. The right hand panel shows minus the ratio of the imaginary and real part of the frequency. This quantity is linked with the ratio of the observational damping time and period. The curves are horizontal lines and thus this quantity is independent of the stratification even for non-linear values of the stratification parameter $\alpha$. It must be stressed that this could change drastically if the stratification parameter is allowed to vary with radial position. As we have no indication that this may be the case in loops, and have no estimate of  $\sigma_{\rm in}/\sigma_{\rm out}$, we did not carry out any calculations for such cases.

	Figure 2: Oscillation frequency $\omega _{\rm r}$ (panel  a)) and damping rate $\omega _{\rm i}$ (panel  b)) as a function of the stratification parameter $\alpha$. The straight line is the prediction from the linear analysis, the curved line is the non-linear numerical result. When $\alpha$ is varied, $\rho _0$ is adjusted so as to keep the top density constant.
Open with DEXTER

The above plots were obtained for varying $\alpha$ and the density at the footpoints kept constant. The density at for example the top of the loop is thus decreased for increasing values of $\alpha$. From an observational point of view, the density is probably easiest to determine at the loop tops and thus it would be better to compare results for different $\alpha$ while keeping the top density constant. Such a plot is shown in Fig. 2. Notice that now the frequency decreases with $\alpha$ contrary to Fig. 1. The linear prediction is modified as a term $\frac{\partial\rho_0}{\partial\varepsilon}$ needs to be included. After some algebra it turns out that in Eq. (8) S_kk needs to be replaced with S_kk-1.

More generally we could keep constant a sort of weighted mean of the density along the loop. The two previous plots then correspond to a weight function represented by a delta function at the footpoints and at the top respectively. It is clear that these two situations are the two extremes: any other weight function will lead to a result in between. As for constant footpoint density the frequency increases with $\alpha$ while it decreases when the top density is kept constant, we may want to search for a weight function so that increasing $\alpha$ while the weighted density is kept constant leaves the frequency invariant (at least linearly). It fact, the search for such a weight function may be inspired by trying to simulate the 2D results by using a mean density as the density in the 1D calculations. Naively, one could try to approximate the 2D result by using the mean density along the loop as the density in the 1D model. As we will show, this can indeed be done, but not simply by using the mean density, but by using a weighted mean density with an appropriate weight function attributing more importance to the loop tops than to the loop footpoints.

For any weight function the weighted density can be expressed as:

$$\rho_w=\rho_0 f(\alpha_n)$$

where f is a different function of $\alpha$ depending on the weight function (f=1 for constant footpoint density, $f=1+\sum \alpha_n$ for constant top density, $f=1+\frac{2}{\pi}\sum_{{\rm odd}n}\frac{\alpha_n}{n}$ for constant mean density). For the linear prediction we then find in general that in formula (8) S_nkk needs to be replaced with $S_{nkk}-\frac{\partialf}{\partial\alpha_n}$. A weight function that would lead to a linearly invariant frequency thus requires that $\frac{\partialf}{\partial\alpha_n}=S_{nkk}$. The definition of S_nkk is exactly the mean of $\sin{\frac{n\pi}{L}z}$ weighted with $\sin^2{\frac{k\pi}{L}z}$ so that it is clear that $\sin^2{\frac{k\pi}{L}z}$ is the weight function we are looking for. Figure 3 shows the frequencies and damping rates as a function of $\alpha$ while keeping that weighted mean of the density constant. The invariance of the frequency with $\alpha$seems to extend well outside the linear $\alpha$ domain. Even for $\alpha=-1$ the deviation is only 10% and for the loop value of $\alpha =-0.78$ the deviation is just a few percent at most.

	Figure 3: Oscillation frequency $\omega _{\rm r}$ (panel  a)) and damping rate $\omega _{\rm i}$ (panel  b)) as a function of the stratification parameter $\alpha$. The straight line is the prediction from the linear analysis, the curved line is the non-linear numerical result. When $\alpha$ is varied, $\rho _0$ is adjusted so as to keep the weighted mean density constant.
Open with DEXTER

We thus conclude that the frequencies and damping times of a stratified loop are more or less the same as those of an unstratified loop with the same weighted mean density (weighted with $\sin^2{\frac{k\pi}{L}z}$). The appearance of the weighting function $\sin^2{\frac{k\pi}{L}z}$ may not come as a surprise as it represents the wave energy density distribution along the loop of the fundamental mode.

The quasi-solutions are clearly influenced by the longitudinal stratification. Figures 4 and 5 show the pressure perturbation and the radial and azimuthal displacement (the longitudinal displacement vanishes as p=0) for a quasi-mode in an equilibrium with $\alpha =-0.5$ and for $\alpha =-0.78$. The jump of $\xi _r$ around the radial boundary position is due to the resonance. The jump in $\xi _\varphi$ is also present in discontinuous models and is just due to the different internal and external media.

It is intriguing that in the plots of the radial and azimuthal displacement no higher order longitudinal harmonics are visible. In the pressure perturbation, however, the higher order sines are visible through the more flattened longitudinal profile for $\alpha=0.5$ and even through the additional bumps for $\alpha =-0.78$.

As the total pressure perturbation is related to compression, these results suggest that the observations should show intensity oscillations somewhere halfway along the loop legs. However, as the oscillations are linear these intensity oscillations are very small and maybe not detectable. The clear detections of loop oscillations (including their period and damping times) are not based on the intensity oscillations due to compression but simply on the displacement of the more intense loop structure. However, in the displacement vector no visible signatures of higher order components are present. This may explain why no higher order harmonic components are visible in the observations although they may be present. In that respect it may be interesting to search for the associated intensity oscillations due to compression in the observations, as those should show signatures of higher order harmonics.

	Figure 4: Quasi-solution for $\alpha =-0.5$. Radial domain: [0.10] with loop radius at 1, longitudinal domain: [0,100] with L=100. Left panels: $p_{\rm T}$, middle panels: $\xi _r$, right panels: $\xi _\varphi$. Upper panels: real part, lower panels: imaginary part.
Open with DEXTER

	Figure 5: Quasi-solution for $\alpha =-0.78$. Radial domain: [0.10] with loop radius at 1, longitudinal domain: [0,100] with L=100. Left panels: $p_{\rm T}$, middle panels: $\xi _r$, right panels: $\xi _\varphi$. Upper panels: real part, lower panels: imaginary part.
Open with DEXTER

6 Summary and conclusions

The straight cylindrically symmetrical flux tube model of a coronal loop is extended by including longitudinal density stratification. In the internal and external regions where the equilibrium quantities are radially invariant, the longitudinal variation of the fast linear motions is decomposed into a linear combination of eigenmodes of the local Alfvén operator, which enables us to write the solution as a sum of separable terms. The matching condition of the pressure perturbation and the radial displacement then yields the dispersion relation in the form of the determinant of an infinite set of linear equations. The dispersion relation is shown to reduce to the familiar impedance-matching expressions when there is no longitudinal stratification.

In the thin boundary approximation the effect of the resonant absorption in the boundary layer can be summarized by means of a connection formula. This jump condition for the boundary layer is easily included in the dispersion relation. Again, the driven impedance-matching expressions are retrieved when there is no longitudinal density stratification.

A linear expansion in the longitudinal stratification parameter is performed. This revealed that the simple impedance-matching expressions remain valid to first order in the stratification. A linear expression for the frequency shift due to longitudinal stratification is obtained. It is found that both period and damping time are affected in the same way so that the ratio of the two (which is the most important observational parameter) remains unaffected. It is pointed out that this does not remain true when the stratification inside and outside the loop is different.

For coronal loops the longitudinal stratification is out of the linear domain and the dispersion relation is solved numerically by cutting down the infinite matrix problem to finite dimensions. When varying the stratification parameter it is important to state whether the footpoint density or the top density or some other weighted mean density is taken to be constant. It is found that the frequency is largely unaffected by the stratification even in the non-linear domain, when the mean density, weighted with $\sin^2(\frac{k\pi}{L}z)$ (k is the harmonic number, k=1 for fundamental loop oscillations), is held constant. Hence, the 2D results for frequencies and damping times may be approximated to a certain level of confidence by using the weighted mean density as the density in 1D computations.

The fast quasi-modes for large values of the stratification parameter show clear signatures of higher order harmonic components. However, these signatures are only visible in the pressure perturbation and not in the radial or azimuthal displacement. This might explain why these components are not visible in the loop motions although they might be present.

Acknowledgements

This work was initiated during a visit of M. Goossens to the Space Science Center of the university of New Hampshire. It is a pleasure for M. Goossens to acknowledge the hospitality of J.V. Hollweg during his visit and the financial support from the NASA Sun-Earth Connection Theory Program under grant NAG5-11797 to the University of New Hampshire.
I. Arregui acknowledges the support from the PLATON research training network: HPRN-CT2000-00152.

References

• Aschwanden, M. J., De Pontieu, B., Schrijver, C. J., & Title, A. M. 2002, Sol. Phys., 206, 99 [NASA ADS] [CrossRef]
• Aschwanden, M. J., Fletcher, L., Schrijver, C. J., & Alexander, D. 1999, ApJ, 520, 880 [NASA ADS] [CrossRef] (In the text)
• Aschwanden, M. J., Nightingale, R. W., Andries, J., Goossens, M., & Van Doorsselaere, T. 2003, ApJ, 598, 1375 [NASA ADS] [CrossRef] (In the text)
• Aschwanden, M. J. 2004, in Proc. of the NATO advanced research workshop on Turbulence, Waves and Instabilities in the Solar Plasma, Lillafured, Hungary, 16-20 september 2002, ed. R. Erdélyi, K. Petrovay, R. B., & M. J. Aschwanden (Dordrecht: The Netherlands: Kluwer Academic Publishers), 215 (In the text)
• Davila, J. M. 1991, in Mechanisms of Chromospheric and Coronal heating, ed. P. Ulmschneider, E. Priest, & R. Rosner (Berlin: Springer), 464 (In the text)
• De Pontieu, B., Martens, P. C. H., & Hudson, H. S. 2001, ApJ, 558, 859 [NASA ADS] [CrossRef] (In the text)
• Díaz, A. J., Oliver, R., & Ballester, J. L. 2002, ApJ, 580, 550 [NASA ADS] [CrossRef]
• Díaz, A. J., Oliver, R., Erdélyi, R., & Ballester, J. L. 2001, A&A, 379, 1083 [EDP Sciences] [NASA ADS] [CrossRef]
• Goossens, M., Hollweg, J. V., & Sakurai, T. 1992, Sol. Phys., 138, 233 [NASA ADS] [CrossRef] (In the text)
• Goossens, M., Andries, J., & Aschwanden, M. J. 2002, A&A, 394, L39 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Hollweg, J. V., & Yang, G. 1988, J. Geophys. Res., 93, 5423 [NASA ADS] (In the text)
• Ionson, J. A. 1978, ApJ, 226, 650 [NASA ADS] [CrossRef] (In the text)
• Nakariakov, V. M., Ofman, L., DeLuca, E. E., Roberts, B., & Davila, J. M. 1999, Science, 285, 862 [NASA ADS] [CrossRef]
• Ofman, L., & Aschwanden, M. J. 2002, ApJ, 576, L153 [NASA ADS] [CrossRef] (In the text)
• Ruderman, M. S., & Roberts, B. 2002, ApJ, 577, 475 [NASA ADS] [CrossRef] (In the text)
• Schrijver, C. J., Aschwanden, M. J., & Title, A. M. 2002, Sol. Phys., 206, 69 [NASA ADS] [CrossRef]
• Thompson, M. J., & Wright, A. N. 1993, J. Geophys. Res., 98, 15541 [NASA ADS] (In the text)
• Van Doorsselaere, T., Andries, J., Poedts, S., & Goossens, M. 2004a, ApJ, 606, 1223 [NASA ADS] [CrossRef] (In the text)
• Van Doorsselaere, T., Debosscher, A., Andries, J., & Poedts, S. 2004b, A&A, 424, 1065 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)

 

  
Online Material

  
Appendix A: Linear expansions of the impedances and the jump contribution

The dispersion relation can be summarized as:

$$Z^{({\rm ex})}-Z^{({\rm in})}=\imath I.$$

Where:

$$\begin{eqnarray*}Z^{({\rm ex})} &=& \frac{\Xi^{({\rm ex},k)}_k}{\Pi^{({\rm ex},k... ...n},k)}_k}{\vert \langle\phi L_{{\rm A}1}\phi\rangle\vert}\phi_k. \end{eqnarray*}$$

The dispersion relation can now be developed to first order in the stratification by requiring:

$$\frac{{\rm d}Z^{({\rm ex})}}{{\rm d}\varepsilon}-\frac{{\rm d... ...d}\varepsilon}=\imath \frac{{\rm d}I}{{\rm d}\varepsilon}\cdot$$ (A.1)

For the impedance terms there are two reasons why the impedance changes as a result of the stratification. Firstly, because the stratification parameter $\varepsilon$ is explicitly present in the expressions, secondly because the frequency of the quasi-mode changes as a result of the change in the stratification:

$$\frac{{\rm d}}{{\rm d}\varepsilon}=\frac{\partial}{\partial\v... ...mega}{\partial\varepsilon}\frac{\partial}{\partial\omega}\cdot$$

Moreover, it is interesting to see that the impedances are structured in the same way viz.:

$$Z=\frac{\kappa}{\lambda}G(\kappa R)$$

where G is a different function depending on the boundary condition. As $\kappa$ is a simple expression involving, except for $\lambda$, only constants, the impedances only change through the changes in $\lambda$:

$$\frac{{\rm d}Z}{{\rm d}\varepsilon}=\frac{{\rm d}\lambda}{{\rm d}\varepsilon}\frac{{\rm d}Z}{{\rm d}\lambda}\cdot$$

We straightforwardly obtain:

$$\frac{{\rm d}\lambda}{{\rm d}\varepsilon}=\rho_0\left(\omega^... ...{k,k}+2\omega\frac{\partial\omega}{\partial\varepsilon}\right)$$

and (using the shorthand $w=\kappa R$):

$$\begin{eqnarray*}\frac{{\rm d}Z}{{\rm d}\lambda} &=& -\frac{\kappa}{\lambda^2}G(... ...t]\\ &=& Z\frac{1}{2\lambda}\left(\frac{wG'(w)}{G(w)}-1\right). \end{eqnarray*}$$

For the change of the resonance contribution I due to the stratification some more care needs to be taken. The resonant position may shift due to the stratification. This shift may be calculated by requiring that the eigenvalue $\lambda$ vanishes at the resonant position:

$$\frac{{\rm d}\lambda}{{\rm d}\varepsilon}=\frac{\partial\lamb... ...tialr}{\partial\varepsilon}\frac{\partial\lambda}{\partialr}=0$$

so that:

$$\frac{\partialr}{\partial\varepsilon}=-\frac{\frac{\partial\l... ...2\omega\frac{\partial\omega}{\partial\varepsilon}\right) \cdot$$ (A.2)

Now let us investigate which factors in I may contribute to a linear change with the stratification parameter. First of all it must be noticed that when one considers the resonant mode that is already interacting in zeroth order the following factor does not yield first order contributions:

$$\begin{eqnarray*}\frac{\left(\overline{\phi}\Pi^{({\rm in})}\right)^{(k)}}{\Pi^{... ...({\rm in},k)}_j}{\Pi^{({\rm in},k)}_k}}_{\mbox{second order}}\\ \end{eqnarray*}$$

where we have used $\phi_k=1$ and for $j\neq k$, both $\Pi^{({\rm in},k)}_j\sim\varepsilon$ and $\phi_j\sim\varepsilon$. The same factor does not cause any first order contributions for additional resonances either. To see this, consider such an additional resonant mode dominant in the lth sine component ($l\neq k$). Then:

$$\begin{eqnarray*}\frac{\left(\overline{\phi}\Pi^{({\rm in})}\right)^{(k)}}{\Pi^{... ...k)}_k}\right)}_{\mbox{first order}}\phi_k}_{\mbox{second order}} \end{eqnarray*}$$

where we have now used $\phi_l=1$, $\Pi^{({\rm in},k)}_l\sim\varepsilon$ and $\phi_j\sim\varepsilon$ for $j \neq l$.

Therefore the linear changes in I can only result from the factor 1/r² and $1/\vert\langle\rangle\vert$(where we have used the shorthand notation $\langle\rangle=\langle\phi L_{{\rm A}1}\phi\rangle$):

 

$$\frac{{\rm d}I}{{\rm d}\varepsilon} \frac{\partialr}{\partial\varepsilon}\frac{\partialI}{\partialr}+... ...le\vert}{{\rm d}\varepsilon}\frac{\partialI}{\partial\vert \langle\rangle\vert}$$ =

$$-2\frac{I}{r}\frac{\partialr}{\partial\varepsilon}-\frac{I}{\vert \langle\rangle\vert}\frac{{\rm d}\vert \langle\rangle\vert}{{\rm d}\varepsilon}$$ =

$$-2\frac{I}{r}\frac{\partialr}{\partial\varepsilon}-\frac{I}{\langle\rangle}\frac{{\rm d}\langle\rangle}{{\rm d}\varepsilon}\cdot$$ = (A.3)

For determining ${\rm d}\langle\rangle/{\rm d}\varepsilon$ it must be noticed that $L_{{\rm A}1}$ is diagonal for $\varepsilon=0$ and that $\langle\phi{\rm d}\phi/{\rm d}\varepsilon\rangle=0$ so that:

$$\frac{{\rm d} \langle\rangle}{{\rm d}\varepsilon}=\frac{{\rm d}\left(L_{{\rm A}1}\right)^{(k)}_k}{{\rm d}\varepsilon}\cdot$$

Remember that we have to take into account the possibility of a shift in the resonant position and thus:

$$\frac{{\rm d}}{{\rm d}\varepsilon}=\frac{\partial}{\partial\v... ...\partialr}{\partial\varepsilon}\frac{\partial}{\partialr}\cdot$$

It is obtained straightforwardly that:

$$\begin{eqnarray*}\frac{{\rm d}\left(L_{{\rm A}1}\right)^{(k)}_k}{{\rm d}\varepsi... ...}{\partialr^2}\omega^2\frac{\partialr}{\partial\varepsilon}\cdot \end{eqnarray*}$$

Inserting this result in Eq. (A.3) and by using expression (A.2) for $\partial r/\partial \varepsilon$ we obtain:

$$\begin{eqnarray*}\frac{{\rm d}I}{{\rm d}\varepsilon}&=&\left(\omega^2\sigma S_{k... ...\omega^2}{\langle\rangle}\frac{\partial\sigma}{\partialr}S_{kk}I \end{eqnarray*}$$

Eq. (A.1) can now be solved with respect to $\frac{\partial\omega}{\partial\varepsilon}$ and yields result (8) with the following definitions:

$$\begin{eqnarray*}\mu_{({\rm in})} &=& -\rho_{0({\rm in})}\frac{{\rm d}Z^{({\rm i... ...ht]\\ \nu_{({\rm res})} &=& I\frac{\rho_0}{\langle\rangle}\cdot \end{eqnarray*}$$