# Linear and nonlinear evolution of current-carrying highly magnetized jets

## Abstract

We investigate the linear and nonlinear evolution of current-carrying jets in a periodic configuration by means of high resolution three-dimensional numerical simulations. The jets under consideration are strongly magnetized with a variable pitch profile and initially in equilibrium under the action of a force-free magnetic field. The growth of current-driven (CDI) and Kelvin-Helmholtz (KHI) instabilities is quantified using three selected cases corresponding to static, Alfvénic and super-Alfvénic jets.

During the early stages, we observe large-scale helical deformations of the jet corresponding to the growth of the initially excited CDI mode. A direct comparison between our simulation results and the analytical growth rates obtained from linear theory reveals good agreement on condition that high-resolution and accurate discretization algorithms are employed.

After the initial linear phase, the jet structure is significantly altered and, while slowly-moving jets show increasing helical deformations, larger velocity shear are violently disrupted on a few Alfvén crossing time leaving a turbulent flow structure. Overall, kinetic and magnetic energies are quickly dissipated into heat and during the saturated regime the jet momentum is redistributed on a larger surface area with most of the jet mass travelling at smaller velocities. The effectiveness of this process is regulated by the onset of KHI instabilities taking place at the jet/ambient interface and can be held responsible for vigorous jet braking and entrainment.

###### keywords:

instabilities - ISM: jets and outflows - stars: jets - MHD - methods: numerical^{1}

^{2}

## 1 Introduction

The investigation of instabilities on the propagation of collimated magnetized flows has been an outstanding research subject during the last 30 years (Cohn, 1983; Ferrari & Trussoni, 1983; Hardee et al., 1992; Todo et al., 1993; Appl, 1996). These instabilities have a substantial importance in characterizing the mechanism of formation and evolution of various observed structures in such flows (Bahcall et al., 1995; Biretta, 1996; Raga & Noriega-Crespo, 1998; Rosado et al., 1999; Nakamura & Meier, 2004; Nakamura et al., 2007; Mignone et al., 2010). Broadly speaking, instabilities in astrophysical jets fall in two classes: external instabilities (caused by the relative motion between jet and external medium) and intrinsic instabilities (related to the toroidal component of magnetic fields).

One of the most important external instabilities is the Kelvin-Helmholtz instability (KHI), responsible for the interaction and mixing between the jet and ambient medium as well as the transfer of linear momentum and jet braking (Bodo et al., 1995; Hardee et al., 1995; Bodo et al., 1998; Rossi et al., 2008). The KHI generally leads to distortions of the interface between the jet and the ambient fluids, eventually producing shocks and turbulent mixing with ambient material. Analyses of the KHI have been extensively accomplished by many researchers either in a linear regime (Turland & Scheuer, 1976; Payne & Cohn, 1985; Hardee et al., 1992; Bodo et al., 1996; Hardee, 2000; Perucho et al., 2004a; Perucho & Lobanov, 2007; Osmanov et al., 2008) or a nonlinear regime (Bodo et al., 1994; Hardee et al., 1997; Bodo et al., 1998; Rosen et al., 1999; Rossi et al., 2004; Perucho et al., 2004b, 2005; Rossi et al., 2008). For example, the linear analysis by Perucho et al. (2004a) showed that the linear growth of KH modes is smaller for faster and colder relativistic jets. Perucho & Lobanov (2007) indicated that the growth rates of KHI modes are significantly reduced by the presence of a thick shear-layer. This important result was confirmed by numerical simulations as well (see, e.g, Perucho et al., 2005, 2007). Some results have been derived for cylindrical relativistic jets that clarified the importance of the nonlinear evolution of KHI in the dichotomy of FRI/FRII extragalactic radio sources (Rossi et al., 2004, 2008). In this respect, the hydrodynamic simulations performed by Rossi et al. (2008) on relativistic light jets have shown that the nonlinear growth of KHI promotes a strong interaction between the jet and the external medium with a consequent mixing and remarkable deceleration.

On the other hand, intrinsic instabilities are generally related to the magnetohydrodynamic (MHD) structure of the jet. In this sense, a key role is played by the relative strength between the poloidal and toroidal (or azimuthal) components of magnetic field. This ratio, referred to as the magnetic pitch parameter, plays an important role in triggering the so-called current-driven instabilities or CDI. According to the mechanisms responsible for jet acceleration and collimation (Blandford & Payne, 1982; Romanova & Lovelace, 1992), the toroidal magnetic field is a significant component that is expected to dominate far from the central engine. Its destabilizing action is at the base of the CDI (Bateman, 1980) and can deeply affect the morphological structure of jets (Ferrari et al., 2011). In this respect, the observations and analytical models presented for the jet in M87 (Sikora et al., 2005; Hardee, 2006, 2011; Walker et al., 2008, 2009) confirm the role of the CDI in conversion Poynting flux to kinetic flux flow within a few hundred gravitational radii. Nevertheless, there it has been suggested over the years that this flux conversion process may be accompanied by other efficient mechanisms such as matter entrainment and jet expansion that can significantly change the nature of jets from magnetically (sub-Alfvénic) to kinetically (super-Alfvénic) dominated (Hardee, 2006, 2011).

Among CDI, the (or kink) mode is the most important one and grows faster than the other modes (Appl et al., 2000; Nakamura & Meier, 2004). In the kink instability, magnetic field lines are compressed on the inner side of a deformed cylindrical flux tube and the magnetic pressure exerted by toroidal component becomes larger than the net magnetic tension. Thus, the jet curvature magnifies so that some of the magnetic energy accumulated by the twisting is released by a kink (Spruit, 1996). Eventually, the hoop stress provided by the toroidal field is reduced (Eichler, 1993) and leads to a whole helical deformation jet (Mizuno et al., 2009, 2011) and/or even jet complete disruption (Nakamura & Meier, 2004).

Linear perturbative analysis of the CDI in MHD jets has been largely performed under different physical conditions by several groups (Cohn, 1983; Appl & Camenzind, 1992; Eichler, 1993; Appl, 1996; Appl et al., 2000; Wanex, 2005; Bonanno & Urpin, 2011). In most of them the force-free limit is considered for helical magnetic field. For instance, Appl et al. (2000) showed that for cold supermagnetosonic jets with dominant azimuthal fields, the CD modes are confined to the jet interior and develop quickly on a time scales of the order of Alfvén crossing time in a frame of reference co-moving with the jet. Wanex (2005) demonstrated that jets with axial sheared flow and sheared magnetic fields reduced the linear growth of CD modes at constant magnetic pitch. In addition, Bonanno & Urpin (2011) considered the case with both azimuthal and axial magnetic fields and showed that the length scale of CD kink modes is strongly sensitive to the magnetic pitch parameter. In the relativistic MHD jets, the linear analysis of CDI by Istomin & Pariev (1994, 1996) demonstrated that cylindrical jets with constant axial (poloidal) magnetic fields are stable against kink CD modes. More recently, the studies by Bodo et al. (2013) on the non-rotating magnetized jets revealed the splitting of the CD kink modes into an inner and outer modes at high magnetization.

On the other side, several numerical studies have concentrated on the nonlinear development of CD instabilities in MHD jets under various physical assumptions (Todo et al., 1993; Lery & Frank, 2000; Lery et al., 2000; Baty & Keppens, 2002; Nakamura & Meier, 2004; Nakamura et al., 2007; Carey & Sovinecy, 2009). In the work by Lery et al. (2000), the non-linear analyses of CDI for cold super fast magneto-sonic jets demonstrates that the CDI play a significant role in re-distributing the current density in the inner parts of the jet. Baty & Keppens (2002) considered the nonlinear interaction between CDI and KHI surface modes on the propagation of supersonic jets showing that CDI prevent the development of KH vortices at the jet surface. According to their results, the magnetic field deformations induced by the nonlinear growth of CDI modes provide a stabilizing factor on the KHI-driven vortical structures leading to a substantial decrease in the mixing process between jets and ambient medium and preventing disruptive effects.

More recently, Mizuno et al. (2009) studied the development of kink instabilities on helically magnetized relativistic static columns showing that, for small pitch values, the growth rate of instability rapidly increases during the linear stages and that the nonlinear evolution features a continually growing helically twisted column. These authors further outstretched this work by considering a sub-Alfvénic velocity shear surface (Mizuno et al., 2011), showing that the temporal evolution of the CDI is dramatically reduced with respect to the static column. Additionally, O’Neill et al. (2012) assumed a similar approximation to that of Mizuno et al. (2009) and studied various local models of co-moving magnetized plasma columns in force-free, rotational, pressure-confined equilibrium configurations and found that the details of initial force balance strongly affect the resulting column morphology.

In this work, we investigate the stability of strongly magnetized jets with initial equilibrium structure described by a force-free helical magnetic field. We adopt a periodic configuration representative of a jet section far from the launching region and consider radially sheared, axial flows with different velocities (Section 2). The configurations are destabilized using an exact eigenfunction corresponding the fastest growing CDI mode of the linearized MHD equations and the evolution is followed through the linear and nonlinear phases using three-dimensional numerical simulations (Section 3). We first assess the impact of grid resolution on the growth of the CDI mode during the linear phase by performing a close comparison with the results from normal mode analysis. Hence we explore the nonlinear behaviour in terms of jet morphology, shear-induced effects triggered by the onset of KHI modes with particular attention to jet braking and momentum transfer from the jet to ambient medium as the system approaches the saturated regime. Contrary to previous studies, our results are relevant to jets with a smaller plasma- () and are based on much larger numerical resolutions ( zones on the jet radius). Finally, our findings are summarized in Section 4.

## 2 Model Setup

### 2.1 Equations and method of solution

In the following, we will investigate the dynamical evolution of current-carrying jets by means of three-dimensional numerical simulations. The simulations are performed by solving the time-dependent ideal MHD equations in Cartesian coordinates

(1) | |||||

(2) | |||||

(3) | |||||

(4) | |||||

(5) |

where and denote the mass density, bulk velocity and magnetic field, respectively. Note that a factor has been absorbed in the definition of . The total (gas + magnetic) pressure is denoted with while the total energy density includes internal, kinetic and magnetic contributions as

(6) |

where the internal energy obeys the perfect gas law with specific heat ratio . Equation (5) sets the evolution of a dynamically passive scalar field that is used to track fluid elements initially residing in the jet ( for and otherwise) or to mark surfaces of constant magnetic flux which, for our case, are initially concentric cylinders ().

Alternatively to equation (3) we also solve, during the initial stages of evolution for , the entropy equation away from shocks

(7) |

where . In a highly magnetized plasma, solving equation (7) instead of equation (3) has the advantage of being more accurate and preventing the occurrence of negative pressure values which could otherwise be triggered by the truncation error of the scheme when retrieving from equation (6).

The numerical computations are carried using the MHD module of the PLUTO code for astrophysical gas dynamics (Mignone et al., 2007). In order to solve equation (1)-(4), we use a third-order algorithm based on a second order Runge-Kutta time stepping, piecewise linear/parabolic reconstruction and the HLL Riemann solver. The employment of a more accurate Riemann solver such the Roe or HLLD Riemann solver (see the discussion in O’Neill et al., 2012) leads, unfortunately, to severe numerical difficulties when dealing with such low-beta plasma configurations.

### 2.2 Initial condition

Our initial condition consists of an infinitely long axisymmetric jet moving in the vertical direction. The jet has an initial equilibrium structure that depends on the cylindrical radius only and characterized by a force-free magnetic field , constant gas pressure and absence of rotations (). In such a way, the toroidal and poloidal magnetic field components obey the time-independent radial component of the momentum equation (equation 2)

(8) |

while density and vertical velocity can be chosen arbitrarily as they do not explicitly appear in the radial balance equation. Here, we set

(9) |

where is the on-axis density, is the jet radius, is the ambient to jet density contrast and is the jet axial velocity. Gas pressure is initially constant and equal to where is the speed of sound.

Following Bodo et al. (2013) we prescribe the azimuthal component of magnetic field to be

(10) |

where is the magnetization radius and determines the maximum field strength. The chosen profile corresponds to a current mainly distributed inside the jet and peaked on the axis. Equation(10) yields a current-free field at large distances as for and a linear profile, , for .

The poloidal component of the field is readily obtained by integrating equation (8) and reads

(11) |

where is the error function while

(12) |

is the magnetic pitch parameter on the axis.

The strength of the magnetic field is controlled by the value of which, without loss of generality, is fixed by the condition that the average Alfvén speed over the jet beam is always unity:

(13) |

Likewise, we take the jet density and radius to be our reference density and length, i.e., and .

Low resolution | High resolution | |||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

Case | ||||||||||||

A0 | 0 | 0 | 0 | |||||||||

A1 | 1 | 10 | 1 | |||||||||

A10 | 10 | 100 | 10 |

Our equilibrium jet model is thus written in terms of the pitch parameter , the jet velocity (in units of the Alfvén speed), the jet density contrast and the sound speed . In the remainder of this paper, we will set and corresponding to highly magnetized jets with and density equal to that of the ambient medium. Also, we take which is close to the lower bound permitted by our equilibrium model, see equation (11).The radial profiles of axial velocity, magnetic field components and magnetic pitch are shown in the three panels of Fig. 1 for the different simulation cases. Note that and (and the pitch) are the same in all cases.

We consider three simulation cases with different values of the axial velocity corresponding, respectively, to the static jet case, a trans-Alfvénic jet and to a super-Alfvénic flow. The three cases will be referred to as A0, A1 and A10 and are listed, together with specific simulation parameters in Table 1. Each simulation case is repeated twice using low and high resolution.

The computational domain is the Cartesian box defined by and , where corresponds to one wavelength of the chosen eigenmode (see §2.3). The total number of zones in the three directions is given by , and , respectively, and is given in Table 1 for low and high resolutions runs. The grid has uniform spacing in the vertical direction and in the region , where (see Table 1). Outside of this region, for , the mesh spacing increases in geometrical progression. In order to obtain cubic cells in the central region of the domain, the number of zones in the and directions is chosen to be while the stretched portions of the grid are discretized symmetrically with and zones on each segment. The number of zones used to resolve the jet radius, , is reported in Table 1 and is, to the extent of our knowledge, the largest one employed so far in simulations of low plasma- in periodic jet configurations. This provides finer resolution in the neighbourhood of the jet where most of the dynamics is expected to take place. Finally, we use periodic boundary conditions in the vertical direction while outflow conditions hold at the remaining boundaries.

### 2.3 Linear theory and choice of perturbation

In order to ease up the comparison with the predictions from global normal mode stability analysis, we perturb the equilibrium configuration using an exact eigenmode of the linearized MHD equations (Appl & Camenzind, 1992). Setting the array of fluid variables, we perturb the initial condition as , where are the equilibrium profiles described in Section 2.2 while

(14) |

is the perturbation. Here denotes the real part, is a complex eigenfunction, is the corresponding (complex) eigenvalue, and are the (real) axial and azimuthal wavenumbers and is a small real number chosen in such a way that the transverse velocity perturbation amplitude equals . The eigenfunction is determined by solving a boundary value problem with appropriate conditions at small and large radii see, for instance, Appl & Camenzind (1992); Appl et al. (2000) or Bodo et al. (2013) for the relativistic treatment. Instability occurs when .

We point out that the choice of an exact eigenmode has revealed to be crucial when comparing the correct growth rate of the desired mode during the linear stages of the numerical computation. This is particularly true when the system can be linearly destabilized by the presence of additional modes in the same range of wavenumbers. This conclusion has been confirmed by a few experiments (not reported here) adopting simpler perturbation forms (e.g. Mizuno et al. (2009) or O’Neill et al. (2012)) and yielding mixed growth rates given by the simultaneous excitation of different modes.

We follow a temporal approach so that and are real numbers while the growth rate of instability is given by the imaginary part of the complex eigenvalue . A plot of the growth rate for as a function of is given in Fig. 2 for different velocities. For (blue curve, first panel), only the CDI mode is present while the appearance of a velocity shear (equation 9) triggers additional KHI modes. When a new KH mode with comparable growth rate appears at larger wavenumbers (red curve in the second panel of Fig. 2) and partially overlaps with the CD mode. The CD and KH growth rates curves intersect through an ’X’ point at and, with increasing jet velocity, the two modes become closer together. Around the lower and upper branches detach from the ’X’ point and the two modes exchange their topological structure. The resulting structure (third panel in Fig. 2, showing the growth rates for ) consists of a mode with larger growth rate and a second mode with smaller amplitude with mixed CD/KH properties. By further increasing the velocity they reach the configuration shown in the right-hand panel of Fig. 2.

In addition we note that reflected KH modes (Bodo et al., 1989) appear as well for . Overall, we expect both KH and CD modes to play a role although it may not be possible to unequivocally isolate their contributions.

We choose the eigenfunction corresponding to the wavenumber at which the growth rate is approximately maximum (vertical dotted line in Fig. 2) and report the precise value in Table 1. Beware that only modes with wavenumbers given by an integer multiple of can actually develop in our computational box.

Finally, since the perturbation in equation (14) is known as a function of the cylindrical radius, we use bilinear interpolation to obtain the values on our discrete Cartesian domain at .

## 3 Results

In the following sections we present our simulation results separately according to their linear and nonlinear evolution. Several diagnostic are computed in a way similar to Mignone et al. (2013) by introducing the horizontal average operator

(15) |

where is any flow quantity, is a weight function and is the position vector. Likewise, we define the volume average operator as

(16) |

When averaging with we use the short-hand notation .

### 3.1 Linear evolution

During the initial stages of the evolution, the perturbation grows and the system slowly departs from equilibrium. As we shall see, this phase is characterized by the dominance of the CDI mode.

#### Comparison with linear theory

As an indicator of the growth of the instability, we use the volume-average of the transverse velocity, defined by

(17) |

where and are the horizontal components of velocity. Equation (17) allows us to perform a direct comparison with the result of linear stability analysis, as shown in Fig. 3 where we plot the time history of the averaged transverse velocity for both low and high resolution computations.

Complementary, we also quantify the growth of the instability by performing a Fourier decomposition of the quantity which, in the limit of small departures from equilibrium, represents essentially the linear radial displacement. In order to compute the discrete Fourier transform in the direction, we consider all the zones lying in a small annular region satisfying and then interpolate the resulting sequence of values on an azimuthal 1D grid using regularly spaced points. The discrete Fourier transform is then computed as

(18) |

where is a fiducial radius, represents the regularly gridded data in the direction and is the azimuthal wave number. The power spectrum is then computed by averaging in the vertical direction: and it is shown in the three panels of Fig 4.

In the static column case most of the main source of energy is magnetic and only the CD mode is present. For (left-hand panel in Fig. 3) the perturbation grows exponentially and both low and high-resolution runs yield growth rates in accordance with linear theory. The maximum is reached around () and () for the low- and high-resolution runs, respectively and the measured growth rate (averaged between ) is and , respectively, to be compared with the exact value reported in Table 1 (). From the left-hand panel in Fig 4, showing the power-spectrum at , the prevalence of the mode is evident.

The evolution of the Alfvénic jet (Case A1) takes place on a time scale similar to the static column case and the duration of the linear phase is approximately the same (). The agreement with linear prediction is excellent and in both cases we measure a growth rate and for () showing convergence for finer mesh spacing (). The power spectrum for this case is shown in the middle panel of Fig. 4 at clearly showing that the kink mode dominates.

In the super-Alfvénic jet, the grid resolution plays a crucial role in determining the numerical value of the growth rate (right-hand panel in Fig. 3). Indeed, from several numerical experiments not reported here, we found this particular case to be the most challenging one owing to the large discretization noise arising from the truncation error of the scheme. The noise triggers grid-sized additional perturbations with sufficiently large amplitudes on top of the initially chosen CDI eigenmode. We note that this spurious effect could be reduced by the combined action of a higher order interpolation scheme (PPM) and by solving for the entropy equation rather than the total energy equation. In addition, the results obtained at higher resolution improve appreciably over the low-resolution computations. A linear phase is observed for and the measured growth rates are and for () at low and high resolutions, respectively, to be compared with the theoretical value of . Similarly to the other two Cases, the dominance of the mode is clear from the right-hand panel of Fig. 4.

#### Three-dimensional structure

The three-dimensional structures of selected jet cases are visible in the panels of Fig. 5 where density and magnetic pressure isosurfaces (top) and magnetic flux surfaces defined by the condition (bottom) are shown at (for case A0 and A1) and (for case A10).

In the static column case (left-hand panel), the initial displacement grows into a twisted helical deformation of the column with density enhancements corresponding to regions of smaller magnetic perturbation and viceversa. This gives rise to a double-stranded spiral structure in which the outward magnetic flux tube stretches and its diameter widen progressively. Owing to the flux-freezing condition, mass and internal energy inside the flux tube are depleted. A similar structure has been also found by other authors (Baty & Keppens, 2002; Mizuno et al., 2009; O’Neill et al., 2012) using different setups. The end of the linear phase is marked by a density build-up on the axis with flow velocities of the order of the Alfvén speed.

In the Alfvénic jet case, the three-dimensional structure at (middle panel in Fig. 5) reveals some similarities with the static column case although we observe the formation of fast magnetosonic disturbances propagating into the ambient medium. At later stages these fronts steepen into weak shock waves. This configuration makes the double helix pattern more difficult to form owing to the presence of the velocity shear. However, similarly to Mizuno et al. (2011), a helical density structure stills persists surrounding the central magnetic helix and propagates along the jet.

For the super-Alfvénic jet (case A10), the linear growth of the perturbation is accompanied by the formation of strong waves propagating obliquely away from the axis and later steepening into shocks (see Section 3.2). The three-dimensional structure, rendered in the right-hand panel of Fig. 5 at , shows a large wavelength non-axisymmetric deformation as well as the presence of small scale surface modes not seen in the previous two cases. This can be inspected from Fig. 6 for . Note also, that during this phase, both low and high resolution computations show similar trends.

### 3.2 Nonlinear evolution

The transition from the linear to the nonlinear phase leads to an overall change of morphology characterized by large scale energy and momentum redistribution. In general, the original equilibrium configuration is destroyed after a few Alfvén crossing time following the end of the linear phase although the different growth of CD or KH modes deeply affects the evolutionary stages of the jets in the three cases considered here.

#### Energetics.

The energy initially carried by the jet is subject to considerable variation as a result of the instability and the ensuing interaction with the external medium. In the top panels of Fig. 6, we plot the relative contributions of the volume-integrated kinetic, thermal and magnetic energies of the jet

(19) |

normalized to the initial total energy of the jet, . Here, is a short-hand notation used to label the different contributions, i.e. , and . Concurrently, to quantify the amount of energy gained or lost by the ambient medium we also plot (bottom panels of Fig. 6) the increment normalized to . Here is defined as in equation (19) with .

Although the three contributions remain approximately constant during the initial evolutionary stages, the end of the linear phase is marked by a rapid transfer of energy from the jet to the ambient medium and, to a lesser degree, by a partial dissipation of magnetic and kinetic energies into heat. In the static column case, a fraction of the initial jet magnetic energy is transferred to the ambient medium thereby accelerating and ultimately heating the surrounding gas. This transition occurs more violently when the jet velocity increases and the beam becomes more kinetically dominated. KH-driven fast magnetosonic disturbances are sheared by the flow and later steepen forming shock waves that provide an efficient dissipation mechanism of mechanical energy. This situation is best depicted in Fig. 7 showing density, magnetic and thermal pressure isosurfaces at different times for the super-Alfvénic jet which becomes forcefully disrupted on a very rapid time scale around . Here the jet loses up to of its initial energy which becomes then available to the ambient medium mostly in the form of heat and, secondly, kinetic energy.

The employment of high resolution introduces modest variations during the onset of the nonlinear phase and results, for sheared jets, in a somewhat more efficient heat generation. Computations performed at higher numerical resolution yield a larger energy gain by the ambient medium and, most remarkably in Case A10, result in an earlier onset of the nonlinear phase. However, its effect is less recognizable at later times.

Note that the total energy is not strictly conserved since the entropy equation is selectively used to evolve the internal energy during the early stages (). However, we have verified that this introduces fluctuations of few percents of the initial value for all simulation cases. In any way, conservation of energy holds until fluid material starts escaping through the lateral boundaries of the computational domain (thin dotted lines in Fig. 6).

#### Morphology

As demonstrated in Section 3.1, the different jet configurations become linearly unstable to the growth of the CDI mode with wavelength equal to the vertical box size. This induces large-scale helical displacements of the jet as already shown in Fig. 5. In order to quantify the amount of jet deformation and distortion we compute the barycenter coordinates as in Mignone et al. (2010, 2013)

(20) |

where is the jet mass density.

In Fig. 8 we plot the radial distance as a function of the vertical coordinate for the selected cases at different simulation times for high resolution (low-resolution computations behave similarly). For the static column (left-hand panel), the initial magnetic surface expands and stretches sideways while maintaining a simple helical structure with constant growing radius up to where . At later times, the helical deformation stretches up but the radius does not remain constant and presents evidence of deformations with smaller vertical wavenumber most likely due to nonlinear wave interactions. A similar behaviour is observed for Case A1 although the departure from a simple constant-radius helix occurs at much earlier times (, middle panel in Fig. 8). Here, for , we see the development of an additional perturbation mode with wavenumber where is the initial perturbation wavenumber (see Section 2.3). This may be explained by the fact that the initial equilibrium for the A1 jet is liable not only to the CDI but also to the KHI mode that has comparable growth rate (second panel of Fig. 2). This mode grows on top of the CDI mode and it reaches a maximum amplitude at while it is later modified by nonlinear interactions. The same effect is also found in the super-fast jet, where large amplitude perturbations with smaller wavelength begin to develop for (see the right-hand panel in Fig. 8).

The final configuration approaching the saturated state is shown in Fig. 9 where we display both the fluid density (red colour) and jet density (green isosurfaces) in the high resolution simulations. At the centre of the domain, a central region with larger magnetic energy and gas density is formed. This region, on the other hand, contains very little of the initial jet material which has been wrapped and twisted inside the helical magnetic flux tubes stretching sideways. These magnetic surfaces enclose a progressively larger volume thereby leading to a substantial dilution of mass and magnetic energy contained therein, see Fig. 9. The same behaviour has been described in the force-free simulation cases presented by O’Neill et al. (2012) with the only difference that, in our case, the nonlinear growth is not solely dominated by the mode.

A different evolution is observed in super-Alfvénic jet (case A10), where the growth of instability goes along with the onset of small-scale surface perturbations developing on top of the CDI mode. While the large-scale column deformation takes place on a few Alfvén time scale, these additional modes are triggered by the velocity shear and grow as small ripples on the jet surface eventually dominating the small-scale structure and leading to a quick violent disruption of the initial cylindrical configuration. This transition is mediated by the onset of KHI acting at the jet/ambient interface and becoming very efficient in promoting entrainment, momentum transfer and mixing. Indeed, by , the jet has completely lost its initial coherent structure and settles into a turbulent chaotic state characterized by a much wider surface area (right-hand panels in Fig. 9).

#### Mass-velocity distribution.

Since the three components of momentum are conserved, it is instructive to compute how the total jet mass is re-distributed as a function of the velocity during the evolution. To this end we partition the velocity value range in small intervals of width and compute the mass-velocity density function as

(21) |

where is a velocity component while when the local zone velocity falls inside the given velocity bin: and otherwise. In other words, the numerator of equation (21) picks out those computational cells having velocity between and where is the width of the velocity bin. Note also that . Using the mass distribution function defined by equation (21), we compute the average and variance as

(22) |

(23) |

respectively.

Fig. 10 shows the distributions of mass as a function of at different times for Cases A0, A1 and A10 in the high resolution simulations. At the distributions are strongly peaked around the initial jet velocity with small dispersion around the mean value meaning that most of the jet mass moves at the same speed. As the system evolves, the net effect of the instabilities taking place inside the jet is that of spreading and re-distributing the initial jet material and momentum over a much wider surface area with consequent reduction of the average jet velocity.

For Case A0, the mass-velocity density function remains symmetric around the initial central value () and spreads progressively forming a sequence of distributions characterized by different modes and dispersion widths (). During the initial linear stages () two secondary lateral peaks corresponding to regions of enhanced vertical velocity perturbation appear (blue line in Fig. 10). These structures are then quickly dissipated leaving an almost-stationary jet configuration (green line). After this phase, we observe, for , the formation of a more chaotic velocity pattern characterized by a scattered distribution with a larger tail (orange line in Fig. 10). Finally, as the system approaches the final state, most of the jet mass has reduced its velocity and the distribution attains a larger peak and smaller dispersion (red line).

For Case A1, the jet gradually slows down as the mass-velocity density function shifts its peak from the initial velocity to smaller values. This transition is particularly effective within the range and gives rise to a sequence of non-symmetric and highly irregular distribution profiles (blue and green lines in Fig. 10). During the final steady state the jet mass is again characterized by a symmetric distribution peaked around with (red line)

A similar shift in velocity can be observed for the super-fast jet (Case A10) although the transition occurs on an even faster time scale. By , in particular, the jet average velocity has already halved and, by , most of the jet mass moves at a much smaller velocity (). As the system approaches the saturated regime, the mass-velocity density function is well described by a symmetric distribution with average value and small dispersion (variance ).

Since the average velocity defined by equation (22) gives basically the volume-integrated jet momentum, it is legitimate to ask what is its relative contribution to the total (i.e. jet+ambient) conserved vertical momentum . In Fig. 11, we plot and as a function of time for Case A1 and A10 in the high resolution simulations where the volume-integrated jet and ambient momenta are, respectively, defined as

(24) |

while . As expected, the end of the linear phase is distinguished by a net transfer of momentum from the jet to the ambient medium. The time scale of this transition closely reflects the changes observed in the mass-velocity distribution profiles and occurs faster in the super-fast jet case. Indeed, for Case A1 and , the jet has lost most of its momentum while for Case A10 the jet has exhausted its momentum already for .

## 4 Summary

In this work we have presented 3D numerical simulations of magnetized jets evolving from an initial cylindrical equilibrium configuration described by uniform density, small magnetic pitch and highly magnetized plasma (). Three cases have been considered corresponding to a static, trans-Alfvénic and a super-Alfvénic jet with Alfvénic Mach number , respectively. Simulations have been performed by perturbing the initial equilibrium state with the exact eigenfunction of linear perturbative theory corresponding to one wavelength of the fastest-growing CDI mode with .

Our results demonstrate that the predicted linear growth rate is well reproduced using a high order reconstruction method and a grid resolution of at least zones per jet radius. Higher resolution may be needed in the case of super-fast jets to reduce grid-induced numerical noise.

Overall, the linear evolution is characterized by large-scale growing helical deformations of the plasma column triggered by the excitation of the CDI mode. The instability breaks the initial axial symmetry and develops on a few tens of Alfvén crossing times while proceeding faster for the super-fast jet. Density and magnetic field tend to form a double-helix pattern featuring regions of alternating enhanced density and magnetic field. The magnetic flux tube stretches sideways and wrap regions of depleted mass and internal energy. As the flow velocity is increased, the jet deformation is accompanied by the propagation of strong fast-magnetosonic waves later steepening into shocks. During this phase, the energy budget remains essentially constant and little energy and momentum exchange are observed.

After the initial transient phase, the equilibrium is considerably altered and the final structure strongly depends on the velocity shear layer which inevitably introduces coupling between CDI and KHI modes. When a velocity-shear is not present, the growth of the CDI mode result in large-scale helical deformations that do not lead to complete disruption. In this sense, our findings favourably compare to the force-free configuration of O’Neill et al. (2012) who accomplished 3D simulations of local comoving jets. However, despite the fact that all magnetized columns under consideration are unstable to CDI modes, the presence of KH surface modes leads to the formation of small-scale distortions that may eventually crumble and destroy the helical structure. This effect becomes more pronounced at larger velocities and, in the most severe case, leads to the complete jet disruption and the formation of a chaotic turbulent flow on a very short time scale. These results are in contrast with the findings of Baty & Keppens (2002) who considered periodic jet configurations threaded by weaker magnetic fields than the ones considered here and found that the presence of CDI modes provides a stabilizing nonlinear interaction mechanism weakening the disruptive effect of KHI perturbations.

The onset of the nonlinear phase is marked by a net dissipation of energy into heat, a process that proceeds gradually for slowly moving jets but becomes violently amplified by the formation of strong magnetosonic shocks for fast jets. Concurrently, as the system evolves towards the saturated regime, the instabilities act so as to spread and re-distribute the initial jet material and momentum over a larger surface area with consequent reduction of the jet average velocity and hence favouring jet braking in few Alfvén time scales. We have verified that, for configurations with non-vanishing initial axial velocity, a large fraction () of the initial jet momentum is transferred to the ambient medium. The re-distribution process is efficiently regulated by the presence of KHI operating at the jet/ambient interface and can thus be held responsible for promoting entrainment, momentum transfer and mixing.

Future extension of this work will enlarge this analysis to global jet simulations in both classical and relativistic regimes.

## Acknowledgements

A.M. wish to thank G. Mamatsashvili for very valuable support on linear perturbative analysis. We acknowledge the CINECA Award no. HP10BCP4GU, 2013 for the availability of high performance computing resources and support.

### Footnotes

- pagerange: Linear and nonlinear evolution of current-carrying highly magnetized jets–References
- pubyear: 2014

### References

- Appl S., & Camenzind M., 1992, A&A, 256, 354.
- Appl S., 1996, A&A, 314, 995.
- Appl S., Lery T., Baty H., 2000, A&A, 355, 818.
- Bahcall J. N., Kirhakos S., Schneider D. P., Davis R. J., Muxlow, T. W. B., Garrington S. T., & Conway R. G., 1995, ApJ, 452, L91.
- Bateman G., 1980, MHD Instabilities, MIT Press, Cambridge.
- Baty H., Keppens R., 2002, ApJ, 580, 800.
- Biretta J. A., 1996, in Hardee P., Bridle A., Zensus J. A., eds, ASP Conf. Ser. Vol. 100, Energy Transport in Radio Galaxies & Quasars. Astron. Soc. Pac., San Francisco, p. 187
- Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883.
- Bodo G., Rosner R., Ferrari A., & Knobloch E., 1989, ApJ, 341, 631
- Bodo G., Massaglia S., Ferrari A., & Trussoni E., 1994, A&A, 283, 655.
- Bodo, G., Massaglia, S., Rossi P., Rosner R., Malagoli A., Ferrari Aet al., 1995, A&A, 303, 281.
- Bodo G., Rosner R., Ferrari A., & Knobloch E., 1996, ApJ, 470, 797.
- Bodo G., Rossi P., Massaglia S., Ferrari A., Malagoli A., Rosner R, 1998, A&A, 333, 1117.
- Bodo G., Mamatsashvili G., Rossi P., & Mignone A. , 2013, MNRAS, 434, 3030.
- Bonanno A., Urpin V., 2011, A&A, 525, A100.
- Carey C. S., & Sovinec C. R., 2009, ApJ, 699, 362.
- Cohn H., 1983, ApJ, 269, 500.
- Eichler D., 1993, ApJ, 419, 111.
- Ferrari A., Trussoni E., 1983, MNRAS, 205, 515.
- Ferrari A., Mignone A., Campigotto M., 2011, in Bonanno A., Kosovichev A., eds, Proc. IAU Symp. 274, Advances in Plasma Astrophysics. Cambridge Univ. Press, Cambridge, p. 410.
- Hardee P. E., 2000, ApJ, 533, 176.
- Hardee P. E., 2006, in Hughes P. A., Bregman J. N., eds, AIP Conf. Proc. Vol.856, Relativistic Jets: The Common Physics of AGN, Microquasars, and Gamma-Ray Bursts. Am. Inst. Phys., New York, p. 57.
- Hardee P. E., 2011, in Romero G. E., Sunyaev R. A., Belloni T., eds, Proc. IAU Symp. 275, Jets at all Scales. Cambridge University Press, Cambridge, p. 41.
- Hardee P. E., Cooper M. A., Norman M. L., Stone J. M., 1992, ApJ, 399, 478.
- Hardee P. E., Clarke D. A., Howell D. A., 1995, ApJ, 441, 644.
- Hardee P. E., Clarke D. A., Rosen A., 1997, ApJ, 485, 533.
- Istomin Y. N., & Pariev V. I., 1994, MNRAS, 267, 629.
- Istomin Y. N., & Pariev V. I., 1996, MNRAS, 281, 1.
- Lery T., Frank A., 2000, ApJ, 533, 897.
- Lery T., Baty H., Appl S., 2000, A&A, 355, 120.
- Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., & Ferrari A., 2007, ApJS, 170, 228.
- Mignone A., Rossi P., Bodo G., Ferrari A., Massaglia S., 2010, MNRAS, 402, 7M.
- Mignone A., Striani E., Tavani M., & Ferrari A. 2013, MNRAS, 436, 1102
- Mizuno Y., Lyubarsky Y., Nishikawa K. I., Hardee P. E., 2009, ApJ, 700, 684.
- Mizuno Y., Hardee P. E., Nishikawa K. I., 2011, ApJ, 734, 19.
- Nakamura M., Meier D. L., 2004,in Bertin G., Farina D., Pozzolieds R., eds, AIP Conf. Proc. Vol.703, Plasma in the Laboratory and in the Universe: New Insights and New Challenges. Am. Inst. Phys., New York, p.308.
- Nakamura M., Li H., Li Sh., 2007, ApJ, 656, 721.
- Osmanov Z., Mignone A., Massaglia S., Bodo G., Ferrari A., 2008, A&A, 490, 493.
- O’Neill S. M., Beckwith K., & Begelman M. C., 2012, MNRAS, 422, 1436.
- Payne D. G., Cohn H., 1985, ApJ, 291, 655.
- Perucho M., Hanasz M., Martí J. M., Sol H., 2004a, A&A, 427, 415.
- Perucho M., Martí J. M., Hanasz M., 2004b, A&A, 427, 431.
- Perucho M., MartÂ´í J. M., Hanasz M., 2005, A&A, 443, 863.
- Perucho M., Lobanov A. P., 2007, A&A, 469, L23.
- Perucho M., Hanasz M., MartÂ´í J. M., & Miralles J. A., 2007, Phys.Rev.E, 75, 6312.
- Raga A.; Noriega-Crespo A., 1998, AJ, 116, 2943.
- Romanova M. M., Lovelace R. V. E., 1992, A&A, 262, 26.
- Rosado M., Raga A. C., & Arias L., 1999, AJ, 117, 462.
- Rosen A., Hardee P. E., Clarke D. A., Johnson A., 1999, ApJ, 510, 136.
- Rossi P., Bodo G., Massaglia S., Ferrari A., & Mignone A., 2004, Ap & SS, 293, 149.
- Rossi P., Mignone A., Bodo G., Massaglia S., & Ferrari A., 2008, A&A, 488, 795.
- Sikora M., Begelman M., Msejski G., & Lasota J.-P., 2005, ApJ, 625, 72.
- Spruit H. C., 1996, Astrophysical Journal, 2022S.
- Todo Y., Uchida Y., Sato T., & Rosner R., 1993, ApJ, 403, 164.
- Turland B. D., Scheuer P. A. G., 1976, MNRAS, 176, 421.
- Walker R. C., Ly C., Junor W., Hardee P. J., 2008, J. Phys.: Conf. Ser., 131a, 2053, 2053.
- Walker R. C., Ly C., Junor W., Hardee P. J., 2009, in Hagiwara Y., Fomalont E., Tsuboi M., Murata Y., eds, ASP Conf. Ser, Vol. 402, Approaching Micro-Arcsecond Resolution with VSOP-2: Astrophysics and Technologies. Astron. Soc. Pac., San Francisco, p. 227 227.
- Wanex Lucas F., 2005, Ap&SS., 298, 337.