Leaky modes in coronal magnetic flux tubes revisited

It is shown that the extensive literature on leaky modes in coronal magnetic flux tubes should be revised since those modes have no physical meaning.


Introduction
To set the scene for the discussion of leaky modes in coronal flux tubes, this section starts with a review of the existing literature, not striving for completeness but just stressing the main issues ( § 1.1).Next, the accepted standard model is reconstructed for the plane slab version of it in order to clarify these issues and even correct some published mistakes ( § 1.2).In the last part of this section, three defects of the standard leaky mode model are pointed out ( § 1.3).The rest of the paper is devoted to the necessary corrections of the model, with numerical examples, and the construction of the modified solutions for the modes in these magnetic configurations ( § § 2-4).They exhibit no damping and associated leakage of energy, so that the final conclusion is the rejection of the theories on leaky modes in coronal magnetic flux tubes ( § 5).

Discussion of the literature on leaky modes
The starting point of this paper is the failed attempt, by adding some corrections of § 10.5 to the Errata of our textbook (Goedbloed, Keppens & Poedts 2019), to construct a physically acceptable presentation of what has been termed 'leaky modes in coronal magnetic flux tubes' in the literature.We show that such modes are mathematical constructs without physical relevance to the dynamics of plasmas in those environments.Consequently, the mentioned § 10.5 of our textbook should be deleted.
In the publications of Meerson, Sasorov & Stepanov (1978), Wilson (1981) and Spruit (1982) an attempt was initiated to classify the enormous variety of waves in coronal magnetic flux tubes for different choices of their geometry and wavenumbers of the perturbations.Amongst the different possibilities also waves were assumed for which 'the wave amplitudes decrease along the tube away from the disturbance as the wave energy leaks away from the tube in the transverse direction' (Wilson 1981) or waves that 'can be damped by radiating an MHD or acoustic wave into the surroundings of the tube' (Spruit 1982).These assumptions were articulated by Cally (1986) who investigated the different classes of waves for rather general cylindrical geometries and wavenumbers and classified the different waves with real frequencies.He also demonstrated that, in fact, 'many other (leaky) modes are found which excite waves in the external medium and thereby lose energy to the surroundings'.The latter require complex values of the frequency ω and the assumption of a sudden transition of standing waves in the main plasma column to exclusively outgoing waves in the external region.The external region was modelled as a homogeneous plasma so that the solutions could be represented as Hankel functions with known asymptotic behaviour at large distances approaching infinity.One of the Hankel functions represents outgoing waves and the other, which was to be rejected, ingoing waves.The result is a dispersion equation with complex eigenfrequencies.
Quite a number of follow-up publications on the subject of leaky waves in coronal flux tubes appeared.For example, investigations on how the damped waves could be described, in analogy with the quasi-modes of the Alfvén subspectrum, by impedance matching of the resonant leaky modes to an external driver (Keppens, Bogdan & Goossens 1994;Keppens 1995;Stenuit, Keppens & Goossens 1998;Stenuit et al. 1999).A similar analysis of the complex dispersion function with the Laplace transform by Ruderman & Roberts (2006b), leading to the elimination of leaky kink modes when the density of the inner plasma is higher than that of the outer plasma (in conflict with an earlier paper by Cally (2003)), was rejected by Cally (2006) using the same method but, in turn, claimed again to be in error by Ruderman & Roberts (2006a).Further analysis was performed by Terradas, Oliver & Ballester (2005) for plane slab geometry and Terradas, Andries & Goossens (2007) for cylindrical geometry, claiming that the results of Cally (2003) were correct after all, supported by an attempt to relate the damped leaky modes to the real continuous spectrum by Andries & Goossens (2007).Finally, § § 5.6 and 6.6 of the textbook on magnetohydrodynamic (MHD) waves by Roberts (2019) contain an extensive analysis of the leaky modes, summarising most of the results obtained so far.
Most of the mentioned publications concern cylindrical geometry and general wavenumbers, but they all centre on the sudden transition to 'outgoing waves only' at the position of the jump of the Alfvén velocity.This leads to a dispersion equation when complex values of ω are admitted.This analysis can be greatly simplified for the analogous plane slab model, presented in § 1.2.That section is constructed according to the same (wrong) assumptions of the mentioned leaky mode analyses.Subsequently, the three major defects of this model, described in § 1.3, are repaired in § § 2-4, leading to the conclusion that all leaky modes of the kind discussed in the citations of this section should be rejected on physical grounds.That conclusion carries over to the cylindrical case.

The standard leaky mode model
The literature discussed in § 1.1 concerns models of coronal flux tubes with quite some details on contributions of the geometry (frequently cylindrical) and wavenumbers (including tangential ones in addition to the radial number) in order to facilitate comparison with observational data.Since there is a grave error in all these papers, we will do quite the opposite and simplify the model as much as possible without losing the essential problem that we address, namely, do these models allow for damping due to leakage of energy at large distances or not?
We consider the simple model of a plasma slab with a dense inner plasma and a tenuous outer plasma, mimicking a cylindrical flux tube of radius r = a surrounded by a tenuous plasma.The plasmas are assumed pressureless (vanishing sound speed, c ≡ γ p/ρ = 0).The magnetic field B is considered constant throughout (otherwise there would be no equilibrium) and the density jumps from ρ to ρ. (Actually, it makes no difference for the theory if the 'dense' and 'tenuous' properties of the two plasmas are interchanged, the only essential feature is a density jump.)Consequently, the Alfvén speed jumps from b ≡ B/ √ ρ to b ≡ B/ ρ at the boundaries of the inner plasma.The hat will be used consistently for variables in the outer region.
We immediately get rid of the three constants a, ρ and b, effectively equating them (and also the constant μ 0 ) to 1.They just provide units of length, mass and time corresponding to the scale independence of the ideal MHD equations; see § 4.1.2 of Goedbloed et al. (2019).The dimensionless 'radial' coordinate will be indicated by x, the jumps occur at |x| = 1, whereas the inner region is indicated by −1 < x < 1 and the outer regions are indicated by x < −1 and x > 1. Apart from the wavenumbers, the only quantity that matters is the Alfvén velocity in the inner region, i.e. b, since b ≡ 1 by normalisation.Consequently, the equilibrium is a one-parameter family described by the parameter δ ≡ b/ b, where δ < 1 for a higher density and a smaller Alfvén velocity in the inner than in the outer region and δ > 1 for the reverse.
The spectrum of MHD modes for this model is obtained by solving the ordinary differential equation (ODE) for the normal component ξ ≡ n • ξ (with n ≡ e x ) of the plasma displacement vector ξ , obtained by Goedbloed (1971a), where we neglect the gravitational terms whereas the slow magnetosonic modes have disappeared by the assumption c = 0: The magnetic field is taken in the z-direction, so that the tangential wavenumbers are k ⊥ ≡ k y and k ≡ k z .The normal modes for this model are then expressed by from which the other components of the vector ξ follow by algebraic expressions in terms of ξ and its first derivative ξ .The wavenumber in the normal direction is indicated by the different symbol q to remind us that, in general, when the plasma is inhomogeneous in x, the solutions cannot be represented by a wavenumber in that direction.We can pull the equilibrium quantities through the differential operator since they are all constant (except of course at |x| = 1).The Alfvén solutions, ω 2 = k 2 b 2 , can then be extracted so that (1.1) simplifies to a description for the fast modes only: (1.3) A similar equation holds ξ in the external region, with Alfvén velocity b.Excepting the boundaries at |x| = 1, these equations have constant coefficients in the two separate regions so that the internal and external solutions may be constructed immediately: Of course, the two dispersion equations refer to the same value of ω so that the wavenumbers q and q are related, The contribution k 2 0 of the tangential wavenumbers may give rise to differences of oscillatory or evanescent behaviour in the two regions.
It remains to apply the proper boundary conditions (BCs) joining the inner and outer solutions.Before we do that, it is expedient to exploit the symmetry about x = 0.This implies that there are two independent sets of solutions, which may be determined separately and, most importantly, only for x ≥ 0 since the dependence for x ≤ 0 follows from symmetry.The antisymmetric (sine-shaped) solutions are subject to the BC the symmetric (cosine-shaped) solutions are subject to the BC In the literature, these modes are called 'sausage modes' and 'kink modes', respectively, because of the analogy with modes in cylindrical geometry.(For that geometry, the solutions close to the origin are, rather confusingly, described by the Bessel function J 1 for the m = 0 sausage modes and by J 0 for the m = 1 kink modes.)At this point, even before applying the BCs at the separating surfaces, we need to impose the crucial conditions that restrict the waves in the outer regions to be outgoing only.It now becomes imperative to admit complex values of q and ω, q ≡ q1 + iq 2 , ω ≡ σ + iν, (1.9) where the subscripts 1 and 2 indicate the real and imaginary component of the complex quantity involved.The normal propagation part of the external waves may then be expressed as ξ(x, t) = αe −(q 2 x−νt) e i(q 1 x−σ t) + βe (q 2 x+νt) e −i(q 1 x+σ t) .
(1.10) Since, by symmetry, we need to consider only x ≥ 0, and time progresses, it is evident that the following constraints should be applied: β = 0 (for outgoing waves), α = 0 (for incoming waves). (1.11) For 'leaky modes', the first restriction applies.The BCs at the interface between the two plasmas express continuity of the normal displacement ξ and of the perturbation Π of the total pressure.The expression for the latter variable may be obtained from the explicit expression given in Eq. (7.88) of Goedbloed et al. (2019), which here reduces to (1.12) The BCs for the outgoing 'leaky waves' now become ξ(1) = ξ(1) ⇒ α(e iq ∓ e −iq ) = αe iq , (1.13) (1.14) where the upper sign applies to the antisymmetric modes and the lower sign applies to the symmetric modes.This yields the following dispersion equation for these waves: 1 .15)and recall that q and q are related through (1.6).
The dispersion equation is a transcendental equation for the determination of the normal wavenumber q, which could be solved by numerical means.However, it serves no purpose to strive for generality here because the essential problem is well exposed by ignoring the dependence on the tangential wavenumbers by setting k 0 = 0, so that (1.16) (One of the square root solutions of (1.6) also yields q/q = −δ, but this case can be covered by flipping the sign of q and assigning a different complex value to α.)The resulting dispersion equation can be solved explicitly to give -antisymmetric modes: (1.17) -symmetric modes: (1.18) The two solutions may be combined into one expression that alternatingly applies to the antisymmetric (even n) and to the symmetric (odd n) solutions: The expressions are given for δ < 1.For δ > 1, an extra phase shift of 1 2 π should be added to q 1 ≡ Re(q), so that the expressions for odd and even modes are interchanged but their damping rates are given by the same logarithmic expression.For inward propagating waves ( α = 0), the sign of the imaginary part of these expressions is flipped so that those modes are unstable with logarithmic growth rates.
This appears to establish that normal fast magnetosonic waves are standing waves in the inner region, characterised by the Alfvén velocity b, but that they propagate outwardly with the Alfvén velocity b.Moreover, all these waves are damped with the same logarithmic factor depending only on the jump δ of those velocities.However, there is a problem.Quite counterintuitively, the damping rates explode logarithmically in the limit δ → 1, both from below and from above, i.e. for any tiny jump.What physical mechanism could possibly cause such a miraculous effect?In addition, it is evident that (1.15) has no solution in that limit.We will have to answer these questions.

Defects of the leaky mode model
In order to address these problems, we note that the standard leaky mode model presented thus far has three defects.
(i) Because of the jump of the Alfvén velocity, the differential equation (1.1), or (1.3) for k ⊥ = 0, has discontinuous coefficients and, hence, discontinuous solutions with a jump q → q of the radial wavenumber and an associated jump of the derivative ξ , or ξ for k ⊥ = 0, at the sharp boundary between the two plasmas.In general, for wave propagation problems in discontinuous media, one has to show that a limiting procedure by smoothing the discontinuity can be designed that leads to smooth solutions.This is usually no problem but in this particular case it is.As we noted previously, the limiting procedure of ever smaller jumps remains singular to the very limit of vanishing jump when the eigenvalue blows up.This is related to the fact that the trivial solution does not even exist because it is impossible that solutions suddenly change from standing waves to propagating waves in a homogeneous medium.The content of our next item is that this remains true for finite jumps of the Alfvén velocity.May be smoothing can cure the problem?(ii) The procedure of imposing the condition of 'outgoing waves only' for x > 1 is not dictated by physical conditions since the jump in the Alfvén velocity does not demand that.It is purely an a priori assumption on what one would like the system as a whole to do, not an a posteriori conclusion on what it actually does.The reasonable assumption on the solution at x = 1 is simply to demand continuity there, as required by the description in terms of a differential equation.(iii) The leaky mode solutions of the simple model become exponentially large at infinity.This corresponds to exponentially large energy outflows at that position.Since this is physically impossible, when solving the differential equation properly, those solutions should be excluded.
We will repair these three defects in the following three sections and, finally, construct the physically acceptable solutions for the elementary model of the flux tube.

Smoothing the jump of the Alfvén velocity profile
We have seen that the solutions found in § 1.2 contain singularities that are related to the discontinuity of the Alfvén velocity.The obvious remedy is to smooth this function.This is easily done, but the drawback is the loss of explicit algebra in favour of numerical analysis.Fortunately, this is readily available in the numerical program ROC, used before to describe the new super-Alfvénic rotational instabilities (SARIs) in accretion discs about black holes by Goedbloed & Keppens (2022).
The mentioned program makes use of the method of the spectral web, extensively described in previous papers on the subject (Goedbloed 2018a,b).The cylindrical geometry used in these papers is here simplified to the plane geometry of the present elementary flux tube model.The relevant differential equation derived by Goedbloed (1971a) will be exploited, eliminating gravity and the contributions of the slow and Alfvén waves, whereas k 0 = 0 is assumed again.The only modification of the model of § 1.2 is the replacement of the constant Alfvén velocity in the inner region by so that the jump is smoothed in the outer part of the inner region.
The spectral web method consists, basically, of solving the differential equations for ξ(x) and ξ(x) in the two regions for arbitrary complex values of ω, join them by applying the BC (1.13), but replacing the BC (1.14) for Π by the computation of the complementary energy over the boundary surface: Here, L y and L z are arbitrary lengths in the y-and z-directions of the tangential plane, having no other function than to provide the proper dimensions to W com .
The procedure replacing the analysis of § 1.2 is as follows.A relevant piece (e.g. a square) of the complex ω-plane is selected and the real and imaginary parts of W com are computed by substituting the solution pairs {ξ, Π} and { ξ, Π } for all values of ω scanned in that piece.The function ξ(x) is obtained numerically by solving the mentioned ODE, or rather the equivalent set of first-order ODEs for ξ and Π , on the inner interval 0 ≤ x ≤ 1 by applying either one of the BCs (1.7) or (1.8) at x = 0.For the outer interval x ≥ 1, the function ξ is the explicitly given solution of (1.5), with β = 0 , whereas the associated expression for Π follows from the equivalent of (1.12).Satisfaction of the BC (1.13) is obtained by simply adapting the amplitude ξ(1) to ξ(1) or vice versa.Satisfaction of the BC (1.14) for [[Π ]] is not obtained for arbitrary values of ω but only for the isolated discrete eigenvalues, where the inner and outer solutions are smoothly joined.Those are found on the intersections of the curves Im(W com ) = 0 (the solution path) and Re(W com ) = 0 (the conjugate path), which, together, constitute the spectral web.
Figure 1 shows the spectral web thus obtained for the antisymmetric modes of the modified leaky mode model with the smoothed Alfvén velocity profile.The red curves represent the solution path and the blue curves represent the conjugate path.The black dots on the intersections are the eigenvalues of the leaky modes.There is a marked difference between the new eigenvalues (for = 0.2), which exhibit increasing damping rates for increasing σ ≡ Re(ω), and the old ones (for = 0), having constant damping rates satisfying (1.19), or rather (1.17).Although there is this striking difference, the first result is that damping is obtained for both models, so that the discrete spectrum of leaky modes appears to survive after our first correction.
The damping rate for the modified equilibrium model becomes larger though, and even logarithmically infinite as σ or q 1 ≡ Re(q) → ∞.This is in agreement with the damping rates obtained for the old model since larger values of q 1 imply more localisation of the modes in the region where small changes of the Alfvén velocity occur, namely, close to |x| = 1.Consequently, the problem that we noted at the end of § 1.2 is now even more prominent: infinite damping rates are present for all values of the variation of the Alfvén velocity, not only for the infinitesimal ones.An inner region with propagation in both directions is necessary because the 'outgoing wave only' BC conflicts with the BCs (1.7) and (1.8) at x = 0. Hence, for this model, there can be no approach of the situation with outgoing waves only for the trivial equilibrium with constant Alfvén velocity everywhere.This already suggests that the logarithmic damping rates are more due to this impossible limit than to a genuine physical effect.
There appears to be another problem with the spectral web method applied thus far.Blindly applying it would also include the 'false eigenvalues' for real ω, indicated by the crosses in figure 1  false ones where the solution on the inner or outer interval accidentally vanishes at the interface at |x| = 1.They are easily tracked by the value of ξ(1) and to be eliminated when it vanishes.Clearly, the 'solution' shown in figure 2(b) vanishes at |x| = 1 and, therefore, fails to connect to an exclusively outwardly propagating mode of non-vanishing amplitude.Nevertheless, these 'solutions' are included in our figures for reasons that will become clear in the next section.
A genuine leaky mode eigenfunctions for the modified equilibrium model is shown in figure 2(c).Note that the derivative ξ rapidly changes in the interval (1 − , 1) to smoothly join the value of ξ (1): the discontinuity of the old model (with = 0) is gone.
The spectral web and associated solutions for the symmetric modes are shown in figures 3 and 4.These modes already exhibit damping from the very first mode onward (n = 0, q 1 = 0), in agreement with (1.18) for the old equilibrium model with = 0. Otherwise, their behaviour is very similar to that of the antisymmetric modes.
Finally, figure 5 shows the spectral web in the opposite case of incoming waves only ( α = 0) for the antisymmetric modes.As we have already noted for the old = 0 equilibrium model, these modes are unstable.For = 0, they exhibit the same peculiarity as the leaky modes, but for incoming waves now: the imaginary part of their frequency tends to infinity as the variation of the Alfvén velocity tends to zero, either from above or from below.This implies infinite growth rates in those limits.Clearly, smoothing the Alfvén velocity profile alone is not enough to produce physically meaningful solutions.
As a parenthetical remark, the discontinuity of ξ while Π is continuous in the sharp-boundary model, that forced us to do a careful study with respect to the limiting process of smooth to discontinuous equilibrium profiles, may evoke the question whether  In (b), the coinciding real and imaginary components of the solutions on the interval (0, 1) have been shifted up or down a few plotting units to distinguish them; they do not match with the complex outgoing wave components for x ≥ 1. all surface mode models (that abound in coronal plasma physics) should be subjected to a similar scrutiny?The answer to that question is: yes, they should.That limit did not invalidate the sharp-boundary model, but it did reveal some important subtleties that could easily have been overlooked otherwise.As a pertinent other example, in the pioneering investigation (Rosenbluth 1959) of the sharp-boundary pinch, it was shown that a singularity ω = ω A (x) ≡ k B = 0 (that was only much later on realised to be the Alfvén continuum extending up to the marginal stability point) has an enormous stabilising effect because of induced skin currents at that position.The result is that, with respect to stability, the wall is effectively moved inward to that singular point.This was proved to be general by Newcomb (1960) investigating the stability of a diffuse linear pinch.That singularity divides the radial interval into independent subintervals that can be studied separately with respect to stability.The consequence of this is that a crude surface mode model of a pinch surrounded by a vacuum magnetic field completely misses this stabilising effect.Moreover, it paved the way for the early success of the Rijnhuizen screw pinch experiment in the Netherlands where replacing the external vacuum field by a tenuous plasma with force-free currents gives another boost to the stability.The theoretical stability analysis for force-free currents j = αB with constant α can be found in Goedbloed (1971b).With respect to the above question on the need to investigate the effect of smoothing the discontinuity in the boundary layer, this shows that this stabilising effect would have been missed completely otherwise.

Admitting reflected waves
Reviewing the text of § 1.2, it appears most logical to question the procedure where it is stated that we need to 'restrict the waves in the outer regions to be outgoing only' and that this implies that 'it becomes imperative to admit complex values of q and ω.' In fact,  In (b), the coinciding real and imaginary components of the solutions on the interval (0, 1) have been shifted up or down a few plotting units to distinguish them; they do not match with the complex outgoing wave components for x ≥ 1. why would one deviate at all from the standard procedure to solve differential equations by propagating a solution from point to point assuming continuity?This requires exploiting both independent solutions, not just one.In addition, even when a wave mainly propagates in one direction, if inhomogeneities or a boundary are present, reflected waves are excited so that soon there is a multitude of outward and inward waves.
Thus, let us redo the leaky mode problem, again assuming k 2 0 = 0, dropping the condition of 'outgoing waves only' ( β = 0) and just assuming continuity of both ξ and Π at |x| = 1.As before, for the modified Alfvén velocity profile, this implies that the solution (1.4) for ξ(x) is replaced by the numerical solution imposing the BC (1.7) or (1.8) at the origin.For the solution in the external region, we can now exploit the general solution (1.5) where the coefficients α1,2 and β1,2 are to be determined by the two BCs at x = 1.This involves some straightforward algebra, giving the following expressions for these coefficients in terms of the inner boundary values ξ 1,2 (1) and Π 1,2 (1), omitting the arguments (1): where the upper sign refers to α1 and β1 and the lower sign to α2 and β2 , whereas the constants η 1 and η 2 are expressions in terms of Π 1 and Π 2 : A more transparent representation of the solution in the external region is obtained by transforming the expression (1.5) for ξ in terms of α and β into one where the exponential factors become unity at x = 1: ξ = Âe iq(x−1) + Be −iq(x−1) , where Â ≡ αe iq and B ≡ βe −iq . (3.3) The coefficients Â and B, expressed in terms of the boundary values ξ 1,2 and η 1,2 of the internal solution, read For real frequencies (ν = q2 = 0), by proper choice of the phase angle of the internal solution, both ξ and ξ will become real functions, so that ξ 2 = 0 and η 1 = 0 whereas 2 ) 1/2 and φ = ψ = arctan (η 2 /ξ 1 ). (3.5) For these frequencies, the outgoing and incoming wave components, represented by the coefficients Â and B, are in perfect balance and there is no freedom to set B = 0 to satisfy the leaky mode condition of 'outgoing waves only'.Hence, assuming continuity of ξ and Π at |x| = 1, there is no restrictive dispersion equation like (1.15) anymore.Equations (3.1) yield a solution for any pair of values {q 1 , q2 }, i.e. for any complex value of ω.In other words, the discrete spectra of leaky modes illustrated in figure 1 or 3, and also of unstable modes illustrated in figure 5, are swallowed by a continuous spectrum covering the whole complex ω-plane.It is schematically indicated in figure 6, obtained by running the numerical code for a coarse 100 × 100 rectangle in the ω-plane (it becomes a black square for finer grids).This result is no surprise since [[Π ]] ≡ 0 and, hence, W com = 0 is a direct consequence of the Spectral Web method for continuous continuation of the solutions.
The consequences of this two-dimensional continuum of solutions, that are no longer constrained by the a priori condition of 'outgoing waves only', are as follows.
(i) For complex q with q2 < 0 (so that ν < 0), the leaky mode eigenfunctions, shown in figures 2(c) and 4(c), will not be modified by the condition of continuity since they are, in fact, already continuous by being smoothly joined to the outgoing wave expression.This is only so because their frequencies satisfy the dispersion equation (1.15).However, with the continuity conditions on ξ and Π , these leaky modes are swallowed by the sea of other complex modes with all kinds of distributions of incoming and outgoing waves.The leaky modes are not special anymore.(ii) For real q with q2 = 0 (so that ν = 0), on the other hand, the inner part of the 'false' solutions shown in figures 2(b) and 4(b) could not be joined smoothly to the outer part when it was restricted to outgoing waves only.Dropping that constraint smooth genuine eigenfunctions are now obtained at the same frequency: the 'false' solutions are 'rehabilitated'.This is illustrated in figure 7(b,c).However, these real modes are not special either, there are numerous other real modes that do not vanish at |x| = 1.(iii) For complex q with q2 > 0 (so that ν > 0), the discrete spectrum of unstable modes shown in figure 5, is also swallowed by the complex continuum.(iv) Actually, as a consequence of the construction, all complex frequencies are part of the new continuum: 'embarrassment of riches'.All solutions illustrated thus far are FIGURE 6. Spectral Web of all unconstrained modes (counterpart of figure 1) admitting both outgoing and ingoing waves.Now the spectral web represents a continuous spectrum occupying the whole complex ω-plane (here plotted on a course grid to restrict the number of pixels).For any value of ω, the mixture of outgoing and ingoing wave components for x ≥ 1 is dictated by continuous continuation of the solution on the interval (0, 1).
not fundamentally different from solutions in other parts of the complex ω-plane, that could be obtained by just modifying their frequencies by some amount.The dispersion equation (1.15) no longer applies and the discrete spectrum is replaced by the continuous spectrum covering the full ω-plane.There is no restriction at all anymore on the frequencies of the modes.Hence, in the dynamics of coronal flux tubes, the leaky modes have lost their distinguishing feature: they are just modes with β = 0 surrounded by infinitely many modes with small values of β.

Eliminating infinite energy flows at infinity
Having demonstrated that there is no mathematical objection against the concept of leaky modes in coronal magnetic flux tubes, just that they form a tiny collection in an infinite sea of other modes with complex frequencies, it remains to demonstrate that there is, actually, no physical mechanism operating in such environments that could preferably (by some resonance) excite them.To that end, we will study the asymptotic behaviour of the modes at large distances, or rather at infinity, being the only location where the outer world is permitted to enter in this model.
In figure 8, the eigenfunctions are shown on a large interval for the two particular frequencies selected for illustration in this paper.In figure 8 Since both outgoing and ingoing wave components are now permitted, the outer parts of these solutions (for x ≥ 1) are smoothly joined to the inner parts (for x < 1).The eigenvalues of these solutions are situated in the real part of the spectral web continuum shown in figure 6.The real and imaginary parts of the solutions now coincide over the whole interval (which implies that the solutions can be renormalised to become purely real).ξ(x, t) = αe −q 2 (x− bt) e iq 1 (x− bt) + βe q2 (x+ bt) e −iq 1 (x+ bt) , with b = 1.(4.1) Clearly, for fixed time t, the amplitude | ξ | explodes exponentially as | α|e −q 2 x for damped modes (q 2 < 0), and as | β|e q2 x for unstable modes (q 2 > 0).It appears that we need to construct a physical mechanism producing such infinities.That mechanism could only be a powerful pumping mechanism at infinity, for x → ∞ producing an exponentially infinite pumping power, i.e. an energy flow per unit time in the normal, or opposite, direction (n ≡ e x ), given by the limit of the expressions x (all leaky modes, q2 < 0), x (all unstable modes, q2 > 0).
Here, 'all leaky modes' is used in the generalised sense of all unconstrained modes with a finite damping rate, with frequencies covering the lower half of the complex ω-plane.
Obviously, that includes the restricted ones that we are criticising here.Likewise for 'all unstable modes' with a finite growth rate, with frequencies covering the upper half of the complex ω-plane.Such miraculous pumping mechanisms do not exist in nature.
Consequently, the continuum of modes in coronal flux tubes with a finite imaginary part (ν = 0) should be rejected.What remains is the continuum of fast magneto-sonic modes (ν = 0) with frequencies along the real axis of the ω-plane.They may be considered as the limit of "wall"-confined modes where the wall position w → ∞.This is illustrated in figure 9 for the axisymmetric modes, and in figure 10 for the symmetric modes, on a large but finite interval.In the top frames the discrete spectra of wall-confined modes are shown to approach the continuum of real unconstrained modes.This justification of the approach of the discrete modes on a finite interval to the continuum of modes on the infinite domain is analogous to the quantum mechanical description of the wave function and associated energy levels of a particle in a box approaching the continuous spectrum of energy levels for a free particle.This has been noted before by Andries & Goossens (2007), however with the wrong conclusion that these real modes could be used in a construction of the damped leaky modes, termed 'quasi-modes' by them.This necessarily violates the above conclusions on energy flows, let alone the absence of any tuning mechanism that could revive the discrete resonances of the leaky modes.
It should be noted that the continuum of real-frequency fast magnetosonic modes constitutes a continuous spectrum in the usual sense in that the corresponding improper eigenfunctions are non-square integrable.Their potential energy W increases linearly with the size L x (= w) in the normal direction because, for a homogeneous outer plasma, the energy per unit length is the same for all sections of unit length.Obviously, this constant linear increase is to be replaced by some other increase if that plasma is not homogeneous but, eventually, W → ∞ as L x → ∞.This infinity is, in a sense, a passive one corresponding to the free outflow condition that was originally envisaged to be provided by the leaky modes.However, the infinity P → ∞ as L x → ∞ of the power required to excite the leaky modes is of a completely different kind.The required exponentially large energy outflow should actively be provided by some dynamical pumping system.Dynamical systems of this kind do not exist in nature.For the present model of a flux tube, with vanishing pressure (c = 0) and exclusively normal wave propagation (k 0 = 0), the slow magnetosonic continuum and the Alfvén continuum, together with their associated cluster sequences of discrete modes, have disappeared (or, rather, condensed) in the origin ω = 0 so that only the fast magnetosonic subspectrum with real frequencies remains.With the two mentioned approximations, this subspectrum is complete: it consists of the usual cluster sequences of discrete fast modes accumulating at ω = ±∞ if the radial interval (1, w) is finite (illustrated in figures 9 and 10) and it condenses into a continuum covering the whole real ω-axis in the limit w → ∞.There are no additional fast discrete modes or quasi-modes elsewhere in the complex ω-plane.Note that, although one could restrict the analysis to finite radial intervals, the idealisation to an infinite size is a useful construct since it makes the whole machinery of distribution theory, in particular developed to justify the use of the improper δ-functions in quantum mechanics, also available for classical MHD.Of course, the present reduction to the fast magnetosonic subspectrum does not eliminate the widely analysed plethora of discrete modes in coronal flux tubes.In particular, dropping the two mentioned approximations, the slow and Alfvén continua and their cluster frequencies reappear in the finite frequency domain.In addition, some further distinctions result from the two possibilities that the relation (1.6) between the wavenumbers q and q permits.When k 0 = 0, depending on the values of the tangential wavenumbers, q is either purely real or purely imaginary.The latter case corresponds to the so-called 'body modes' which are evanescent in the outer region and, hence, have finite energies, also for the infinite radial sizes.(Ironically, for that reason, the exponentially increasing component is routinely discarded in the construction of those discrete modes.)Which of the different kinds of spectra applies depends on the values of k ⊥ and k chosen.
Another spin-off from the enormous efforts on quantum mechanical systems is the extensive literature on the mathematical methods used to analyse scattering and absorption phenomena in atomic and subatomic physics.This has been summarised and applied to MHD in the thesis of Keppens (1994).The methods are quite general and provide a convenient way of distinguishing genuine discrete modes from quasi-modes.The latter should be considered as a superposition of the genuine discrete or continuum modes of the spectrum proper.They apply, for example, to the well-known Alfvén quasi-modes which admit absorption of the energy supplied by an external driver on a certain time scale set by the imaginary part of the frequency of the quasi-mode.In terms of the scattering matrix S, which is the ratio of the outgoing to the incoming wave amplitudes, this happens when |S| < 1.As shown by (3.5), the fast magnetosonic modes along the real ω-axis are represented by points in the complex S-plane on the 'elastic scattering circle', |S| ≡ |B|/|A| = 1.Of course, they do not have absorption, but the question is whether they admit absorption by superposition of all the modes of the fast subspectrum.The present analysis shows that such a superposition of real-frequency fast modes cannot lead to the existence of damped quasi-modes: energies that linearly increase with the size of the plasma cannot produce energies that exponentially increase with that size.
What about the extensive analysis of fast leaky modes by Cally (2003Cally ( , 2006) ) and of the initial value problem (IVP) for those modes by Ruderman & Roberts (2006b,a)?The latter appear to prove that damped fast quasi-modes come about from the Laplace transform of the time-dependent MHD equations, with appropriate approximations, folding in reasonable initial data, constructing the resolvent operator, and subsequently applying the inverse Laplace transform to construct the desired time dependence of the perturbations.This involves deforming the Laplace contour down into the complex ω-plane in such a way that most of the integral becomes negligible and only the contributions of the encircled poles remain.Hence, the only question to be addressed is: are there such poles (away from the real axis) or not?The IVP for the MHD equations, appropriate to the present problem, is discussed extensively, and applied to Alfvén wave damping in § § 10.1-10.4 and to absorption and heating by that wave in Chapter 11 of our textbook (Goedbloed et al. 2019).(Those sections remain valid; as stated previously, only § 10.5 on leaky modes should be deleted.)The crucial expression for the Laplace transformed variable ξ (indicated by the tilde, note the difference with the hat) is the integral over the radial coordinate x involving the Green's function given by (10.17) of our textbook, that we reproduce here: Here, R represents the initial data, whereas the constituent functions Γ and (indicated by the symbols T and D in the quoted paper) of the Green's function G are constructed from the two elementary solutions U 1 and U 2 of the homogenous ODE that satisfy either the left or the right BC.For the spectral problem, the dispersion function , defined by is of particular significance since its zeros determine the eigenvalues of the problem.At this point, the fatal error is made, in the analysis by Cally (2006) by assigning physical significance to these eigenvalues and in the analysis by Ruderman & Roberts (2006b) by exploiting them in the solution of the IVP.In constructing the solutions U 1 and U 2 , they omit the Hankel function H (2)  1 (corresponding to the incoming wave contribution β exp(−iqx) of the plane slab model) in order to satisfy the 'outgoing waves only' BC.(In the analysis by Ruderman & Roberts (2006b) this is justified differently, namely, by requiring that the solution tends to zero in the limit x → ∞ in the upper part of the complex ω-plane, but the end result is the same: one of the fundamental solutions of the differential equation is discarded.)In this way, they obtain what they have put in, namely, the spurious discrete leaky modes (illustrated in figures 1-4 for the plane slab model).However, at the same time, they lose the actual complete fast mode subspectrum (illustrated in figures 9 and 10) since all of those modes consist of perfectly balanced incoming and outgoing wave components; they are on the 'elastic scattering circle', |S| = 1, as shown by (3.5).
Consequently, the preconceived notion of leaky modes that are supposed to be generated automatically by open boundaries, and that imply the sudden transition of standing towards propagating waves at the rather arbitrary position of what one assigns to be the boundary of the actual flux tube, is self-contradictory.The attempt to construct them by analytic continuation from the eigenfunctions of the problem fails completely since the two kinds of modes refer to different assumptions on the BCs of the system.The leaky modes obtained from the 'outgoing waves only' BC require the existence of powerful pumps at infinity that do not approximate any known physical system.On the other hand, for the fast eigenmodes of the conservative MHD model of the flux tube no such unphysical assumptions are necessary.The only approximation that could be questioned there is the transition of the cluster spectrum of discrete fast modes to a continuous spectrum when the size of the system is imagined to become infinite.However, this Gedankenexperiment has good papers thanks to the analogy with the quantum mechanical theory of free electrons.
It will not have escaped the attentive reader that the corrected initial value approach described here is identical to that presented by Goedbloed (1998) once more rejecting the conjecture by Grad (1973) that there are additional (real) magnetosonic continua associated with the apparent singularities of the denominator expression in the ODE for ξ .This is no accident: the actual fast discrete cluster spectrum (for w < ∞) or continuous spectrum (for w → ∞) is complete, additional real or complex modes cannot be found by solving the IVP.

Conclusion
By three modifications of the standard plane slab model for coronal magnetic flux tubes, we have shown the following.
(i) A smooth Alfvén velocity profile, instead of a discontinuous one, cures the shortcoming of the standard model with discontinuous derivatives in the solutions, but it also makes the unphysical divergence for small variations of that profile more manifest.(ii) Replacing the unfounded restriction to 'outgoing waves only' at the arbitrary position of the plasma interface between the two plasmas with different densities by the regular condition of continuity of the normal plasma displacement ξ and of the perturbation Π of the total pressure causes the leaky modes to lose their identity by immersion in a continuous spectrum of mathematically permitted modes that cover the complete complex ω-plane.(iii) The complex continuum thus constructed has no physical meaning, except on the real axis, since all corresponding modes require an exponentially large pumping power at infinity.
One could object that most of the references cited in § 1.1 refer to cylindrical magnetic flux tubes having much more complicated dependencies on the solutions than for the plane slab.To counter that argument, one would have to replace the ODE (1.1) with the constant Alfvén velocities, or the smoothed profile (2.1) or any of the more general equilibria permitted by the differential equation derived by Goedbloed (1971a), by its cylindrical counterpart (Goedbloed 1971b), where the restriction to isothermal equilibria of the equation of Hain & Lüst (1958) was dropped.This would involve rather cumbersome numerical calculations, replacing the ordinary Bessel functions for the internal region by their numerical counterparts, but keeping the Hankel functions for the external region.However, such a laborious investigation is not necessary since it would only obscure the central result already obtained for the sharp-boundary cylindrical case, extensively studied in the literature cited in § 1.1.That result is a dispersion equation governing the discrete spectrum of leaky modes with complex frequencies which is fully analogous to the plane slab counterpart (1.15).For the cylindrical case, that dispersion equation is an equality with ordinary Bessel functions on the left-hand side and Hankel functions on the right-hand side, see Eq. (A 8) of Meerson et al. (1978) which is probably the first time this equation explicitly appears in the literature.(We owe this remark to a careful referee.)Replacing the ordinary Bessel functions by numerical solutions for the smoothed Alfvén velocity would be a subordinate detail.The essential point of all leaky mode calculations is the choice to restrict the possible solutions in the external region to one of the two Hankel functions only, namely, the one that represents the outgoing wave component.That restriction leads to a complex discrete spectrum, which disappears in the complex continuum when all combinations of the two Hankel functions are admitted.In turn, that mathematical continuum is forbidden by the physical condition that the associated power at infinity remains finite.The asymptotic form of the Hankel functions for large argument is actually of the same form as the exponential functions in the plane slab case, except for a multiplier r −1/2 .Replacing the length L y by the integral L θ = 2πr over the angle θ , that factor is canceled and the resulting expressions are identical to (4.2).The crux of § 3 is that there is no compelling reason to make that restriction, and of § 4 that restricting all possible combinations of the elementary solutions in the external region to those that correspond to a finite power at infinity eliminates all non-real values of ω.Hence, there is no dispersion equation for complex frequencies, and there are no associated resonances either.
In conclusion, exit leaky modes in coronal magnetic flux tubes.
& Goossens (2007), Terradas et al. (2005) and Wang et al. (2023)) is: yes, the real fast continuum modes provide the physical solution of the EVP, but the complex leaky modes are quasi-modes consisting of a superposition of the singular continuum modes appearing in the solution of the IVP.To validate this, the analogies with Landau damping in the IVP in the velocity space description of plasmas (Landau 1946), dissipationless damping of plasma oscillations in a cold plasma (Sedláček 1971) or the excitation of Alfvén quasi-modes in the MHD description (Tataronis & Grossmann 1973;Chen & Hasegawa 1974;Poedts & Kerner 1991) are cited.However, these analogies are not valid.
The mentioned examples with damping concern resonances at singular points of the continuum that are located right in the middle of the physical interval (in ordinary or velocity space), whereas in the present case the resonances (if existing at all) would have to be located at infinity.Moreover, in those examples, both the basic singular continuum modes and the damped quasi-modes are solutions of the pertinent boundary value problem (BVP) applying the same BCs to both sets of modes.As already demonstrated in § 4, the present fast continuum modes and the leaky modes are solutions of different BVPs.At this point, a note on the nomenclature of the modes is in order.The adjective 'leaky' would be appropriate for the regular fast continuum modes (infinite L ≡ L x ) or the discrete modes (finite L) since the energy carried by these modes can leak out of the ends of the flux tube at x = ±L if the waves are travelling.This requires excitation somewhere, since the eigenmodes themselves (improper or proper) are standing waves, not propagating energy.The crucial observation is another one though: whereas the total energy of the fast continuum modes is infinite (∼ L), leakage at the ends only concerns the finite energy content of the travelling waves on the outermost part of the flux tube that is passively transported out of the system at the Alfvén speed.Quite the contrary in the case of the leaky modes: There, the adjective 'leaky' is a very misleading misnomer since the process of excitation of these modes is not passive at all but requires active, exponentially large, pumping power on the outside of the flux tube (the adjective 'power-driven' would have been more descriptive).For the construction of these modes, the infinite energy of the fast continuum modes is irrelevant, since the main result of § 4 does not change at all if the external interval is limited to a finite size, L = 1000 is already big enough to demonstrate that the pumping power required for the existence of leaky modes in this system, given by (4.2), cannot be provided by any physical system in the universe.
Finally, it remains to address the question on why the signatures of leaky modes (the real frequencies and approximate damping rates of the most global ones) nevertheless appear to enter in the calculation of the initial phase of the temporal evolution.This feature is clearly seen in the illustrations of the papers by Terradas et al. (2005, figures 4 and 6), by Andries & Goossens (2007, figures 2 and 3) and in the recent paper by Wang et al. (2023, figures 2-5).(We are indebted to one of the referees to bring this paper to our attention; it provides an extensive list of relevant references and includes a discussion of the implications of the theory for the interpretation of the abundant observations in solar coronal seismology.)In evident contradiction to our discussion, in the quoted paper by Andries and Goossens the time-dependent calculations of the initial evolution of the modes is claimed to provide the connection between the complex leaky modes and the continuous spectrum of real fast magnetosonic waves in an infinite medium.This is done by an explicit construction of the solution of the IVP for the complete spectrum of MHD waves (including the slow and Alfvén subspectra, quite a monumental expression), exploiting what is called the spectral measure to fold-in the spatial and temporal distribution of the initial data.However, this spectral measure is constructed by exploiting the dispersion function of the leaky modes, i.e. by folding in the desired answer of the IVP, similar to the procedure exploited in the analysis of Ruderman & Roberts (2006b) that we rejected in § 4. That does not take away that in the initial responses of the time-dependent studies by Terradas et al. (2005) and by Wang et al. (2023), peaks appear at the leaky mode frequencies.There is no contradiction here: the system is perturbed by waves that propagate outwardly, as if leaky modes are excited.Inevitably, they excite regular fast waves at positions down the loop.These reflected waves arrive at the original region of excitation at some later time.For example, at t ∼ 10/v A reflected waves from x ∼ 10 arrive carrying the implicit message 'The requested amplitudes for leaky modes cannot be produced here due to the excessive power requirements'.Consequently, the leaky mode response is quenched and the evolution continues with the excitation of regular fast magnetosonic continuum or discrete waves, nor requiring any external power source.Hence, a leaky mode response over the full length of the loop, of the kind discussed and rejected in the present paper, cannot possibly be excited.In addition, remember that the quoted examples of genuine continuum damping do not refer to initial perturbations but to the long-time asymptotic behaviour of the system.The failed attempt to construct a leaky mode counterpart may explain why these modes were actually never observed in solar coronal seismology, as noted by Wang et al. (2023).
FIGURE 1. Spectral web of the antisymmetric leaky modes (dots) and false solutions (crosses) for the smoothed Alfvén velocity profile b(x) ≡ B/ √ ρ(x).In this paper, that velocity is chosen to vary from b(0) = δ to b(1) = 1.0 over the interval 1.0 − ≤ x ≤ 1.0, where δ = 0.4 and = 0.2.The marked modes are constrained to have outgoing wave components only.The dots represent the discrete spectrum of the antisymmetric leaky modes thus obtained.The eigenvalues for the original model, with a jump in the Alfvén velocity ( = 0), are represented by the symbol •.

FIGURE 2 .
FIGURE 2. Two antisymmetric solutions corresponding to the spectral web depicted in figure 1.(a) Profile of the smoothed Alfvén velocity b(x).(b) False solution corresponding to one of the crosses in figure 1 (σ = 2.5734, ν = 0.0).(c) Leaky mode solution corresponding to one of the dots in figure 1 (σ = 3.6168, ν = −0.21830).In (b), the coinciding real and imaginary components of the solutions on the interval (0, 1) have been shifted up or down a few plotting units to distinguish them; they do not match with the complex outgoing wave components for x ≥ 1.

FIGURE 3 .
FIGURE 3. Spectral web of the symmetric leaky modes (dots) and false solutions (crosses) for the smoothed Alfvén velocity profile b(x) ≡ B/ √ ρ(x).The marked modes are constrained to have outgoing wave components only.The dots represent the discrete spectrum of the symmetric leaky modes.The eigenvalues for the model with a jump in the Alfvén velocity are represented by the symbol •.

FIGURE 4 .
FIGURE 4. Two symmetric solutions corresponding to the spectral web depicted in figure 3. (a) Profile of the smoothed Alfvén velocity b(x).(b) False solution corresponding to one of the crosses of figure 3 (σ = 3.2446, ν = 0.0).(c) Leaky mode solution corresponding to one of the dots of figure 3 (σ = 4.3381, ν = −0.22939).In (b), the coinciding real and imaginary components of the solutions on the interval (0, 1) have been shifted up or down a few plotting units to distinguish them; they do not match with the complex outgoing wave components for x ≥ 1.

FIGURE 5 .
FIGURE 5. Spectral web of the discrete antisymmetric overstable (= unstable) modes (dots) and false solutions (crosses) for the smoothed Alfvén velocity profile b(x) ≡ B/ρ(x).The marked modes are constrained to have ingoing wave components only.The dots represent the discrete spectrum of the antisymmetric ingoing waves.The eigenvalues for the model with a jump in the Alfvén velocity are represented by the symbol •.The corresponding solutions are omitted since they are identical to figures 2 and 4 when ω is replaced by ω * and the real and imaginary components of both ξ and Π are interchanged.

FIGURE 7 .
FIGURE 7. The unconstrained solutions for the real frequencies corresponding to the constrained 'false' solutions of figures 2(b) and 4(b).(a) Profile of the smoothed Alfvén velocity b(x).(b) Antisymmetric solution (σ = 2.5731, ν = 0.0).(c) Symmetric solution (σ = 3.2446, ν = 0.0).Since both outgoing and ingoing wave components are now permitted, the outer parts of these solutions (for x ≥ 1) are smoothly joined to the inner parts (for x < 1).The eigenvalues of these solutions are situated in the real part of the spectral web continuum shown in figure6.The real and imaginary parts of the solutions now coincide over the whole interval (which implies that the solutions can be renormalised to become purely real).

FIGURE 8 .
FIGURE 8. Unconstrained antisymmetric modes continued on a large interval.(a) Profile of the Alfvén velocity b(x) on a large interval.(b) The 'false' solution of figure 2(b), or rather the solution of figure 7(b), on the large interval (σ = 2.5731, ν = 0.0); that solution now corresponds to a genuine eigenvalue.(c) The leaky mode of figure 2(c) on the large interval (σ = 3.6168, ν = −0.21830).The eigenvalues of these solutions are situated in, respectively, (b) the real (ν = 0) and (c) the complex (ν = 0) parts of the spectral web continuum shown in figure 6.
FIGURE 9. (a) Discrete spectrum of wall-confined real antisymmetric modes for a large interval.(b) Profile of the Alfvén velocity b(x).(c) The 31st eigenfunction, indicated by the arrow in panel (a) (σ = 2.6018, ν = 0.0).This discrete eigenfunction and eigenvalue approaches the continuum ones of the unconstrained solution shown in figure 8(b).Clearly, the influence of the wall position becomes ever smaller for ever larger values of the wavenumber.
FIGURE 10.(a) Discrete spectrum of wall-confined real symmetric modes for a large interval.(b) Profile of the Alfvén velocity b(x).(c) The 39th eigenfunction, indicated by the arrow in panel (a) (σ = 3.2320, ν = 0.0).This discrete eigenfunction and eigenvalue approaches the continuum ones of the unconstrained solution shown in figure 7(c) continued on the larger interval.For these highly oscillatory solutions the position of the wall becomes insignificant.