Impulsively Generated Wave Trains in Coronal Structures. I. Effects of Transverse Structuring on Sausage Waves in Pressureless Tubes

Appropriate for the solar corona, we work in the framework of cold MHD, for which the gas pressure is neglected and hence slow modes become irrelevant. We consider sausage waves in a structured corona modeled by a density-enhanced straight tube with mean radius R embedded in a uniform magnetic field ${\boldsymbol{B}}$. Both ${\boldsymbol{B}}$ and the tube are directed in the z-direction in a cylindrical coordinate system $(r,\theta ,z)$. The left column of Figure 1 offers an illustration of this setup. The equilibrium density is assumed to be a function of r only and of the form

$$\begin{eqnarray}&&\rho (r)={\rho }_{{\rm{e}}}+({\rho }_{{\rm{i}}}-{\rho }_{{\rm{e}}})f(r),\end{eqnarray} \tag{ 1 }$$

where the function $f(r)$ decreases from unity at r = 0 to zero when r approaches infinity. In other words, the equilibrium density decreases from ${\rho }_{{\rm{i}}}$ at tube axis to ${\rho }_{{\rm{e}}}$ far from the tube. The Alfvén speed is given by ${v}_{{\rm{A}}}(r)=B/\sqrt{4\pi \rho (r)}$, which increases from ${v}_{\mathrm{Ai}}$ to ${v}_{\mathrm{Ae}}$ with ${v}_{\mathrm{Ai},{\rm{e}}}=B/\sqrt{4\pi {\rho }_{{\rm{i}},{\rm{e}}}}$.

Zoom In Zoom Out Reset image size

While our methodology applies for arbitrary choices of $f(r)$, we nonetheless focus only on the following three families of profiles. The first one is called "μ power" and given by

$$\begin{eqnarray}&&f(r)=\displaystyle \frac{1}{1+{(r/R)}^{\mu }}.\end{eqnarray} \tag{ 2 }$$

The next one is called "outer μ" and given by

$$\begin{eqnarray}f(r)=\left\{\begin{array}{ll}1, & 0\leqslant r\leqslant R,\\ {(r/R)}^{-\mu }, & r\geqslant R.\end{array}\right.\end{eqnarray} \tag{ 3 }$$

The third one is called "inner μ" and in the form

$$\begin{eqnarray}f(r)=\left\{\begin{array}{ll}1-{\left(\displaystyle \frac{r}{R}\right)}^{\mu }, & 0\leqslant r\leqslant R,\\ 0, & r\geqslant R.\end{array}\right.\end{eqnarray} \tag{ 4 }$$

These profiles are illustrated in the right column of Figure 1, where ${\rho }_{{\rm{i}}}$/${\rho }_{{\rm{e}}}$ is arbitrarily chosen to be 5. Evidently, μ is a measure of the profile steepness, and all profiles converge to a top-hat one when $\mu \to \infty$. We will focus on coronal tubes with parameters typical of EUV active region loops, namely with density contrasts between 2 and 10 (e.g., Aschwanden et al. 2004). In addition, we will consider only $\mu \geqslant 1$ because otherwise the "μ power" and "inner μ" profiles become cusped at r = 0.

Now let $\delta \rho ,\delta {\boldsymbol{v}}$ and $\delta {\boldsymbol{b}}$ represent the perturbations to the density, velocity, and magnetic field, respectively. Working in the ideal cold MHD framework, and specializing to the axisymmetric case ($\partial /\partial \theta \equiv 0$), one finds that $\delta {b}_{\theta },\delta {v}_{\theta }$, and $\delta {v}_{z}$ vanish. Eliminating $\delta {b}_{r}$ and $\delta {b}_{z}$, one derives a single equation for $\delta {v}_{r}(r,z;t)$ (e.g., Nakariakov et al. 2012; Chen et al. 2015a, 2015b)

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}\delta {v}_{r}}{\partial {t}^{2}}={v}_{{\rm{A}}}^{2}(r)\left(\displaystyle \frac{{\partial }^{2}}{\partial {z}^{2}}+\displaystyle \frac{{\partial }^{2}}{\partial {r}^{2}}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}-\displaystyle \frac{1}{{r}^{2}}\right)\delta {v}_{r}.\end{eqnarray} \tag{ 5 }$$

In addition, the density perturbation $\delta \rho$ is governed by

$$\begin{eqnarray}&&\displaystyle \frac{\partial \delta \rho }{\partial t}=-\left[\delta {v}_{r}\displaystyle \frac{\partial \rho }{\partial r}+\displaystyle \frac{\rho }{r}\displaystyle \frac{\partial (r\delta {v}_{r})}{\partial r}\right].\end{eqnarray} \tag{ 6 }$$

For each profile with a given combination of $[{\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}},\mu ]$, we will always start with an eigenmode analysis to establish the dispersive behavior of the trapped modes, the behavior of the group speed curves in particular. Fourier-decomposing any perturbation $\delta g(r,z;t)$ as

$$\begin{eqnarray}&&\delta g(r,z;t)=\mathrm{Re}\{\tilde{g}(r)\exp [-i(\omega t-{kz})]\},\end{eqnarray} \tag{ 7 }$$

one finds from Equation (5) that

$$\begin{eqnarray}&&\displaystyle \frac{1}{r}\displaystyle \frac{d}{{dr}}\left(r\displaystyle \frac{d\tilde{\xi }}{{dr}}\right)+\left(\displaystyle \frac{{\omega }^{2}}{{v}_{{\rm{A}}}^{2}}-{k}^{2}-\displaystyle \frac{1}{{r}^{2}}\right)\tilde{\xi }=0,\end{eqnarray} \tag{ 8 }$$

where $\tilde{\xi }=i{\tilde{v}}_{r}/\omega$ is the Fourier amplitude of the radial Lagrangian displacement. Equation (8) can be formulated into a standard eigenvalue problem (EVP) when supplemented with the following boundary conditions (BCs)

$$\begin{eqnarray}&&\tilde{\xi }(r=0)=0,\tilde{\xi }(r\to \infty )\to 0.\end{eqnarray} \tag{ 9 }$$

To solve this EVP, we employ a MATLAB boundary-value-problem solver BVPSuite in its eigenvalue mode (see Kitzhofer et al. 2009, for a detailed description of the code). In the appendix of Li et al. (2014), we performed an extensive validation study of this code using available analytical solutions in the context of coronal seismology. The end result is that the dimensionless angular frequency $\omega R/{v}_{\mathrm{Ai}}$ can be formally expressed as

$$\begin{eqnarray}&&\displaystyle \frac{\omega R}{{v}_{\mathrm{Ai}}}={ \mathcal G }[{kR};{\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}},f(r)].\end{eqnarray} \tag{ 10 }$$

Note that both ω and the axial wavenumber k are real-valued. The axial phase and group speeds simply follow from the definitions ${v}_{\mathrm{ph}}=\omega /k$ and ${v}_{\mathrm{gr}}=d\omega /{dk}$, respectively.

With the group speed curves at hand, we will then examine impulsively generated sausage waves by asking how a coronal tube responds to a localized initial perturbation. For this purpose, we choose to solve Equation (5), which is equivalent to, but simpler than, the linearized ideal cold MHD equations. However, to initiate our simulations, an initial condition (IC) is needed for $\partial \delta {v}_{r}/\partial t$ in addition to the IC for $\delta {v}_{r}$. This can be found by noting that if we initially perturb only $\delta {v}_{r}$, then

$$\begin{eqnarray}&&\displaystyle \frac{\partial \delta {v}_{r}}{\partial t}(r,z;t=0)=0,\end{eqnarray} \tag{ 11 }$$

given the absence of $\delta {b}_{r}$ and $\delta {b}_{z}$ at t = 0. On the other hand, the following form is adopted as the initial radial velocity perturbation,

$$\begin{eqnarray}\delta {v}_{r}(r,z;t=0) & = & A{{\rm{e}}}^{1/2}{v}_{\mathrm{Ai}}\left(\displaystyle \frac{r}{{\sigma }_{r}}\right)\\ & & \times \exp \left(-\displaystyle \frac{{r}^{2}}{2{\sigma }_{r}^{2}}\right)\exp \left(-\displaystyle \frac{{z}^{2}}{2{\sigma }_{z}^{2}}\right),\end{eqnarray} \tag{ 12 }$$

which ensures the parity of the generated wave trains by not displacing the tube axis. A constant $\exp (1/2)$ is introduced such that initially $\delta {v}_{r}$ in units of ${v}_{\mathrm{Ai}}$ attains a maximum of A, which is arbitrarily chosen to be unity. Furthermore, ${\sigma }_{r}$ and ${\sigma }_{z}$ determine the extent to which the initial perturbation spans in the radial and axial directions, respectively. Throughout this study, we choose ${\sigma }_{r}={\sigma }_{z}=R/\sqrt{2}$ to ensure that the initial perturbation is neither too localized nor too extended. This perturbation is shown by the red arrows in the left column of Figure 1.

Having specified the ICs, we then develop a simple finite-difference (FD) scheme to numerically solve Equation (5) in a computational domain extending from 0 to L_r ($-{L}_{z}/2$ to ${L}_{z}/2$) in the r-(z-) direction. A uniform grid spacing ${\rm{\Delta }}z=0.08\,R$ is adopted in the z-direction for simplicity. However, to speed up our computations, the grid points in the r-direction are chosen to be nonuniformly distributed. We fix the spacing at ${\rm{\Delta }}{r}_{1}=0.02\,R$ when $r\leqslant 3\,R$. We then allow it to increase as ${\rm{\Delta }}{r}_{j+1}=1.025{\rm{\Delta }}{r}_{j}$ with j being the grid index until the spacing reaches $\sqrt{{\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}}}{\rm{\Delta }}{r}_{1}$, beyond which it remains a constant again. Despite this complication, we make sure that the FD approximations to the spatial derivatives in Equation (5) are second order accurate. Time integration is implemented with the leap-frog method, for which a uniform time step ${\rm{\Delta }}t=0.6{{\rm{\Delta }}}_{\min }/{v}_{{\rm{A}},\max }$ is adopted to comply with the Courant condition to ensure numerical stability. Here ${{\rm{\Delta }}}_{\min }$ (${v}_{{\rm{A}},\max }$) represents the smallest (largest) value that the grid spacing (Alfvén speed) attains in the computational domain. At the left boundary r = 0, we fix $\delta {v}_{r}$ at zero in view of the symmetric properties of sausage modes. We choose both L_r and L_z to be sufficiently large such that the perturbations will not have sufficient time to reach the rest of the boundaries (at $r={L}_{r}$ and $z=\pm {L}_{z}/2$). In this regard, the BCs there are irrelevant. Nonetheless, we fix $\delta {v}_{r}$ at 0 for simplicity in practice. Furthermore, L_r and L_z are chosen to be 250 R and 500 R, respectively.

Given that available measurements attributed to impulsively generated sausage waves are primarily made with imaging rather than spectroscopic instruments, density perturbations will be more relevant compared with velocity perturbations. In fact, $\delta \rho$ can be readily evaluated with Equation (6), which is also solved with a scheme second order accurate in both time and space. To initiate this evaluation, we note that $\delta \rho (r,z;t=0)=0$ because no density perturbation is applied. The computed density perturbations are then sampled at a point along the tube axis and some distance h away from the impulsive source. A value of $h=75\,R$ is chosen such that only trapped modes play a role in determining $\delta \rho (r=0,z=h,t)$.

To summarize at this point, for each family of density profiles, we will adopt the following procedure to examine the signatures of impulsively generated wave trains.

• 1.  
  The ${v}_{\mathrm{gr}}$–k curves will be found by solving the EVP (Equations (8) and (9)) with BVPSuite. We will also present analytical expressions for ${v}_{\mathrm{gr}}$ in some limiting cases to help better understand the behavior of the group speed curves. In particular, we will experiment with different combinations of $[{\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}},\mu ]$ to see whether cutoff wavenumbers exist and whether the ${v}_{\mathrm{gr}}$–k curves behave in a non-monotonical fashion.
• 2.  
  We will solve Equation (5) with the initial conditions (11) and (12), for which we fix $[A,{\sigma }_{r},{\sigma }_{z}]$ at $[1,R/\sqrt{2},R/\sqrt{2}]$. Equation (6) is simultaneously solved at each time step when Equation (5) is advanced. A time series corresponding to the density perturbation sampled at a specified location, $\delta \rho (r=0,z=75\,R,t)$, will be subsequently analyzed with wavelet analysis, which is placed in the context of the ${v}_{\mathrm{gr}}$–k curves.

Before proceeding, we note that whether cutoff wavenumbers exist is solely determined simply by the asymptotic behavior of $f(r)$ when r is large. With the aid of Kneser's oscillation theorem, LN15 established a general result that cutoff wavenumbers are present (absent) if $f(r)$ is steeper (less steep) than ${r}^{-2}$ at large distances. However, solving the EVP is necessary to find out whether ${v}_{\mathrm{gr}}$ possesses a non-monotonical dependence on k.

For top-hat profiles, the solution to Equation (8) in the uniform interior (exterior) is proportional to ${J}_{1}({nr})$ (${K}_{1}({mr})$), where

$$\begin{eqnarray}&&{n}^{2}=\displaystyle \frac{{\omega }^{2}}{{v}_{\mathrm{Ai}}^{2}}-{k}^{2},\hspace{0.5cm}{m}^{2}={k}^{2}-\displaystyle \frac{{\omega }^{2}}{{v}_{\mathrm{Ae}}^{2}},\end{eqnarray} \tag{ 13 }$$

and J_j and K_j represent Bessel functions of the first kind and modified Bessel function of the second kind, respectively (here j = 1). Both n² and m² are non-negative for trapped modes. The dispersion relation (DR) reads (e.g., Edwin & Roberts 1983; Kopylova et al. 2007)

$$\begin{eqnarray}&&n\displaystyle \frac{{J}_{0}({nR})}{{J}_{1}({nR})}=-m\displaystyle \frac{{K}_{0}({mR})}{{K}_{1}({mR})}.\end{eqnarray} \tag{ 14 }$$

There are an infinite number of branches of solutions to Equation (14). Ordered by the radial harmonic number l, all of these branches possess a low wavenumber cutoff given by

$$\begin{eqnarray}&&{k}_{{\rm{c}},l}=\displaystyle \frac{{j}_{0,l}}{R\sqrt{{\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}}-1}},\end{eqnarray} \tag{ 15 }$$

where ${j}_{0,l}$ is the lth zero of J₀ ($l=1,2,\ldots$). Evidently, ${k}_{{\rm{c}},l}$ increases with l. Furthermore, it is no surprise to see that a cutoff wavenumber exists for all radial harmonic numbers, given that $f(r)$ is identically zero when $r\gt R$ and consequently steeper than ${r}^{-2}$.

It will also be informative to find out how ${v}_{\mathrm{gr}}$ behaves at large kR. It turns out that (e.g., Li et al. 2014; Oliver et al. 2015; Yu et al. 2016a)

$$\begin{eqnarray}&&\displaystyle \frac{{v}_{\mathrm{gr}}^{2}}{{v}_{\mathrm{Ai}}^{2}}\approx \displaystyle \frac{1}{1+{j}_{1,l}^{2}/{({kR})}^{2}},\end{eqnarray} \tag{ 16 }$$

where ${j}_{1,l}$ is the lth zero of J₁ ($l=1,2,\cdots$). This means that ${v}_{\mathrm{gr}}$ should eventually approach ${v}_{\mathrm{Ai}}$ from below when kR increases. Note that Equation (16) comes directly from the approximate behavior of ${v}_{\mathrm{ph}}$ for sufficiently large kR

$$\begin{eqnarray}&&\displaystyle \frac{{v}_{\mathrm{ph}}^{2}}{{v}_{\mathrm{Ai}}^{2}}\approx 1+\displaystyle \frac{{j}_{1,l}^{2}}{{k}^{2}{R}^{2}}.\end{eqnarray} \tag{ 17 }$$

Figure 2 presents the dependence on the axial wavenumber k of the axial phase (the upper row) and group (lower) speeds for both a density contrast of 3 (the left column) and of 10 (right). The solid and dashed curves correspond to radial harmonic numbers of 1 and 2, respectively. These numerical results are found with BVPSuite and agree exactly with the numerical solutions to the analytical DR (14). One sees that cutoff wavenumbers exist for both density contrasts and for both branches, and ${k}_{{\rm{c}},l}$ increases with l but decreases with ${\rho }_{{\rm{i}}}$/${\rho }_{{\rm{e}}}$ as expected from Equation (15). Furthermore, all of the ${v}_{\mathrm{gr}}$–k curves show a non-monotonical dependence on k. Letting + (−) denote the tendency for ${v}_{\mathrm{gr}}$ to increase (decrease) with k, this "−/+" behavior can be partly explained by the fact that, asymptotically, ${v}_{\mathrm{gr}}$ approaches ${v}_{\mathrm{Ai}}$ from below for arbitrary l and ${\rho }_{{\rm{i}}}$/${\rho }_{{\rm{e}}}$.

Zoom In Zoom Out Reset image size

Figure 3 displays the temporal evolution (the upper row) and the corresponding Morlet spectra (lower) of the density perturbations $\delta \rho$ sampled at a distance $h=75\,R$ along the tube axis for both a density contrast ${\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}}=3$ (the left column) and 10 (right). In the lower row, the left vertical axis represents the angular frequency ω, and the right one represents the period P. The Morlet spectra are created by using the standard wavelet toolkit devised by Torrence & Compo (1998), and the dashed contours represent the $95 \%$ confidence level computed by assuming a white-noise process for the mean background spectrum. The dotted vertical lines correspond to the arrival times of wavepackets traveling at the internal (${v}_{\mathrm{Ai}}$) and external (${v}_{\mathrm{Ae}}$) Alfvén speeds. Furthermore, the yellow curves represent ω as a function of h/${v}_{\mathrm{gr}}$ with the radial harmonic number increasing from bottom to top. These curves are essentially a replot of the curves in the lower row of Figure 2.

Zoom In Zoom Out Reset image size

Consider the lower row first. Placing the group speed curves on the Morlet spectrum indicates that the energy contained in the initial perturbation predominantly goes to the fundamental radial mode. Furthermore, the fact that the Morlet spectrum looks like a "crazy tadpole" is in close agreement with the reasoning by Edwin & Roberts (1986). The narrow tail of the tadpole is threaded by the portion of the $\omega \mbox{--}h/{v}_{\mathrm{gr}}$ curve where ω varies little. On the other hand, the broad head occurs in the portion where ${v}_{\mathrm{gr}}$ is close to its local minimum such that multiple wavepackets arrive simultaneously to yield some stronger Morlet power. Comparing Figures 3(b) and (d), one sees that this behavior is clearer for stronger density contrasts. Now examining the upper row, one finds that wavepackets corresponding to higher radial harmonics also show up (see, e.g., the interval $t\gtrsim 110R/{v}_{\mathrm{Ai}}$ in Figure 3(c)). They are simply too weak to discern in the Morlet spectra.

LN15 showed that $\mu =2$ is a dividing line separating the cases where cutoff wavenumbers do not exist ($\mu \lt 2$) from those where they do exist ($\mu \gt 2$). Care needs to be exercised when $\mu =2$, in which case LN15 showed that the cutoff wavenumber for radial harmonic number l is given by

$$\begin{eqnarray}&&{k}_{{\rm{c}},l}=\displaystyle \frac{1}{R\sqrt{{\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}}-1}}.\end{eqnarray} \tag{ 18 }$$

Note that ${k}_{{\rm{c}},l}$ does not depend on l. Our numerical computations with BVPSuite demonstrate that the expected behavior of cutoff wavenumbers indeed holds.

Some expectations on whether the ${v}_{\mathrm{gr}}$–k curves behave in a non-monotonical manner can be made by examining the asymptotic behavior of the group speeds at large axial wavenumbers. LN15 showed that while, in general, the DR cannot be analytically established, it is possible to find an approximate one with the WKB approach when ${\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}}\gg 1$. This approximate DR at large kR yields that

$$\begin{eqnarray}&&\displaystyle \frac{{v}_{\mathrm{ph}}^{2}}{{v}_{\mathrm{Ai}}^{2}}\approx 1+{\left(\displaystyle \frac{c}{{kR}}\right)}^{\beta },\end{eqnarray} \tag{ 19 }$$

where c is a constant that depends on l and μ, and $\beta =2\mu /(\mu +2)$.³ It then follows that

$$\begin{eqnarray}&&\displaystyle \frac{{v}_{\mathrm{gr}}^{2}}{{v}_{\mathrm{Ai}}^{2}}\approx 1+(1-\beta ){\left(\displaystyle \frac{c}{{kR}}\right)}^{\beta }.\end{eqnarray} \tag{ 20 }$$

Given that $\beta \lt 1$ ($\beta \gt 1$) when $\mu \lt 2$ ($\mu \gt 2$), one finds that ${v}_{\mathrm{gr}}$ approaches ${v}_{\mathrm{Ai}}$ from above (below) asymptotically. Using the notations of "−" and "+," this means that asymptotically the ${v}_{\mathrm{gr}}$–k curves are of the "−" ("+") type when $\mu \lt 2$ ($\mu \gt 2$). The case with $\mu =2$ is special because, while Equation (20) suggests that ${v}_{\mathrm{gr}}^{2}/{v}_{\mathrm{Ai}}^{2}-1\to 0$, it does not say whether ${v}_{\mathrm{gr}}$ approaches ${v}_{\mathrm{Ai}}$ from above or below. One way to tell this is to extend Equations (19) and (20) by incorporating terms of even higher order in $1/({kR})$. Alternatively, we may employ full numerical solutions to the EVP by using BVPSuite. These numerical solutions indicate that Equations (19) and (20) hold for arbitrary combinations of $[{\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}},\mu ]$. In particular, they are true as long as ${\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}}\gt 1$, way beyond the nominal range of applicability of the WKB analysis that assumes ${\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}}\gg 1$. It is just that, in addition to l and μ, the constant c depends on the density contrast as well. When μ increases, this dependence gets weaker in that c approaches ${j}_{1,l}$ as happens for top-hat profiles (see Equation (17)).

Figure 4 largely follows Figure 2 in presenting the dependence on the axial wavenumber k of the axial phase and group speeds. The difference is that now the parameter μ is relevant, and a number of different values for μ are examined as given by the curves in different colors. Note that all of the curves are found with BVPSuite. Consider cutoff wavenumbers first. From both columns one sees that when $\mu \lt 2$ (the red and blue curves), both the first and second branches start at k = 0. In fact, our numerical solutions indicate that this holds for sausage modes of arbitrary radial harmonic number l. On the contrary, when $\mu \gt 2$ (the yellow and black curves), one sees that cutoff wavenumbers exist for sausage modes of any radial harmonic number l, and ${k}_{{\rm{c}},l}$ increases with l. In the particular case where $\mu =2$ (the green curves), ${k}_{{\rm{c}},l}$ is non-zero but does not depend on l, as expected from Equation (18). Now examine the wavenumber dependence of the group speeds. One sees that ${v}_{\mathrm{gr}}$ behaves in a monotonical "−" (non-monotonical "−/+") manner when $\mu \lt 2$ ($\mu \gt 2$), which can be partly explained by the asymptotic behavior of ${v}_{\mathrm{gr}}$ as given by Equation (20). When $\mu =2$, the green curves indicate that ${v}_{\mathrm{gr}}$ approaches ${v}_{\mathrm{Ai}}$ from above at large wavenumbers. In other words ${v}_{\mathrm{gr}}^{2}/{v}_{\mathrm{Ai}}^{2}-1\to {0}^{+}$, which cannot be directly deduced from Equation (20).

Zoom In Zoom Out Reset image size

To sum up here, the behavior of the group speed curves is determined solely by μ. Cutoff wavenumbers are absent (present) when $\mu \lt 2$ ($\mu \geqslant 2$), and the ${v}_{\mathrm{gr}}$–k curves behave in a "−" ("−/+") manner when $\mu \leqslant 2$ ($\mu \gt 2$).

Figure 5 presents the temporal evolution and Morlet spectra for density perturbations sampled at a distance of $h=75\,R$ along the tube axis pertinent to a μ of 1. Note that the third branch corresponding to a radial harmonic number l of 3 is also plotted in the lower row. For the present choice of μ, trapped modes exist for arbitrary wavenumbers and their group speeds decrease monotonically with k from ${v}_{\mathrm{Ae}}$ to ${v}_{\mathrm{Ai}}$. Consequently, two qualitative differences from the top-hat case arise. First, the signals in the upper row do not extend beyond h/${v}_{\mathrm{Ai}}$. Second, now the radial harmonics can be more readily excited because they can receive a substantial fraction of the energy contained in the initial perturbation. This is particularly clear in Figure 5(d), where not only modes of l = 2 but those of l = 3 can be clearly discerned. Instead of tadpoles, the Morlet spectra look more like carps, for which the fins derive from modes of radial harmonic numbers $l\geqslant 2$ and the bodies are obliquely directed due to the frequency drift of the modes of l = 1.

Zoom In Zoom Out Reset image size

Moving on to the cases with $\mu =3$ as displayed in Figure 6, one sees that overall the Morlet spectra look similar to the top-hat cases. However, in this case associating the strongest power with the superposition of multiple wavepackets is less clear, because the strongest power is not as closely related to the portion of the group speed curves embedding a local minimum (see, e.g., Figure 6(b)). In addition, modes of l = 2 can now be more readily discerned: when ${\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}}=10$, these modes show up not only in the Morlet spectrum (Figure 6(d)) but also in the temporal evolution (e.g., the interval between 20 and $50\,R/{v}_{\mathrm{Ai}}$ in Figure 6(c)). Nonetheless, because cutoff wavenumbers exist and ${k}_{{\rm{c}},l}$ increases with l, sausage modes of l = 1 remain dominant and account for the appearance of the narrow tails in the Morlet spectra. The broad heads, however, turn out to be a result of the gradual frequency increase in the group speed curves beyond the speed minimum.

Zoom In Zoom Out Reset image size

Let us start with the case where $\mu =1$. The solution to Equation (8) in the outer portion can be shown to have the form

$$\begin{eqnarray*}&&\tilde{\xi }\propto {x}^{-1/2}{W}_{\nu ,1}(x),\end{eqnarray*}$$

where $x=2{mr}$, and ${W}_{\nu ,1}(x)$ is Whittaker's W function (DLMF 2016, Section 13.14). In addition,

$$\begin{eqnarray}&&\nu =\displaystyle \frac{{kR}}{2}\displaystyle \frac{{v}_{\mathrm{ph}}^{2}(1-{\rho }_{{\rm{e}}}/{\rho }_{{\rm{i}}})}{{v}_{\mathrm{Ai}}^{2}\sqrt{1-{v}_{\mathrm{ph}}^{2}/{v}_{\mathrm{Ae}}^{2}}}.\end{eqnarray} \tag{ 21 }$$

With $\tilde{\xi }$ expressible in terms of ${J}_{1}({nr})$ in the inner portion, the DR reads

$$\begin{eqnarray}&&{nR}\displaystyle \frac{{J}_{0}({nR})}{{J}_{1}({nR})}={mR}-\nu +\displaystyle \frac{1}{2}-\displaystyle \frac{{W}_{\nu +\mathrm{1,1}}(2{mR})}{{W}_{\nu ,1}(2{mR})}.\end{eqnarray} \tag{ 22 }$$

We have used the fact that

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dx}}{W}_{\nu ,1}(x)=\left(\displaystyle \frac{1}{2}-\displaystyle \frac{\nu }{x}\right){W}_{\nu ,1}(x)-\displaystyle \frac{{W}_{\nu +\mathrm{1,1}}(x)}{x}.\end{eqnarray} \tag{ 23 }$$

Some approximate results can be found in the limiting cases, where $k\to 0$ or $k\to \infty$. When $k\to 0$, it turns out that $\nu \to (2l+1)/2$ for radial harmonic number l with $l=1,2,3,\cdots$. It then follows from the definition of ν that

$$\begin{eqnarray}&&{v}_{\mathrm{ph}}\approx {v}_{\mathrm{Ae}}\sqrt{1-{\left(\displaystyle \frac{{\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}}-1}{2l+1}\right)}^{2}{({kR})}^{2}}.\end{eqnarray} \tag{ 24 }$$

Consequently,

$$\begin{eqnarray}&&{v}_{\mathrm{gr}}\approx {v}_{\mathrm{Ae}}\left[1-\displaystyle \frac{3}{2}{\left(\displaystyle \frac{{\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}}-1}{2l+1}\right)}^{2}{({kR})}^{2}\right].\end{eqnarray} \tag{ 25 }$$

While the non-existence of cutoff wavenumbers is readily expected since $f(r)$ is asymptotically less steep than ${r}^{-2}$, Equations (24) and (25) offer explicit expressions for what happens when kR is small. When kR approaches infinity, the right-hand side (rhs) of Equation (22) tends to infinity. As happens in the top-hat case, the group and phase speeds for trapped modes of radial harmonic number l can still be approximated by Equations (16) and (17), respectively.

Now consider the case where $\mu =2$. The solution to Equation (8) in the outer portion can be shown to have the form

$$\begin{eqnarray*}&&\tilde{\xi }\propto {K}_{\nu }({mr}),\end{eqnarray*}$$

where K_ν is modified Bessel's function of the second kind with

$$\begin{eqnarray}&&\nu =\sqrt{1-\displaystyle \frac{{v}_{\mathrm{ph}}^{2}}{{v}_{\mathrm{Ai}}^{2}}\left(1-\displaystyle \frac{{\rho }_{{\rm{e}}}}{{\rho }_{{\rm{i}}}}\right){({kR})}^{2}}.\end{eqnarray} \tag{ 26 }$$

Now that the cutoff wavenumbers are given by Equation (18), one sees that $\nu \to 0$ at the cutoff because ${v}_{\mathrm{ph}}={v}_{\mathrm{Ae}}={v}_{\mathrm{Ai}}\sqrt{{\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}}}$. Moving away from this cutoff, ${v}_{\mathrm{ph}}k$ increases with k, meaning that ν is purely imaginary for trapped sausage modes (see DLMF 2016, Section 10.45, for a discussion of K_ν of imaginary order). With $\tilde{\xi }$ expressible in terms of ${J}_{1}({nr})$ in the inner portion, the DR reads

$$\begin{eqnarray}&&({nR})\displaystyle \frac{{J}_{0}({nR})}{{J}_{1}({nR})}=1-\nu -({mR})\displaystyle \frac{{K}_{\nu -1}({mR})}{{K}_{\nu }({mR})}.\end{eqnarray} \tag{ 27 }$$

As is the case for $\mu =1$, the rhs of Equation (27) grows unbounded when kR approaches infinity. This means that once again ${v}_{\mathrm{ph}}$ and ${v}_{\mathrm{gr}}$ can be approximated by Equations (17) and (16), respectively.

Figure 7 presents, in a form identical to Figure 4, the wavenumber dependence of the phase and group speeds as found with BVPSuite. We have also numerically solved the analytical DRs pertinent to $\mu =1$ and 2 (Equations (22) and 27), and found that these solutions agree exactly with the red and green curves. Note that these analytical results yield an asymptotic behavior for ${v}_{\mathrm{ph}}$ and ${v}_{\mathrm{gr}}$ to depend on k in the ways described by Equations (17) and (16), respectively. In fact, this asymptotic behavior holds for arbitrary combinations of $[{\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}},\mu ]$ as found from an extensive set of computations, and partly explains why the ${v}_{\mathrm{gr}}$–k curves are exclusively of the −/+ type for "outer μ" profiles. Figure 7 is a subset of our computations and illustrates that with increasing wavenumber, ${v}_{\mathrm{gr}}$ always decreases first to some minimum before eventually approaching ${v}_{\mathrm{Ai}}$ from below.

Zoom In Zoom Out Reset image size

When it comes to cutoff wavenumbers, one sees exactly the same behavior as in the "μ power" cases. Cutoff wavenumbers do not exist for any radial harmonic number l when $\mu \lt 2$, but do exit and increase with l when $\mu \gt 2$. In the particular case, where $\mu =2$, cutoff wavenumbers exist but do not depend on the radial harmonic number. They are actually exactly given by Equation (18). The behavior of cutoff wavenumbers is expected, because it is determined solely by the asymptotic r-dependence of $f(r)$, which agrees exactly with the "μ power" case: $1/[1+{(r/R)}^{\mu }]\to {(r/R)}^{-\mu }$ when $r\to \infty$.

To sum up here, cutoff wavenumbers are absent (present) when $\mu \lt 2$ ($\mu \geqslant 2$) as happens for the "μ power" cases. However, we note that in this case the ${v}_{\mathrm{gr}}$–k curves are unanimously of the "−/+" type, which is different from the "μ power" cases where $\mu =2$ also separates the "−" from the "−/+" behavior.

Figure 8 presents the temporal evolution and Morlet spectra of the sampled density perturbations for "outer μ" profiles with $\mu =1$. One sees that the Morlet spectra bear similarities to those in both Figures 5 and 6: they look like crazy tadpoles with fins. This is understandable because the characteristics of the ${v}_{\mathrm{gr}}$–k curves in the present case is a mixture of those for $\mu =1$ and $\mu =3$ in the "μ power" cases. To be specific, no cutoff wavenumbers exist for trapped modes of any radial harmonic number l, the consequence being that modes of $l\geqslant 2$ can be excited to account for the fins. On the other hand, the non-monotonical wavenumber dependence of the group speeds explains why the main bodies of the wavelet spectra are in the shape of "crazy tadpoles."

Zoom In Zoom Out Reset image size

Moving on to the cases with $\mu =3$ as given in Figure 9, one sees the classical crazy tadpoles in the Morlet spectra. This behavior agrees closely with the top-hat case, and we need only to mention that the superposition of multiple wavepackets plays a key role in shaping the broad head of the tadpoles.

Zoom In Zoom Out Reset image size

An analytical treatment for this particular μ was initially given by Pneuman (1965) who assumed the ambient to be a vacuum (${\rho }_{{\rm{e}}}=0$). Our Paper I offered a generalization by presenting a DR valid for arbitrary finite ${\rho }_{{\rm{e}}}$ but focused only on the fundamental radial modes. We now further exploit that DR to examine trapped modes of arbitrary radial harmonic number l ($l=1,2,3,\ldots$).

The following definitions are necessary,

$$\begin{eqnarray*}p & = & \displaystyle \frac{\omega R}{{v}_{\mathrm{Ai}}}\sqrt{1-\displaystyle \frac{{\rho }_{{\rm{e}}}}{{\rho }_{{\rm{i}}}}},\\ \alpha & = & 1-\displaystyle \frac{{(\omega R/{v}_{\mathrm{Ai}})}^{2}-{({kR})}^{2}}{4p}.\end{eqnarray*}$$

The solution to Equation (8) in the inner portion can be shown to have the form

$$\begin{eqnarray*}&&\tilde{\xi }\propto {x}^{1/2}\exp (-x/2)M(\alpha ,2,x),\end{eqnarray*}$$

where $x=p{(r/R)}^{2}$, and $M(a,b,x)$ is Kummer's M function (DLMF 2016, Section 13.2). With $\tilde{\xi }$ expressible in terms of ${K}_{1}({mr})$ in the outer portion, the DR then reads

$$\begin{eqnarray}&&-({mR})\displaystyle \frac{{K}_{0}({mR})}{{K}_{1}({mR})}=2-p+\alpha p\displaystyle \frac{M(\alpha +1,3,p)}{M(\alpha ,2,p)}.\end{eqnarray} \tag{ 28 }$$

Approximate expressions can be found for both the cutoff wavenumbers and the asymptotic behavior of the phase and group speeds at large wavenumbers. In fact, both are related to the fact that the DR for any radial harmonic number l can be approximated by $\alpha =1-l$ for nearly the entire range of wavenumbers. Given that ${v}_{\mathrm{ph}}={v}_{\mathrm{Ae}}$ at the cutoff, one sees that

$$\begin{eqnarray*}p & = & {k}_{{\rm{c}}}R\sqrt{\displaystyle \frac{{\rho }_{{\rm{i}}}}{{\rho }_{{\rm{e}}}}-1},\\ \alpha & = & 1-\displaystyle \frac{{k}_{{\rm{c}}}^{2}{R}^{2}({\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}}-1)}{4p}=1-\displaystyle \frac{p}{4}.\end{eqnarray*}$$

This means that the cutoff can be approximated by $p\approx 4l$, or equivalently

$$\begin{eqnarray}&&{k}_{{\rm{c}},l}R\approx \displaystyle \frac{4l}{\sqrt{{\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}}-1}}.\end{eqnarray} \tag{ 29 }$$

It turns out this approximation is increasingly accurate with l, and is accurate to within $13.7 \%$ for l = 1.⁵ On the other hand, α at large wavenumbers is almost exactly $1-l$, yielding a simple equation quadratic in ω. Solving this equation for ω, one finds that

$$\begin{eqnarray*}&&\displaystyle \frac{\omega R}{{v}_{\mathrm{Ai}}}=2l\sqrt{1-\displaystyle \frac{{\rho }_{{\rm{e}}}}{{\rho }_{{\rm{i}}}}}+\sqrt{4{l}^{2}\left(1-\displaystyle \frac{{\rho }_{{\rm{e}}}}{{\rho }_{{\rm{i}}}}\right)+{k}^{2}{R}^{2}}.\end{eqnarray*}$$

Consequently,

$$\begin{eqnarray}&&\displaystyle \frac{{v}_{\mathrm{ph}}}{{v}_{\mathrm{Ai}}}\approx 1+\displaystyle \frac{2l\sqrt{1-{\rho }_{{\rm{e}}}/{\rho }_{{\rm{i}}}}}{{kR}}+\displaystyle \frac{2{l}^{2}(1-{\rho }_{{\rm{e}}}/{\rho }_{{\rm{i}}})}{{({kR})}^{2}},\end{eqnarray} \tag{ 30 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{v}_{\mathrm{gr}}}{{v}_{\mathrm{Ai}}}\approx 1-\displaystyle \frac{2{l}^{2}(1-{\rho }_{{\rm{e}}}/{\rho }_{{\rm{i}}})}{{({kR})}^{2}}.\end{eqnarray} \tag{ 31 }$$

In other words, one expects to see that ${v}_{\mathrm{gr}}$ approaches ${v}_{\mathrm{Ai}}$ from below.

With BVPSuite, we have computed the phase and group speed curves for an extensive set of combinations of $[{\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}},\mu ]$. We find that cutoff wavenumbers (${k}_{{\rm{c}},l}$) exist for arbitrary radial harmonic number l and ${k}_{{\rm{c}},l}$ increases with l. The existence of cutoff wavenumbers is expected because $f(r)$ is identically zero for $r\gt R$ and therefore steeper than ${r}^{-2}$. However, it is interesting to note that ${v}_{\mathrm{gr}}$ at large wavenumbers follows exactly the same k dependence as given by Equation (20) for "μ power" profiles. Not only the exponent β but also the constant c are found to be identical. In particular, for the particular choice $\mu =2$ (and hence $\beta =1$), we find that ${v}_{\mathrm{gr}}^{2}/{v}_{\mathrm{Ai}}^{2}-1\to 0$ as expected from Equation (20). It is just that in the present case ${v}_{\mathrm{gr}}^{2}/{v}_{\mathrm{Ai}}^{2}-1\to {0}^{-}$ (see Equation (31)), in contrast to the "μ power" case where ${v}_{\mathrm{gr}}^{2}/{v}_{\mathrm{Ai}}^{2}-1\to {0}^{+}$. Then why is the asymptotic group speed behavior identical to what happens for "μ power" profiles, barring this minor difference happening for $\mu =2$? To understand this, we note that at large density contrasts, LN15 showed that the dispersive properties of sausage modes are determined by the behavior of the turning points of the function $Q(r)$ as given by Equation (44) therein. For large wavenumbers, the locations of both turning points correspond to an r that is substantially smaller than R. In this portion, however, the "inner μ" profiles are not too different from the "μ power" ones because $1/[1+{(r/R)}^{\mu }]\approx 1-{(r/R)}^{\mu }$. While the analytical analysis in LN15 was conducted in the WKB framework by assuming ${\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}}\gg 1$, our numerical results indicate that actually the asymptotic behavior holds for arbitrary values of ${\rho }_{{\rm{i}}}$/${\rho }_{{\rm{e}}}$ as long as ${\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}}\gt 1$.

Figure 10 presents, in the same format as Figure 4, some examples of the phase (the upper row) and group (lower) speed curves taken from the extensive set of numerical results. Regarding cutoff wavenumbers, one sees that they exist for all the examined μ values and for both density contrasts. In addition, one finds for both density contrasts that $\mu =2$ is a dividing line that separates the cases where ${v}_{\mathrm{gr}}$ eventually approaches ${v}_{\mathrm{Ai}}$ from above from those where ${v}_{\mathrm{gr}}$ approaches ${v}_{\mathrm{Ai}}$ from below. As to the particular case, where $\mu =2$, one sees that ${v}_{\mathrm{gr}}$ tends to ${v}_{\mathrm{Ai}}$ from below, as expected from Equation (31). However, while overall the ${v}_{\mathrm{gr}}$–k curves behave in a "−/+" manner when $\mu \geqslant 2$, some interesting complications arise when $\mu \lt 2$. For a combination $[{\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}},\mu ]$ of $[3,1]$ (the red curves in the left column), ${v}_{\mathrm{gr}}$ is found to be a monotonically decreasing function of k. However, when $[{\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}},\mu ]=[10,1]$, the red curves in the right column indicate that ${v}_{\mathrm{gr}}$ decreases with k first and then increases before eventually decreasing toward ${v}_{\mathrm{Ai}}$. In other words, the ${v}_{\mathrm{gr}}$–k curves are of the "−/+/−" type, possessing not only a local minimum but also a local maximum. If examining the cases where $\mu =1.5$, the blue curves in both columns suggest that the ${v}_{\mathrm{gr}}$–k curves both belong to this "−/+/−" category.

Zoom In Zoom Out Reset image size

Specializing to the density contrasts in the range of $[2,10]$ and μ in the range of $[1,3]$, Figure 11 summarizes the behavior of the group speed curves for "inner μ" profiles. One sees that in this case, cutoff wavenumbers exist regardless of ${\rho }_{{\rm{i}}}$/${\rho }_{{\rm{e}}}$ or μ. While the ${v}_{\mathrm{gr}}$–k curves are unanimously of the "−/+" type for $\mu \geqslant 2$, the area $\mu \lt 2$ is further subdivided into two regions separated by the left dash–dotted curve. For combinations of $[{\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}},\mu ]$ that lie below (above) this curve, the ${v}_{\mathrm{gr}}$–k curves behave in a "−" ("−/+/−") manner.

Zoom In Zoom Out Reset image size

Figure 12 presents the temporal evolution and Morlet spectra of density perturbations for two "inner μ" profiles, both with $\mu =1$. Common to both columns is that the fundamental radial modes dominate the signals. Nonetheless, some substantial difference exists in that the Morlet spectrum for ${\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}}=3$ (the left column) looks like a submarine, whereas it recovers the classical shape of crazy tadpoles when ${\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}}=10$ (right). The former appearance is understandable because the energy contained in the initial perturbation is primarily distributed to the starting portion of the ${v}_{\mathrm{gr}}$–k curves, where ω shows little variation. The appearance of a Morlet spectrum in a tadpole shape in the right column requires some explanation, because for this combination of $[{\rho }_{{\rm{i}}}/{\rho }_{{\rm{e}}},\mu ]$, the ${v}_{\mathrm{gr}}$–k curves are actually of the "−/+/−" type, which is different from the top-hat case. In fact, the broad head can still be attributed to the superposition of multiple wavepackets with group speeds close to the local minimum. Nonetheless, different from the tadpoles in Figures 3(d) and 9(d), the strongest Morlet power in the present tadpole does not extend beyond h/${v}_{\mathrm{Ai}}$. This occurs because the local minimum in the ${v}_{\mathrm{gr}}$–k curve is larger than the internal Alfvén speed.

Zoom In Zoom Out Reset image size

Figure 13 displays what happens when $\mu =3$. One sees the classical crazy tadpoles in the Morlet spectra, which is especially true when the density contrast is chosen to be 10. In fact, these spectra are very similar to those in Figure 3, which is readily understandable given the similarities of the group speed curves in both sets of figures.

Zoom In Zoom Out Reset image size