Sensitivity of Coronal Loop Sausage Mode Frequencies and Decay Rates to Radial and Longitudinal Density Inhomogeneities:
A Spectral Approach
Abstract
Fast sausage modes in solar magnetic coronal loops are only fully contained in unrealistically short dense loops. Otherwise they are leaky, losing energy to their surrounds as outgoing waves. This causes any oscillation to decay exponentially in time. Simultaneous observations of both period and decay rate therefore reveal the eigenfrequency of the observed mode, and potentially insight into the tubes’ nonuniform internal structure. In this article, a global spectral description of the oscillations is presented that results in an implicit matrix eigenvalue equation where the eigenvalues are associated predominantly with the diagonal terms of the matrix. The offdiagonal terms vanish identically if the tube is uniform. A linearized perturbation approach, applied with respect to a uniform reference model, is developed that makes the eigenvalues explicit. The implicit eigenvalue problem is easily solved numerically though, and it is shown that knowledge of the real and imaginary parts of the eigenfrequency is sufficient to determine the width and density contrast of a boundary layer over which the tubes’ enhanced internal densities drop to ambient values. Linearized density kernels are developed that show sensitivity only to the extreme outside of the loops for radial fundamental modes, especially for small density enhancements, with no sensitivity to the core. Higher radial harmonics do show some internal sensitivity, but these will be more difficult to observe. Only kink modes are sensitive to the tube centres. Variation in internal and external Alfvén speed along the loop is shown to have little effect on the fundamental dimensionless eigenfrequency, though the associated eigenfunction becomes more compact at the loop apex as stratification increases, or may even displace from the apex.
1 Introduction
Coronal magnetoseismology seeks to infer the physical characteristics of solar coronal structures, particularly coronal loops, from observations of low frequency magnetohydrodynamic (MHD) waves and oscillations.
Currently unmeasurable physical parameters inside coronal loops such as magnetic field strength, fieldaligned flow magnitude, plasma temperature, and effective adiabatic index can be probed by various coronal seismology inversion schemes. Collective wave modes supported by magnetized tubes include kink and sausage modes. Both kink and sausage modes abound in the lower solar atmosphere.
Although much attention has been devoted to kink modes (, assuming an dependence on cylindrical angle ), there is a long history of interest in fast sausage modes () as well [1, 2, 3, 4, 5, 6]. Sausage modes are better suited to probing the radial structure of loop filaments, whereas kink modes are more suited to exploring their cores. They are widely thought to be implicated in quasiperiodic pulsations (QPPs) [7] in flare loops.
Slow sausage modes have independently received considerable attention [8, 9, 10]. These “anklebiter modes” are essentially fieldguided acoustic waves largely restricted to the lower reaches of loops. They are not present under the pressureless assumption adopted here, and will not be discussed further. Fast sausage waves on the other hand are fast magnetohydrodynamic waves partially (leaky) or totally (nonleaky) confined to overdense magnetic fluxtubes acting as wave guides.
Nonideal MHD mechanisms such as electron heat conduction, ion viscosity, and finite plasma are are not included at this stage as they are believed to be too weak to cause the temporal damping observed in QPP events [2, 3]. A detailed survey of the more extensive literature on coronal sausage waves and QPPs is given in [5].
Except for very short fat loops, (fast) sausage modes in coronal conditions are normally leaky. That is, they couple to radially outgoing oscillations in the surrounding plasma that extract energy from the loop oscillations [11, 12]. This causes them to exhibit complex eigenfrequencies, with negative imaginary parts. Observations of both period and decay rate therefore have the potential to usefully constrain loop properties.
Recently, [5] and [6] explored sausage oscillations in pressureless straight tube models using numerical solution of the governing ordinary differential equation (ODE) in radius , implementing a variety of density crosssections. Finite plasma is addressed in [13] but will not be considered here. These studies aimed to constrain the internal to external density contrast and the transverse Alfvén crossing time , where is the loop radius and the internal Alfvén speed, as well as obtain some information on the steepness of the radial density profile in the boundary layer between loop and external corona.
In this article, we address essentially the same model, though optionally with stratification along the loop allowed as well (the effects of stratification on kink waves have been wellstudied: [14, 15], etc.). However, rather than using a direct numerical ODE solver, we adopt a spectral decomposition that has several desirable properties. In particular, the resulting matrix eigenvalue problem is highly diagonally dominant, with the offdiagonal terms representing the radial inhomogeneity. This makes it easy to linearize the problem (in essentially boundary layer width), and thereby to identify density inversion kernels that explicitly show which modes are sensitive to which parts of the loop. This is carried out for the unstratified case, where a consequent perturbation analysis yields an explicit formula for frequency variance in terms of radial inhomogeneity parameters. Accuracy of the linearized perturbation formula is tested against full numerical solution.
Finally, the method is extended to allow for stratification along the loop, both internally and externally (independently). It is found that even mild stratification can significantly affect damping rate in particular, especially for an unstratified (hot) loop in a stratified atmosphere.
2 Mathematical Development
2.1 General Equations
Consider an axisymmetric Alfvén speed distribution in a pressureless straight “loop” of length and radius , surrounded by an external atmosphere with Alfvén speed dependent only on , (). In the pressureless plasma approximation, sound speed is neglected compared to Alfvén speed, freezing out slow modes. (Loops of nonuniform cross section are considered by [4], and nonzero plasma is addressed by [13].)
As is conventional, the loop is straightened for mathematical convenience [16], as depicted in Fig. 1, thereby neglecting curvature effects [17, 18, 19]. The coordinate is distance along the loop from one footpoint. The terms “loop” and “tube” will be used interchangeably throughout, despite the rectification.
If wave perturbations on this system are also assumed axisymmetric – the sausage modes – the linearized wave equation for the radial displacement assuming an time dependence is
(1) 
where represents the partial derivative with respect to , etc. (see for example [5], Equation (6)). The Alfvén speed is related to the magnetic field strength and plasma density via (in SI units)
(2) 
where henry is the permeability of free space [20].
For this sausage mode, there is no coupling to the Alfvén wave, unlike the case of kink waves, where coupling with the Alfvén wave results in resonant absorption and consequent kink mode decay. The wave under consideration therefore is a pure fast sausage wave.
The adopted boundary conditions are on , , and , and that and match continuously to an outgoing or evanescent external solution at .
Using separation of variables, the external solution may be written as
(3) 
for arbitrary , where is the Hankel function of the first kind and order 1 and is the eigenvalue/eigenfunction of the regular SturmLiouville equation
(4) 
with boundary conditions . The will be alternately even and odd about since is assumed even. By general SturmLiouville theory, orthogonality (adopting a particular normalization) is guaranteed provided , making the convenient expansion functions.
Internally, a Bessel/SturmLiouville expansion of the radial displacement is adopted,
(5) 
where the orthogonal expansion functions are
(6) 
and the are the internal radial wavenumbers belonging to a radially uniform reference tube. Here
(7) 
normalizes the kernel :
(8) 
where is the Kronecker delta. The orthogonality (8) is valid whether is real or complex.
If the external medium is not stratified in , i.e., , then , where are the longitudinal wavenumbers. The expansion is then Bessel/Fourier. In that case, the are given simply by
(9) 
. Otherwise, they are the eigenvalues of Equation (4), as above.
The internal radial wavenumbers are determined by the condition that each of the expansion functions individually couples to an evanescent or outgoing radiation solution in , i.e., that both and are continuous at . That is,
(10) 
If and then is pure imaginary (), the Hankel functions become modified Bessel functions of the second kind (the evanescent exterior), and the are real.
It is convenient to define . Multiplying Equation (1) by and integrating yields
(11) 
where
(12) 
When truncated ( and ), Equation (11) has nontrivial solutions if and only if the matrix
(13) 
is singular, where has components . Here , and Equation (13) defines . In practice, the eigenvalue problem is solved iteratively to make singular.
2.2 Uniform Density Tube
In the unstratified case, and Equation (9) applies, it is convenient to define , where is the fractional density increase inside the tube.
Equation (10) is also the appropriate matching equation when is uniform and there is an Alfvén speed discontinuity at . In that case though, the are not just expansion functions; they are independent modes with their own eigenfrequencies
(14) 
The condition for a sausage mode to be trapped (nonleaky) is that be imaginary, i.e., , hence , or . With , trapped sausage modes are only to be expected for very large density enhancements , where as is the zero of the Bessel function. For typical long thin loops this is highly implausible for the longitudinal fundamental or low harmonics, even for , so only leaky sausage modes need be considered. (An example of a short dense uniform loop that supports a nonleaky sausage mode is , , where for , for which as anticipated.)
The eigenfrequencies set out in Equation (14) provide useful reference states for perturbative solution of the nonuniform case. Table 1 lists the first few uniformdensity eigenfrequencies (all leaky) for , 1, and 4 for fat () and thin () loops. The corresponds to the longitudinal wavenumber , and the refers to the radial order. These are the socalled “Trig Modes” analysed by [11, 12]. Without loss of generality the loop length and external Alfvén speed have both been scaled to unity.
Thick Loop  Thin Loop  

Clearly, dependence of frequency on is weak, and on is strong, indicating the relative dependence of frequency on the longitudinal and radial variations in . Sausage modes, by their nature, are sustained by radial total pressure variations, not longitudinal tension (which is the primary mechanism for kink waves).
2.3 Nondimensionalization
For numerical purposes it is convenient to nondimensionalize throughout. Adopting as the fiducial length (the loop length) and velocity (the external Alfvén speed at the loop apex), we also have the time scale . Hence, any dimensionless frequency calculated below should be understood as corresponding to the dimensional frequency . The radii are of course all relative to .
3 Unstratified Loop with Boundary Layer
In this section, it is assumed that both the loop and background are unstratified in . Then , where is the fractional density increase inside the tube with respect to the exterior.
3.1 Sensitivity of Eigenfrequencies to Density Inhomogeneities
Assume that the spectral expansion is truncated to and . Therefore, Equation (11) represents equations in unknowns, . Let be a single vector rearrangement: , , …, , , …, where . Then the truncated Equation (11) may be written in matrix form
(15) 
where and . Here the diagonal matrix , and is the symmetric matrix with element , where , , , and , i.e.,
(16) 
By Gerschgorin’s Circle Theorem [21], all eigenvalues lie within the union of the Gerschgorin disks in the complex plane centred at
(17) 
with radius
(18) 
Furthermore, any set of contiguous disks that is disjoint from the other disks contains exactly eigenvalues. The disk centres are in exact accord with Equation (14), the uniformtube eigenfrequencies.
Several points should be made about the eigenvalue equation (15).

The offdiagonal terms in vanish for uniform , while the diagonal terms are just .

Increasing the truncation limits and does not alter the Gerschgorin centre , but does (slightly) increase the radius . Typically, for the fundamental mode , introducing a boundary layer (with as in Section 3.3 below) and increasing and from 3 to 4 results in a change in the numerically determined frequency only in the sixth significant figure. For the , first radial overtone, the change is in the fifth significant figure. In either case, is adequate in practice.
With that understanding in place, a linearized perturbation approach offers useful insights.
3.2 Perturbative Approach
The Gerschgorin disk centres , as defined by Equation (17), provide estimates of the true eigenvalues, becoming exact as the offdiagonal terms vanish with diminishing inhomogeneity. However, since actually depends on , the estimate is implicit rather than explicit. Assuming is small though, where the tilde refers to the uniform reference state, a linear expansion may be made. Note that this does not require to be pointwisesmall, only that for all , , , . In particular, this is valid for a thin enough boundary layer in which drops to zero on the edge, despite the total drop being large.
Recalling that is a function of , and expanding , , and , Equation (15) becomes to first order in small quantities. If say, then the row (i.e., the row) of is identically zero. Hence, since the associated eigenvector is , where the “1” is the element, it follows that
(19) 
where is the diagonal element. In other words, the perturbation to any eigenvalue is just the perturbation to the corresponding diagonal element of , to leading order. That means that the perturbed eigenvalue is simply the centre of the Gerschgorin disk, to first order.
With that result, it follows from Equation (17) that
(20) 
to first order. The departure of from is entirely due to the direct change in and not the change in , since the normalization (8) remains valid for any . Hence, may be calculated using the radial wavenumbers. The additional term proportional to results from the additional dependence of the radial wavenumbers on via the matching condition (10).
Now,
(21) 
defining the dimensionless quantity
(22) 
which is obtained by differentiating Equation (10). Ultimately, the fractional perturbation in the eigenfrequency of the mode is
(23) 
to first order, where
(24) 
is the kernelaveraged variation in density from the uniform state. Note that has not been assumed small, but all the terms have.
Typically, the corrections are not negligible. For example, with and , for ; for ; and for .
An arbitrary density perturbation cannot uniquely be reconstructed from even a complete and perfect knowledge of all the sausage mode frequencies. Most obviously, is insensitive to any part of that is odd about , because is even about that point for all . Consequently, the odd part of depends only, and very weakly, on the offdiagonal terms that enter at higher order.
The density sensitivity kernels for the thick loop are displayed in Fig. 2 for and 2 (the radial fundamental and first overtone), (longitudinal fundamental), and the three values of listed in Table 1.
The radial fundamental clearly gives good sensitivity to density at the loop edge. The first harmonic is more sensitive to the interior, especially for larger . Because of the factor in , there is no sensitivity at all to the tube centre. This is to be expected, since is a node for sausage modes. Indeed, it is a node for all higher azimuthal order modes except for kink modes, so only kink mode frequencies will be sensitive to the centre. Similarly, the loop apex is a node for (not shown), or any even . The modes are most sensitive at and .
Corresponding plots of the kernels for the thin loop are very similar (not shown).
3.3 Results: Density Independent of
Consider a twoparameter density model independent of ,
(25) 
that has a central plateau peaking at at the centre, and a boundary layer falling continuously to zero at the edge . The boundary layer becomes thinner as diminishes (see Fig. 3), with its thickness at halfheight being . The uniform reference model is recovered for and .
3.3.1 Linear Perturbation Solutions
3.3.2 Nonlinear Solution
Although more expensive, the eigenfrequencies for a specific model may be determined numerically without recourse to linearization. This is achieved using a standard rootfinder to iteratively make singular (where ; recall that also depends on in a complicated manner) given and . A starting guess is supplied by linear theory, allowing the selection of and .
Figure 5 shows contours of the real and imaginary parts of and for the thick tube over a wide range of and . It was calculated using spectral resolution. Note that there is now no need for a reference model . Figure 6 similarly depicts for the thin tube .
The fact that the contours of the real and imaginary parts of in Figures 5 and 6 for the most part cross at angles that are not too fine suggests that observations of period and decay rate are sufficient to determine and with some accuracy, assuming of course that the model (25) is valid. Simultaneous inversions using two modes, say the fundamental and first radial overtone , would lend validity to the model if the resulting were consistent.
4 Formulation with a Stratified Loop and Atmosphere
We return to the gravitationally stratified atmosphere. For concreteness, assume that the external density diminishes exponentially with height in a semicircular loop, so that
(26) 
where is the external Alfvén speed at the apex . It is assumed that , making Alfvén speed continuous at the tube boundary. No restriction is made on at the loop ends, and . Without loss of generality, time may be scaled by setting throughout, so the external atmosphere is characterized by one parameter only, .
Although more complex internal densities may be accommodated by the expansion method, for simplicity it will be assumed here that takes the form (with )
(27) 
which again represents a continuous Alfvén speed at . The density scale heights inside and outside the loop match if , in which case
(28) 
Equation (27) corresponds to a loop centre Alfvén speed specified by , representing a hot loop (relative to the external atmosphere) if , and a cool loop if .
With Equation (26) in place, the eigenvalue problem (4) with zero endpoint conditions is solved numerically for the and . For simplicity, attention will be restricted to even modes in by selecting only expansion functions that are symmetric about . With a symmetric Alfvén speed profile, as chosen here, the odd and even modes are decoupled in any case.
4.1 Stratified Atmosphere Results
For the most part, the fundamental mode is adequately represented with spectral resolution , (used throughout unless otherwise noted). That is, there is very little energy in the and last retained expansion functions, as judged by the squared magnitudes of the expansion coefficients . Indeed, by Parseval’s theorem, , so the are just the fractional energies in each expansion mode. The are calculated by solving once the frequency has been iteratively adjusted to make singular (see Eqns. (11) and (13)). The longitudinal expansion modes (solutions of Eqn. (4)) are ordered by decreasing , which assures that the longitudinal external eigenfunctions become increasingly oscillatory as increases, so that is identifiably the fundamental.
Since stratification is relevant chiefly to long high loops, attention will be restricted to the “thin loop” case, (radius divided by length). Figure 7 displays the eigenfrequencies for the case , , and varying stratification parameters and . Clearly, these dimensionless frequencies dependend only weakly on stratification.
However, dimensional frequency scales as . So, considering a semicircular loop of length , fixed base Alfvén speed , and Alfvén (i.e., density) dimensional scale height , the dimensional frequency scales with the dimensionless frequency according to . Hence, for a highly stratified loop, , the dimensional frequency (and decay rate) increase exponentially with loop length.
Although the dimensionless eigenfrequencies are almost insensitive to stratification, the displacement eigenfunctions evolve considerably as increases (Fig. 8), becoming progressively more compactly situated at the loop apex . This corresponds to the fundamental sausage mode being trapped at the apex by stratification rather than the boundary conditions at and 1. The energies in the constituent expansion functions for each of these eigenfunctions are displayed in Fig. 9, showing that dominates throughout, as one might expect from Fig. 8.
On the other hand, when the loop itself is unstratified () but the external region is stratified (quantified by ), some higher overtones become prominent (Fig. 10), ultimately shifting the peak displacement away from the apex, though again it is concentrated by increasing . Figure 11 confirms that the immediate longitudinal overtones are indeed quite prominent. Note that in this case, it was necessary to use enhanced resolution or 12 for small .
Figures 9 and 11 also demonstrate that there is very little contribution of to the fundamental. This is hardly surprising given the wide disparity in different eigenfrequencies of the uniform tube exhibited in Table 1.
Finally, very similar behaviour in the eigenfunctions is observed when is held fixed (the external atmosphere is unstratified) and increased (the loop is stratified internally; Fig. 12). In this case the expansion functions are just the , so the compacting of the lobe with increasing is produced by Fourier superposition, and not a change in , as is the predominant mechanism in the other cases.
There are two surprising results from this section.

Figure 7 reveals a remarkable insensitivity of the fundamental complex eigenfrequency to stratification, if the Alfvén speed at the apex is kept fixed. This suggests that is strongly determined by that local value of .

The extent to which the fundamental eigenfunctions are modified by even weak stratification is remarkable. This is a case where the eigenfunction is much more sensitive than the eigenfrequency. This may be understood qualitatively by recourse to eigenvector sensitivity theory in linear algebra, where it is known that eigenvectors depend more sensitively on parameters when the eigenvalues are closely spaced [22]. Reference to Table 1 confirms that indeed the eigenvalues of the uniform tube are indeed very closely packed in . This makes nearby harmonics more accessible.
5 Conclusions
One conceptual and practical advantage of the spectral approach to the unstratified loop presented here is that it places the dominant contribution to the eigenvalues on a matrix diagonal, making them almost explicit, with the offdiagonal corrections being moments of the density inhomogeneity, vanishing if the tube is uniform. This makes the spectral formulation the most natural for exploring sensitivity of sausage mode eigenfrequencies to density inhomogeneities.
A particular benefit is that the method allows the identification of the linear density kernel , that directly shows the sensitivity of the various modes to density variations in different parts of the loop. Figure 2 shows clearly that the fundamental is strongly dependent upon the extreme loop edge, and that would be needed to more accurately probe the interior. Observing this mode presents a challenge, but would be worthwhile.
Related to this, it is important to understand that the results presented in Section 3.3 in terms of and do not actually probe the loop centre, despite being ostensibly the density there. That is an artefact of the model (25). In fact, a density anomaly near is all but invisible to the fundamental mode, so and should be more properly thought of as parameterizing the boundary layer when used with . Higher radial order modes can ameliorate this deficiency somewhat, but the kink mode is a more natural discriminator of loop centre density.
Finally, in Section 4 the general SturmLiouville expansion formalism is applied to allow stratification along the loop, both internally and externally. Numerical results carried out for several cases indicate that stratification leaves dimensionless eigenfrequencies almost unchanged, as measured in a dimensionless system based on the apex external Alfvén speed. Of course, stratification will alter this apex Alfvén speed and hence the actual dimensional frequency and decay rate. On the other hand, stratification has significant consequences for the fundamental eigenfunction, by either making it more compact around the loop apex as stratification increases, or possibly adding significant contributions of higher longitudinal order (). The main lobe(s) of the fundamental may even be displaced from the apex.
References
 [1] Aschwanden M J, Nakariakov V M and Melnikov V F 2004 ApJ 600 458–463 (Preprint astroph/0309493)
 [2] Kopylova Y G, Melnikov A V, Stepanov A V, Tsap Y T and Goldvarg T B 2007 Astronomy Letters 33 706–713
 [3] Inglis A R, van Doorsselaere T, Brady C S and Nakariakov V M 2009 A&A 503 569–575 (Preprint 1303.6301)
 [4] Pascoe D J, Nakariakov V M, Arber T D and Murawski K 2009 A&A 494 1119–1125
 [5] Chen S X, Li B, Xiong M, Yu H and Guo M Z 2015 ApJ 812 22 (Preprint 1509.01442)
 [6] Guo M Z, Chen S X, Li B, Xia L D and Yu H 2016 Sol. Phys. 291 877–896 (Preprint 1512.03692)
 [7] Nakariakov V M and Melnikov V F 2009 Space Sci. Rev. 149 119–151
 [8] De Moortel I, Ireland J, Walsh R W and Hood A W 2002 Sol. Phys. 209 61–88
 [9] Nakariakov V M and Verwichte E 2005 Living Rev. Solar Phys. 2 URL http://www.livingreviews.org/lrsp20053
 [10] De Moortel I 2009 Space Sci. Rev. 149 65–81
 [11] Cally P S 1985 Australian Journal of Physics 38 825–837
 [12] Cally P S 1986 Sol. Phys. 103 277–298
 [13] Chen S X, Li B, Xiong M, Yu H and Guo M Z 2016 ApJ 833 114 (Preprint 1610.03254)
 [14] Andries J, Goossens M, Hollweg J V, Arregui I and Van Doorsselaere T 2005 A&A 430 1109–1118
 [15] Erdélyi R and Verth G 2007 A&A 462 743–751
 [16] Ruderman M S and Roberts B 2002 ApJ 577 475–486
 [17] Smith J M, Roberts B and Oliver R 1997 A&A 317 752–760
 [18] Van Doorsselaere T, Debosscher A, Andries J and Poedts S 2004 A&A 424 1065–1074
 [19] Verwichte E, Foullon C and Nakariakov V M 2006 A&A 446 1139–1149
 [20] Priest E R 1982 Solar Magnetohydrodynamics (Dordrecht: D. Reidel)
 [21] Gerschgorin S 1931 Izv. Akad. Nauk. USSR Otd. Fiz.Mat. Nauk 6 749–754
 [22] Golub G H and van Loan C F 1996 Matrix computations 3rd ed (The Johns Hopkins University Press)