Standing Sausage Perturbations in Solar Coronal Loops with Diffuse Boundaries: An Initial Value Problem Perspective

There have been abundant observational instances of low-frequency waves and oscillations in the structured solar corona (see, e.g., the reviews by Banerjee et al. 2007; De Moortel & Nakariakov 2012; Nakariakov et al. 2016; Wang 2016). Combined with magnetohydrodynamic (MHD) wave theory, these observations can help deduce the atmospheric parameters that prove difficult to directly measure, thereby constituting the field of "coronal seismology" (see, e.g., the reviews by Nakariakov & Verwichte 2005; Nakariakov & Kolotkov 2020; also the textbook by Roberts 2019). Evidently, a physical interpretation needs to be assigned to an observed oscillatory signal for it to be seismologically exploited. For this purpose, it has been customary to contrast observations with the theoretical expectations for waves in field-aligned cylinders that are structured only in the radial direction and in a step form (developed by, e.g., Wentzel 1979; Spruit 1982; Edwin & Roberts 1983, hereafter ER83; also Zajtsev & Stepanov 1975; Cally 1986). It turns out that this apparently simple equilibrium supports a rich variety of waves, and we restrict ourselves to the fast family (see Wang et al. 2021, for the most recent review of the slow family). Indeed, radial fundamental kink modes in the sense of Goossens et al. (2009, 2012) have been amply identified and put to seismological practice (see the review by Nakariakov et al. 2021). As an outcome, the spatial variations of the magnetic field strength were deduced not only for individual active regions (ARs; Anfinogentov & Nakariakov 2019) but also over a substantial fraction of the lower corona (Yang et al. 2020).

Candidate fast sausage modes (FSMs), however, have only been sporadically reported in coronal observations (see Li et al. 2020, for the most recent review). As detailed therein, this rarity is intimately connected to the cutoff axial wavenumbers k_cutoff, to explain which it suffices to consider the pressureless MHD. In fact, we will adopt pressureless MHD throughout and, additionally, restrict ourselves to flare and AR loops as waveguiding inhomogeneities. ¹ Let R denote the cylinder radius and ρ_i (ρ_e ) the internal (external) density with ρ_i > ρ_e . Likewise, let v_Ai (v_Ae) represent the internal (external) Alfvén speed. Standard analysis of the eigenvalue problem (EVP) on a laterally open domain then yields that FSMs in an ER83 equilibrium possess a series of ${k}_{\mathrm{cutoff},m}={j}_{0,m}/(R\sqrt{{\rho }_{i}/{\rho }_{e}-1})$, where j_0,m is the mth zero of Bessel J₀ with m = 1, 2, ⋯ (e.g., Roberts et al. 1984; Vasheghani Farahani et al. 2014). Now consider an open system unbounded in both the axial (z) and radial (r) directions and suppose that a cylinder is perturbed by an axisymmetric perturbation localized both radially and axially. The pertinent two-dimensional initial value problem (2D IVP) has been examined rather extensively, with the majority of solutions found by directly evolving the MHD equations (e.g., Shestov et al. 2015; Yu et al. 2016a, 2017). Physical insights, on the other hand, can also be gleaned from a modal approach, which was discussed heuristically by Edwin & Roberts (1986) and made more formal by Berghmans et al. (1996). Our study makes frequent reference to Oliver et al. (2015, hereafter ORT15), who were the first to offer an explicit expression for the solution to the 2D IVP. Noting that a continuous range of axial wavenumbers (k) is involved, Equation (25) in ORT15 expresses the solution as the summation of the contributions associated with an individual k, which in turn were written as the superposition of eigenmodes with individual angular frequencies (ω, hereafter "frequency" for brevity). A finite number of discrete ω pertaining to proper eigenmodes (or "trapped modes" in physical terms) are relevant only when k > k_cutoff,1, whereby the periodicity is consistently $\lesssim 2\pi /({k}_{\mathrm{cutoff},1}{v}_{\mathrm{Ae}})\approx 2.6\sqrt{1-{\rho }_{e}/{\rho }_{i}}(R/{v}_{\mathrm{Ai}})\lt 2.6R/{v}_{\mathrm{Ai}}$. Regardless of k, however, a continuum of improper eigenmodes is always involved, the associated ω extending from kv_Ae out to infinity. The point is, only proper modes survive in the sausage wave trains sampled sufficiently far from the exciter, the characteristic periodicities therefore being similar to the transverse Alfvén time R/v_Ai. If an individual k is examined, as happens for standing modes, then one finds by directly evolving the MHD equations that R/v_Ai consistently characterizes FSMs regardless of k (e.g., Terradas et al. 2007; Nakariakov et al. 2012; Guo et al. 2016; Lim et al. 2020). While this result is much expected for k > k_cutoff,1, its physical understanding for the opposite situation is a bit involved, given the likely contributions due to improper modes with ω not far exceeding kv_Ae (ORT15; see also our Appendix A). The quick answer is that the interference of the improper modes tends to make their superposition favor a discrete set of periods P_leaky that pertain to the so-called "leaky modes" (Andries & Goossens 2007, and references therein), and P_leaky is well known to be either similar to or substantially shorter than 2π/(k_cutoff,1 v_Ae) (e.g., Meerson et al. 1978; Cally 1986; Kopylova et al. 2007). The damping time of the discrete leaky modes (τ_leaky) is also known to offer a shortcut estimate for the timescale characterizing the wave attenuation, the result being that τ_leaky/P_leaky ≈ (ρ_i /ρ_e )/π² (e.g., Kopylova et al. 2007). Two primary reasons are now clear to account for the rarity of candidate coronal FSMs. First, R/v_Ai typically evaluates to at most a couple of tens of seconds, thereby demanding a high instrumental cadence and ruling out the possibility for typical (extreme) ultraviolet instruments to resolve an FSM (see Su et al. 2012; Tian et al. 2016, for exceptions). Second, there tends to be a stringent requirement on instrumental sensitivity as well. For AR loops, that they are thin and tenuous means that FSMs tend to be detectable only as wave trains, which are indeed compatible with a number of high-cadence ground-based measurements in visible light during total eclipses (e.g., Williams et al. 2002; Katsiyannis et al. 2003; Samanta et al. 2016). For flare loops, that they are thick and dense means that FSMs have a better chance to be detected as both wave trains and standing modes, provided once again that the instrumental cadence is sufficient (see Li et al. 2020, and references therein). This explains why the candidate coronal FSMs reported so far have been primarily connected to radio measurements of short-period quasiperiodic pulsations (QPPs; see the recent reviews by McLaughlin et al. 2018; Zimovets et al. 2021).

With the ER83 equilibrium apparently idealized, one may argue that cutoff wavenumbers are not inherent to coronal FSMs in reality. Indeed, there have been a considerable number of theoretical studies that extend ER83 by incorporating various aspects of reality (see Li et al. 2020, and references therein). Among these, we focus on the equilibria that differ from ER83 only by replacing the step density profile with a continuous one, the reason being that a generic guiding principle can be established to tell when cutoff wavenumbers exist (Lopin & Nagorny 2015a, hereafter LN15; also Lopin & Nagorny 2015b). Let R now refer to some mean cylinder radius, and let the subscript i (e) refer to the equilibrium quantities at the cylinder axis (infinitely far from the cylinder). The radial profile for the equilibrium density ρ₀(r) can then be described in a generic form, ρ₀(r) = ρ_e + (ρ_i − ρ_e )f(r), where the function f(r) evaluates to unity (zero) when r = 0 (r → ∞). Restrict ourselves to the case where f(r) is monotonic. With Kneser's oscillation theorem, LN15 were the first to point out that cutoff wavenumbers exist only when r² f(r) does not diverge when r approaches infinity. This expectation was then verified numerically by Li et al. (2018), one example being for the so-named "outer-μ" profile where f(r) is identically unity for r < R but of the form (r/R)^−μ otherwise. Figure 8 therein indicates that for the mth radial harmonic (m = 1, 2, ⋯), no cutoff wavenumber exists (or, equivalently, k_cutoff,m = 0) when μ < 2, whereas the combination $({k}_{\mathrm{cutoff},m}R)\sqrt{{\rho }_{i}/{\rho }_{e}-1}$ increases monotonically from unity for μ = 2 to j_0,m for a step profile (μ = ∞).

Some important consequences arise for FSMs when cutoff wavenumbers are absent. Theoretically, the dispersive properties of FSMs in this situation are distinct from FSMs in ER83 in two aspects, one being that trapped modes are allowed regardless of the axial wavenumber k or the density contrast ρ_i /ρ_e , and the other being that FSMs tend to be weakly dispersive for small k with axial phase speeds only marginally smaller than v_Ae (e.g., Figure 3 in LN15 and Figure 7 in Yu et al. 2017). With the former distinction evident, we note that trapped FSMs in ER83 are highly dispersive, at least when k is not far larger than a cutoff (e.g., Edwin & Roberts 1983; Roberts et al. 1983). Observationally, these two distinctions offer a richer possibility for interpreting oscillatory signals, to illustrate which we consider a spatially resolved QPP measured with the Nobeyama RadioHeliograph as reported by Kupriyanova et al. (2013). As detailed therein, multiple periodicities were simultaneously found, with the associated spatial distributions of the spectral power strongly indicating an axial fundamental together with its overtones in the involved flare loop. Contrasting the observations with the canonical ER83 theory, the authors deduced that these standing modes belong to the kink family, and FSMs were ruled out because of their dispersive properties. However, adopting density profiles similar to the outer-μ one with μ < 2, both LN15 and Lopin & Nagorny (2019) suggested that the observations may be interpreted as FSMs as well. Put to seismology, this interpretation returned values for the internal Alfvén speed v_Ai that may differ considerably from those returned with the interpretation in terms of kink modes, for which purpose we quote ∼1100 km s⁻¹ from Lopin & Nagorny (2019, Section 6) and ∼1750 km s⁻¹ from Kupriyanova et al. (2013, Section 5.1). Strictly speaking, a comparison between the two sets of v_Ai is not straightforward because Lopin & Nagorny (2019) adopted pressureless MHD, whereas a finite gas pressure is considered in Kupriyanova et al. (2013). ² Our point is that the disappearance of cutoff wavenumbers for FSMs as a result of some straightforward departure of the equilibrium from ER83 can offer more physical interpretations for observations in the first place and enable more seismological possibilities afterward.

Some further consequence arises if we now focus on standing FSMs in AR loops. For ease of description, let us recall that we consistently work in pressureless MHD and adopt the customary assumption that sees AR loops as straight, density-enhanced, field-aligned cylinders. We further assume that the radial density distribution is of the outer-μ type, which is reasonable but admittedly difficult to prove or disprove (e.g., Aschwanden et al. 2003; Goddard et al. 2017). In addition, we assume that lower coronal eruptions (LCEs), the primary exciter for the much-observed large-amplitude radial fundamental kink modes (Zimovets & Nakariakov 2015; Nechaeva et al. 2019), can deposit a nonnegligible amount of energy as axisymmetric perturbations to AR loops as well. Note that this assumption is not that bold but has been implied in the interpretation of rapidly propagating waves as sausage wave trains (see the review by Roberts 2008, and references therein). For our purposes, it suffices to consider only axial fundamentals. With dimensionless cutoff wavenumbers ${k}_{\mathrm{cutoff}}R\geqslant 1/\sqrt{{\rho }_{i}/{\rho }_{e}-1}$ when μ ≥ 2, one finds k_cutoff R ≥ 1/3, given the typical range of [2, 10] quoted for the density contrast ρ_i /ρ_e (Aschwanden et al. 2004). One further finds kR = π R/L ≲ π/15 in view of the measurements of widths and lengths for AR loops imaged in the extreme ultraviolet (EUV; Schrijver 2007, Figure 1). Standing FSMs, at least the axial fundamentals, are therefore unlikely to be observable for two reasons. One is that their periodicities will be ∼R/v_Ai and therefore short, and the other is that they tend to experience rapid attenuation as well. Let us stress that these two signatures have not been explicitly shown for this particular outer-μ profile but are expected with the studies on 1D IVPs addressing standing FSMs in the leaky regime for an ER83 equilibrium (e.g., Terradas et al. 2007; Nakariakov et al. 2012). Now consider those AR loops with μ < 2. Given the absence of cutoff wavenumbers, the system is expected to settle into a trapped mode or some combination of trapped modes, the quality of the oscillatory signals therefore being sufficiently high. Likewise, the periodicities will eventually be characterized by the longitudinal Alfvén time L/v_Ae, which can be readily resolved with, say, the majority of available UV/EUV instruments. In fact, both expectations have already been invoked in seismological applications, albeit in the context of flare loops (LN15; Lopin & Nagorny 2019). Supposing that AR loops with μ < 2 are not uncommon, one further deduces that a substantial fraction of kink oscillations will be mixed with standing FSMs when LCEs occur. As advocated by Chen et al. (2015) and Guo et al. (2016), the simultaneous observations of multiple modes of distinct nature will then considerably mitigate the nonuniqueness issue inherent to coronal seismology (see Arregui & Goossens 2019, for dedicated remarks). ³ However, standing FSMs have not been reported or even implicated in observations of oscillating AR loops to our knowledge. An obvious excuse is that observers have nearly exclusively adopted the ER83 framework and therefore dismissed the possibility that FSMs may possess periodicities ≳2L/v_Ae altogether. ⁴ Our point, however, is that this possibility is expected solely on the basis of EVP analyses on an open domain, and one has yet to demonstrate that FSMs with periodicities ≳2L/v_Ae do exist as solutions to the pertinent IVP.

Zoom In Zoom Out Reset image size

Focusing on sausage oscillations in AR loops with the outer-μ family of density profiles, we intend to address the question "does the absence of cutoff wavenumbers guarantee a temporal behavior that is distinct from the situation where cutoff wavenumbers are present?" We decompose this question into two interconnected aspects. One, how does the value of μ influence the timescale that characterizes the energy attenuation? Two, does the transverse or longitudinal Alfvén time characterize the periodicity when a wave signal is sufficiently strong? This manuscript is structured as follows. Section 2 formulates the IVP for a radially open system, which is then solved with a direct finite-difference (FD) approach in Section 3. While the answer to our question is already clear in the FD solutions, Section 4 moves on to solve the IVP by superposing eigensolutions to the relevant EVP on a closed domain. These modal solutions are presented for more than just cross-validation purposes. Rather, they help quantify the specific contributions from individual frequencies. By experimenting with various domain sizes, we will better connect the solutions to our IVP with the theoretical expectations from the analyses of EVPs on an open domain. Section 5 summarizes this study, ending with some concluding remarks.

We adopt pressureless ideal MHD as our theoretical framework, in which the primitive variables are the mass density ρ, velocity v , and magnetic field B . The equilibrium quantities are denoted with a subscript 0, and the equilibrium is taken to be static ( v ₀ = 0). Working in a cylindrical coordinate system (r, θ, z), we take the equilibrium magnetic field to be uniform and directed in the z-direction ( B ₀ = B₀ e _z ). Seeing AR loops as density-enhanced cylinders with some mean radius R, we assume that the equilibrium density (ρ₀) depends only on r and decreases from ρ_i at the cylinder axis (r = 0) to ρ_e infinitely far from the cylinder (r → ∞). The Alfvén speed is defined by ${v}_{{\rm{A}}}^{2}={B}_{0}^{2}/({\mu }_{0}{\rho }_{0})$, with μ₀ being the magnetic permeability of free space. From here onward, by "internal" (subscript i) and "external" (subscript e), we consistently refer to the equilibrium quantities evaluated at r = 0 and r → ∞, respectively. The internal (external) Alfvén speed is therefore denoted by v_Ai (v_Ae).

2.1. Preliminary Formulation of the IVP

We now formulate the preliminary version of the IVP in a radially open system. Let the subscript 1 denote small-amplitude perturbations to the equilibrium. Specializing to sausage perturbations (∂/∂θ = 0), the linearized, pressureless, ideal MHD equations read

$$\begin{eqnarray}&&{\rho }_{0}\displaystyle \frac{\partial {v}_{1r}}{\partial t}=\displaystyle \frac{{B}_{0}}{{\mu }_{0}}\left(\displaystyle \frac{\partial {B}_{1r}}{\partial z}-\displaystyle \frac{\partial {B}_{1z}}{\partial r}\right),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {B}_{1r}}{\partial t}={B}_{0}\displaystyle \frac{\partial {v}_{1r}}{\partial z},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {B}_{1z}}{\partial t}=-{B}_{0}\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}\left({{rv}}_{1r}\right)\,.\end{eqnarray} \tag{ 3 }$$

With coronal cylinders bounded by the planes z = 0 and L in mind, the following ansatz

$$\begin{eqnarray}\begin{array}{rcl}{v}_{1r}(r,z;t) & = & \hat{v}(r;t)\sin ({kz}),\\ {B}_{1r}(r,z;t) & = & {\hat{B}}_{r}(r;t)\cos ({kz}),\\ {B}_{1z}(r,z;t) & = & {\hat{B}}_{z}(r;t)\sin ({kz}),\end{array}\end{eqnarray} \tag{ 4 }$$

is appropriate for axial standing modes, with k = n π/L being the quantized axial wavenumber (n = 1, 2, ⋯ ). Equations (1)–(3) then become

$$\begin{eqnarray}&&{\rho }_{0}\displaystyle \frac{\partial \hat{v}}{\partial t}=-\displaystyle \frac{{B}_{0}}{{\mu }_{0}}\left(k{\hat{B}}_{r}+\displaystyle \frac{\partial {\hat{B}}_{z}}{\partial r}\right),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\hat{B}}_{r}}{\partial t}={{kB}}_{0}\hat{v},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\hat{B}}_{z}}{\partial t}=-{B}_{0}\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}\left(r\hat{v}\right)\,.\end{eqnarray} \tag{ 7 }$$

Without a loss of generality, the initial conditions (ICs) are specified as

$$\begin{eqnarray}&&\hat{v}(r,t=0)=u(r),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\hat{B}}_{r}(r,t=0)={\hat{B}}_{z}(r,t=0)=0.\end{eqnarray} \tag{ 9 }$$

The boundary condition (BC) at the cylinder axis (r = 0) reads

$$\begin{eqnarray}&&\hat{v}={\hat{B}}_{r}=\partial {\hat{B}}_{z}/\partial r=0,\end{eqnarray} \tag{ 10 }$$

whereas the BC at r → ∞ is irrelevant.

It proves necessary to examine the energetics associated with the IVP as well. Let V refer to a volume bounded laterally by a cylindrical surface with radius r and horizontally by the planes z = 0 and L. One then finds from Equations (5)–(7) that

$$\begin{eqnarray}&&{E}_{\mathrm{tot}}(r,t)-{E}_{\mathrm{tot}}(r,t=0)=-F(r,t),\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{E}_{\mathrm{tot}}(r,t) & = & \pi L{\displaystyle \int }_{0}^{r}(r^{\prime} {dr}^{\prime} )\left\{\displaystyle \frac{1}{2}{\rho }_{0}(r^{\prime} ){\hat{v}}^{2}(r^{\prime} ,t)\right.\\ & & \left.+\,\displaystyle \frac{1}{2{\mu }_{0}}\left[{\hat{B}}_{r}^{2}(r^{\prime} ,t)+{\hat{B}}_{z}^{2}(r^{\prime} ,t)\right]\right\}\,,\end{array}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&F(r,t)=\pi L{\int }_{0}^{t}{dt}^{\prime} \left[r{\hat{p}}_{\mathrm{tot}}(r,t^{\prime} )\hat{v}(r,t^{\prime} )\right]\,.\end{eqnarray} \tag{ 13 }$$

Here a common factor π L is retained to ensure that E_tot(r, t) represents the instantaneous total energy in V, while F(r, t) represents the cumulative energy loss from V. Furthermore, ${\hat{p}}_{\mathrm{tot}}={B}_{0}{\hat{B}}_{z}/{\mu }_{0}$ is connected to the Eulerian perturbation of total pressure. Evidently, the terms in the brackets in Equation (13) stem from the radial component of the Poynting vector.

2.2. Reformulation of the IVP and Parameter Specification

For mathematical convenience, Equations (5)–(10) are reformulated to the following form.

IVP 1. Solutions are sought for the following equation:

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}\hat{v}}{\partial {t}^{2}}={v}_{{\rm{A}}}^{2}(r)\left(\displaystyle \frac{{\partial }^{2}\hat{v}}{\partial {r}^{2}}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial \hat{v}}{\partial r}-\displaystyle \frac{\hat{v}}{{r}^{2}}-{k}^{2}\hat{v}\right),\end{eqnarray} \tag{ 14 }$$

subjected to the ICs

$$\begin{eqnarray}&&\hat{v}(r,t=0)=u(r),\quad \displaystyle \frac{\partial \hat{v}}{\partial t}(r,t=0)=0,\end{eqnarray} \tag{ 15 }$$

together with the BC

$$\begin{eqnarray}&&\hat{v}(r=0,t)=0,\end{eqnarray} \tag{ 16 }$$

on a domain spanning from r = 0 to ∞.

Necessary for energetics considerations, ${\hat{B}}_{r}(r,t)$ and ${\hat{B}}_{z}(r,t)$ can be found with $\hat{v}(r,t)$ by integrating Equations (6) and (7) from the initial state (Equation (9)).

We proceed to make IVP 1 more specific. The equilibrium density distribution is chosen to be the outer-μ profile in Yu et al. (2017), namely,

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{0}(r) & = & {\rho }_{e}+({\rho }_{i}-{\rho }_{e})f(r),\\ f(r) & = & \left\{\begin{array}{ll}1, & 0\leqslant r\leqslant R,\\ {\left(r/R\right)}^{-\mu }, & r\geqslant R.\end{array}\right.\end{array}\end{eqnarray} \tag{ 17 }$$

Here μ ≥ 1 measures the steepness of ρ₀(r) outside the cylinder. We focus on axial fundamentals (k = π/L). In addition, we specify the initial perturbation in Equation (15) as

$$\begin{eqnarray}\displaystyle \frac{u(r)}{{v}_{\mathrm{Ai}}}=\left\{\begin{array}{ll}{\sin }^{3}(\pi r/{\rm{\Lambda }}), & 0\leqslant r\leqslant {\rm{\Lambda }},\\ 0, & \quad r\geqslant {\rm{\Lambda }},\end{array}\right.\end{eqnarray} \tag{ 18 }$$

which is localized within r = Λ and prescribed to be sufficiently smooth with a magnitude arbitrarily set to be the internal Alfvén speed (v_Ai).

The solution to IVP 1 is fully determined by the dimensionless parameters [ρ_i /ρ_e , μ; L/R; Λ/R], among which we see μ as the primary adjustable one. The density contrast and loop length-to-radius ratio are fixed at [ρ_i /ρ_e , L/R] = [2.25, 15], both close to the lower end of but nonetheless within the accepted range for AR loops (e.g., Aschwanden et al. 2004; Schrijver 2007). We take Λ = 4 R unless otherwise specified. Figure 1(a) illustrates our equilibrium, and the blue arrows represent the initial velocity field in any cut through the cylinder axis as appropriate for an axial fundamental. Specializing to IVP 1, Figure 1(b) shows the radial profiles for $u(r)=\hat{v}(r,t=0)$ (blue dashed curve) and the equilibrium density (ρ₀; solid curves). Two values are adopted for μ, one being 1.5 (black curve) and the other being 5 (red curve). As already stressed, FSMs do not suffer from cutoff wavenumbers k_cutoff for μ < 2. When μ ≥ 2, k_cutoff R always exceeds $1/\sqrt{{\rho }_{i}/{\rho }_{e}-1}$ and therefore ≳0.89 with the chosen ρ_i /ρ_e , making trapped modes irrelevant for the chosen kR = π/15 ≈ 0.21.

We choose to solve IVP 1 with two independent methods, one being an FD approach and the other being a modal approach involving eigenmodes for the EVPs on either a closed or an open domain. In practice, the FD approach turns out to be orders of magnitude less time-consuming and is therefore more suitable for parametric studies.

3.1. Method

3.2. Numerical Results

Figure 2 presents our FD solutions by showing the distribution of the radial speed $\hat{v}$ in the r − t plane. We contrast two cases with the steepness parameter (a) μ = 1.5 and (b) μ = 5. One sees from Figure 2 that the temporal evolution in both cases is characterized by some dispersive propagation of the axisymmetric disturbance. To proceed, we note that the disturbances belong to the fast family, and the local fast speed is equivalent to the local Alfvén speed in pressureless MHD. For ease of description, the most prominent wavefronts are singled out and labeled such that the plus (minus) sign pertains to outward (inward) propagation. When necessary, the same label will be used to denote the associated wake as well. The following features can then be readily told in both cases, to describe which it suffices to consider only Figure 2(a). First, the initial perturbation splits into two wave fronts manifested as the two bright stripes labeled 1₊ and 1₋. While not that evident, wave front 1₋ is actually accompanied by a wake appearing as a narrow dark stripe. At least a substantial fraction of both wave front 1₋ and its wake then make it into the cylinder (r < R), as evidenced by the change of slope of the stripe. Here, by "a substantial fraction," we mean that some reflection is expected but difficult to identify. Second, once reaching the cylinder axis (r = 0), wave front 1₋ is reflected to form wave front 2₊, which then propagates outward as a dark stripe. The change from a bright to a dark stripe for essentially the same wave front is simply because $\hat{v}$ necessarily reverses sign at the cylinder axis, which acts as a rigid wall in the present context (see also Berghmans et al. 1996, hereafter BDBG96). In this sense, the bright stripe following wave front 2₊, namely, wake 2₊, is actually the reflected wake 1₋. Third, the partial reflection of wake 2₊ around the cylinder boundary (r = R) then leads to wave front 3₋, part of which is guided by the dashed curve. Wave front 4₊ then results from the reflection of wave front 3₋ at the cylinder axis. Having described these common features, we note that some differences nonetheless exist between the two cases. For instance, wave front 3₋ is easier to tell for μ = 5 than for μ = 1.5. This is understandable because the case with μ = 5 corresponds to a steeper v_A profile around the cylinder boundary and hence a stronger partial reflection there.

Zoom In Zoom Out Reset image size

The slight differences notwithstanding, in both cases, one expects a continuous decrease for the wave energy in the cylindrical volume V bounded by r = Λ, namely, where the initial perturbation is applied. This expectation is indeed true, ⁵ and to demonstrate it, we display the temporal variations of the total energy in V (E_tot; dashed curves) and the cumulative energy loss from V (F; dashed–dotted curves) in Figure 3. The sum E_tot + F is further given by the solid curves. We discriminate between the cases with μ = 1.5 and 5 with different colors. From the solid curves, one sees that energy conservation is maintained remarkably well, to quantify which we quote an accuracy of better than 0.05% for all FD computations. Examining the dashed curves, one sees that E_tot shows a couple of plateaus, with the behavior of E_tot for μ = 5 around t ∼ 5 R/v_Ai being an example. In view of Equations (11) and (13), these plateaus appear simply because the radial component of the Poynting flux tends to vanish therein. More importantly, E_tot rapidly decreases with time, with the two most prominent intervals readily accounted for by the passage of wave fronts 1₊ and 2₊ (see Figure 2). Virtually no energy is left in V when t ≳ 7 R/v_Ai, which is true for μ = 5 and 1.5 alike.

Zoom In Zoom Out Reset image size

The rapid attenuation of wave energy can be told, in a more straightforward way, by directly showing the temporal evolution of the radial speed itself $\hat{v}$. This is done in Figure 4, where we plot $\hat{v}$ at r = R and use different colors to discriminate different values of μ. We note that the FD solutions are shown by the solid curves, labeled "FD open" to reflect that they are found on an open domain. The asterisks, labeled "modal closed," represent the solutions found by superposing the eigenfunctions for the EVP on a closed domain. The details of the modal solutions are not important for now. What matters is that they agree with the FD ones exactly, thereby suggesting the reliability of both sets of solutions. Consider now only the solid curves. With the aid of Figure 2, one readily identifies the first three extrema with wave fronts 1₋ and 2₊, as well as wake 2₊. One further sees one more extremum (two more extrema) in the black (red) curve, the corresponding wave fronts/wakes also identifiable in Figure 2. When discernible, any extremum in the black curve appears later than its counterpart in the red curve. The explanation for this behavior is actually straightforward because the extrema in $\hat{v}(R,t)$ ultimately derive from wave front 1₋. One sees from Figure 1 that the local fast speed v_A(r) at any r > R is larger for μ = 5 than for μ = 1.5. It therefore takes more time for wave front 1₋ to enter the cylinder (r ≤ R) in the case with μ = 1.5, thereby making the relevant extrema in $\hat{v}(R,t)$ appear later. With this understanding, the spacing between two consecutive prominent extrema or, equivalently, the periodicity when the signal is strong is intimately connected to the passage of the relevant wave fronts/wakes before they appear in $\hat{v}(R,t)$. The end result is that the periodicity is expected to depend on the details of both the equilibrium density profile and the initial perturbation. In our setup, this translates into the dependence on μ and Λ.

Zoom In Zoom Out Reset image size

Figure 5 quantifies the dependence of the wave behavior on the steepness parameter (μ) for a number of values of the spatial extent of the initial perturbation (Λ) as labeled. Two quantities are examined, one being the time that it takes for E_tot(Λ, t) to drop from the initial value by a factor of e⁴ ≈ 55 (τ_ener; Figure 5(a)) and the other being the temporal spacing between the first two extrema in the $\hat{v}(R,t)$ profile (Δ₁; Figure 5(b)). Let us examine Figure 5(a) first and start by noting that we deliberately choose a rather large factor (e⁴) to determine τ_ener. We note further that the initial perturbation peaks at r = Λ/2 (see Equation (18)). Now two prominent features are evident. First, τ_ener at any given μ increases with Λ. Second, the μ-dependence of τ_ener tends to be weak when Λ ≲ 2 R. Let V still denote the cylindrical volume laterally bounded by r = Λ. It turns out that the departure of wave fronts 1₊ and 2₊ from V is the primary reason for E_tot(Λ, t) to decrease to the designated threshold (see Figure 2). In particular, τ_ener tends not to be much longer than the time t₂₊ at which wave front 2₊ arrives at r = Λ. In turn, this transit time t₂₊ comprises two components,

$$\begin{eqnarray}&&{t}_{2+}={t}_{1-}({\rm{\Lambda }}/2\to 0)+{t}_{2+}(0\to {\rm{\Lambda }}),\end{eqnarray} \tag{ 19 }$$

where the symbols on the right-hand side (RHS) are such that t₁₋(Λ/2 → 0) represents the time that wave front 1₋ spends when traveling from r = Λ/2 to r = 0. For a given μ, both terms on the RHS increase with Λ, meaning that t₂₊ and hence τ_ener increase monotonically with Λ. Now move on to the μ-dependence for a given Λ. When Λ ≲ 2 R, the cylinder exterior (r > R) is relevant for determining t₂₊ only by being partially involved in t₂₊(0 → Λ). When Λ > 2 R, however, it is involved in both terms on the RHS of Equation (19). The end result is that t₂₊ and hence τ_ener are insensitive to μ when Λ ≲ 2 R but tend to decrease with μ when the opposite is true, which is understandable, given that the local fast speed v_A(r) at any r > R increases with μ.

Zoom In Zoom Out Reset image size

Now move on to Figure 5(b). One sees that Δ₁ possesses a considerably more complicated behavior, by which we mean the features difficult to understand with the simple kinematic considerations that were applied to Figure 5(a). Take the cases where Λ = 4 R and 8 R. In both cases, the first and second extrema in the $\hat{v}(R,t)$ profile correspond to wave fronts 1₋ and 2₊, respectively. In kinematic terms, the temporal spacing between the two then comprises t₁₋(R → 0) and t₂₊(0 → R), neither of which is supposed to involve Λ or μ. Consequently, the blue and maroon curves are expected to overlap, an expectation evidently at variance with the numerical results. One therefore deduces that Δ₁ embeds some subtleties that the kinematic arguments cannot address, which will become evident in the modal solutions to be presented shortly. The quick answer is the common sense that μ is relevant for determining the eigenstructures, while Λ determines how the energy in the initial perturbation is distributed among the eigenmodes (see Equation (23)). Important for now is that Figure 5 has already answered the questions we laid out in the Introduction. First, Figure 5(a) indicates that the energy imparted by the initial perturbation is attenuated within a characteristic timescale ${\tau }_{\mathrm{ener}}\sim { \mathcal O }({\rm{\Lambda }}/{v}_{\mathrm{Ai}})$. If Λ is not far different from R, then axisymmetric perturbations will rapidly become too weak to detect. Second, even if some instrument happens to capture a perturbation immediately after its excitation, Figure 5(b) indicates that Δ₁ is consistently $\sim { \mathcal O }(R/{v}_{\mathrm{Ai}})$, thereby placing rather stringent demands on the instrumental cadence. Third, for any examined Λ, no abrupt change is seen in the behavior of τ_ener or Δ₁ when μ crosses the nominally critical value of 2. Put another way, trapped modes are not discernible even though they suddenly appear when μ drops below 2 in EVPs on an open domain. We address why in what follows.

4.1. Method

Our modal approach starts with specifying the following EVP on a closed domain.

EVP 1. Nontrivial solutions are sought for the following equation:

$$\begin{eqnarray}&&-{\omega }^{2}\breve{v}={v}_{{\rm{A}}}^{2}(r)\left(\displaystyle \frac{{d}^{2}}{{{dr}}^{2}}\breve{v}+\displaystyle \frac{d}{{rdr}}\breve{v}-\displaystyle \frac{\breve{v}}{{r}^{2}}-{k}^{2}\breve{v}\right),\end{eqnarray} \tag{ 20 }$$

defined on a domain of [0, d] and subjected to the BCs

$$\begin{eqnarray}&&\breve{v}(r=0)=\breve{v}(r=d)=0.\end{eqnarray} \tag{ 21 }$$

Equation (20) is found by replacing $\hat{v}$ with $\mathrm{Re}[\breve{v}(r)\exp (-i\omega t)]$ in Equation (14).

It is known that EVP 1 possesses the following Sturm–Liouville properties (see BDBG96 for details, even though a step profile is examined therein). First of all, the eigenvalues (ω²) are positive and form an infinite, discrete, monotonically increasing sequence {ω_l } with respect to the mode number l = 1, 2, ⋯ . The associated eigenfunction ${\breve{v}}_{l}(r)$ can be made and will be seen as real-valued. It then follows that ${\breve{v}}_{l}(r)$ possesses l − 1 nodes inside the domain. In addition, the set $\{{\breve{v}}_{l}(r)\}$ is complete and satisfies the orthogonality condition

$$\begin{eqnarray}&&{\int }_{0}^{d}{\breve{v}}_{l}(r){\breve{v}}_{m}(r){\rho }_{0}(r){rdr}=0,\end{eqnarray} \tag{ 22 }$$

provided l ≠ m. Eventually, the solution to IVP 1 can be written as

$$\begin{eqnarray}&&\begin{array}{l}{\hat{v}}^{(d)}(r,t)=\displaystyle \sum _{l=1}^{\infty }{c}_{l}{\breve{v}}_{l}(r)\cos ({\omega }_{l}t),\\ 0\leqslant r\leqslant d,\quad t\leqslant {\displaystyle \int }_{{\rm{\Lambda }}}^{d}\displaystyle \frac{{dr}}{{v}_{{\rm{A}}}(r)},\end{array}\end{eqnarray} \tag{ 23 }$$

where the coefficient c_l measures the contribution from the lth mode,

$$\begin{eqnarray}&&{c}_{l}=\displaystyle \frac{{\int }_{0}^{d}u(r){\breve{v}}_{l}(r){\rho }_{0}(r){rdr}}{{\int }_{0}^{d}{\breve{v}}_{l}^{2}(r){\rho }_{0}(r){rdr}}.\end{eqnarray} \tag{ 24 }$$

The superscript (d) in Equation (23) is meant to indicate that the modal structure depends on the domain size d. Here and hereafter, by "modal structure," we further mean the l-dependence of ω_l . Expressed formally, ω_l can be written as

$$\begin{eqnarray}&&\displaystyle \frac{{\omega }_{l}R}{{v}_{\mathrm{Ai}}}={{ \mathcal F }}_{l}({\rho }_{i}/{\rho }_{e},\mu ;{kR};d/R).\end{eqnarray} \tag{ 25 }$$

However, we stress that the modal solution ${\hat{v}}^{(d)}(r,t)$ itself does not depend on d in the time frame of validity explicitly given in Equation (23). Physically, this time frame of validity simply represents the interval when the outermost edge of the perturbation is within the domain of EVP 1.

We now describe some details involved in the evaluation of Equation (23). To start, [ρ_i /ρ_e , kR, Λ/R] is fixed at [2.25, π/15, 4] for all modal solutions. In the main text, we contrast only two steepness parameters (μ = 1.5 and 5) but experiment with a substantial number of dimensionless domain sizes (d/R). The step profile (μ = ∞) is also of interest, the discussions nonetheless collected in Appendix A. Regardless, we consistently formulate and solve EVP 1 with the general-purpose finite-element code PDE2D (Sewell 1988), which was first introduced into the solar context by Terradas et al. (2006). A uniform grid is adopted if d/R ≤ 150; otherwise, we employ a grid whereby the spacing is uniform for r ≤ 5 R but increases by a constant factor afterward. We make sure that different grid setups yield consistent solutions to EVP 1. Likewise, we make sure that the number of modes (${l}_{\max }$) incorporated in the summation in Equation (23) is sufficiently large, meaning that increasing ${l}_{\max }$ does not influence the modal solutions to be analyzed.

Some insights into EVP 1 are further necessary to address the roles of μ and d/R. These are made transparent when Equation (20) is transformed into a Schrödinger form with standard techniques employed by, e.g., BDBG96 and LN15, the result being

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{d}^{2}{\rm{\Phi }}}{{{dr}}^{2}}+Q(r){\rm{\Phi }}=0,\\ Q(r)=\displaystyle \frac{{\omega }^{2}-V(r)}{{v}_{{\rm{A}}}^{2}(r)}.\end{array}\end{eqnarray} \tag{ 26 }$$

Here ${\rm{\Phi }}(r)={r}^{1/2}\breve{v}(r)$ defines some "wave function," while

$$\begin{eqnarray}&&V(r)={v}_{{\rm{A}}}^{2}(r)\left({k}^{2}+\displaystyle \frac{3}{4{r}^{2}}\right)\end{eqnarray} \tag{ 27 }$$

defines the potential. Three properties then ensue. One, ω_l , is bound to exceed kv_Ai regardless of μ or d/R (see Appendix B). Two, high-frequency modes with ω_l ≫ kv_Ae are permitted regardless of μ or d/R, and they follow the relation ω_l ≈ l π v_Ae/d when d/R is sufficiently large (also see Appendix B). Three, the spatial behavior of mode functions $\breve{v}(r)$ is determined by the sign of Q(r). We therefore classify the modes into two categories, labeling those with Q(d) > 0 (Q(d) < 0) as "oscillatory" ("evanescent"). For ease of description, we see this classification scheme as applicable only to closed systems but adopt the viewpoint that such terms as "trapped" versus "leaky" (e.g., ER83; Cally 1986) or "proper" versus "improper" (ORT15; also Oliver et al. 2014) apply to open systems. Evidently, an evanescent mode becomes a trapped mode when d goes infinite. Its frequency ω_l is therefore d/R-independent for sufficiently large d/R and bound to be lower than the "critical frequency" ω_crit = kv_Ae because $V(r)\to {k}^{2}{v}_{\mathrm{Ae}}^{2}$ when r → ∞. Now specialize to our chosen set, [ρ_i /ρ_e , kR] = [2.25, π/15]. Evanescent modes are possible only when μ < 2. All modes are oscillatory for μ ≥ 2, the frequency ω_l for any l always higher than ω_crit and d/R-dependent.

4.2. Numerical Results

4.2.1. Frequency Distribution of Modal Contributions

We fix the domain size to be d = 50 R in this subsection. As shown by Figure 4, the modal solutions thus constructed agree with their FD counterparts. The dependencies on the mode frequency (ω_l ) of the contributions of individual modes are then shown in Figure 6, where the modes are represented by asterisks, and different colors are adopted to discriminate between the two steepness parameters. The critical frequency ω_crit = kv_Ae is plotted by the vertical dashed–dotted lines for reference. By examining ∣c_l ∣, Figure 6(a) overviews the gross contributions from individual modes, where by "gross" we mean that c_l is position-independent (Equation (24)). For μ = 1.5 (black) and 5 (red) alike, the ω_l -dependence of ∣c_l ∣ features a number of peaks whose magnitude weakens with ω_l . The modal contributions at the specific location r = R are plotted in Figure 6(b), from which one sees that the spatial dependence of ${\breve{v}}_{l}(r)$ makes the contributions from modes with ω_l ≳ 3.5 v_Ai/R less pronounced than expected with Figure 6(a). Regardless, the point is that ${\hat{v}}^{(d)}(R,t)$ is dominated by modes with ω_l higher but not far higher than v_Ai/R, in agreement with the periodicities found in Figure 4. More importantly, the lowest mode frequency ω₁ exceeds ω_crit, meaning that all modes are oscillatory. This is true not only for μ = 5 but also for μ = 1.5, despite the fact that trapped modes are bound to appear in the latter case for a truly open system.

Zoom In Zoom Out Reset image size

4.2.2. Dependence of Modal Structure on Domain Size

This subsection examines how the modal structure depends on the dimensionless domain size d/R; the reasons for doing this are twofold. First, trapped modes appear for an open system when μ < 2 was found on solid mathematical grounds by LN15. One naturally argues that d/R in Figure 6 is not large enough for an evanescent mode to appear. However, the modal solution (Equation (23)) does not depend on d/R within the time frame explicitly given there. Figure 6 then indicates that the contribution from evanescent modes is negligible, even if d/R is larger still. We will quantify how negligible this contribution is. Second, to our knowledge, the only study involving EVP 1 was conducted by BDBG96 for a step profile. However, the d/R-dependence was not of interest and hence only briefly mentioned in Figure 2 therein. We analytically examine this dependence in some detail for step and continuous profiles in Appendices A and B, respectively. The analytical results in turn help better quantify the modal behavior that we find numerically and present in the main text.

We start by examining the d/R-dependence of the modal contributions to ${\hat{v}}^{(d)}(r,t)$. For this purpose, we rewrite the modal solution (Equation (23)) as

$$\begin{eqnarray}&&{\hat{v}}^{(d)}(r,t)=\sum _{{\omega }_{l}\lt {{kv}}_{\mathrm{Ae}}}{c}_{l}{\breve{v}}_{l}(r)\cos ({\omega }_{l}t)+\sum _{{\omega }_{l}\gt {{kv}}_{\mathrm{Ae}}}{S}_{l}{\breve{v}}_{l}(r)\cos ({\omega }_{l}t){\rm{\Delta }}{\omega }_{l},\end{eqnarray} \tag{ 28 }$$

where the second summation accommodates all oscillatory modes, whereas the first collects evanescent ones. The frequency spacing

$$\begin{eqnarray}&&{\rm{\Delta }}{\omega }_{l}={\omega }_{l+1}-{\omega }_{l}\end{eqnarray} \tag{ 29 }$$

is relevant only for oscillatory modes in Equation (28) but actually defined for all l. We further view the combination ${S}_{l}{\breve{v}}_{l}(r)$ as some local "spectral density," with S_l defined by

$$\begin{eqnarray}&&{S}_{l}=\displaystyle \frac{{c}_{l}}{{\omega }_{l+1}-{\omega }_{l}}.\end{eqnarray} \tag{ 30 }$$

Why is the modal solution decomposed in such a simple but cumbersome way? We choose to leave a detailed answer until later and for now stress only that it does not make sense to compare the frequency dependencies of ${c}_{l}{\breve{v}}_{l}$ at a specific location between different values of d/R because of the d/R-dependence of ω_l . This is made more specific by Figure 7, where we specialize to r = R and display the ω_l -dependencies of the local spectral densities ${S}_{l}{\breve{v}}_{l}(R)$ for both μ = 1.5 (black symbols) and μ = 5 (red symbols). Two domain sizes are examined, with the result for d/R = 50 (100) represented by asterisks (plus signs). For both values of μ, one sees that the spectral densities for the two values of d/R outline exactly the same curve, even though the mode frequencies are more closely spaced for d/R = 100 than for d/R = 50. One further sees that all modes remain oscillatory for μ = 1.5, even when d/R = 100. In fact, an evanescent mode appears for this μ only when d/R ≳ 4000, a remarkably large value that makes it numerically formidable to compute the necessary set of modes to yield a further ω_l -dependence of ${S}_{l}{\breve{v}}_{l}(R)$. We rather arbitrarily choose d/R = 12,800 and compute only two small subsets for each μ, the first (second) comprising those five modes with ω_l just exceeding 1.5 v_Ai/R (3 v_Ai/R). The corresponding spectral densities are plotted by diamonds. For each μ, one then expects to see 10 but actually can discern only two diamonds because the mode frequencies in each subset are too close to tell apart. Regardless, the diamonds lie exactly on the curve outlined by the result for any smaller d/R, reinforcing the insignificance of the evanescent modes despite that one such mode does exist. To be precise, ${c}_{1}{\breve{v}}_{1}(R)$ for μ = 1.5 evaluates to 3.51 × 10⁻¹¹ in units of v_Ai, or rather in units of the magnitude of the initial perturbation.

Zoom In Zoom Out Reset image size

We focus on how the modal structure varies when the domain size varies. For a sequence of d/R as labeled, Figure 8(a) shows the frequencies (ω_l ) of the first 20 modes as horizontal ticks stacked vertically at a given d/R, with the results for μ = 1.5 and 5 differentiated by black and red. Note that in this "level scheme," ω_l is measured in units of the critical frequency (ω_crit), and the horizontal dashed–dotted line marks ω_l = ω_crit. Two features then follow. By examining the case with d/R = 50 for either value of μ, one sees that the frequency spacing Δω_l becomes increasingly uniform with increasing l. As detailed in Appendices A.2 and B, this feature derives from the fact that for any μ at a sufficiently large d/R, the mode frequency ω_l for a large enough l can be approximated by l π v_Ae/d (Equation (A17)). Slightly surprising is that this approximation is accurate to better than 5% for both μ when l merely exceeds 11, despite the fact that l is nominally required to be ≫(d/R)/π ≈ 15.9 for Equation (A17) to hold. While not shown, we find that ω_l may be approximated by l π v_Ae/d for l beyond its nominal range of validity at other values of d/R as well. Regardless, with "feature 1," we refer to the μ-independent fact that ω_l ∝ lv_Ae/d, and hence Δω_l ∝ v_Ae/d at large l and large d/R. Feature 2, on the other hand, concerns the modes with ω_l that differ little from ω_crit. This turns out to be difficult to examine with Figure 8(a) because the mode frequencies become increasingly packed when d/R increases. In fact, Δω_l eventually becomes so small that we choose to exaggerate the fractional difference δ_l = ω_l /ω_crit − 1 by a factor of 10⁵ in Figure 8(b), where the horizontal dashed–dotted line again marks ω_l = ω_crit. Now one sees that ω_l consistently exceeds ω_crit for μ = 5, meaning that the modes are consistently oscillatory. When μ = 1.5, however, the first mode shows up as an evanescent mode for d/R = 50 × 3⁵ = 12, 150, and so does the second mode at the even larger d/R.

Zoom In Zoom Out Reset image size

Figure 9 further examines the modes with ω_l close to ω_crit for (a) μ = 1.5 and (b) μ = 5. Here the modulus of δ_l = ω_l /ω_crit − 1 for a given l is displayed as a function of d/R by the solid (dashed) curves when δ_l is positive (negative). Among the 50 modes examined, one mode out of five is plotted when l ranges from 10 to 50, whereas all of the first five are presented. For reference, the eigenfunctions ${\breve{v}}_{l}$ of the first three modes are given by Figure 10 for μ = 1.5 (left column) and μ = 5 (right column). A number of d/R are examined and can be directly read from the figure. Examine the case with μ = 5 first. Figure 9(b) indicates that δ_l is positive for all modes, the oscillatory nature of which is made clearer by the spatial behavior of ${\breve{v}}_{l}$ in Figure 10. Furthermore, δ_l for all l follows a 1/d²-dependence, shown by the blue dashed–dotted curve in Figure 9(b). We note that this 1/d²-dependence is not empirically found but rather inspired by the analytical behavior of δ_l for step density profiles when δ_l ≪ 1 (see Equation (A11) for details). That this dependence applies to the case μ = 5 reinforces the notion that the mode behavior is qualitatively similar when μ > 2. Now move on to the more interesting case where μ = 1.5. For any of the first three modes, Figure 9(a) indicates a transition from an oscillatory to an evanescent mode, as evidenced by the change of the sign of δ_l at some critical (d/R)_cri,l . When d/R becomes larger still, δ_l and hence ω_l become independent of d/R. In addition, (d/R)_cri,l is seen to increase with l, a feature that can be readily understood with the left column of Figure 10. Let D_l denote the spatial extent of the eigenfunction of an evanescent mode, meaning mathematically that Q(r) becomes negative when r > D_l (see Equations (26) and (27)). When multiple evanescent modes exist on a sufficiently large domain, their frequencies are necessarily such that ω₁ < ω₂ < ⋯ because the entire set {ω_l } is a monotonically increasing sequence with respect to l. On the other hand, it can be readily shown that the potential V(r) eventually approaches ${k}^{2}{v}_{\mathrm{Ae}}^{2}$ from below (above) when μ < 2 (μ > 2). It then follows that D₁ < D₂ < ⋯ . Consequently, (d/R)_cri,1 is necessarily smaller than (d/R)_cri,2 for the domain to accommodate the diminishing portion of the eigenfunction of the second mode. In fact, the sequence (d/R)_cri,l necessarily increases monotonically with respect to l. Figure 9(a) further indicates that at sufficiently large d/R, the oscillatory modes for μ = 1.5 remain characterized by δ_l ∝ 1/d² unless δ_l is extremely small.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

With the aid of Equation (28), we now offer some general remarks on the dependencies on the steepness parameter μ of both the modal behavior on a closed domain with large d/R and the solution to IVP 1. However, we choose to focus on the chosen [ρ_i /ρ_e , kR, Λ/R] = [2.25, π/15, 4] to avoid this manuscript becoming even longer. When μ < 2, more and more evanescent modes appear when d/R increases. With the exception of the first several, the oscillatory modes are such that their frequency spacing Δω_l starts with a 1/d²-dependence before eventually settling to a 1/d-dependence. All modes are oscillatory when μ > 2, and the corresponding Δω_l simply transitions from a 1/d²-dependence for small l to a 1/d-dependence for large l. Now focus on the two values of μ that we have adopted. Recall that evanescent modes are irrelevant when μ = 5 and make no contribution to the time-dependent solution when μ = 1.5. One therefore recognizes that only the second summation in Equation (28) matters. What results from Equation (28) when d/R increases is then an increasingly refined discretization of some Fourier integral over a continuum of ω extending from kv_Ae to infinity. The relevant terms of this integral, applicable to a truly open system (d/R → ∞), was explicitly worked out for μ = ∞ by ORT15 (see Appendix A.1 for details). Evidently, one by-product of our modal approach on a closed domain is the numerical distribution in the ω − r space of the terms in the Fourier integral, which cannot be expressed in closed form for general μ to our knowledge.

Supposing Λ/R is adjustable in view of Figure 5, we move on to demonstrate a generic condition for evanescent modes to be negligible. It suffices to adopt a truly open system and consider the signal at a specific location, such as r = R. We start with the assumption that evanescent modes do not contribute and deduce the condition that ensures this assumption. Let τ denote the extent of the duration of interest, by which we mean that the signal becomes too weak to discern when t > τ. Evidently, the outermost edge travels to a distance of D_τ . Recall that the spatial extent of the eigenfunctions of the evanescent modes D_l increases with l. One therefore deduces that evanescent modes are bound to be negligible when D_τ < D₁. On the other hand, D_τ is evidently lower than Λ + v_Ae τ because the speed at which the outermost edge travels (v_A(r)) is consistently lower than v_Ae. It then follows that evanescent modes can be neglected, provided

$$\begin{eqnarray}&&{\rm{\Lambda }}+{v}_{\mathrm{Ae}}\tau \lt {D}_{1}.\end{eqnarray} \tag{ 31 }$$

Equating τ to τ_ener in Figure 5, one readily finds that the inequality holds for all values of Λ examined therein, thereby explaining why the wave behavior is solely characterized by dispersive propagation but shows no sign of wave trapping. In fact, we can slightly generalize Equation (31) by supposing Λ ≫ R and adopting the worst-case scenario that τ is given by t₂₊ in Equation (19). We further neglect the deviation of v_A(r) from v_Ae, meaning that ${t}_{{1}^{-}}({\rm{\Lambda }}/2\to 0)\approx {\rm{\Lambda }}/2{v}_{\mathrm{Ae}}$ and ${t}_{{2}^{+}}(0\to {\rm{\Lambda }})\approx {\rm{\Lambda }}/{v}_{\mathrm{Ae}}$. A rather safe estimate for t₂₊ and hence τ is then 3Λ/2v_Ae. The inequality (Equation (31)) therefore becomes

$$\begin{eqnarray}&&{\rm{\Lambda }}\lt \displaystyle \frac{2{D}_{1}}{5}.\end{eqnarray} \tag{ 32 }$$

Note that the RHS evaluates to ∼4000 R in view of Figure 10 and further evaluates to ∼8 × 10⁶ km if we quote R ∼ 2000 km from Figure 1 in Schrijver (2007), in keeping with the adopted L/R = 15. One therefore deduces that trapped modes are unlikely to be relevant in the temporal evolution of axial fundamental sausage modes in AR loops, at least for the [ρ_i /ρ_e , L/R] examined here. We stress that trapped modes are allowed as eigensolutions on an open system when μ < 2, as pointed out by LN15. Likewise, we stress that their distinct dispersive behavior relative to trapped modes in the canonical ER83 equilibrium is relevant when large values of k are involved, such as happens for impulsive sausage wave trains in AR loops. This latter point is clear if one contrasts Figure 8 with Figure 3 in Yu et al. (2017). It is just that the existence of trapped modes for μ < 2 on an open system does not guarantee that they contribute to the temporal evolution of axial fundamentals.

Focusing on fast sausage modes (FSMs) in solar coronal loops, this study aimed at examining the consequences of some peculiar dispersive properties that may arise in an equilibrium differing from ER83 only by replacing the step with a continuous density profile (ρ₀(r)). By "peculiar," we mean that FSMs are not subject to cutoff axial wavenumbers when ρ₀(r) outside the cylinder possesses a sufficiently shallow r-dependence, which was first recognized on firm mathematical grounds by LN15 when analyzing the relevant eigenvalue problem (EVP) on a radially open system. Two effects follow. First, FSMs may be trapped regardless of the axial wavenumber k and density contrast ρ_i /ρ_e . Second, long-wavelength trapped FSMs are nearly dispersionless, with their axial phase speeds differing little from the external Alfvén speed. These two effects then led LN15 and Lopin & Nagorny (2019) to deduce that fast sausage perturbations of observable quality may exist in AR and flare loops alike, with their periodicities characterized by the longitudinal rather than the canonical transverse Alfvén time. If true, this deduction may substantially broaden the range of observed periodicities that FSMs can account for and therefore offer more seismological possibilities.

We made an effort to make our scope as narrow as possible by addressing the question "does the existence of trapped modes in EVPs on an open system guarantee that they play a role in determining the temporal behavior of sausage perturbations"? To be specific, we chose to work in the framework of linearized, pressureless, ideal MHD and specialize to an outer-μ density profile (Equation (17)). The solution to the relevant IVP on an open system (IVP 1) is then determined by the dimensionless parameters [ρ_i /ρ_e , μ; kR; Λ/R], where μ characterizes the steepness of the density profile outside the nominal radius R, and Λ represents the spatial extent of the initial perturbation. We focus on axial fundamentals in AR loops by taking [ρ_i /ρ_e , kR] = [2.25, π/15] in view of the observational constraints from EUV measurements (Aschwanden et al. 2004; Schrijver 2007). We distinguish between the cases with μ < 2 and μ > 2 because trapped modes are present (absent) in the former (latter) for the chosen [ρ_i /ρ_e , kR]. The IVP 1 is solved with both a direct FD approach and a modal approach, whereby the solution is expressed as the superposition of eigenmodes for the pertinent EVP on a closed domain (EVP 1). The dimensionless domain size d/R is involved in the latter approach, the evanescent modes of which are the counterparts of trapped modes on a truly open system. Our findings can be summarized as follows.

The answer to the question we laid out is "no." We came to this conclusion primarily because the FD solutions for a substantial range of μ and Λ/R are consistently characterized by some dispersive propagation but show no sign of wave trapping (Figure 5). In particular, the solutions show a smooth transition when μ crosses the nominally critical value of 2. With the modal approach, we showed that more and more evanescent modes appear when the domain size increases, thereby lending further support to the recognition of LN15. However, even the shortest spatial extent of the evanescent eigenfunctions is well beyond the observationally reasonable range of the spatial extent of the initial perturbations. ⁶ Consequently, the initial perturbations cannot impart a discernible fraction of energy to the evanescent modes, which in turn means that these modes do not contribute to the temporal evolution of the system.

Before closing, we offer some further remarks on the influence on coronal FSMs due to the deviation of the equilibrium from ER83. We start by noting that the formulation of the transverse structuring actually offers a mixed message for FSMs in terms of their observational applications. On the one hand, sausage-like perturbations are robust in the sense that they are permitted even when waveguides are not strictly axisymmetric, as happens for waveguides with, say, elliptic (e.g., Erdélyi & Morton 2009; Aldhafeeri et al. 2021) or even irregular cross sections (Aldhafeeri 2021; Guo et al. 2021). On the other hand, that the dispersive behavior of FSMs in non-ER83 equilibria may be qualitatively different does not necessarily mean that FSMs can be invoked to interpret a broader range of periodicities. That said, one cannot rule out coronal FSMs as an interpretation for oscillations with periodicities on the order of the longitudinal Alfvén time either. Let us name only one possible equilibrium configuration where essentially, the only difference from ER83 is the introduction of a magnetically twisted boundary layer. Radial fundamental FSMs may be trapped regardless of the axial wavenumber (e.g., Khongorova et al. 2012; Lim et al. 2018), the corresponding eigenfunctions possessing spatial scales that do not far exceed the cylinder radius (e.g., Mikhalyaev & Bembitov 2014; Lopin 2021). One readily deduces that these radial fundamental FSMs may indeed show up in reality, even though a definitive answer relies on a detailed study from the IVP perspective.

We thank the referee for constructive comments. This research was supported by the National Natural Science Foundation of China (41974200 and 11761141002). We gratefully acknowledge ISSI-BJ for supporting the workshop on "Oscillatory Processes in Solar and Stellar Coronae," during which this study was partly initiated.

A.1. Solution to IVP 1 in Terms of Eigenmodes for an Open System

This subsection presents the formal solution to IVP 1 expressed as the superposition of eigenmodes for a laterally unbounded system. We closely follow the Fourier integral–based approach that ORT15 adopted to examine the 2D propagation of sausage wave trains in a system that is unbounded in the axial direction as well. The modal solution to our IVP 1 is actually part of their 2D solution. Consequently, only slight revisions to Equation (25) in ORT15 are needed to ensure dimensional and notational consistency. We choose to provide a minimal set of equations leading to the modal solution for two reasons. One, some equations will find immediate applications to the pertinent EVP on a closed domain. Two, the conceptual understanding embedded in the formal solution to IVP 1 should be informative for future studies on standing sausage perturbations in generic coronal structures. To this end, we see the axial wavenumber (k) as given and arbitrary.

The following notations are necessary. We use J_n and Y_n to represent Bessel functions of the first and second kind, respectively. Likewise, I_n and K_n denote modified Bessel functions of the first and second kind, respectively. Only orders n = 0 and 1 are relevant. In particular, j_n,m denotes the mth zero of J_n with m = 1, 2, ⋯ . The cutoff wavenumbers are then given by

$$\begin{eqnarray}&&{k}_{\mathrm{cutoff},{\rm{m}}}R=\displaystyle \frac{{j}_{0,m}}{\sqrt{{\rho }_{i}/{\rho }_{e}-1}}.\end{eqnarray} \tag{ A1 }$$

Defining

$$\begin{eqnarray}\begin{array}{rcl}{k}_{i}^{2} & = & \displaystyle \frac{{\omega }^{2}-{k}^{2}{v}_{\mathrm{Ai}}^{2}}{{v}_{\mathrm{Ai}}^{2}},\\ {k}_{e}^{2} & = & \displaystyle \frac{{\omega }^{2}-{k}^{2}{v}_{\mathrm{Ae}}^{2}}{{v}_{\mathrm{Ae}}^{2}},\\ {\kappa }_{e}^{2} & = & -\displaystyle \frac{{\omega }^{2}-{k}^{2}{v}_{\mathrm{Ae}}^{2}}{{v}_{\mathrm{Ae}}^{2}}=-{k}_{e}^{2},\end{array}\end{eqnarray} \tag{ A2 }$$

we note that ω always exceeds kv_Ai; hence, ${k}_{i}^{2}$ are always positive.

The modal solution to IVP 1 involves both proper and improper modes, which are discriminated by the sign of ${k}_{e}^{2}$. Proper modes (${k}_{e}^{2}\lt 0$) are relevant when k exceeds k_cutoff,1 and correspond to a discrete set of frequencies. Let j label a proper mode. Its eigenfunction reads

$$\begin{eqnarray}{\breve{v}}_{j}(r)=\left\{\begin{array}{ll}-\displaystyle \frac{{v}_{\mathrm{Ai}}}{{k}_{i}R}{K}_{0}({\kappa }_{e}R){J}_{1}({k}_{i}r), & 0\leqslant r\leqslant R,\\ \displaystyle \frac{{v}_{\mathrm{Ai}}}{{\kappa }_{e}R}{J}_{0}({k}_{i}R){K}_{1}({\kappa }_{e}r), & r\gt R.\end{array}\right.\end{eqnarray} \tag{ A3 }$$

Written in this form, Equation (A3) ensures that $d{\breve{v}}_{j}(r)/{dr}$ is continuous. The mode frequency is then dictated by the continuity of ${\breve{v}}_{j}$ itself across r = R, which leads to the well-known dispersion relation (DR; e.g., ER83; see also Meerson et al. 1978; Spruit 1982; Cally 1986)

$$\begin{eqnarray}&&{k}_{i}\displaystyle \frac{{J}_{0}({k}_{i}R)}{{J}_{1}({k}_{i}R)}+{\kappa }_{e}\displaystyle \frac{{K}_{0}({\kappa }_{e}R)}{{K}_{1}({\kappa }_{e}R)}=0\,.\end{eqnarray} \tag{ A4 }$$

Improper modes (${k}_{e}^{2}\gt 0$) are relevant regardless of k, their frequencies continuously spanning the range (kv_Ae, ∞). The eigenfunction reads

$$\begin{eqnarray}{\breve{v}}_{\omega }(r)=\left\{\begin{array}{ll}-{v}_{\mathrm{Ai}}\displaystyle \frac{{k}_{e}^{2}R}{{k}_{i}}{J}_{1}({k}_{i}r), & 0\leqslant r\leqslant R,\\ -{v}_{\mathrm{Ai}}({k}_{e}R)\left[{C}_{J}{J}_{1}({k}_{e}r)+{C}_{Y}{Y}_{1}({k}_{e}r)\right], & r\gt R,\end{array}\right.\end{eqnarray} \tag{ A5 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{C}_{J} & = & \displaystyle \frac{\pi {k}_{e}R}{2{k}_{i}}\left[-{k}_{i}{J}_{0}({k}_{i}R){Y}_{1}({k}_{e}R)+{k}_{e}{J}_{1}({k}_{i}R){Y}_{0}({k}_{e}R)\right],\\ {C}_{Y} & = & \displaystyle \frac{\pi {k}_{e}R}{2{k}_{i}}\left[-{k}_{e}{J}_{1}({k}_{i}R){J}_{0}({k}_{e}R)+{k}_{i}{J}_{0}({k}_{i}R){J}_{1}({k}_{e}R)\right].\end{array}\end{eqnarray} \tag{ A6 }$$

Both ${\breve{v}}_{\omega }$ and $d{\breve{v}}_{\omega }(r)/{dr}$ are continuous. The modal solution eventually reads

$$\begin{eqnarray}&&\begin{array}{l}\hat{v}(r,t)=\displaystyle \sum _{j=1}^{J}{c}_{j}{\breve{v}}_{j}(r)\cos ({\omega }_{j}t)+{\displaystyle \int }_{{{kv}}_{\mathrm{Ae}}}^{\infty }{S}_{\omega }{\breve{v}}_{\omega }(r)\cos (\omega t)d\omega ,\\ 0\lt r\lt \infty ,\quad 0\lt t\lt \infty ,\end{array}\end{eqnarray} \tag{ A7 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{c}_{j} & = & \displaystyle \frac{{\displaystyle \int }_{0}^{\infty }u(r){\breve{v}}_{j}(r){\rho }_{0}(r){rdr}}{{\displaystyle \int }_{0}^{\infty }{\breve{v}}_{j}^{2}(r){\rho }_{0}(r){rdr}},\\ {S}_{\omega } & = & \displaystyle \frac{\omega {\displaystyle \int }_{0}^{\infty }u(r){\breve{v}}_{\omega }(r){\rho }_{0}(r){rdr}}{({\rho }_{e}{v}_{\mathrm{Ai}}^{2}{R}^{2}){\left({k}_{e}{v}_{\mathrm{Ae}}\right)}^{2}({C}_{J}^{2}+{C}_{Y}^{2})}.\end{array}\end{eqnarray} \tag{ A8 }$$

Let us summarize the steps to solve IVP 1 given ρ_i /ρ_e , k, and u(r). First of all, with Equation (A1), one counts J, the number of cutoff wavenumbers that are smaller than k. Evidently, J = 0 if k < k_cutoff,1, making proper modes irrelevant. Second, if J ≥ 1, then for each allowed j, one evaluates its eigenfrequency ω_j and then its eigenfunction ${\breve{v}}_{j}$ with Equations (A4) and (A3), respectively. The contribution of the proper mode, c_j , can then be found with Equation (A8). Third, with Equations (A5) and (A6), one evaluates the improper eigenfunction for any ω > kv_Ae. The contribution from the improper mode, S_ω , is then readily found with Equation (A8).

A.2. Eigenmodes for a Closed System

This subsection examines some properties of the eigenmodes for a closed system, namely, the solutions to EVP 1 specialized to a step density profile. We start with a concrete example found for [ρ_i /ρ_e , kR] = [2.25, π/15] and an initial perturbation given by Equation (18) with Λ/R = 4. We solve IVP 1 with three independent methods. The first, to be called "modal open," is based on eigenmodes on an open system, whereas the second ("modal closed") is based on Equation (23) for a domain size d = 50 R. The two sets of solutions agree with each other exactly and further agree with the solution found with the simpler FD approach. We choose not to present the comparison among the time-dependent solutions because a comparison between the frequency dependencies of the modal contributions seems more informative but is unavailable as far as we know. Note that these contributions for the chosen [ρ_i /ρ_e , kR] are solely due to improper (oscillatory) modes in the modal open (modal closed) approach. With Equations (28) and (A8) in mind, Figure 11 then specializes to r = R and compares the local spectral density ${S}_{l}{\breve{v}}_{l}(R)$ from the modal closed approach (asterisks) with ${S}_{\omega }{\breve{v}}_{\omega }(R)$ from the modal open approach (blue solid curve). It is reassuring to see that the solid curve exactly threads the symbols, meaning that the continuum of improper modes is adequately resolved by the discrete oscillatory modes despite the rather modest domain size.

Zoom In Zoom Out Reset image size

We now focus on the discrete modes themselves by recalling the discussions immediately following Equation (26). First, ${k}_{i}^{2}\gt 0$ is guaranteed because ω_l necessarily exceeds kv_Ai. Second, the mode classification is eventually determined by how ω_l compares with kv_Ae, meaning that ${k}_{e}^{2}$ is positive (negative) when a mode is oscillatory (evanescent). Regardless, the eigenfunction ${\breve{v}}_{l}(r)$ can be consistently described by Equation (A5), the reason being that J₁(k_e r) and Y₁(k_e r) always form a numerically satisfactory pair in the outer region (R < r < d). ⁷ The requirement ${\breve{v}}_{l}(r=d)=0$ then gives a DR that governs the mode frequency ω_l ,

$$\begin{eqnarray}&&{C}_{J}{J}_{1}({k}_{e}d)+{C}_{Y}{Y}_{1}({k}_{e}d)=0.\end{eqnarray} \tag{ A9 }$$

We now specialize to the situation where k < k_cutoff,1 to better connect with the main text, the associated analytical progress being new to our knowledge.

Let us examine the analytical behavior of the modes with frequencies ω_l just above ω_crit = kv_Ae. Expressing ω_l as kv_Ae(1 + δ) with 0 < δ ≪ 1, one finds ${k}_{i}^{2}\approx {k}^{2}({\rho }_{i}/{\rho }_{e}-1)$ and ${k}_{e}^{2}\approx {k}^{2}(2\delta )$ (see Equation (A2)). Now suppose that k_e R ≪ 1. The approximate expressions of J_n and Y_n for small arguments then indicate that Y₁(k_e R) ∼ 1/(k_e R) is the most singular term in the coefficients C_J and C_Y (see Equation (A6)). The left-hand side of the DR (Equation (A9)) is therefore dominated by the first term, meaning that k_e d ≈ j_1,l . In other words,

$$\begin{eqnarray}&&\delta =\displaystyle \frac{{\omega }_{l}}{{{kv}}_{\mathrm{Ae}}}-1\approx \displaystyle \frac{{j}_{1,l}^{2}}{2{\left({kd}\right)}^{2}}.\end{eqnarray} \tag{ A10 }$$

Given the assumptions δ ≪ 1 and k_e R ≪ 1, the modes in question are characterized by

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\omega }_{l}}{{{kv}}_{\mathrm{Ae}}}\approx 1+\displaystyle \frac{{j}_{1,l}^{2}}{2{\left({kd}\right)}^{2}},\\ \ {\rm{provided}}\,d/R\gg {j}_{1,l},\quad {\left({kd}\right)}^{2}\gg {j}_{1,l}^{2}.\end{array}\end{eqnarray} \tag{ A11 }$$

Now consider high-frequency modes in a system with d/R ≫ 1. By "high," we assume that

$$\begin{eqnarray}&&{\omega }_{l}\gg {v}_{\mathrm{Ae}}/R,\quad {\omega }_{l}\gg {{kv}}_{\mathrm{Ae}}.\end{eqnarray} \tag{ A12 }$$

It then follows from Equation (A2) that k_i R, k_e R ≫ 1 and hence k_i d, k_e d ≫ 1. With the expressions for J_n and Y_n at large arguments, some algebra indicates that the DR (Equation (A9)) approximates to

$$\begin{eqnarray}&&\sin [{k}_{e}(d-R)+{\phi }_{e}]=0,\end{eqnarray} \tag{ A13 }$$

where ϕ_e satisfies the relation

$$\begin{eqnarray}&&\tan {\phi }_{e}=\sqrt{{\rho }_{e}/{\rho }_{i}}\tan \left[\sqrt{{\rho }_{i}/{\rho }_{e}}({k}_{e}R)-\displaystyle \frac{\pi }{4}\right].\end{eqnarray} \tag{ A14 }$$

Note that Equation (A13) is implicit in ω_l because ω_l is involved in Equation (A14). Note further that the range of ϕ_e is not restricted by Equation (A14) per se. However, as can be verified a posteriori, ϕ_e is negligible to leading order, meaning that k_e ≈ l π/(d − R) ≈ l π/d. If desired, this solution can be plugged into Equation (A14) to yield the first-order correction due to ϕ_e , resulting in

$$\begin{eqnarray}&&{\omega }_{l}\approx \left(l\pi -{\phi }_{l}\right)\displaystyle \frac{{v}_{\mathrm{Ae}}}{d-R}.\end{eqnarray} \tag{ A15 }$$

Here ϕ_l is given by

$$\begin{eqnarray}\begin{array}{rcl}{\phi }_{l} & = & \arctan \left[\sqrt{{\rho }_{e}/{\rho }_{i}}\tan \left(\sqrt{{\rho }_{i}/{\rho }_{e}}\displaystyle \frac{l\pi }{d/R-1}-\displaystyle \frac{\pi }{4}\right)\right]\\ & & +\unicode{x0230A}\sqrt{{\rho }_{i}/{\rho }_{e}}\displaystyle \frac{l}{d/R-1}+\displaystyle \frac{1}{4}\unicode{x0230B}\pi ,\end{array}\end{eqnarray} \tag{ A16 }$$

with the floor function ⌊ · ⌋ introduced to make ϕ_l continuous with respect to l. Restricting oneself to those modes with l ≫ d/(π R) and l ≫ kd/π in view of the assumptions (Equation (A12)), one recognizes that the RHS of Equation (A16) is necessarily far smaller than l π. Overall, it suffices for our purposes to summarize the properties of the high-frequency modes as

$$\begin{eqnarray}&&{\omega }_{l}\approx l\pi ({v}_{\mathrm{Ae}}/d),\quad {\rm{for}}\,l\gg \displaystyle \frac{d/R}{\pi },\quad l\gg \displaystyle \frac{{kd}}{\pi }.\end{eqnarray} \tag{ A17 }$$

This section examines some further properties of EVP 1 by capitalizing on the Schrödinger equation (Equation (26)). We recall that the equilibrium density is of the outer-μ type (see Equation (17)). This subsection extends Section A.2 in that μ is no longer restricted to being infinite.

We start by showing that all solutions to EVP 1 necessarily possess an ω_l that exceeds kv_Ai, regardless of μ or d/R. What we offer is only a slight generalization of the arguments given by BDBG96 for a step density profile. This generalization is possible because the arguments therein rely on only two conditions, one being that the potential V(r) consistently exceeds ${k}^{2}{v}_{\mathrm{Ai}}^{2}$ and the other being that the eigenfunction ${\breve{v}}_{l}(r)$ vanishes at both r = 0 and r = d. Neither condition is restricted to the particular μ = ∞. Now suppose that ω_l < kv_Ai for some mode, meaning that ${\omega }_{l}^{2}\lt V(r)$. The wave function Φ(r) and hence the eigenfunction ${\breve{v}}_{l}(r)$ then necessarily peak somewhere in the domain, diminishing toward both smaller and larger r. One then deduces that $d{\breve{v}}_{l}(r)/{dr}$ is discontinuous, thereby violating the continuity requirement for the Eulerian perturbation of total pressure.

We now address high-frequency modes in a system with d/R ≫ 1, "high" in the same sense as in the assumption (Equation (A12)). Our approach is essentially the classical WKB one detailed in Bender & Orszag (1999, Chapter 10). Somehow different is that we avoid the complication associated with the turning points (namely, where Q(r) = 0), which are bound to occur at small r for high-frequency modes because V(r) diverges at r = 0. This is done by handling the Schrödinger equation (Equation (26)) in the interior (r < R) and exterior (r > R) separately. First consider the exterior, where the condition ω² ≫ V(r) is ensured by Equation (A12). The leading order WKB solution to Equation (26) reads

$$\begin{eqnarray}&&{\rm{\Phi }}(r)\approx {A}_{e}{v}_{{\rm{A}}}^{1/2}(r)\sin [{\rm{\Theta }}(r)],\end{eqnarray} \tag{ B1 }$$

where

$$\begin{eqnarray*}&&{\rm{\Theta }}(r)=\omega {\int }_{R}^{r}\displaystyle \frac{{dr}^{\prime} }{{v}_{{\rm{A}}}(r^{\prime} )}+\phi ,\end{eqnarray*}$$

with A_e and ϕ being constants. Requiring Φ(d) = 0 then leads to

$$\begin{eqnarray}&&{\omega }_{l}\approx \displaystyle \frac{l\pi -{\phi }_{l}}{{\int }_{R}^{d}{dr}^{\prime} /{v}_{{\rm{A}}}(r^{\prime} )}\end{eqnarray} \tag{ B2 }$$

for mode l. Evidently, the integral in Equation (B2) is well approximated by d/v_Ae for large d/R. The high-frequency modes are therefore still characterized by Equation (A17), provided ∣ϕ_l ∣ ≪ l π.

We proceed to show that the inequality ∣ϕ_l ∣ ≪ l π indeed holds via three mutually complementary methods. First, we numerically solve EVP 1 for a substantial number of combinations of [ρ_i /ρ_e , μ; k; d/R]. Comparing the computed ω_l with the RHS of Equation (B2) then indicates that ∣ϕ_l ∣ ≪ l π. Second, we make some analytical progress to estimate ϕ_l . Now the interior (r < R) needs to be examined, for which purpose Equation (A5) indicates that the exact solution to Equation (26) reads Φ(r) = A_i r^1/2 J₁(k_i r), with A_i being an arbitrary constant. In addition, the condition k_i R ≫ 1 is ensured by the assumption (Equation (A12)), enabling one to employ the approximate expressions for Bessel functions at large arguments to find

$$\begin{eqnarray}&&{\left.\displaystyle \frac{d{\rm{\Phi }}/{dr}}{{\rm{\Phi }}}\right|}_{r={R}^{-}}\approx {k}_{i}\cot \left({k}_{i}R-\pi /4\right)-\displaystyle \frac{1}{2R}.\end{eqnarray} \tag{ B3 }$$

One further finds with the external WKB solution (Equation (B1)) that

$$\begin{eqnarray}&&{\left.\displaystyle \frac{d{\rm{\Phi }}/{dr}}{{\rm{\Phi }}}\right|}_{r={R}^{+}}\approx {k}_{i}\cot {\phi }_{l}+{\left.\displaystyle \frac{1}{2}\displaystyle \frac{d\mathrm{ln}{v}_{{\rm{A}}}}{{dr}}\right|}_{r={R}^{+}}.\end{eqnarray} \tag{ B4 }$$

Equating Equation (B3) to Equation (B4) then leads to

$$\begin{eqnarray}&&\cot {\phi }_{l}\approx \cot \left({k}_{i}R-\pi /4\right)-\displaystyle \frac{1}{2}\left(\displaystyle \frac{1}{{k}_{i}R}+{\left.\displaystyle \frac{1}{{k}_{i}}\displaystyle \frac{d\mathrm{ln}{v}_{{\rm{A}}}}{{dr}}\right|}_{r={R}^{+}}\right).\end{eqnarray} \tag{ B5 }$$

Now suppose that the second term on the RHS of Equation (B5) is negligible, implying that μ is not too large. Let us further see the RHS of Equation (B5) as known from Equation (A17). It turns out that ϕ_l can be approximated to leading order by ${\phi }_{l}\approx {k}_{i}R-\pi /4\approx l\pi \sqrt{{\rho }_{i}/{\rho }_{e}}/(d/R)$, meaning that ∣ϕ_l ∣ ≪ l π for sufficiently large d/R. Third, we offer some heuristic arguments to estimate ∣ϕ_l ∣ for arbitrary μ > 1. Let N_int (N_ext) be the number of extrema in the eigenfunction ${\breve{v}}_{l}(r)$ in the interior r < R (the exterior R < r < d). Evidently, ϕ_l in Equation (B2) stems from the influence of the interior, making ∣ϕ_l /π∣/l essentially a measure of N_int/N_ext. With d/R ≫ 1, one readily deduces that N_int ≪ N_ext and hence ∣ϕ_l ∣ ≪ l π.

This section examines some subtlety associated with our outer-μ formulation (Equation (17)) for the equilibrium density ρ₀(r). Let us recall our argument that, for the fixed pair [ρ_i /ρ_e , L/R] = [2.25, 15] and some chosen μ < 2, the dimensionless spatial extent of the initial perturbation (Λ/R) needs to be unrealistically large for evanescent modes to be nonnegligible in the time-dependent solutions to IVP 1. However, Equation (17) indicates that the spatial range containing the density enhancement broadens when μ decreases, making the nominal radius R less and less ideal for characterizing the spatial variation of ρ₀(r). Let the spatial scale of ρ₀(r) be measured by some effective radius R_eff. One may question whether our argument still holds because there may be a regime where evanescent modes are visible for not-so-extreme values of Λ measured in units of R_eff rather than R. Somehow it is nontrivial to quantify this aspect in an exhaustive manner, to explain which we note that we will exclusively adopt the FD approach to solve IVP 1 here for computational convenience. We will additionally fix the steepness parameter at μ = 1.5 but experiment with two different values for the density contrast ρ_i /ρ_e . Let τ denote the duration to be examined in an FD solution and D_τ the distance that the outer edge of the perturbation reaches when t = τ. For now, consider the modal structure for EVP 1 on a sufficiently large domain of size d, where, by "sufficiently large," we mean that d ≫ D_τ such that a multitude of evanescent modes exist. From Section 4, we know that oscillatory modes are always permitted, and the frequencies of evanescent modes (ω < ω_crit = kv_Ae) may not differ much from those of the low-frequency oscillatory modes (ω ≳ ω_crit). Now that low-frequency oscillatory modes are excited in general, it may take a considerable amount of time for them to interfere such that their contribution to a time-dependent solution eventually becomes sufficiently weak to make evanescent modes visible. In practice, however, the value of τ cannot be made infinite. Three regimes then arise in the signal behavior within a large but nonetheless finite time frame, where evanescent modes are not discernible, somehow discernible but weak, and stronger than some threshold, respectively. We see evanescent modes as observationally relevant only when the last regime occurs.

Some definitions and remarks are necessary here. We define R_eff as the radial distance where the function f(r) in Equation (17) attains 1/10, a factor that is meant to be small but admittedly arbitrary. One nonetheless finds that R_eff/R = 10^1/μ and evaluates to 4.64 for μ = 1.5. We see R_eff rather than R as being observationally relevant and consistently use R_eff to measure Λ and the loop length L. Likewise, time will be measured in units of the longitudinal Alfvén time τ_long = 2π/ω_crit = 2L/v_Ae, which is more relevant for the large-time behavior. We examine only the time frame t ≲ 40τ_long, which is seen as sufficiently long. Overall, the time-dependent solutions to IVP 1 are determined by the set of dimensionless parameters [ρ_i /ρ_e , L/R_eff, Λ/R_eff] given a fixed μ = 1.5. We deem the range Λ/R_eff ≤ 20 as observationally realistic for initial perturbations. On the other hand, we will examine the following two quantities to assess the significance of evanescent modes. The first is Γ(t) = E_tot(R_eff, t)/E₀, with E₀ = E_tot(Λ, t = 0) being the total energy initially deposited to the entire system. The second is simply the instantaneous radial speed at the effective radius $\hat{v}({R}_{\mathrm{eff}},t)$. While both Γ(t) and $\hat{v}({R}_{\mathrm{eff}},t)$ measure the signal strength in the volume r ≤ R_eff, we find that the former can better bring out the differences when L/R_eff or Λ/R_eff varies.

We start by examining an AR loop with [ρ_i /ρ_e , L/R_eff] = [10, 10], which is only marginally realistic because in general, ρ_i /ρ_e (L/R_eff) is lower (larger) in observations (e.g., Aschwanden et al. 2004; Schrijver 2007). Figure 12 displays the temporal profiles of (a) Γ(t) and (b) $\hat{v}({R}_{\mathrm{eff}},t)$ for a number of Λ/R_eff as labeled. Examining any Λ/R_eff in the chosen time frame, one sees that a periodic behavior develops at large t for Γ(t) and $\hat{v}({R}_{\mathrm{eff}},t)$ alike. We will focus on this periodic stage here and hereafter. A slight difference between Figures 12(a) and (b) is that the period in Γ(t) is half that in $\hat{v}({R}_{\mathrm{eff}},t)$, which arises simply because E_tot involves the perturbations as squared terms by definition (see Equation (12)). More importantly, the signal strengthens when Λ/R_eff increases, as can be discerned in Figure 12(b) and seen more clearly in Figure 12(a). Regardless, the signal for any Λ/R_eff weakens monotonically with time, eventually resulting in extremely small values for both Γ(t) and $\hat{v}({R}_{\mathrm{eff}},t)$. Note that this is true even for Λ/R_eff = 40, which exceeds the range that we deem observationally realistic. Note further that the signal in the periodic stage tends to weaken when L/R_eff increases or ρ_i /ρ_e decreases from the value we choose. Figure 12 therefore means that evanescent modes are not discernible for realistic combinations of [ρ_i /ρ_e , L/R_eff], thereby strengthening our conclusion that the existence of evanescent modes in the pertinent EVP analysis does not guarantee their relevance in the temporal evolution of the system. It then follows that one needs to invoke, say, kink modes to account for a periodicity on the order of the longitudinal Alfvén time when analyzing oscillating AR loops, even given our outer-μ formulation. Furthermore, whether an interpretation in terms of kink modes is justifiable can be readily assessed by looking for the telltale signature of transverse displacements because AR loops tend to be imaged with high spatial resolution on a routine basis.

Zoom In Zoom Out Reset image size

That said, evanescent modes may indeed be discernible or visible if one experiments with, say, drastically different values of ρ_i /ρ_e . We proceed with a fixed combination [ρ_i /ρ_e , Λ/R_eff] = [100, 20], the chosen density contrast being relevant for flare loops (e.g., Aschwanden et al. 2004). Figure 13 presents the temporal profiles of (a) Γ(t) and (b) $\hat{v}({R}_{\mathrm{eff}},t)$ for a variety of dimensionless loop lengths L/R_eff as labeled. Three features are evident regarding the periodic stage. First, Figure 13(b) indicates that the signal for any L/R_eff further settles to a stage where $\hat{v}({R}_{\mathrm{eff}},t)$ actually possesses two periodicities, the dominant one (P_domi) being close to but nonetheless longer than τ_long. In addition, the signal amplitude is modulated by a second period (P_env ≫ P_domi), ⁸ as can be seen more clearly in Figure 13(a). Second, when L/R_eff increases, the signal weakens, and P_domi (P_env) in units of τ_long slightly decreases (increases). This feature can be seen in both Figure 13(a) and (b) but is clearer in the latter. Third, and more importantly, the signal envelope fluctuates about a time-independent level, which can be seen by examining the profiles of the maxima/minima of either Γ(t) or $\hat{v}({R}_{\mathrm{eff}},t)$. Given this feature, we take the relevance of evanescent modes as self-evident and focus on the more interesting L/R_eff-dependencies of P_domi/τ_long and P_env/τ_long. Evidently, the envelope modulation stems from the beat among evanescent modes, which themselves possess frequencies that are marginally lower than ω_crit = 2π/τ_long. We recall that a larger L/R_eff means a larger L/R and hence a smaller dimensionless axial wavenumber kR. We proceed to consider evanescent eigenmodes of some given radial harmonic number l (say, l = 1 or 2, as in Figure 7 of Yu et al. 2017). Let ω denote the eigenfrequency, P = 2π/ω the eigenperiod, and v_ph = ω/k the phase speed. One readily finds P/τ_long = ω_crit/ω = v_Ae/v_ph. The reason for P_domi/τ_long to decrease with L/R_eff is that the phase speeds of evanescent modes increase toward v_Ae when kR decreases (see Figure 7 in Yu et al. 2017). On the other hand, that P_env increases with L/R_eff is because a reduction in kR makes the difference in the values of v_ph/v_Ae for modes with adjacent values of l smaller.

Zoom In Zoom Out Reset image size

That evanescent modes can be seen in Figure 13 does not necessarily mean that they are observationally relevant. Consider only a fixed ρ_i /ρ_e = 100 and only the periodic stage. Let A denote the maximum amplitude of $\hat{v}({R}_{\mathrm{eff}},t)$ in units of the magnitude of the initial perturbation. Experimenting with a substantial number of combinations [L/R_eff, Λ/R_eff], we find that A always decreases (increases) with increasing L/R_eff (Λ/R_eff) when the other parameter is fixed, at least when 8 ≤ L/R_eff ≤ 15 and 0 < Λ/R_eff ≤ 20. Suppose that evanescent modes can be seen as observationally relevant only when A exceeds some threshold A_c and let ${(L/{R}_{\mathrm{eff}})}_{{\rm{c}}}$ denote the critical L/R_eff beyond which A < A_c. Seeing A_c as variable, this then means that ${(L/{R}_{\mathrm{eff}})}_{{\rm{c}}}$ is a function of A_c. Figure 14 displays the A_c-dependence of ${(L/{R}_{\mathrm{eff}})}_{{\rm{c}}}$ for a number of Λ/R_eff as labeled; before describing them, we need to mention some technical subtlety. We consistently adopt the time frame t ≤ 40τ_long to determine A and eventually construct a symbol in Figure 14. However, Figure 13(a) has already shown that P_env becomes longer when L/R_eff increases, making it possible that the amplitude maximum is not captured for t ≤ 40τ_long if P_env is too long. This turns out not to be a real concern for two reasons, one being that one needs to specify a time frame in any case and the other being that the envelope modulation is already clear for t ≤ 40τ_long. We will return to this point later. Now let "allowed range" refer specifically to those L/R_eff for which evanescent modes are visible. One sees from Figure 14 that ${(L/{R}_{\mathrm{eff}})}_{{\rm{c}}}$ decreases monotonically with A_c for a given Λ/R_eff, meaning that the allowed range narrows when the relevant instrumental sensitivity weakens. Likewise, ${(L/{R}_{\mathrm{eff}})}_{{\rm{c}}}$ increases when Λ/R_eff increases, meaning that the allowed range broadens with increasing Λ/R_eff for a given sensitivity and hence a given A_c.

Zoom In Zoom Out Reset image size

It should be informative to place Figure 14 in the ER83 context. Consider an "ER" loop, by which we mean a loop with an equilibrium density ρ₀(r) that equals ρ_i for r < R_eff but ρ_e otherwise, and restrict to axial fundamentals. A critical ${(L/{R}_{\mathrm{eff}})}_{\mathrm{ER}}$ then follows from Equation (A1),

$$\begin{eqnarray}&&{(L/{R}_{\mathrm{eff}})}_{\mathrm{ER}}=\displaystyle \frac{\pi \sqrt{{\rho }_{i}/{\rho }_{e}-1}}{{j}_{\mathrm{0,1}}},\end{eqnarray} \tag{ C1 }$$

only below which evanescent modes are possible in an ER loop. With Equation (C1), one finds that ${(L/{R}_{\mathrm{eff}})}_{\mathrm{ER}}$ evaluates to 13 for ρ_i /ρ_e = 100, and this value is represented by the horizontal dashed–dotted line in Figure 14. Two situations arise regarding the relevance of ${(L/{R}_{\mathrm{eff}})}_{\mathrm{ER}}$; to describe them, we assume that L/R_eff is known observationally and that a time series is available for the pertinent observable over tens of longitudinal Alfvén times. First, consider the situation where L/R_eff is observed to exceed ${(L/{R}_{\mathrm{eff}})}_{\mathrm{ER}}$. Suppose that the observed signal eventually possesses a stage where the amplitude does not diminish with time. It then follows from the left corner of Figure 14 that this stage cannot be interpreted in the ER83 framework because no evanescent modes can be excited regardless of Λ. In this regard, a diffuse loop model such as that formulated by the outer-μ profile may indeed broaden the range of flare-associated QPPs that FSMs may account for. However, this broadening is rather limited because it happens only when the relevant instrument is highly sensitive and Λ/R_eff is nearly beyond the range that we deem observationally realistic. Now consider the situation where L/R_eff is observed to be smaller than ${(L/{R}_{\mathrm{eff}})}_{\mathrm{ER}}$. Evanescent modes are allowed by ER loops in this case and visible for the majority of combinations [A_c, Λ/R_eff] explored in Figure 14 for our outer-μ loops. One may then question whether the temporal profile of the pertinent observable can help categorize the involved flare loop as an ER loop or an outer-μ one. This is indeed possible; to illustrate this, we consider the case where L/R_eff = 11. Suppose that the flare loop is describable as an ER loop. Equation (C1) then yields a value of 5.66 if one replaces j_0,1 = 2.41 therein with j_0,2 = 5.52. This means that the evanescent mode with a radial harmonic number l = 2 is not relevant, and the pertinent signal eventually settles to a monochromatic variation with a constant amplitude. Now suppose that the flare loop agrees with an outer-μ loop for which μ is presumed to be 1.5, and suppose further that Λ/R_eff = 10. It turns out that the $\hat{v}({R}_{\mathrm{eff}},t)$ profile is very similar to the black curve in Figure 13(b), even though the amplitude maximum A in the periodic stage attains a smaller value of ∼0.1. Nonetheless, if an amplitude of this magnitude can be resolved, then the relevant observable will be seen to experience some envelope modulation, which is distinct from what is expected for an ER loop. In addition, this amplitude modulation may also be useful to distinguish between kink modes and FSMs in flare loops if the flare loops are not well spatially resolved to tell whether they are experiencing transverse displacements. We refrain from discussing this aspect further for two reasons, one being that kink modes in our outer-μ setup have yet to be examined, and the other being that a dedicated forward modeling approach seems necessary to establish the detailed observational signatures. Rather, with Figure 14, we conclude that whether evanescent modes in outer-μ loops can be observed depends critically on instrumental sensitivity. We conclude further that the ER83 cutoff ${(L/{R}_{\mathrm{eff}})}_{\mathrm{ER}}$ serves as a useful reference in that evanescent modes in outer-μ loops are hardly observable when $L/{R}_{\mathrm{eff}}\geqslant {(L/{R}_{\mathrm{eff}})}_{\mathrm{ER}}$ unless the pertinent instrument is extremely sensitive.