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 (Reference Goedbloed and Keppens2022).

The mentioned program makes use of the method of the spectral web, extensively described in previous papers on the subject (Goedbloed Reference Goedbloed2018a,Reference Goedbloedb). 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 (Reference Goedbloed1971a) 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 Reference Goedbloed2018a,Reference Goedbloedb; Goedbloed & Keppens Reference Goedbloed and Keppens2022): 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 Reference Rosenbluth1959) 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 (Reference Newcomb1960) 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 (Reference Goedbloed1971b). 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 (Reference Andries and Goossens2007), 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 (Reference Keppens1994). 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 (Reference Cally2003, Reference Cally2006) and of the initial value problem (IVP) for those modes by Ruderman & Roberts (Reference Ruderman and Roberts2006b,Reference Ruderman and Robertsa)? 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. Reference Goedbloed, Keppens and Poedts2019). (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 (Reference Cally2006) by assigning physical significance to these eigenvalues and in the analysis by Ruderman & Roberts (Reference Ruderman and Roberts2006b) 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 (Reference Ruderman and Roberts2006b) 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 (Reference Goedbloed1998) once more rejecting the conjecture by Grad (Reference Grad1973) 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.