The Kinked Jet of the Crab Nebula

Modelling the Kinked Jet of the Crab Nebula


We investigate the dynamical propagation of the South-East jet from the Crab pulsar interacting with supernova ejecta by means of three-dimensional relativistic MHD numerical simulations with the PLUTO code. The initial jet structure is set up from the inner regions of the Crab Nebula. We study the evolution of hot, relativistic hollow outflows initially carrying a purely azimuthal magnetic field. Our jet models are characterized by different choices of the outflow magnetization ( parameter) and the bulk Lorentz factor ().

We show that the jet is heavily affected by the growth of current-driven kink instabilities causing considerable deflection throughout its propagation length. This behavior is partially stabilized by the combined action of larger flow velocities and/or reduced magnetic field strengths. We find that our best jet models are characterized by relatively large values of () and small values of . Our results are in good agreement with the recent X-ray (Chandra) data of the Crab Nebula South-East jet indicating that the jet changes direction of propagation on a time scale of the order of few years. The 3D models presented here may have important implications in the investigation of particle acceleration in relativistic outflows.

MHD - pulsars: individual: Crab Nebula - ISM: jets and outflows - instabilities - shock waves

1 Introduction

Pulsars loose their rotational energy through a relativistic wind of waves and particles. The interaction of these outflows with the surrounding ambient produces Pulsar Wind Nebulae (PWNe), observable from radio to -rays. Pulsar Wind Nebulae often show a torus-jet structure (see, e.g., Kargaltsev & Pavlov (2008) for a review, and the Chandra images of the Crab, Vela and B1509-58 Nebulae). Several theoretical (Lyubarsky & Kirk, 2001; Lyubarsky, 2002; Pétri & Lyubarsky, 2007) and numerical (Komissarov & Lyubarsky, 2003; Del Zanna, Amato, & Bucciantini, 2004; Del Zanna et al., 2006) studies attempted to explain and reproduce this structure.

The Crab Nebula is surely the most popular and studied PWN. It’s powered by a pulsar with a very large spin-down luminosity, erg/s, that is carried away by a relativistic and highly magnetized wind. According to estimates of the number of pairs emitted by the Crab Pulsar, the energy flux carried by the wind is dominated by the Poynting flux while close to the star, and dominated by the particles when close to the termination shock (Kennel Coroniti 1984). First theoretical (Kennel & Coroniti, 1984) and numerical (Komissarov & Lyubarsky, 2003; Del Zanna, Amato, & Bucciantini, 2004) studies suggested that the magnetization parameter, defined as the ratio between the Poynting flux and the kinetic energy of the particles, should lie in the range at the termination shock.

Low values are required to avoid the excessive axial compression of high 1- and 2-D models, that push the pulsar wind termination shock too close to the pulsar (Komissarov, 2013; Lyubarsky, 2012). Nevertheless, as recently pointed out by Porth et al. (2013), 3D high- models of PWN have the same morphology of two-dimensional axisymmetric low- models owing to the presence of kink instabilities (Begelman, 1998) that reduce the axial compression and lead to uniform pressure within the Nebula. In addition, a lateral dependence of the wind magnetization, that increases towards the axis yielding close to the poles, was proposed by recent investigations (Lyubarsky, 2012; Komissarov, 2013). In these models, the Nebula and the jet would be therefore injected with a highly magnetized plasma, and they would be regions of strong magnetic dissipation (Komissarov, 2013; Porth et al., 2013).

Unexpectedly, the Crab Nebula produces strong and day long -ray flares (Tavani et al., 2011; Abdo et al., 2011; Striani et al., 2011; Buehler, 2012; Striani et al., 2013). The average Nebula spectrum, that shows a cutoff around MeV, is interpreted as synchrotron emission of electrons and positrons accelerated at the termination shock. The -ray flares show, instead, an energy peak at MeV. Magnetic fields of mG3 and particles’ Lorentz factors (Striani et al., 2013) are required to match the observable data for the flares. A very fast and efficient acceleration process is required to explain this very short time-scale and the hard energy peak (Tavani et al., 2011).

Figure 1: Summed image for 7 Chandra ACIS observations of the Crab Nebula’s southern jet. The data were taken between 2010, September and 2011, April. The color bar gives the summed counts per ACIS pixel over a total effective exposure of about 4.2 ksec. One ACIS pixel is 0.492 seconds of arc. North is up, East is to the left, and the pulsar is at (0, 0) in the figure. These data are from an exhaustive study of all of the Chandra observations (Weisskopf et al. in preparation).

The Chandra images of the Crab Nebula show a strongly bent jet (Weisskopf et al., 2000), see Fig. 1. Interestingly, the comparison of the 2001 and 2010 Chandra images show that the deflection changed its orientation in this lapse of time (Weisskopf, 2011)4. This phenomenon can be the consequence of kink instabilities taking place in the jet. Previous MHD simulations of the Crab Nebula (see, for instance Komissarov & Lyubarsky, 2003, 2004; Del Zanna, Amato, & Bucciantini, 2004; Del Zanna et al., 2006) assumed a 2-D axisymmetric model, that does not allow deformations of the jet induced by current-driven (CD) modes with azimuthal wave number .

In this paper, for the first time, we present 3-D relativistic MHD simulations of the Crab Nebula jet and investigate for which values of the magnetization and Lorentz factor our simulations reproduce the observed jet behavior. In this context, previous numerical simulations by Mignone et al. (2010) have demonstrated that relativistic jets possessing an axial current can be significantly affected by the growth of non-axial symmetric perturbations leading to large deflection off the main longitudinal axis. Similar results have been confirmed by other investigators, e.g., Moll et al. (2008) (for Newtonian MHD) and Porth (2013) (for relativistic MHD) who have studied the growth of current driven modes in magnetocentrifugally driven jets through three dimensional simulations. Likewise, Mizuno et al. (2011) have demonstrated the importance of three-dimensional dynamics on the stability of a hot plasma column threaded by a toroidal magnetic field.

Besides, kink instabilities in the jet could trigger magnetic reconnection episodes that could eventually be responsible for the observed -ray flares. Magnetic reconnection in current sheets created in the torus and/or the jet is, indeed, considered as one of the possible mechanism for the acceleration of particles at PeV energies (e.g., Cerutti et al. (2012); Uzdensky et al. (2011); Cerutti et al. (2013)).

The paper is organized as follows: in section 2 we describe the initial model setup and its connection to the 2D axisymmetric models of (Del Zanna, Amato, & Bucciantini (2004), dZAB04 henceforth) and describe the relevant simulation parameters. In section 3 we present and discuss our results while conclusions are drawn in section 4.

2 Model Setup

2.1 Equations and Method of Solution

In what follows we consider an ideal relativistic magnetized fluid with rest-mass density , bulk velocity , magnetic field and thermal (gas) pressure . Numerical simulations are carried out by solving the equations of special relativistic MHD (RMHD) in conservation form (Anile, 1990):

where is the Lorentz factor, is the momentum density, denotes the electric field and is the gas enthalpy which relates to and via the ideal gas law:


We adopt here appropriate for a hot relativistic plasma. Total pressure and energy include thermal and magnetic contributions and can be written as


Note that velocities are naturally expressed in units of the speed of light .

An additional equation, describing the advection of a passive scalar (or tracer) , is included to discriminate between jet material (where ) and the environment (where ):


The system of Equations (2.1) is solved in the laboratory frame using the PLUTO code for astrophysical gasdynamics (Mignone et al., 2007) with linear reconstruction and a second-order TVD Runge-Kutta time stepping with a Courant number of . PLUTO is a shock-capturing code targeting supersonic fluid in presence of discontinuous waves and largely based on the usage of Riemann solvers. For the present purpose we employ the HLLD Riemann solver of Mignone, Ugliano, & Bodo (2009) and revert to the HLL solver in presence of strong shocks following the hybrid approach described in the appendix of Mignone et al. (2012). The divergence-free condition of the magnetic field is accurately treated using the constrained transport method.

2.2 Initial and Boundary Conditions

Figure 2: Two-dimensional schematic representation of the initial condition. The jet enters from the nozzle at the bottom boundary into a hot cavity region confined by the supernova remnants (SNR) expanding into the outer interstellar medium (ISM). Tick labels are given in units of the jet radius .

Our initial conditions draw upon the same configuration used by dZAB04 where a freely expanding supernova remnant initially fills the region where is the spherical radius and , see Fig 2. The supernova ejecta is unmagnetized with total mass given by


and a radially increasing velocity profile where is fixed by the condition


Here we take and in accordance with Eq (5) and (6) of Del Zanna, Amato, & Bucciantini (2004). Further out, for , the fluid is uniform and static with density and pressure values representative of the interstellar medium (ISM), i.e., and , where is the proton mass.

Differently from dZAB04, however, we do not include the pulsar wind nebula for but consider, instead, a static hot plasma region where jet acceleration is assumed to take place. Density and pressure in this region are denoted with and while the magnetic field is absent, see also Fig 2.

As the acceleration mechanism cannot be consistently described within our model, we assume that the jet has already formed as the result of the magnetic hoop stress and collimation processes taking place around the polar axis. For this reason, the jet is modeled as a continuous injection of mass, momentum, magnetic field and energy from the lower -boundary inside the circular nozzle where is the cylindrical radius and is the jet radius. Here inflow values are prescribed in terms of constant density and axial velocity


where is the bulk Lorentz factor.

Since present axisymmetric models of PWN adopt a purely toroidal field we initialize the magnetic field to be azimuthal with the following radial profile (Komissarov, 1999; Mignone, Ugliano, & Bodo, 2009):

where encloses a cylinder carrying a constant current and sets the magnetic field strength. We take and fix the value of from the jet magnetization parameter :


where is the average magnetic energy inside the beam.

The gas pressure is recovered from the radial momentum balance across the jet which, in absence of rotations, takes the form


where, for simplicity, a constant value of is used. The previous equation has solution

where is the ambient pressure determined by flow Mach number . Note that in this configuration the pressure is maximum on the axis and monotonically decreases until where it matches the ambient pressure. The maximum value is obtained from Eq. (2.2) with and it takes the value . Consequently, highly magnetized jets also possess larger internal energies, see Table 1.

As axisymmetric models of PWN predict hot and hollow jets, we prescribe the jet mass density to be with and set the sonic Mach number . In such a way the initial density contrast between the supernova remnant and the jet is . Conversely, we leave and as free parameters and perform computations with two different values of and three different magnetizations for a total of six cases, as shown in Table 1. This choice of parameters is consistently based on the results obtained from the 2D axisymmetric simulations of dZAB04 (also repeated by our group with good agreement) from which comparable values of density, pressure and magnetic fields could be inferred.

In the injection nozzle, we perturb the transverse velocities by introducing pinching, helical and fluting modes with corresponding wave numbers (Rossi et al., 2008):


where are randomly chosen phase shifts while high () and low-frequency () modes are given by and . The amplitude of the perturbation is chosen in such a way that the fractional change in the Lorentz factor is :


The computational domain is defined by the Cartesian box and with and covered by computational cells. The grid resolution is uniform in the direction and inside the region where zones are used. The grid spacing increases geometrically outside this region up to the lateral sides of the domain.

We employ outflow (i.e. zero-gradient) boundary conditions on the and sides of the computational box as well as on the top boundary. In the ghost zones at the bottom boundary and outside the injection nozzle, we set , and to be anti-symmetric with respect to the plane while the remaining quantities are symmetrized. Fluid variables inside and outside the nozzle are then smoothly joined using a profile function:


where is a fluid variable value inside the nozzle, is the symmetric (or anti-symmetric) value with respect to the plane. Finally, we use for density, velocity, pressure and vertical component of the field while is used for and .

A1 2 0.1
A2 2 1
A3 2 10
B1 4 0.1
B2 4 1
B3 4 10
Table 1: Simulation cases describing our initial jet configuration. While and are the primary parameters used in our model we also give the derived values of magnetic field strength (, 4-th column), ambient pressure ( 5-th column), on-axis pressure (, 6-th column) and the ratio between average thermal and magnetic pressures (plasma , last column).

3 Results

Figure 3: Three-dimensional contour surfaces of the gas pressure (blue) and parameter (orange) for the six jet configurations. Tick labels are given in units of the jet radius, cm. Low-speed jets with (cases A1, A2 and A3) are shown, from left to right, in the top panel, while high-speed jets with (cases B1, B2 and B3) are shown in the bottom panel. Snapshots are taken at different times when the jet has approximately reached the end of the computational domain.

The six simulation cases introduced in Table 1 show remarkable differences in several aspects such as the propagation velocity, large-scale morphology, interaction and mixing with the environment. These will be discussed in the following.

As the computations produced a significant amount of data ( TB) an efficient post-processing analysis is crucial in order to reduce and extract relevant quantitative results. Our experience has shown that several morphological and dynamical aspects can be readily interpreted by means of horizontally-averaged quantities defined as


where integration is performed over horizontal planes at constant-, can be any flow quantity, is the position vector and is a weight (or filter) function used to include or exclude certain regions of the flow according to specific criteria. We select, for instance, computational zones containing more than of the jet material and moving at least at of the speed of light using

where is the passive scalar obeying equation (3).

3.1 Overall Features.

During the very early phases of evolution, the jet propagates almost undisturbed until it impacts the backward dense layers of the supernova remnant. From this time on (typically year), we observe a drastic deceleration as the jet pushes against the much heavier material of the remnant.

After a few tens of years, the typical structure consists of a large over-pressurized turbulent cocoon enclosing a collimated magnetized central spine moving at mildly relativistic velocities. This is illustrated in Fig. 3 where a volume rendering of thermal pressure and magnetization parameter is shown for each of the six cases at the end of the simulation5.

The cocoon appears to be weakly magnetized and the field remains mainly concentrated in the beam (a similar structure is also observed in the 2D simulations of dZBL04) preserving the initial toroidal structure. Along the jet spine the flow experiences a series of acceleration and deceleration phases owing to the presence of conical recollimation shock waves (or working surfaces) corresponding to regions of jet pinching.

The large-scale morphology is characterized by elongated curved structures which, depending on the case, may considerably departs from axial symmetry. The amount of bending and twisting varies according to the combination the flow Lorentz factor and magnetization () and it is described in more detail in Section 3.3.

Jet Position.

Figure 4: Jet head position as function of time for the six simulation cases A1, A2, A3 (black, green and red solid lines) and B1, B2, B3 (corresponding dashed lines). The thin dotted line on the left represents the position of the outer supernova shock.

Fig. 4 shows the jet head position as a function of time, measured as the maximum height at which jet material has propagated.

Low-speed jets advance slowly () owing to the large density contrast and evolve entirely within the remnant confined by the outer SN shock (see Fig 2). For increasing magnetization the propagation is driven by the additional magnetic pressure support while the mechanism of instability tends to saturate.

Conversely, jets with larger (cases B1, B2 and B3) advance more rapidly () and cross the outer SN shock (dotted line) at earlier times ( years) where they suddenly accelerate because of the reduced density contrast. In particular, owing to its low magnetization and relatively large kinetic energy, the B1 jet is very little affected by the growth of CD modes and its trajectory remains essentially parallel to the axis. With increasing magnetization, the jets in case B2 and B3 are slowed down by appreciable bendings of the flow direction and show comparable propagation speeds.

3.2 Jet Internal Structure.

Figure 5: Volume rendering of the magnitude of the current density for the A3 jet, at showing the formation of helical structures in the front-end regions.

The jet structure does not remain homogeneous during its propagation but, rather, shows substantial variations of several fluid quantities all along its length. Broadly speaking, we are able to identify two regions with different properties. In the back-end region, close to the injection nozzle, the jet has reached a quasi-stationary structure with a number of standing conical shocks. In the front-end region, the dynamics is characterized by a rapid variability and strong interaction between the jet and the remnant. In these regions and for large magnetizations, the beam takes the shape of a twisted helical structure, see Fig. 5 showing a three-dimensional view of the current density for the A3 jet after years.


Figure 6: Averages profiles of the electromagnetic (top) and kinetic (middle) energies normalized to their initial value at as functions of the vertical distance at different times (reported in the legend) for the six simulation cases. The bottom panel shows the maximum Lorentz factor. Regions of strong compression are evident by the quasi-periodic oscillations. The Lorentz factor grows immediately upstream of the shocked flow where magnetic and kinetic energies are smaller and drops discontinuously in the post-shock regions.

A common feature that can be identified all along the jet is the presence of pinching regions corresponding to the formation of magnetized shock waves. These can be distinguished, for instance, by looking at the horizontally-averaged electromagnetic and matter kinetic energies



where the weight function selects only material that is mainly composed by jet particles, see Eq. (3).

In Fig. 6 we plot , and the maximum Lorentz factor (taken on planes) as functions of just before the jet has exited the computational domain or encountered the outer supernova shock. Average magnetic and kinetic energies exhibit quasi-periodic oscillations along the beam due to jet pinching with the corresponding formation of internal shocks with large compression factors. These cycles are more evident in the slowly-moving jets that reveal shocks with larger strengths. Here the frequency of oscillations increases with and magnetic fields tend to dissipate more rapidly.

Indeed, as discussed in Mignone et al. (2010), the presence of a dominant azimuthal magnetic field component prevents the inner jet core from interacting with the surrounding thus substantially reducing the loss and transfer of momentum. The net effect of this shielding mechanism is to sustain the kinetic energy at the expenses of magnetic energy thus leading to a significant decrease of along the beam.

Figure 7: Two-dimensional color distribution map of the plane-averaged parameter normalized to the initial injection value as function of time (abscissa, in years) and vertical height (ordinate, in light-years) for the six cases. Note the substantial decrease of in the more magnetized cases (A3 and B3).

This is best illustrated in Fig. 7, where we show a 2D color distribution map of the horizontally-averaged parameter normalized to its initial value.


As the jet advances into the remnant, the propagation is accompanied by the formation of highly intermittent unstable structures during which jet fragmentation is frequently observed. These events take place on a short time-scale (typically less than a year) and in correspondence of large kinked deflection where the jet beam temporarily breaks down forming strong intermediate shock waves resembling the main termination shock. A typical example is illustrated in Fig. 8 where we show, from left to right, a volume rendering of the parameter, the current density and the Lorentz factor for the A2 jet during a short temporal evolution.

We point out that the features produced in our jet fragmentation (hot-spots) do not have an equivalent observed optical or X-ray signature, although some hints of such structures are present in the X-ray maps of the jet terminal part (see left panel in Fig. 14). On the other hand, the locations of our jets that correspond to higher values of, e.g, the parameter, might have a relation with the synchrotron processes responsible for the optical and X-ray emission. This issue will be further investigated in forthcoming studies.

Figure 8: Jet fragmentation for the A2 case illustrated through a sequence of close-by frames at (left column), (central column) and (right column) years. From top to bottom the figure shows, respectively, the volume rendering of the parameter, thermal pressure, current density magnitude and bulk flow Lorentz factor.

Magnetic Field Topology.

Figure 9: Average angle between the magnetic field vector and the mean flow direction as a function of for cases A1, A2 and A3 (solid lines). The dotted lines show the (approximate) position of the terminal reverse shock for the three cases. The evolution times corresponds to those indicated in Fig. 3
Figure 10: Magnetic field lines and pressure volume map for the A2 jet at years.

In order to gain some insight on the topology of magnetic field, we first compute the average flow direction by integrating, for each , the velocity vector on horizontal planes:


with defined by Eq. (3). We then decompose the magnetic field into components that are parallel and perpendicular to the direction given by :


where is the magnetic field component parallel to the horizontally-averaged velocity vector and is the component of the field perpendicular to it. The average cosine between magnetic field and mean flow direction is then simply obtained from:


Fig. 9 shows that the magnetic field in the jet remains essentially perpendicular to the (average) flow trajectory for most of the jet length while it acquires a poloidal component immediately after the terminal reverse shock. This is confirmed by the direct three-dimensional visualization of the magnetic field lines in proximity of the beam (shown in Fig. 10 for the A2 jet) which reveals that the field retains the initial toroidal shape that progressively evolves into a bent helical structure with a small pitch. Indeed, as mentioned in section 3.2, the magnetic field acts as an effective screening sheath that hampers the development of small scale perturbations at the interface between the jet and the surrounding.

3.3 Jet Deflections.

Figure 11: Maximum radial deflection as a function of time for the different cases discusses in the text. Solid lines refer to case A1 (black), case A2 (green) and A3 (red) while dashed lines refer to case B1, B2 and B3, respectively.

A remarkable feature observed in several runs is a curved trajectory featuring large time-dependent deflection of the jet beam away from the main longitudinal axis (see Fig. 3). As discussed by Mignone et al. (2010), this behavior may be ascribed to the onset of CD instabilities triggered by the presence of a toroidal magnetic field component (see also Moll et al., 2008; Porth, 2013). This result is confirmed by numerical investigations of infinitely-long periodic jets adopting the same initial structure (Mignone et al., in preparation) revealing the presence of CD instabilities with a rapid growth of the (or kink) mode.

We introduce a measure of the deflection radius and the deflection angle by computing at each time and vertical height the centroids and :


and choosing the weight function as in Eq. (3).

Fig. 11 plots the maximum of taken over the spatial coordinate as a function of time. Case A2 and A3 show the largest bendings reaching values in excess of jet radii. Although to a less degree, cases B2 and B3 are also prone to appreciable wiggling () whereas case B1 propagates almost parallel to the main longitudinal axis with very weak bending of the beam.

A more detailed investigation is given in Fig. 12 and Fig. 13 showing 2D colored distribution maps of the deflection radius and angle as functions of time and vertical coordinate. Jets become increasingly more flexed as they propagate from the injection region up to the head where they attain the largest deformations (Porth, 2013). The shape of these deformed structures roughly resembles the morphological features that can be inferred from observational data. This is shown in Fig. 14 where we present a qualitative comparison between the X-ray Chandra observation of the terminal part of the jet in 2010 (top panel) and our simulated jet (case A2, bottom panel).

Figure 12: Two-dimensional colored distribution maps of the average deflection radius (in units of the jet radius) as function of time (abscissa, in years) and vertical height (ordinate, in light-years) for the six simulated cases. The amount of deflection at a given time and height is quantified by a different color given by the varied shades of blue (), green (), orange () and red ().
Figure 13: Two-dimensional colored distribution maps (similar to Fig. 12) of the average deflection angle as function of time (abscissa, in years) and vertical height (ordinate, in light-years) for the six simulated cases.
Figure 14: An example of qualitative agreement between data (top showing Chandra X-ray map of the terminal jet part) and simulations (bottom, showing a volume rendering of for the A2 jet).

For moderate and strong magnetizations in low-velocity jets (A2 and A3), the deflection keeps growing indefinitely during the course of the simulation. In these cases we observe, for , a change in the propagation direction on a time scale of years while, afterwards, jet material keeps flowing along a particular curved direction established by the large-scale backflow circulation pattern and changes direction very slightly. Here the inertia of the flexed jet becomes so effectively large that any restoring mechanism is unable to push the beam back along the main longitudinal axis. In the corresponding higher velocity cases (B2 and B3) the amount of bending is reduced owing to the increased inertia which acts as a stabilizing factor. In these cases the jet changes its propagation angle by a large amount ( or more) around the main longitudinal axis over a time period years (see the right panels in Fig. 13).

We point out, however, that previous simulations (not shown here) have demonstrated that the evolution of the trajectory and the corresponding amount of deflection may be noticeably affected by the shape of the initial perturbation (given by Eq. 9).

We thus conclude that the magnetization parameter plays a crucial role in destabilizing the jet: weakly magnetized configurations (A1 and B1) are less affected by the growth of instability than the corresponding moderately (A2 and B2) and highly (A3 and B3) magnetized models which show comparable growth rates. The effect of the Lorentz factor, on the other hand, is that of reducing the growth of instability, in accordance with the results obtained from the linear stability analysis of Bodo et al. (2013).

Flow Direction.

Figure 15: Same as Fig. 12 but for the average propagation angle (in degrees) in the positive direction. The color map spans from (white) to (dark red).
Figure 16: Same as Fig. 15 but for the negative direction. In this case, ranges from (white) to (dark red).

The change in the jet trajectory is associated with a corresponding variation of the average propagation velocity. We compute the average angle between the velocity vector and the axial direction by integrating, on constant -planes, the projected velocity:


where and are used to discriminate between jet material moving in the positive or negative vertical direction, respectively. In other words, gives a measure of the forward flow while is associated with the backflow motion. The filter function is defined as usual (Eq. 3).

Fig. 15 shows a colored distribution map of as a function of time and vertical distance. In all cases, sudden changes in the trajectory occur at the jet head where magnetic field are amplified and the flow is abruptly decelerated through the termination shock. However, low-speed jets tend to assume a large-scale curved structured since the average propagation angle gradually changes from (straight propagation) close to the launching region up to at the jet head. Conversely, high speed jets are stabilized by the larger Lorentz factor and propagate more parallel to the longitudinal axis building large kicks mainly in proximity of the jet head. Here the velocity is drastically reduced and the magnetic field is increased thus de-stabilizing the motion.

Backflow Motion.

The departure from axial symmetry produced by kinked deformations has the side effect of promoting a predominantly one-sided backflow motions along the negative direction. This is shown in Fig. 16 where we plot a 2D color map of (computed from Eq. 19) as function of time and vertical height. Backflows are strongest in case A1 and progressively lessen at larger magnetizations or for larger Lorentz factor cases (e.g., they are almost absent in the B1 jet). A three-dimensional view is given in Fig. 17 for the small Lorentz factor jets. In the most prominent cases, the backflow is able to survive over a distance of and for several years before its kinetic energy is dissipated in the form of turbulent motion during the interaction with the ambient medium.

Figure 17: Three-dimensional rendering of the passive scalar distribution for the A1, A2 and A3 jets at the time reported in the legend. The size of the backflow region reduces with increasing magnetization. Regions in red mark fluid elements composed by at least of the jet material () whereas blue regions mark fluids elements that have mixed with remnant and are composed by less than of the jet material ().

3.4 Dissipation.

A crucial aspect in the modeling of magnetically driven outflows is the dissipation of magnetic fields at both large and small scales. This is discussed in the following sections.

Large-scale Dissipation

Energy dissipation at large scale may be quantified by comparing the electromagnetic to thermal energy ratios inside the jet with that of the turbulent cocoon. We distinguish the two regions by taking advantage of the passive scalar and choosing different weight functions. Inside the jet we compute


with defined by Eq. (3). This choice ensures that only high-velocity regions containing more than of the jet material are included. On the other hand, in the environment we define


using when , which includes moving material that has already mixed. The ratio is plotted for the jet (solid lines) and the environment (dash-dot lines) in Fig. 18 at the end of each simulation case, just before reaching the end of the computational domain. Inside the jet, the ratio between magnetic and thermal energies present gradually decreasing quasi-periodic oscillations (see §3.2) that drop sharply at the jet termination shock. Conversely, the energy distribution proportions inside the cocoon are more homogeneous and settle down to approximately the same values () independently of the value of jet . This demonstrates that both Poynting and kinetic energy fluxes are efficiently diverted at the termination shock and thereafter scattered and dispersed via the backflow to feed the cocoon. During this process the field becomes significantly randomized and the energy distributions reach a state far from equipartition, independently of the initial magnetization. From this perspective, jets appear to be a very efficient way to dissipate magnetic energy through the interaction of the head of the jet with the surrounding remnant.

Figure 18: Ratio between average electromagnetic and thermal energies inside the jet (solid lines) and outside in the cocoon (dot-dash). Plots in the top and bottom panels refers to the slower and faster jet cases, respectively.

Small-scale Dissipation

Figure 19: Maximum current density taken over planes as a function of for the A2 jet (top) and A3 jet (bottom). The solid, dashed and dotted lines mark, respectively, three different close-by simulation times reported in the corresponding legends. The largest peaks are shown by the dashed lines.

Although our simulations do not include dissipation mechanisms other than numerical that acts at the cell size, it is instructive to localize strong current sheets that may host regions where particle acceleration is likely to take place. In Fig. 19 we plot the maximum value of the current density, as a function of the vertical height at three successive simulation times for the A2 and A3 jets. Current peaks take place in regions of strong pinching where the flow is shocked. According to the broad distinction of back-end and front-end regions given in Section 3.2, these structures may have considerably different lifetimes. While the inner regions of the jet have reached a quasi-steady periodic structure with little time-variability, the outer regions reveals the formation of short-lived current peaks that diffuse on a time scale which is of the same order or less than our temporal resolution, approximately months. The 3D spatial distribution of the current density, corresponding to the dashed line in the bottom panel of Fig. 19, is shown at yrs in Fig. 20 for the A3 jet. Note that the formation of the strong current peak at occurs concurrently with the development of a jet kink and the abrupt change of flow direction. The sudden rise of these localized current peaks favors the formation of reconnection layers where induced electric fields may lead to efficient particle acceleration. This possibility will be addressed in future works.

Figure 20: Three-dimensional contour plot of the current density after years for the A3 jet. Six surfaces of constant current density are displayed using different colors, see the attached legend. The current peak at corresponds to the one shown by the dashed line in the bottom panel of Fig. 19.

4 Summary

In this work, we have presented numerical simulations of three-dimensional relativistic magnetized jets propagating into a heavier supernova remnant. The proposed jet models aim at investigating the dynamics and morphology of the bipolar outflows observed in the Crab Nebula. Our initial configuration stems from the results of previous 2D axisymmetric numerical models of Pulsar wind nebulae (Del Zanna, Amato, & Bucciantini, 2004) that predict the formation of hot under-dense jets as the result of the magnetic hoop stress collimation process. Using these models to constrain our input parameters, we have performed numerical simulations of jets initially carrying a purely azimuthal field by varying the bulk flow Lorentz factor and the ratio between Poynting and kinetic energy fluxes () at the injection region.

Our results show that jets with moderately/high magnetic fields () are prone to large scale non-axisymmetric current-driven instabilities leading to prominent deflections of the jet beam away from the axis. This effect is enhanced by the high density contrast () between the jet and the supernova remnant which leads to the formation of a strongly perturbed beam and largely over-pressurized cocoons. The typical timescale for the formation of these curved patterns is, for the representative parameters adopted in our model, of the order of a few years and thus compatible with observations. While the kink mode is mainly responsible for the formation of non axisymmetric morphologies, the presence of axial pinch modes engenders a knotty structure with a chain of strong intermediate shocks with large compression factors. High variability is observed in the outer regions where the dynamics is strongly influenced by the interaction of the jet material with the ambient medium. Here rapid variations of the flow properties are characterized by the formation of intermediate magnetized shocks in proximity of sudden kinked deflection of the flow trajectory. These are short-lived episodes leading to the formation of intermittent unstable structures such as sporadic jet fragmentation or strong reconnection layers on a time-scale of a few months.

Our computations demonstrate that the development of these unstable patterns is more pronounced in relatively low- jets () which show the largest deflections and evolve entirely inside the remnant. Conversely, jets with larger Lorentz factor () propagate with larger inertia, are less affected by the growth of pinch or kink modes and drill out of the remnant in less than years. Furthermore, jets with low magnetic fields () are weakly affected by the onset of current-driven modes and tend to propagate more parallel to the longitudinal axis.

Based on our relativistic 3D MHD simulations, for the first time we can conclude that moderately to high- jets with relatively small Lorentz factors (case A2 and A3) are the most likely candidates to account for the dynamical behavior observed in the Crab Nebula jet. In particular, the change in orientation of the jet recently noticed (Weisskopf et al. in preparation) can be reproduced by our simulations. Our findings are complementary to 3D numerical calculations (Mizuno et al., 2011; Porth et al., 2013) studying the effects of large values of in pulsar winds.

Future extensions of this work will take into account the full three-dimensional structure of the PWN allowing us to model the jet launching region and the collimation process. Furthermore, our results can be used as a starting point for a detailed analysis of the impulsive particle acceleration mechanism in the Crab Nebula as due to induced electric fields at the localized reconnection layers in the termination zone of the jet.


We thank Martin Weisskopf for providing Fig. 1 and for many stimulating discussions about the Crab Nebula morphology and the characteristics of the South East jet. A.M. wishes to thank Luca Del Zanna and Elena Amato for valuable comments on the initial jet configuration and G. Bodo for helpful discussions on current-driven instabilities. We acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support. This work was partially supported by ASI grants n. I/042/10/0, I/028/12/0.


  1. pagerange: Modelling the Kinked Jet of the Crab NebulaReferences
  2. pubyear: 2013
  3. The Crab Nebula average magnetic field is B mG.
  4. A bent jet seems to be a common feature of jets in PWN, as shown, e.g., by the jets of the Vela PWN (Durant et al., 2013)
  5. A collection of movies for the simulation cases presented here can be found at


  1. Abdo, A.A., et al. 2011, Science, 331, 739.
  2. Anile, A. M. 1990, Relativistic Fluids and Magneto-fluids, by A. M. Anile, pp. 348. ISBN 0521304067. Cambridge, UK: Cambridge University Press, February 1990.,
  3. Begelman, M. C. 1998, ApJ, 493, 291
  4. Bodo, G., Mamatsashvili, G., Rossi, P., and Mignone, A. , 2013, submitted to MNRAS
  5. Buehler, R., Scargle, J. D., Blandford, R. D., et al. 2012, ApJ, 749, 26
  6. Cerutti, B., Uzdensky, D. A., & Begelman, M. C. 2012, ApJ, 746, 148
  7. Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2013, arXiv:1302.6247
  8. Del Zanna L., Amato E., Bucciantini N., 2004, A&A, 421, 1063
  9. Del Zanna L., Volpi D., Amato E., Bucciantini N., 2006, A&A, 453, 621
  10. Durant, M., Kargaltsev, O., Pavlov, G. G., Kropotina, J., & Levenfish, K. 2013, Astrophysical Journal, 763, 72
  11. Kargaltsev, O., & Pavlov, G. G. 2008, 40 Years of Pulsars: Millisecond Pulsars, Magnetars and More, 983, 171
  12. Kennel, C. F., & Coroniti, F. V. 1984, Astrophysical Journal, 283, 710
  13. Komissarov, S. S. 1999, MNRAS, 308, 1069
  14. Komissarov S. S., Lyubarsky Y. E., 2003, MNRAS, 344, L93
  15. Komissarov S. S., Lyubarsky Y. E., 2004, MNRAS, 349, 779
  16. Komissarov, S. S. 2012, MNRAS, 422, 326
  17. Komissarov S. S., 2013, MNRAS, 428, 2459
  18. Lyubarsky, Y., & Kirk, J. G. 2001, Astrophysical Journal, 547, 437
  19. Lyubarsky, Y. E. 2002, MNRAS, 329, L34
  20. Lyubarsky, Y. E. 2012, MNRAS, 427, 1497
  21. Mignone, A., Bodo, G., Massaglia, S., Matsakos, T., Tesileanu, O., Zanni, C., & Ferrari, A. 2007, ApJS, 170, 228
  22. Mignone A., Ugliano M., Bodo G., 2009, MNRAS, 393, 1141
  23. Mignone A., Rossi P., Bodo G., Ferrari A., Massaglia S., 2010, MNRAS, 402, 7
  24. Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7
  25. Mizuno, Y., Lyubarsky, Y., Nishikawa, K.-I., & Hardee, P. E. 2011, ApJ, 728, 90
  26. Moll, R., Spruit, H. C., & Obergaulinger, M. 2008, Astronomy & Astrophysics, 492, 621
  27. Pétri, J., & Lyubarsky, Y. 2007, Astronomy & Astrophysics, 473, 683
  28. Porth, O. 2013, MNRAS, 429, 2482
  29. Porth, O., Komissarov, S. S., & Keppens, R. 2013, MNRAS, 431, L48
  30. Rossi, P., Mignone, A., Bodo, G., Massaglia, S., & Ferrari, A. 2008, Astronomy & Astrophysics, 488, 795
  31. Striani, E., Tavani, M., Piano, G., et al. 2011, ApJL, 741, L5
  32. Striani, E., Tavani, M., Vittorini, V., et al. accepted for pubblication on Astrophysical Journal
  33. Tavani, M., et al. 2011, Science, 331, 736
  34. Uzdensky, D. A., Cerutti, B., & Begelman, M. C. 2011, ApJL, 737, L40
  35. Vittorini, V., et al. 2011, ApJ, 732, L22
  36. Weisskopf, M. C., Hester, J. J., Tennant, A. F., et al. 2000, ApJL, 536, L81
  37. Weisskopf, M. 2011, Extreme and Variable High Energy Sky (Extremesky 2011)
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description