© The Authors 2024

1. Introduction

The past two decades have seen considerable progress in solar coronal seismology (e.g., the reviews by De Moortel & Nakariakov 2012; Nakariakov & Kolotkov 2020). One key factor has been the increasingly refined theories for magnetohydrodynamic (MHD)¹ waves in structured media (see recent topical collections by Nakariakov et al. 2022; Kolotkov et al. 2023). This can be said even for waves in such canonical equilibria as straight cylinders where the physical parameters are only structured radially (see the textbooks by Roberts 2019 and Goedbloed et al. 2019). Indeed, kink motions have been better understood theoretically, facilitating fuller seismological applications of, say, the cyclic transverse displacements extensively observed in active region (AR) loops (e.g., see the dedicated review by Nakariakov et al. 2021). Likewise, substantial theoretical progress has been made for sausage motions (see the review by Li et al. 2020), enabling one to better exploit the abundantly measured quasi-periodic pulsations (QPPs) in flares, for instance (see Van Doorsselaere et al. 2016 and Zimovets et al. 2021 for reviews).

Whichever observable is used for seismology, an unambiguous theoretical understanding needs to be offered, which is not always possible. It suffices to regard “observables” as such timescales as periods or damping times. It also suffices to consider kink motions in a z-aligned cylindrical setup and only consider those temporal intervals where one may invoke the expectations in classic mode analyses. By “classic”, we mean the customary procedure that starts by Fourier-decomposing any perturbation as exp[−i(Ωt − kz)], where k is the axial wave number, and the possibly complex-valued Ω is found by seeking nontrivial solutions to some boundary value problem (BVP; see, e.g., Chapter 6 in Roberts 2019). We consider only small k given that AR loops are of primary interest throughout. We additionally use “long” (“short”) to refer to timescales on the order of the axial (transverse) Alfvén time. The periodicities at large times for kink motions are accepted as long and pertain to the kink frequency ℜΩ_kink (see, e.g., Eq. (1) in Edwin & Roberts 1983)², as has been extensively implemented (e.g., Nakariakov et al. 1999; Nakariakov & Ofman 2001). That said, some controversy exists as to whether the accompanying damping is solely due to the resonant absorption in the Alfvén continuum (see Goossens et al. 2011 for conceptual clarifications). The “principal fast leaky kink” mode (PFLK) with ℜΩ_PFLK ≈ ℜΩ_kink was additionally found in the mode analysis by Cally (1986; see also Meerson et al. 1978; Spruit 1982). Given a non-vanishing ℑΩ_PFLK, it was proposed by Cally (2003) that PFLK may occasionally play a role in the observed damping. However, Ruderman & Roberts (2006b) argued otherwise, thereby leading to some extensive discussions regarding whether the timescales of PFLKs can make it into the system evolution (Cally 2006; Ruderman & Roberts 2006a; Andries & Goossens 2007). Evidently, this controversy may be partially settled by examining kink motions as an initial value problem (IVP). Terradas et al. (2007) offered such a study, reporting no evidence of PFLKs in their 1D solutions found with a finite-element code. Instead, temporal signatures were clearly demonstrated for the “trig” modes that solve the same dispersion relation (DR) that governs PFLKs in Cally (1986). We note that some key features of the expected spatial signatures can also manifest themselves, as evidenced by a recent 3D numerical study based on the finite-volume approach (Shi et al. 2023). We further note that trig modes are not of concern in the above-mentioned controversy, even though they are subject to some further debate.

Some clarifications on the nomenclature prove necessary. The term “classic BVPs” refers to those that arise in classic mode analyses; the associated governing equations are ordinary differential equations (ODEs), given that only 1D structuring is of interest. By “modes”, we specifically mean the nontrivial solutions to classic BVPs, and a mode is characterized by a “mode frequency” (Ω) and “mode function” pair. “Spectrum” refers to the collection of Ω. Two factors are typically responsible for ℑΩ not vanishing in our context: the physical relevance of the Alfvén continuum (e.g., Goedbloed 1998; Goossens et al. 2011) and the lateral domain being open. We focus on the latter. The boundary condition (BC) at infinity in classic analyses then amounts to the condition that no ingoing waves are allowed (e.g., Andries & Goossens 2007; Wang et al. 2023; Goedbloed et al. 2023). This BC yields two distinct types of mode behavior at large lateral distances, enabling the spectrum to be divided as such. One sub-spectrum comprises the well known “trapped modes”, the mode functions being evanescent at large distances and the mode frequencies being purely real (ℑΩ = 0, e.g., Edwin & Roberts 1982, 1983). In contrast, the modes in the other sub-spectrum are characterized by a non-vanishing ℑΩ, and the mode functions are oscillatory. A rich variety of modes have been reported in this sub-spectrum, the PFLK and trig modes being examples (e.g., Fig. 1 in Cally 2003). We concentrate on the trig modes, a feature of which is that the elements in the set {ℜΩ} for small axial wave numbers all correspond to short periodicities and are roughly equally spaced (Cally 1986; Terradas et al. 2005). However, we follow Wang et al. (2023) in calling them discrete leaky modes (DLMs). We choose not to use the customary term “leaky modes” in order to avoid confusions; in the literature, this term may include PFLKs (e.g., Cally 2003; Ruderman & Roberts 2006b) or may actually refer to the entire ℑΩ ≠ 0 sub-spectrum (Cally 1986; also Zajtsev & Stepanov 1975; Spruit 1982).

Our study will address kink motions in slab equilibria from an IVP perspective, paying special attention to the relevance of the DLM expectations to our time-dependent solutions. We stress that the slab geometry, although in less frequent use, may be relevant for interpreting the oscillatory behavior measured in, for example, post-flare supra-arcades (Verwichte et al. 2005), streamer stalks (e.g, Chen et al. 2010; Feng et al. 2011; Decraemer et al. 2020), and flare current sheets (e.g., Jelínek & Karlický 2012; Karlický et al. 2013; Yu et al. 2016). More importantly, key to our study is the wave behavior far from wave-guiding inhomogeneities, meaning that our qualitative findings are applicable to more general geometries. That said, DLMs in slab geometry are indeed controversial. The classic mode analysis in Roberts (2019, Sect. 5.6) indicates that no kink DLM exists that solves the classic BVP when the slab is density-enhanced relatively to its ambient corona, in contrast to earlier mode analyses that suggest the opposite (e.g., Terradas et al. 2005; Yu et al. 2015; Chen et al. 2018). Accepting DLMs as mathematical solutions, one further subtlety is that DLMs may have no physical meaning as argued on energetics grounds by Goedbloed et al. (2023). We note that this argument holds regardless of geometries and, hence, apparently contradicts direct numerical simulations that quantitatively display some expected DLM signatures (e.g., Terradas et al. 2005, 2007; Shi et al. 2023). Our study is partly motivated by these inconsistencies.

Our essentially analytical approach builds on the standard spectral theory of differential operators (e.g., Richtmyer 1978) or more intuitively on the Fourier-integral-based method (e.g., Whitham 1974). What is key is that one formulates an eigenvalue problem (EVP) consistent with the IVP in both the spatial operator and the BCs. Given localized initial exciters, what happens at infinity should not affect the perturbations at any finite distance. Consequently, the BCs at infinity in both our IVP and EVP are only nominal; no definitive requirement is specified. Formally speaking, this condition is the only difference between our EVP and the pertinent classic BVP. This difference turns out to be crucial; to clarify this, we use “eigen” in any term that arises in our EVP whenever possible. Eigenmode refers to a nontrivial solution jointly characterized by an eigenfrequency (ω) and an eigenfunction. Eigenspectrum refers to the collection of ω. We emphasize eigen because the nontrivial solutions qualify as such, the defining features being that all eigenfrequencies are real-valued and the set of eigenfunctions is orthogonal and complete. The absence of a definitive BC at infinity introduces a continuous eigenfrequency distribution in situations where one expects DLMs in classic BVPs. These details notwithstanding, the wave field at any instant can be decomposed into the eigenfunctions, any coefficient depending only on time as a simple harmonic oscillator. The initial conditions (ICs) for the IVP then fully determine all coefficients and, hence, a time-dependent solution. We choose to call this approach the method of eigenfunction expansion to avoid confusion with the Fourier-decomposition in classic mode analyses. Although well established, this approach has only occasionally been adopted in solar contexts (Oliver et al. 2014, 2015; Li et al. 2022; Wang et al. 2023; see also Cally 1991; Soler & Terradas 2015; Ebrahimi et al. 2020 for related studies). The outlined steps were considered self-evident, but they are outlined here to clarify matters.

This manuscript examines the response of symmetric slab setups to localized kink exciters. We allow for continuous nonuniformities, but restrict ourselves to 2D motions to make the Alfvén continuum irrelevant (see Goossens et al. 2011 for reasons). We largely focus on conceptual understandings behind how the periodicities in the t-dependent solutions connect to classic mode analysis and what role may be played for these periodicities by the equilibrium quantities and initial exciters. Our study is new in the following aspects, relatively to the most relevant studies in the literature. Firstly, an analytical study dedicated to our purposes with our approach is not available. Similar questions were addressed by Ruderman & Roberts (2006b) and Andries & Goossens (2007) for step profiles. However, the Laplace transform approach therein is quite involved, and such complications as the multi-valued-ness of the Green function may potentially obscure the relevance of DLMs. On the other hand, our approach bears close resemblance to a few available studies, which nonetheless also adopted step profiles and were devoted to different contexts (e.g., Oliver et al. 2014; Li et al. 2022). By addressing continuous profiles, our study will not only demonstrate the generality of the eigenfunction expansion approach, it will also shed more light on how the profile steepness impacts the system evolution. Secondly, a considerable number of numerical studies exist in various contexts that partially overlap with what we examine here (e.g., Terradas et al. 2005, 2007; Nakariakov et al. 2012; Guo et al. 2016; Lim et al. 2020, to name only a few). However, the numerical schemes therein are exclusively grid-based, and hence they do not explicitly involve such concepts as modes or eigenmodes. Our EVP-based approach makes it easier to explore the connection between a t-dependent solution and modal expectations.

This manuscript is structured as follows. Section 2 presents the steps that eventually lead to a 1D IVP, for which the analytical solution is then formulated in Sect. 3 with the eigenfunction expansion method. Section 4 presents a rather systematic set of example solutions, thereby helping to better visualize our conceptual understandings. Our study is summarized in Sect. 5, where some concluding remarks are also offered.

2. Problem formulation

2.1. Overall description

We adopt pressureless, ideal, MHD throughout, in which the primitive dependents are mass density ρ, velocity v, and magnetic field B. Let (x, y, z) denote a Cartesian coordinate system, and let subscript 0 denote the equilibrium quantities. We consider only static equilibria (v₀ = 0), assuming that the equilibrium magnetic field is uniform and z-directed (B₀ = B₀e_z). The equilibrium density ρ₀ is assumed to be a function of x only, and it is symmetric with regard to x = 0. Let ρ_i and ρ_e denote the internal and external densities, respectively. We model density-enhanced coronal slabs (ρ_i > ρ_e) of half-width d by prescribing

$$\begin{array}{c}\hfill {\rho }_0(x)={\rho }_{\mathrm{i}}+({\rho }_{\mathrm{e}}-{\rho }_{\mathrm{i}})f(x),\end{array}$$(1)

where the function f(x) attains zero at the slab axis (x = 0), but unity outside the slab (|x|> d). The Alfvén speed is defined by v_A²(x) = B₀²/μ₀ρ₀(x), with μ₀ being the magnetic permeability of free space. By v_Ai(v_Ae), we mean the internal (external) Alfvén speed evaluated with ρ_i (ρ_e). We further place two bounding planes at z = 0 and z = L to mimic magnetically closed structures anchored in the dense photosphere.

We examine axially standing linear kink motions in a strictly 2D fashion. Let subscript 1 denote small-amplitude perturbations. By 2D we mean ∂/∂y = 0, and we end up with an initial-boundary-value problem (IBVP) involving only v_1x, B_1x, and B_1z. The equilibria are perturbed in velocity only as implemented by the following ICs:

$$\begin{array}{cc}\hfill v_{1x}(x,z,t=0)& =u(x)sin(kz),\hfill \end{array}$$(2a)

$$\begin{array}{cc}\hfill B_{1x}(x,z,t=0)& =B_{1z}(x,z,t=0)=0.\hfill \end{array}$$(2b)

Axially standing motions are ensured by prescribing quantized axial wave numbers k = nπ/L (n = 1, 2, ⋯), and we arbitrarily focus on axial fundamentals (n = 1). Kink motions, on the other hand, are ensured by assuming the function u(x) in Eq. (2) to be even. It then follows that the 2D IBVP can be examined in the domain where 0 ≤ x < ∞ and 0 ≤ z ≤ L. We require that v_1x, ∂B_1x/∂z, and B_1z vanish at the lower (z = 0) and upper boundaries (z = L). All dependents are taken to be bounded at x → ∞, while the BCs at the slab axis (x = 0) write

$$\begin{array}{c}\hfill \frac{\partial v_{1x}}{\partial x}=\frac{\partial B_{1x}}{\partial x}=B_{1z}=0.\end{array}$$(3)

2.2. Formulation of the 1D initial-value problem

Our 2D IBVP simplifies to a simpler 1D version given the ICs and BCs. Formally, one may adopt the following ansatz:

$$\begin{array}{cc}\hfill v_{1x}(x,z,t)& =\widehat{v}(x,t)sin(kz),\hfill \end{array}$$(4a)

$$\begin{array}{cc}\hfill B_{1x}(x,z,t)& ={\widehat{B}}_x(x,t)cos(kz),\hfill \end{array}$$(4b)

$$\begin{array}{cc}\hfill B_{1z}(x,z,t)& ={\widehat{B}}_z(x,t)sin(kz),\hfill \end{array}$$(4c)

such that a single equation results for $\widehat{v}$:

$$\begin{array}{cc}\hfill \frac{{\partial }^2\widehat{v}}{\partial t^2}& ={v_{\mathrm{A}}}^2(x)(\frac{{\partial }^2\widehat{v}}{\partial x^2}-k^2\widehat{v}).\hfill \end{array}$$(5)

The ICs are written as

$$\begin{array}{c}\hfill \widehat{v}(x,t=0)=u(x),\frac{\partial \widehat{v}}{\partial t}(x,t=0)=0,\end{array}$$(6)

which follow from Eq. (2). We note that the 1D problem is in fact an IVP, despite being defined on 0 ≤ x < ∞ and subjected to the following BCs:

$$\begin{array}{c}\hfill \frac{\partial \widehat{v}}{\partial x}(x=0,t)=0,\widehat{v}(x\to \infty ,t)<\infty .\end{array}$$(7)

The BC at the slab axis (x = 0) follows from parity considerations, while the BC at infinity is only nominal.

2.3. Parameter specification

We now specify the necessary parameters. To start, we follow Li et al. (2018) to prescribe an “inner-μ” profile for the density distribution by taking f(x) in Eq. (1) to be

$$\begin{array}{c}\hfill f(x)=\{\begin{array}{cc}{(x/d)}^{\mu },\hfill & 0\le x\le d,\hfill \\ 1,\hfill & x>d,\hfill \end{array}\end{array}$$(8)

where μ > 0. Evidently, the density profile becomes steeper with μ, with the much-studied step distribution recovered for μ = ∞. We further specify the initial exciter in Eq. (2) as³

$$\begin{array}{c}\hfill \frac{u(x)}{v_{\mathrm{Ai}}}=\{\begin{array}{cc}{cos}^3(\pi x/2\mathrm{\Lambda }),\hfill & 0\le x\le \mathrm{\Lambda },\hfill \\ 0,\hfill & x\ge \mathrm{\Lambda }.\hfill \end{array}\end{array}$$(9)

This u(x) is rather arbitrarily chosen to be sufficiently smooth, with its spatial extent being characterized by the parameter Λ. Figure 1a illustrates our 2D IBVP by displaying the equilibrium density ρ₀ (the filled contours) and the initial velocity field (arrows). The x dependences of ρ₀ are shown in Fig. 1b for a number of steepness parameters, μ, as labeled, the density contrast being fixed at ρ_i/ρ_e = 5. Likewise, the function u(x) is displayed in Fig. 1c for an arbitrarily chosen set of Λ.

Fig. 1. Equilibrium setup and initial exciters. (a) Representation of 2D initial boundary value problem (IBVP). The x − z distribution of the equilibrium density is shown by the filled contours, superimposed on which is the initial velocity field (the arrows). Axial fundamentals are ensured by the z dependence of the initial perturbation. (b) Transverse profiles of the equilibrium density ρ0 as prescribed by Eq. (8). A number of steepness parameters μ are considered as labeled, whereas the density contrast is fixed at ρi/ρe = 5. (c) Transverse profiles of the initial perturbation u(x) as prescribed by Eq. (9). Labeled here are a number of values of Λ, which characterizes the spatial extent of the initial exciter.

Technical details aside, linear kink motions are fully dictated by two sets of parameters. The dimensional set is taken to be {ρ_i, d, v_Ai}, which merely provides normalizing constants. The behavior of linear motions then hinges on the dimensionless set:

$$\begin{array}{c}\hfill \{{\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}},\mu ;kd=\pi d/L;\mathrm{\Lambda }/d\}.\end{array}$$(10)

The effects of μ and Λ/d are of primary interest, and we also adjust the density contrast ρ_i/ρ_e when necessary. With AR loops in mind, we fix the axial wave number at kd = π/20. For axial fundamentals, this corresponds to an axial length of L = 20d, which lies toward the lower end of the range for AR loops typically imaged in the EUV (e.g., Aschwanden et al. 2004; Schrijver 2007). The density contrast ρ_i/ρ_e is taken to be in the range [2, 10], which is representative of AR loops as well (e.g., Aschwanden et al. 2004, and references therein).

3. EVP-based solutions: Formalism

This section solves the 1D IVP formulated in Sect. 2.2 by the method of eigenfunction expansion.

3.1. EVP-based solutions for generic transverse structuring

As described in the introduction, our approach is based on the general spectral theory of differential operators. In principle, there is no need for the initial exciter u(x) to follow Eq. (9) or for the equilibrium density ρ₀(x) to follow Eq. (8). This subsection specializes to the inner-μ prescription of ρ₀(x) for definitiveness; nonetheless, the steepness parameter μ is allowed to be arbitrary. The corresponding formulations generalize those in Sect. 3.1 of Wang et al. (2023), where we handled t-dependent sausage motions in coronal slabs with step profiles.

Key to our approach is the EVP that pertains to the spatial operator, ℒ, of the 1D IVP (see Eq. (5)). Specifically, nontrivial solutions are sought for the following equation:

$$\begin{array}{c}\hfill \mathcal{L}\breve{v}-{v_{\mathrm{A}}}^2(x)(\frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}\breve{v}-k^2\breve{v})={\omega }^2\breve{v},\end{array}$$(11)

which is defined on [0, ∞) and subject to the BCs (see Eq. (7))

$$\begin{array}{c}\hfill \mathrm{d}\breve{v}/\mathrm{d}x(x=0)=0\text{and}\breve{v}(x\to \infty )<\infty .\end{array}$$(12)

We note that k > 0 serves as an arbitrary parameter here.

One recognizes that the operator ℒ is self-adjoint under the definition of the scalar product

$$\begin{array}{c}\hfill ⟨U|V⟩{\displaystyle {\int }_0^{\infty }}{U}^{\ast }(x)V(x){\rho }_0(x)\mathrm{d}x,\end{array}$$(13)

where the asterisk represents a complex conjugate. Some properties then readily follow from general theory (see, in particular, Andries & Goossens 2007, for in-depth discussions on the self-adjoint aspect):

• All eigenvalues ω² are positive, and we see the eigenfrequency ω as positive without loss of generality.

• All eigenfunctions can be made and are seen as real-valued.

• The eigenspectrum, namely the collection of {ω}, comprises a point sub-spectrum and a continuum sub-spectrum.

• The eigenmodes in the point sub-spectrum are jointly characterized by an eigenfrequency ω_j and eigenfunction ${\breve{v}}_j$, with j = 1, 2, ⋯ being the mode label. The eigenfrequency ω_j is dictated by a DR, satisfying kv_Ai < ω_j < kv_Ae ≔ ω_crit. The eigenfunction ${\breve{v}}_j$ is evanescent for x > d, making the modes “proper” in that ${\breve{v}}_j$ is square integrable (i.e., $⟨{\breve{v}}_j|{\breve{v}}_j⟩$ converges in the classic sense).

• Proper eigenmodes are subject to cutoff wave numbers 0 < k_cut, 2 < k_cut, 3 < ⋯ with the exception of the lowest mode j = 1. By cutoff we mean that only one proper eigenmode arises when 0 < k < k_cut, 2, while a total of J proper eigenmodes are allowed when k_cut, J < k < k_{cut, J + 1} for J ≥ 2.

• The eigenmodes in the continuum subspectrum are jointly characterized by some eigenfrequency ω and eigenfunction ${\breve{v}}_{\omega }$. The eigenfrequency ω continuously spans (ω_crit, ∞), meaning the irrelevance of the concept of DRs. The eigenfunction ${\breve{v}}_{\omega }$ is oscillatory for x > d, making the modes “improper” in that ${\breve{v}}_{\omega }$ is not square integrable (i.e., $⟨{\breve{v}}_{\omega }|{\breve{v}}_{\omega }⟩$ is not defined in the classic sense).

• The eigensolutions in the proper set and the improper continuum are complete, satisfying the following orthogonality condition:

  $$\begin{array}{cc}\hfill ⟨{\breve{v}}_j|{\breve{v}}_{j^{\prime }}⟩& =⟨{\breve{v}}_j|{\breve{v}}_j⟩{\delta }_{j,j^{\prime }},\hfill \end{array}$$(14a)

  $$\begin{array}{cc}\hfill ⟨{\breve{v}}_j|{\breve{v}}_{\omega }⟩& =0,\hfill \end{array}$$(14b)

  $$\begin{array}{cc}\hfill ⟨{\breve{v}}_{\omega }|{\breve{v}}_{{\omega }^{\prime }}⟩& =q(\omega )\delta (\omega -{\omega }^{\prime }),\hfill \end{array}$$(14c)

  where δ_j, j′ is the Kronecker delta.

The EVP-based solution to our 1D IVP eventually gives

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \widehat{v}(x,t)=& {\displaystyle \sum _{j=1}^J}{c}_j{\breve{v}}_j(x)cos\left({\omega }_{j}t\right)+{\displaystyle {\int }_{k{v}_{\mathrm{Ae}}}^{\infty }}{S}_{\omega }{\breve{v}}_{\omega }(x)cos(\omega t)\mathrm{d}\omega ,\hfill \\ \hfill & 0\le x<\infty ,0\le t<\infty .\hfill \end{array}\end{array}$$(15)

We have made it explicit that the solution (15) is valid for arbitrary x and t. The summation on the right hand side (RHS) collects all (up to J) proper eigenmodes that may arise for a given k, while the integration accounts for the contribution from improper eigenmodes. The coefficients, c_j and S_ω, are given by

$$\begin{array}{c}\hfill c_j=\frac{⟨u|{\breve{v}}_j⟩}{⟨{\breve{v}}_j|{\breve{v}}_j⟩},S_{\omega }=\frac{⟨u|{\breve{v}}_{\omega }⟩}{q(\omega )}·\end{array}$$(16)

We refer to $S_{\omega }{\breve{v}}_{\omega }(x)$ as a “local spectral density”.

Some general remarks can be made with Eq. (15) regarding the possible relevance of DLMs. We recall that the concept of DRs only applies to proper eigensolutions in our EVP. However, the proper and improper eigenmodes are not determined beforehand. Rather, this distinction arises naturally to comply with the nominal BC at infinity, where no restriction is imposed. We further recall that the BC at infinity in “classic BVPs” is that no ingoing waves are allowed therein. Our Appendix B collects some necessary examinations on classic BVPs, and the concept of DRs is shown to always apply. Two types of nontrivial solutions arise, with the trapped modes being identical to our proper eigenmodes. DLMs, on the other hand, are demonstrated to exist for all the examined density profiles, thereby corroborating and generalizing available findings by, say, Goedbloed et al. (2023). We now consider the oscillation frequencies of DLMs {ℜΩ^DLM}. Given Eq. (15), there is clearly no guarantee for {ℜΩ^DLM} to stand out in a t-dependent solution; DLMs pertain to the classic BVP rather than our EVP. In this regard, we agree with Goedbloed et al. (2023) that DLMs do not have any physical meaning per se. However, we argue that the theories concerning DLMs may still prove seismologically relevant if such timescales as {ℜΩ^DLM} can be identified in some t-dependent solutions. Evidently, for DLMs to be relevant, the set {ℜΩ^DLM} needs to play a special role in the frequency dependence of S_ω in Eq. (15). Equally evident is that the details of the initial exciter must be crucial in determining whether this relevance materializes (see Eq. (16)).

Some further remarks are now necessary, given their immediate relevance to the inner-μ prescription in Eq. (8). We start with the following definitions:

$$\begin{array}{cc}\hfill {k_{\mathrm{i}}}^2& \frac{{\omega }^2-k^2{v_{\mathrm{Ai}}}^2}{{v_{\mathrm{Ai}}}^2}=\frac{{\omega }^2}{{v_{\mathrm{Ai}}}^2}-k^2,\hfill \end{array}$$(17a)

$$\begin{array}{cc}\hfill {k_{\mathrm{e}}}^2& \frac{{\omega }^2-k^2{v}_{\mathrm{Ae}}^2}{v_{\mathrm{Ae}}^2}=\frac{{\omega }^2}{v_{\mathrm{Ae}}^2}-k^2,\hfill \end{array}$$(17b)

$$\begin{array}{cc}\hfill {{\kappa }_{\mathrm{e}}}^2& -{k_{\mathrm{e}}}^2=k^2-\frac{{\omega }^2}{v_{\mathrm{Ae}}^2},\hfill \end{array}$$(17c)

$$\begin{array}{cc}\hfill D& {k_{\mathrm{i}}}^2+{{\kappa }_{\mathrm{e}}}^2=\frac{{\omega }^2}{{v_{\mathrm{Ai}}}^2}-\frac{{\omega }^2}{v_{\mathrm{Ae}}^2}·\hfill \end{array}$$(17d)

We note that k_i² and D are positive definite, and we see k_i as positive. Some intricacy arises for k_e². We chose to work with k_e² and see k_e as positive when handling improper eigenmodes (k_e² > 0). For proper eigenmodes (k_e² < 0), however, we always opt for κ_e² and see κ_e as positive. Regardless, the following dimensionless quantities are further defined:

$$\begin{array}{cc}\hfill {\overline{k}}_{\mathrm{i}}& k_{\mathrm{i}}d,{\overline{k}}_{\mathrm{e}}{k}_{\mathrm{e}}d,\hfill \end{array}$$(18a)

$$\begin{array}{cc}\hfill {\overline{\kappa }}_{\mathrm{e}}& {\kappa }_{\mathrm{e}}d,\overline{D}D{d}^2;\hfill \end{array}$$(18b)

these are used only when absolutely necessary. We now reformulate Eq. (11) into a Schrödinger form:

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2\breve{v}}{\mathrm{d}{x}^2}+[\frac{{\omega }^2}{{v_{\mathrm{A}}}^2(x)}-k^2]\breve{v}=0.\end{array}$$(19)

Some general results can be obtained for the uniform exterior (vA(x) = v_Ae for x > d). The external solution always gives

$$\begin{array}{c}\hfill {\breve{v}}_j(x)\propto {\mathrm{e}}^{-{\kappa }_{\mathrm{e}}x}\end{array}$$(20)

for proper eigenmodes, and it is always expressible as

$$\begin{array}{c}\hfill {\breve{v}}_{\omega }(x)=v_{\mathrm{Ai}}[A_{\mathrm{c}}cos(k_{\mathrm{e}}x)+A_{\mathrm{s}}sin(k_{\mathrm{e}}x)]\end{array}$$(21)

for improper eigenmodes. Here, A_c and A_s are some constants of integration. The function q(ω) in Eq. (16) becomes

$$\begin{array}{c}\hfill q(\omega )=({\rho }_{\mathrm{e}}{v_{\mathrm{Ai}}}^2)\frac{k_{\mathrm{e}}{v}_{\mathrm{Ae}}^2}{\omega }\frac{\pi ({A_{\mathrm{c}}}^2+{A_{\mathrm{s}}}^2)}{2}·\end{array}$$(22)

Equation (22) was first given in our previous sausage study (Wang et al. 2023, Eq. (34)), the derivation benefiting substantially from Appendix B of Oliver et al. (2015). Sausage and kink eigensolutions in our setup only formally differ in the interior (x < d). The reason for Eq. (22) involving only the external solution is that the integral over the interior is regular and hence does not contribute (see Eqs. (13) and (14)). We also note that Eqs. (21) and (22) apply to all steepness parameters (μ), the effects of which show up only indirectly via A_c and A_s.

We now consider the interior (x < d). We only examine a representative set of μ, for which compact analytical solutions exist for Eq. (19). Regardless, the following aspects always hold:

• An internal eigenfunction is found by solving Eq. (19) in conjunction with the BC at x = 0 (see Eq. (12)).

• The magnitude of an eigenfunction is irrelevant in principle. In practice, we nonetheless scale any internal improper eigenfunction ${\breve{v}}_{\omega }$ in such a way that it remains regular for vanishingly small k_e, as is the case when ω → ω_crit = kv_Ae.

• The internal solution is connected to the external one by demanding the continuity of both $\breve{v}$ and $\text{d}\breve{v}/\text{d}x$ at x = d.

• We always start with improper eigenmodes. Proper eigenmodes are described afterwards. The expression for a proper eigenfunction ${\breve{v}}_j$ is written to ensure the continuity of ${\breve{v}}_j$ itself, and the continuity of $\text{d}{\breve{v}}_j/\text{d}x$ then yields a DR. Expressions for $⟨{\breve{v}}_j|{\breve{v}}_j⟩$ are provided when available. We also analytically examine cutoff axial wave numbers, k_cut, j, but refrain from fully analyzing the DRs.

3.2. The case with μ = ∞

This subsection examines the much-studied step profile (μ = ∞). The internal solution to Eq. (19) is simply ∝cos(k_ix) given that the interior is uniform. An improper eigenfunction ${\breve{v}}_{\omega }(x)$ is expressed as

$$\begin{array}{c}\hfill \frac{{\breve{v}}_{\omega }(x)}{v_{\mathrm{Ai}}}=\{\begin{array}{cc}\frac{k_{\mathrm{e}}}{k_i}cos(k_{\mathrm{i}}x),\hfill & 0\le x\le d,\hfill \\ A_{\mathrm{c}}cos(k_{\mathrm{e}}x)+A_{\mathrm{s}}sin(k_{\mathrm{e}}x),\hfill & x>d,\hfill \end{array}\end{array}$$(23)

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill A_{\mathrm{c}}& =\frac{k_{\mathrm{e}}}{k_{\mathrm{i}}}cos(k_{\mathrm{i}}d)cos(k_{\mathrm{e}}d)+sin(k_{\mathrm{i}}d)sin(k_{\mathrm{e}}d),\hfill \\ \hfill A_{\mathrm{s}}& =\frac{k_{\mathrm{e}}}{k_{\mathrm{i}}}cos(k_{\mathrm{i}}d)sin(k_{\mathrm{e}}d)-sin(k_{\mathrm{i}}d)cos(k_{\mathrm{e}}d).\hfill \end{array}\end{array}$$(24)

We now consider proper eigenmodes. An eigenfunction is given as

$$\begin{array}{c}\hfill \frac{{\breve{v}}_j(x)}{v_{\mathrm{Ai}}}=\{\begin{array}{cc}{\mathrm{e}}^{-{\kappa }_{\mathrm{e}}d}cos(k_{\mathrm{i}}x),\hfill & 0\le x\le d,\hfill \\ cos(k_{\mathrm{i}}d){\mathrm{e}}^{-{\kappa }_{\mathrm{e}}x},\hfill & x>d.\hfill \end{array}\end{array}$$(25)

A rather concise expression is available for $⟨{\breve{v}}_j|{\breve{v}}_j⟩$, reading

$$\begin{array}{c}\hfill \frac{⟨{\breve{v}}_j|{\breve{v}}_j⟩}{{\rho }_{\mathrm{i}}{v_{\mathrm{Ai}}}^{2}d}=\frac{{\mathrm{e}}^{-2{\kappa }_{\mathrm{e}}d}}{2}[1+\frac{sin(2k_{\mathrm{i}}d)}{2k_{\mathrm{i}}d}+\frac{{cos}^2(k_{\mathrm{i}}d)}{{\kappa }_{\mathrm{e}}d}\frac{{\rho }_{\mathrm{e}}}{{\rho }_{\mathrm{i}}}]·\end{array}$$(26)

The DR governing the eigenfrequency ω_j is

$$\begin{array}{c}\hfill k_{\mathrm{i}}tan(k_{\mathrm{i}}d)={\kappa }_{\mathrm{e}}.\end{array}$$(27)

Cutoff wave numbers (k_cut, j) arise when the axial phase speed ω/k approaches v_Ae in response to a varying k. Now that κ_e → 0⁺, one recognizes from Eq. (27) that

$$\begin{array}{c}\hfill k_{\mathrm{cut},j}d=\frac{(j-1)\pi }{\sqrt{{\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}}-1}},(j=2,3,\cdots ).\end{array}$$(28)

The expressions for both the DR (Eq. (27)) and cutoff wave numbers (Eq. (28)) are standard textbook material (e.g., Roberts 2019, Chapter 5), albeit largely in the context of classic BVPs.

3.3. The case with μ = 2

This subsection addresses the case where μ = 2, for which purpose the following definitions are necessary:

$$\begin{array}{cc}\hfill p& \sqrt{\overline{D}}=\frac{\omega d}{v_{\mathrm{Ai}}}\sqrt{1-\frac{{\rho }_{\mathrm{e}}}{{\rho }_{\mathrm{i}}}},\hfill \end{array}$$(29a)

$$\begin{array}{cc}\hfill \alpha & \frac{1}{4}-\frac{{{\overline{k}}_{\mathrm{i}}}^2}{4p}=\frac{1}{4}-\frac{{(\omega d/v_{\mathrm{Ai}})}^2-{(kd)}^2}{4p},\hfill \end{array}$$(29b)

$$\begin{array}{cc}\hfill X& p{(x/d)}^2.\hfill \end{array}$$(29c)

The internal solution is always written e^−X/2M(α, 1/2, X) to respect the BC at x = 0, with M(⋅, ⋅, ⋅) being Kummer’s M function (e.g., DLMF 2016, Chapter 13). An improper eigenfunction is written

$$\begin{array}{c}\hfill \frac{{\breve{v}}_{\omega }(x)}{v_{\mathrm{Ai}}}=\{\begin{array}{cc}{\overline{k}}_{\mathrm{e}}{\mathrm{e}}^{-X/2}M(\alpha ,1/2,X),\hfill & 0\le x\le d,\hfill \\ A_{\mathrm{c}}cos(k_{\mathrm{e}}x)+A_{\mathrm{s}}sin(k_{\mathrm{e}}x),\hfill & x>d,\hfill \end{array}\end{array}$$(30)

where

$$\begin{array}{cc}\hfill A_{\mathrm{c}}& ={\overline{k}}_{\mathrm{e}}{\mathrm{e}}^{-X_1/2}M(\alpha ,1/2,X_1)cos(k_{\mathrm{e}}d)-Qsin(k_{\mathrm{e}}d),\hfill \end{array}$$(31a)

$$\begin{array}{cc}\hfill A_{\mathrm{s}}& ={\overline{k}}_{\mathrm{e}}{\mathrm{e}}^{-X_1/2}M(\alpha ,1/2,X_1)sin(k_{\mathrm{e}}d)+Qcos(k_{\mathrm{e}}d).\hfill \end{array}$$(31b)

Furthermore, X₁ = p is the value of X at x = d, and the symbol Q is defined by

$$\begin{array}{c}\hfill Q=-p{\mathrm{e}}^{-X_1/2}M(\alpha ,1/2,X_1)+4p\alpha {\mathrm{e}}^{-X_1/2}M(\alpha +1,3/2,X_1).\end{array}$$(32)

Here, M(α + 1, 3/2, X₁) appears as a result of the identify

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}\mathfrak{z}}M(a,b,\mathfrak{z})=\frac{a}{b}M(a+1,b+1,\mathfrak{z}).\end{array}$$(33)

We now consider proper eigenmodes. An eigenfunction is written as

$$\begin{array}{c}\hfill \frac{{\breve{v}}_j(x)}{v_{\mathrm{Ai}}}=\{\begin{array}{cc}{\mathrm{e}}^{-X/2}M(\alpha ,1/2,X){\mathrm{e}}^{-{\kappa }_{\mathrm{e}}d},\hfill & 0\le x\le d,\hfill \\ {\mathrm{e}}^{-X_1/2}M(\alpha ,1/2,X_1){\mathrm{e}}^{-{\kappa }_{\mathrm{e}}x},\hfill & x>d.\hfill \end{array}\end{array}$$(34)

A compact expression is not available for $⟨{\breve{v}}_j|{\breve{v}}_j⟩$. However, a DR may be readily derived, reading

$$\begin{array}{c}\hfill -{\overline{\kappa }}_{\mathrm{e}}=-p+4p\alpha \frac{M(\alpha +1,3/2,p)}{M(\alpha ,1/2,p)}·\end{array}$$(35)

The DR for trapped sausage modes in the same equilibrium was first given by Edwin & Roberts (1988). However, the DR for kink motions (Eq. (35)) is new as far as we are aware.

We now derive the approximate expressions for cutoff wave numbers k_cut, j. Interestingly, this derivation is connected to the behavior at large axial wave numbers (kd → ∞), in which case ω/k → v_Ai. One readily deduces that ${\overline{\kappa }}_{\text{e}}\to \infty$, p → ∞, and ${\overline{\kappa }}_{\text{e}}/p\to 1^{-}$, given their definitions (see Eqs. (17), (18), and (29)). The DR (35) then dictates that

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill E& \alpha \frac{M(\alpha +1,3/2,p)}{M(\alpha ,1/2,p)}\hfill \\ \hfill & =\alpha \frac{1+\frac{\alpha +1}{3/2}p+\frac{(\alpha +1)(\alpha +2)}{(3/2)(3/2+1)2!}{p}^2+\cdots }{1+\frac{\alpha }{1/2}p+\frac{(\alpha )(\alpha +1)}{(1/2)(1/2+1)2!}{p}^2+\cdots }\hfill \\ \hfill & \to 0,\hfill \end{array}\end{array}$$(36)

where the second equal sign follows from the definition of Kummer’s M function. Evidently, Eq. (36) is guaranteed provided α → 0, −1, −2, ⋯. One may then formally write

$$\begin{array}{c}\hfill {\alpha }_j\to 1-j,\text{for}kd\to \infty ,\end{array}$$(37)

where j = 1, 2, ⋯ is recalled to be the mode label. It turns out that Eq. (37), while derived for kd → ∞, holds approximately, even for relatively small values of kd. We note that α is evaluated as $(1-k_{\text{cut}}\sqrt{{\rho }_{\text{i}}/{\rho }_{\text{e}}-1})/4$ at cutoffs given that ω/k = v_Ae (see Eq. (29)). Approximating this α with Eq. (37) then yields that

$$\begin{array}{c}\hfill k_{\mathrm{cut},j}d\approx \frac{4j-3}{\sqrt{{\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}}-1}},(j=2,3,\cdots ),\end{array}$$(38)

where the absence of k_cut, 1 is accounted for. This approximate expression is increasingly accurate when j increases, overestimating the exact value by 16.6%, 8.3%, and 5.5% for j = 2, 3, and 4, respectively.

3.4. The case with μ = 1

This subsection addresses the case where μ = 1, for which the following definition is necessary:

$$\begin{array}{c}\hfill X\frac{-{{\overline{k}}_{\mathrm{i}}}^2+\overline{D}(x/d)}{{\overline{D}}^{2/3}}·\end{array}$$(39)

The internal solution is then a linear combination of Airy’s functions Ai and Bi, reading

$$\begin{array}{c}\hfill \breve{v}(x)\propto \frac{\mathrm{Ai}(X)}{{\mathrm{Ai}}^{\prime }(X_0)}-\frac{\mathrm{Bi}(X)}{{\mathrm{Bi}}^{\prime }(X_0)},\end{array}$$(40)

with Ai′ and Bi′ being Airy’s prime functions such that the BC at x = 0 is respected (e.g., DLMF 2016, Chapter 9). Furthermore,

$$\begin{array}{c}\hfill X_0=\frac{-{{\overline{k}}_{\mathrm{i}}}^2}{{\overline{D}}^{2/3}}\end{array}$$(41)

evaluates X at x = 0. An improper eigenfunction gives

$$\begin{array}{c}\hfill \frac{{\breve{v}}_{\omega }(x)}{v_{\mathrm{Ai}}}=\{\begin{array}{cc}{\overline{k}}_{\mathrm{e}}[\frac{\mathrm{Ai}(X)}{{\mathrm{Ai}}^{\prime }(X_0)}-\frac{\mathrm{Bi}(X)}{{\mathrm{Bi}}^{\prime }(X_0)}],\hfill & 0\le x\le d,\hfill \\ A_{\mathrm{c}}cos(k_{\mathrm{e}}x)+A_{\mathrm{s}}sin(k_{\mathrm{e}}x),\hfill & x>d,\hfill \end{array}\end{array}$$(42)

where

$$\begin{array}{cc}\hfill A_{\mathrm{c}}& ={\overline{k}}_{\mathrm{e}}{Y}_{1}cos(k_{\mathrm{e}}d)-{\overline{D}}^{1/3}{Y}_1^{\prime }sin(k_{\mathrm{e}}d),\hfill \end{array}$$(43a)

$$\begin{array}{cc}\hfill A_{\mathrm{s}}& ={\overline{k}}_{\mathrm{e}}{Y}_{1}sin(k_{\mathrm{e}}d)+{\overline{D}}^{1/3}{Y}_1^{\prime }cos(k_{\mathrm{e}}d).\hfill \end{array}$$(43b)

Here, Y₁ and Y₁′ are some coefficients given by

$$\begin{array}{cc}\hfill Y_0& =\frac{\mathrm{Ai}(X_0)}{{\mathrm{Ai}}^{\prime }(X_0)}-\frac{\mathrm{Bi}(X_0)}{{\mathrm{Bi}}^{\prime }(X_0)},\hfill \end{array}$$(44a)

$$\begin{array}{cc}\hfill Y_1& =\frac{\mathrm{Ai}(X_1)}{{\mathrm{Ai}}^{\prime }(X_0)}-\frac{\mathrm{Bi}(X_1)}{{\mathrm{Bi}}^{\prime }(X_0)},\hfill \end{array}$$(44b)

$$\begin{array}{cc}\hfill Y_1^{\prime }& =\frac{{\mathrm{Ai}}^{\prime }(X_1)}{{\mathrm{Ai}}^{\prime }(X_0)}-\frac{{\mathrm{Bi}}^{\prime }(X_1)}{{\mathrm{Bi}}^{\prime }(X_0)},\hfill \end{array}$$(44c)

with Y₀ defined for immediate future use, and

$$\begin{array}{c}\hfill X_1=\frac{{{\overline{\kappa }}_{\mathrm{e}}}^2}{{\overline{D}}^{2/3}}\end{array}$$(45)

being the value of X evaluated at x = d.

We now consider proper eigenmodes. An eigenfunction is

$$\begin{array}{c}\hfill \frac{{\breve{v}}_j(x)}{v_{\mathrm{Ai}}}=\{\begin{array}{cc}[\frac{\mathrm{Ai}(X)}{{\mathrm{Ai}}^{\prime }(X_0)}-\frac{\mathrm{Bi}(X)}{{\mathrm{Bi}}^{\prime }(X_0)}]{\mathrm{e}}^{-{\kappa }_{\mathrm{e}}d},\hfill & 0\le x\le d,\hfill \\ [12pt][\frac{\mathrm{Ai}(X_1)}{{\mathrm{Ai}}^{\prime }(X_0)}-\frac{\mathrm{Bi}(X_1)}{{\mathrm{Bi}}^{\prime }(X_0)}]{\mathrm{e}}^{-{\kappa }_{\mathrm{e}}x},\hfill & x>d.\hfill \end{array}\end{array}$$(46)

Some algebra leads to an expression for $⟨{\widehat{v}}_j|{\widehat{v}}_j⟩$, which reads

$$\begin{array}{cc}\hfill \frac{⟨{\breve{v}}_j|{\breve{v}}_j⟩}{{\rho }_{\mathrm{i}}{v_{\mathrm{Ai}}}^{2}d}=& \frac{{\mathrm{e}}^{-2{\overline{\kappa }}_{\mathrm{e}}}}{{\overline{D}}^{1/3}}[1-(1-\frac{{\rho }_{\mathrm{e}}}{{\rho }_{\mathrm{i}}})\frac{{{\overline{k}}_{\mathrm{i}}}^2}{\overline{D}}][X_1{Y}_1^2-{(Y_1^{\prime })}^2-X_0{Y_0}^2]\hfill \\ \hfill & -\frac{{\mathrm{e}}^{-2{\overline{\kappa }}_{\mathrm{e}}}}{{\overline{D}}^{2/3}}\frac{1-{\rho }_{\mathrm{e}}/{\rho }_{\mathrm{i}}}{3}[Y_1^{\prime }{Y}_1-X_1{(Y_1^{\prime })}^2+X_1^2{Y}_1^2-X_0^2{Y}_0^2]\hfill \\ \hfill & +Y_1^2\frac{{\rho }_{\mathrm{e}}}{{\rho }_{\mathrm{i}}}\frac{{\mathrm{e}}^{-2{\overline{\kappa }}_{\mathrm{e}}}}{2{\overline{\kappa }}_{\mathrm{e}}}·\hfill \end{array}$$(47)

When deriving Eq. (47), we used some identities for the indefinite integrals of Airy’s functions (see Sect. 9.11 (iv) in DLMF 2016). The DR for proper kink eigenmodes is further written as

$$\begin{array}{c}\hfill \frac{{\mathrm{Ai}}^{\prime }(X_1){\mathrm{Bi}}^{\prime }(X_0)-{\mathrm{Ai}}^{\prime }(X_0){\mathrm{Bi}}^{\prime }(X_1)}{\mathrm{Ai}(X_1){\mathrm{Bi}}^{\prime }(X_0)-{\mathrm{Ai}}^{\prime }(X_0)\mathrm{Bi}(X_1)}=-\frac{{\overline{\kappa }}_{\mathrm{e}}}{{\overline{D}}^{1/3}}·\end{array}$$(48)

We note that this DR is not available in the solar literature per se, despite the well-known applications of Airy’s functions in wave contexts (e.g., Snyder & Love 1983; Li et al. 2018).

We now derive the approximate expressions for cutoff wave numbers, k_cut, j, which arise when ω/k = v_Ae and, hence, ${\overline{\kappa }}_{\text{e}}=0$ (see Eqs. (17) and (18)). It follows from Eq. (45) that X₁ = 0, meaning that the DR (48) becomes

$$\begin{array}{c}\hfill \frac{{\mathrm{Ai}}^{\prime }(X_0)}{{\mathrm{Bi}}^{\prime }(X_0)}=\frac{{\mathrm{Ai}}^{\prime }(0)}{{\mathrm{Bi}}^{\prime }(0)}=-\frac{1}{\sqrt{3}}=-cot\frac{\pi }{3}·\end{array}$$(49)

Defining

$$\begin{array}{c}\hfill \zeta \frac{2}{3}{(-X_0)}^{3/2},\end{array}$$(50)

one finds that Ai′(X₀)/Bi′(X₀) is well approximated by −cot(ζ + π/4) (see Sect. 9.7 (ii) in DLMF 2016). We note that $\overline{D}={{\overline{k}}_{\text{i}}}^2={k_{\text{cut}}}^2({\rho }_{\text{i}}/{\rho }_{\text{e}}-1)$ at cutoffs (see Eqs. (17) and (18)). With X₀ now being $-{\overline{D}}^{1/3}$, Eq. (49) then leads to

$$\begin{array}{c}\hfill k_{\mathrm{cut},j}d\approx \frac{\frac{3}{2}(j-\frac{11}{12})\pi }{\sqrt{{\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}}-1}},(j=2,3,\cdots ),\end{array}$$(51)

where the absence of k_cut, 1 is also made explicit. This approximate expression is increasingly accurate when j increases, overestimating the exact value by merely 0.7%, even for j = 2.

3.5. Further remarks

This subsection presents some further remarks that help visualize the mathematical developments so far.

Fixing the density contrast at ρ_i/ρ_e = 5, Fig. 2 presents how the axial phase speed (ω/k) depends on the axial wave number (k) for the first several branches of proper kink eigenmodes. Different panels pertain to different values of the steepness parameter μ, with the two horizontal dashed lines representing the internal and external Alfvén speeds (v_Ai and v_Ae). We consistently distinguish between different branches with different line styles, and specifically illustrate our mode labeling convention for the step case (μ = ∞, Fig. 2a). Also annotated therein are cutoff wave numbers, which are labeled so as to be consistent with the mode labels. We note that cutoffs are absent for the first branch. We further note that the cutoff wave number k_cut, j tends to increase when μ decreases for a given ρ_i/ρ_e and a given label j ≥ 2. In fact, k_cut, 4 is so large when μ = 2 or μ = 1 that the j = 4 branch is outside the range for plotting Figs. 2b and 2c. This behavior of cutoffs k_cut, j can be adequately explained by their exact or approximate expressions (e.g., Eq. (28) for μ = ∞).

Fig. 2. Dependence of dimensionless axial phase speed (ω/kvAi) on dimensionless axial wave number (kd) for the first several branches of proper kink eigenmodes. The density contrast is fixed at ρi/ρe = 5, whereas a number of steepness parameters μ are examined in different panels. The two horizontal dashed lines represent the internal and external Alfvén speeds (vAi and vAe). Different branches are distinguished by the line styles as illustrated in panel a, where our convention for labeling the eigenmodes is also shown. Cutoff axial wave numbers arise for all branches except the first one (j = 1), their labels chosen to be consistent with our mode labeling convention. Panel a singles out the dimensionless axial wave number kd = π/20 to illustrate the spectrum for the EVP presented in Sect. 3.1. The continuum sub-spectrum is represented by the vertical arrow, the eigenfrequencies ω continuously extending from kvAe to infinity. Only one element is present in the point sub-spectrum for this chosen kd, and it is represented by the asterisk.

We now describe the implementation of the EVP-based solution (15). To start, we recall that the density contrast ρ_i/ρ_e is allowed to vary between 2 and 10, whereas the axial wave number will be fixed at kd = π/20. It follows from Eq. (28) that the inequality kd = π/20 < k_cut, 2 consistently holds for the step profile (μ = ∞), and therefore holds for other values of μ as well. Consequently, only one proper eigenmode is present in the point sub-spectrum associated with our EVP (see Sect. 3.1). Taking the step profile as an example, Fig. a then displays this proper eigenmode via an asterisk. Furthermore, the arrow is intended to represent the continuum sub-spectrum, with the ∞ symbol indicating that the eigenfrequency ω extends out to infinity. The following steps are then adopted to implement Eq. (15) for a given combination [ρ_i/ρ_e, kd, Λ/d]:

• Given a value of μ, we solve the corresponding DR for the eigenfrequency ω₁ of the only relevant proper eigenmode, and we then evaluate the associated eigenfunction ${\breve{v}}_1$. The inner products $⟨{\breve{v}}_1|{\breve{v}}_1⟩$ and $⟨u|{\breve{v}}_1⟩$ are evaluated afterwards, thereby enabling the evaluation of the coefficient c₁ (see Eq. (16)). The proper contribution is computed with the first term on the RHS of Eq. (15), taking J = 1.

• We evaluate the eigenfunctions ${\breve{v}}_{\omega }$ for the continuum eigenmodes, obtaining the coefficient A_c and A_s as byproducts. The function q(ω) is then evaluated with Eq. (22), making it straightforward to evaluate the coefficient S_ω with Eq. (16). The improper contribution is computed with the second term on the RHS of Eq. (15).

Despite being essentially analytical, the above steps nonetheless involve some numerical evaluations, with the integration over the continuum in Eq. (15) being an example. Convergence studies are therefore conducted to ensure that no difference can be discerned in our time-dependent solutions when, say, a different grid is employed to discretize the continuum.

4. EVP-based solutions: Specific examples

4.1. General spatio-temporal patterns

This subsection presents some generic features of our time-dependent solutions. We start with a computation pertaining to [ρ_i/ρ_e = 5, μ = ∞, kd = π/20, Λ/d = 1], the lateral speed for which is displayed as a function of x and t in Fig. 3. One readily sees that the system evolution comprises two qualitatively different stages. The short-time response (t ≲ 12d/v_Ai here) to the kink exciter features a series of oblique ridges in the slab exterior (x > d), signifying partial reflections and partial transmissions at the slab-ambient interface. By construction, these partial reflections and transmissions are associated with the improper contribution in view of the trigonometric x dependence of the improper eigenfunctions (${\breve{v}}_{\omega }$, Eq. (21)). The x − t diagram at large times is then characterized by a series of horizontal stripes, meaning that the external perturbations eventually become laterally standing. Evidently, this standing behavior derives from the dominance of the proper contribution (see Eq. (15) and note that only one proper eigenmode is involved in the summation). The qualitative behavior in Fig. 3 is common to our t-dependent solutions. Partial reflections or transmissions may not be readily recognizable for some combinations of ρ_i/ρ_e, μ, and Λ/d.

Fig. 3. Time-dependent solution obtained with eigenfunction expansion method for the combination [ρi/ρe = 5, μ = ∞, kd = π/20, Λ/d = 1]. Shown is the distribution of the lateral speed $\widehat{v}$ in the x − t plane, with the filled contours chosen to better visualize the wave dynamics at short times.

We focus on the time sequences of the lateral speed at the slab axis $\widehat{v}(x=0,t)$ hereafter. Adopting a fixed combination of [ρ_i/ρ_e = 5, kd = π/20, Λ/d = 1], Fig. 4 displays the $\widehat{v}(x=0,t)$ profiles for different steepness parameters, μ, in different panels as labeled. In addition to $\widehat{v}(x=0,t)$ itself (the solid curve in each panel), the contributions from proper and improper eigenmodes are further displayed by the dashed and dash-triple-dotted curves, respectively. Once again, one readily discerns the two-stage behavior. We first consider the proper contribution, namely a simple sinusoid with period P₁ = 2π/ω₁ (see Eq. (15)). We recall that ω₁ is only slightly smaller than ω_crit = kv_Ae (see Fig. 2), meaning that P₁ is only slightly longer than the axial Alfvén time 2π/ω_crit = 2L/v_Ae(≈17.9d/v_Ai here) for all the examined steepness parameters. We further recall that the proper contribution is stationary in magnitude (see Eq. (15)), meaning that any overall two-stage behavior in $\widehat{v}(x=0,t)$ is carried by the improper contribution. For any dash-triple-dotted curve, one then sees a transition from some rapid attenuation at short times to a later stage where the temporal attenuation is substantially slower. A comparison between different panels further indicates that this transition tends to occur later for larger μ. One also sees that some short periodicities on the order of the lateral Alfvén time (d/v_Ai) may be present in the initial stage (e.g., Fig. 4a), and these short periodicities tend to persist longer for large μ (compare, e.g., Fig. 4a with Fig. 4c). Regardless, the improper contribution at large times consistently features a longer periodicity that is only marginally below 2π/ω_crit = 2L/v_Ae.

Fig. 4. Temporal evolution of lateral speed at slab axis $\widehat{v}(x=0,t)$. A fixed combination [ρi/ρe = 5, kd = π/20, Λ/d = 1] is adopted, whereas a number of steepness parameters μ are examined in different panels as labeled. The solid curve in each panel shows the lateral speed itself. The contributions from proper and improper eigenmodes are presented by the dashed and dash-triple-dotted curves, respectively.

The temporal behavior in the improper contribution can be understood from two complementary perspectives. The first perspective, based on partial reflection and transmission, is visually more intuitive for interpreting the short periodicities (∼d/v_Ai) at early times (see also Fig. 3). Energetically, this perspective is also more intuitive in understanding the overall tendency for the improper contribution to diminish with time, because this tendency can only take place via the incessant transmission of fast perturbations into the ambient fluid. Our second perspective, built directly on Eq. (15), is to attribute all the temporal features to the interference among the continuum eigenmodes. For instance, the temporal attenuation is connected to the effect whereby any adjacent monochromatic components of the continuum will become increasingly out of phase as time proceeds. Physically, this effect is nothing but destructive interference; the increasingly out-of-phase cosines in Eq. (15) tend to cancel out. One then expects that the low-frequency portion with ω ≳ ω_crit = kv_Ae will eventually stand out, because it takes longer for the improper eigenmodes in this portion to become out of phase. This expectation is reproduced by Fig. 4, where any dash-triple-dotted curve is eventually characterized by slow attenuation and periodicity ≲2π/ω_crit. One also expects that the contribution from the high-frequency portion will attenuate more rapidly, which is indeed seen at early times. That this initial stage involves high-frequency improper eigenmodes is also corroborated by the short periodicities ∼d/v_Ai therein.

The interference perspective is applicable to broader solar contexts. We only consider studies where the eigenfunction expansion method is adopted. We further restrict ourselves to time intervals where some perturbation at a given position decays with time as a result of the destructive interference of continuum eigenmodes. Evidently, this interference perspective holds regardless of the origin of the continuum. We then discriminate two situations, in one of which the continuum results from the system being laterally open, and in the other the continuum arises because of continuous transverse structuring inside the domain. Our study is an example for the former. The physics is identical to what happens to wave motions in equilibrium setups with step transverse profiles for cylindrical (Oliver et al. 2014, 2015; Li et al. 2022) and slab geometries alike (Wang et al. 2023; see also Andries & Goossens 2007). What our study demonstrates is that a further physically relevant continuum does not necessarily arise when the transverse structuring is continuous. When indeed physically relevant, however, the pertinent continuum may play a vital role in explaining, say, the intricate dependence on the equilibrium quantities of damping envelopes of kink motions in coronal cylinders (Soler & Terradas 2015; see also Cally 1991). On this aspect, we note that the study by Soler & Terradas (2015) focused on the resonant damping in the Alfvén continuum, and hence discarded the continuum that is introduced by the system being unbounded. To our knowledge, the eigenfunction expansion method remains to be applied to the situation that physically involves continua of distinct origins.

4.2. Contribution from proper eigenmodes

This subsection focuses on the proper contribution, the purpose being to illustrate subtleties in assessing the capabilities for the examined equilibria to confine the energy imparted by the initial exciter. An intuitive indicator will be the dimensionless axial cutoff wave number k_cut, jd for any given label, j, and we take k_cut, 2d for ease of description. The intuition behind this is that a smaller k_cut, 2d represents a stronger confinement capability. We note that $k_{\text{cut},j}d\sqrt{{\rho }_{\text{i}}/{\rho }_{\text{e}}-1}$, for a given j, depends solely on the function f(x) for any equilibrium density distribution describable by Eq. (1), as explained by Li et al. (2018), even though proper sausage eigenmodes were examined therein. Focusing on our inner-μ profile (Eq. (8)), one recognizes that $k_{\mathrm{cut},2}d\sqrt{{\rho }_{\text{i}}/{\rho }_{\text{e}}-1}$ depends exclusively on the steepness parameter μ. It then follows that k_cut, 2d decreases with ρ_i/ρ_e or μ, at least for the examined steepness parameters (see Eqs. (28), (38), and (51)). This is indeed intuitive given that a slab deviates more strongly from its ambient fluid when ρ_i/ρ_e or μ increases, and hence one expects a stronger proper contribution inside the slab. We focus on the slab axis, and take $c_1{\breve{v}}_1(x=0)$ to illustrate that this expectation does not necessarily hold.

Figure 5 displays $c_1{\breve{v}}_1(x=0)$ as a function of the density contrast ρ_i/ρ_e for several Λ/d, as distinguished by the different line styles. The specific values of Λ/d are further labeled. In addition, the results for the step case (μ = ∞, the black curves) are compared with those for μ = 2 (red; see Fig. 5a) and μ = 1 (blue, Fig. 5b). Examining any curve, one sees that $c_1{\breve{v}}_1(x=0)$ consistently increases with ρ_i/ρ_e, as is expected for stronger structuring. The μ dependence, on the other hand, is more subtle if one inspects a pair of curves with the same line style but different colors. The calculation for a larger μ may indeed yield a larger $c_1{\breve{v}}_1(x=0)$ when the initial exciter is more spatially extended (e.g., the dash-triple-dotted curves in Fig. 5a labeled Λ/d = 2). However, the opposite may also occur as indicated by, say, the rightmost portions of the curves labeled Λ/d = 1 in Fig. 5a.

Fig. 5. Dependence of $c_1{\breve{v}}_1(x=0)$ on density contrast ρi/ρe for several values of Λ/d as labeled and highlighted by line styles. The results for the step case (μ = ∞; black curves) are compared with those for (a) μ = 2 (red) and (b) μ = 1 (blue). The axial wave number is fixed at kd = π/20. We note that $c_1{\breve{v}}_1(x=0)$ measures the magnitude of the asymptotic variation at the slab axis (see text for details).

One may wonder why the somehow counterintuitive μ dependence occurs occasionally. We address this by capitalizing on Eqs. (16) and (13) to rewrite $c_1{\breve{v}}_1(x=0)$ as

$$\begin{array}{c}\hfill c_1{\breve{v}}_1(x=0)={\displaystyle \underset{0}{\overset{\mathrm{\Lambda }}{\int }}}u(x)\frac{G(x)}{d}\mathrm{d}x,\end{array}$$(52)

where

$$\begin{array}{c}\hfill G(x)\frac{[{\rho }_0(x){\breve{v}}_1(x)]{\breve{v}}_1(x=0)}{{\displaystyle {\int }_0^{\infty }}{\rho }_0(x){\breve{v}}_1^2(x)\mathrm{d}x}d.\end{array}$$(53)

The introduction of the slab half-width d into Eqs. (52) and (53) is immaterial, the purpose being simply to make G(x) dimensionless. Defined this way, G(x) allows the proper eigenfunction to be arbitrarily scaled. One may therefore interpret G(x) in Eq. (52) as a response function of a system that is fully determined by [ρ_i/ρ_e, μ, kd]. The specific output $c_1{\breve{v}}_1(x=0)$ is then determined by how the input u(x) is distributed over the G(x) profile. We arbitrarily specialize this in a combination: [ρ_i/ρ_e = 10, kd = π/20]. The response function G(x) then depends solely on μ, and the associated profiles are plotted in Fig. 6 for two different values of μ, one being μ = ∞ (the black solid curves) and the other being μ = 2 (red). We note that Figs. 6a and 6b differ only in the range of the horizontal axis as far as the solid curves are concerned. The u(x) profile is additionally plotted by the dotted curve for (a) Λ/d = 1 and (b) Λ/d = 2, and we note that u(x) is independent of G(x). A comparison between the two G(x) profiles shows that they differ primarily in that for μ = 2, G(x) exceeds its μ = ∞ counterpart in the x ≲ 0.4d interval. The integral in Eq. (52) then means that this portion weighs more when the initial exciter u(x) is more spatially localized, making $c_v{\breve{v}}_1(x=0)$ larger for μ = 2 rather than for μ = ∞. This pertains to Fig. 6a, where Λ/d = 1. However, the x ≳ 0.4d portion plays a more important role when u(x) is sufficiently extended, thereby explaining why $c_v{\breve{v}}_1(x=0)$ for μ = ∞ is larger, as happens when Λ/d = 2 (see Fig. 6b). One may speculate that the subtle μ dependence in Fig. 5 only holds when u(x) takes the specific form of Eq. (9), which may indeed be true. However, our point is that the interpretation of G(x) as a response function provides a sufficiently general framework for understanding the system responses to different implementations of u(x). Furthermore, with the subtle μ dependence, we actually highlight the importance of the details of the initial exciter u(x) for determining how the system responds.

Fig. 6. Function G(x) plotted as solid curves over ranges (a) [0, d] and (b) [0, 2d] for two different steepness parameters μ, as highlighted by different colors. Plotted by the dotted curves are two initial perturbations u(x) with (a) Λ/d = 1 and (b) Λ/d = 2. A fixed combination [ρi/ρe = 10, kd = π/20] is employed (see text for details).

4.3. Contribution from improper eigenmodes

This subsection examines the improper contribution, again specializing to the slab axis (x = 0). Let ${\widehat{v}}_{\mathrm{improper}}(x,t)$ denote the integral in Eq. (15). We now focus specifically on the connection between the short periodicities at early times in ${\widehat{v}}_{\mathrm{improper}}(x=0,t)$ and the pertinent DLM expectations. An examination of DLMs is therefore necessary in the framework of classic BVPs, for which we recall that the concept of DRs always applies. We present such an examination in Appendix B, where the DRs for continuous profiles are new in solar contexts. A countable infinity of modes arises for any axial wave number kd when given a combination of [ρ_i/ρ_e, μ]. Let Ω_j be the frequency of the j-th mode, and see its real part ℜΩ_j as positive without loss of generality. By DLMs, we mean the infinite subset of modes with the defining features that ℜΩ_j > kv_Ae and ℑΩ_j < 0. In our study, all modes with j ≥ 2 qualify as DLMs for all [ρ_i/ρ_e, μ] given the chosen kd = π/20. We nonetheless denote any such frequency as ${\mathrm{\Omega }}_j^{\mathrm{DLM}}$ to emphasize the distinction between the oscillation frequency $\mathfrak{R}{\mathrm{\Omega }}_j^{\mathrm{DLM}}$ and any (real-valued) eigenfrequency ω in our improper continuum.

Some general expectations follow from Eq. (15). We start by reformulating $S_{\omega }{\breve{v}}_{\omega }(x=0)$ as

$$\begin{array}{c}\hfill S_{\omega }{\breve{v}}_{\omega }(x=0)={\displaystyle \underset{0}{\overset{\mathrm{\Lambda }}{\int }}}u(x)G_{\omega }(x)\mathrm{d}x,\end{array}$$(54)

where

$$\begin{array}{c}\hfill G_{\omega }(x)\frac{\left[{\rho }_0(x){\breve{v}}_{\omega }(x)\right]{\breve{v}}_{\omega }(x=0)}{q(\omega )};\end{array}$$(55)

we used Eqs. (13) and (16). Equation (54) closely resembles Eq. (52), meaning that G_ω(x) can also be interpreted as a response function. We note that neither G_ω(x) nor $\mathfrak{R}{\mathrm{\Omega }}_j^{\mathrm{DLM}}$ depends on u(x), and we suppose for now that the combination [ρ_i/ρ_e, μ, kd] is given. One then expects that some elements in the set $\{\mathfrak{R}{\mathrm{\Omega }}_j^{\mathrm{DLM}}\}$ are bound to appear for certain u(x), provided that this set is special regarding the ω dependence of G_ω(x). Whether this expectation holds depends on the details of u(x), which is solely characterized by the spatial extent Λ here.

Figure 7 takes a fixed [ρ_i/ρ_e = 5, μ = 1, kd = π/20] to illustrate what happens to ${\widehat{v}}_{\mathrm{improper}}(x=0,t)$ when Λ/d varies as labeled. Two features are evident. Firstly, some long periodicities of ∼2π/kv_Ae can always be identified at large times, and this tends to be more prominent for larger Λ. Secondly, and more importantly, short periodicities of ∼d/v_Ai do show up at early times provided that Λ/d is sufficiently small, one example being the solid line that pertains to Λ/d = 0.5. Figure 8c then address why these features arise by displaying the local spectral density $S_{\omega }{\breve{v}}_{\omega }(x=0)$ as a function of ω, distinguishing between different Λ/d with different line styles. The vertical dash-dotted line in magenta represents the critical frequency ω_crit = kv_Ae. We note that $S_{\omega }{\breve{v}}_{\omega }(x=0)$ approaches zero when ω approaches ω_crit from above, or equivalently when ${\overline{k}}_{\text{e}}\to 0^{+}$ (see Eq. (17)). As can be readily verified, the denominator q(ω) on the RHS of Eq. (55) scales as ${\overline{k}}_{\text{e}},$ whereas the numerator scales as ${{\overline{k}}_{\text{e}}}^2$ in this situation. Figure 8c then indicates that a peak always stands out in $S_{\omega }{\breve{v}}_{\omega }(x=0)$ (or equivalently its modulus) at a low frequency close to ω_crit, and this peak strengthens with Λ/d, thereby explaining both the persistence of the long periodicity and its Λ dependence. On the other hand, more peaks at frequencies substantially higher than ω_crit become increasingly prominent when Λ/d decreases. This naturally accounts for the behavior of the short periodicities in Fig. 7. On top of that, the high-frequency peaks increase close to the oscillation frequencies expected for the DLMs, among which the first two ($\mathfrak{R}{\mathrm{\Omega }}_2^{\mathrm{DLM}}$ and $\mathfrak{R}{\mathrm{\Omega }}_3^{\mathrm{DLM}}$) are marked by the magenta arrows. We now move on to the results for (a) μ = ∞ and (b) μ = 2. All qualitative features in Fig. 8c are seen to hold for these cases as well. Some quantitative differences arise nonetheless, one example being that the short periodicities are easier to develop for a larger μ (compare the solid curve in Fig. 8c with that in, say, Fig. 8a).

Fig. 7. Temporal evolution of improper contribution as given by the integral in Eq. (15). A fixed combination [ρi/ρe = 5, μ = 1, kd = π/20] is adopted, whereas a number of values are examined for Λ/d as labeled.

Fig. 8. Frequency dependences of local spectral density evaluated at the slab axis, $S_{\omega }{\breve{v}}_{\omega }(x=0)$. A fixed combination [ρi/ρe = 5, kd = π/20] is adopted, whereas a number of values are examined for the steepness parameter μ in different panels. For each μ, a number of values of Λ/d are tested and shown according to the line styles. The magenta arrows in each panel mark the oscillation frequencies of the first two discrete leaky modes ($\mathfrak{R}{\mathrm{\Omega }}_2^{\mathrm{DLM}}$ and $\mathfrak{R}{\mathrm{\Omega }}_3^{\mathrm{DLM}}$) computed for the same set of [ρi/ρe, μ, kd]. All vertical dash-dotted lines correspond to the critical frequency ωcrit = kvAe (see text for more details).

On may question what makes $\mathfrak{R}{\mathrm{\Omega }}_j^{\mathrm{DLM}}$ special in the improper contribution. In view of Eq. (54), the most likely reason for the systematic Λ dependence of the high-frequency peaks is that the set $\{\mathfrak{R}{\mathrm{\Omega }}_j^{\mathrm{DLM}},j\ge 2\}$ is special for the response function G_ω(x). However, why this happens remains puzzling because the DLMs pertain to a classic BVP, which observes a different BC from our IVP. It turns out that the function q(ω) plays a decisive role in determining the ω dependence of G_ω(x) (see Eq. (55)). We illustrate this point by considering the simplest case of μ = ∞, noting that similar arguments can be offered for other values of μ. To start, we note that the inequality $|{\mathrm{\Omega }}_j^{\mathrm{DLM}}{|}^2\gg k^2{v_{\text{Ae}}}^2$ consistently holds in our study. When μ = ∞, the DR pertinent to DLMs (Eq. (B.5)) dictates that

$$\begin{array}{c}\hfill \mathfrak{R}{\mathrm{\Omega }}_j^{\mathrm{DLM}}d/v_{\mathrm{Ai}}\approx (j-1)\pi ,j=2,3,\cdots ,\end{array}$$(56)

as was first shown for the limiting case kd → 0 by Terradas et al. (2005). We now consider the portion of the improper continuum that satisfies ω² ≫ k²v_Ae² > k²v_Ai² and, hence, k_i² ≈ ω²/v_Ai², k_e² ≈ ω²/v_Ae² (see Eq. (17)). The expressions for the coefficients A_c and A_s (see Eq. (24)) lead to

$$\begin{array}{c}\hfill {A_{\mathrm{c}}}^2+{A_{\mathrm{s}}}^2\approx \frac{{\rho }_{\mathrm{e}}}{{\rho }_{\mathrm{i}}}+(1-\frac{{\rho }_{\mathrm{e}}}{{\rho }_{\mathrm{i}}}){sin}^2\left(\frac{\omega d}{v_{\mathrm{Ai}}}\right)·\end{array}$$(57)

Evidently, the values of $\mathfrak{R}{\mathrm{\Omega }}_j^{\mathrm{DLM}}$ minimize A_c² + A_s² and hence tend to minimize q(ω) (see Eq. (22)), meaning that these values tend to stand out as extrema in the response function G_ω(x) and, hence, in the output $S_{\omega }{\breve{v}}_{\omega }(x=0)$ (see Eqs. (54) and (55)). Equation (57) is reminiscent of Eq. (83) in Andries & Goossens (2007), which addressed a similar IVP with the more involved Laplace transform approach. Equation (57) bears some even closer resemblance to Eq. (39) in Wang et al. (2023), where we examined the temporal responses to sausage exciters with the same eigenfunction expansion method.

At this point, we reiterate that IVP studies in general offer a fuller picture for the wave dynamics than classic mode analyses. As concrete examples, our Figs. 7 and 8 illustrate that the DLM expectations do not necessarily materialize. That said, the concept of DLMs may still prove useful as a shortcut approach for interpreting the complicated system evolution in appropriate spatio-temporal domains. The shortcut refers to the general fact that classic mode analyses are computationally much cheaper than IVP studies. Furthermore, the theoretical results for DLMs are independent of initial exciters, and hence they are much easier to implement than IVP studies from a seismological standpoint. It should therefore be informative to explore the requirements on the initial exciter (Λ/d here) for the concept of DLMs to make practical sense.

Figure 9 takes [ρ_i/ρ_e = 5, μ = 1, kd = π/20] to illustrate how we quantify the range of Λ/d in which the concept of DLMs is useful. Basically, we focus on the ω dependence of $|S_{\omega }{\breve{v}}_{\omega }(x=0)|$, and we systematically reduce Λ/d to track the variation of the position of the second peak (ω_peak, 2). We note that $|S_{\omega }{\breve{v}}_{\omega }(x=0)|$ always features a peak at some frequency ≳ω_crit = kv_Ae, and we consider this peak as the first one. Figure 9 shows a discontinuity around some sufficiently small Λ/d, which is to be denoted (Λ/d)₁ and reads ∼1.5 in this particular case. This behavior is actually common to all of our computations, and it results from our convention for numbering the peaks. Specifically, it always holds that the second peak initially moves toward higher frequencies when Λ/d decreases toward (Λ/d)₁. When Λ/d decreases further, however, an additional peak becomes identifiable at a lower frequency and therefore qualifies as our second peak. The value ω_peak, 2 then increases with decreasing Λ/d and eventually approaches $\mathfrak{R}{\mathrm{\Omega }}_2^{\mathrm{DLM}}$. Let (Λ/d)_crit denote the value of Λ/d where ω_peak, 2 equals $0.9\mathfrak{R}{\mathrm{\Omega }}_2^{\mathrm{DLM}}$. Equivalently, (Λ/d)_crit is where the Λ − ω_peak, 2 curve in Fig. 9 enters the hatched area, for which the lower edge is set to be $0.9\mathfrak{R}{\mathrm{\Omega }}_2^{\mathrm{DLM}}$ and the upper edge is arbitrarily taken to be $1.1\mathfrak{R}{\mathrm{\Omega }}_2^{\mathrm{DLM}}$. We deem the portion Λ/d ≤ (Λ/d)_crit to be where the concept of DLMs is useful. Here, the specific threshold factor 0.9 is not that important; taking (Λ/d)₁ as (Λ/d)_crit serves our purposes equally well. One may question the fact that the frequency range delineated by the hatched area is also crossed by the Λ − ω_peak, 2 curve in the portion Λ/d > (Λ/d)₁. The reason for discarding this large-Λ portion is that additional peaks emerge when Λ/d decreases from (Λ/d)₁. Once identifiable, the positions of these peaks become almost instantly close to the expectations for additional DLMs ($\mathfrak{R}{\mathrm{\Omega }}_j^{\mathrm{DLM}}$ with j ≥ 3). This can be readily seen by comparing, say, the solid and dashed curves in Fig. 8. In other words, the occasional proximity of ω_peak, 2 to $\mathfrak{R}{\mathrm{\Omega }}_2^{\mathrm{DLM}}$ for Λ/d > (Λ/d)₁ is actually irrelevant.

Fig. 9. Dependence of frequency ωpeak, 2 on Λ/d, where $|S_{\omega }{\breve{v}}_{\omega }(x=0)|$ attains the second peak. A fixed combination [ρi/ρe = 5, μ = 1, kd = π/20] is adopted. The hatched portion corresponds to the frequency range $[0.9,1.1]\times \mathfrak{R}{\mathrm{\Omega }}_2^{\mathrm{DLM}}$ (see text for details).

Figure 10 displays (Λ/d)_crit as a function of the density contrast ρ_i/ρ_e for a number of steepness parameters μ as labeled. The axial wave number is fixed at kd = π/20. Evidently, the curve for each μ serves as a dividing line in the ρ_i/ρ_e − Λ/d plane, with the concept of DLMs only being helpful in the portion below the curve. We now see (Λ/d)_crit as a function of ρ_i/ρ_e and μ. Figure 10 then indicates that (Λ/d)_crit increases monotonically with ρ_i/ρ_e or μ when the other parameter is fixed. Put together, this agrees with our intuition that the concept of DLMs may help interpret the system evolution for a broader range of initial exciters, provided that a slab makes a stronger distinction from its surroundings. However, (Λ/d)_crit is consistently smaller than ∼3.3 for density contrasts representative of AR loops. The DLM theories therefore only seem to be practically useful for situations where the exciters laterally span no more than several 10³ km, given typical widths of AR loops. We argue that this spatial extent is not that small when placed in the context of, say, energy-release processes in solar flares (see, e.g., the review by Shibata & Magara 2011). In our opinion, it is more of an issue that the short periodicities need to be resolved with a cadence higher than typically implemented for EUV instruments.

Fig. 10. Dividing curves in ρi/ρe − Λ/d plane, separating the area where the concept of discrete leaky modes helps (the portion below a curve) from where it does not (above). The axial wave number is fixed at kd = π/20, whereas a number of steepness parameters μ are examined as discriminated by the different colors.

5. Summary

This study was largely intended to address, in a physically transparent manner, how classic mode analyses connect to the time-dependent wave behavior in structured media. We chose to work in linear ideal MHD for simplicity, and we examined axial fundamental kink motions that ensue when localized velocity perturbations are introduced in some symmetric, straight, field-aligned slab equilibria. Only 2D motions were of interest, and the equilibria were taken to be structured only transversely. Continuous structuring was allowed for inside the slab, while the slab exterior was assumed to be uniform (Eqs. (1) and (8)). A 1D IVP was formulated for laterally open systems, with no definitive requirement imposed at infinity. An EVP was constructed correspondingly, enabling the IVP to be analytically solved in terms of the complete set of eigenfunctions. This t-dependent solution, given by Eq. (15), allows a clear distinction between the contribution from proper eigenmodes and that due to improper eigenmodes. The wave dynamics depends on two subgroups of parameters. One subgroup characterizes the initial exciter (the dimensionless spatial extent Λ/d here, Eq. (9)). The other subgroup characterizes the equilibria involving the density contrast ρ_i/ρ_e, profile steepness μ, and the half-width-to-length ratio d/L. A systematic set of example solutions was offered for parameters representative of AR loops.

Our results are summarized as follows; by “long” (“short”), we refer to periodicities on the order of the axial (lateral) Alfvén time. Overall, all spatio-temporal patterns consistently involve the improper contribution, which nonetheless attenuates with time such that the evolution at a given location is eventually dominated by the long periodicities carried by the proper contribution. Physically, this attenuation is attributed to the destructive interference among the improper continuum. Specializing in the slab axis, we demonstrate that the proper contribution is strengthened with the density contrast, but it may occasionally be stronger for less steep density profiles. We find that short periodicities can only be clearly identified in the improper contribution for sufficiently localized exciters. When identifiable, these periodicities tend to agree closely with the oscillation frequencies expected for the discrete leaky modes (DLMs) in classic analyses, despite the boundary conditions therein being different from those in our IVP. We demonstrate that the eigenfunction expansion approach allows the system response to be interpreted as the interplay between the initial exciter and some response function, the latter depending only on the equilibrium quantities (Eqs. (52) and (54)). All qualitative features can be explained as such, for proper and improper contributions alike.

Our results allow for some general remarks on the seismological applicability of the extensively studied DLMs. Above all, the examples in our Appendix B help clarify that DLMs are mathematically allowed as nontrivial solutions in classic mode analyses. Conceptually, this does not contradict the improper continuum in any example eigenspectrum in Fig. 2a; our eigensolutions observe different boundary conditions from those in classic analyses. It is therefore justifiable to say that modes that do not exist cannot be excited (Goedbloed et al. 2023, page 22), even though we prefer to take this as meaning the inadequacy for DLMs to account for the system evolution in the entire spatio-temporal volume. However, this inadequacy does not invalidate previous attempts that invoke the oscillation frequencies or even the exponential damping rates of DLMs to diagnose, say, the physical parameters of flaring loops in the context of flare QPPs (e.g., Kopylova et al. 2007; Nakariakov et al. 2012; Chen et al. 2015). Rather, the interference among the improper continuum eigenmodes can indeed make the DLM expectations visible. It is just that whether the DLM expectations materialize depends sensitively on the somehow intricate interplay between the initial exciter and the equilibrium setup. We take this intricacy as encouraging rather than discouraging; to illustrate this point, we note that DLMs are the only modes that are accepted to yield short periodicities in classic mode analyses. Supposing that a short periodicity is observed together with some long periodicity, which is admittedly rare but not impossible (see, e.g., Kolotkov et al. 2015, for a specific observation), the simultaneous use of both periodicities helps alleviate the non-uniqueness issue for inversion problems in coronal seismology (see Arregui & Goossens 2019 for more on this issue). If we now suppose that no short periodicity can be identified in some oscillatory signal measured with adequate temporal cadence, our Fig. 10 allows one to deduce the minimal lateral extent of the initial exciter, thereby complementing the customary practice that focuses on diagnosing the equilibrium quantities.