TY - JOUR
AU - Begelman, Mitchell, C
AB - ABSTRACT We carry out an extensive linear stability analysis of magnetized cylindrical jets in a global framework. Foregoing the commonly invoked force-free limit, we focus on the small-scale, internal instabilities triggered in regions of the jet dominated by a toroidal magnetic field, with a weak vertical field and finite thermal pressure gradient. Such regions are likely to occur far from the jet source and boundaries, and are potential sites of magnetic energy dissipation that is essential to explain the particle acceleration and radiation observed from astrophysical jets. We validate the local stability analysis of Begelman by verifying that the eigenfunctions of the most unstable modes are radially localized. This finding allows us to propose a generic stability criterion in the presence of a weak vertical field. A stronger vertical field with a radial gradient complicates the stability criterion, due to the competition between the destabilizing thermal pressure gradient and stabilizing magnetic pressure gradients. Nevertheless, we argue that the jet interiors generically should be subject to rapidly growing small-scale instabilities, capable of producing current sheets that lead to dissipation. We identify some new instabilities, not predicted by the local analysis, which are sensitive to the background radial profiles but have smaller growth rates than the local instabilities, and discuss the relevance of our work to the findings of recent numerical jet simulations. instabilities, MHD, methods: numerical, stars: jets, galaxies: jets 1 INTRODUCTION Astrophysical jets are bright, collimated, and fast-moving outflows that emanate from a wide variety of systems, including young stellar objects (YSOs), active galactic nuclei (AGNs), microquasars, and gamma-ray bursts (GRBs). Depending upon the source, these jets can be relativistic (e.g. AGNs, microquasars, GRBs) or non-relativistic (e.g. YSOs). They traverse distances much larger than their initial radius, which is typically the size of their source, all the while maintaining a large-scale coherent structure. For instance, the jets from YSOs traverse ≲ 1pc, which is ∼105–107 times the size of the central star, and those from AGNs traverse ≳ 106 pc, which is ∼109 times the size of the central supermassive black hole (see e.g. de Gouveia Dal Pino 2005 and Beskin 2010 for an observational overview). Magnetic fields are believed to play an important role in governing the structure, dynamics, and origin of these jets (e.g. Blandford & Znajek 1977). In spite of the differences in their environments and scales, the jets display very similar overall properties, not all of which are well understood. One of the major open areas of inquiry in this field is the question of jet stability (see e.g. Hardee 2011, 2013 for reviews on this subject). It has long been known from basic plasma physics theory and laboratory experiments that cylindrical, magnetized plasmas are unstable to various kinds of instabilities (see e.g. Bateman 1978 and Freidberg 1982, especially in the context of controlled fusion in tokamaks). In the case of a jet, these instabilities can be either kinetically or magnetohydrodynamically driven. The Kelvin–Helmholtz instability epitomizes the former category, and arises mainly due to the velocity shear created by the bulk motion of the jet with respect to the ambient medium (see Hardee 2013). The latter category is sub-classified as (see e.g. Longaretti 2008) (i) current-driven instabilities, which are driven by the current parallel (j∥) to the total magnetic field B and are most easily studied in the context of cold, pressureless jets, and (ii) pressure-driven instabilities, which are driven by the gradient of the background plasma pressure and the electric current perpendicular (j⊥) to B. The curvature of the magnetic field lines plays an important role in confining the plasma in this case, and instability occurs when the destabilizing plasma pressure gradient is strong enough to push the plasma out of this curvature. If both components of the current are present, along with a finite plasma pressure gradient, then the instabilities will be of a mixed nature (see e.g. Striani et al. 2016). It can be roughly shown that pressure-driven effects dominate if Bz ≪ Bϕ, which yields j⊥ > j∥, and current-driven effects dominate if Bz ≳ Bϕ, which yields j⊥ < j∥, where Bϕ and Bz are the toroidal and vertical magnetic fields, respectively (see e.g. Bonanno & Urpin 2011b). The two most commonly studied magnetohydrodynamic (MHD) instabilities in magnetized plasmas are the m = 0 pinch mode and the m = 1 kink mode, where m is the azimuthal wavenumber. The former leads to a pinching distortion of the jet and triggers the so-called sausage instability, while the latter induces helical perturbations giving rise to the ‘kink instability’ (see e.g. Hardee 1979). Depending on the length-scales at which they operate, the above instabilities can be further classified as internal or external. Internal instabilities are small-scale modes compared to the jet radius, which are confined deep within the jet interior as the perturbations have no displacement at the jet boundaries. By definition these induce local changes (in morphology and energetics) without destroying the jet as a whole. External instabilities, on the other hand, are large-scale modes of the order of the jet radius, such that the perturbations have a non-zero displacement at the jet boundaries. These can potentially disrupt the whole jet. For instance, the external kink instability can bend and wiggle the entire jet off-axis, the external pinch instability can constrict the jet radius at various heights, and the Kelvin–Helmholtz instability due to the velocity shear at the jet surface can lead to mixing and interaction with the ambient medium. The classical approach to understanding instabilities in jets is to perform a linear stability analysis of the ideal MHD equations. Many studies have been carried out using this approach, while invoking different magnetic field profiles, boundary conditions (impenetrable, free, outflow, etc.), methods of solution, and physical assumptions (inclusion/exclusion of rotation, velocity shear, plasma pressure, etc.). Most of the analyses, however, have been carried out under the force-free assumption, where the system is magnetically dominated and the effect of plasma pressure is neglected. These include both non-relativistic (Appl 1996; Appl, Lery & Baty 2000) and relativistic (Istomin & Pariev 1996; Lyubarskii 1999; Tomimatsu, Matsuoka & Takahashi 2001; Narayan, Li & Tchekhovskoy 2009) analyses. The key developments in the linear stability analysis of cylindrical, relativistic, rotating, and force-free jets can be summarized as follows. Istomin & Pariev (1996) showed that stability is achieved if Bz is radially constant, while Lyubarskii (1999) showed that a Bz decreasing radially outwards induces instability. Tomimatsu et al. (2001), however, found a completely different criterion governed by the angular velocity of the magnetic field lines. Narayan et al. (2009) studied the effect of poloidal field curvature along with rigid rotation, with the aim of establishing a link between the various stability criteria. They found that their criterion is a generalization of the Kruskal–Shafranov stability criterion (which, however, does not account for rotation) for the m = 1 kink mode (most relevant for tokamaks; e.g. Freidberg 1982), which states that instability ensues when Bϕ/Bz < (2|$\pi$|Rl)/L, where L is the length of the plasma column and Rl is its cylindrical radius. Recent numerical investigations in this context have been carried out by Mizuno et al. (2012), Bodo et al. (2013), Mizuno, Hardee & Nishikawa (2014), Bodo et al. (2016), Sobacchi, Lyubarsky & Sormani (2017), and Sobacchi & Lyubarsky (2018). Begelman (1998; hereafter B98) was the first to consider the effect of the plasma pressure gradient balancing the magnetic gradients in a regime dominated by a toroidal magnetic field (later works accounting for the plasma pressure include Baty & Keppens 2002; Bonanno & Urpin 2011b; Nalewajko & Begelman 2012; O’Neill, Beckwith & Begelman 2012; Kim et al. 2015). He performed a local study of the linear properties of short-wavelength internal instabilities and proposed a corresponding stability criterion. One of our aims in this work is to compare and validate the solutions from B98’s purely local calculations with those from the linear eigenvalue analysis in our global framework. With the advancement of computational resources, numerical simulations of both non-relativistic and relativistic jets are now possible. These simulations are either local, where a small section of the 3D jet that is comoving with the bulk flow is modelled using periodic boundary conditions along the jet axis (Mizuno et al. 2009; Mizuno et al. 2012; O’Neill et al. 2012; Porth & Komissarov 2015); semi-global, where a pre-existing jet is simulated in a non-periodic computational domain to better study the spatial development of various instabilities (Mizuno et al. 2014; Singh, Mizuno & de Gouveia Dal Pino 2016); or fully global, where the whole 3D jet is simulated starting from its injection/launch to large-scale propagation (McKinney & Blandford 2009; Moll 2009; Mignone et al. 2010; Bromberg & Tchekhovskoy 2016). Numerical simulations are important in order to study the non-linear evolution and saturation of both the internal and the external instabilities, as well as to visualize the morphological changes brought about by them in order to compare with observations. However, they still have limitations, as the results may depend significantly on the numerical solvers employed, resolution, launching mechanisms, etc. (see e.g. O’Neill et al. 2012). Thus, linear analyses are still relevant and essential for interpreting the complex numerical results and for making useful predictions. Stability is a fundamental issue for jet studies in two complementary ways. First, since we know observationally that astrophysical jets survive and propagate for large distances, we can infer that they are somehow stable to external instabilities. Such stability may be achieved due to (i) the jet speed becoming sufficiently supersonic or relativistic that external modes are suppressed or do not have enough time to grow; (ii) rapid sideways expansion of the jet that breaks causal connection across it; (iii) the stabilizing influence of rotation and velocity shear; or (iv) a high density contrast or other properties of the ambient medium that promote stability (see e.g. Hardee 2011 for a discussion and, also, McKinney & Blandford 2009; Porth & Komissarov 2015; Bromberg & Tchekhovskoy 2016). Secondly, some kind of internal instability is needed in order to explain several jet observations, although overall jet stability and coherence are still desirable. In spite of the various different mechanisms proposed (see references in Bromberg & Tchekhovskoy 2016), there is significant debate over what is responsible for the multiwavelength emission from relativistic jets in AGNs, microquasars, and GRBs. It has been shown that internal kink-type instabilities are capable of triggering the dissipation of magnetic energy into thermal energy and radiation (see e.g. Mizuno et al. 2011 and Porth, Komissarov & Keppens 2014 in the context of pulsar wind nebulae and Lyutikov & Blandford 2003 and Bromberg & Tchekhovskoy 2016 in the context of GRBs). This most likely occurs due to the formation of thin current sheets at small scales that facilitate reconnection of oppositely directed magnetic field lines being pushed closer together (see e.g. Kagan et al. 2015 for a review) – as demonstrated by MHD simulations (e.g. Singh et al. 2016) as well as fully kinetic particle-in-cell (PIC) simulations (e.g. Sironi, Petropoulou & Giannios 2015). Also, 3D resistive MHD simulations have shown that magnetized plasma columns can become unstable to both pressure-driven and current-driven instabilities, which leads to the fragmentation of current sheets into filaments, followed by magnetic reconnection at an enhanced rate (Striani et al. 2016). The local dissipation zones created by internal instabilities can also serve as sites for particle acceleration, which is required to explain the observed non-thermal synchrotron and inverse Compton emission from relativistic jets (Sikora et al. 2005). This idea has gained considerable support from test particle acceleration studies that have been carried out in magnetic reconnection domains inside relativistic jets (see e.g. de Gouveia Dal Pino & Kowal 2015 for a review). More recently, efficient particle acceleration has also been demonstrated in both reconnection and turbulence studies using PIC simulations, and the results are in general consistent with jet observations (see e.g. Guo et al. 2016; Petropoulou, Giannios & Sironi 2016; Zhdankin et al. 2017; Werner et al. 2018, and references therein). Furthermore, internal pinch and kink instabilities together can explain the rich small-scale radial structures observed in jets, such as equispaced emission knots and wiggles (de Gouveia Dal Pino 2005; Hardee 2013). A lot of work has been done in trying to understand what perpetuates large-scale jet integrity. However, a fundamental understanding of internal instabilities is equally important, in order to answer basic questions such as how jets radiate. In this work, we perform a comprehensive study of the linear development of internal instabilities and the associated stability criteria, by solving the global eigenvalue problem for a cylindrical magnetized plasma. Most linear studies of internal instabilities are carried out in the force-free limit and/or in conjunction with studies of large-scale external instabilities that require the presence of an ambient medium. In our work, we drop the force-free assumption and focus exclusively on the internal instabilities that are confined deep within the jet interior. We primarily model regions having a dominant toroidal magnetic field, a sub-dominant vertical magnetic field, and a non-zero thermal pressure gradient, which is more prone to trigger pressure-driven instabilities. Most of the magnetic energy dissipation is supposed to take place in such a regime, which is likely to occur significantly away from the jet source. Our work thus expands on B98’s local, analytical calculations by self-consistently including the effects of radial gradients, geometric and magnetic curvatures, and a more complex toroidal field (B98 considered a radial power law). We find that the maximum growth rates of the unstable modes are in excellent accord with the local prediction. More importantly, the B98 instability criterion – which essentially states that when the vertical field is very weak, the thermal pressure gradient is sufficiently strong, and Bϕ falls slower than 1/r, the underlying region is unstable to both axisymmetric and higher-order non-axisymmetric modes – holds true even in a global framework. This validation implies that dissipation in the jet interior is likely dominated by myriad local instabilities, i.e. radially localized modes, which tend to occur at large vertical wavenumbers and have the highest growth rates. Thus, we predict that fast-growing, internal instabilities will always be triggered inside a jet as long as the B98 instability criterion is met locally, irrespective of the exact nature of the background radial profiles. Note that once triggered, the instabilities should have enough time to grow before the jet fluid moves to a significant distance from the source. B98 showed that these instabilities have growth times of the order of Alfvén crossing times, which are short enough in a toroidal-field-dominated region to develop non-linearly. Since in this work we recover growth rates very similar to those predicted by B98, the same conclusion applies. Not surprisingly, we also identify some intrinsically global instabilities, whose characteristics differ somewhat from the local predictions. These usually occur at smaller vertical wavenumbers and have smaller growth rates than the local instabilities. We clarify here that both the local and the global instabilities in this work refer to small-scale, internal instabilities. In our context, ‘local’ implies that the instabilities are insensitive to the background radial curvature and gradients, whereas the global instabilities are sensitive to these factors (i.e. global here is not to be confused with large-scale, external instabilities that we do not talk about in this work). We additionally test the effect of introducing a stronger vertical field and a corresponding radial gradient. Such conditions tend to introduce aforementioned current-driven effects, and the stability criterion becomes more complex due to the interplay between the thermal and magnetic pressure gradients. We also compare our results with recently performed numerical jet simulations, in order to obtain more insight into the workings of real jets. This paper is organized as follows. In Section 2, we lay out the basic properties of the jet model we study, along with the assumptions invoked and the underlying linearized MHD equations governing the problem. In Section 3, we discuss previously carried out analytical stability analyses most relevant to our work. We recall the local dispersion relation and stability criterion proposed by B98, the local stability criterion of Suydam (1958), and the magnetic resonance condition, all of which shall be tested in our global framework. In Section 4, we describe the numerical set-up and normalization scheme used for our global eigenvalue analysis. In Section 5, we present detailed global solutions for different background conditions appropriate for jet interiors, and compare them with the local predictions wherever applicable. In Section 6, we compare our results with those of existing 3D jet simulations, with the aim of better understanding what triggers as well as quenches internal instabilities in jets. We conclude in Section 7, by summarizing our main results and their implications for observed astrophysical jets. 2 THEORY The literature studying MHD instabilities in the context of astrophysical jets is vast and rather confusing. It is indeed a complex topic and only a glimpse of this is given in the Introduction. Keeping this in mind, in this section we outline the main rationale behind our jet model, justify our assumptions, and lay out the basic equations underlying the problem. 2.1 Jet model and assumptions We carry out a linear stability analysis for the following jet model: We assume the jet to have a cylindrical geometry (r, ϕ, z), threaded by an axisymmetric magnetic field B = (0, B0ϕ(r), B0z(r)), which is often referred to as a screw-pinch in the plasma physics literature. Note that we are particularly interested in studying the internal instabilities that may arise in a toroidal magnetic field dominated regime, which are predominantly pressure driven in nature. Such a region may arise due to the following reasons. For a non-relativistic jet, as it propagates away from the source along the z-axis, the poloidal field is expected to fall off faster (∼1/r2) than the toroidal field (∼1/r) in the radial direction due to magnetic flux conservation during lateral expansion (Begelman, Blandford & Rees 1984). In a relativistic jet, a toroidally dominated region can arise if the outer flux surfaces diverge away from the jet core, as demonstrated by Begelman & Li (1994), which can further explain the puzzle behind converting Poynting flux into kinetic energy in a magnetically dominated jet. We shall also briefly discuss the effect of a stronger poloidal field, likely to occur near the jet core, that makes current-driven effects important. In a real jet, it is likely that the complex magnetic field structure, bulk motion, and velocity gradients make it quite hard to distinguish between the various instabilities triggered except at suitable limits (Baty & Keppens 2002). In this work, however, our goal is to isolate the internal instabilities to the extent possible. Our global framework thus comprises an annular region of a static jet bounded by rigid walls (see Section 4.1 for details), which eliminates contamination due to the Kelvin–Helmholtz instability as well as the external kink acting at the jet boundaries (the vertical direction is still assumed to be local). Note that our inner radial boundary does not extend all the way up to the r = 0 singularity (i.e. we do not study the instabilities that may arise in the jet core itself). This restriction does not interfere with our objective of studying the internal instabilities or testing the local calculations. In fact, it allows us to solve the global eigenvalue problem more accurately, i.e. without making any assumptions about the asymptotic behaviour of the eigenfunctions at small radii (e.g. Bodo et al. 2013; Kim et al. 2015). The more commonly invoked force-free equilibrium is quite restrictive and, hence, we invoke the presence of a non-zero thermal pressure gradient in the jet that balances the magnetic forces (following B98). Also, as internal instabilities are triggered in the jet, they are likely to dissipate magnetic energy into thermal energy, which further enhances the role of such a thermal pressure gradient. We perform our analysis in the non-relativistic limit, and assume the jet to be non-rotating as well as static (i.e. the fluid velocity v = 0). In other words, we work in the jet fluid frame for simplicity. Note that we also ignore the effect of a comoving, poloidal velocity shear that may be present in the jet interior (see e.g. Nalewajko & Begelman 2012). 2.2 Basic equations We model our jet by invoking the ideal, non-relativistic MHD equations describing a magnetized, compressible fluid: \begin{eqnarray*} \frac{\partial \rho }{\partial t} + \nabla \cdot (\rho \boldsymbol{ v}) = 0 , \end{eqnarray*} (1) \begin{eqnarray*} \frac{\partial \mathbf {\boldsymbol{ v}}}{\partial t} + (\mathbf {\boldsymbol{ v}} \cdot \nabla) \mathbf {\boldsymbol{ v}} + \frac{1}{\rho }\nabla \left(P + \frac{\mathbf {\boldsymbol{ B}}^2}{8\pi }\right) - \frac{1}{4\pi \rho }(\mathbf {\boldsymbol{ B}}\cdot \nabla)\mathbf {\boldsymbol{ B}} = 0 , \end{eqnarray*} (2) \begin{eqnarray*} \frac{\partial \mathbf {\boldsymbol{ B}}}{\partial t} - \nabla \times (\mathbf {\boldsymbol{ v}}\times \mathbf {\boldsymbol{ B}}) = 0 , \end{eqnarray*} (3) \begin{eqnarray*} \nabla \cdot {\mathbf {\boldsymbol{ B}}} = 0 \quad {\rm or}\quad {\mathbf {\boldsymbol{ B}}} = \nabla \times {\mathbf {\boldsymbol{ A}}} , \end{eqnarray*} (4) where ρ is the density, P is the plasma pressure, and A is the magnetic vector potential. We substitute equation (4) in equation (3) to obtain \begin{eqnarray*} \frac{\partial \mathbf {\boldsymbol{ A}}}{\partial t} - (\mathbf {\boldsymbol{ v}}\times \mathbf {\boldsymbol{ B}}) = 0 , \end{eqnarray*} (5) which is a particularly useful form of the magnetic induction equation when solving for non-axisymmetric perturbations. Recall that B is invariant under the gauge transformation A → A + ∇ψ, where ψ is any scalar function. We also define the pitch |${\cal P}$| of the magnetic field, often invoked during the stability analysis of magnetized jet columns: \begin{eqnarray*} {\cal P} = r \frac{B_z}{B_\phi }. \end{eqnarray*} (6) The magnetic configurations studied in this work involve a wide range of pitch parameters. Note that the vertical equilibrium is trivial as we do not consider any vertical gradients of the background quantities. The radial equilibrium is derived from equation (2), which on assuming an isothermal equation of state |$P_0(r)= \rho _0 (r) c_s^2$|⁠, can be written as \begin{eqnarray*} \partial _r \rho _0 = - \frac{B_{0\phi }}{4\pi c_\mathrm{ s}^2} \left[ \frac{B_{0\phi }}{r} + \partial _r B_{0\phi } \right] - \frac{B_{0z}}{4\pi c_\mathrm{ s}^2}\partial _r B_{0z} , \end{eqnarray*} (7) where cs is the sound speed and the subscript ‘0’ denotes unperturbed background variables (see Appendix A). Note that because of the isothermal approximation, the density and plasma pressure gradient are synonymous in this work. The complete set of linearized MHD equations, obtained by perturbing equations (1), (2), (4), and (5), is given in Appendix A. 3 ANALYTICAL STABILITY ANALYSES In this work, we compare the outcome of our global analysis with the solutions from B98’s local dispersion relation, which was derived for a magnetized, relativistic, non-rotating, compressible, cylindrically symmetric plasma (see equation 3.32 of B98). In this section, we present only the relevant equations and final expression for this dispersion relation, and refer the readers to B98 and Das, Begelman & Lesur (2018) for a detailed derivation. We also discuss the stability criteria proposed by B98 and Suydam (1958) and the magnetic resonance condition. 3.1 B98 local dispersion relation B98 considered a toroidally dominated configuration, having a weak, uniform vertical field, B0z ≪ B0ϕ, and a power-law toroidal field, B0ϕ ∝ rα. The perturbed variables were assumed to have the form \begin{eqnarray*} \exp i(l r + m \phi + k_z z - \omega t), \end{eqnarray*} (8) where ω is the modal frequency, t is the time, and l, m, and kz are the radial, azimuthal, and vertical wavenumbers, respectively. The above form of the perturbations corresponds to the short-wavelength approximation, i.e. lr, kzr ≫ 1, where the amplitudes are treated as constant coefficients. Note that m can take only integer values such that m = 0, 1, 2, etc., of which the m = 0 (pinch) and m = 1 (kink) modes are of particular interest to jet stability, as discussed in the Introduction (the m = 2 mode causes elliptical distortions). m can be positive or negative, which indicates the sense of the spatial twist of the non-axisymmetric mode. The radial equilibrium given by equation (7) for the case considered by B98 simplifies to \begin{eqnarray*} \partial _r \rho _0 = - (1+\alpha)\frac{B_{0\phi }^2}{4\pi c_\mathrm{ s}^2 r}. \end{eqnarray*} (9) Finally, the non-relativistic version of the local dispersion relation obtained by B98 can be written as (see also appendix B of Das et al. 2018) \begin{eqnarray*} &&{\frac{1}{2} \biggl (1 + \frac{l^2}{k_z^2} \biggr) \biggl (1 + \frac{v_{A\phi }^2}{c_\mathrm{ s}^2} \biggr) \omega ^4 - \biggl [ \biggl (1 + \frac{l^2}{k_z^2} \biggr) (m+ \eta)^2 \biggl (1 + \frac{v_{A\phi }^2}{2c_\mathrm{ s}^2} \biggr) }\nonumber\\ &&{\quad \times\,\frac{v_{A\phi }^2}{r^2}- \alpha \biggl (1 + \frac{v_{A\phi }^2}{c_\mathrm{ s}^2} \biggr) \frac{v_{A\phi }^2}{r^2} + \biggl (1 - \frac{v_{A\phi }^2}{c_\mathrm{ s}^2} \biggr) \frac{v_{A\phi }^2}{r^2} \biggr ] \omega ^2} \nonumber \\ &&{\quad+ \, (m + \eta)^2 \frac{v_{A\phi }^2}{r^2} \biggl [ \frac{1}{2} \biggl (1 + \frac{l^2}{k_z^2} \biggr) (m + \eta)^2 \frac{v_{A\phi }^2}{r^2} - (1 + \alpha) \frac{v_{A\phi }^2}{r^2} \biggr ] \!=\! 0,}\nonumber\\ \end{eqnarray*} (10) where all the new variables are defined in Table 1. This is a fourth-degree dispersion relation in ω, obtained after neglecting the fast magnetosonic modes for simplicity. To derive this relation we have assumed, following B98, that the plasma-beta |$\beta \approx 2c_\mathrm{ s}^2/v_{A\phi }^2$|⁠, since vAϕ ≫ vAz . Hence, vAz only appears in the dispersion relation through η. We non-dimensionalize this equation such that all lengths are scaled by the local radial coordinate r0, all wavenumbers by 1/r0, all velocities by cs, all frequencies by cs/r0, and all densities by the local background density ρ0 at r0, which reduces equation (10) to the dimensionless form (see Table 1 for the new notation): \begin{eqnarray*} &&{\frac{1}{2} \biggl (1 + \frac{l^2}{k_z^2} \biggr) d(1 + d) {\tilde{\omega }}^4 - \biggl [ \biggl (1 + \frac{l^2}{k_z^2} \biggr) (m+ \eta)^2 (1 + 2d) }\nonumber\\ &&{\quad- \, \alpha (1 + d) + (d - 1) \biggr ]{\tilde{\omega }}^2}\nonumber\\ &&{\qquad + \, (m + \eta)^2 \biggl [ \frac{1}{2} \biggl (1 + \frac{l^2}{k_z^2} \biggr) (m + \eta)^2 - (1 + \alpha) \biggr ] = 0 .} \end{eqnarray*} (11) Table 1. Summary of the parameters introduced in the local analysis of Section 3.1 and Section 3.2. Parameters Notes |$\lbrace \tilde{k}_z, \tilde{l} \rbrace = r_0 \lbrace k_z, l\rbrace$| Dimensionless wavenumbers |${\tilde{\omega }} = \frac{r_0}{c_\mathrm{ s}} \omega$| Dimensionless modal frequency |$\eta = \tilde{k}_z \frac{B_{0z}}{B_{0\phi }}$| – |$v_{A\phi } = \frac{B_{0\phi }}{\sqrt{4\pi \rho _0}}$| Toroidal Alfvén velocity |$v_{Az} = \frac{B_{0z}}{\sqrt{4\pi \rho _0}}$| Vertical Alfvén velocity |$d = \frac{c_\mathrm{ s}^2}{v_{A\phi }^2}$| – |$\beta = \frac{2 c_\mathrm{ s}^2}{v_{A\phi }^2 + v_{Az}^2} \approx 2d$| Plasma-beta for vAz ≪ vAϕ |$\alpha =\frac{\partial \ln B_{0\phi }}{\partial \ln r}$| – Parameters Notes |$\lbrace \tilde{k}_z, \tilde{l} \rbrace = r_0 \lbrace k_z, l\rbrace$| Dimensionless wavenumbers |${\tilde{\omega }} = \frac{r_0}{c_\mathrm{ s}} \omega$| Dimensionless modal frequency |$\eta = \tilde{k}_z \frac{B_{0z}}{B_{0\phi }}$| – |$v_{A\phi } = \frac{B_{0\phi }}{\sqrt{4\pi \rho _0}}$| Toroidal Alfvén velocity |$v_{Az} = \frac{B_{0z}}{\sqrt{4\pi \rho _0}}$| Vertical Alfvén velocity |$d = \frac{c_\mathrm{ s}^2}{v_{A\phi }^2}$| – |$\beta = \frac{2 c_\mathrm{ s}^2}{v_{A\phi }^2 + v_{Az}^2} \approx 2d$| Plasma-beta for vAz ≪ vAϕ |$\alpha =\frac{\partial \ln B_{0\phi }}{\partial \ln r}$| – View Large Table 1. Summary of the parameters introduced in the local analysis of Section 3.1 and Section 3.2. Parameters Notes |$\lbrace \tilde{k}_z, \tilde{l} \rbrace = r_0 \lbrace k_z, l\rbrace$| Dimensionless wavenumbers |${\tilde{\omega }} = \frac{r_0}{c_\mathrm{ s}} \omega$| Dimensionless modal frequency |$\eta = \tilde{k}_z \frac{B_{0z}}{B_{0\phi }}$| – |$v_{A\phi } = \frac{B_{0\phi }}{\sqrt{4\pi \rho _0}}$| Toroidal Alfvén velocity |$v_{Az} = \frac{B_{0z}}{\sqrt{4\pi \rho _0}}$| Vertical Alfvén velocity |$d = \frac{c_\mathrm{ s}^2}{v_{A\phi }^2}$| – |$\beta = \frac{2 c_\mathrm{ s}^2}{v_{A\phi }^2 + v_{Az}^2} \approx 2d$| Plasma-beta for vAz ≪ vAϕ |$\alpha =\frac{\partial \ln B_{0\phi }}{\partial \ln r}$| – Parameters Notes |$\lbrace \tilde{k}_z, \tilde{l} \rbrace = r_0 \lbrace k_z, l\rbrace$| Dimensionless wavenumbers |${\tilde{\omega }} = \frac{r_0}{c_\mathrm{ s}} \omega$| Dimensionless modal frequency |$\eta = \tilde{k}_z \frac{B_{0z}}{B_{0\phi }}$| – |$v_{A\phi } = \frac{B_{0\phi }}{\sqrt{4\pi \rho _0}}$| Toroidal Alfvén velocity |$v_{Az} = \frac{B_{0z}}{\sqrt{4\pi \rho _0}}$| Vertical Alfvén velocity |$d = \frac{c_\mathrm{ s}^2}{v_{A\phi }^2}$| – |$\beta = \frac{2 c_\mathrm{ s}^2}{v_{A\phi }^2 + v_{Az}^2} \approx 2d$| Plasma-beta for vAz ≪ vAϕ |$\alpha =\frac{\partial \ln B_{0\phi }}{\partial \ln r}$| – View Large Table 2. Summary of the parameters introduced in the global analysis of Section 4 and Section 5. Parameters Definitions |$r \equiv \frac{r}{R_{\rm in}}$| Dimensionless radius kz ≡ kzRin Dimensionless vertical wavenumber |$\rho \equiv \frac{\rho }{\rho _{\rm in}}$| Dimensionless density |$\omega _{I,R} \equiv \omega _{I,R} \frac{R_{\rm in}}{c_\mathrm{ s}}$| Dimensionless eigenvalues (I ≡ imaginary; R ≡ real) |$v_{1r,1\phi ,1z} \equiv \frac{v_{1r,1\phi ,1z}}{c_\mathrm{ s}}$| Dimensionless perturbed velocities |$v_{A\phi ,Az} \equiv \frac{B_\phi ,z}{\sqrt{4\pi \rho _{\rm in}}} \frac{1}{c_\mathrm{ s}}$| Dimensionless background Alfvén velocities |$v_{A\phi0 ,Az0}$| Background Alfvén velocities at r = Rin |$\alpha =\frac{\partial \ln v_{A\phi }}{\partial \ln r}$| Dimensionless parameter in vAϕ-PL models |$a_1$| Dimensionless parameter in vAϕ-Gen models ϵ Dimensionless parameter in vAz-Var models; |$\frac{\mathrm{ d}P}{\mathrm{ d}r} = \frac{\epsilon }{1-\epsilon } \frac{\mathrm{ d}}{\mathrm{ d}r} \biggl (\frac{B_z^2}{8\pi }\biggr)$| Parameters Definitions |$r \equiv \frac{r}{R_{\rm in}}$| Dimensionless radius kz ≡ kzRin Dimensionless vertical wavenumber |$\rho \equiv \frac{\rho }{\rho _{\rm in}}$| Dimensionless density |$\omega _{I,R} \equiv \omega _{I,R} \frac{R_{\rm in}}{c_\mathrm{ s}}$| Dimensionless eigenvalues (I ≡ imaginary; R ≡ real) |$v_{1r,1\phi ,1z} \equiv \frac{v_{1r,1\phi ,1z}}{c_\mathrm{ s}}$| Dimensionless perturbed velocities |$v_{A\phi ,Az} \equiv \frac{B_\phi ,z}{\sqrt{4\pi \rho _{\rm in}}} \frac{1}{c_\mathrm{ s}}$| Dimensionless background Alfvén velocities |$v_{A\phi0 ,Az0}$| Background Alfvén velocities at r = Rin |$\alpha =\frac{\partial \ln v_{A\phi }}{\partial \ln r}$| Dimensionless parameter in vAϕ-PL models |$a_1$| Dimensionless parameter in vAϕ-Gen models ϵ Dimensionless parameter in vAz-Var models; |$\frac{\mathrm{ d}P}{\mathrm{ d}r} = \frac{\epsilon }{1-\epsilon } \frac{\mathrm{ d}}{\mathrm{ d}r} \biggl (\frac{B_z^2}{8\pi }\biggr)$| View Large Table 2. Summary of the parameters introduced in the global analysis of Section 4 and Section 5. Parameters Definitions |$r \equiv \frac{r}{R_{\rm in}}$| Dimensionless radius kz ≡ kzRin Dimensionless vertical wavenumber |$\rho \equiv \frac{\rho }{\rho _{\rm in}}$| Dimensionless density |$\omega _{I,R} \equiv \omega _{I,R} \frac{R_{\rm in}}{c_\mathrm{ s}}$| Dimensionless eigenvalues (I ≡ imaginary; R ≡ real) |$v_{1r,1\phi ,1z} \equiv \frac{v_{1r,1\phi ,1z}}{c_\mathrm{ s}}$| Dimensionless perturbed velocities |$v_{A\phi ,Az} \equiv \frac{B_\phi ,z}{\sqrt{4\pi \rho _{\rm in}}} \frac{1}{c_\mathrm{ s}}$| Dimensionless background Alfvén velocities |$v_{A\phi0 ,Az0}$| Background Alfvén velocities at r = Rin |$\alpha =\frac{\partial \ln v_{A\phi }}{\partial \ln r}$| Dimensionless parameter in vAϕ-PL models |$a_1$| Dimensionless parameter in vAϕ-Gen models ϵ Dimensionless parameter in vAz-Var models; |$\frac{\mathrm{ d}P}{\mathrm{ d}r} = \frac{\epsilon }{1-\epsilon } \frac{\mathrm{ d}}{\mathrm{ d}r} \biggl (\frac{B_z^2}{8\pi }\biggr)$| Parameters Definitions |$r \equiv \frac{r}{R_{\rm in}}$| Dimensionless radius kz ≡ kzRin Dimensionless vertical wavenumber |$\rho \equiv \frac{\rho }{\rho _{\rm in}}$| Dimensionless density |$\omega _{I,R} \equiv \omega _{I,R} \frac{R_{\rm in}}{c_\mathrm{ s}}$| Dimensionless eigenvalues (I ≡ imaginary; R ≡ real) |$v_{1r,1\phi ,1z} \equiv \frac{v_{1r,1\phi ,1z}}{c_\mathrm{ s}}$| Dimensionless perturbed velocities |$v_{A\phi ,Az} \equiv \frac{B_\phi ,z}{\sqrt{4\pi \rho _{\rm in}}} \frac{1}{c_\mathrm{ s}}$| Dimensionless background Alfvén velocities |$v_{A\phi0 ,Az0}$| Background Alfvén velocities at r = Rin |$\alpha =\frac{\partial \ln v_{A\phi }}{\partial \ln r}$| Dimensionless parameter in vAϕ-PL models |$a_1$| Dimensionless parameter in vAϕ-Gen models ϵ Dimensionless parameter in vAz-Var models; |$\frac{\mathrm{ d}P}{\mathrm{ d}r} = \frac{\epsilon }{1-\epsilon } \frac{\mathrm{ d}}{\mathrm{ d}r} \biggl (\frac{B_z^2}{8\pi }\biggr)$| View Large 3.2 B98 local stability criterion B98 found that relativistic effects have little effect on the growth rate of the unstable modes and, hence, his findings are applicable in the non-relativistic formalism as well. Seeking unstable modes such that ω2 < 0 and ωI > 0 (as per the convention in equation 8), B98 obtained a necessary and sufficient condition for instability to occur (see equation 4.2 of B98) such that \begin{eqnarray*} \frac{\partial \ln B_{0\phi }}{\partial \ln r} \equiv \alpha \gt \frac{1}{2} (m + \eta)^2 - 1 \end{eqnarray*} (12) or, equivalently, unstable modes exist for \begin{eqnarray*} 0\lt (m+\eta)^2\lt 2(1+\alpha). \end{eqnarray*} (13) The above inequality clearly indicates that α = −1 is the case of marginal stability for all m or, equivalently, any local patch of the jet in which B0ϕ decreases faster than 1/r should be stable with respect to both axisymmetric and higher-order non-axisymmetric modes. From equation (12) we can also say that for a given m and α > −1, instability occurs in the range \begin{eqnarray*} -\sqrt{2(1+\alpha)} - m \lt \eta \lt \sqrt{2(1+\alpha)} - m \end{eqnarray*} (14) or \begin{eqnarray*} (-\sqrt{2(1+\alpha)} - m)\frac{B_{0\phi }}{B_{0z}} \lt \tilde{k}_{z}|_{\rm unst} \lt (\sqrt{2(1+\alpha)} - m) \frac{B_{0\phi }}{B_{0z}}. \end{eqnarray*} (15) Apart from the above instability criterion, B98 made certain estimates from his analytical calculations, which we quote here in the non-relativistic limit. First, B98 showed that the fastest-growing modes (i.e. the modes having the most negative ω2) occur in the limit |$l^2/k_z^2 \ll 1$|⁠. Thus, putting |$l^2/k_z^2 \approx 0$|⁠, followed by minimizing ω2 with respect to (m +η)2 in the dispersion relation given by equation (11), one obtains \begin{eqnarray*} {\tilde{\omega }}^2_{\rm min} = - 8 d \biggl [ \biggl (1 - \frac{1+\alpha }{4d}\biggr) - \biggl (1 - \frac{1+\alpha }{2d} \biggr)^{1/2} \biggr ] , \end{eqnarray*} (16) which is equivalent to equation (4.8) of B98, when the relativity parameter e therein is assumed to be e ≪ 1. Thus, the maximum growth rate of the system is given by \begin{eqnarray*} {\tilde{\omega }}_{g}|_{\rm max}= \sqrt{8 d \biggl [ \biggl (1 - \frac{1+\alpha }{4d}\biggr) - \biggl (1 - \frac{1+\alpha }{2d} \biggr)^{1/2} \biggr ]} , \end{eqnarray*} (17) which occurs at \begin{eqnarray*} (m + \eta)^2_{\rm min} = 1 + \alpha + \frac{1}{2}(1 + 2d) {\tilde{\omega }}^2_{\rm min} , \end{eqnarray*} (18) which in turn is equivalent to equation (4.9) of B98 for e ≪ 1. Alternatively, the vertical wavenumbers corresponding to the most unstable modes are given by \begin{eqnarray*} \tilde{k}_{z}|_{\rm max} = \biggl (\pm \sqrt{1 + \alpha + \frac{1}{2}(1 + 2d) {\tilde{\omega }}^2_{\rm min}} -m \biggr) \frac{B_{0\phi }}{B_{0z}}. \end{eqnarray*} (19) We note from equation (17) that the maximum growth rate is independent of m, which is a characteristic of pressure-driven instabilities. This is in contrast to current-driven instabilities, whose maximum growth rate is inversely proportional to m (e.g. Appl et al. 2000). 3.3 Magnetic resonance condition For non-axisymmetric perturbations in a cylindrical plasma column having the form given by equation (8), the magnetic resonance condition is defined as k · B = 0, where k is the wave-vector with components (l, m, kz). Consequently, one can define a resonant surface r = rres, where this condition is satisfied such that \begin{eqnarray*} k_z B_z (r_{\rm res}) + \frac{m}{r_{\rm res}} B_\phi = 0. \end{eqnarray*} (20) The restoring force due to the magnetic tension is supposed to be minimum at these surfaces, which are hence more likely to trigger instabilities (see e.g. Longaretti 2008). By studying the global eigenfunctions of the most unstable modes, we shall test whether the above condition influences the instabilities in this work (see Section 5). 3.4 Suydam’s local stability criterion Suydam (1958) proposed a necessary condition for stability for the non-axisymmetric (m ≠ 0) modes of a screw-pinch such that \begin{eqnarray*} \frac{r}{4} \biggl (\frac{1}{{\cal P}}\frac{\mathrm{ d} {\cal P}}{\mathrm{ d}r} \biggr)^2 + \frac{8\pi }{B_z^2} \frac{\mathrm{ d}P}{\mathrm{ d}r} \ge 0 , \end{eqnarray*} (21) where |${\cal P}$| is the magnetic pitch given by equation (6). Note that this criterion was derived from a variational energy principle, which states that an equilibrium is unstable if the radial displacement of a mode makes the change in potential energy negative. Furthermore, this is a local stability criterion as the displacement of the concerned (m, kz) modes is assumed to be localized only about the magnetic resonance surface (see equation 20), for reasons cited in Section 3.3. The thermal pressure gradient (second term in the inequality 21) has a negative sign and is known to drive the system to instability. Thus, Suydam’s criterion basically states that the magnetic shear (first term in the inequality 21) must be strong enough to counteract this effect to maintain stability. We shall revisit the applicability of this criterion when discussing our global solutions in Section 5. We note that although Sudyam’s criterion is local, it conveys the crucial point that a competition between the thermal and magnetic pressure gradients plays an important role in determining the stability of the system. 4 GLOBAL EIGENVALUE ANALYSIS 4.1 Numerical set-up In order to construct global solutions from the linearized axisymmetric system of equations (A6)–(A14), we consider the jet to be extended in radius r ∈ [Rin, Rout], to account for the curvature terms. As mentioned before, since we are only interested in the internal instabilities characterizing the jet, we assume this domain to be located far from the influence of the ambient medium as well as from the central axis (i.e. Rin ≠ 0). We solve equations (A6)–(A14) as a linear eigenvalue problem using the pseudo-spectral code dedalus (Burns et al., in preparation).1 The solutions are decomposed on the basis of Chebyshev polynomials along the radial grid and on a Fourier basis in the vertical direction such that the perturbed variables have a form \begin{eqnarray*} f(r)\exp i(k_z z + m\phi -\omega t) , \end{eqnarray*} (22) and the linear eigenvalue problem becomes \begin{eqnarray*} {\cal M}(r) {\boldsymbol{\xi }}(r) = \omega {\cal I} {\boldsymbol{\xi }}(r) , \end{eqnarray*} (23) where ω = ωR + iωI is the complex eigenvalue; |${\cal I}$| is the identity operator; |${\cal M}(r)$| is the MHD linear differential operator; and |${\boldsymbol{\xi }}(r)$| is the eigenfunction constituted of the perturbed variables in the system. The growth rate of an unstable mode is given by ωI > 0, whereas ωR is an indicator of the overstability of the mode. We impose rigid/impenetrable boundary conditions at our inner and outer jet boundaries, such that (also, see Das et al. 2018) \begin{eqnarray*} v_{1r} = 0 \quad{\rm at} \quad r = R_{\rm in}, R_{\rm out} \end{eqnarray*} \begin{eqnarray*} \partial _r B_{1z} = 0 \quad {\rm at}\quad r = R_{\rm in}, R_{\rm out} , \end{eqnarray*} (24) where the subscript ‘1’ indicates perturbed variables (see Appendix A). These four boundary conditions suffice for the case with the Weyl gauge; however, for the Coulomb gauge we additionally impose \begin{eqnarray*} A_{1r} = 0 \quad {\rm at}\quad r = R_{\rm in}, R_{\rm out} \end{eqnarray*} (25) to close the system of equations. We also impose boundary conditions on the background pressure/density, which we shall describe in Section 5. We choose our fiducial resolution to be either Nr = 200 or 256, where Nr is the number of radial grid points. The fiducial radial extent of the jet is set to be r ∈ [1, 5]Rin, unless otherwise mentioned. 4.2 Normalization scheme We non-dimensionalize our linearized equations (A6)–(A14) for the global analysis by using quantities at the inner radial jet boundary Rin (unless otherwise mentioned). We scale all length-scales by Rin, all wavenumbers by 1/Rin, all velocities by the sound speed cs, all frequencies by cs/Rin and all densities by the inner radial density ρin. In order to solve the linearized equations, we replace |$\partial _z \rightarrow ik_z$| and |$\partial _t \rightarrow -i\omega$|⁠, and define dimensionless perturbed Alfvén velocities as \begin{eqnarray*} v_{A1 r,A1\phi ,A1z} \equiv \frac{B_{1r,1\phi ,1z}}{\sqrt{4 \pi \rho _{\rm in}}} \biggl (\frac{1}{c_\mathrm{ s}} \biggr) , \end{eqnarray*} (26) dimensionless perturbed vector potential components as \begin{eqnarray*} A_{1 r,\phi ,z} \equiv \frac{A_{1r,\phi ,z}}{\sqrt{4 \pi \rho _{\rm in}}} \biggl (\frac{1}{c_\mathrm{ s} R_{\rm in}} \biggr) , \end{eqnarray*} (27) and dimensionless background Alfvén velocities as \begin{eqnarray*} v_{A \phi , Az} \equiv \frac{B_{\phi ,z}}{\sqrt{4 \pi \rho _{\rm in}}} \biggl (\frac{1}{c_\mathrm{ s}} \biggr) . \end{eqnarray*} (28) Note that from here onwards we will describe the magnetic field in terms of the definition given by equation (28) (see also Table 2 for a summary of the definitions to be used for the global analysis). Finally, depending on the gauge condition imposed, we solve for the set of eigenfunctions |${\boldsymbol{\xi }}(r)|_{\mathrm{ Weyl}}= \lbrace \rho _1, v_{1r}, v_{1\phi }, v_{1z}, v_{A1\phi }, v_{A1z}, A_{1r}, A_{1\phi }, A_{1z}\rbrace$| or |${\boldsymbol{\xi }}(r)|_{\mathrm{ Coulomb}} =\lbrace \rho _1, v_{1r}, v_{1\phi }, v_{1z}, v_{A1\phi }, v_{A1z}, A_{1r}, A_{1\phi }, A_{1z}, \psi \rbrace$|⁠. 5 GLOBAL SOLUTIONS In this section, we present the results of our global stability analysis. We divide the results into three categories depending upon the transverse magnetic field profiles adopted for the jet interior. There is a lot of uncertainty regarding the nature of magnetic fields inside jets, which has led different groups to adopt a wide range of field prescriptions (see references in the Introduction). However, our analysis reveals certain fundamental properties of the internal instabilities, which stand out irrespective of the complex nature of the magnetic field imposed. 5.1 Results for a power-law toroidal field and constant vertical field: |$\boldsymbol{ v}$|Aϕ-PL−|$\boldsymbol{ v}$|Az-Cons We begin by choosing a power law (PL) for the toroidal field, vAϕ = vAϕ0rα, and a constant (Cons) vertical field, vAz = vAz0, where vAϕ0,Az0 are the respective quantities at the inner jet radius (r = 1). This case provides a direct comparison with the local analysis of B98 and allows us to test the validity of the stability criterion proposed therein. The magnetic pitch for this case follows from equation (6) such that \begin{eqnarray*} {\cal P} = r^{1-\alpha } \frac{v_{Az0}}{v_{A\phi 0}} . \end{eqnarray*} (29) Thus, the pitch is an increasing function of radius if α < 1, a decreasing function if α > 1, and a constant if α = 1. All the runs carried out for this case are summarized in Table 3. Table 3. Summary of runs for the case vAϕ-PL–vAz-Cons presented in Section 5.1. m α vAϕ0 vAz0/vAϕ0 Nr Boundary Conditions (B.C.s) imposed on ρ(r) Other notes 0,1,2 −0.3 0.926 0.005 256 ρ(1) = 1; ρ(r → ∞) = 0 Fiducial α < 0 model; |$v_{A\phi 0}= \sqrt{\frac{-2\alpha }{(1+\alpha)}}$| 0 −0.1 0.6 0.008 300 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.3 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.7 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.9 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.5 0.2 0.025 300 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.5 0.6 0.008 300 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.5 1.0 0.005 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.5 1.25 0.004 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.5 1.5 0.003 300 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.1 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.25 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.3 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.4 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.5 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.6 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.7 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.9 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.5 |$\sqrt{2}$| 0.005 256 ρ(1) = 1 ρ(r → ∞) = 0 Pressure equipartition; |$v_{A\phi 0}= \sqrt{\frac{-2\alpha }{(1+\alpha)}}$| 1 −0.3 0.926 0.001 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 0.01 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 0.1 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 0.5 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 0.75 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 1.0 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 0.005 256 ρ(r → ∞) = 10 Constant density background 1 −0.3 0.926 0.005 256 ρ(r → ∞) = 100 Constant density background m α vAϕ0 vAz0/vAϕ0 Nr Boundary Conditions (B.C.s) imposed on ρ(r) Other notes 0,1,2 −0.3 0.926 0.005 256 ρ(1) = 1; ρ(r → ∞) = 0 Fiducial α < 0 model; |$v_{A\phi 0}= \sqrt{\frac{-2\alpha }{(1+\alpha)}}$| 0 −0.1 0.6 0.008 300 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.3 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.7 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.9 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.5 0.2 0.025 300 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.5 0.6 0.008 300 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.5 1.0 0.005 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.5 1.25 0.004 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.5 1.5 0.003 300 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.1 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.25 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.3 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.4 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.5 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.6 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.7 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.9 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.5 |$\sqrt{2}$| 0.005 256 ρ(1) = 1 ρ(r → ∞) = 0 Pressure equipartition; |$v_{A\phi 0}= \sqrt{\frac{-2\alpha }{(1+\alpha)}}$| 1 −0.3 0.926 0.001 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 0.01 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 0.1 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 0.5 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 0.75 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 1.0 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 0.005 256 ρ(r → ∞) = 10 Constant density background 1 −0.3 0.926 0.005 256 ρ(r → ∞) = 100 Constant density background View Large Table 3. Summary of runs for the case vAϕ-PL–vAz-Cons presented in Section 5.1. m α vAϕ0 vAz0/vAϕ0 Nr Boundary Conditions (B.C.s) imposed on ρ(r) Other notes 0,1,2 −0.3 0.926 0.005 256 ρ(1) = 1; ρ(r → ∞) = 0 Fiducial α < 0 model; |$v_{A\phi 0}= \sqrt{\frac{-2\alpha }{(1+\alpha)}}$| 0 −0.1 0.6 0.008 300 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.3 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.7 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.9 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.5 0.2 0.025 300 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.5 0.6 0.008 300 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.5 1.0 0.005 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.5 1.25 0.004 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.5 1.5 0.003 300 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.1 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.25 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.3 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.4 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.5 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.6 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.7 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.9 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.5 |$\sqrt{2}$| 0.005 256 ρ(1) = 1 ρ(r → ∞) = 0 Pressure equipartition; |$v_{A\phi 0}= \sqrt{\frac{-2\alpha }{(1+\alpha)}}$| 1 −0.3 0.926 0.001 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 0.01 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 0.1 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 0.5 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 0.75 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 1.0 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 0.005 256 ρ(r → ∞) = 10 Constant density background 1 −0.3 0.926 0.005 256 ρ(r → ∞) = 100 Constant density background m α vAϕ0 vAz0/vAϕ0 Nr Boundary Conditions (B.C.s) imposed on ρ(r) Other notes 0,1,2 −0.3 0.926 0.005 256 ρ(1) = 1; ρ(r → ∞) = 0 Fiducial α < 0 model; |$v_{A\phi 0}= \sqrt{\frac{-2\alpha }{(1+\alpha)}}$| 0 −0.1 0.6 0.008 300 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.3 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.7 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.9 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.5 0.2 0.025 300 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.5 0.6 0.008 300 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.5 1.0 0.005 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.5 1.25 0.004 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 0 −0.5 1.5 0.003 300 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.1 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.25 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.3 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.4 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.5 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.6 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.7 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.9 0.6 0.008 200 ρ(1) = 1 ρ > 0 ∀ r ∈ [1, 5] 1 −0.5 |$\sqrt{2}$| 0.005 256 ρ(1) = 1 ρ(r → ∞) = 0 Pressure equipartition; |$v_{A\phi 0}= \sqrt{\frac{-2\alpha }{(1+\alpha)}}$| 1 −0.3 0.926 0.001 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 0.01 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 0.1 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 0.5 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 0.75 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 1.0 200 ρ(1) = 1; ρ(r → ∞) = 0 – 1 −0.3 0.926 0.005 256 ρ(r → ∞) = 10 Constant density background 1 −0.3 0.926 0.005 256 ρ(r → ∞) = 100 Constant density background View Large The background density for this case is obtained by integrating equation (7). Keeping our normalization scheme in mind (see Table 2), this yields \begin{eqnarray*} \rho (r) = - \frac{(1+\alpha)v_{A\phi 0}^2}{2\alpha }r^{2\alpha } + \rho _{\rm cons} , \end{eqnarray*} (30) where ρcons is the integration constant. Now, in order to determine ρcons, we impose the inner boundary condition \begin{eqnarray*} \rho =1 \quad{\rm at}\quad r=1 , \end{eqnarray*} (31) which yields \begin{eqnarray*} \rho (r) = 1 + \frac{(1+\alpha)v_{A\phi 0}^2}{2\alpha }(1-r^{2\alpha }) . \end{eqnarray*} (32) The plasma-beta can be written as \begin{eqnarray*} \beta = \frac{2\rho (r)}{v_{A\phi 0}^2 r^{2\alpha } + v_{Az0}^2} . \end{eqnarray*} (33) However, one needs to keep in mind that if vAϕ0 becomes too large, and/or α assumes a large positive value, then the density may become negative beyond a certain radius. This puts a restriction on the computational domain for certain parameter values. If α < 0 on the other hand, then an additional outer boundary condition can be imposed such that \begin{eqnarray*} \rho \rightarrow 0 \quad {\rm as}\quad r \rightarrow \infty . \end{eqnarray*} (34) The above condition sets ρcons = 0 and self-consistently assures ρ > 0 throughout the domain, yielding \begin{eqnarray*} \rho (r) = r^{2\alpha } \quad {\rm and} \quad v_{A\phi 0} = \sqrt{\frac{-2\alpha }{(1+\alpha)}} . \end{eqnarray*} (35) Furthermore, on assuming vAz0 ≪ vAϕ, equation (33) reduces to \begin{eqnarray*} \beta \approx -\frac{1+\alpha }{\alpha } = {\rm constant} . \end{eqnarray*} (36) Our fiducial model corresponds to α = −0.3, vAϕ0 = 0.926, vAz0 = 0.005, and β ≈ 2.33. For this model, we impose both conditions (31) and (34) on the background density. However, in this case, vAϕ0 depends on α and is no longer a free parameter (see equation 35). Thus, for some cases (e.g. when varying α for a given vAϕ0, when α > 0, etc.), we impose only condition (31). We use r ∈ [1, 5] for all the cases listed in Table 3, which also assures ρ > 0 throughout the domain. We shall now discuss the properties of the axisymmetric m = 0, and the non-axisymmetric m = 1 and 2 modes, corresponding to the cases listed in Table 3, separately. 5.1.1 Axisymmetric m = 0 modes The main findings of this subsection, listed below, are summarized by Figs 1, 2, and 3: The left-hand panel of Fig. 1 shows the imaginary part ωI of the eigenvalues, or growth rate, as a function of the vertical wavenumber kz for the fiducial α = −0.3 case (see Table 3). The right-hand panel shows the real part ωR of the corresponding eigenvalues. We also demarcate the fast magnetosonic, Alfvén, and slow magnetosonic modes in the right-hand panel, indicated by the letters F, A, and S, respectively. Note that the total number of modes in the global problem is proportional to the radial resolution. The cyan (ωI ≠ 0) and magenta (ωI = 0) lines in this figure represent the corresponding local solutions predicted by the B98 dispersion relation, which are given by equation (11) with |$l^2/k_z^2=0$|⁠. From the left-hand panel of Fig. 1, we see that the m = 0 modes are symmetric about the kz-axis (see also equation 11), with ωI = 0 at kz = 0. We also find that the ωI > 0 cyan lines are in very good agreement with the growth rates of the most unstable modes of the global eigenvalue solution. Some of the measurable quantities predicted by the local dispersion relation for this case are the maximum growth rate ωg|max ≈ 0.35 (see equation 17), the corresponding vertical wavenumbers kz|max ≈ ±129 (see equation 19), and the range of instability −219 ≲ kz|unst ≲ 219 (see relation 15). All the corresponding global quantities have comparable but slightly lower values as seen from Fig. 1: ωg|max ≈ 0.33, kz|max ≈ ±124, and −206 ≲ kz|unst ≲ 206. The growth rates of the most unstable modes in the range 0 ≤ kz ≲ 60 are also slightly larger than the local prediction [see point (v) below]. On the right-hand panel of Fig. 1, note that the cyan line indicates ωR corresponding to the most unstable modes of the left-hand panel (ωI > 0), while the magenta lines indicate the ωR corresponding to the stable modes (ωI = 0 and ωI < 0). From the global solutions we find that the most unstable modes are purely growing with zero phase velocity ωR = 0, while the stable modes are purely oscillatory, corroborating the local predictions. In fact, the unstable modes appear to be destabilized slow modes. Note that the (stable) fast modes do not have a local counterpart as they are excluded from the local analysis (which yields a fourth-degree dispersion relation instead of sixth; see equation 11). A subtle difference between the local and global solutions is as follows: In the global solutions, as kz → 0, ωR → 0 for all the Alfvén modes, unlike the local solutions, which have ωR ≠ 0 for kz = 0. Also, for |kz| > 160 and |ωR| > 1.25 (where the black lines intersect and then diverge from the magenta lines), additional modes (possibly Alfvénic) appear in the global solutions that do not have a local counterpart. We next discuss Fig. 2, which shows the eigenfunction v1r as a function of radius for the most unstable mode at different kz corresponding to Fig. 1. Note that the eigenfunctions are normalized with respect to their respective maximum amplitudes. Let us define two characteristics of the radial eigenfunctions, which will help us better understand the nature of the instabilities in this work. These are Δr = r2 − r1, which is the full-width at half maximum of the eigenfunction, and Rmax, which gives the peak position of the eigenfunction as measured from the inner boundary r = 1 (demonstrated with kz = 8 in the top left-hand panel of Fig. 2). Based on these, we can broadly classify any instability in this work either as (1) Local or (2) Global. A local instability has very narrow eigenfunctions (i.e. Δr ≲ 0.25), which are highly concentrated (or localized) towards the inner boundary (i.e. Rmax → 1). These properties indicate that the instability does not care about the radial variation of the background quantities like magnetic field and density and, hence, can be fully described by a local calculation (i.e. radially localized). A global instability, on the other hand, is affected by the background radial profiles and, hence, can only be captured in a global framework (i.e. radially extended). These can be of three kinds, having eigenfunctions that are either wide and peak significantly away from the inner boundary (Global I, i.e. Δr > 0.25 and Rmax↛1); or wide and peak towards the inner boundary (Global II, i.e. Δr > 0.25 and Rmax → 1); or narrow and peak significantly away from the inner boundary (Global III, i.e. Δr ≲ 0.25 and Rmax↛1). Table 4 summarizes the characteristics of the various instabilities. We emphasize that although we distinguish between local and global instabilities here, they are both still internal, small-scale instabilities, which do not affect the large-scale properties of the jet. Keeping the above characteristics in mind, we can now look at the different eigenfunctions in Fig. 2. The top left-hand panel shows the eigenfunctions for 8 ≤ kz ≤ 55, the top right-hand panel for 85 ≤ kz ≤ 205, the bottom left-hand panel for −55 ≤ kz ≤ −8, and the bottom right-hand panel for −205 ≤ kz ≤ −85 (the bottom panels thus representing the kz < 0 counterparts of the respective top panels). Note that the radial axes of the right-hand panels of Fig. 2 have been zoomed close to the inner boundary r ∈ [1, 1.5] for clarity, whereas the left-hand panels show the full radial domain r ∈ [1, 5]. We see that the eigenfunctions in the bottom panels are identical to their top counterparts (except for a 180° phase shift displayed by some), which again reflects the symmetry of the m = 0 modes about the kz-axis. The most unstable modes are found to be nodeless in general, across the entire range of unstable kz. The modes shown in the left-hand panels of Fig. 2 have smaller |kz| than those in the right-hand panel. They represent Global I and II instabilities. The eigenfunctions span a large radial extent, and as |kz| → 0, the modes seem to peak towards the outer boundary (Global I), gradually shifting towards the inner boundary as |kz| increases (Global II). This explains why the global solution in Fig. 1 deviates from the local prediction as kz → 0. On the other hand, the larger |kz| modes in the right-hand panels are Local instabilities, having narrow eigenfunctions that are highly localized close to the inner boundary. In general, the modes become more localized towards r = 1 as|kz| increases. In Fig. 3, we study some interesting trends displayed by the m = 0 modes. The left-hand panel shows the evolution of the most unstable m = 0 modes for a fixed vAϕ0 = 0.6 (and |$\beta _0 \approx 2/v_{A\phi 0}^2 = 5.55$|⁠), when α is varied. The dotted and dashed lines represent the global and local solutions, respectively, which are in fairly good agreement. We find that as α decreases, the growth rates of the modes (in both local and global solutions) keep decreasing till they are completely stabilized at α = −1, thus validating the B98 local stability criterion (see Section 3.2). The right-hand panel of Fig. 3, on the other hand, shows the corresponding evolution for a fixed α = −0.5, when vAϕ0 is varied. We find that as vAϕ0 increases, the growth rates of the most unstable modes increase. We also see significant deviation from the local solution for vAϕ0 = 1.5 (β0 = 0.89) at low |kz| values. This is possibly a consequence of the stronger toroidal field, which makes the radial curvature effects more important than the corresponding weaker field cases and, hence, necessitates a global analysis. Figure 1. View largeDownload slide Global eigenvalue solutions for the fiducial m = 0 case of vAϕ-PL−vAz-Cons, with α = −0.3, vAϕ0 = 0.926, and vAz = vAz0 = 0.005. Left-hand panel: the imaginary part ωI of the modal frequency or, the growth rate, as a function of kz. Right-hand panel: the real part ωR of the modal frequency as a function of kz. The global eigenvalue solutions are shown in black. The cyan (unstable modes with ωI > 0 and stable modes with ωI < 0 in the left-hand panel, and the corresponding ωR in the right-hand panel) and magenta (stable modes with ωI = 0 in the left-hand panel and the corresponding ωR in the right-hand panel) lines represent the solutions of the local dispersion relation, equation (11) with |$l^2/k_z^2=0$|⁠. The letters F, S, and A denote the fast, slow, and Alfvén modes, respectively. The global problem is solved on a radial grid r ∈ [1, 5] with resolution Nr = 256. Figure 1. View largeDownload slide Global eigenvalue solutions for the fiducial m = 0 case of vAϕ-PL−vAz-Cons, with α = −0.3, vAϕ0 = 0.926, and vAz = vAz0 = 0.005. Left-hand panel: the imaginary part ωI of the modal frequency or, the growth rate, as a function of kz. Right-hand panel: the real part ωR of the modal frequency as a function of kz. The global eigenvalue solutions are shown in black. The cyan (unstable modes with ωI > 0 and stable modes with ωI < 0 in the left-hand panel, and the corresponding ωR in the right-hand panel) and magenta (stable modes with ωI = 0 in the left-hand panel and the corresponding ωR in the right-hand panel) lines represent the solutions of the local dispersion relation, equation (11) with |$l^2/k_z^2=0$|⁠. The letters F, S, and A denote the fast, slow, and Alfvén modes, respectively. The global problem is solved on a radial grid r ∈ [1, 5] with resolution Nr = 256. Figure 2. View largeDownload slide Normalized radial eigenfunction v1r for the most unstable modes as a function of radius r, for the case discussed in Fig. 1. Top panels: kz > 0 modes. Bottom panels: kz < 0 counterparts of the respective top panels (see coloured legends for kz values; note that |kz| increases from right to left in all panels). Δr = r2 − r1 is the full-width at half maximum and Rmax is the radial peak location of the eigenfunctions. The radial axes in the right panels have been zoomed close to the inner boundary for clarity. Figure 2. View largeDownload slide Normalized radial eigenfunction v1r for the most unstable modes as a function of radius r, for the case discussed in Fig. 1. Top panels: kz > 0 modes. Bottom panels: kz < 0 counterparts of the respective top panels (see coloured legends for kz values; note that |kz| increases from right to left in all panels). Δr = r2 − r1 is the full-width at half maximum and Rmax is the radial peak location of the eigenfunctions. The radial axes in the right panels have been zoomed close to the inner boundary for clarity. Figure 3. View largeDownload slide Growth rate ωI as a function of kz for the m = 0 case of vAϕ-PL−vAz-Cons with vAz0 = 0.005. Left-hand panel: cases when vAϕ0 = 0.6 is fixed and α is varied (see coloured legends; note that α decreases from top to bottom). Right-hand panel: cases when α = −0.5 is fixed and vAϕ0 is varied (see coloured legends; note that vAϕ0 decreases from top to bottom). The dotted and dashed lines represent the global eigenvalue and local (equation 11) solutions, respectively. Figure 3. View largeDownload slide Growth rate ωI as a function of kz for the m = 0 case of vAϕ-PL−vAz-Cons with vAz0 = 0.005. Left-hand panel: cases when vAϕ0 = 0.6 is fixed and α is varied (see coloured legends; note that α decreases from top to bottom). Right-hand panel: cases when α = −0.5 is fixed and vAϕ0 is varied (see coloured legends; note that vAϕ0 decreases from top to bottom). The dotted and dashed lines represent the global eigenvalue and local (equation 11) solutions, respectively. Table 4. Classification of the internal instabilities based on the properties of their radial eigenfunctions. Full-width at half maximum Radial peak location Type of instability Δr ≲ 0.25 Rmax → 1 Local Δr > 0.25 Rmax↛ 1 Global I Δr > 0.25 Rmax → 1 Global II Δr ≲ 0.25 Rmax↛ 1 Global III Full-width at half maximum Radial peak location Type of instability Δr ≲ 0.25 Rmax → 1 Local Δr > 0.25 Rmax↛ 1 Global I Δr > 0.25 Rmax → 1 Global II Δr ≲ 0.25 Rmax↛ 1 Global III View Large Table 4. Classification of the internal instabilities based on the properties of their radial eigenfunctions. Full-width at half maximum Radial peak location Type of instability Δr ≲ 0.25 Rmax → 1 Local Δr > 0.25 Rmax↛ 1 Global I Δr > 0.25 Rmax → 1 Global II Δr ≲ 0.25 Rmax↛ 1 Global III Full-width at half maximum Radial peak location Type of instability Δr ≲ 0.25 Rmax → 1 Local Δr > 0.25 Rmax↛ 1 Global I Δr > 0.25 Rmax → 1 Global II Δr ≲ 0.25 Rmax↛ 1 Global III View Large 5.1.2 Non-axisymmetric m = 1 modes The main findings of this subsection, described below, are summarized by Figs 4–7: Fig. 4 is similar to Fig. 1, except that it represents the fiducial α = −0.3 case for the m = 1 modes. All the symbols and colours have the same meaning in both figures. We note that the m = 1 modes are no longer symmetric about the kz-axis, as shown by both the local and the global solutions. Also, non-axisymmetry seems to induce a somewhat greater deviation from the local prediction as compared to the axisymmetric case.  First, we discuss the growth rates of the modes with kz > 0. From the left-hand panel of Fig. 4 we see that these modes have a much smaller growth rate than the local prediction. The local solution predicts a non-zero growth rate ωI ≈ 0.28 at kz = 0 and instability in the range 0 ≤ kz|unst ≲ 34 . The growth rate in the global solution, on the other hand, has ωg|max ≈ 0.1 at kz|max ≈ 9 and instability for 1 ≲ kz|unst ≲ 20. Next, we discuss the growth rates of the modes with kz < 0. Note that, as already mentioned at the beginning of Section 3.1, the sign of m determines the direction of the helical twist. However, if m is assumed to be positive, as we do in this work, then a negative kz in conjunction with a positive m can be interpreted as the helical perturbation, or kink, twisting the magnetic field in a sense opposite to itself. On the other hand, a positive kz with a positive m twists the magnetic field in the same direction as the kink. This is because the local dispersion relation (equation 10) is invariant under the transformation (−m, kz) → (m, −kz) or (−m,η) → (m, −η) .  The local solution predicts two peaks with identical maximum growth rates ωg|max ≈ 0.35 (see cyan lines in the left-hand panel of Fig. 4). The first instability peak occurs at kz|max ≈ −55 and the second at kz|max ≈ −315, with the instability in the range −404 ≲ kz|unst ≤ 0. The double-peaked structure is due to the local magnetic resonance occurring at r = rres = 1 for kz|res = −mvAϕ0/vAz0 ≈ −185 (see equation 20), which segregates the two kz < 0 instability peaks. One can thus interpret the instabilities to be originating from rres = 1 (see Section 3.3) in the local scenario. However, an important point to note is that magnetic resonance is an intrinsically global condition. Thus, in order to truly understand if it plays any role, we need to check if the unstable eigenfunctions have significant amplitude at the corresponding resonant radius. This can be verified only in a global framework, as demonstrated in point (vi) below.  The global solution also exhibits two peaks in the growth rate, which are, however, not identical unlike the local solution (left-hand panel of Fig. 4). The peak values and kz locations differ slightly from the local values, matching more closely for the second peak: ωg|max ≈ 0.32 at kz|max ≈ −63 and ωg|max ≈ 0.34 at kz|max ≈ −304 and instability in the range −392 ≲ kz|unst ≤ 0. The behaviour at kz ≈ −185, where local magnetic resonance is satisfied, is also very different. In the global solution, although the first instability peak stabilizes at kz ≈ −185, the second peak does not originate from kz ≈ −185, but from kz = 0. This leads to an interaction between the unstable modes in the global scenario, which appears as a band between the two instability peaks (these are slow-mode–slow-mode interactions). On the right-hand panel of Fig. 4, we find that the most unstable modes have zero phase velocity (ωR = 0) like the m = 0 case, and are likely to be destabilized slow modes. For kz > 0, the local solutions seem to be tracing the slow and Alfven modes of the global solutions fairly well. For kz < 0, we see the appearance of possibly another branch of slow modes (marked as S2) that does not have a local counterpart, which also seems to be interacting with the Alfvén modes. For a detailed understanding of the nature of the m = 1 instabilities, we need to look at the corresponding eigenfunctions. Fig. 5 shows the eigenfunction v1r as a function of radius for the range of unstable kz corresponding to Fig. 4. The eigenfunction shown here is that of the most unstable mode at a given kz. Starting from the left, the first panel shows the eigenfunctions for 3 ≤ kz ≤ 18, the second for −135 ≤ kz ≤ −10, the third for −250 ≤ kz ≤ −140, and the last for −390 ≤ kz ≤ −270, plotted on radial axes, r ∈ [1, 5], r ∈ [1, 3], r ∈ [1, 3], and r ∈ [1, 1.2], respectively, for clarity. We now classify the instabilities according to Table 4.  The first panel of Fig. 5 shows that the eigenfunctions for kz > 0 represent Global I (kz < 5) and Global II (5 < kz < 18) modes. This explains the deviation of the global solution from the local prediction in this kz-range, as seen in the left-hand panel of Fig. 4. The decrease of vAϕ with r suppresses the growth rate of the global kz > 0 modes compared to the local value.  The modes with kz < 0 exhibit mixed behaviour. The second panel shows that the modes with kz ≥ −45 are Global II instabilities, which progress to becoming Local instabilities for −135 ≤ kz ≤ −63. This is reflected in the left-hand panel of Fig. 4, where we see that the global eigenvalues approach the local prediction (cyan line) as kz → −135.  The eigenfunctions in the third panel represent Global III instabilities and are unique in the sense that they are shifted away from both the radial boundaries. These eigenfunctions represent the interaction between the two sets of kz < 0 instabilities mentioned in point (ii) above (also see the left-hand panel of Fig. 4). As kz → −250, the eigenfunctions tend to get localized close to the inner boundary again.  Finally, the eigenfunctions in the last panel represent Local instabilities, and hence the global eigenvalue solution for this range of kz matches well with the local prediction, as seen in the left-hand panel of Fig. 4. Overall, we find that the eigenfunctions become more localized towards r = 1 as |kz| increases, similar to the m = 0 case. Fig. 6 shows the effect of varying α on the growth rates of the m = 1 modes, for a fixed vAϕ0 = 0.6 and vAz0 = 0.005 (compare with the left panel of Fig. 3 for the equivalent m = 0 case). As α decreases, the growth rates of all the modes decrease for both the local and global solutions. The kz > 0 unstable modes stabilize first. In the global eigenvalue solution these modes stabilize at α = −0.4, slightly higher than the local prediction of stabilization at α = −0.5 (see equation 15). This is probably due to the fact that the kz > 0 modes are Global instabilities, as discussed above and, hence, they stabilize earlier due to the rapid fall in vAϕ with radius as α decreases. Nevertheless, we may still consider α = −0.5 as an absolute stability criterion for the kz> 0 modes, when m = 1. As α decreases beyond −0.5, instability exists only for kz < 0. We see that the onset of instabilities as per the local solution keeps shifting to lower kz. The instabilities in the global solution also seem to pull away from kz = 0 slightly as α → −1, however at a much slower rate than the local solutions. The small |kz| modes thus correspond to Global instabilities [I or II; see point (iv) above], which deviate from the local prediction. Finally, at α = −1 there is complete stabilization of all the m = 1 modes, again validating B98’s local stability criterion. We now discuss the resonant character of the m = 1 modes from a global perspective. For the power-law profile in this section, the m = 1 modes satisfying the magnetic resonance condition (equation 20) are given by kz = −vAϕ0rα/vAz0. Thus, for our fiducial domain 1 ≤ r ≤ 5, the most unstable modes with −185 ≲ kz ≲ −114 are supposed to be resonant and the ones outside this range are non-resonant. We now look at the eigenfunctions of these modes.  Note that the most unstable modes in the range −185 ≲ kz ≲ −140, shown in the third panel of Fig. 5, are the interacting slow modes mentioned in point (ii) above, and are unlikely to be resonant. The most unstable modes in the range −135 ≲ kz ≲ −114, shown in the second panel, could in principle be resonant. In order to verify this, we need to compare the peak locations of these eigenfunctions with the corresponding resonant radius rres from equation (20).  Consider the example of kz = −135 in the second panel of Fig. 5. It is highly localized close to the inner boundary, peaking at r ∼ 1.02, i.e. well inside the resonant surface for this mode at rres = (− kzvAz0/vAϕ0)1/α ∼ 2.86 (green dashed vertical line in the second panel of Fig. 5). Thus, the amplitude of the eigenfunction is effectively zero at rres, which implies that magnetic resonance is not actually triggering these modes. Note that the role of magnetic resonance might be sensitive to the background as well as boundary conditions. For example, Bodo et al. (2013) indeed found unstable modes that peak at the corresponding resonant radius in their stability analyses carried out for relativistic, force-free, pressureless jets with surface velocity shear. Alternatively, Kim et al. (2015) showed that the presence of a finite pressure dilutes the importance of the magnetic resonance condition in triggering instabilities. With reference to the above discussion on magnetic resonance, we can further conclude that Suydam’s criterion is not applicable to check for instability in this case, as the unstable modes are not localized about rres (see Section 3.4). We do, however, find an interesting result for α = −1. The left-hand side of Suydam’s stability criterion (inequality 21) reduces, in our units and notation, to \begin{eqnarray*} \frac{1}{4r}(1-\alpha)^2 - \frac{2}{r}(1+\alpha)\frac{v_{A\phi }^2}{v_{Az}^2} , \end{eqnarray*} (37) which is ≥0 for α ≤ −1. This happens to be consistent with the stability condition for the m = 1 modes established above (i.e. B98 stability criterion). The m = 1 instabilities are likely to be responsible for the dissipation of magnetic energy into thermal energy in the jet due to the formation of current sheets (as mentioned in the Introduction). This may lead to an enhanced background pressure, which in turn can affect the instabilities. In order to study this, we add a constant background pressure (which in our case corresponds to a constant density) to our fiducial model. We test two cases: ρcons = ρ(r → ∞) = 10 and ρcons = ρ(r → ∞) = 100 (see equation 30 and Table 3). Note that the above densities are normalized using the same inner density ρin as our fiducial model, i.e. when ρcons = ρ(r → ∞) = 0, in order to maintain the same magnetic field profile for all three cases. We find that an enhancement in the background pressure leads to a suppression of the maximum growth rates, which can be roughly approximated as |$\propto \omega _{\rm fidu}\rho _{\rm cons}^{-0.5}$|⁠, where ωfidu is the maximum growth rate of the fiducial case. However, the stability criterion remains unchanged, and we observe the same range of unstable kz, as well as the same |$k_{z\rm max}$|⁠, for all the cases. We also verified that even in the presence of an enhanced background pressure, the system gets stabilized as α → −1. So far, we have restricted our discussion to the limit vAz ≪ vAϕ (see Introduction and Section 2.1 for the motivation behind this) and, hence, to predominantly pressure-driven instabilities. However, if the vertical field becomes strong enough, it can have a stabilizing effect by counteracting the hoop stress due to the toroidal field. Also, as mentioned towards the beginning of the Introduction, a stronger vertical field is likely to introduce current-driven effects, which may induce a change in the behaviour of the instabilities.  Fig. 7 shows the growth rates of the m = 1 modes as vAz0 is increased with respect to vAϕ0 for our fiducial α = −0.3 case (only global eigenvalue solutions are shown). For 0.001vAϕ0 ≤ vAz0 ≤ 0.01vAϕ0, the maximum growth rate of the m = 1 modes remains fairly similar and this is the limit up to which the local calculations are valid. As vAz0 increases further, we see that the maximum growth rate starts decreasing, although it is still appreciable for vAz0 = 0.1vAϕ0. There is a sharp drop in the maximum growth rate for vAz0 ≥ 0.5vAϕ0 and complete stabilization for vAz0 = vAϕ0, where the vertical field dominates over the toroidal field everywhere in the domain (since vAz0 is a constant and vAϕ0 decreases with radius). Nevertheless, it has been pointed out that even a strong vertical field may lead to instability if m is sufficiently high (see Bonanno & Urpin 2011a for a discussion). Note that when vAz0 = 0.1vAϕ0, the kz > 0 modes stabilize for α = −0.3 in contrast to the case shown in Fig. 6 (vAz0 = 0.005vAϕ0), where they stabilize for α = −0.4. Thus, with an increase in the vertical field, the kz > 0 modes are marginally stabilized at higher α. Similarly, there is complete stabilization for all kz modes at a higher value of α, and α ≤ −1 reduces to a sufficient (and no longer necessary) stability criterion. Figure 4. View largeDownload slide Global eigenvalue solutions for the fiducial m = 1 case of vAϕ-PL−vAz-Cons, with everything else the same as Fig. 1. In the right-hand panel, S2 denotes a new branch of slow modes that appears in the global eigenvalue solutions. Figure 4. View largeDownload slide Global eigenvalue solutions for the fiducial m = 1 case of vAϕ-PL−vAz-Cons, with everything else the same as Fig. 1. In the right-hand panel, S2 denotes a new branch of slow modes that appears in the global eigenvalue solutions. Figure 5. View largeDownload slide Normalized radial eigenfunction v1r for the most unstable modes as a function of radius r, for the case discussed in Fig. 4 (see coloured legends for kz values; note that |kz| increases from right to left in all panels). The green dashed vertical line, in the second panel from the left, indicates the resonant radius rres = 2.86 for kz = −135 (green dashed curve), given by equation (20). Note that the radial range shown differs from panel to panel for clarity, although the global problem is solved on a grid r ∈ [1, 5]. Figure 5. View largeDownload slide Normalized radial eigenfunction v1r for the most unstable modes as a function of radius r, for the case discussed in Fig. 4 (see coloured legends for kz values; note that |kz| increases from right to left in all panels). The green dashed vertical line, in the second panel from the left, indicates the resonant radius rres = 2.86 for kz = −135 (green dashed curve), given by equation (20). Note that the radial range shown differs from panel to panel for clarity, although the global problem is solved on a grid r ∈ [1, 5]. Figure 6. View largeDownload slide Growth rate ωI as a function of kz for the m = 1 case of vAϕ-PL−vAz-Cons, when vAϕ0 = 0.6 and vAz0 = 0.005 are fixed, and α is varied. The magenta vertical line indicates the kz = 0 line in each panel. The global eigenvalue solutions are shown in black. The cyan lines are the corresponding local solutions from the dispersion relation given by equation (11) with |$l^2/k_z^2=0$|⁠. Figure 6. View largeDownload slide Growth rate ωI as a function of kz for the m = 1 case of vAϕ-PL−vAz-Cons, when vAϕ0 = 0.6 and vAz0 = 0.005 are fixed, and α is varied. The magenta vertical line indicates the kz = 0 line in each panel. The global eigenvalue solutions are shown in black. The cyan lines are the corresponding local solutions from the dispersion relation given by equation (11) with |$l^2/k_z^2=0$|⁠. Figure 7. View largeDownload slide Growth rate ωI as a function of kz for the m = 1 case of vAϕ-PL−vAz-Cons, when vAϕ0 = 0.926 and α = −0.3 are fixed, and vAz0 is varied (see coloured legends; note that the maximum growth rate decreases as vAz0 increases). Only global eigenvalue solutions are shown in this figure. The kz-axis is in log scale. Figure 7. View largeDownload slide Growth rate ωI as a function of kz for the m = 1 case of vAϕ-PL−vAz-Cons, when vAϕ0 = 0.926 and α = −0.3 are fixed, and vAz0 is varied (see coloured legends; note that the maximum growth rate decreases as vAz0 increases). Only global eigenvalue solutions are shown in this figure. The kz-axis is in log scale. 5.1.3 Non-axisymmetric m = 2 modes The main findings of this subsection, discussed below, are summarized by Figs 8 and 9: Fig. 8 is similar to Figs 1 and 4, with the same colour scheme and symbolization, except that it represents the m = 2 modes. Like the m = 1 modes, the m = 2 modes are not symmetric about the kz-axis. There is also significant deviation from the local prediction. First, we discuss the growth rates of the m = 2 modes in the left-hand panel of Fig. 8. We see that there is no instability for kz > 0, as shown by both the local and the global solutions. This is true for a wide range of α, namely α ≤ 1, as seen from equation (15). The local solution predicts two peaks with the same maximum growth rate ωg|max ≈ 0.35. From the cyan lines we see that the first instability peak occurs at kz|max ≈ −240 and the second at kz|max ≈ −500, with the instability occupying a total range −590 ≲ kz|unst ≤ −151. The local magnetic resonance corresponds to kz|res = −mvAϕ0/vAz0 ≈ −370 at r = 1, which again separates the two peaks in the local solution.  The global solution also exhibits two peaks in the growth rate, which match fairly well with the local prediction (the first peak has a slightly smaller ωI ≈ 0.34 than the second peak). However, the instability onsets at a much lower |kz| in the global solution, yielding a total range −580 ≲ kz|unst ≤ −20 . Also, both the instability peaks originate from kz ≈ −20, contrary to the local solution.  Examining the right-hand panel of Fig. 8, the global solutions have similar characteristics to the m = 1 case (see right-hand panel of Fig. 4). The most unstable modes have zero phase velocity like both the m = 0 and m = 1 cases, and are likely to be destabilized slow modes. There is significant deviation from the local solutions, with a new branch of slow modes (marked S2) appearing and interacting with the Alfvén modes, just like in the m = 1 case. Fig. 9 shows the eigenfunctions v1r for a range of unstable kz. Note that we only plot the eigenfunction corresponding to the most unstable mode at a given kz to avoid confusion. Starting from the left, the first panel shows the eigenfunctions for −160 ≤ kz ≤ −25, the second for −305 ≤ kz ≤ −190, the third for −450 ≤ kz ≤ −310, and the last for −570 ≤ kz ≤ −487, with different radial axes for clarity, namely r ∈ [1, 5], r ∈ [1, 2], r ∈ [1, 1.8], and r ∈ [1, 1.1], respectively. We now classify the instabilities according to Table 4.  From the first panel of Fig. 9, we see that the eigenfunction for kz = −25 peaks towards the outer boundary. We have verified that all modes in the range −25 < kz < 0 exhibit a similar behaviour and, hence, represent Global I instabilities. These then transform to Global III instabilities in the range −160 < kz < −25, which is an indication of slow-mode–slow-mode interaction. The eigenfunctions in the second panel clearly represent Local instabilities. The eigenfunctions in the third panel represent Global III instabilities, signifying the slow-mode–slow-mode interaction in the range −450 < kz < −310. The eigenfunctions in the last panel again represent Local instabilities. Note that there are two sets of interacting instabilities for m = 2 (first and third panel of Fig. 9), unlike the m = 1 case, which has only one (third panel of Fig. 5). The m = 2 modes satisfying the magnetic resonance condition (equation 20) are given by kz = −2vAϕ0rα/vAz0. Thus, for our fiducial domain 1 ≤ r ≤ 5, the modes with −370 ≲ kz ≲ −228 are supposed to be resonant. The most unstable modes in the range −370 ≲ kz ≲ −310 are interacting in nature and, hence, non-resonant (third panel of Fig. 9). The modes in the range −305 ≲ kz ≲ −228, in the second panel, are supposed to be resonant, but they are localized far away from the corresponding resonant radii. Thus, again magnetic resonance does not play any role in triggering the m = 2 instabilities. The inapplicability of Suydam’s criterion remains the same as for the m = 1 case, since inequality (21) holds true for all m ≠ 0. Figure 8. View largeDownload slide Global eigenvalue solutions for the fiducial m = 2 case of vAϕ-PL−vAz-Cons, with everything else the same as Figs 1 and 4. Figure 8. View largeDownload slide Global eigenvalue solutions for the fiducial m = 2 case of vAϕ-PL−vAz-Cons, with everything else the same as Figs 1 and 4. Figure 9. View largeDownload slide Normalized radial eigenfunction v1r for the most unstable modes as a function of radius r, for the case discussed in Fig. 8 (see coloured legends for kz values; note that |kz| increases from right to left in all panels). The red dashed vertical line, in the second panel from the left, indicates the resonant radius rres = 1.91 for kz = −305 (red dashed curve), given by equation (20). Note that the radial axis has been plotted for different extents in the four panels for clarity, although the global problem is solved on a grid r ∈ [1, 5]. Figure 9. View largeDownload slide Normalized radial eigenfunction v1r for the most unstable modes as a function of radius r, for the case discussed in Fig. 8 (see coloured legends for kz values; note that |kz| increases from right to left in all panels). The red dashed vertical line, in the second panel from the left, indicates the resonant radius rres = 1.91 for kz = −305 (red dashed curve), given by equation (20). Note that the radial axis has been plotted for different extents in the four panels for clarity, although the global problem is solved on a grid r ∈ [1, 5]. The global solutions yield very similar maximum growth rates for the m = 0, 1, and 2 cases studied in this section (see left-hand panels of Figs 1, 4, and 8). This is in accordance with the local prediction for pressure-driven instabilities (see equation 17). Although not shown here, we have verified that the m = 2 modes also get completely stabilized as α → −1. Physically this condition means that the destabilizing thermal pressure gradient vanishes when α = −1 (see equation 9). Interestingly, the background equilibrium becomes force free but not pressureless when α = −1. Thus, α ≤ −1 can be interpreted as a global stability criterion for all axisymmetric and non-axisymmetric modes, for a nearly constant vertical field such that vAz ≪ vAϕ. The reason the local B98 stability criterion holds true even for the global solutions is that the fastest-growing modes are purely Local instabilities (occurring at large |kz|) irrespective of m. 5.2 Results for a power-law toroidal field and varying vertical field: |$\boldsymbol{ v}$|Aϕ-PL−|$\boldsymbol{ v}$|Az-Var In this section, we discuss the effect of a radially varying (Var) vertical field on the instabilities (for a power-law toroidal field). In order to do so, we make the following assumption about the relation between the thermal pressure gradient and the vertical field gradient: \begin{eqnarray*} \frac{\mathrm{ d}}{\mathrm{ d}r} \biggl (\frac{B_z^2}{8\pi }\biggr) = \biggl (\frac{1-\epsilon }{\epsilon } \biggr) \frac{\mathrm{ d}P}{\mathrm{ d}r} , \end{eqnarray*} (38) where 0 < ϵ < 1. Thus, a larger value of ϵ indicates a stronger thermal pressure gradient and vice versa. Note that equation (38) is an idealization and may not hold true inside a jet column. Nevertheless, this formalism helps us understand the basic contribution of a vertical field gradient towards the stability of the system. Using the above equation and assuming vAϕ = vAϕ0rα, with α < 0, we can integrate equation (7) to obtain the background density (or equivalently, thermal pressure) as a function of radius. We impose the boundary conditions (31) and (34) to obtain \begin{eqnarray*} \rho (r) = r^{2\alpha } \quad {\rm and} \quad v_{Az}(r) = v_{Az0}r^{\alpha } , \end{eqnarray*} (39) such that \begin{eqnarray*} v_{Az0} = \sqrt{\frac{2(1-\epsilon)}{\epsilon }}\quad {\rm and}\quad v_{A\phi 0} = \sqrt{\frac{-2\alpha }{(1+\alpha)\epsilon }} . \end{eqnarray*} (40) Note that we also assumed vAz → 0 as r → ∞. The plasma-beta for this case is a constant with respect to radius (for a given α and ϵ) such that \begin{eqnarray*} \beta = \frac{2}{(v_{A\phi 0}^2 + v_{Az0} ^2)} = \frac{\epsilon (1+\alpha)}{1 - \epsilon (1+\alpha)} . \end{eqnarray*} (41) The magnetic pitch for this case can be written using equation (6) as \begin{eqnarray*} {\cal P} = r\frac{v_{Az0}}{v_{A\phi 0}} = r \sqrt{\frac{(1-\epsilon)(1+\alpha)}{-\alpha }} , \end{eqnarray*} (42) which always increases with radius. We shall now study the effect of varying ϵ on the stability of the m = 0 and m = 1 modes for α = −0.3. The cases are listed in Table 5. We use r ∈ [1, 5] for all cases. Table 5. Summary of runs for the case vAϕ-PL−vAz-Var with α = −0.3 presented in Section 5.2. m ϵ |$v_{A\phi 0} = \sqrt{\frac{-2\alpha }{\epsilon (1+\alpha)}}$| |$v_{Az0} = \sqrt{\frac{2(1-\epsilon)}{\epsilon }}$| vAz0/vAϕ0 Nr B.C.s imposed on ρ(r) 0 0.999 0.926 0.045 0.048 200 ρ(1) = 1; ρ(r → ∞) = 0 0 0.99 0.930 0.142 0.153 200 ρ(1) = 1; ρ(r → ∞) = 0 0 0.9 0.976 0.471 0.483 200 ρ(1) = 1; ρ(r → ∞) = 0 0 0.8 1.035 0.707 0.683 200 ρ(1) = 1; ρ(r → ∞) = 0 1 0.999 0.926 0.045 0.048 256 ρ(1) = 1; ρ(r → ∞) = 0 1 0.99 0.930 0.142 0.153 256 ρ(1) = 1; ρ(r → ∞) = 0 1 0.9 0.976 0.471 0.483 256 ρ(1) = 1; ρ(r → ∞) = 0 1 0.8 1.035 0.707 0.683 256 ρ(1) = 1; ρ(r → ∞) = 0 1 0.7 1.106 0.926 0.837 200 ρ(1) = 1; ρ(r → ∞) = 0 1 0.4 1.464 1.732 1.18 300 ρ(1) = 1; ρ(r → ∞) = 0 m ϵ |$v_{A\phi 0} = \sqrt{\frac{-2\alpha }{\epsilon (1+\alpha)}}$| |$v_{Az0} = \sqrt{\frac{2(1-\epsilon)}{\epsilon }}$| vAz0/vAϕ0 Nr B.C.s imposed on ρ(r) 0 0.999 0.926 0.045 0.048 200 ρ(1) = 1; ρ(r → ∞) = 0 0 0.99 0.930 0.142 0.153 200 ρ(1) = 1; ρ(r → ∞) = 0 0 0.9 0.976 0.471 0.483 200 ρ(1) = 1; ρ(r → ∞) = 0 0 0.8 1.035 0.707 0.683 200 ρ(1) = 1; ρ(r → ∞) = 0 1 0.999 0.926 0.045 0.048 256 ρ(1) = 1; ρ(r → ∞) = 0 1 0.99 0.930 0.142 0.153 256 ρ(1) = 1; ρ(r → ∞) = 0 1 0.9 0.976 0.471 0.483 256 ρ(1) = 1; ρ(r → ∞) = 0 1 0.8 1.035 0.707 0.683 256 ρ(1) = 1; ρ(r → ∞) = 0 1 0.7 1.106 0.926 0.837 200 ρ(1) = 1; ρ(r → ∞) = 0 1 0.4 1.464 1.732 1.18 300 ρ(1) = 1; ρ(r → ∞) = 0 View Large Table 5. Summary of runs for the case vAϕ-PL−vAz-Var with α = −0.3 presented in Section 5.2. m ϵ |$v_{A\phi 0} = \sqrt{\frac{-2\alpha }{\epsilon (1+\alpha)}}$| |$v_{Az0} = \sqrt{\frac{2(1-\epsilon)}{\epsilon }}$| vAz0/vAϕ0 Nr B.C.s imposed on ρ(r) 0 0.999 0.926 0.045 0.048 200 ρ(1) = 1; ρ(r → ∞) = 0 0 0.99 0.930 0.142 0.153 200 ρ(1) = 1; ρ(r → ∞) = 0 0 0.9 0.976 0.471 0.483 200 ρ(1) = 1; ρ(r → ∞) = 0 0 0.8 1.035 0.707 0.683 200 ρ(1) = 1; ρ(r → ∞) = 0 1 0.999 0.926 0.045 0.048 256 ρ(1) = 1; ρ(r → ∞) = 0 1 0.99 0.930 0.142 0.153 256 ρ(1) = 1; ρ(r → ∞) = 0 1 0.9 0.976 0.471 0.483 256 ρ(1) = 1; ρ(r → ∞) = 0 1 0.8 1.035 0.707 0.683 256 ρ(1) = 1; ρ(r → ∞) = 0 1 0.7 1.106 0.926 0.837 200 ρ(1) = 1; ρ(r → ∞) = 0 1 0.4 1.464 1.732 1.18 300 ρ(1) = 1; ρ(r → ∞) = 0 m ϵ |$v_{A\phi 0} = \sqrt{\frac{-2\alpha }{\epsilon (1+\alpha)}}$| |$v_{Az0} = \sqrt{\frac{2(1-\epsilon)}{\epsilon }}$| vAz0/vAϕ0 Nr B.C.s imposed on ρ(r) 0 0.999 0.926 0.045 0.048 200 ρ(1) = 1; ρ(r → ∞) = 0 0 0.99 0.930 0.142 0.153 200 ρ(1) = 1; ρ(r → ∞) = 0 0 0.9 0.976 0.471 0.483 200 ρ(1) = 1; ρ(r → ∞) = 0 0 0.8 1.035 0.707 0.683 200 ρ(1) = 1; ρ(r → ∞) = 0 1 0.999 0.926 0.045 0.048 256 ρ(1) = 1; ρ(r → ∞) = 0 1 0.99 0.930 0.142 0.153 256 ρ(1) = 1; ρ(r → ∞) = 0 1 0.9 0.976 0.471 0.483 256 ρ(1) = 1; ρ(r → ∞) = 0 1 0.8 1.035 0.707 0.683 256 ρ(1) = 1; ρ(r → ∞) = 0 1 0.7 1.106 0.926 0.837 200 ρ(1) = 1; ρ(r → ∞) = 0 1 0.4 1.464 1.732 1.18 300 ρ(1) = 1; ρ(r → ∞) = 0 View Large 5.2.1 Effect on m = 0 and m = 1 modes The left-hand panels of Fig. 10 show the profiles for the toroidal (solid) and vertical (dashed) fields as ϵ is varied, keeping α = −0.3 fixed (top left: ϵ = 0.999, 0.99, 0.9; bottom left: ϵ = 0.8, 0.7, 0.4). The background density (and pressure) is independent of ϵεand is shown by the dot–dashed cyan line in the left-hand panels (see equation 39). We see from the figure that for ϵ = 0.999, vAz is very sub-dominant compared to vAϕ. Note that although vAϕ for ϵ = 0.999 and ϵ = 0.99 coincides (black and red solid lines in the top left-hand panel of Fig. 10), vAz is different (black and red dashed lines in the same figure). As ϵ decreases, vAz approaches vAϕ quite steeply, and finally dominates over vAϕ when ϵ = 0.4. Note that since vAz follows the same power law as vAϕ, the two profiles are identical when vAϕ0 = vAz0. This happens when ϵ = (1 + 2α)/(1 + α), giving ϵ = 0.57 for α = −0.3 and β = 0.667. Thus, for ϵ < 0.57, the vertical field is stronger than the toroidal field. Figure 10. View largeDownload slide Global eigenvalue solutions for the case vAϕ-PL−vAz-Var with α = −0.3. Left-hand panels: background radial profiles of vAϕ (solid lines), vAz (dashed lines), and ρ (cyan dot–dashed line) for different ϵ given by equation (38) (see coloured legends; note that ϵ increases from top to bottom for the dashed curves). Right-hand panels: growth rate ωI as a function of kz for the m = 0 (top) and m = 1 (bottom) cases, corresponding to the background profiles in the left-hand panels. Note that the maximum growth rate as well as the maximum unstable |kz| decreases as ϵ decreases. Figure 10. View largeDownload slide Global eigenvalue solutions for the case vAϕ-PL−vAz-Var with α = −0.3. Left-hand panels: background radial profiles of vAϕ (solid lines), vAz (dashed lines), and ρ (cyan dot–dashed line) for different ϵ given by equation (38) (see coloured legends; note that ϵ increases from top to bottom for the dashed curves). Right-hand panels: growth rate ωI as a function of kz for the m = 0 (top) and m = 1 (bottom) cases, corresponding to the background profiles in the left-hand panels. Note that the maximum growth rate as well as the maximum unstable |kz| decreases as ϵ decreases. The top right-hand panel of Fig. 10 shows the growth rates as a function of kz (from the global eigenvalue analysis) for the m = 0 modes, when ϵ = 0.999, 0.99, 0.9, and 0.8. The bottom right-hand panel shows the same for the m = 1 modes, when ϵ = 0.999, 0.99, 0.9, 0.8, 0.7, and 0.4. We see that for both the axisymmetric and non-axisymmetric cases, as the vertical field gradient and strength become stronger compared to the thermal pressure gradient (which is independent of ϵ), the growth rates decrease drastically. The m = 0 modes appear to stabilize completely at a critical value, ϵ = ϵcr ∼ 0.8, contrary to the m = 1 modes, which seem to persist under stronger vertical fields, and stabilize completely when ϵcr ∼ 0.4. The kz > 0 modes for the m = 1 case, however, stabilize for ϵ = 0.999. This is because the vertical field strength for this case is almost a factor of 10 higher at the inner boundary, i.e. vAz0 = 0.045, compared to the fiducial vAϕ-PL−vAz-Cons case shown in Fig. 4 with vAz0 = 0.005 (although the toroidal field profile remains similar for both the cases). This behaviour is very similar to the stabilizing effect of a strong (constant) vertical field shown in Fig. 7, and is likely due to current-driven effects. Note that as the vertical field gradient, and hence strength, gets stronger, the instabilities shift to smaller |kz| and become increasingly Global in nature (and thus can no longer be captured in a local analysis). This makes determining a global stability criterion for this case more difficult, as the B98 stability criterion is not applicable in this limit (unlike in Section 5.1). Absolute stability is attained for a critical combination of the magnetic pressure gradients and the thermal pressure gradient in the jet, as indicated by ϵcr (for a given m and α). In the case of m ≠ 0 modes, this is similar to the statement of Suydam’s stability criterion, which however is not directly applicable in this case because of its local nature (see Section 3.4). Alternatively, it has been suggested by Bromberg & Tchekhovskoy (2016) that the internal kink instability (equivalent to the m = 1 modes in this work) stabilizes when there is equipartition between the thermal and magnetic pressures. However, that seems to be neither a necessary nor a sufficient condition from our analysis. For instance, equipartition is attained when (in our units) \begin{eqnarray*} \rho = \frac{1}{2}\left(v_{A\phi }^2 + v_{Az}^2\right) \end{eqnarray*} (43) or, equivalently, β = 1 in equation (41), which yields ϵ = 0.71 for α = −0.3. However, from Fig. 10 we see that the m = 0 modes attain stability for ϵ > 0.71, when the thermal pressure dominates over the total magnetic pressure, while the m = 1 modes attain stability for ϵ < 0.71, when the total magnetic pressure dominates over the thermal pressure. Although ϵ is just a parametrization, our results indicate that pressure balance alone may not guarantee stability, and one needs to consider the interplay between the magnetic and thermal pressure gradients as well (see also Section 6 below). 5.3 Results for a generic toroidal field and constant vertical field: |$\boldsymbol{ v}$|Aϕ-Gen−|$\boldsymbol{ v}$|Az-Cons In this section, we consider a more generic profile for the background toroidal field such that \begin{eqnarray*} v_{A\phi } = v_{A\phi 0}\frac{(1+a_1)r}{1+ a_1r^2} , \end{eqnarray*} (44) where a1 is a constant parameter governing the radial location of the peak of vAϕ. The above equation mimics the magnetic field profiles obtained by Begelman & Li (1992) for their Crab Nebula models (see also Appl et al. 2000; O’Neill et al. 2012). We consider a constant vertical field vAz = vAz0, and obtain the density by integrating equation (7) subject to the boundary condition (31), \begin{eqnarray*} \rho (r) = 1 + \frac{v_{A\phi 0}^2(1+a_1)^2}{2 a_1} \biggl (\frac{1}{(1 + a_1 r^2)^2} - \frac{1}{(1+a_1)^2} \biggr) . \end{eqnarray*} (45) For our fiducial model, we also impose the boundary condition (34) to obtain \begin{eqnarray*} \rho (r) = \frac{(1 +a_1)^2}{(1 + a_1 r^2)^2} \quad {\rm and} \quad v_{A\phi 0}= \sqrt{2a_1} . \end{eqnarray*} (46) The plasma-beta for the fiducial model, on assuming vAz0 ≪ vAϕ, then becomes \begin{eqnarray*} \beta \approx \frac{1}{a_1 r^2} . \end{eqnarray*} (47) The pitch for this magnetic configuration follows from equation (6): \begin{eqnarray*} {\cal P} = \frac{v_{Az0}(1+a_1r^2)}{v_{A\phi 0}(1+a_1)} , \end{eqnarray*} (48) which increases with radius. Fig. 11 shows the different background profiles as a function of r for our fiducial case having a1 = 0.1, vAϕ0 = 0.447, and vAz0 = 0.01 ≪ vAϕ0. We use r ∈ [1, 5] for all the cases studied in this section, which are listed in Table 6. Fig. 12 shows ωI (left-hand panels) and ωR (right-hand panels) as a function of kz for the m = 0 (top), m = 1 (middle), and m = 2 (bottom) modes described by the fiducial model shown in Fig. 11. The global eigenvalue solutions are in black, while the cyan (ωI ≠ 0) and magenta (ωI = 0) lines indicate local solutions derived using the B98 dispersion relation with |$l^2/k_z^2=0$| (equation 11) for a power-law toroidal field vAϕ = 0.447r0.8. Note that there is no fixed α corresponding to the generic profile; however, a slope of α = 0.8 is tangent to the generic profile at the inner radius (see Fig. 11). This gives us the opportunity to directly compare the solutions of the generic and power-law toroidal fields. Note that all colours and symbols have the same meaning as in Fig. 1. The main findings of this section, discussed below, are summarized by Figs 12, 13, and 14: The growth rates of the m = 0 modes are symmetric about the kz-axis, as expected. The global solution matches well with the local solution (representing the most unstable modes) for |kz| ≳ 30. The local power-law solution predicts a maximum growth rate ωg|max ≈ 0.42 (see equation 17), the corresponding vertical wavenumbers kz|max ≈ ±40 (see equation 19), and the range of instability −85 ≲ kz|unst ≲ 85 (see relation 15). All the corresponding global quantities are in very good agreement but with slightly lower values. However, the growth rates are somewhat higher than the local prediction for |kz| < 30. Note that unlike the m = 0 case for α = −0.3 shown in Fig. 1, the most unstable modes do not have ωI → 0 as kz → 0. This can be explained by the fact that the present solution is roughly approximated by a α > 0 power-law toroidal field, which tends to give higher growth rates because vAϕ is increasing with radius. Thus, the solutions behave differently from our fiducial α < 0 case (see Section 5.1.1). The envelope of the global solutions for both the m = 1 and m = 2 cases, i.e. the most unstable modes, seem to match reasonably well with the local prediction, with some differences. The local solutions for both m = 1 and 2 have a double-peaked structure. This is again due to the occurrence of local magnetic resonance at r = 1 for m = 1 with kz|res = −vAϕ0/vAz0 ≈ −45, and for m = 2 with kz|res = −2vAϕ0/vAz0 ≈ −89. However, ωI ≠ 0 at kz|res, unlike the behaviour of the local resonance discussed in Section 5.1.2. The two peaks in the global solution are barely distinguishable due to slow-mode–slow-mode interaction, leading to a band of higher growth rate about kz ≈ −45 for m = 1 and about kz ≈ −89 for m = 2. For m = 1, the local solution slightly overestimates the growth rates for the kz > 0 modes, as well as that of the first peak. For m = 2, the kz > 0 modes are completely stabilized; however, the kz < 0 modes start slightly further away from the axis than the local prediction. In both cases, the local solution predicts the total range of instability extremely well. From the right-hand panels of Fig. 12, we find that the most unstable modes are destabilized slow modes that are purely growing, for all m = 0, 1, and 2 cases. In fact, the local predictions match quite well with the global solutions across all kz and m, with no new branches of slow-Alfvén interaction appearing in the ωR–kz plane, unlike those observed in Section 5.1. In Fig. 13, we illustrate the normalized, v1r eigenfunctions corresponding to the most unstable m = 1 modes. The kz > 0 modes shown in the first panel (from left) are Global II stabilities, which tend to get more localized to the inner boundary as kz increases. The unstable modes in the range −20 ≤ kz ≤ −5 shown in the second panel also represent Global II instabilities, however with a larger Δr than the first panel. The third panel showing −70 ≤ kz ≤ −30 represents Global III instabilities, i.e. the interacting slow modes. Finally, the last panel showing −125 < kz < −85 represents Local instabilities.  Note that the m = 1 modes satisfying the magnetic resonance condition for the present field profile are given by kz = −(vAϕ0/vAz0)(1 + a1)/(1 + a1r2). Thus, for r ∈ [1, 5], the modes in −45 ≲ kz ≲ −14 are supposed to be resonant. Again, the most unstable modes in −45 ≲ kz ≲ −30 are interacting and non-resonant (third panel), while the eigenfunctions of the modes in −20 ≲ kz ≲ −14 peak far away from the resonant radius (second panel). For example, the eigenfunction for kz = −20 peaks at r ∼ 1.35, which is inside the resonant surface for this mode at |$r_{\rm res}= \sqrt{(1/a_1)(-1 - v_{A\phi 0}(1+a_1)/(k_z v_{Az0}))} \sim 3.82$| (vertical red dashed line in the second panel of Fig. 13). This is similar to the conclusion drawn regarding the unstable modes discussed in Section 5.1.2.  As to the eigenfunctions of the most unstable m = 0 modes, they represent Global I instabilities for |kz| < 15, Global II instabilities for 15 < |kz| < 3, and Local instabilities for |kz| > 30 – which corroborates the deviation of the global solution from the local prediction at small |kz|. The eigenfunctions of the most unstable m = 2 modes represent Global II instabilities for −40 < kz < −10, Local instabilities for −70 < kz < −40, Global III instabilities for −115 < kz < −70 (interacting modes), and again Local instabilities for −170 < kz < −115. Finally, in Fig. 14, we compare the global solutions for the m = 1 modes corresponding to multiple generic profiles having different a1. The toroidal magnetic fields (coloured solid lines) are shown as a function of radius in the top left-hand panel. They are constructed such that they all have the same vAϕ0 = 0.447 and vAz = constant = 0.01 as the fiducial model (see Table 6 for details). We note that as a1 increases, the maximum magnetic field strength decreases as well as shifting towards the inner boundary. The dotted lines are ∝ r−0.5, which become tangential to the different magnetic field curves at different radii. The aim here is to establish a correlation with the solutions of the power-law profiles discussed in Section 5.1. The bottom panel shows the growth rate as a function of kz corresponding to the cases in the top left-hand panel (global eigenvalue solutions only). We observe that as a1 increases, the growth rates of the unstable modes decrease and the range of instability shrinks. Furthermore, the two growth rate peaks (for kz < 0) become more distinguishable and the modes start resembling the left-hand panel of Fig. 4. Interestingly, the kz > 0 modes stabilize for a1 > 0.4, which we now try to explain using the top two panels (note that some of the incomplete kz > 0 curves are due to inadequate resolution). In Section 5.1.2 and Fig. 6, we saw that the kz > 0 modes are completely stabilized when α = −0.5 as per the B98 prediction. Keeping that in mind, in the top left-hand panel of Fig. 14 we observe that the radius at which the r−0.5 curve is tangential to the respective field profiles shifts radially inwards as a1 increases. Now, the most unstable kz > 0 mode for each a1 roughly corresponds to kz = 9 (see bottom panel). The top right-hand panel shows the v1r eigenfunctions for kz = 9 when a1 = 0.1, 0.2, 0.3, and 0.4. All four eigenfunctions represent Global II instabilities, which implies that they are influenced by the background toroidal field profile. More interestingly, the eigenfunctions also shift radially inwards as a1 increases. Thus, we can conclude that if the toroidal field transitions to a r−0.5 at sufficiently small radii to influence the eigenfunctions of the kz > 0 unstable modes, they get stabilized. For the present scenario, this happens when a1 = 0.5. Similarly, we can conclude that there will be complete stabilization (i.e. of all kz modes) with a further increase in a1, as the background toroidal field transitions to a r−1 profile. Thus, we recover the B98 stability criterion even for a more complex magnetic background, precisely due to the local nature of the fastest-growing internal (pressure-driven) instabilities (provided vAz ≪ vAϕ). Figure 11. View largeDownload slide Background radial profiles for the fiducial vAϕ–Gen–vAz–Cons case, with vAϕ (solid line) given by equation (44), ρ (dotted line) given by equation (46), vAϕ0 = 0.447, a1 = 0.1 , and vAz = constant = 0.01 (dashed line). The dot–dashed line shows a power-law profile for vAϕ = 0.447r0.8. Figure 11. View largeDownload slide Background radial profiles for the fiducial vAϕ–Gen–vAz–Cons case, with vAϕ (solid line) given by equation (44), ρ (dotted line) given by equation (46), vAϕ0 = 0.447, a1 = 0.1 , and vAz = constant = 0.01 (dashed line). The dot–dashed line shows a power-law profile for vAϕ = 0.447r0.8. Figure 12. View largeDownload slide Global eigenvalue solutions for the fiducial vAϕ-Gen−vAz-Cons case, with vAϕ0 = 0.447, a1 = 0.1 , and vAz = 0.01. The top, middle, and bottom panels represent the m = 0, 1, and 2 cases, respectively. Left-hand panels: growth rate ωI as a function of kz. Right-hand panels: the real part of the eigenvalue ωR as a function of kz. The cyan (unstable modes with ωI > 0 and stable modes with ωI < 0 in the left-hand panels, and the corresponding ωR in the right-hand panels) and magenta (stable modes with ωI = 0 in the left-hand panels and the corresponding ωR in the right-hand panels) lines represent the solutions of the local dispersion relation, equation (11), for vAϕ = 0.447r0.8, vAz = 0.01, and |$l^2/k_z^2=0$|⁠. The global problem is solved on a radial grid r ∈ [1, 5] with resolution Nr = 200. All symbols have the same meaning as in Fig. 1. Figure 12. View largeDownload slide Global eigenvalue solutions for the fiducial vAϕ-Gen−vAz-Cons case, with vAϕ0 = 0.447, a1 = 0.1 , and vAz = 0.01. The top, middle, and bottom panels represent the m = 0, 1, and 2 cases, respectively. Left-hand panels: growth rate ωI as a function of kz. Right-hand panels: the real part of the eigenvalue ωR as a function of kz. The cyan (unstable modes with ωI > 0 and stable modes with ωI < 0 in the left-hand panels, and the corresponding ωR in the right-hand panels) and magenta (stable modes with ωI = 0 in the left-hand panels and the corresponding ωR in the right-hand panels) lines represent the solutions of the local dispersion relation, equation (11), for vAϕ = 0.447r0.8, vAz = 0.01, and |$l^2/k_z^2=0$|⁠. The global problem is solved on a radial grid r ∈ [1, 5] with resolution Nr = 200. All symbols have the same meaning as in Fig. 1. Figure 13. View largeDownload slide Normalized radial eigenfunction v1r for the most unstable modes as a function of radius r, for the m = 1 case shown in Fig. 12 (see coloured legends for kz values; note that |kz| increases from right to left in all panels). The red dashed vertical line, in the second panel from the left, indicates the resonant radius rres = 3.82 for kz = −20 (red dashed curve), given by equation (20). Note that the radial axis in the last panel has been zoomed close to the inner boundary for clarity. Figure 13. View largeDownload slide Normalized radial eigenfunction v1r for the most unstable modes as a function of radius r, for the m = 1 case shown in Fig. 12 (see coloured legends for kz values; note that |kz| increases from right to left in all panels). The red dashed vertical line, in the second panel from the left, indicates the resonant radius rres = 3.82 for kz = −20 (red dashed curve), given by equation (20). Note that the radial axis in the last panel has been zoomed close to the inner boundary for clarity. Figure 14. View largeDownload slide Global eigenvalue solutions for the m = 1 case of vAϕ-Gen−vAz-Cons, when vAϕ0 = 0.447 and vAz = 0.01 are fixed, and a1 is varied. Top left-hand panel: background radial profiles of vAϕ given by equation (44), with different a1 (see coloured solid curves; note that a1 increases from top to bottom). The dotted lines represent the r−1/2 profiles tangential to each vAϕ. Bottom panel: growth rate as a function kz for different a1 (the colour scheme is the same as in the top left-hand panel; as a1 increases, the maximum growth rate decreases). Top right-hand panel: normalized radial eigenfunction v1r as a function of r, corresponding to the most unstable mode at kz = 9 for different a1 (note that a1 increases from right to left). Figure 14. View largeDownload slide Global eigenvalue solutions for the m = 1 case of vAϕ-Gen−vAz-Cons, when vAϕ0 = 0.447 and vAz = 0.01 are fixed, and a1 is varied. Top left-hand panel: background radial profiles of vAϕ given by equation (44), with different a1 (see coloured solid curves; note that a1 increases from top to bottom). The dotted lines represent the r−1/2 profiles tangential to each vAϕ. Bottom panel: growth rate as a function kz for different a1 (the colour scheme is the same as in the top left-hand panel; as a1 increases, the maximum growth rate decreases). Top right-hand panel: normalized radial eigenfunction v1r as a function of r, corresponding to the most unstable mode at kz = 9 for different a1 (note that a1 increases from right to left). Table 6. Summary of runs for case vAϕ-Gen−vAz-Cons presented in Section 5.3. m a1 |$v_{A\phi 0}$| vAϕ0(1 + a1) |$v_{Az0}$| Nr B.C.s imposed on ρ(r) Other notes 0,1,2 0.1 0.447 0.492 0.01 200 ρ(1) = 1; ρ(r → ∞) = 0 Fiducial model; |$v_{A\phi 0}= \sqrt{2a_1}$| 1 0.2 0.447 0.537 0.01 200 ρ(1) = 1 ρ ≥ 0 ∀ r 1 0.3 0.447 0.581 0.01 200 ρ(1) = 1 ρ ≥ 0 ∀ r 1 0.4 0.447 0.626 0.01 200 ρ(1) = 1 ρ ≥ 0 ∀ r 1 0.5 0.447 0.671 0.01 200 ρ(1) = 1 ρ ≥ 0 ∀ r m a1 |$v_{A\phi 0}$| vAϕ0(1 + a1) |$v_{Az0}$| Nr B.C.s imposed on ρ(r) Other notes 0,1,2 0.1 0.447 0.492 0.01 200 ρ(1) = 1; ρ(r → ∞) = 0 Fiducial model; |$v_{A\phi 0}= \sqrt{2a_1}$| 1 0.2 0.447 0.537 0.01 200 ρ(1) = 1 ρ ≥ 0 ∀ r 1 0.3 0.447 0.581 0.01 200 ρ(1) = 1 ρ ≥ 0 ∀ r 1 0.4 0.447 0.626 0.01 200 ρ(1) = 1 ρ ≥ 0 ∀ r 1 0.5 0.447 0.671 0.01 200 ρ(1) = 1 ρ ≥ 0 ∀ r View Large Table 6. Summary of runs for case vAϕ-Gen−vAz-Cons presented in Section 5.3. m a1 |$v_{A\phi 0}$| vAϕ0(1 + a1) |$v_{Az0}$| Nr B.C.s imposed on ρ(r) Other notes 0,1,2 0.1 0.447 0.492 0.01 200 ρ(1) = 1; ρ(r → ∞) = 0 Fiducial model; |$v_{A\phi 0}= \sqrt{2a_1}$| 1 0.2 0.447 0.537 0.01 200 ρ(1) = 1 ρ ≥ 0 ∀ r 1 0.3 0.447 0.581 0.01 200 ρ(1) = 1 ρ ≥ 0 ∀ r 1 0.4 0.447 0.626 0.01 200 ρ(1) = 1 ρ ≥ 0 ∀ r 1 0.5 0.447 0.671 0.01 200 ρ(1) = 1 ρ ≥ 0 ∀ r m a1 |$v_{A\phi 0}$| vAϕ0(1 + a1) |$v_{Az0}$| Nr B.C.s imposed on ρ(r) Other notes 0,1,2 0.1 0.447 0.492 0.01 200 ρ(1) = 1; ρ(r → ∞) = 0 Fiducial model; |$v_{A\phi 0}= \sqrt{2a_1}$| 1 0.2 0.447 0.537 0.01 200 ρ(1) = 1 ρ ≥ 0 ∀ r 1 0.3 0.447 0.581 0.01 200 ρ(1) = 1 ρ ≥ 0 ∀ r 1 0.4 0.447 0.626 0.01 200 ρ(1) = 1 ρ ≥ 0 ∀ r 1 0.5 0.447 0.671 0.01 200 ρ(1) = 1 ρ ≥ 0 ∀ r View Large 6 DISCUSSION In this section, we discuss some of our results in the context of recent numerical simulations. Bromberg & Tchekhovskoy (2016) carried out global, 3D relativistic MHD simulations of relativistic, Poynting-flux-dominated jets. They simulated two kinds of jets – (i) headless jets that propagate in a vacuum funnel and experience no direct influence from the ambient medium and (ii) headed jets that propagate through an ambient medium. They found that the ambient medium is responsible for collimating the headed jets, and that this compression effect triggers the internal kink instability (and, also, the external kink instability that we do not discuss here). Since these instabilities lead to the dissipation of magnetic energy, the simulations showed a build-up of background thermal pressure, which can no longer be ignored under the force-free approximation. One of the main findings of Bromberg & Tchekhovskoy (2016) was that the internal kink instabilities saturated once equipartition was reached between the thermal and magnetic energies, about which we commented briefly in Section 5.2 above. In Fig. 15, we reconstruct a situation similar (but not identical) to that in Fig. 14 of Bromberg & Tchekhovskoy (2016), to analyse the conditions under which the internal kink stabilizes. Figure 15. View largeDownload slide Testing the role of equipartition between thermal and magnetic pressures on the stability of m = 1 modes. Top left-hand panel: background radial profiles of a toroidal Alfvén velocity vAϕ = 1.41r−0.5 (blue solid line), constant vertical Alfvén velocity vAz = 0.14 (green dashed line), total Alfvén velocity |$v_{\mathrm{ Atot}} = \sqrt{v_{A\phi }^2 + v_{Az}^2}$| (solid grey line), and |$\sqrt{2\rho }$|⁠, which is a quantity proportional to the thermal pressure (red dot–dashed line), so defined to compare with the Alfvén velocities. Bottom panel: growth rate as a function of kz from the global eigenvalue analysis of the case shown in the top left-hand panel for m = 1, and radial grid resolutions Nr = 200 (black) and Nr = 256 (red). Top right-hand panel: normalized radial eigenfunction v1r as a function of r, corresponding to the most unstable modes at kz = −2, −4, and −12 (see coloured legends for kz values; note that |kz| increases from right to left). Figure 15. View largeDownload slide Testing the role of equipartition between thermal and magnetic pressures on the stability of m = 1 modes. Top left-hand panel: background radial profiles of a toroidal Alfvén velocity vAϕ = 1.41r−0.5 (blue solid line), constant vertical Alfvén velocity vAz = 0.14 (green dashed line), total Alfvén velocity |$v_{\mathrm{ Atot}} = \sqrt{v_{A\phi }^2 + v_{Az}^2}$| (solid grey line), and |$\sqrt{2\rho }$|⁠, which is a quantity proportional to the thermal pressure (red dot–dashed line), so defined to compare with the Alfvén velocities. Bottom panel: growth rate as a function of kz from the global eigenvalue analysis of the case shown in the top left-hand panel for m = 1, and radial grid resolutions Nr = 200 (black) and Nr = 256 (red). Top right-hand panel: normalized radial eigenfunction v1r as a function of r, corresponding to the most unstable modes at kz = −2, −4, and −12 (see coloured legends for kz values; note that |kz| increases from right to left). In their fig. 14, Bromberg & Tchekhovskoy (2016) showed the azimuthally averaged magnetic field and thermal pressure profiles inside their 3D simulated jet, at a height where there was no more dissipation due to the internal kink instability. The top left-hand panel of our Fig. 15 shows a power-law profile for the toroidal field (solid blue line) such that α = −0.5 and vAϕ0 = 1.41 (see equation 35), in order to match the scaling found by Bromberg & Tchekhovskoy (2016). The ratio of the vertical to toroidal field strength is roughly similar to that in Fig. 14 of Bromberg & Tchekhovskoy (2016), such that we choose vAz = 0.14 (dashed green line). Note that our vertical field is constant, unlike that of Bromberg & Tchekhovskoy (2016) (we shall comment on this later). In order to compare the thermal pressure with the magnetic fields, we plot the equivalent quantity given by equation (43), i.e. |$\sqrt{2\rho }$| (red dot–dashed line). We see that for this case, the thermal pressure is in exact equipartition with the toroidal field, and very nearly in equipartition with the total magnetic field, |$v_{\mathrm{ Atot}} = \sqrt{v_{A\phi }^2 + v_{Az}^2}$| (grey solid line). Since our computational domain is located away from the jet boundaries and core, this plot can be compared with the toroidally dominated regime in fig. 14 of Bromberg & Tchekhovskoy (2016), i.e. the outer jet region from ∼20 to 60RL, where RL is the light cylinder radius in their work [we have kept the same colour and line-style scheme in the top left-hand panel of our figure as fig. 14 of Bromberg & Tchekhovskoy (2016) for easier comparison]. The main point of this exercise is to demonstrate that despite equipartition, as well as the r−0.5 scaling for vAϕ, we find unstable m = 1 modes in our analysis (with a non-negligible maximum growth rate ∼0.3), unlike Bromberg & Tchekhovskoy (2016). Note that the bottom panel of Fig. 15 presents solutions for two resolutions Nr = 200 and Nr = 256. This is to show that the apparent gaps in the solutions towards higher |kz| are not real, and they close as the resolution increases. The eigenfunctions of the most unstable modes at kz = −2, −4, and −12 are plotted in the top right-hand panel. We shall now discuss a few possible explanations for this apparent discrepancy. First, we check whether the above eigenfunctions can be resolved in the set-up employed by Bromberg & Tchekhovskoy (2016), who use a logarithmic radial grid with uniform (dimensionless) spacing rgrid = (rbr/rin)1/N − 1, where N is their radial grid resolution and rin = 0.8RL and rbr = 800RL are the inner and outer radial boundaries (see their table 1). They use two resolutions, namely N = 128 and N = 256, which correspond to rgrid ∼ 0.06 and 0.03, respectively. The dimensionless quantity describing our eigenfunctions, equivalent to rgrid, is therefore Δr/r1 (see the top left-hand panel of Fig. 2). We see that as |kz| increases, the eigenfunctions in the top right-hand panel of Fig. 15 become more localized towards the inner boundary having Δr/r1 ∼ 0.7, 0.4, and 0.24 for kz = −2, −4, and −12, respectively. We next define a quality factor |${\cal Q} = (\Delta r/r_1)/r_{\mathrm{ grid}}$|⁠, which determines how well resolved an eigenfunction is. Although fig. 14 of Bromberg & Tchekhovskoy (2016) is for the N = 128 case, we compare our eigenfunctions with their higher-resolution case N = 256. We find that for rgrid ∼ 0.03 the eigenfunctions appear to be well resolved, having |${\cal Q} =23,13, \mathrm{ and}\, 8$| for kz = −2, −4, and −12, respectively. Thus, in principle, the unstable modes we find in our Fig. 15 should have led to the internal kink instability in the 3D jets simulated by Bromberg & Tchekhovskoy (2016). Next, we take a closer look at fig. 14 of Bromberg & Tchekhovskoy (2016) and point out its differences with our Fig. 15. We observe from fig. 14 of Bromberg & Tchekhovskoy (2016) that although the very inner core of the jet (which we exclude from our computation) is indeed in equipartition, the outer jet region, ∼20−60RL, is not. In fact, the thermal pressure dominates over the magnetic pressure in the outer region. Interestingly, the thermal pressure profile is nearly flat in this region, which in turn implies a very weak thermal gradient. This is probably a consequence of the rapid thermal pressure build-up due to magnetic energy dissipation. Moreover, the vertical field has a non-negligible gradient, and its strength, although less than that of the toroidal field, cannot be considered as sub-dominant (hence the B98 stability criterion cannot be applied in this case). Thus, we believe that the main cause for stability seen in Bromberg & Tchekhovskoy (2016) is the negligible thermal pressure gradient in the jet region above the dissipation zone, which is easily stabilized by the combined vertical and toroidal field gradients present therein (as discussed in Section 5.2). Our stability analysis is also consistent with the ‘core–envelope’ jet model of Komissarov (1999), which was used to carry out periodic box, 3D simulations by Porth & Komissarov (2015). This model consists of a Z-pinched2 inner jet core, where the toroidal field balances the plasma pressure, and a force-free, purely toroidal outer envelope such that Bϕ ∝ 1/r. Note that the main objective of Porth & Komissarov (2015) was to investigate the stability of their jet to external, large-scale instabilities by exploring different pressure distributions for the ambient medium. They found that for a critical ambient pressure (Pamb ∝ z−2) the overall jet becomes stable to the external instabilities, the outer envelope is featureless, and the internal instabilities are confined to a narrow jet core. We point out that irrespective of the nature of the ambient medium and the jet core (both of which are excluded from our analysis), the outer envelope would be stable to all internal instabilities according to the B98 stability criterion (see Section 3.2), which has been extensively verified in this work (see Section 5.1). Thus, if an unstable jet core additionally develops, the magnetic energy would be dissipated to heat by the internal instabilities therein. This would lead to a bright core or spine, which would in turn be surrounded by a stable and darker envelope or sheath, giving rise to the so-called spine–sheath structure, commonly invoked to explain observations of jets in AGNs (e.g. Mimica et al. 2015) and blazars (e.g. Sikora, Rutkowski & Begelman 2016). In an earlier study, Mizuno et al. (2012) also carried out 3D relativistic, periodic box MHD simulations of rotating jets, having a poloidal-field-dominated core and Bϕ → 1/r in the envelope. They considered a low gas density and pressure background decreasing with radius such that their jets were still force free. They found that the jets were helically distorted due to the development of the internal kink instability, whose growth rate depended strongly on the gradient of the poloidal field. Their findings were in accord with the linear stability analyses mentioned in our Introduction (Istomin & Pariev 1996; Lyubarskii 1999), which in turn emphasizes the relevance of linear studies. Mizuno et al. (2012) further found that multiple unstable wavelengths were excited during the non-linear evolution, which coupled to give rise to an external kink instability that disrupted the entire jet. We consider fig. 8 of Mizuno et al. (2012), which shows the time evolution of 2D density slices in the xz plane of the jet (z being the jet axis), for different angular velocities. If we interpret the loop-like high-density structures in this figure (denoted by orange and red) as locations where the internal kink instability is triggered, we observe that they are localized in the interior of the jet (for all the panels). However, there is no instability seen towards the edge of the jet where Bϕ ∝ 1/r, which is again compatible with the B98 stability criterion. Finally, we discuss the 3D resistive MHD simulations by Striani et al. (2016) to highlight the significance of thermal pressure gradients. These were performed to study reconnection in magnetically confined cylindrical plasma columns. Striani et al. (2016) modelled both a force-free equilibrium and one with a thermal pressure gradient. For the latter case, they found small-scale finger-like structures in the current density triggered by the pressure-driven internal kink instability (see their fig. 7). These features were shown to cause efficient magnetic dissipation via reconnection of multiple fragmented current sheets, which in turn reinforces the premise of our work. For the force-free case, Striani et al. (2016) observed a bending of the plasma column due to the current-driven kink instability as well as the formation of plume-like structures (see their fig. 9). More interestingly, they found that the currents through several plumes were perpendicular to the magnetic field, indicative of a secondary pressure-driven instability (as the Lorentz force is required to be balanced by a non-zero thermal pressure gradient in these regions; see their fig. 10), which led to the formation of mixed modes. Fragmentation of current sheets and an enhanced magnetic dissipation rate were also noted in this case. This finding leads to the important inference that initially Poynting-flux-dominated jets can also develop thermal pressure gradients, which in turn are essential to facilitate magnetic reconnection and dissipation. 7 SUMMARY AND CONCLUSIONS Below we summarize the main findings of our analysis: We have carried out a global eigenvalue analysis of the linearized, ideal MHD equations, in order to study the stability of cylindrical jets to both axisymmetric and non-axisymmetric perturbations. Our calculations have been performed in the jet fluid frame, and we have ignored velocity shear (both at the surface and in the interior), rotation, and relativistic effects for simplicity. We have focused on characterizing the small-scale internal instabilities that are confined to the interiors of jets, away from the effects of both the jet core and the boundaries. We have further analysed the importance of a thermal pressure gradient for triggering instabilities in a region of the jet dominated by a toroidal magnetic field (with a sub-dominant vertical field), as is likely to occur at some distance from the jet source. This gives rise to predominantly pressure-driven instabilities, which are often overlooked in the literature compared to their purely current-driven counterparts arising in force-free jets. Since the exact transverse structure of magnetic fields inside jets is not known, we have adopted three different background configurations for our stability analysis, yielding the results summarized below.  Before we proceed with the summary, note that we have further classified the internal instabilities as Local and Global, based on the properties of their radial eigenfunctions (see Table 4 for details). The Local instabilities have narrow eigenfunctions and peak very close to the inner radial boundary. Hence, these agree very well with local (i.e. radially localized) analytic calculations. The Global instabilities fall into several different types, any of which can be captured accurately only in a global (radially extended) framework as they are influenced by the background radial gradients and curvature (again, not to be confused with large-scale external instabilities). Note that all the instabilities found in this work are purely growing slow magnetosonic modes. First, we consider a power-law profile for the toroidal magnetic field, such that vAϕ ∝ rα, and a constant vertical magnetic field vAz (Section 5.1; case vAϕ-PL−vAz-Cons). We do a parameter scan for different values of α, and find internal instabilities for a range of vertical wavenumbers kz (both kz > 0 and kz < 0; the sign indicates the direction of the helical twist with respect to the magnetic field if m > 0), for m = 0, 1, and 2 (see Figs 1, 4, and 8). The maximum growth rate (which is nearly independent of m, as is characteristic of pressure-driven instabilities) for a given α matches very well with the predictions of B98, who performed a local, linear stability analysis of a similar jet model. Interestingly, we find that in a region where vAz ≪ vAϕ, all (m, kz) modes are stabilized for α ≤ −1 (see e.g. Figs 3 and 6). This global stability criterion is in excellent agreement with the local criterion of B98. This criterion is better understood by looking at the radial equilibrium equation (9), which shows that the thermal pressure gradient (or, equivalently, density gradient in our case) has a negative sign in general (for α > −1), and is oppositely directed to the restoring magnetic tension force. However, when α = −1, the thermal pressure gradient vanishes, the toroidal magnetic pressure gradient is balanced by its own tension force, and the system becomes stable irrespective of m and kz. Note that for m = 1, the kz > 0 modes get stabilized at a shallower toroidal field profile, i.e. α ≲ −0.5, which is also in direct corroboration with B98. In a region where vAz is comparable to vAϕ, however, the modes get stabilized even for α > −1 (Fig. 7).  The m = 0 unstable modes can be classified as Local instabilities across almost their entire kz range (Fig. 2). The global solutions for the m ≠ 0 cases, however, exhibit a mixture of Local and Global instabilities. The Global instabilities tend to occur at smaller |kz|, and also include those arising due to the interaction between distinct slow modes occupying the same range of kz (Figs 5 and 9). Nevertheless, the fastest-growing modes for all m are Local instabilities, which occur predominantly at large negative values of kz. This is the reason for the excellent agreement between B98’s local analysis and our global solutions. The important consequence of this finding is as follows: In a toroidally dominated region of a jet with a negative thermal pressure gradient strong enough to trigger instability (and very weak vertical field), internal instabilities can occupy any local patch of the jet, irrespective of the exact configuration of the magnetic field in the jet interior. According to B98, the growth time-scales of these instabilities are inversely proportional to the toroidal field strength; thus, they grow fast enough for a substantial region of the jet to become internally unstable. This in turn can serve as a likely site for dissipation of magnetic energy to heat (and non-thermally accelerated particles) via current sheet formation. The Global instabilities, which are subject to the details of the radial profiles, would further enhance this effect. Although internal instabilities can develop close to the base of Poynting-flux-dominated jets (e.g. Bromberg & Tchekhovskoy 2016), we propose that a thermal pressure gradient is crucial to sustain them farther from the source, where much of the particle acceleration leading to radiation takes place. Next, we consider the same toroidal field profile as above but with a radially varying vertical magnetic field (Section 5.2; case vAϕ-PL−vAz-Var). In order to self-consistently determine the equilibrium, we assume the vertical magnetic pressure to be proportional to the thermal pressure, characterized by the parameter ϵ (see equation 38) – the smaller the value of ϵ, the stronger the vertical magnetic field. As the vertical field gradient increases (for a given α), there is a consequent increase in the vertical field strength relative to the toroidal magnetic field, although the thermal pressure gradient remains unchanged (as it is independent of ϵ in this formalism). A decrease in ϵ induces stabilization of both the m = 0 and m = 1 modes for a given α (Fig. 10). Complete stabilization for the m = 0 and m = 1 modes is, however, attained at different values of ϵ, indicating that the pinch and kink modes react differently to the background gradients. Moreover, our results indicate that stability in such a situation is not determined solely by whether equipartition is achieved between the thermal and magnetic pressures (see the last point below). Finally, we consider a more complex toroidal field profile (see equation 44) with a constant (weak) vertical magnetic field (Section 5.3; case vAϕ-Gen−vAz-Cons). Interestingly, we find that for all the m = 0, 1, and 2 cases, the growth rate and the kz range of instability of the most unstable global solutions show excellent agreement with the corresponding local solutions obtained for a power-law toroidal field with α = 0.8 and the same vAϕ0 (see Figs 11 and 12). We also establish a direct correlation between the stability criterion for the complex toroidal field configuration and the simple power-law configuration. We find that as the background toroidal field becomes proportional to r−1/2, i.e. α = −0.5, the kz > 0 modes get completely stabilized for m = 1, as in the power-law case (see Fig. 14). A similar line of reasoning allows us to conclude that the B98 stability criterion for all kz modes (i.e. α ≤ −1) would also be valid be in this case. The m = 0 modes again represent mostly Local instabilities, whereas the m ≠ 0 modes exhibit both Global and Local instabilities, occupying different kz. As kz assumes increasingly negative values, however, the instabilities become more and more localized close to the inner radial boundary for the non-axisymmetric case (see Fig. 13). This further validates our claim that highly localized internal instabilities are triggered as long as the B98 instability criterion is satisfied locally inside the jet, irrespective of the background magnetic field profile. We have also established that magnetic resonance does not play any role in triggering the m ≠ 0 instabilities (irrespective of the toroidal field profiles), as their eigenfunctions peak far away from the corresponding resonant radii (see Section 3.3, and second panels of Figs 5, 9, and 13). We have also compared our results with those of recent 3D numerical simulations by Mizuno et al. (2012), Porth & Komissarov (2015), Bromberg & Tchekhovskoy (2016), and Striani et al. (2016). We conclude that thermal pressure gradients are indeed key to incite internal instabilities, which are fundamental for causing magnetic energy dissipation via reconnection of current sheets. With regard to the question of what causes internal instabilities to saturate, our results indicate two possible answers depending upon the region of the jet in question: (1) In a toroidally dominated region with a very weak vertical field, complete stability for all m modes is attained when the toroidal field behaves as vAϕ ∝ 1/r, as proposed by B98. This supports the ‘spine–sheath’ model of jets, with a Z-pinched jet core and an envelope stable to internal instabilities, observed in several AGN jets; (2) in a region with an appreciable toroidal field but a non-negligible vertical field strength and radial gradient, stability ultimately ensues from the complex interplay between the destabilizing thermal pressure gradient and the stabilizing magnetic pressure gradients, or equivalently, magnetic shear. Contrary to the claim by Bromberg & Tchekhovskoy (2016), we conclude that equipartition alone may not guarantee stability (see Fig. 15). Both of these situations are possible in a real jet (at different distances from the source), depending upon how fast the vertical field falls off relative to the toroidal field due to either the lateral expansion of the jet (Begelman et al. 1984) or the divergence of outer flux surfaces (Begelman & Li 1994). ACKNOWLEDGEMENTS This work was supported in part by NASA (National Aeronautics and Space Administration) Astrophysics Theory Program grant NNX14AB37G and NSF (National Science Foundation) grant AST-1411879. UD is immensely grateful to Ben Brown for allocating time on the LCD machines, where all the calculations using dedalus were performed. UD specially thanks Phil Armitage, Bhupendra Mishra, and Vladimir Zhdankin for general discussions. The authors thank the anonymous referee for their constructive comments. Footnotes 1 Dedalus is available at http://dedalus-project.org. 2 A Z-pinch is a cylindrical geometry in plasma physics, where an electric current flows in the z-direction that generates a magnetic field in the ϕ-direction. REFERENCES Appl S. , 1996 , A&A , 314 , 995 Appl S. , Lery T. , Baty H. , 2000 , A&A , 355 , 818 Bateman G. , 1978 , MHD Instabilities . MIT Press , Cambridge, Mass Baty H. , Keppens R. , 2002 , ApJ , 580 , 800 https://doi.org/10.1086/343893 Crossref Search ADS Begelman M. C. , 1998 , ApJ , 493 , 291 https://doi.org/10.1086/305119 (B98) Crossref Search ADS Begelman M. C. , Li Z.-Y. , 1992 , ApJ , 397 , 187 https://doi.org/10.1086/171778 Crossref Search ADS Begelman M. C. , Li Z.-Y. , 1994 , ApJ , 426 , 269 https://doi.org/10.1086/174061 Crossref Search ADS Begelman M. C. , Blandford R. D. , Rees M. J. , 1984 , Rev. Mod. Phys. , 56 , 255 https://doi.org/10.1103/RevModPhys.56.255 Crossref Search ADS Beskin V. S. , 2010 , Phys.-Usp. , 53 , 1199 https://doi.org/10.3367/UFNe.0180.201012b.1241 Crossref Search ADS Blandford R. D. , Znajek R. L. , 1977 , MNRAS , 179 , 433 https://doi.org/10.1093/mnras/179.3.433 Crossref Search ADS Bodo G. , Mamatsashvili G. , Rossi P. , Mignone A. , 2013 , MNRAS , 434 , 3030 https://doi.org/10.1093/mnras/stt1225 Crossref Search ADS Bodo G. , Mamatsashvili G. , Rossi P. , Mignone A. , 2016 , MNRAS , 462 , 3031 https://doi.org/10.1093/mnras/stw1650 Crossref Search ADS Bonanno A. , Urpin V. , 2011a , Phys. Rev. E , 84 , 056310 https://doi.org/10.1103/PhysRevE.84.056310 Crossref Search ADS Bonanno A. , Urpin V. , 2011b , A&A , 525 , A100 https://doi.org/10.1051/0004-6361/200913836 Crossref Search ADS Bromberg O. , Tchekhovskoy A. , 2016 , MNRAS , 456 , 1739 https://doi.org/10.1093/mnras/stv2591 Crossref Search ADS Das U. , Begelman M. C. , Lesur G. , 2018 , MNRAS , 473 , 2791 https://doi.org/10.1093/mnras/stx2518 Crossref Search ADS de Gouveia Dal Pino E. M. , 2005 , Adv. Space Res. , 35 , 908 https://doi.org/10.1016/j.asr.2005.03.145 Crossref Search ADS de Gouveia Dal Pino E. M. , Kowal G. , 2015 , in Lazarian A. , de Gouveia Dal Pino E. M. , Melioli C. , eds, Magnetic Fields in Diffuse Media, Vol. 407 . Springer-Verlag , Berlin Heidelberg , p. 373 Freidberg J. P. , 1982 , Rev. Mod. Phys. , 54 , 801 https://doi.org/10.1103/RevModPhys.54.801 Crossref Search ADS Guo F. et al. , 2016 , ApJ , 818 , L9 https://doi.org/10.3847/2041-8205/818/1/L9 Crossref Search ADS Hardee P. E. , 1979 , ApJ , 234 , 47 https://doi.org/10.1086/157471 Crossref Search ADS Hardee P. E. , 2011 , in Romero G. E. , Sunyaev R. A. , Belloni T. , eds, Proc. IAU Symp. 275, Jets at All Scales . Cambridge Univ. Press , Cambridge , p. 41 Hardee P. E. , 2013 , EPJ Web Conf. , 61 , 02001 Crossref Search ADS Istomin Y. N. , Pariev V. I. , 1996 , MNRAS , 281 , 1 https://doi.org/10.1093/mnras/281.1.1 Crossref Search ADS Kagan D. , Sironi L. , Cerutti B. , Giannios D. , 2015 , Space Sci. Rev. , 191 , 545 https://doi.org/10.1007/s11214-014-0132-9 Crossref Search ADS Kim J. , Balsara D. S. , Lyutikov M. , Komissarov S. S. , George D. , Siddireddy P. K. , 2015 , MNRAS , 450 , 982 https://doi.org/10.1093/mnras/stv606 Crossref Search ADS Komissarov S. S. , 1999 , MNRAS , 308 , 1069 https://doi.org/10.1046/j.1365-8711.1999.02783.x Crossref Search ADS Longaretti P.-Y. , 2008 , in Massaglia S. , Bodo G. , Mignone A. , Rossi P. , eds, Lecture Notes in Physics, Vol. 754, Jets From Young Stars III . Springer Verlag , Berlin , p. 131 Lyubarskii Y. E. , 1999 , MNRAS , 308 , 1006 https://doi.org/10.1046/j.1365-8711.1999.02763.x Crossref Search ADS Lyutikov M. , Blandford R. , 2003 , preprint (arXiv:e-print) arXiv:astro-ph/0312347 McKinney J. C. , Blandford R. D. , 2009 , MNRAS , 394 , L126 https://doi.org/10.1111/j.1745-3933.2009.00625.x Crossref Search ADS Mignone A. , Rossi P. , Bodo G. , Ferrari A. , Massaglia S. , 2010 , MNRAS , 402 , 7 https://doi.org/10.1111/j.1365-2966.2009.15642.x Crossref Search ADS Mimica P. , Giannios D. , Metzger B. D. , Aloy M. A. , 2015 , MNRAS , 450 , 2824 https://doi.org/10.1093/mnras/stv825 Crossref Search ADS Mizuno Y. , Lyubarsky Y. , Nishikawa K.-I. , Hardee P. E. , 2009 , ApJ , 700 , 684 https://doi.org/10.1088/0004-637X/700/1/684 Crossref Search ADS Mizuno Y. , Lyubarsky Y. , Nishikawa K.-I. , Hardee P. E. , 2011 , ApJ , 728 , 90 https://doi.org/10.1088/0004-637X/728/2/90 Crossref Search ADS Mizuno Y. , Lyubarsky Y. , Nishikawa K.-I. , Hardee P. E. , 2012 , ApJ , 757 , 16 https://doi.org/10.1088/0004-637X/757/1/16 Crossref Search ADS Mizuno Y. , Hardee P. E. , Nishikawa K.-I. , 2014 , ApJ , 784 , 167 https://doi.org/10.1088/0004-637X/784/2/167 Crossref Search ADS Moll R. , 2009 , A&A , 507 , 1203 https://doi.org/10.1051/0004-6361/200912266 Crossref Search ADS Nalewajko K. , Begelman M. C. , 2012 , MNRAS , 427 , 2480 https://doi.org/10.1111/j.1365-2966.2012.22117.x Crossref Search ADS Narayan R. , Li J. , Tchekhovskoy A. , 2009 , ApJ , 697 , 1681 https://doi.org/10.1088/0004-637X/697/2/1681 Crossref Search ADS O’Neill S. M. , Beckwith K. , Begelman M. C. , 2012 , MNRAS , 422 , 1436 https://doi.org/10.1111/j.1365-2966.2012.20721.x Crossref Search ADS Petropoulou M. , Giannios D. , Sironi L. , 2016 , MNRAS , 462 , 3325 https://doi.org/10.1093/mnras/stw1832 Crossref Search ADS Porth O. , Komissarov S. S. , 2015 , MNRAS , 452 , 1089 https://doi.org/10.1093/mnras/stv1295 Crossref Search ADS Porth O. , Komissarov S. S. , Keppens R. , 2014 , MNRAS , 438 , 278 https://doi.org/10.1093/mnras/stt2176 Crossref Search ADS Sikora M. , Begelman M. C. , Madejski G. M. , Lasota J.-P. , 2005 , ApJ , 625 , 72 https://doi.org/10.1086/429314 Crossref Search ADS Sikora M. , Rutkowski M. , Begelman M. C. , 2016 , MNRAS , 457 , 1352 https://doi.org/10.1093/mnras/stw107 Crossref Search ADS Singh C. B. , Mizuno Y. , de Gouveia Dal Pino E. M. , 2016 , ApJ , 824 , 48 https://doi.org/10.3847/0004-637X/824/1/48 Crossref Search ADS Sironi L. , Petropoulou M. , Giannios D. , 2015 , MNRAS , 450 , 183 https://doi.org/10.1093/mnras/stv641 Crossref Search ADS Sobacchi E. , Lyubarsky Y. E. , 2018 , MNRAS , 473 , 2813 https://doi.org/10.1093/mnras/stx2592 Crossref Search ADS Sobacchi E. , Lyubarsky Y. E. , Sormani M. C. , 2017 , MNRAS , 468 , 4635 https://doi.org/10.1093/mnras/stx807 Crossref Search ADS Striani E. , Mignone A. , Vaidya B. , Bodo G. , Ferrari A. , 2016 , MNRAS , 462 , 2970 https://doi.org/10.1093/mnras/stw1848 Crossref Search ADS Suydam B. R. , 1958 , in Proceedings of the Second International Conference on the Peaceful Uses of Atomic Energy, Vol. 31 . United Nations , Geneva , p. 157 Tomimatsu A. , Matsuoka T. , Takahashi M. , 2001 , Phys. Rev. D , 64 , 123003 https://doi.org/10.1103/PhysRevD.64.123003 Crossref Search ADS Werner G. R. , Uzdensky D. A. , Begelman M. C. , Cerutti B. , Nalewajko K. , 2018 , MNRAS , 473 , 4840 https://doi.org/10.1093/mnras/stx2530 Crossref Search ADS Zhdankin V. , Werner G. R. , Uzdensky D. A. , Begelman M. C. , 2017 , Phys. Rev. Lett. , 118 , 055103 https://doi.org/10.1103/PhysRevLett.118.055103 Crossref Search ADS PubMed APPENDIX A: LINEARIZED MHD EQUATIONS In order to carry out a linear stability analysis, we perturb the MHD equations (1), (2), (4), and (5), and retain only the terms of linear order such that \begin{eqnarray*} \rho = \rho _0 + \rho _1 , \end{eqnarray*} (A1) \begin{eqnarray*} P = P_0 + P_1 , \end{eqnarray*} (A2) \begin{eqnarray*} \boldsymbol{ v} = {\boldsymbol{v}_1} , \end{eqnarray*} (A3) \begin{eqnarray*} {\boldsymbol{B}} = {\boldsymbol{B}_0} + {\boldsymbol{B}_1} , \end{eqnarray*} (A4) \begin{eqnarray*} {\boldsymbol{A}} = {\boldsymbol{A}_0} + {\boldsymbol{A}_1} . \end{eqnarray*} (A5) We can now write the complete set of linearized MHD equations for a generic case involving non-axisymmetric perturbations, background radial gradients, compressibility, and radial curvature: \begin{eqnarray*} \partial _t \rho _1 + \rho _0 \biggl [ \frac{1}{r} \partial _r (r v_{1r}) + \frac{1}{r} \partial _\phi v_{1\phi } + \partial _z v_{1z} \biggr ] + v_{1r} \partial _r \rho _0 = 0 , \end{eqnarray*} (A6) \begin{eqnarray*} &&{\rho _0 \partial _t v_{1r} + \partial _r \rho _1 + \frac{1}{4 \pi } \biggl [ B_{0\phi } \partial _r B_{1 \phi } + B_{1\phi }\partial _r B_{0\phi } - \frac{B_{0\phi }}{r} \partial _\phi B_{1r}} &&{ - B_{0z} \partial _z B_{1r} + B_{0z} \partial _r B_{1z} + \frac{2 B_{0\phi } B_{1\phi }}{r} + B_{1z} \partial _r B_{0z} \biggr ] = 0 , } \end{eqnarray*} (A7) \begin{eqnarray*} &&{\rho _0 \partial _t v_{1\phi } + \frac{1}{r} \partial _\phi \rho _1 - \frac{1}{4\pi } \biggl [ B_{1r}\partial _r B_{0\phi } + B_{0z} \partial _z B_{1\phi } - \frac{B_{0z}}{r} \partial _\phi B_{1z}}\nonumber\\ &&{\qquad + \frac{B_{0\phi } B_{1r}}{r} \biggr ] = 0 , } \end{eqnarray*} (A8) \begin{eqnarray*} \rho _0 \partial _t v_{1z} + \partial _z \rho _1 + \frac{1}{4\pi } \biggl [ B_{0\phi } \partial _z B_{1\phi } - \frac{B_{0\phi }}{r}\partial _\phi B_{1z} - B_{1r} \partial _r B_{0z} \biggr ] = 0 , \end{eqnarray*} (A9) \begin{eqnarray*} B_{1\phi } - \partial _z A_{1r} + \partial _r A_{1z} = 0 , \end{eqnarray*} (A10) \begin{eqnarray*} B_{1z} + \partial _r A_{1\phi } - \frac{1}{r} \partial _r (r A_{1r}) = 0 , \end{eqnarray*} (A11) \begin{eqnarray*} \partial _t A_{1r} - B_{0z} v_{1\phi } + B_{0\phi } v_{1z} = 0 , \end{eqnarray*} (A12) \begin{eqnarray*} \partial _t A_{1\phi } + B_{0z} v_{1r} = 0 , \end{eqnarray*} (A13) \begin{eqnarray*} \partial _t A_{1z} - B_{0\phi } v_{1r} = 0 , \end{eqnarray*} (A14) where the relation |$P_1 = c_s^2 \rho _1$| has been used to eliminate P1. Note that B1r is not an independent quantity as it can be defined as |$B_{1r} = \partial _\phi A_{1z}/r - \partial _z A_{1\phi }$|⁠. This set of nine equations for the nine perturbed variables is obtained by imposing the Weyl gauge, ψ = 0. However, we have also solved these equations for a non-zero ψ by imposing the Coulomb gauge, ∇ · A = 0 (which serves as the 10th equation for the 10th variable ψ). Both gauges lead to identical results as expected. © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Internal instabilities in magnetized jets
JF - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/sty2675
DA - 2019-01-11
UR - https://www.deepdyve.com/lp/oxford-university-press/internal-instabilities-in-magnetized-jets-iKEfCcsN8s
SP - 2107
VL - 482
IS - 2
DP - DeepDyve
ER -