Leaky modes in coronal magnetic flux tubes revisited

Abstract

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.

1. 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).

1.1. 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 $\omega$ 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 wave-numbers, 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 $\omega$ 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.

1.2. 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 \equiv \gamma p/\rho = 0$). The magnetic field $B$ is considered constant throughout (otherwise there would be no equilibrium) and the density jumps from $\rho$ to $\hat {\rho }$. (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 \equiv B/\sqrt {\rho }$ to $\hat {b} \equiv B/\sqrt {\hat {\rho }}$ 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$, $\hat {\rho }$ and $\hat {b}$, effectively equating them (and also the constant $\mu _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 $\hat {b} \equiv 1$ by normalisation. Consequently, the equilibrium is a one-parameter family described by the parameter $\delta \equiv b/\hat {b}$, where $\delta < 1$ for a higher density and a smaller Alfvén velocity in the inner than in the outer region and $\delta > 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 $\xi \equiv {\boldsymbol {n}} \boldsymbol {\cdot } \boldsymbol {\xi }$ (with ${\boldsymbol {n}} \equiv {\boldsymbol {e}}_x$) of the plasma displacement vector $\boldsymbol {\xi }$, obtained by Goedbloed (1971a), where we neglect the gravitational terms whereas the slow magnetosonic modes have disappeared by the assumption $c = 0$:

(1.1)\begin{equation} b^2 \frac{{\rm d}}{{\rm d}x} \left(\frac{\omega^2 - k_\parallel^2 b^2}{\omega^2 - k_0^2 b^2} \,\frac{{\rm d}\xi}{{\rm d}x}\right) + (\omega^2 - k_\parallel^2 b^2) \xi = 0 ,\quad k_0^2 \equiv k_\perp^2 + k_\parallel^2. \end{equation}

The magnetic field is taken in the $z$-direction, so that the tangential wavenumbers are $k_\perp \equiv k_y$ and $k_\parallel \equiv k_z$. The normal modes for this model are then expressed by

(1.2)\begin{equation} \xi({\boldsymbol{r}};t) = A \exp [{\rm i}(q x + k_\perp y + k_\parallel z - \omega t)], \end{equation}

from which the other components of the vector $\boldsymbol {\xi }$ follow by algebraic expressions in terms of $\xi$ and its first derivative $\xi '$. 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, $\omega ^2 = k_\parallel ^2 b^2$, can then be extracted so that (1.1) simplifies to a description for the fast modes only:

(1.3)\begin{equation} b^2 \frac{{\rm d}^2 \xi}{{\rm d}x^2} + (\omega^2 - k_0^2 b^2) \xi = 0 ,\quad k_0^2 \equiv k_\parallel^2 + k_\perp^2 . \end{equation}

A similar equation holds $\hat {\xi }$ in the external region, with Alfvén velocity $\hat {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:

(1.4)\begin{gather} \xi = \alpha{\rm e}^{{\rm i} q x} + \beta{\rm e}^{-{\rm i} q x} ,\quad \omega^2 = (q^2 + k_0^2){}b^2, \end{gather}

(1.5)\begin{gather}\hat{\xi} = \hat{\alpha}{\rm e}^{{\rm i} \hat{q} x} + \hat{\beta}{\rm e}^{-{\rm i} \hat{q} x} ,\quad \omega^2 = (\hat{q}^2 + k_0^2){}\hat{b}^2. \end{gather}

Of course, the two dispersion equations refer to the same value of $\omega$ so that the wavenumbers $\hat {q}$ and $q$ are related,

(1.6)\begin{equation} \hat{q}^2 = \delta^2 {}q^2 - (1 - \delta^2) {}k_0^2. \end{equation}

The contribution $k_0^2$ 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 \ge 0$ since the dependence for $x \le 0$ follows from symmetry. The antisymmetric (sine-shaped) solutions are subject to the BC

(1.7)\begin{equation} \xi(0) = 0 \quad\Rightarrow\quad \alpha ={-} \beta, \end{equation}

the symmetric (cosine-shaped) solutions are subject to the BC

(1.8)\begin{equation} \xi'(0) = 0 \quad\Rightarrow\quad \alpha = \beta. \end{equation}

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 ${\rm J}_1$ for the $m = 0$ sausage modes and by ${\rm 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 $\hat {q}$ and $\omega$,

(1.9)\begin{equation} \hat{q} \equiv \hat{q}_1 + {\rm i}\hat{q}_2 ,\quad \omega \equiv \sigma + {\rm i}\nu , \end{equation}

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

(1.10)\begin{equation} \hat{\xi}(x,t) = \hat{\alpha}{\rm e}^{- (\hat{q}_2 x - \nu t)}{\rm e}^{{\rm i}(\hat{q}_1 x - \sigma t)} + \hat{\beta}{\rm e}^{(\hat{q}_2 x + \nu t)}{\rm e}^{{- {\rm i}}(\hat{q}_1 x + \sigma t)} . \end{equation}

Since, by symmetry, we need to consider only $x \ge 0$, and time progresses, it is evident that the following constraints should be applied:

(1.11)\begin{equation} \hat{\beta} = 0 \quad (\text{for outgoing waves}),\quad \hat{\alpha} = 0 \quad (\text{for incoming waves}). \end{equation}

For ‘leaky modes’, the first restriction applies.

The BCs at the interface between the two plasmas express continuity of the normal displacement $\xi$ and of the perturbation $\varPi$ 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)\begin{equation} \varPi ={-} \frac{B^2 (\omega^2 - k_\parallel^2 b^2)}{\omega^2 - k_0^2 b^2} \xi' ={-} B^2 (1 + k_\perp^2/q^2) \xi'. \end{equation}

The BCs for the outgoing ‘leaky waves’ now become

(1.13)\begin{gather} \xi(1) = \hat{\xi}(1) \quad \Rightarrow \quad \alpha ({\rm e}^{{\rm i}q} \mp {\rm e}^{-{\rm i}q}) = \hat{\alpha} {\rm e}^{{\rm i}\hat{q}}, \end{gather}

(1.14)\begin{gather}\varPi(1) = \hat{\varPi}(1) \quad \Rightarrow \quad \alpha q (1 + k_\perp^2/q^2) ({\rm e}^{{\rm i}q} \pm {\rm e}^{-{\rm i}q}) = \hat{\alpha}\hat{q} (1 + k_\perp^2/\hat{q}^2) {\rm e}^{{\rm i}\hat{q}} , \end{gather}

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)\begin{equation} \frac{1 \pm {\rm e}^{{-}2{\rm i}q}}{1 \mp {\rm e}^{{-}2{\rm i}q}} = d \quad\Rightarrow\quad {\rm e}^{{-}2{\rm i}q} ={\mp} \frac{1 - d}{1 + d}, \quad \text{where}\ d \equiv \frac{\hat{q}(1 + k_\perp^2/\hat{q}^2)}{ q (1 + k_\perp^2/q^2)}, \end{equation}

and recall that $\hat {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)\begin{equation} d \rightarrow d(k_0 = 0) \equiv \hat{q}/q = b/\hat{b} \equiv \delta. \end{equation}

(One of the square root solutions of (1.6) also yields $\hat {q}/q = - \delta$, but this case can be covered by flipping the sign of $q$ and assigning a different complex value to $\alpha$.) The resulting dispersion equation can be solved explicitly to give

(1.17)\begin{align} & \text{- antisymmetric modes:} \quad q = \left(n' + \displaystyle{\displaystyle\frac{1}{2}}\right) {\rm \pi}+ {\displaystyle\frac{1}{2}} {\rm i} \ln \left|\frac{1 - \delta}{1 + \delta}\right| \quad(n' = 0, \pm 1, \pm 2, \ldots), \end{align}

(1.18)\begin{align} & \text{- symmetric modes:} \quad q = n'' {\rm \pi}+ {\displaystyle\frac{1}{2}} {\rm i} \ln \left|\frac{1 - \delta }{1 + \delta}\right| \quad(n'' = 0, \pm 1, \pm 2, \ldots). \end{align}

The two solutions may be combined into one expression that alternatingly applies to the antisymmetric (even $n$) and to the symmetric (odd $n$) solutions:

(1.19)\begin{equation} \omega = q{}b ,\quad q = {\displaystyle\frac{1}{2}} n {\rm \pi}+ {\displaystyle\frac{1}{2}} {\rm i} \ln \left|\frac{1 - \delta }{1 + \delta}\right| \quad(n = 0, \pm 1, \pm 2, \ldots). \end{equation}

The expressions are given for $\delta < 1$. For $\delta > 1$, an extra phase shift of ${\textstyle \frac {1}{2}} {\rm \pi}$ should be added to $q_1 \equiv {\rm 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 ($\hat {\alpha } = 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 $\hat {b}$. Moreover, all these waves are damped with the same logarithmic factor depending only on the jump $\delta$ of those velocities. However, there is a problem. Quite counterintuitively, the damping rates explode logarithmically in the limit $\delta \rightarrow 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.

1.3. 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.

1. (i) Because of the jump of the Alfvén velocity, the differential equation (1.1), or (1.3) for $k_\perp = 0$, has discontinuous coefficients and, hence, discontinuous solutions with a jump $q \rightarrow \hat {q}$ of the radial wavenumber and an associated jump of the derivative $\xi '$, or $\xi ''$ for $k_\perp = 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?

2. (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.

3. (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.

2. 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

(2.1)\begin{equation} b(x) = \begin{cases} \delta & (\text{for }0 \le x \le 1 - \epsilon), \\ {\textstyle\frac{1}{2}}(1 + \delta) + {\textstyle\frac{1}{2}} (1 - \delta) \sin[(x - 1 + {\textstyle\frac{1}{2}}\epsilon){\rm \pi}/\epsilon] & (\text{for } 1 - \epsilon \le x \le 1), \end{cases} \end{equation}

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 $\xi (x)$ and $\hat {\xi }(x)$ in the two regions for arbitrary complex values of $\omega$, join them by applying the BC (1.13), but replacing the BC (1.14) for $\varPi$ by the computation of the complementary energy over the boundary surface:

(2.2)\begin{equation} W_{\rm com} ={-} {\displaystyle\frac{1}{2}} \int \xi^* \lbrack\!\lbrack\varPi\rbrack\!\rbrack \,{\rm d}y\,{\rm d}z ={-} {\displaystyle\frac{1}{2}} L_y L_z \,\xi^* \lbrack\!\lbrack\varPi\rbrack\!\rbrack , \quad\text{where}\ \lbrack\!\lbrack\varPi\rbrack\!\rbrack \equiv \hat{\varPi}(1) - \varPi(1). \end{equation}

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_{\rm com}$.

The procedure replacing the analysis of § 1.2 is as follows. A relevant piece (e.g. a square) of the complex $\omega$-plane is selected and the real and imaginary parts of $W_{\rm com}$ are computed by substituting the solution pairs $\{\xi,\varPi \}$ and $\{\hat {\xi },\hat {\varPi }\}$ for all values of $\omega$ scanned in that piece. The function $\xi (x)$ is obtained numerically by solving the mentioned ODE, or rather the equivalent set of first-order ODEs for $\xi$ and $\varPi$, on the inner interval $0 \le x \le 1$ by applying either one of the BCs (1.7) or (1.8) at $x = 0$. For the outer interval $x \ge 1$, the function $\hat {\xi }$ is the explicitly given solution of (1.5), with $\hat {\beta } = 0$ , whereas the associated expression for $\hat {\varPi }$ follows from the equivalent of (1.12). Satisfaction of the BC (1.13) is obtained by simply adapting the amplitude $\xi (1)$ to $\hat {\xi }(1)$ or vice versa. Satisfaction of the BC (1.14) for $\lbrack \!\lbrack \varPi \rbrack \!\rbrack$ is not obtained for arbitrary values of $\omega$ 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 ${\rm Im}(W_{\rm com}) = 0$ (the solution path) and ${\rm Re}(W_{\rm 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 $\epsilon = 0.2$), which exhibit increasing damping rates for increasing $\sigma \equiv {\rm Re}(\omega )$, and the old ones (for $\epsilon = 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.

Figure 1. Spectral web of the antisymmetric leaky modes (dots) and false solutions (crosses) for the smoothed Alfvén velocity profile $b(x) \equiv B/\sqrt {\rho (x)}$. In this paper, that velocity is chosen to vary from $b(0) = \delta$ to $b(1) = 1.0$ over the interval $1.0 - \epsilon \le x \le 1.0$, where $\delta = 0.4$ and $\epsilon = 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 ($\epsilon = 0$), are represented by the symbol $\circ$.

The damping rate for the modified equilibrium model becomes larger though, and even logarithmically infinite as $\sigma$ or $q_1 \equiv {\rm Re}(q) \rightarrow \infty$. 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 $\omega$, indicated by the crosses in figure 1 and associated with the ‘eigenfunctions’ depicted in figure 2(b), in the discrete spectrum of modes. This is a known feature of this method, noted before (Goedbloed 2018a,b; Goedbloed & Keppens 2022): along the solution path, the integrand $\xi ^* \lbrack \!\lbrack \varPi \rbrack \!\rbrack$ of $W_{\rm com}$ alternates between genuine eigenvalues where $\lbrack \!\lbrack \varPi \rbrack \!\rbrack = 0$ and 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 $\xi (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.

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 ($\sigma = 2.5734$, $\nu = 0.0$). (c) Leaky mode solution corresponding to one of the dots in figure 1 ($\sigma = 3.6168$, $\nu = -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 \ge 1$.

A genuine leaky mode eigenfunctions for the modified equilibrium model is shown in figure 2(c). Note that the derivative $\xi '$ rapidly changes in the interval $(1 - \epsilon, 1)$ to smoothly join the value of $\hat {\xi }'(1)$: the discontinuity of the old model (with $\epsilon = 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 $\epsilon = 0$. Otherwise, their behaviour is very similar to that of the antisymmetric modes.

Figure 3. Spectral web of the symmetric leaky modes (dots) and false solutions (crosses) for the smoothed Alfvén velocity profile $b(x) \equiv B/\sqrt {\rho (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 $\circ$.

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 ($\sigma = 3.2446$, $\nu = 0.0$). (c) Leaky mode solution corresponding to one of the dots of figure 3 ($\sigma = 4.3381$, $\nu = - 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 \ge 1$.

Finally, figure 5 shows the spectral web in the opposite case of incoming waves only ($\hat {\alpha } = 0$) for the antisymmetric modes. As we have already noted for the old $\epsilon = 0$ equilibrium model, these modes are unstable. For $\epsilon \ne 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.

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) \equiv B/\rho (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 $\circ$. The corresponding solutions are omitted since they are identical to figures 2 and 4 when $\omega$ is replaced by $\omega ^*$ and the real and imaginary components of both $\xi$ and $\varPi$ are interchanged.

As a parenthetical remark, the discontinuity of $\xi$ while $\varPi$ 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 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 $\omega = \omega _{\rm A}(x) \equiv k_\parallel 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 ${\boldsymbol {j}} = \alpha {\boldsymbol {B}}$ with constant $\alpha$ 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.

3. 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 $\hat {q}$ and $\omega$.’ In fact, 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_0^2 = 0$, dropping the condition of ‘outgoing waves only’ ($\hat {\beta } = 0$) and just assuming continuity of both $\xi$ and $\varPi$ at $|x| = 1$. As before, for the modified Alfvén velocity profile, this implies that the solution (1.4) for $\xi (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 $\hat {\alpha }_{1,2}$ and $\hat {\beta }_{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 $\xi _{1,2}(1)$ and $\varPi _{1,2}(1)$, omitting the arguments (1):

(3.1)\begin{equation} \left. \begin{aligned} \hat{{\alpha}} _{1,2} & = {\textstyle\frac{1}{2}} {\rm e}^{\hat{q}_2} [(\xi_{1,2} + \eta_{1,2}) \cos\hat{q}_1 \pm (\xi_{2,1} + \eta_{2,1}) \sin\hat{q}_1 ],\\ \hat{{\beta}} _{1,2} & ={\textstyle\frac{1}{2}} {\rm e}^{-\hat{q}_2} [ (\xi_{1,2} - \eta_{1,2}) \cos\hat{q}_1 \mp (\xi_{2,1} - \eta_{2,1}) \sin\hat{q}_1 ], \end{aligned} \right\} \end{equation}

where the upper sign refers to $\hat {\alpha }_1$ and $\hat {\beta }_1$ and the lower sign to $\hat {\alpha }_2$ and $\hat {\beta }_2$, whereas the constants $\eta _1$ and $\eta _2$ are expressions in terms of $\varPi _1$ and $\varPi _2$:

(3.2)\begin{equation} \eta_1 \equiv (\hat{q}_2 \varPi_1 - \hat{q}_1 \varPi_2)/|\hat{q}|^2 \quad\text{and}\quad \eta_2 \equiv (\hat{q}_1 \varPi_1 + \hat{q}_2 \varPi_2)/|\hat{q}|^2 . \end{equation}

A more transparent representation of the solution in the external region is obtained by transforming the expression (1.5) for $\hat {\xi }$ in terms of $\hat {\alpha }$ and $\hat {\beta }$ into one where the exponential factors become unity at $x = 1$:

(3.3)\begin{equation} \hat{\xi} = \hat{A}{\rm e}^{{\rm i} \hat{q} (x - 1)} + \hat{B}{\rm e}^{-{\rm i} \hat{q} (x - 1)} ,\quad\text{where}\ \hat{A} \equiv \hat{\alpha}{\rm e}^{{\rm i} \hat{q}} \ \text{and}\ \hat{B} \equiv \hat{\beta}{\rm e}^{-{\rm i} \hat{q}}. \end{equation}

The coefficients $\hat {A}$ and $\hat {B}$, expressed in terms of the boundary values $\xi _{1,2}$ and $\eta _{1,2}$ of the internal solution, read

(3.4)\begin{equation} \left. \begin{aligned} \hat{A} = |\hat{A}| {\rm e}^{{\rm i}\hat{\phi}} , \quad |\hat{A}| & = {\textstyle\frac{1}{2}} \sqrt{(\xi_1 + \eta_1)^2 + (\xi_2 + \eta_2)^2},\quad \hat{\phi} = \arctan\frac{\xi_2 + \eta_2}{\xi_1 + \eta_1},\\ \hat{B} = |\hat{B}| {\rm e}^{-{\rm i}\hat{\psi}} ,\quad |\hat{B}| & = {\textstyle\frac{1}{2}} \sqrt{(\xi_1 - \eta_1)^2 + (\xi_2 - \eta_2)^2}, \quad \hat{\psi} ={-} \arctan\frac{\xi_2 - \eta_2}{\xi_1 - \eta_1}. \end{aligned} \right\} \end{equation}

For real frequencies ($\nu = \hat {q}_2 = 0$), by proper choice of the phase angle of the internal solution, both $\xi$ and $\hat {\xi }$ will become real functions, so that $\xi _2 = 0$ and $\eta _1 = 0$ whereas

(3.5)\begin{equation} |\hat{A}| = |\hat{B}| = {\textstyle\frac{1}{2}} (\xi_1^2 + \eta_2^2)^{1/2} \quad \text{and}\quad \hat{\phi} = \hat{\psi} = \arctan{(\eta_2/\xi_1)}. \end{equation}

For these frequencies, the outgoing and incoming wave components, represented by the coefficients $\hat {A}$ and $\hat {B}$, are in perfect balance and there is no freedom to set $\hat {B} = 0$ to satisfy the leaky mode condition of ‘outgoing waves only’.

Hence, assuming continuity of $\xi$ and $\varPi$ at $|x|= 1$, there is no restrictive dispersion equation like (1.15) anymore. Equations (3.1) yield a solution for any pair of values $\{\hat {q}_1, \hat {q}_2\}$, i.e. for any complex value of $\omega$. 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 $\omega$-plane. It is schematically indicated in figure 6, obtained by running the numerical code for a coarse $100 \times 100$ rectangle in the $\omega$-plane (it becomes a black square for finer grids). This result is no surprise since $\lbrack \!\lbrack \varPi \rbrack \!\rbrack \equiv 0$ and, hence, $W_\textrm {com} = 0$ is a direct consequence of the Spectral Web method for continuous continuation of the solutions.

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 $\omega$-plane (here plotted on a course grid to restrict the number of pixels). For any value of $\omega$, the mixture of outgoing and ingoing wave components for $x \ge 1$ is dictated by continuous continuation of the solution on the interval $(0,1)$.

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.

1. (i) For complex $\hat {q}$ with $\hat {q}_2 < 0$ (so that $\nu < 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 $\xi$ and $\varPi$, 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.

2. (ii) For real $\hat {q}$ with $\hat {q}_2 = 0$ (so that $\nu = 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$.

3. (iii) For complex $\hat {q}$ with $\hat {q}_2 > 0$ (so that $\nu > 0$), the discrete spectrum of unstable modes shown in figure 5, is also swallowed by the complex continuum.

4. (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 not fundamentally different from solutions in other parts of the complex $\omega$-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 $\omega$-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 $\hat {\beta } = 0$ surrounded by infinitely many modes with small values of $\hat {\beta }$.

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 ($\sigma = 2.5731$, $\nu = 0.0$). (c) Symmetric solution ($\sigma = 3.2446$, $\nu = 0.0$). Since both outgoing and ingoing wave components are now permitted, the outer parts of these solutions (for $x \ge 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).

4. 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(b), the real ‘false’ solution of figure 2(b), turned genuine in figure 7(b), is shown to keep oscillating with a constant amplitude. On the other hand, in figure 8(c), the complex leaky mode solution of figure 2(c) is shown to keep increasing, suggesting that, eventually, it will explode exponentially as $x \rightarrow \infty$. This is actually what happens, in agreement with the expression (1.10) of the plasma displacement for normal propagation ($k_0 = 0$), where $\sigma = \hat {q}_1$ and $\nu = \hat {q}_2$:

(4.1)\begin{equation} \hat{\xi}(x,t) = \hat{\alpha}{\rm e}^{- \hat{q}_2 (x - \hat{b}t)}{\rm e}^{{\rm i}\hat{q}_1 (x - \hat{b}t)} + \hat{\beta}{\rm e}^{\hat{q}_2 (x + \hat{b}t)}{\rm e}^{{- {\rm i}}\hat{q}_1 (x + \hat{b}t)} ,\quad \text{with}\ \hat{b} = 1. \end{equation}

Clearly, for fixed time $t$, the amplitude $|\hat {\xi }|$ explodes exponentially as $|\hat {\alpha }|\textrm {e}^{- \hat {q}_2 x}$ for damped modes ($\hat {q}_2 < 0$), and as $|\hat {\beta }|\textrm {e}^{\hat {q}_2 x}$ for unstable modes ($\hat {q}_2 > 0$). It appears that we need to construct a physical mechanism producing such infinities.

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 ($\sigma = 2.5731$, $\nu = 0.0$); that solution now corresponds to a genuine eigenvalue. (c) The leaky mode of figure 2(c) on the large interval ($\sigma = 3.6168$, $\nu = -0.21830$). The eigenvalues of these solutions are situated in, respectively, (b) the real ($\nu = 0$) and (c) the complex ($\nu \ne 0$) parts of the spectral web continuum shown in figure 6.

That mechanism could only be a powerful pumping mechanism at infinity, for $x \rightarrow \infty$ producing an exponentially infinite pumping power, i.e. an energy flow per unit time in the normal, or opposite, direction (${\boldsymbol {n}} \equiv {\boldsymbol {e}}_x$), given by the limit of the expressions

(4.2)\begin{equation} \left. \begin{aligned} P^{\rm out} & \equiv{-}\int \left(\frac{\partial\hat{\xi}}{\partial t} \varPi^* \right)^{\rm out} \,{\rm d}S \approx L_y L_z |\omega|^2 |\hat{\alpha}|^2 {\rm e}^{{-}2\hat{q}_2 x} \quad(\text{all leaky modes}, \hat{q}_2 < 0),\\ P^{\rm in} & \equiv \int \left(\frac{\partial\hat{\xi}}{\partial t} \varPi^* \right)^{\rm in} \,{\rm d}S \approx{-} L_y L_z |\omega|^2 |\hat{\beta}|^2 {\rm e}^{2\hat{q}_2 x} \quad(\text{all unstable modes}, \hat{q}_2 > 0). \end{aligned} \right\} \end{equation}

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 $\omega$-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 $\omega$-plane. Such miraculous pumping mechanisms do not exist in nature. Consequently, the continuum of modes in coronal flux tubes with a finite imaginary part ($\nu \ne 0$) should be rejected.

What remains is the continuum of fast magneto-sonic modes ($\nu = 0$) with frequencies along the real axis of the $\omega$-plane. They may be considered as the limit of “wall”-confined modes where the wall position $w \rightarrow \infty$. 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.

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) ($\sigma = 2.6018$, $\nu = 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) ($\sigma = 3.2320$, $\nu = 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.

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 \rightarrow \infty$ as $L_x \rightarrow \infty$. 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 \rightarrow \infty$ as $L_x \rightarrow \infty$ 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 $\omega = 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 $\omega = \pm \infty$ if the radial interval $(1,w)$ is finite (illustrated in figures 9 and 10) and it condenses into a continuum covering the whole real $\omega$-axis in the limit $w \rightarrow \infty$. There are no additional fast discrete modes or quasi-modes elsewhere in the complex $\omega$-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 $\delta$-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 $\hat {q}$ permits. When $k_0 \ne 0$, depending on the values of the tangential wavenumbers, $\hat {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_\perp$ and $k_\parallel$ 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 $\omega$-axis are represented by points in the complex $S$-plane on the ‘elastic scattering circle’, $|S| \equiv |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 (2003, 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 $\omega$-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 $\tilde {\xi }$ (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:

(4.3)\begin{equation} \tilde{\xi}(x;\omega^2) = \lim_{w \rightarrow \infty}\int_{0}^{w} G(x, x';\omega^2) R(x') \,{{\rm d}x}' ,\quad G(x, x'; \omega^2) \equiv \frac{\varGamma(x, x'; \omega^2)}{\Delta(\omega^2)} . \end{equation}

Here, $R$ represents the initial data, whereas the constituent functions $\varGamma$ and $\Delta$ (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 $\Delta$, defined by

(4.4)\begin{equation} \Delta(\omega^2) \equiv b^2(x) \left[ U_1(x; \omega^2) U_2'(x; \omega^2) - U_1'(x; \omega^2) U_2(x; \omega^2) \right], \end{equation}

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_1^{(2)}$ (corresponding to the incoming wave contribution $\hat {\beta }\exp (-\textrm {i}\hat {q}x)$ 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 \rightarrow \infty$ in the upper part of the complex $\omega$-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 $\xi$. This is no accident: the actual fast discrete cluster spectrum (for $w < \infty$) or continuous spectrum (for $w \rightarrow \infty$) is complete, additional real or complex modes cannot be found by solving the IVP.

5. Conclusion

By three modifications of the standard plane slab model for coronal magnetic flux tubes, we have shown the following.

1. (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.

2. (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 $\xi$ and of the perturbation $\varPi$ 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 $\omega$-plane.

3. (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_\theta = 2{\rm \pi} r$ over the angle $\theta$, 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 $\omega$. 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.