© The Authors 2026

1. Introduction

The past two decades have witnessed a rapid accumulation of observations on low-frequency waves and oscillations in the solar corona (see, e.g., Banerjee et al. 2007; De Moortel & Nakariakov 2012; Shen et al. 2022, for reviews). These observations have largely been placed in two contexts, one being the problem of coronal heating (see the reviews by, e.g., Parnell & De Moortel 2012; Arregui 2015; Van Doorsselaere et al. 2020) and the other being solar coronal seismology (SCS; see, e.g., Nakariakov & Kolotkov 2020; Nakariakov et al. 2024, for reviews). Whichever the application, theoretical investigations into linear magnetohydrodynamic (MHD) waves in structured media often prove indispensable. The transverse structuring is of particular interest given the filamentary nature of the observed wave hosts (see, e.g., Erdélyi & Goossens 2011; Nakariakov et al. 2022; Kolotkov et al. 2023, for topical issues).

We focus on how the transverse structuring influences wave dispersion, by which we specifically mean the distinction between the phase and group velocities. This distinction has long been of theoretical interest in SCS, and tends to be largely addressed for the response of some density-enhanced equilibria to impulsive localized exciters (see, e.g., the reviews by Nakariakov & Kolotkov 2020; Li et al. 2020). We consider one-dimensional (1D) equilibria, for which the density structuring is restricted to only one transverse direction. We let 2D motions refer to those that prohibit wave propagation in the extra transverse direction. It has been customary to examine 2D fast sausage motions to bring out the dispersive effects in slab or cylindrical geometry (e.g., Roberts et al. 1983, 1984; Edwin & Roberts 1988). Impulsively excited sausage motions were predicted to possess a telltale three-phase signature in their time sequences sampled sufficiently far from the exciters (Roberts et al. 1983). An almost monochromatic periodic phase appears first, followed by a stronger quasi-periodic phase and eventually by a monochromatic decay phase (see Edwin & Roberts 1986, for theoretical reasons). This three-phase behavior turns out to be compatible with 2D time-dependent simulations (e.g., Murawski & Roberts 1993, 1994). Likewise, the predicted periodicities were compatible with the observed pulsating behavior in radio bursts (e.g., Zhao et al. 1990; Fu et al. 1990) and visible forbidden line intensities (e.g., Pasachoff & Landman 1984; Pasachoff & Ladd 1987). Seismology was therefore enabled, yielding such key coronal parameters as the transverse Alfvén time (e.g., Roberts et al. 1984; Fu et al. 1990).

Considerable progress has been made for the dispersive evolution of impulsively excited waves in the new century. For instance, Nakariakov et al. (2004) showed that the three-phase signature in time sequences translates into “crazy tadpoles” in the associated Morlet spectra. Some subtleties notwithstanding, the expected temporal or Morlet features were shown to largely persist regardless of the details of the initial exciter (e.g., Goddard et al. 2019; Kolotkov et al. 2021) or the wave host (e.g., Nakariakov et al. 2005; Jelinek & Karlicky 2010; Pascoe et al. 2013; Yu et al. 2016, 2017; Guo et al. 2022; Shi et al. 2025, 2026). That these features tend to be robust is not surprising, given the robust applicability to large-time signals of the method of stationary phase (MSP; Edwin & Roberts 1986; Li et al. 2023). The scenario of impulsive waves has therefore been broadly invoked for or implicated in understanding the oscillatory behavior observed in the radio (e.g., Mészárosová et al. 2009; Kaneda et al. 2018) and visible passbands (e.g., Williams et al. 2001; Katsiyannis et al. 2003; Samanta et al. 2016). Likewise, instances in line with this scenario have been found in imaging data; two examples are the cyclic transverse displacements of streamer stalks (streamer waves; e.g., Chen et al. 2010; Kwon et al. 2013; Decraemer et al. 2020), and the quasi-periodic fast propagating waves (QFPs; Liu et al. 2010, 2011; Shen & Liu 2012).

This study is intended to offer an in-depth mode analysis of 3D kink motions in coronal slabs. Of particular interest are two questions, namely how wave dispersion is influenced by the inclusion of the third dimension, and the spatio-temporal features that can be expected for impulsive 3D kink motions. To our knowledge, the most relevant study is the one by Li et al. (2023, hereafter paper I), where we focus on demonstrating the peculiar propagation features of impulsively excited 3D kink motions via a representative time-dependent simulation. This manuscript is distinct from paper I in both scope and objective. First, mode analysis was performed in paper I, the purpose being primarily to support the interpretation of the specific time-dependent simulation results. Accordingly, the mode analysis was conducted largely in a numerical manner and was restricted to the transverse fundamental. In contrast, this manuscript is dedicated to a systematic investigation into the group velocity behavior of 3D kink modes trapped in a slab configuration. Numerical solutions to the pertinent dispersion relation are presented, and we make substantial analytical progress in multiple asymptotic regimes. Furthermore, both the transverse fundamental and its overtones are addressed. We note that transverse overtones have received little attention to date, despite the long-lasting interest in oblique kink modes (e.g., Ionson 1978; Wentzel 1979; Hollweg & Yang 1988). Second, we capitalized on our mode analysis to come up with a scheme for categorizing 3D kink modes from the group velocity perspective, thereby complementing their restoring-force-based characterization (e.g., Goossens et al. 2009; Bahari & Khalvandi 2017). Third, our theoretical results are connected to the dispersive evolution of impulsively excited 3D kink motions. We specifically examined what morphological features were expected before performing a time-dependent simulation, in contrast to paper I where the simulation results were interpreted largely a posteriori.

This manuscript is structured as follows. Section 2 formulates our problem, presenting the relevant dispersion relations (DRs) and collecting some necessary definitions. Section 3 then details the behavior of the group velocities of 3D kink modes by numerically solving the DR. Also presented are a set of approximate analytical solutions that prove valuable for understanding the numerical results. Section 4 places our findings in the context of the large-time features of impulsively excited kink motions. Our study is summarized in Sect. 5, where some concluding remarks are also offered.

2. Problem formulation

2.1. General formulation

We adopt ideal, gravity-free, pressureless MHD throughout, in which the primitive quantities are the mass density ρ, velocity v, and magnetic field B. We let the subscript 0 denote equilibrium quantities, and consider only static equilibria (v₀ = 0). We let (x, y, z) be a Cartesian coordinate system. The equilibrium magnetic field is taken to be uniform and z-directed (B₀ = B₀e_z). We further assume the equilibrium density ρ₀ to be a function of x only. The Alfvén speed is defined by v_A² = B₀²/(μ₀ρ₀), with μ₀ being the magnetic permeability of free space.

We let the subscript 1 denote small-amplitude perturbations. The linearized, time-dependent, ideal MHD equations write

$$\begin{array}{cc}\hfill {\rho }_0\frac{\partial {\mathit{v}}_1}{\partial t}& =\frac{(\mathrm{\nabla }\times {\mathit{B}}_1)\times {\mathit{B}}_0}{{\mu }_0}=-\frac{\mathrm{\nabla }({\mathit{B}}_0·{\mathit{B}}_1)}{{\mu }_0}+\frac{{\mathit{B}}_0·\mathrm{\nabla }{\mathit{B}}_1}{{\mu }_0},\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \frac{\partial {\mathit{B}}_1}{\partial t}& =\mathrm{\nabla }\times ({\mathit{v}}_1\times {\mathit{B}}_0).\hfill \end{array}$$(2)

Evidently, the right-hand side (RHS) of Eq. (1) represents the linearized Ampére force, the only restoring force (f) in pressureless MHD. By the second equal sign we show the nominal decomposition of f into the magnetic pressure gradient force (f^PG) and magnetic tension force (f^T). These further write

$$\begin{array}{cc}\hfill {\mathit{f}}^{\mathrm{PG}}& =f_x^{\mathrm{PG}}{\mathit{e}}_x+f_y^{\mathrm{PG}}{\mathit{e}}_y=-\frac{B_0}{{\mu }_0}(\frac{\partial B_{1z}}{\partial x}{\mathit{e}}_x+\frac{\partial B_{1z}}{\partial y}{\mathit{e}}_y),\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill {\mathit{f}}^{\mathrm{T}}& =f_x^{\mathrm{T}}{\mathit{e}}_x+f_y^{\mathrm{T}}{\mathit{e}}_y=\frac{B_0}{{\mu }_0}(\frac{\partial B_{1x}}{\partial z}{\mathit{e}}_x+\frac{\partial B_{1y}}{\partial z}{\mathit{e}}_y),\hfill \end{array}$$(4)

where the non-contributing field-aligned components are excluded (see the cylindrical study by Goossens et al. 2009).

We perform classic mode analysis for 3D motions hereafter, kicking off with the Fourier ansatz,

$$\begin{array}{c}\hfill g_1(x,y,z;t)=\mathfrak{R}\{\stackrel{\sim }{g}(x)exp[-ı(\omega t-k_{y}y-k_{z}z)]\},\end{array}$$(5)

where g₁ represents any linear perturbation, with k_y (k_z) being real-valued out-of-plane (axial) wavenumbers. We see the angular frequency ω as real-valued as well. In component form, Equations (1) and (2) then write

$$\begin{array}{cc}\hfill \omega {\stackrel{\sim }{v}}_x& =-\frac{B_0}{{\mu }_0{\rho }_0}(k_z{\stackrel{\sim }{B}}_x+ı{\stackrel{\sim }{B\prime }}_z),\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill \omega {\stackrel{\sim }{v}}_y& =-\frac{B_0}{{\mu }_0{\rho }_0}(k_z{\stackrel{\sim }{B}}_y-k_y{\stackrel{\sim }{B}}_z),\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill \omega {\stackrel{\sim }{B}}_x& =-B_0{k}_z{\stackrel{\sim }{v}}_x,\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill \omega {\stackrel{\sim }{B}}_y& =-B_0{k}_z{\stackrel{\sim }{v}}_y,\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill \omega {\stackrel{\sim }{B}}_z& =-ıB_0({\stackrel{\sim }{v}}_x^{\prime }+ık_y{\stackrel{\sim }{v}}_y),\hfill \end{array}$$(10)

where we adopt the shorthand notation ′ ≔ d/dx. By mode we refer to a nontrivial solution, which is jointly characterized by a mode frequency ω and an eigenvector $\{{\stackrel{\sim }{v}}_x,{\stackrel{\sim }{v}}_y,{\stackrel{\sim }{B}}_x,{\stackrel{\sim }{B}}_y,{\stackrel{\sim }{B}}_z\}$. We see the wavenumbers k_y and k_z as independents, taking equilibrium quantities to be parameters only.

2.2. Fast waves in unbounded uniform media

This subsection examines an unbounded uniform medium for future reference. Now that ρ₀ = const, Fourier decomposition is allowed for the x-direction as well, enabling the quantity $\stackrel{\sim }{g}$ in Eq. (5) to be expressible as

$$\begin{array}{c}\hfill \stackrel{\sim }{g}(x)=\breve{g}{\mathrm{e}}^{ık_{x}x}.\end{array}$$(11)

Here ğ is a constant. We focus on fast modes, or equivalently compressional Alfvén waves, by assuming k_y ≠ 0. Some explicit expressions for the restoring forces readily follow,

$$\begin{array}{c}\hfill {\breve{\mathit{f}}}^{\mathrm{PG}}=\frac{B_0^2}{{\mu }_0}\frac{{\omega }^2/v_{\mathrm{A}}^2-k_z^2}{ı\omega }{\breve{\mathit{v}}}_{\perp },{\breve{\mathit{f}}}^{\mathrm{T}}=\frac{B_0^2}{{\mu }_0}\frac{k_z^2}{ı\omega }{\breve{\mathit{v}}}_{\perp }\frac{{\breve{\mathit{f}}}^{\mathrm{PG}}}{\mathrm{\Lambda }},\end{array}$$(12)

where ${\breve{\mathit{v}}}_{\perp }={\breve{v}}_x{\mathit{e}}_x+{\breve{v}}_y{\mathit{e}}_y$. It then follows that the parameter

$$\begin{array}{c}\hfill \mathrm{\Lambda }=\frac{{\omega }^2}{k_z^2{v}_{\mathrm{A}}^2}-1\end{array}$$(13)

adequately quantifies how one force is related to the other.

Some general properties for fast modes can be readily deduced. To start, the DR reads

$$\begin{array}{c}\hfill {\omega }^2=(k_x^2+k_y^2+k_z^2)v_{\mathrm{A}}^2{\mathsf{k}}^2{v}_{\mathrm{A}}^2,\end{array}$$(14)

where a 3D wavevector k = k_xe_x + k_ye_y + k_ze_z is introduced. Plugging Eq. (14) into Eq. (13) yields that

$$\begin{array}{c}\hfill \mathrm{\Lambda }=\frac{k_x^2+k_y^2}{k_z^2}·\end{array}$$(15)

This means that the magnetic pressure gradient force and the tension force are always in-phase, thereby consistently complementing each other to drive fast motions. Some subtlety arises for near-parallel propagation (k_x² + k_y² ≪ k_z²), in which case the tension force dominates such that fast modes become nearly degenerate with shear Alfvén waves. We proceed to define the phase and group velocities as

$$\begin{array}{c}\hfill {\mathit{v}}_{\mathrm{ph}}\frac{\omega }{\mathsf{k}}{\mathit{e}}_{\mathsf{k}},{\mathit{v}}_{\mathrm{gr}}\frac{\partial \omega }{\partial \mathsf{k}}=\frac{\partial \omega }{\partial k_x}{\mathit{e}}_x+\frac{\partial \omega }{\partial k_y}{\mathit{e}}_y+\frac{\partial \omega }{\partial k_z}{\mathit{e}}_z,\end{array}$$(16)

where e_k = k/k. We assume {k_x, k_y, k_z, ω} to be non-negative without loss of generality. The DR (Eq. (14)) then dictates that

$$\begin{array}{c}\hfill {\mathit{v}}_{\mathrm{ph}}={\mathit{v}}_{\mathrm{gr}}=v_{\mathrm{A}}{\mathit{e}}_{\mathsf{k}},\end{array}$$(17)

namely fast modes are dispersionless. That fast modes are Alfvén-like for near-parallel propagation is reflected in the fact that v_gr is largely aligned with the equilibrium field B₀ = B₀e_z.

2.3. Kink modes in slab equilibria

This subsection proceeds to examine a slab equilibrium,

$$\begin{array}{c}\hfill {\rho }_0(x)=\{\begin{array}{cc}{\rho }_{\mathrm{i}},\hfill & |x|<d,\hfill \\ [0.2cm]{\rho }_{\mathrm{e}},\hfill & |x|>d,\hfill \end{array}\end{array}$$(18)

where d is the slab half-width. By internal (subscript i) and external (subscript e), we consistently refer to the equilibrium quantities inside and outside the slab, respectively. In particular, the internal (external) Alfvén speed, v_Ai (v_Ae), is evaluated with the internal (external) density. We note that a step profile (Eq. (18)) is adopted to avoid the resonant absorption of 3D kink motions in the Alfvén continuum (see Goossens et al. 2011 for conceptual clarifications). We note further that pressureless MHD does not allow slow motions per se. However, we avoid the single use of fast to describe the collective behavior of compressible motions (see Goossens et al. 2009 for more; see also Goossens et al. 2020, 2021). Rather, such terms as B₀-straddling or B₀-same-side are used to characterize 3D kink motions from the standpoint of group velocities.

We focus on trapped kink modes by supplementing Eqs. (6)–(10) with appropriate boundary conditions (BCs). Only the half volume x ≥ 0 needs to be considered for the resulting boundary value problem (BVP). The BC at the slab axis (x = 0) is specified as ${\stackrel{\sim }{v}}_x\prime ={\stackrel{\sim }{v}}_y={\stackrel{\sim }{B}}_x\prime ={\stackrel{\sim }{B}}_y={\stackrel{\sim }{B}}_z=0$, while all Fourier amplitudes are required to vanish when x → ∞. All mode frequencies are real-valued. Now that the BCs are homogeneous, some straightforward dimensional analysis yields that the mode frequencies can be formally written as

$$\begin{array}{c}\hfill \frac{{\omega }_{j}d}{v_{\mathrm{Ai}}}=W_j(k_{y}d,k_{z}d|{\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}})=W_j(k_{y}d,k_{z}d).\end{array}$$(19)

By the second equal sign we recall that the density contrast ρ_i/ρ_e is seen as known. The transverse order (j = 1, 2, …) numbers the mode frequencies at a given pair [k_y, k_z] by increasing order. We follow the convention that the transverse fundamental corresponds to j = 1, while its overtones correspond to j ≥ 2. We detail only the cases j = 1 and j = 2, noting that our analysis readily generalizes to any j. The subscript j is dropped for brevity hereafter, unless confusions may arise.

Some generic properties ensue without solving the BVP. One readily verifies that if ω is a mode frequency for a given pair [k_y, k_z], then so is −ω. Likewise, if ω is a mode frequency for [k_y, k_z], then it remains so for the pairs [ − k_y, k_z], [k_y, −k_z], and [ − k_y, −k_z]. It therefore suffices to see ω as positive, and consider only the quadrant k_y > 0, k_z > 0. Furthermore, Equations (6)–(10) allow the Fourier amplitudes of the restoring forces to be expressed as (see Eqs. (3) and (4))

$$\begin{array}{c}\hfill {\stackrel{\sim }{\mathit{f}}}^{\mathrm{PG}}=\frac{B_0^2}{{\mu }_0}\frac{{\omega }^2/v_{\mathrm{A}}^2-k_z^2}{ı\omega }{\stackrel{\sim }{\mathit{v}}}_{\perp },{\stackrel{\sim }{\mathit{f}}}^{\mathrm{T}}=\frac{B_0^2}{{\mu }_0}\frac{k_z^2}{ı\omega }{\stackrel{\sim }{\mathit{v}}}_{\perp }\frac{{\stackrel{\sim }{\mathit{f}}}^{\mathrm{PG}}}{\mathrm{\Lambda }},\end{array}$$(20)

where ${\stackrel{\sim }{\mathit{v}}}_{\perp }={\stackrel{\sim }{v}}_x{\mathit{e}}_x+{\stackrel{\sim }{v}}_y{\mathit{e}}_x$. Equation (20) is formally identical to its counterpart in uniform media (Eq. (12)), the key difference being that the Alfvén speed v_A and the quantity Λ need to be discriminated between the interior and exterior. We focus on the interior, supplementing Λ with the subscript i. Evidently,

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{\mathrm{i}}=\frac{{\omega }^2}{k_z^2{v}_{\mathrm{Ai}}^2}-1.\end{array}$$(21)

Some specific results follow from explicit solutions to the BVP. We proceed by defining

$$\begin{array}{cc}\hfill {\kappa }_{\mathrm{i},\mathrm{e}}^2& k_z^2-\frac{{\omega }^2}{v_{\mathrm{Ai},\mathrm{e}}^2},\hfill \end{array}$$(22)

$$\begin{array}{cc}\hfill m_{\mathrm{i},\mathrm{e}}^2& k_y^2+{\kappa }_{\mathrm{i},\mathrm{e}}^2=k_y^2+k_z^2-\frac{{\omega }^2}{v_{\mathrm{Ai},\mathrm{e}}^2},\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill n_{\mathrm{i}}^2& -m_{\mathrm{i}}^2=\frac{{\omega }^2}{v_{\mathrm{Ai}}^2}-(k_y^2+k_z^2).\hfill \end{array}$$(24)

We see both m_e² and m_e as positive. The external mode functions therefore write ∝e^−m_ex, meaning that m_e can be taken as a measure of the capability for slabs to trap kink motions. The signs of κ_i, e² and m_i² are unknown at this point. We see kink modes as belonging to the surface (body) subfamily when m_i² > 0 (−n_i² = m_i² < 0), taking m_i > 0 (n_i > 0) without loss of generality. Trapped 3D kink modes in our slab configuration are known to obey the DR (e.g., Arregui et al. 2007; Yu et al. 2021),

$$\begin{array}{c}\hfill coth(m_{\mathrm{i}}d)=\left(\frac{m_{\mathrm{e}}}{m_{\mathrm{i}}}\right)\left(\frac{-{\kappa }_{\mathrm{i}}^2}{{\kappa }_{\mathrm{e}}^2}\right),\end{array}$$(25)

or equivalently

$$\begin{array}{c}\hfill cot(n_{\mathrm{i}}d)=\left(\frac{m_{\mathrm{e}}}{n_{\mathrm{i}}}\right)\left(\frac{-{\kappa }_{\mathrm{i}}^2}{{\kappa }_{\mathrm{e}}^2}\right).\end{array}$$(26)

We choose to work with Eq. (25) (Eq. (26)) when handling surface (body) modes. Let k = k_ye_y + k_ze_z. This study pays special attention to both the phase and group velocities as defined by

$$\begin{array}{c}\hfill {\mathit{v}}_{\mathrm{ph}}\frac{\omega }{k}{\mathit{e}}_k,{\mathit{v}}_{\mathrm{gr}}=v_{\mathrm{gr},y}{\mathit{e}}_y+v_{\mathrm{gr},z}{\mathit{e}}_z\frac{\partial \omega }{\partial k_y}{\mathit{e}}_y+\frac{\partial \omega }{\partial k_z}{\mathit{e}}_z,\end{array}$$(27)

where e_k = k/k.

Of future use is the situation where k_y = 0. The DR of 2D kink modes is now textbook material, writing (see Sect. 5.5 in Roberts 2019, and references therein)

$$\begin{array}{c}\hfill cot(n_{\mathrm{i}}d)=\frac{n_{\mathrm{i}}}{m_{\mathrm{e}}}·\end{array}$$(28)

Mathematically, Equation (28) can be derived from its 3D counterpart by letting k_y = 0. Equation (26) is more appropriate than Eq. (25) for this purpose, because 2D kink modes belong exclusively to the body family. A series of cutoff axial wavenumbers k_z^(J) are well known to arise,

$$\begin{array}{c}\hfill k_z^{(J)}d=\frac{J\pi }{\sqrt{{\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}}-1}},\end{array}$$(29)

where J = 1, 2, … in the 2D context. Transverse fundamentals (j = 1) exist for any k_z > 0. Taking the first overtone (j = 2) as example, by cutoff we then mean that trapped modes are allowed only when k_z > k_z^{(j − 1)} = k_z⁽¹⁾ (e.g., Nakariakov & Roberts 1995; Li et al. 2013). The solutions to Eq. (28) are denoted by ω^2D = ω^2D(k_z) for clarity. Equation (29) is also involved in our examination of 3D kink modes, which nonetheless require that half integers (J = 1/2, 3/2, …) be addressed.

3. Results

This study is intended to demonstrate how classic mode analysis can be better placed in the context of impulsively excited waves. Regarding the DR (Eq. (25) or equivalently Eq. (26)), it should be ideal that one deduces those generic properties that are insensitive to the density contrast ρ_i/ρ_e. However, the transcendental nature of the DR does not allow a completely analytical treatment. We adopt a fixed ρ_i/ρ_e = 3, a value relevant for active region loops (e.g., Aschwanden et al. 2004), polar plumes (e.g., Wilhelm et al. 2011), and streamer stalks (e.g., Chen et al. 2011). We solve the DR with standard root-finders for this specific ρ_i/ρ_e = 3, and build a more generic picture by proceeding analytically in various limiting cases.

3.1. Generic characterization

This subsection characterizes kink modes by inspecting the form of the DR (Eqs. (25) and (26)). The mode frequencies ω for a given transverse order j then map to a dispersion sheet in the 3D space [k_y, k_z, ω]. We choose to present a dispersion sheet by showing some of its cuts through constant values of k_z. A dispersion curve then results, expressing ω as a function of k_y for a given k_z. One readily recognizes from Eq. (26) that the surfaces κ_i² = 0, κ_e² = 0, m_e = 0, and n_i = Jπ (J = 0, 1/2, 1, 3/2, …) are important in organizing the dispersion sheets; these surfaces are where the LHS or RHS of Eq. (26) changes sign. It follows that the intersections of these surfaces are likely to be important as well. In particular, the intersection between the surfaces n_i = Jπ and m_e = 0 satisfies

$$\begin{array}{c}\hfill (k_y^2+k_z^2)d^2=\frac{{(J\pi )}^2}{{\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}}-1},(J=1/2,1,3/2,\dots )\end{array}$$(30)

which defines a set of circles in the k_y − k_z plane. The reason for k_z^(J) (see Eq. (29)) to be relevant is then that k_z^(J) locates where these circles intersect k_y = 0. The quantity k_z^(J) evaluates to

$$\begin{array}{c}\hfill k_z^{(1/2)}d=1.11,k_z^{(1)}d=2.22,k_z^{(3/2)}d=3.33,\dots \end{array}$$(31)

for the chosen ρ_i/ρ_e = 3.

Figure 1 shows a series of k_y − ω planes for an increasing sequence of dimensionless axial wavenumbers (k_zd) as indicated. The associated mode frequency ω of the transverse fundamental (first overtone) is shown by the black solid (dash-dotted) curve for each k_z. A number of curves are further displayed for reference, with the color and linestyle being consistent across all panels. These are also labeled (see Fig. 1c1) to ease our description. Specifically,

Fig. 1. Characterization of trapped 3D kink modes in solar coronal slabs. A density contrast of ρi/ρe = 3 is chosen for illustrative purposes. Plotted in each panel are the dependences on the out-of-plane wavenumber (ky) of the mode frequency of the transverse fundamental (with transverse order j = 1, the black solid curve) and the first overtone (j = 2, black dash-dotted). The hatched portions in any ky − ω plane represent where trapped modes are forbidden, with four curves serving as the relevant borders (labeled in panel c1). Curves 0, 1, 20, and 3 correspond to where ${\kappa }_{\mathrm{i}}^2:=k_z^2-{\omega }^2/v_{\mathrm{Ai}}^2=0$, κe2 := kz2 − ω2/vAe2 = 0, ni2 := ω2/vAi2 − (ky2 + kz2) = 0, and me2 := ky2 + kz2 − ω2/vAe2 = 0, respectively. The curves labeled 2Jπ (J = 1/2, 1, 3/2) further correspond to where ni = Jπ. These curves pertain to the characterization of the j = 1 and j = 2 modes, whose dispersive behavior is typified by different panels where the dimensionless axial wavenumber (kzd) increases sequentially (see text for more details).

• Line 0, the blue solid line, corresponds to κ_i² = 0.

• Line 1, the blue dashed line, corresponds to κ_e² = 0.

• Curve 3, the upper red solid curve, corresponds to m_e = 0.

• Curve 2₀, the lower red solid curve, corresponds to n_i = 0.

• Curves 2_Jπ with J = [1/2, 1, 3/2], correspond to n_i = Jπ and are shown by the red dotted, dashed, and dash-dotted curves, respectively.

Above all, these curves are important for showing that trapped kink modes are allowed only in the non-hatched portion in a k_y − ω plane. The portion bounded by curve 3 from below does not allow trapped modes by definition, given that m_e² < 0 therein. That trapped modes are prohibited in the other two hatched portions is because there is always a mismatch between the signs of the LHS and RHS of Eq. (25).

The specific value of k_z does not qualitatively impact the transverse fundamental in the following aspects. First, the transverse fundamental is present for any k_z > 0 and k_y ≥ 0. Second, the fundamental always transitions from a body to a surface type as k_y increases, the dividing line being n_i = 0 (or equivalently m_i² = 0). Third, the dispersion curve of the fundamental always lies in the horizontal stripe bounded by lines 0 and 1, and is further bounded from above by curve 2_π/2 when the fundamental is of the body type. The reason is simply that the LHS of Eq. (25) or Eq. (26) needs to be positive definite.

The behavior of the first overtone may be qualitatively different for different values of k_z. This is so despite that the first overtone is always of the body type. Three regimes need to be discriminated regarding how the 2_Jπ (with J = 1/2, 1, 3/2) curves are positioned with respect to the stripe between lines 0 and 1.

• Regime I with 0 < k_z < k_z^(1/2) as typified by Fig. 1a. All three 2_Jπ curves lie above the horizontal stripe. One readily verifies that the dispersion curve starts from the intersection between curves 3 and 2_π/2, running always between curves 2_π and 2_π/2. Of relevance is then some cutoff value of k_y, only beyond which is the first overtone allowed.

• Regime II with k_z^(1/2) < k_z < k_z⁽¹⁾ as typified by Fig. 1b. Only the 2_π/2 curve extends into the horizontal stripe. One readily verifies that the dispersion curve starts from the intersection between curve 3 and line 1. The pair [k_y, ω] reads [0, k_zv_Ae] at this intersection, where the transverse overtone is not allowed per se because m_e = 0. The dispersion curve therefore lies entirely outside the horizontal stripe, running between curves 2_π and 2_π/2.

• Regime III with k_z > k_z⁽¹⁾. Figures 1c1 and 1c2 further typify what happens when k_z is below and above k_z^(3/2), respectively. Curve 2_π extends into the horizontal stripe in the former, while curve 2_3π/2 also does so in the latter. However, the dispersion curve is qualitatively the same in the following sense. One readily verifies that the intersection between curve 2_π and line 1 always lies on the dispersion curve, which therefore can be divided into two portions. The portion in the horizontal stripe is always bounded by curve 2_π from below, and never touches curve 2_3π/2 or line 1. On the other hand, the portion outside the horizontal stripe always runs between curves 2_π and 2_π/2.

All these properties follow from the inspection of the signs of the LHS and RHS of Eq. (26). Similar conclusions can therefore be inferred for overtones of higher transverse order (j). All overtones belong exclusively to the body type, even though some quantitative difference from the first overtone does arise for j ≥ 3 regarding how the associated dispersion curve is positioned in the k_y − ω plane. We take the second overtone j = 3 for instance, and by difference we mean that curves 2_Jπ with J = 3/2, 2, 5/2 are relevant simply because the cot function is π-periodic.

3.2. Transverse fundamental

This subsection is devoted to the transverse fundamental. Figure 2 displays the distribution in the k_y − k_z plane of the phase speed as equally spaced contours. Also plotted are the filled contours of the parameter n_i, to evaluate which we now see ω in Eq. (24) as the mode frequency. We note that n_i makes sense only for body modes (n_i² > 0), meaning that the portion not occupied by the color map corresponds to surface modes. Figure 3 further plots how the y-component (v_gr, y, black contours) and z-component (v_gr, z, red) of the group velocity depend on [k_y, k_z].

Fig. 2. Transverse fundamental kink mode in a coronal slab with a density contrast ρi/ρe = 3. The dependence on [ky, kz] of the phase speed $v_{\mathrm{ph}}=\omega /\sqrt{k_y^2+k_z^2}$ is shown by the equally spaced contours. Trapped modes are allowed for all [ky ≥ 0, kz > 0], with both body (ni2 > 0) and surface (ni2 < 0) modes being relevant. The color map represents the parameter ni, and is restricted to body modes.

Fig. 3. Transverse fundamental kink mode in a coronal slab with a density contrast ρi/ρe = 3. Shown are the dependences on [ky, kz] of the y-component (the black contours) and the z-component (red) of the group velocity. All contours are equally spaced.

Some analytical progress can be made to help digest Figs. 2 and 3. Consider near-parallel propagation first (k_y² ≪ k_z²), and start with the particular case k_y = 0. It is well documented for these 2D modes (e.g., Roberts 2019, Chapter 5) that the phase speed ω^2D(k_z)/k_z decreases monotonically with k_z, approaching v_Ae (v_Ai) when k_zd → 0 (k_zd → ∞). This behavior is what one sees along k_y = 0 in Fig. 2. One further deduces that

$$\begin{array}{cc}\hfill {\omega }^{2\mathrm{D}}(k_z)& \approx k_z{v}_{\mathrm{Ae}}\sqrt{1-{({\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}}-1)}^2{(k_{z}d)}^2},k_{z}d\ll 1,\hfill \end{array}$$(32)

$$\begin{array}{cc}\hfill {\omega }^{2\mathrm{D}}(k_z)& \approx k_z{v}_{\mathrm{Ai}}\sqrt{1+{(\pi /2)}^2{(k_{z}d)}^{-2}},k_{z}d\gg 1,\hfill \end{array}$$(33)

by specializing Eqs. (18) and (16) in Li et al. (2013) to pressureless MHD. We note that n_i evaluates to $\sqrt{{({\omega }^{2\mathrm{D}}/v_{\mathrm{Ai}})}^2-k_z^2}$ for 2D modes, and is itself a monotonically increasing function of k_z (see Fig. 2). Equations (32) and (33) then help show that n_id → 0 for k_zd → 0 and n_id → π/2 when k_zd → ∞.

We now consider the situation 0 < k_y²/k_z² ≪ 1. It follows from Li et al. (2023, Eq. (9)) that

$$\begin{array}{c}\hfill \omega (k_y,k_z)\approx {\omega }^{2\mathrm{D}}[1+\frac{1}{2}\frac{m_{\mathrm{e}}^{2\mathrm{D}}d-1}{{(k_{z}d)}^2+(m_{\mathrm{e}}^{2\mathrm{D}}d){({\omega }^{2\mathrm{D}}d/v_{\mathrm{Ai}})}^2}{(k_{y}d)}^2],\end{array}$$(34)

where ω^2D = ω^2D(k_z) and ${(m_{\mathrm{e}}^{2\mathrm{D}})}^2:=k_z^2-{({\omega }^{2\mathrm{D}})}^2/v_{\mathrm{Ae}}^2$. Equation (34) explicitly shows that the correction to ω^2D is only quadratic in k_y, meaning that v_gr, y = ∂ω/∂k_y → 0 when k_yd → 0 (see Fig. 3). More importantly, one notices from Eqs. (32) and (33) that $m_{\mathrm{e}}^{2\mathrm{D}}\propto k_z^2$ for k_zd ≪ 1 and $m_{\mathrm{e}}^{2\mathrm{D}}\propto k_z$ for k_zd ≫ 1. Overall, $m_{\mathrm{e}}^{2\mathrm{D}}$ turns out to be a monotonically increasing function of k_z. One then recognizes the existence of some critical axial wavenumber k_z^2D, c, which renders $m_{\mathrm{e}}^{2\mathrm{D}}-1$ negative (positive) when k_z < k_z^2D, c (k_z > k_z^2D, c). We let 𝒞 denote the coefficient in front of (k_yd)² in Eq. (34). It then follows that 𝒞 reverses its sign when k_zd passes through k_z^2D, cd. Somehow Fig. 2 indicates that both v_ph and n_i tend to decrease with k_y at any given k_z, demonstrating no signature of k_z^2D, cd. This behavior can be explained by Eqs. (27) and (24) where v_ph and n_i are defined; the direct k_y-dependences tend to dominate the indirect k_y-dependences through ω. However, k_z^2D, cd plays an important role for v_gr, y and v_gr, z. One sees that v_gr, y tends to decrease to negative values with increasing k_y for small k_z, whereas v_gr, y increases to positive values when k_y increases from zero for large k_z. This is a direct consequence of Eq. (34), which actually helps one to deduce k_z^2D, cd ≈ 1.38 from Fig. 3 by locating where v_gr, y changes sign for k_y → 0. The quantity k_z^2D, c is relevant for v_gr, z as well, to show which we focus on the range k_z < k_z^2D, c. One sees that v_gr, z tends to decrease (increase) with k_y when k_z is fixed at some small (large) value. This is readily understandable if one notices the relevance of d𝒞/dk_z when evaluating ∂ω/∂k_z with Eq. (34).

We now address the situation k_y²/k_z² ≫ 1 where sufficiently oblique modes can be safely taken as surface modes. It is known that (Tatsuno & Wakatani 1998; Yu et al. 2021)

$$\begin{array}{c}\hfill \omega (k_y,k_z)\approx k_{z}C(k_y),\text{with}\frac{C^2(k_y)}{v_{\mathrm{Ai}}^2}=\frac{1+tanh(k_{y}d)}{{\rho }_{\mathrm{e}}/{\rho }_{\mathrm{i}}+tanh(k_{y}d)}\end{array}$$(35)

approximately solves Eq. (25). The more restrictive situation k_yd → ∞ is well studied, yielding $C\to C_{\mathrm{kink}}=v_{\mathrm{Ai}}\sqrt{2/({\rho }_{\mathrm{e}}/{\rho }_{\mathrm{i}}+1)}$ (e.g., Ionson 1978; Goossens et al. 1992; Ruderman et al. 1995). One sees from Fig. 2 that the contours of the phase speed v_ph largely become radially directed in the k_y − k_z plane. This behavior can be explained by the simplest approximation ω ≈ k_zC_kink, which yields that $v_{\mathrm{ph}}\approx C_{\mathrm{kink}}/\sqrt{1+{(k_y/k_z)}^2}$. Furthermore, this simplest approximation also helps understand the overall tendencies for v_gr, y to be small and for v_gr, z to approach C_kink when k_y²/k_z² ≫ 1 (see Fig. 3). However, some subtlety remains regarding the k_y-dependence of v_gr, y or v_gr, z even if the improved approximation (Eq. (35)) is employed. For instance, one expects with Eq. (35) that v_gr, y approaches zero from below given that C(k_y) decreases monotonically with k_y. Figure 3, however, indicates that v_gr, y > 0 for sufficiently large k_y. Likewise, Eq. (35) suggests that v_gr, z should approach C_kink from above. This contradicts the behavior of v_gr, z at large k_y. It turns out that Eq. (35) needs to be improved to address the k_y-dependence for those ranges of k_y where tanh(k_yd) varies only slowly. Writing

$$\begin{array}{c}\hfill \omega \approx k_{z}C(k_y)(1-Q\frac{k_z^2}{k_y^2}),\end{array}$$(36)

we proceed to Taylor-expand all terms on both sides of the DR (Eq. (25)) to the lowest order in k_z²/k_y². The coefficient Q then results from some algebra, reading

$$\begin{array}{c}\hfill Q=\frac{1}{4}\frac{{(1-{\rho }_{\mathrm{e}}/{\rho }_{\mathrm{i}})}^2}{{[{\rho }_{\mathrm{e}}/{\rho }_{\mathrm{i}}+tanh(k_{y}d)]}^2}\{tanh(k_{y}d)-(k_{y}d)[1-tanh(k_{y}d)]\}.\end{array}$$(37)

For the chosen ρ_i/ρ_e = 3, Equation (36) proves to be remarkably accurate when k_y/k_z ≳ 2. Consequently, one readily explains the large-k_y behavior of v_gr, y or v_gr, z by noting that

$$\begin{array}{cc}\hfill Q& \approx \frac{1}{4}\frac{{(1-{\rho }_{\mathrm{e}}/{\rho }_{\mathrm{i}})}^2}{{(1+{\rho }_{\mathrm{e}}/{\rho }_{\mathrm{i}})}^2},\omega \approx k_z{C}_{\mathrm{kink}}(1-Q\frac{k_z^2}{k_y^2}),\hfill \end{array}$$(38a)

$$\begin{array}{cc}\hfill \frac{v_{\mathrm{gr},y}}{C_{\mathrm{kink}}}& \approx 2Q\frac{k_z^3}{k_y^3},\frac{v_{\mathrm{gr},z}}{C_{\mathrm{kink}}}\approx 1-3Q\frac{k_z^2}{k_y^2}\hfill \end{array}$$(38b)

for sufficiently large k_yd.

Figure 4 gathers the directional information of the group (v_gr) and phase velocities (v_ph), aiming to categorize transverse fundamental kink modes on the k_y − k_z plane. Plotted are the [k_y, k_z]-dependences of (a) the angle (∠(v_gr, e_z)) between v_gr and the z-direction, and (b) the angle (∠(v_gr, v_ph) between the two vector velocities themselves. The angles are measured counterclockwise from v_gr, the readings being in degrees. All contours are equally spaced. Also shown by the filled contours is the parameter Λ_i, which is recalled to be the ratio of the magnetic pressure gradient force to the magnetic tension force (see Eq. (21)).

Fig. 4. Transverse fundamental kink mode in a coronal slab with a density contrast ρi/ρe = 3. Shown are the dependences on [ky, kz] of (a) the angle between the group velocity vgr and the z-direction, and (b) the angle between the group velocity vgr and the phase velocity vph. The angles are in degrees and measured counterclockwise from vgr to ${\widehat{\mathit{e}}}_z$ or vph. Superimposed are the filled contours of Λi, the ratio of the magnetic pressure gradient force to the magnetic tension force. All contours are equally spaced. The angle information allows the fundamental to be classified into two regimes as labeled in panel c (see text for details).

We consider the quantity Λ_i first. By definition, Λ_i is positive definite for the transverse fundamental and overtones alike, given that all trapped modes ensure κ_i² = k_z² − ω²/v_Ai² < 0 (see Fig. 1). Figure 4a further indicates that Λ_i tends to be smaller than unity for the majority of the [k_y, k_z] combinations, with the portion Λ_i > 1 restricted to the lower left corner. This can be understood with the approximate expressions that were presented previously. The expressions for 2D modes (k_y = 0, Eqs. (32) and (33)) suggest that Λ_i → ρ_i/ρ_e − 1 when k_zd ≪ 1 and Λ_i → 0 when when k_zd → ∞. Addressing a small k_y²/k_z², Equation (34) further suggests that Λ_i decreases (increases) with k_y for relatively small (large) values of k_z that ensure a negative (positive) k_y²-correction. Likewise, one deduces from Eq. (36) that

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{\mathrm{i}}\to \frac{C_{\mathrm{kink}}^2}{v_{\mathrm{Ai}}^2}-1=\frac{1-{\rho }_{\mathrm{e}}/{\rho }_{\mathrm{i}}}{1+{\rho }_{\mathrm{e}}/{\rho }_{\mathrm{i}}}<1,\end{array}$$(39)

when k_y²/k_z² → ∞. Actually, one sees from Fig. 4a that Λ_i is nearly uniform for sufficiently large k_y; the subtle k_y− (and k_z-) dependences of C(k_y) and the Q term play only a marginal role in determining the gross distribution of Λ_i.

We now attempt to place the angle distributions in the context of Λ_i. For immediate future reference, we classify trapped 3D kink modes into the following three regimes, noting that v_gr, z is positive definite in our context.

• B₀-straddling, by which we mean the situation where v_gr, y < 0 (i.e., v_gr and v_ph straddle B₀).

• B₀-same-side A, by which we mean the situation where v_gr, y > 0 and |∠(v_gr, e_z)| < 15°.

• B₀-same-side F, by which we refer to the situation where v_gr, y > 0, |∠(v_gr, e_z|> 15°, and |∠(v_gr, v_ph)| < 15°.

We note that by B₀-same-side we mean v_gr and v_ph sit on the same side of B₀. This classification scheme largely comes from the analogy with 3D waves (k_y > 0) in unbounded uniform media (see Sect. 2.2), with the threshold angle arbitrarily chosen to be some small value (15° here). We note that the waves examined in Sect. 2.2, broadly called fast therein, need to be categorized into the B₀-same-side A and B₀-same-side F subgroups here. Equation (15) then offers an unambiguous association of this categorization with the force ratio; the B₀-same-side A subgroup arises when the force ratio Λ < tan²(15° ) ≈ 0.07, whereas the B₀-same-side F subgroup arises when the opposite is true. The letters A and F are meant to distinguish the Alfvén-like modes and the genuinely fast-like modes. The B₀-straddling subgroup is absent in Sect. 2.2, but bears substantial resemblance to slow waves in unbounded uniform media with non-vanishing puressue (e.g., Roberts 2019, Sect. 2.7). Figures 4a and 4b indicate that transverse fundamental kink modes qualify as B₀-same-side A in the majority of the k_y − k_z plane. This is true even for highly oblique propagation (k_y²/k_z² ≫ 1), despite that the force ratio Λ_i is only marginally small (1/2 for the chosen ρ_i/ρ_e = 3, see Eq. (39)). In this regard, our classification scheme agrees with the force-ratio-based argument by Goossens et al. (2009) in that transverse fundamental kink modes are largely Alfvén-like, despite that cylindrical geometry was examined therein. More striking is the B₀-straddling subgroup, which shows up at the lower left corner in the k_y − k_z plane. By striking we mean that pressureless MHD does not allow slow waves in textbook sense as far as unbounded uniform media are concerned. However, these motions do make it into Fig. 4, thereby highlighting the intricacies that transverse structuring brings into wave studies. For clarity, we collect our classification in Fig. 4c and label the subgroups accordingly. One sees that the force ratio Λ_i ≳ 1 in the B₀-straddling corner.

3.3. First transverse overtone

This subsection is devoted to the first transverse overtone. Figure 5 follows the format of Fig. 2 to present the [k_y, k_z]-dependences of the phase speed v_ph and the quantity n_i. We note that trapped 3D modes are forbidden within the circle $d\sqrt{k_y^2+k_z^2}=k_z^{(1/2)}d=(\pi /2)/\sqrt{{\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}}-1}$, which reads 1.11 given a ρ_i/ρ_e = 3. We note also that trapped 2D modes (k_y = 0) are further prohibited for k_z^(1/2) < k < k_z⁽¹⁾ = 2k_z^(1/2). When allowed, trapped modes belong exclusively to the body type (n_i² > 0), meaning that n_i makes sense. Figure 6 presents, in a form identical to Fig. 3, the distributions in the k_y − k_z plane of the y component (v_gr, y) and z-component (v_gr, z) of the group velocity.

Fig. 5. Similar to Fig. 2, but for the first transverse overtone. Trapped 3D modes are forbidden within the circle $d\sqrt{k_y^2+k_z^2}=k_z^{(1/2)}d=(\pi /2)/\sqrt{{\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}}-1}$, which reads 1.11 for the chosen ρi/ρe = 3. Trapped 2D modes are further prohibited for kz(1/2) < k < kz(1) = 2kz(1/2) along ky = 0. Trapped modes are exclusively of the body type (ni2 > 0).

It proves convenient to make analytical progress in various limiting cases to help understand Figs. 5 and 6. First consider highly oblique propagation k_y²/k_z² ≫ 1. One readily sees from the DR (Eq. (26)) that n_i → π⁻ when k_y²/k_z² → ∞, meaning that the leading order solution $\widehat{\omega }$ reads

Fig. 6. Similar to Fig. 3, but for the first transverse overtone.

$$\begin{array}{c}\hfill \widehat{\omega }(k_y,k_z)\approx \frac{v_{\mathrm{Ai}}}{d}\sqrt{{\pi }^2+{(k_{y}d)}^2+{(k_{z}d)}^2}.\end{array}$$(40)

It then follows that $v_{\mathrm{ph}}=\omega /\sqrt{k_y^2+k_z^2}$ is approximately proportional to $\sqrt{1+{\pi }^2/[{(k_{y}d)}^2+{(k_{z}d)}^2]}$, meaning that the contours of v_ph eventually become a series of circles centered at the origin on the k_y − k_z plane. This is what one sees from Fig. 5. Furthermore, Equation (40) indicates that

$$\begin{array}{c}\hfill \frac{v_{\mathrm{gr},y}}{v_{\mathrm{Ai}}}\approx \frac{1}{\sqrt{1+[{\pi }^2+{(k_{z}d)}^2]/{(k_{y}d)}^2}}·\end{array}$$(41)

At a given k_z, this means that v_gr, y approaches the internal Alfvén speed v_Ai as a monotonically increasing function of k_y, in close agreement with Fig. 6. Likewise, one sees from Eq. (40) that

$$\begin{array}{c}\hfill \frac{v_{\mathrm{gr},z}}{v_{\mathrm{Ai}}}\approx \frac{1}{\sqrt{1+[{\pi }^2+{(k_{y}d)}^2]/{(k_{z}d)}^2}},\end{array}$$(42)

which further simplifies to $v_{\mathrm{gr},z}\propto 1/\sqrt{1+{(k_y/k_z)}^2}$ when (k_yd)² ≫ π². This means that the contours of v_gr, z eventually become a series of rays radially directed from the origin on the k_y − k_z plane, thereby explaining the red contours for k_y²/k_z² ≫ 1 in Fig. 6. For completeness, we note that an improved version of Eq. (40) can be readily derived by writing

$$\begin{array}{c}\hfill \omega =\widehat{\omega }(1+C_3{ϵ}^3+C_4{ϵ}^4+\dots ),\end{array}$$(43)

with $ϵ=k_z{v}_{\mathrm{Ai}}/\widehat{\omega }$ seen as a small parameter. Expanding both sides of Eq. (26), one readily verifies that the ϵ-corrections start with third order in the parentheses of Eq. (43). An approximate solution that corrects $\widehat{\omega }$ to forth order reads

$$\begin{array}{c}\hfill \omega (k_y,k_z)\approx \widehat{\omega }[1-\frac{Q_3}{{(\widehat{\omega }d/v_{\mathrm{Ai}})}^3}+\frac{Q_4}{{(\widehat{\omega }d/v_{\mathrm{Ai}})}^4}],\end{array}$$(44)

where

$$\begin{array}{cc}\hfill Q_3& =\frac{{\pi }^2}{\sqrt{({\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}})({\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}}-1)}},\hfill \end{array}$$(45a)

$$\begin{array}{cc}\hfill Q_4& =\frac{3}{2}\frac{{\pi }^2}{({\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}})({\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}}-1)}·\hfill \end{array}$$(45b)

We now consider small values of k_yd. Some subtlety arises regarding what small means given the zone where trapped modes are forbidden. Three cases need to be discriminated in terms of k_z. We start with Case I where $0<k_z<k_z^{(1/2)}$. As indicated by Fig. 1a, the dispersion curve in this case starts from the intersection (P) between the curves m_e = 0 and n_i = π/2 therein. See k_z as given, and let the coordinates of P be denoted by [k_y^cut, ω^cut]. One finds with the definitions of m_e and n_i (Eqs. (23) and (24)) that

$$\begin{array}{cc}\hfill k_y^{\mathrm{cut}}& =\sqrt{{[k_z^{(1/2)}]}^2-k_z^2},\hfill \end{array}$$(46a)

$$\begin{array}{cc}\hfill {\omega }^{\mathrm{cut}}& =\frac{\pi /2}{\sqrt{1-{\rho }_{\mathrm{e}}/{\rho }_{\mathrm{i}}}}\frac{v_{\mathrm{Ai}}}{d}=k_z^{(1/2)}{v}_{\mathrm{Ae}}.\hfill \end{array}$$(46b)

Evidently, k_y^cut is some critical out-of-plane wavenumber, only beyond which trapped 3D modes are allowed for a given k_z. Defining a small parameter 0 < ϵ ≪ 1 as

$$\begin{array}{c}\hfill k_y=k_y^{\mathrm{cut}}(1+ϵ),\text{i.e.,}ϵ=k_y/k_y^{\mathrm{cut}}-1,\end{array}$$(47)

we proceed to examine how trapped modes behave in the immediate vicinity of the forbidden zone by writing

$$\begin{array}{c}\hfill \omega (k_y,k_z)={\omega }^{\mathrm{cut}}[1+Q_1ϵ+Q_2{ϵ}^2+\dots ].\end{array}$$(48)

Some tedious algebra yields that

$$\begin{array}{cc}\hfill Q_1& =1-\frac{k_z^2}{{(k_z^{(1/2)})}^2},\hfill \end{array}$$(49a)

$$\begin{array}{cc}\hfill Q_2& =\frac{1}{2{(k_z^{(1/2)}d)}^2}[1-\frac{k_z^2}{{(k_z^{(1/2)})}^2}]\hfill \\ \hfill & \{{(k_{z}d)}^2-{\left(\frac{\pi }{2}\right)}^4\frac{{[1-k_z^2/{(k_z^{(1/2)})}^2]}^3}{{[{\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}}-k_z^2/{(k_z^{(1/2)})}^2]}^2}\}.\hfill \end{array}$$(49b)

Evidently, both Q₁ and Q₂ approach zero when k_z → k_z^(1/2), meaning that the expansion (48) is of little use therein. Furthermore, Q₂ is a bit intricate in that it switches from negative to positive values when k_z increases through some critical value. Regardless, both Q₁ and Q₂ are involved in the approximate expression for m_e,

$$\begin{array}{c}\hfill m_{\mathrm{e}}\approx ϵ\sqrt{{(k_y^{\mathrm{cut}})}^2+{(k_z^{(1/2)})}^2(Q_1^2+2Q_2)},\end{array}$$(50)

thereby providing a direct measure of the trapping capability of the slab when k_y deviates slightly from k_y^cut.

Equation (48) proves largely adequate for explaining Figs. 5 and 6, provided that k_z is not that close to k_z^(1/2). We take the n_i distribution in Fig. 5 for instance. One finds from Eq. (48) that

$$\begin{array}{c}\hfill n_{\mathrm{i}}\approx \frac{\pi }{2}[1+\frac{{(k_y^{\mathrm{cut}})}^2}{{(k_z^{(1/2)})}^2}ϵ],\end{array}$$(51)

suggesting that n_i − π/2 scales largely as ∝ϵ when k_y deviates slightly from k_y^cut for a given k_z. One further finds that

$$\begin{array}{c}\hfill v_{\mathrm{ph}}\approx v_{\mathrm{Ae}}[1+\frac{Q_1^2+2Q_2-Q_1}{2}{ϵ}^2].\end{array}$$(52)

The phase speed v_ph therefore tends to decrease from the external Alfvén speed at a rather slow rate (∝ϵ²). Now consider Fig. 6. It follows from Eq. (48) that

$$\begin{array}{c}\hfill v_{\mathrm{gr},y}\approx \frac{{\omega }^{\mathrm{cut}}}{k_y^{\mathrm{cut}}}(Q_1+2Q_2ϵ)=v_{\mathrm{Ae}}\sqrt{1-\frac{k_z^2}{{(k_z^{(1/2)})}^2}}+2Q_2\frac{{\omega }^{\mathrm{cut}}}{k_y^{\mathrm{cut}}}ϵ.\end{array}$$(53)

It immediately follows that v_gr, y decreases from v_Ae toward zero when k_z increases from zero toward k_z^(1/2) along the outer edge of the forbidden zone. Furthermore, that Q₂ changes sign means that v_gr, y decreases (increases) from the edge value when k_y increases from k_y^cut for a given k_z that is below (above) some critical value. Equation (48) can also be used to derive an approximate expression for v_gr, z. This derivation is somehow complicated by the k_z-dependences of k_y^cut and ϵ (see Eqs. (46) and (47)). Some algebra yields that

$$\begin{array}{c}\hfill v_{\mathrm{gr},z}\approx v_{\mathrm{Ae}}\frac{k_z}{k_z^{(1/2)}}[1+(-1+2Q_2\frac{{(k_z^{(1/2)})}^2}{{(k_y^{\mathrm{cut}})}^2})ϵ].\end{array}$$(54)

Taking ϵ → 0, one deduces from Eq. (54) that v_gr, z varies from zero toward v_Ae when k_z is surveyed along the edge of the forbidden zone from zero toward k_z^(1/2). The ϵ-correction, one the other hand, immediately suggests that v_gr, z decreases from its edge value when k_y increases from k_y^cut for a k_z ensuring Q₂ < 0. This is in agreement with Fig. 6. Somehow puzzling is that the same k_y-dependence actually persists for those k_z that yield positive Q₂. This can be addressed by expanding the Q₂-term in Eq. (54) with the aid of Eq. (49), the result being

$$\begin{array}{c}\hfill 2Q_2\frac{{(k_z^{(1/2)})}^2}{{(k_y^{\mathrm{cut}})}^2}=\frac{1}{{(k_z^{(1/2)}d)}^2}\{{(k_{z}d)}^2-{\left(\frac{\pi }{2}\right)}^4\frac{{[1-k_z^2/{(k_z^{(1/2)})}^2]}^3}{{[{\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}}-k_z^2/{(k_z^{(1/2)})}^2]}^2}\}.\end{array}$$(55)

Evidently, this Q₂-term is consistently smaller than unity. It then follows that the coefficient of ϵ in Eq. (54) is always negative.

We now move on to Case II where $k_z^{(1/2)}<k_z<k_z^{(1)}$. This case is somehow peculiar in that there exist no trapped 2D modes (k_y = 0), whereas trapped 3D modes suddenly arise when k_y becomes finite, regardless of how small k_yd is. One readily verifies that m_e → 0 when k_y → 0, meaning that the leading-order solution to Eq. (26) reads ω ≈ k_zv_Ae. Defining ϵ = k_y²/k_z² ≪ 1, we look for improved approximate solutions in the form ω = k_zv_Ae(1 + Q₁ϵ − Q₂ϵ² + …), or equivalently

$$\begin{array}{c}\hfill \omega (k_y,k_z)\approx k_z{v}_{\mathrm{Ae}}[1+Q_1\frac{k_y^2}{k_z^2}-Q_2\frac{k_y^4}{k_z^4}+\dots ].\end{array}$$(56)

Some algebra yields that

$$\begin{array}{c}\hfill Q_1=\frac{1}{2},Q_2=\frac{1}{2}[\frac{1}{4}+\frac{{cot}^2(k_{z}d\sqrt{{\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}}-1})}{{\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}}-1}].\end{array}$$(57)

The reason for us to retain the Q₂ term is that it is involved in m_e², which approximates to

$$\begin{array}{c}\hfill m_{\mathrm{e}}^2\approx -k_z^2(Q_1^2-2Q_2){ϵ}^2.\end{array}$$

More specifically, one finds that

$$\begin{array}{c}\hfill m_{\mathrm{e}}\approx \frac{cot(k_{z}d\sqrt{{\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}}-1})}{\sqrt{{\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}}-1}}\frac{k_y^2}{k_z}·\end{array}$$(58)

For a fixed k_z, the slab therefore possesses some poorer trapping capability (m_e ∝ k_y²) than in Case I where $m_{\mathrm{e}}\propto (k_y-k_y^{(\mathrm{cut})})$ (Eq. (50)). More importantly, that m_e possesses a continuous k_y-dependence means that the sudden appearance of trapped 3D modes for k_y > 0 does not result in any abrupt change in the the temporal response of the slab to impulsive, localized, 3D drivers.

Equation (56) proves adequate for explaining some key features of Figs. 5 and 6. One deduces from Eq. (56) that

$$\begin{array}{c}\hfill n_{\mathrm{i}}\approx \sqrt{{\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}}-1}\sqrt{k_y^2+k_z^2}.\end{array}$$(59)

When k_y → 0, the quantity n_i therefore increases monotonically from π/2 to π as k_z increases from k_z^(1/2) to k_z⁽¹⁾. The contours of n_i for small k_y, on the other hand, are expected to be concentric circles as indeed seen in Fig. 5. One further deduces that

$$\begin{array}{c}\hfill v_{\mathrm{ph}}\approx v_{\mathrm{Ae}}[1-\frac{1}{2}\frac{{cot}^2(k_{z}d\sqrt{{\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}}-1})}{({\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}}-1)}\frac{k_y^4}{k_z^4}],\end{array}$$(60)

meaning that v_ph decreases from v_Ae at a very slow rate (∝k_y⁴) when k_y increases for a given k_z. We now move on to Fig. 6, to explain which it suffices to retain only the Q₁ term in Eq. (56). One finds that

$$\begin{array}{c}\hfill v_{\mathrm{gr},y}\approx v_{\mathrm{Ae}}\frac{k_y}{k_z},\end{array}$$(61)

thereby explaining why the v_gr, y contours are largely radially directed. One further finds that

$$\begin{array}{c}\hfill v_{\mathrm{gr},z}\approx v_{\mathrm{Ae}}(1-\frac{k_y^2}{2k_z^2}),\end{array}$$(62)

meaning that v_gr, z for a given k_z tends to decrease rather slowly with k_y. Equation (62) also explains why the v_gr, z contours are largely radially aligned when k_z is not that close to k_z⁽¹⁾.

We now address Case III where k_z > k_z⁽¹⁾. It is textbook material that trapped 2D modes (k_y = 0) are now allowed, their phase speeds ω^2D(k_z)/k_z decreasing monotonically from v_Ae to v_Ai when k_z increases from k_z⁽¹⁾ to large values (e.g., Roberts 2019, Chapter 5). One further finds that

$$\begin{array}{c}\hfill {\omega }^{2\mathrm{D}}(k_z)\approx k_z{v}_{\mathrm{Ae}}[1-\frac{({\rho }_{\mathrm{i}}/{\rho }_{\mathrm{e}}-1){\pi }^2}{2}{(\frac{k_z}{k_z^{(1)}}-1)}^2]\end{array}$$(63)

when k_z exceeds k_z⁽¹⁾ only marginally. A similar expression is available for trapped 2D sausage modes (Li et al. 2018, Eq. (22)), and Eq. (63) can be derived with the same procedure therein. Likewise, one finds that

$$\begin{array}{c}\hfill {\omega }^{2\mathrm{D}}(k_z)\approx k_z{v}_{\mathrm{Ai}}\sqrt{1+{(3\pi /2)}^2{(k_{z}d)}^{-2}},k_{z}d\gg 1,\end{array}$$(64)

by specializing Eq. (16) in Li et al. (2013) to the pressureless limit. We note that the quantity $n_{\mathrm{i}}=\sqrt{{({\omega }^{2\mathrm{D}}/v_{\mathrm{Ai}})}^2-k_z^2}$ for trapped 2D overtones behaves similarly to transverse fundamentals in that n_i also increases monotonically with k_z (see Fig. 5). Equations (63) and (64) help explicitly show that n_id → π when k_z → k_z⁽¹⁾ and n_id → 3π/2 when k_zd → ∞.

We now consider what happens when 0 < k_y²/k_z² ≪ 1. One readily verifies that Eq. (34) holds for transverse overtones as well, as long as the quantities ω^2D and $m_{\mathrm{e}}^{2\mathrm{D}}$ are interpreted appropriately. We rewrite Eq. (34) as

$$\begin{array}{c}\hfill \omega (k_y,k_z)\approx {\omega }^{2\mathrm{D}}(k_z)[1+\mathcal{C}(k_z){(k_{y}d)}^2]\end{array}$$(65)

for the ease of description. It follows from Eqs. (63) and (64) that $m_{\mathrm{e}}^{2\mathrm{D}}=\sqrt{k_z^2-{({\omega }^{2\mathrm{D}})}^2/v_{\mathrm{Ae}}^2}\propto (k_z/k_z^{(1)}-1)$ when k_z barely exceeds k_z⁽¹⁾, while $m_{\mathrm{e}}^{2\mathrm{D}}\propto k_z$ when k_zd ≫ 1. The quantity $m_{\mathrm{e}}^{2\mathrm{D}}$ itself turns out to be a monotonically increasing function of k_z. Equation (34) then indicates that there must exist some critical axial wavenumber, k_z^2D, c, such that 𝒞 < 0 (𝒞 > 0) when k_z⁽¹⁾ < k_z < k_z^2D, c (k_z > k_z^2D, c). One therefore deduces from Eq. (65) that v_gr, y = ∂ω/∂k_y varies from zero when k_y → 0 toward negative (positive) values for those values of k_z that ensure 𝒞 < 0 (𝒞 > 0). This is what one sees from Fig. 6, and the v_gr, y = 0 contour therein actually identifies k_z^2D, cd ≈ 2.74 for the chosen ρ_i/ρ_e = 3. We restrict ourselves to k_z⁽¹⁾ < k_z < k_z^2D, c. Figure 6 shows the subtle feature that v_gr, z tends to decrease (increase) from its edge value (k_z → 0) when k_y increases for relatively small (large) values of k_z. This can be adequately addressed by taking the k_z-derivative of Eq. (65), even though the cumbersome d𝒞/dk_z term somehow complicates the matter.

Figure 7 applies the same format of Fig. 4 and the associated classification scheme to the first transverse overtone. We note that the filled contours display only a limited range for the force ratio Λ_i, despite that Λ_i can be evaluated wherever trapped modes are permitted. Consider trapped 2D modes first (k_y = 0), which are recalled to arise for k_z > k_z⁽¹⁾. One readily finds from Eqs. (63) and (64) that Λ_i → ρ_i/ρ_e − 1 when k_z approaches k_z⁽¹⁾ from above, and Λ_i → 0 when k_zd → ∞. For highly oblique propagation (k_y²/k_z² ≫ 1), however, Equation (44) suggests that Λ_i = ω²/k_z²v_Ai² − 1 approximates to [π² + (k_yd)²)]/(k_zd)² and may therefore take extremely large values.

Fig. 7. Similar to Fig. 4, but for the first transverse overtone. Only a limited range is plotted for the quantity Λi because Λi takes up too diverse a range. Panel c overviews how the overtone behaves on the ky − kz plane. Trapped 3D modes are forbidden where $\sqrt{k_y^2+k_z^2}\le k_z^{(1/2)}$. When allowed, the overtone is further classified into three regimes based on the angle information (see text for more details).

Figure 7c overviews how the first overtone behaves on the k_y − k_z plane. The disc $\sqrt{k_y^2+k_z^2}\le k_z^{(1/2)}$, labeled 0, is recalled to be where trapped 3D modes are forbidden. Our classification scheme is then applied, and all the three subgroups are now relevant. The B₀-straddling subgroup (with v_gr, y < 0, labeled 1) remains associated with force ratios Λ_i ≳ 1, despite that this subgroup is restricted to a substantially smaller portion when compared with the fundamental. The B₀-same-side F subgroup corresponds to the vast portion to the right of the |∠(v_gr, e_z)| = 15° contour. Our additional criterion |∠(v_gr, v_ph)| < 15° does not provide further constraints per se. In fact, Figure 7b indicates a small ∠(v_gr, v_ph) almost everywhere, which is particularly true for oblique propagation (k_y²/k_z² ≫ 1). This feature, in conjunction with the overall tendency for Λ_i to be (relatively) large, means that the B₀-same-side F subgroup somehow closely resembles what happens in unbounded uniform media. However, the B₀-same-side A subgroup deviates considerably from its counterpart in an unbounded uniform medium, in which case one recalls a Λ ≲ 0.07. The force ratio Λ_i for the B₀-same-side A subgroup here, in contrast, may readily exceed ρ_i/ρ_e − 1. We note that |∠(v_gr, e_z)| evaluates to 45° when Λ = 1 in an unbounded uniform medium (see Eq. (15)). Our specific threshold of 15° for mode categorization is therefore only for illustrative purposes. A larger threshold (up to 45°) does not impact the categorization of the transverse fundamental, given that |∠(v_gr, e_z)| does not exceed ∼11° (Fig. 4). Nonetheless, a larger threshold will broaden the k_y − k_z portion where the first overtone qualifies as B₀-same-side A in Fig. 7c. The rest of the B₀-same-side portion always qualifies as B₀-same-side F modes; |∠(v_gr, v_ph)| is consistently smaller than ∼15° and therefore not a discriminating factor. Importantly, the threshold angle is not involved hereafter; only the distinction between B₀-straddling and B₀-same-side is employed.

4. Implications for impulsively excited 3D kink motions

The question arises of how the above-presented mode analysis connects to the spatio-temporal evolution of a density-enhanced slab in response to localized initial exciters. We start by noting that the problem at hand is an initial value problem (IVP) governed by the time-dependent MHD Equations (1) and (2) over the infinite [x, y, z]-volume. Conceptually, the nominal boundaries x, y, z → ±∞ need to be such that no boundary conditions (BCs) are specified or equivalently that the dependent variables do not diverge (see, e.g., Gao et al. 2024, and references therein). The most unambiguous way to specify initial exciters, on the other hand, is to perturb the velocities only (v₁(x, y, z, t = 0)≠0, B₁(x, y, z, t = 0) = 0; see Sect. 6.1.2 in Goedbloed et al. 2019). We proceed with the following assumptions for definitiveness. First, the only non-vanishing independent quantity at t = 0 is v_1x. Second, v_1x(x, y, z, t = 0) is even in x to ensure kink parity, and is even in both y and z as well. Third, the implementation of v_1x(x, y, z, t = 0) is such that one needs only to consider the range of [k_y, k_z] examined in Figs. 3 and 6.

Quantitative progress can be made by presenting the formal solution to the IVP (Li et al. 2023),

$$\begin{array}{cc}& v_{1x}(x,y,z,t)={\displaystyle {\int }_{-\infty }^{\infty }}\mathrm{d}{k}_y{\displaystyle {\int }_{-\infty }^{\infty }}\mathrm{d}{k}_z\hfill \\ \hfill & \{{\displaystyle \sum _j}\left[{\mathcal{F}}_j(x;k_y,k_z){\mathrm{e}}^{ı({\omega }_{j}t-k_{y}y-k_{z}z)}\right]+\text{improper}\}.\hfill \end{array}$$(66)

Equation (66) is so written by following the approach of Fourier integrals (e.g., Oliver et al. 2014, 2015) or equivalently the idea of decomposing the time-dependent wave fields into the eigenmodes of the MHD force operator (e.g., Li et al. 2022; Wang et al. 2023). We adopt the nomenclature of the latter. The decomposition in the y- or z-direction is self-evident, and one then ends up with an eigenvalue problem with x being the only independent variable and [k_y, k_z] serving as parameters. At any given pair [k_y, k_z], the eigenspectrum of the MHD force operator for our slab equilibrium comprises two subspectra. The point subspectrum is populated by a finite number of eigensolutions, which are labeled by j and are equivalent to our trapped modes (m_e² > 0). By construction, ω_j = ω_j(k_y, k_z) is dictated by the DR. The other subspectrum (with m_e² < 0), called improper continuum, comprises those eigensolutions whose eigenfrequencies continuously cover the range from some critical frequency out to infinity. This improper continuum arises in our context as a result of the physical requirement that no BCs should be specified at x → ±∞. Regardless, the contribution from the continuum eigenmodes (labeled “improper” in Eq. (66)) attenuates rapidly with time and hence are not of further interest here.

What survives at large times is the contribution from trapped modes. We focus on the slab axis x = 0, supposing further that only the transverse fundamental and its first overtone are relevant. These simplifications suffice to illustrate both the usefulness of and intricacies in the MSP that we adopt to digest the morphological features of wave propagation. Suppose that the MSP applies to some [y, z, t] or more appropriately to some [y/t, z/t] with [y, z, t] seen as given. Equation (66), while nominally involving all [k_y, k_z], is actually dominated only by a series of narrow regions (or wavepackets in physical terms) on the k_y − k_z plane. Let wavepackets (WPs) be numbered by n and represented by the central wavenumbers [K_n, y, K_n, z]. It follows from the general theory of the MSP that (e.g., Whitham 1974, Chapter 11)