THREE-DIMENSIONAL RELATIVISTIC MAGNETOHYDRODYNAMIC SIMULATIONS OF CURRENT-DRIVEN INSTABILITY. I. INSTABILITY OF A STATIC COLUMN

Relativistic jets occur in active galactic nuclei (AGNs; e.g., Urry & Padovani 1995; Ferrari 1998; Meier et al. 2001), in microquasars (e.g., Mirabel & Rodriguez 1999), and are thought responsible for the gamma-ray bursts (GRBs; e.g., Zhang & Mészáros 2004; Piran 2005; Mészáros 2006). The most promising mechanism for producing relativistic jets involves magnetohydrodynamic (MHD) acceleration from an accretion disk around a black hole (e.g., Blandford 1976; Lovelace 1976; Blandford & Payne 1982; Fukue 1990; Meier 2005; Narayan et al. 2007), and/or involves the extraction of energy from a rotating black hole (Penrose 1969; Blandford & Znajek 1977).

General relativistic MHD (GRMHD) codes have been used to study the extraction of rotational energy from a spinning black hole, i.e., Blandford-Znajek mechanism (Koide 2003; Komissarov 2005; McKinney 2005; Komissarov & McKinney 2007; Komissarov & Barkov 2009) and from an accretion disk, i.e., Blandford-Payne mechanism (Blandford & Payne 1982). GRMHD codes also probe the formation of GRB jets in collapsars (e.g., Mizuno et al. 2004a, 2004b; Liu et al. 2007; Barkov & Komissarov 2008; Stephens et al. 2008; Nagataki 2009; Komissarov & Barkov 2009) and are used to study the universal magnetic jet acceleration mechanism by the propagation of jets injected from an unresolved Keplerian disk (e.g., Komissarov et al. 2007, 2009; Tchekhovskoy et al. 2008, 2009). GRMHD simulations with a spinning black hole indicate jet production with a magnetically dominated high Lorentz factor spine with v ∼  c, and a matter-dominated sheath with v ≳  c/2 possibly embedded in a lower speed, v ≪ c, disk/coronal wind (e.g., De Villiers et al. 2003, 2005; Hawley & Krolik 2006; McKinney & Gammie 2004; McKinney 2006; McKinney & Narayan 2007a, 2007b; McKinney & Blandford 2009; Beckwith et al. 2008; Hardee et al. 2007; Komissarov & Barkov 2009). Circumstantial evidence such as the requirement for large Lorentz factors suggested by the TeV BL Lacs when contrasted with much slower observed motions (Ghisellini et al. 2005) suggests such a spine-sheath morphology, although alternative interpretations are also possible (Georganopoulos & Kazanas 2003; Bromberg & Levinson 2009; Stern & Poutanen 2008; Giannios et al. 2009).

Numerical simulations indicate that the jet spine and sheath are accelerated and, in part, collimated by strong magnetic fields twisted in the rotating black hole ergosphere and in the accretion disk, respectively. The helically twisted magnetic fields are expected to be dominated by the toroidal component in the far zone as a result of jet expansion because the poloidal component falls off faster with expansion and distance. In configurations with strong toroidal magnetic field, the current-driven (CD) kink mode is unstable. This instability excites large-scale helical motions that can strongly distort or even disrupt the system. For static cylindrical force-free equilibria, the well-known Kruskal–Shafranov criterion states that the instability develops if the length of the column, ℓ, is long enough for the field lines to go around the cylinder at least once (e.g., Bateman 1978): |B_p/B_ϕ| < ℓ/2πR. This criterion suggests that jets are unstable beyond the Alfvén surface because of B_p ≲  B_ϕ and ℓ ≳  R_j at the Alfvén surface, where R_j is the cylindrical radius of the jet. However, rotation and shear motions could significantly affect the instability criterion. For relativistic force-free configurations, the linear instability criteria were studied by Istomin & Pariev (1994, 1996), Begelman (1998), Lyubarskii (1999), Tomimatsu et al. (2001), and Narayan et al. (2009).

The linear mode analysis provides conditions for the instability but says little about the impact the instability has on the system. The instability of the potentially disruptive kink mode found from a linear analysis must be followed into the nonlinear regime. The nonlinear development of the CD kink instability is different for external and internal kink modes (e.g., Bateman 1978). External disruptive instability develops in a current-carrying plasma column surrounded by a vacuum magnetic field. In astrophysical systems, some amount of plasma is present everywhere. In this case, an internal kink instability develops inside a resonant surface on which the helicity of the magnetic field lines matches that of the eigenmode. In the dissipationless MHD, the internal instability is known to saturate. Finite resistivity makes the internal kink mode disruptive, but the timescale significantly exceeds the Alfvén timescale and may be too long to disrupt a fast outflow. In this case, helical structures may be formed in the flow. Recently, helical structures have been found in nonrelativistic/relativistic simulations of magnetized jets (e.g., Lery et al. 2000; Ouyed et al. 2003; Nakamura & Meier 2004; Nakamura et al. 2007; Moll et al. 2008; McKinney & Blandford 2009; Carey & Sovinec 2009).

Twisted structures are observed in many AGN jets on subparsec, parsec, and kiloparsec scales (e.g., Gómez et al. 2001; Lobanov & Zensus 2001). Dissipation processes may result in relaxation to helical equilibrium as evoked by Königl & Choudhuri (1985) for helical structures in AGN jets. In the absence of CD kink instability and resistive relaxation, helical structures may be attributed to the Kelvin–Helmholtz (KH) helical instability driven by velocity shear at the boundary between the jet and the surrounding medium (e.g., Hardee 2004, 2007) or triggered by precession of the jet ejection axis (Begelman et al. 1980). It is still not clear whether current driven, velocity shear driven, or jet precession is responsible for the observed structures, or whether these different processes are responsible for the observed twisted structures at different spatial scales.

This is the first in a series of papers in which we study kink instability in relativistic systems. By relativistic, we mean not only relativistically moving systems but any with magnetic energy density comparable to or greater than the plasma energy density, including the rest mass energy. In this paper, we present three-dimensional results of the CD kink instability of a static plasma column. We start from static configurations because in the case of interest, the source of the free energy is the magnetic field, not the kinetic energy. Therefore, static configurations (or more generally rigidly moving flows considered in the proper reference frame) are the simplest ones for studying the basic properties of the kink instability. At the next stage, we will investigate the influence of shear motions and rotation on the stability and nonlinear behavior of jets. This paper is organized as follows. We describe the numerical method and setup used for our simulations in Section 2, present our results in Section 3, discuss the astrophysical implications in Section 4, and present some numerical tests in the Appendix.

In order to study time evolution of the CD kink instability in the relativistic MHD (RMHD) regime, we use the three-dimensional GRMHD code "RAISHIN" in Cartesian coordinates. RAISHIN is based on a 3 + 1 formalism of the general relativistic conservation laws of particle number and energy momentum, Maxwell's equations, and Ohm's law with no electrical resistance (ideal MHD condition) in a curved spacetime (Mizuno et al. 2006). In the RAISHIN code, a conservative, high-resolution shock-capturing scheme is employed. The numerical fluxes are calculated using the Harten, Lax, & van Leer (HLL) approximate Riemann solver, and flux-interpolated, constrained transport is used to maintain a divergence-free magnetic field (Harten et al. 1983). The RAISHIN code has proven to be accurate to second order and has passed a number of numerical tests including highly relativistic cases and highly magnetized cases in both special and general relativity (Mizuno et al. 2006). The RAISHIN code performs special relativistic calculations in Minkowski spacetime by changing the metric.

For our simulations, we will choose a force-free helical magnetic field as the initial configuration. A force-free configuration is a reasonable choice for the strong magnetic field cases that we study here. In general, the force-free equilibrium of a static cylinder is described by the equation

$$\begin{equation} B_{z} {d B_{z} \over dR} + {B_{\phi } \over R} {d B_{\phi } \over dR} =0. \end{equation} \tag{ 1 }$$

In particular, we choose a poloidal magnetic field component of the form

$$\begin{equation} B_{z}= {B_{0} \over [1+ (R/a)^{2}]^{\alpha }}, \end{equation} \tag{ 2 }$$

for which one finds a toroidal magnetic field component of the form

$$\begin{equation} B_{\phi }= {B_{0} \over (R/a)[1+ (R/a)^{2}]^{\alpha }} \sqrt{ { [1 + (R/a)^{2}]^{2 \alpha } -1 - 2 \alpha (R/a)^2 \over 2 \alpha -1}}, \end{equation} \tag{ 3 }$$

where R is the cylindrical radial position normalized by a simulation scale unit L ≡ 1, B₀ parameterizes the magnetic field amplitude, a is the characteristic radius of the column, and α is the pitch profile parameter. The pitch profile parameter determines the radial profile of the magnetic pitch P = RB_z/B_ϕ and provides a measure of the twist of the magnetic field lines. With our choice for the force-free field, the magnetic pitch can be written as

$$\begin{equation} P = (R/a)^{2} \sqrt{ { 2 \alpha -1 \over [1 + (R/a)^{2} ]^{2 \alpha } -1 -2 \alpha (R/a)^{2}} }. \end{equation} \tag{ 4 }$$

If the pitch profile parameter 0.5 < α < 1, the magnetic pitch increases with radius. If α>1, the magnetic pitch decreases. When α = 1, the magnetic pitch is constant. This configuration is the same as that used in previous nonrelativistic work (Appl et al. 2000; Baty 2005).

The simulation grid is periodic along the axial direction. As a consequence, the allowed axial wavelengths are restricted to λ = L_z/n ⩽ L_z, with n a positive integer and L_z the grid length. The grid is a Cartesian (x, y, z) box of size 4L × 4L × L_z with grid resolution of ΔL = L/40. The grid resolution is the same in all directions. In simulations, we choose two different column radii, a, relative to the column length, L_z: (cases A) a = (1/16)L_z = (1/8)L and (cases B) a = (1/12)L_z = (1/4)L. For case A L_z = 2L, and for case B L_z = 3L. In terms of a, the simulation box size is 32a × 32a × 16a for case A and 16a × 16a × 12a for case B. We impose outflow boundary conditions on the transverse boundaries at x = y = ±2L (x = y = ±16a for case A and x = y = ±8a for case B).

We consider a low gas pressure medium with constant p = p₀ = 0.02ρ₀c² for the equilibrium state, and two different density profiles: (case u) uniform density ρ = 1.6ρ₀ and (case n) nonuniform density decreasing proportional to the magnetic field strength, ρ = ρ₁B² with ρ₁ = 10.0ρ₀. The equation of state is that of an ideal gas with p = (Γ − 1)ρe, where e is the specific internal energy density and the adiabatic index Γ = 5/3⁵. The specific enthalpy is h ≡ 1 + e/c² + p/ρc². The magnetic field amplitude is $B_{0} =0.4\sqrt{4\pi \rho _{0}c^{2},}$ leading to a low plasma-β near the axis. The sound speed is c_s/c = (Γp/ρh)^1/2, and the Alfvén speed is v_A/c = [B²/(ρh + B²)]^1/2.

Zoom In Zoom Out Reset image size

In order to investigate different radial pitch profiles, we perform simulations with: (case 1) constant pitch α = 1, (case 2) increasing pitch α = 0.75, and (case 3) decreasing pitch α = 2.0. The radial profiles of the magnetic field components, the magnetic pitch, and the sound and Alfvén speeds for the different cases are shown in Figure 1. The radial profile of the magnetic field components in all cases are the same when normalized by a.

The initial MHD equilibrium configuration is perturbed by a small radial velocity with profile given by

$$\begin{equation} v_{R} = \delta v \exp \left(- { R \over R_{a} } \right) \cos (m \theta) \sin \left({2 \pi n z \over L_{z}}\right). \end{equation} \tag{ 5 }$$

The amplitude of the perturbation is taken to be δv = 0.01c with radial width R_a = 0.5L (4a for case A and 2a for case B). We choose m = 1 and n = 1 in the above formula. This is identical to imposing (m, n) = (−1, − 1), because of the symmetry between (m, n) and (−m, − n) pairs. The various different cases that we have considered are listed in Table 1.

3.1. Time Evolution of the CD Kink Instability

As an indicator of the growth of the CD kink instability, we use the volume-averaged kinetic energy transverse to the z-axis within a cylinder of radius R/L ⩽ 1.0 written as

$$\begin{equation} E_{{\rm kin},xy}={1 \over V_{b}} \int _{V_b} {\rho v^{2}_{x} + \rho v^{2}_{y} \over 2} dx dy dz \end{equation} \tag{ 6 }$$

and normalized by the initial volume-averaged transverse magnetic energy. The quantity E_kin,xy allows determination of different evolutionary stages (initial exponential (linear growth phase) growth, and subsequent nonlinear evolution). In all cases, the initial growth regime is characterized by an exponential increase in E_kin,xy by several orders of magnitude to a maximum amplitude followed by a slow decline in the nonlinear regime. The time evolution of the maximum radial velocity evaluated over half the grid length L_z/2 shows a growth trend similar to the growth of E_kin,xy.

Figure 2 shows the time evolution of E_kin,xy and maximum radial velocity for cases A and B with constant pitch for uniform and decreasing radial density profiles. The effect of changing the characteristic radius a relative to the grid length L_z is equivalent to changing the wavelength. The wavelength in case A is λ = 16a, and in case B is λ = 12a. Note that according to the Kruskal–Shafranov criterion, the instability develops at λ>2πa. The instability growth rate reaches a maximum at λ_max ≈ 10a, the exact coefficient being dependent on the transverse distribution of the density and magnetic pitch. Specifically in the case of constant pitch and uniform density, Appl et al. (2000) found λ_max = 8.43a and a corresponding growth rate of Γ_max = 0.133v_A0/a. In general, one can use the estimate Γ_max ≈ 0.1v_A0/a.

Zoom In Zoom Out Reset image size

In uniform density (solid and dashed lines in Figure 2) cases, E_kin,xy grows exponentially (linear growth phase), reaches maximum amplitude at t_S ∼  100t_Au for a = (1/16)L_z and at t_S ∼  90t_Au for a = (1/12)L_z, and slowly declines in the nonlinear regime. The timescale t is in units of the Alfvén crossing time of the characteristic radius a, i.e., t_Au ≡ a/v_A0, where the initial Alfvén velocity on the axis is v_A0 ∼  eq0.3c. The difference in growth is a result of the different wavelengths, longer in case A and shorter in case B. Both cases show similar amplitude of E_kin,xy at transition to the nonlinear stage. Here, we find that the shorter wavelength, λ = 12a, grows slightly more rapidly than the longer wavelength, λ = 16a, as would be expected if λ_max < 12a.

In decreasing density cases (dotted lines and dash-dotted lines in Figure 2), the initial growth in the linear growth phase is more rapid and reaches maximum amplitude at t_S ∼  50t_An for a = (1/16)L_z and for a = (1/12)L_z, and slowly declines in the nonlinear regime. This more rapid growth is a result of dependence of the growth rate on the Alfvén velocity. The density decline in the decreasing case leads to a more gradual radial decline in the Alfvén velocity (see Figure 1(g)) than for the uniform density cases, and on average a higher Alfvén velocity. The maximum amplitude of E_kin,xy relative to the initial volume-averaged magnetic energy for the decreasing density cases is almost the same as the uniform density cases. The maximum radial velocities normalized by the initial axial Alfvén velocity are higher for the decreasing density cases, but in absolute terms are almost the same. We note that the normalized maximum radial velocity in the decreasing density cases appears larger than in the uniform density cases because the normalizing axial Alfvén velocity is smaller for the decreasing density cases.

The effect of different radial pitch on the time evolution of E_kin,xy and the maximum radial velocity for cases with different radial pitch are shown in Figure 3. The three uniform density and three decreasing density cases show similar linear and nonlinear evolution but with different timescales. Growth in the linear stage is more rapid for the decreasing density cases than for the uniform density cases. This can be attributed to the higher average Alfvén velocity in the decreasing density case. For both uniform and decreasing density cases, the increasing pitch case (α = 0.75) grows more slowly and reaches maximum at a later time with a smaller value for E_kin,xy than the constant pitch case (α = 1.0). On the other hand, the decreasing pitch case (α = 2.0) grows more rapidly and reaches maximum at an earlier time with a larger value for E_kin,xy than the constant pitch case. Although the transition time from linear to nonlinear evolution is different for each pitch case, the maximum radial velocity is almost the same at transition. The different growth rates as a function of the radial pitch profile are consistent with the nonrelativistic linear analysis in Appl et al. (2000).

Zoom In Zoom Out Reset image size

3.2. Three-dimensional Structure of the CD Kink Instability

The results of the previous subsection show that the instability excites relatively slow motions, v ≪ v_A, and that the overall energy of the system does not vary too much. Nevertheless, the three-dimensional structure of the system is strongly distorted, as will be shown in this subsection. We note that results for case B with the shorter wavelength are essentially the same as for case A. Moreover, the evolution is similar in the uniform and decreasing density cases and only occurs faster in the decreasing density cases. Therefore, we show the results of case A with decreasing density only.

Figure 4 shows the time evolution of a density isosurface for the constant pitch case A1n. Displacement of the initial force-free helical magnetic field leads to a helically twisted magnetic filament around the density isosurface. At longer times, the radial displacement of the high-density region, best seen in transverse density slices at the grid midplane z = L_z/2 (see the right panels in Figure 4), slows significantly. Continuing outward radial motion is confined to a lower density sheath around the high-density core.

Zoom In Zoom Out Reset image size

The transverse growth of helical twisting is illustrated by the time evolution of magnetic field line displacement seen from the pole-on view shown in Figure 5. Here, we follow field lines anchored initially at R/L = 0.2 at the z = 0 surface. Displacement of the magnetic field rapidly increases with time in the linear stage, and the displacement continues to gradually increase throughout the nonlinear phase. Thus, the magnetic field displacement indicates continued growth, even though the high-density region has ceased significant outward motion by the terminal simulation time. Evolution of the helical perturbation in the magnetic field is accompanied by an appropriate evolution of the current distribution. Initially, a large axial current is located near the axis. In the linear growth phase, the axial current is displaced and twists helically. The axial current twist is cospatial with the high-density twist. In the linear growth phase the axial current gradually decreases, and in the nonlinear phase remains constant in magnitude.

Zoom In Zoom Out Reset image size

Figure 6 shows density isosurfaces and transverse density slices at the grid midplane, z = L_z/2, for the cases A2n (increasing pitch) and A3n (decreasing pitch) that can be compared to Figure 4 for case A1n with constant pitch. The three-dimensional density structure from A2n looks similar to that from A1n. The CD kink instability grows exponentially initially. The transverse density slice at the grid midplane (Figure 6(c)) shows little outward motion of the high-density region in the nonlinear stage similar to the constant pitch case A1n, however, there is little outward motion in the low-density sheath surrounding the high-density region. This result is somewhat different from the constant pitch case A1n and suggests a significant reduction in kink amplitude growth. Results from the decreasing pitch case A3n are very different. Figure 6(b) shows a more slender helical density structure wrapped by the magnetic field. Recall that transition from linear to nonlinear evolution was reached at a much earlier time, t_S ∼  50t_An (see Figure 3), than is shown here. While the density cross section (Figure 6(d)) is similar to that of the constant pitch case A1n (Figure 4(f)), radial motion continues in the nonlinear stage to the end of the simulation.

Zoom In Zoom Out Reset image size

Figure 7 shows the axial current isosurface and the magnetic field line displacement seen from the pole-on view for cases A2n and A3n. Like the constant pitch decreasing density case A1n, the axial current twist is cospatial with the high-density twist. In the nonlinear stage, the helically twisted axial current does not change significantly in case A2n. A pole-on view of field lines anchored initially at R/L = 0.2 at the z = 0 surface (see Figure 7(c)) shows that the magnetic field displacement does not increase significantly in the nonlinear stage. This along with the density results shown in Figures 6(a) and (c) suggest significant nonlinear stabilization of the kink mode in the presence of increasing pitch. On the other hand, the axial current isosurface for the decreasing pitch case A3n seen in Figure 7(b) makes a slender helical structure coincident with the high-density structure. Large displacement in the magnetic field can be seen in the pole-on view (Figure 7(d)).

Zoom In Zoom Out Reset image size

We have investigated the development of the CD kink instability of force-free helical magnetic equilibria in order to study the influence of the initial configuration on the linear and nonlinear evolution. We followed the time evolution where the plasma column was perturbed to excite the m = ±1 azimuthal mode. The kink developed initially as predicted by linear stability analysis. The rate at which the kink developed and, in particular, the nonlinear behavior depended on the density and magnetic pitch radial profile. For a constant magnetic pitch profile, a density decline with radius led to faster linear growth and made a transition to the nonlinear stage sooner than occurred for a uniform density as a result of the more gradual decrease in the Alfvén velocity with radius. While there were some differences in the kink density structure in these two different cases, amplitude growth of the kink continued to the terminal simulation time in both cases.

More profound differences in development accompanied different magnetic pitch profiles. An increase in the magnetic pitch with radius led to a reduced growth rate in the linear (exponential) growth stage and made a transition to nonlinear behavior at a later time when compared to the cases with constant magnetic pitch. In the nonlinear stage, amplitude growth of the kink appeared to have ceased by the terminal simulation time. On the other hand, a decrease in the magnetic pitch led to more rapid growth in the linear growth stage and made a transition to nonlinear behavior at an earlier time when compared to the cases with constant magnetic pitch. In the nonlinear stage, a slender helical density and current-carrying plasma column wrapped by magnetic field developed and amplitude growth continued to increase throughout the simulation. Thus, the linear growth and nonlinear evolution of the CD kink instability depended on the radial density profile and strongly depended on the magnetic pitch profile. Increasing magnetic pitch with radius was clearly stabilizing.

In all cases, the initial axisymmetric structure is strongly distorted by the kink instability, even though not disrupted. It is important to note that the instability develops relatively slowly, that is, the characteristic time for the instability to affect strongly the initial structure is roughly τ ∼  100v_A/a. The growth rate of the kink instability is roughly κ ∼  0.1v_A/c, the exact coefficient being dependent on the structure of the undisturbed state. Therefore, it is not accidental that the full development of the instability takes a long time. In a jet context, our simulations correspond to a perturbation that remains at rest in the flow frame. Therefore, in order to check whether the instability would affect a jet flow, one has to compare τ with the propagation time. In the relativistically moving case, time dilation slows the instability further by the jet Lorentz factor so that the condition for the instability to affect the jet structure may be written as

$$\begin{equation} \frac{z}{c\gamma }> 100\frac{a}{v_A}. \end{equation} \tag{ 7 }$$

In relativistic jets, one could take v_A ≈ c. In order to find the final criterion, one has to know how the jet radius, a, and Lorentz factor, γ, vary with the distance, z. If the jet is narrow enough so that Ωa²/c < z, where Ω is the angular velocity at the base of the jet, one can use the scaling (Tchekhovskoy et al. 2008; Komissarov et al. 2009; Lyubarsky 2009) γ ∼  Ωa. In this case, one finds that the criterion for the kink instability is

$$\begin{equation} zc>100\Omega a^2. \end{equation} \tag{ 8 }$$

This means that the instability could affect the jet structure only if the jet expands slowly enough. Assuming the parabolic shape for the jet, Ωa/c = ξ(Ωz/c)^k, where k < 1 and ξ ∼  1 are dimensionless numbers, one finds that the instability develops only if k < 1/2. Even in this case, the characteristic scale for the development of the instability is very large,

$$\begin{equation} \Omega z/c\,{\sim}\, (10\xi)^{2/(1-2k)}. \end{equation} \tag{ 9 }$$

The above estimates should be considered as preliminary because we assumed that the jet moves as a whole so that there is a frame of reference in which the plasma is at rest. Clearly, the criterion for the CD kink instability can be modified by the effects of relativistic rotation, gradual shear, a surrounding sheath, and sideways expansion. For example, in relativistic Poynting dominated jets the growth rate of the instability strongly depends on the transverse profile of the poloidal magnetic field so that τ goes to infinity when the poloidal field is uniform (Istomin & Pariev 1994, 1996; Lyubarskii 1999). On the other hand, Poynting dominated jets are accelerated such that the acceleration zone spans a large range of scales (Tchekhovskoy et al. 2008; Komissarov et al. 2009; Lyubarsky 2009). Therefore, a small kink growth rate does not necessarily preclude significant growth of helical perturbations. We plan to perform three-dimensional RMHD simulations including some of the missing physical effects such as jet flow and rotation in future work. In particular, this will allow us to investigate the interaction between CD and KH-driven structure in the relativistic regime.

Y.M., K.I.N., and P.H. acknowledge partial support by NSF AST-0506719, AST-0506666, NASA NNG05GK73G, NNX07AJ88G, and NNX08AG83G. This work is supported in part by US-Israeli Binational Science Foundation grant 2006170. The simulations were performed on the Columbia Supercomputer at NAS Division at NASA Ames Research Center and the Altix3700 BX2 at YITP in Kyoto University.

In the simulations presented here, numerical viscosity and resistivity are associated with the finite grid resolution. It is generally expected that numerical viscosity affects the quantitative growth or damping of instabilities, and numerical resistivity allows unexpected magnetic reconnection even when we assume the ideal MHD condition. Therefore, we perform several numerical tests to check whether saturation is adequately captured and investigate how the results in the nonlinear evolution stage depend on numerical resolution.

First, we repeat case A1u (uniform density with constant pitch) for four grid resolutions from 20 to 60 computational zones per simulation length unit L = 8a. Figure 8 shows the time evolution of E_kin,xy and the maximum radial velocity for different grid resolutions. The results depend on grid resolution. Higher grid resolution shows faster growth and an earlier transition time to the nonlinear stage at a larger amplitude. This result is the same as that found in the resolution study of Baty & Keppens (2002). The lowest grid resolution with ΔL = L/20 is clearly inadequate. Although growth does depend on grid resolution, the difference between our choice of ΔL = L/40 and the highest grid resolution of ΔL = L/60 is not significant. Therefore, our simulation results adequately capture the CD kink instability.

Zoom In Zoom Out Reset image size

Second, we check the influence of the transverse boundary. Figure 9 shows the time evolution of E_kin,xy and the maximum radial velocity within a box of length L_z/2 for a simulation box with transverse boundaries at x = y = ±2L (±16a) and a simulation box with transverse boundaries at x = y = ±3L (±24a). The results for the time evolution of E_kin,xy and the maximum radial velocity are the same in both cases. We conclude that behavior of the CD kink instability is not affected by our choice for the transverse boundary.

Zoom In Zoom Out Reset image size

Third, we check the influence of the initial gas pressure. In the force-free model, we ignore the inertia and pressure of the plasma. However, in the simulation the relativistic ideal MHD equations are solved numerically not the force-free equations. Figure 10 shows the time evolution of E_kin,xy and the maximum radial velocity in a box of length L_z/2 for an initial pressure of p₀ = 0.02ρ₀c² and p₀ = 0.005ρ₀c². The results of the time evolution of E_kin,xy and the maximum radial velocity show the same linear growth rate in both cases. However, the lower initial pressure case reaches saturation at a slightly later time with higher amplitude but with a larger decline in the nonlinear stage. Thus, our required nonzero initial gas pressure does affect the saturation level and nonlinear development of the CD kink instability somewhat.

Zoom In Zoom Out Reset image size

Fourth, we check the influence of different adiabatic index. In the nonuniform density case (case n) with Γ = 5/3, the sound speed becomes greater than the relativistic limit of $c/ \sqrt{3}$ at large radial position (see Figure 1(f)). In that region, we should use Γ ∼  eq4/3. However, Γ ∼  eq4/3 is not appropriate in the axial region. We have performed the nonuniform density case (n) with adiabatic indices Γ = 5/3, 13/9, and 4/3. Figure 11 shows results for the time evolution of E_kin,xy and the maximum radial velocity in a box of length L_z/2. The results are almost identical for the different adiabatic indices. Therefore, using a nonrelativistic adiabatic index Γ = 5/3 is appropriate for our simulations.

Zoom In Zoom Out Reset image size