1 Introduction
Pulsar Wind Nebulae (PWNe) represent the best laboratory for relativistic plasma astrophysics and high-energy astrophysics, since the very extreme conditions which can be found in these objects cannot be obtained artificially in a laboratory. They are a particular class of supernova remnants, also known as plerions due to their centre filled emission, which makes them easily distinguishable from shell-type remnants, which on the contrary show a limb-brightened morphology.
At present, approximately 100 PWNe are known, and they all show a broad band spectrum, extending from radio to X-ray or even gamma-ray frequencies (Gaensler & Slane Reference Gaensler and Slane2006; Kargaltsev, Rangelov & Pavlov Reference Kargaltsev, Rangelov and Pavlov2013). Spectra are dominated by non-thermal synchrotron emission, due to relativistic particles moving in the nebular magnetic field, from radio up to a few MeV, and by inverse Compton scattering (IC) at higher energies. The radio component is characterized by a flat spectral index, typically in the range
$-0.3\lesssim \unicode[STIX]{x1D6FC}_{r}\lesssim 0$
, while the high-energy component has a much steeper spectrum, with typical photon indices
$\unicode[STIX]{x1D6E4}\sim 1{-}2$
: this implies the presence of one or more spectral breaks between the radio and X-ray frequencies (Gaensler & Slane Reference Gaensler and Slane2006).
PWNe are basically bubbles of relativistic particles and electromagnetic fields, born in the core collapse supernova explosion of massive stars (
$M\gtrsim 8M_{\odot }$
). The nucleus of the progenitor star collapses to a highly magnetized and rapidly rotating neutron star, often observed as a pulsar. When the supernova event terminates, the majority of the rotational energy lost by the neutron star fills the remnant with a magnetized, cold, relativistic plasma wind, mainly made of electron–positron pairs and possibly a small amount of hadrons (ions) (Gallant & Arons Reference Gallant and Arons1994). The pairs are copiously produced in the pulsar magnetosphere where the primary electrons, extracted from the star surface by the extremely strong induced electric field, interact with radiation and give rise to an electromagnetic cascade. The exact amount of pair production in the pulsar magnetosphere and its dependence on the pulsar parameters is one of the missing ingredients in our understanding of pulsar physics and a subject of active investigation (Hibschman & Arons Reference Hibschman and Arons2001; Timokhin & Harding Reference Timokhin and Harding2015). Constraints on this unknown quantity deriving from PWNe modelling will be extensively discussed in the § 2.2.
In the case of a young pulsar (usually Crab-like PWNe are observed for pulsars younger than
${\sim}20\,000$
yrs), the emitted wind is still surrounded by the denser and slower debris of the supernova explosion, which are expanding outward with non-relativistic velocity. The interaction with this material generates a reverse shock, that propagates back towards the pulsar. If the magnetic energy in the wind does not exceed by far the particle kinetic energy, then the reverse shock will propagate back only down to a distance where the wind ram pressure upstream equals the nebular pressure downstream. From this point on, this termination shock (TS hereafter) will slowly move outward, as the nebula expands, as long as the pulsar input can be considered constant in time. At the TS, the pulsar wind is slowed down and heated: its bulk energy is efficiently dissipated and it goes in accelerated particles, which will be responsible for the observed nebular synchrotron emission across the electromagnetic spectrum (Pacini & Salvati Reference Pacini and Salvati1973; Rees & Gunn Reference Rees and Gunn1974; Kennel & Coroniti Reference Kennel and Coroniti1984a
). Between
${\sim}10\,\%$
–30 % of the energy lost by the pulsar shows up as synchrotron emission of the daughter PWN (Kennel & Coroniti Reference Kennel and Coroniti1984b
), to be compared with the ten times smaller fraction that typically appears as gamma-ray pulsed emission (and in radio it is even much smaller). It is appropriate to notice that in the case of a strongly magnetized wind (Poynting flux much larger than the particle kinetic energy flux) this picture must be modified somewhat: a reverse shock hitting such a wind will likely produce little dissipation and possibly never reach a self-similar stage of expansion. In any case, in PWNe effective dissipation is usually inferred to take place in a very thin layer of plasma, which either means a strong magnetohydrodynamic (MHD) shock develops in a moderately magnetized medium or some non-MHD effects (likely magnetic reconnection) guarantee the same result.
Thanks to its proximity (
${\sim}2$
kpc) and to the high emitted power, the Crab nebula is unique among PWNe, and it is naturally considered as the prototype of the whole class. It is indeed one of the best studied objects in the sky, and the one from which we have learnt most of what we currently know about PWNe (Hester Reference Hester2008). The Crab nebula is almost certainly associated with a supernova explosion that happened in 1054 AD. Chinese astronomers reported that a ‘guest’ star appeared in the sky and that it was visible during daytime for three weeks and for 22 months at night (Stephenson & Green Reference Stephenson and Green2002, and references therein). The nebular remnant was then discovered in 1731 by John Bevis, and observed a few years later by Charles Messier, becoming the first item of his famous catalogue of non-cometary objects, with the name M1.
As inferred from energetic arguments by Pacini (Reference Pacini1967) prior to the actual discovery of pulsars, the central engine of the Crab nebula was indeed identified as a pulsar rotating with period
$P=33$
ms, and slowly increasing with time (
${\dot{P}}=4.21\times 10^{-13}~\text{s}~\text{s}^{-1}$
). This observed spin-down implies that the rotational energy is dissipated at a rate of
$L\sim 5\times 10^{38}~\text{erg}~\text{s}^{-1}$
, equivalent to
$1.3\times 10^{5}~L_{\odot }$
, a value similar to the inferred rate at which the energy is supplied to the nebula (Gold Reference Gold1969).
After the launch of the Chandra X-ray observatory, the high energy morphology of the Crab nebula and of many other PWNe became accessible, and many different variable and bright features were observed in the inner regions with great detail. A puzzling and unexpected jet–torus structure was identified in many objects, including the Crab nebula (Hester et al.
Reference Hester1995; Helfand, Gotthelf & Halpern Reference Helfand, Gotthelf and Halpern2001; Pavlov et al.
Reference Pavlov, Teter, Kargaltsev and Sanwal2003; Gaensler, Pivovaroff & Garmire Reference Gaensler, Pivovaroff and Garmire2001; Gaensler et al.
Reference Gaensler, Arons, Kaspi, Pivovaroff, Kawai and Tamura2002; Lu et al.
Reference Lu, Wang, Aschenbach, Durouchoux and Song2002). The Crab torus, in particular, appeared to be composed of several rings, the brightest of which is located at a distance of
${\sim}0.1$
pc from the pulsar. This ring, usually referred to as ‘inner ring’, is likely associated with the TS location (Weisskopf et al.
Reference Weisskopf2000).
In the optical band, Scargle (Reference Scargle1969) had previously identified many bright, arc-like, variable structures, that he named ‘wisps’. These are currently known to fluctuate in brightness on very short time scales (of the order of ks (Hester et al.
Reference Hester, Mori, Burrows, Gallagher, Graham, Halverson, Kader, Michel and Scowen2002; Melatos et al.
Reference Melatos, Scheltus, Whiting, Eikenberry, Romani, Rigaut, Spitkovsky, Arons and Payne2005)). They are also seen to move outward with mildly relativistic velocities (
${\lesssim}0.5c$
, with
$c$
the speed of light), and progressively fade. The typical period between their appearance close to the inner ring and disappearance in the bulk of the nebula can be of the order of hours, days, months or even years (Hester et al.
Reference Hester, Mori, Burrows, Gallagher, Graham, Halverson, Kader, Michel and Scowen2002; Bietenholz et al.
Reference Bietenholz, Hester, Frail and Bartel2004; Melatos et al.
Reference Melatos, Scheltus, Whiting, Eikenberry, Romani, Rigaut, Spitkovsky, Arons and Payne2005), depending on the wavelength of observation (infra-red (IR), X-rays, optical and radio respectively). As we will see in the following, since wisps originate so close to the supposed TS location, where the flow is highly turbulent and swirling, they represent the best features to investigate the physical conditions of the plasma and of the emitting particles near to the shock surface.
Other variable features, the so-called ‘knots’, had been later observed by van den Bergh & Pritchet (Reference van den Bergh and Pritchet1989) as point-like structures, that always appear and fade out at the same locations. These are visible at basically all wavelengths. In the optical band the most luminous one (called knot-1) is located at
$0.65^{\prime \prime }$
from the pulsar, while the second brightest one is at a distance of
$3.8^{\prime \prime }$
. All knots appear to be aligned with the polar jet in the south-east (SE) direction.
The presence of polar X-ray jets is certainly one of the most intriguing features in the high-energy morphology of PWNe. They appear to originate so close to the pulsar location that, when first discovered, it was natural to think that any collimation mechanism had to operate upstream of the TS, in the pulsar wind region. The best candidate among the possible models was soon recognized to be collimation by the hoop stresses due to the dominant toroidal magnetic field component in the wind, but this can be shown to be very inefficient inside an ultra-relativistic MHD flow, such as the pulsar wind before the TS. Later Lyubarsky (Reference Lyubarsky2002) suggested that an anisotropic distribution of the energy flux in the wind could cause the shock to be highly non-spherical, rather developing a characteristic oblate shape (see also Bogovalov & Khangoulyan Reference Bogovalov and Khangoulyan2002). The TS is more expanded in the equatorial direction, where the energy flux is higher, but is much closer to the pulsar along the polar axis. Upstream collimation of the flow is then only apparent, while in reality it is the downstream flow that is diverted towards the axis by hoop stresses and forms a jet. Thus, the collimation mechanism is actually operating in the downstream PWN plasma, not in the pulsar wind, and theoretical difficulties are avoided.
This scenario was soon successfully confirmed by means of axisymmetric relativistic MHD numerical simulations of PWNe (Del Zanna, Amato & Bucciantini Reference Del Zanna, Amato and Bucciantini2004; Komissarov & Lyubarsky Reference Komissarov and Lyubarsky2004; Bogovalov et al. Reference Bogovalov, Chechetkin, Koldoba and Ustyugova2005). MHD numerical modelling has progressively become more and more sophisticated, first with the introduction of increasingly more accurate emission diagnostics, and more recently with the upgrade to full three-dimensional (3-D) simulations. The questions that one has attempted to answer have also become deeper. The description of the most recent progresses and of the physical implications of the numerical findings will be the subject of the present paper.
2 The big open problems in PWN physics
The final purpose of the increasingly refined methods of analysis that have been developed to study the relativistic plasma dynamics in PWNe is to answer a number of old and new questions, the most important of which, in our view, are shortly reviewed in this section.
2.1 The pulsar wind magnetization
The fundamental parameter that characterizes the physical behaviour of the (cold) pulsar wind, and hence also the properties of the PWN arising from its interaction with the surrounding supernova remnant material, is the magnetization
$\unicode[STIX]{x1D70E}$
. This is defined as the ratio between the Poynting and the particle kinetic energy fluxes in the wind:
$$\begin{eqnarray}\unicode[STIX]{x1D70E}=\frac{B^{2}}{4\unicode[STIX]{x03C0}nm_{e}\unicode[STIX]{x1D6FE}^{2}c^{2}},\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FE},B$
and
$n$
are the Lorentz factor, the magnetic field and the rest-frame particle number density of the flow, and
$m_{e}$
is the electron mass (here we consider a wind with kinetic energy flux dominated by pairs). This quantity is expected to vary with distance from the central pulsar, but also with latitude from the equatorial plane of the pulsar rotation. According to theoretical models of pulsar magnetospheres, the wind is expected to be highly magnetized near the pulsar light cylinder, where clearly
$\unicode[STIX]{x1D70E}\gg 1$
. On the other hand, the very first attempts at modelling the Crab nebula within the framework of ideal MHD, showed that effective deceleration of the flow at the TS is only possible if its magnetization is not too high. In particular, steady-state analytical MHD models of the Crab nebula (Rees & Gunn Reference Rees and Gunn1974; Kennel & Coroniti Reference Kennel and Coroniti1984a
; Begelman & Li Reference Begelman and Li1992) require a value of
$\unicode[STIX]{x1D70E}\sim v_{n}/c$
, with
$v_{n}$
the expansion velocity of the PWN, always much less than one (which in the case of the Crab nebula is
$\unicode[STIX]{x1D70E}\sim 3\times 10^{-3}$
). Assuming that the downstream flow is well described by ideal MHD, strong magnetic dissipation is then needed between the light cylinder and the wind TS.
Since pulsars are oblique rotators, the wind has a complex structure: the current sheet that separates magnetic field lines of opposite polarities twists and tangles, cutting the equatorial plane along twin spirals. These oscillations take place within a region around the equator of angular extent equal to the inclination angle between the rotational and magnetic axes of the pulsar. Neighbouring stripes of opposite polarities offer a perfect location for magnetic dissipation to occur, and hence
$\unicode[STIX]{x1D70E}$
can be substantially decreased in that region (Coroniti Reference Coroniti1990).
Nevertheless, in the case of the Crab nebula, the difference between the value of
$\unicode[STIX]{x1D70E}$
expected at the light cylinder and that estimated at the TS in a model similar to that of Kennel & Coroniti requires an efficiency of conversion of magnetic field energy into particle kinetic energy which appears too high to be explained by magnetic reconnection in the striped wind (Lyubarsky & Kirk Reference Lyubarsky and Kirk2001). In the last decade, 2-D axisymmetric MHD simulations actually pointed out that at the TS the magnetization is likely higher than estimated based on steady-state 1-D MHD, and in particular, for the Crab PWN
$\unicode[STIX]{x1D70E}\gtrsim 0.01$
is required at the TS in order to generate polar jets (Del Zanna et al.
Reference Del Zanna, Volpi, Amato and Bucciantini2006). However, even in the framework of 2-D MHD, the estimated values of
$\unicode[STIX]{x1D70E}$
at the TS are still much smaller than unity, leaving the problem of efficient magnetic dissipation along the way to the shock.
Recently Porth, Komissarov & Keppens (Reference Porth, Komissarov and Keppens2014b
) presented the first 3-D numerical simulation of the Crab nebula, showing that values of the magnetization greater than unity can be worked out by just allowing for non-axisymmetric models, by increasing the dimensionality of the simulations. Downstream of the shock, the onset of kink instabilities leads to an efficient conversion of a large fraction of the toroidal magnetic field in the nebula into a poloidal component, reducing the magnetic tension, and also producing a larger amount of dissipation with respect to that observed in two dimensions. These simulations, however, have been run only for approximately
$1/10$
of the Crab lifetime (
$\simeq 1000$
yr), and end when the system is still far from the self-similar expansion phase, since the ratio between the nebula radius and the TS radius is not yet saturated. In particular, it is not clear yet whether the average nebular magnetic field strength will be close to that inferred from observations (
${\sim}200~\unicode[STIX]{x03BC}\text{G}$
) at the final stage of the evolution, and if the synthetic non-thermal emission will match observational data in all spectral bands, thus confirming the model assumptions on the pulsar wind and emitting particles properties. Further investigation is certainly needed in this respect.
2.2 Origin of radio emitting particles and implications for pulsar multiplicity
Although radio emitting particles are dominant by number in the Crab nebula, their origin is still a mystery. Since pulsars are believed to be primary antimatter factories in the Galaxy, reasonably good knowledge of the amount of particles generated in their magnetospheres has fundamental implications, not only for pulsar physics, but also for quantifying their contribution to the measured cosmic ray positron excess (Adriani et al.
Reference Adriani2009; Aguilar et al.
Reference Aguilar2013), and to assess the role of other possible sources, as e.g. dark matter annihilation. The pair production in pulsar magnetospheres is still not completely understood, but there is a general consensus that each primary electron extracted from the star surface generates a high number
$\unicode[STIX]{x1D705}$
of pairs (the so called ‘pulsar pair multiplicity’) in the cascading process. For young, energetic pulsars, the expected value of the multiplicity is in the range
$10^{3}\lesssim \unicode[STIX]{x1D705}\lesssim 10^{7}$
.
Since secondary particles become part of the pulsar wind, a direct estimate of
$\unicode[STIX]{x1D705}$
can be obtained by accurate radiation modelling of the nebula. An important thing to notice is that ions can only be primary, i.e. extracted from the star surface, since they cannot be generated in electromagnetic cascades. As a consequence they can only be at most a fraction
${\sim}1/\unicode[STIX]{x1D705}$
of the total number of the wind pairs. Nevertheless, if
$\unicode[STIX]{x1D705}<m_{i}/m_{e}\approx 10^{3}-10^{4}$
, with
$m_{i}$
the ion mass, the hadronic component of the wind would be energetically dominant even if numerically negligible in the wind.
One of the most important physics issues connected with the origin of the radio emitting particles is the fact that this origin is directly connected to the value of the pair multiplicity
$\unicode[STIX]{x1D705}$
. Estimates of
$\unicode[STIX]{x1D705}$
change by approximately two orders of magnitude depending on whether radio emitting particles are considered as part of the pulsar outflow or not. The first MHD models and high-energy observations of the Crab suggested that
$\unicode[STIX]{x1D705}\sim 10^{4}$
, assuming that only particles responsible for high-energy emission would be part of the pulsar outflow, while radio emitting particles are considered to be of different origin, possibly generated in a primordial outburst of the pulsar or elsewhere in the nebula (Kennel & Coroniti Reference Kennel and Coroniti1984b
; Atoyan & Aharonian Reference Atoyan and Aharonian1996; Gaensler et al.
Reference Gaensler, Arons, Kaspi, Pivovaroff, Kawai and Tamura2002). It is also important to notice that in this scenario the hadronic component of the wind, if present, is expected to be energetically dominant. Interestingly in the case of Crab these would be PeV ions.
On the contrary, if radio particles are considered as part of the pulsar outflow, the expected pulsar pair multiplicity is much higher, namely
$\unicode[STIX]{x1D705}\sim 10^{6}$
, in which case the energy flux carried by a possible hadronic component would be completely irrelevant (Bucciantini, Arons & Amato Reference Bucciantini, Arons and Amato2011).
It is thus clear that discriminating between the two scenarios is of fundamental importance for a correct modelling of the pulsar wind composition.
2.3 Particle acceleration mechanism at the wind TS
It is a matter of fact that the Crab nebula is an efficient particle accelerator, with evidence of particles accelerated up to PeV energies (Arons Reference Arons2012). The puzzling issue with this is the fact that relativistic magnetized shocks, as the TS is, are very hostile environments for particle acceleration, and the mechanism at the basis of the observed particle spectra is still a mystery.
The observed spectrum of the Crab nebula, similar to other PWNe, seems to suggest different acceleration mechanisms at work for low- and high-energy emitting particles, since it can be roughly modelled by assuming a broken power-law distribution function of the injected ultra-relativistic emitting particles,
$N(E)\propto E^{-p}$
, with
$p\sim 1.5$
at low energies (for particles with
$E\lesssim 100$
GeV) and
$p\sim 2.2$
at higher ones (up to a few TeV).
The very flat slope of the radio emitting component is compatible with driven magnetic reconnection at the TS. This mechanism was investigated by Sironi & Spitkovsky (Reference Sironi and Spitkovsky2011), who performed particle-in-cell (PIC) simulations of a striped relativistic wind impinging on a shock, a structure similar to what is expected to occur at the pulsar wind TS. They assume that the striped morphology of the pulsar wind is maintained all the way to the TS, where compression leads the field to reconnect. Many reconnection islands form in the flow, where the unscreened electric field is able to accelerate particles. The properties of the resulting spectrum depend on the flow magnetization
$\unicode[STIX]{x1D70E}$
and on the ratio between the wavelength of the stripes and the particles Larmor radii. This ratio is connected to the value of the pulsar pair multiplicity
$\unicode[STIX]{x1D705}$
, and the result in the case of the Crab nebula is that in order to reproduce the observed radio spectrum, a magnetization
$\unicode[STIX]{x1D70E}\gtrsim 30$
and a multiplicity of
$\unicode[STIX]{x1D705}\gtrsim 10^{8}$
are needed. Such large multiplicity is well above the current predictions of pulsar magnetosphere models (Timokhin & Harding Reference Timokhin and Harding2015), and in any case it makes it difficult to draw a self-consistent picture of the system, since a wind with the required high number of pairs is expected to undergo efficient dissipation well before the TS (Kirk & Skjæraasen Reference Kirk and Skjæraasen2003), destroying the stripes before the forced reconnection might start to work (Amato Reference Amato2014).
One possibility to alleviate the requirements on
$\unicode[STIX]{x1D705}$
is to accelerate particles near the polar cusps of the TS, where the shock front is closer to the pulsar, and as a consequence, the particle density is much higher (it decreases as
${\sim}1/r^{2}$
with distance from the pulsar). However, on the other hand, no stripes are expected to be present at these locations unless the pulsar is an almost orthogonal rotator, i.e. with the angle between the magnetic and rotation axes close to
$90^{\circ }$
.
A different mechanism must be invoked to account for the much steeper high-energy component of the spectrum. A particle injection with slope
$p\simeq 2.2$
is in fact what one expects from diffusive shock acceleration (DSA), i.e. Fermi I-like acceleration at relativistic shocks. Contrary to driven magnetic reconnection, DSA is effective only when the magnetization is sufficiently low, in particular
$\unicode[STIX]{x1D70E}\lesssim 10^{-3}$
(Sironi & Spitkovsky Reference Sironi and Spitkovsky2009). Since MHD simulations of the Crab nebula lead to values of
$\unicode[STIX]{x1D70E}$
of at least one order of magnitude greater, DSA can only be effective in a few sectors of the shock, where magnetic dissipation is sufficiently strong to ensure that locally the plasma is almost unmagnetized. It is not clear however if a sufficiently large fraction of the flow satisfies this condition. This again depends on the width of the striped region and on the efficiency of magnetic dissipation. According to current results of the 2-D MHD simulations that best reproduce the Crab nebula X-ray morphology, only few % of the wind energy would be carried by a sufficiently low
$\unicode[STIX]{x1D70E}$
flow (Amato Reference Amato2014).
An alternative acceleration mechanism that received some attention in the literature (Hoshino et al. Reference Hoshino, Arons, Gallant and Langdon1992; Amato & Arons Reference Amato and Arons2006) is resonant absorption of ion-cyclotron waves in an ion-doped plasma (RCA hereafter). Its viability poses no special requirements on the wind magnetization, but does require that most of the energy of the wind is carried by ions. The idea is that pairs would be accelerated by resonant absorption of the cyclotron radiation emitted by ions in the wind, when these are set into gyration by the enhanced magnetic field at the shock crossing. The outcome of this acceleration process in terms of efficiency, maximum achievable energy and accelerated particle spectrum depends on the relative amount of wind energy that is carried by ions. In order for RCA to provide any acceleration, the ions must be energetically dominant and therefore the mechanism is not viable if the multiplicity is larger than the pair to ion mass ratio, as would typically be the case if radio emitting electrons are part of the pulsar outflow (Amato Reference Amato2014). On the other hand, it is not straightforward to make a connection with the flow properties that are relevant for MHD modelling and single locations, along the shock front, where RCA is more likely to work. This is why in the following we will mostly consider the first two processes discussed above when addressing the question of what MHD simulations can tell us about the acceleration sites of particles in different energy ranges.
2.4 The gamma-ray flares
For a very long time the Crab nebula was believed to be a very stable source at both X-ray and gamma-ray energies, and hence it was considered as the standard candle and calibrator for high-energy astrophysics instrumentation. This was one of the reasons why the evidence of its gamma-ray flaring activity, first observed in September 2010 by the AGILE satellite, and then confirmed by FERMI, was so astonishing to the community (Abdo et al.
Reference Abdo2011; Tavani et al.
Reference Tavani2011). Additional similar flaring events with a sudden release of an isotropic energy of
$E_{\unicode[STIX]{x1D6FE}}\sim 10^{41}$
erg have then been observed in the following years, confirming that the Crab is not a stable source at all (Striani et al.
Reference Striani2011, Reference Striani2013; Buehler et al.
Reference Buehler2012; Becerra, Buehler & Hays Reference Becerra, Buehler and Hays2014).
During the flaring phase, the gamma-ray emission is enhanced dramatically in the spectral range between 50–100 MeV and a few GeV, basically the region where the steady-state synchrotron radiation drops and the IC contribution rises. The newly discovered non-thermal emission seems best modelled as due to a quasi-mono-energetic particle distribution, or a relativistic Maxwellian with Lorentz factor of order a few
$10^{9}$
(Buehler et al.
Reference Buehler2012). We are thus in the presence of an extremely powerful cosmic accelerator, reaching energies above 1 PeV!
The timing of these flares is quite impressive. The increase in the gamma-ray light curve can be as fast as just a few hours, while the flares typically last from a few days up to a couple of weeks. Together with the extreme energetics involved (luminosities beyond 100 MeV up to
$L_{\unicode[STIX]{x1D6FE}}\sim 10^{36}~\text{erg}~\text{s}^{-1}$
), this clearly indicates that the source of the emission must be located in the innermost region of the Crab nebula, most probably around the TS. Unfortunately, in spite of several attempts, no clear indication of a parent candidate location has been found yet, among the variable features observed in the optical and X-ray bands (e.g. the knots): none of these shows any special variation, during the gamma-ray flares (Weisskopf et al.
Reference Weisskopf2013; Bietenholz et al.
Reference Bietenholz, Yuan, Buehler, Lobanov and Blandford2015; Rudy et al.
2015).
The current favourite physical interpretation of the gamma-ray flares relies on acceleration of particles in an electric field beyond the ideal limit, since the constraint
$E\leqslant B$
of relativistic MHD leads to an emission cutoff at approximately
$200$
MeV. The obvious choice among the non-ideal effects is certainly that of sudden and powerful reconnection events in a magnetically dominated plasma, as appropriate for some regions in the post-shock nebular flow, in current sheets with sizes corresponding to the variable features observed in optical and X-ray bands in the innermost regions of the Crab nebula.
As far as particle acceleration is concerned, state-of-the-art 2-D and 3-D particle-in-cell numerical simulations of reconnection in relativistic pair plasmas, including radiative feedback, confirmed the possibility of acceleration of particles, with the correct energetics, higher magnetic fields (
${\sim}1$
mG) and slope, in the narrow reconnection layers (current sheets) (Cerutti et al.
Reference Cerutti, Werner, Uzdensky and Begelman2014). Acceleration is allowed beyond the radiation reaction limit, with a predicted flux in the gamma-rays above 100 MeV exceeding the quiescence level by factors 2–3, as required to reproduce the flare events in the Crab nebula (Cerutti et al.
Reference Cerutti, Werner, Uzdensky and Begelman2013, Reference Cerutti, Werner, Uzdensky and Begelman2014). The dynamics is dominated by driven tearing and by drift-kink instabilities, and the flare event is seen when the beam of accelerated particles is close to the assumed line-of-sight in the periodic simulation box (Sironi & Spitkovsky Reference Sironi and Spitkovsky2014; Kagan et al.
Reference Kagan, Sironi, Cerutti and Giannios2015; Lyutikov et al.
Reference Lyutikov, Sironi, Komissarov and Porth2016; Sironi, Giannios & Petropoulou Reference Sironi, Giannios and Petropoulou2016; Yuan et al.
Reference Yuan, Nalewajko, Zrake, East and Blandford2016).
3 Results from 2-D simulations
The observed complex X-ray morphology of PWNe, characterized by the presence of an equatorial torus and polar jets, was accounted for by different groups by means of 2-D MHD axisymmetric models, as a natural consequence of the pulsar wind anisotropy and striped structure (Del Zanna et al. Reference Del Zanna, Amato and Bucciantini2004; Komissarov & Lyubarsky Reference Komissarov and Lyubarsky2004; Bogovalov et al. Reference Bogovalov, Chechetkin, Koldoba and Ustyugova2005). Moreover the high-energy emission features were reproduced with great accuracy down to very fine detail (Del Zanna et al. Reference Del Zanna, Volpi, Amato and Bucciantini2006; Camus et al. Reference Camus, Komissarov, Bucciantini and Hughes2009), as well as polarization (Bucciantini et al. Reference Bucciantini, del Zanna, Amato and Volpi2005) and the complete emission spectrum (Volpi et al. Reference Volpi, Del Zanna, Amato and Bucciantini2008). In this latter work the spectrum was computed including the gamma-ray emission component, due to the IC scattering process. The IC spectrum was computed by considering the extended Klein–Nishina differential cross-section and various photon targets for relativistic particles, namely synchrotron emitted photons, thermal photons coming from dust emission and photons from the cosmic microwave background. Particles responsible for low-energy emission (radio particles) and high-energy emission (optical/X-ray particles) were considered as part of two distinct families, the first extending up to an energy of approximately 100 GeV and the second starting from an energy approximately 2–3 times larger. This particle spectrum automatically ensures the appearance of a spectral break in the emission spectrum between the infrared and the optical bands. A second break, at ultra-violet (UV) frequencies, also arose naturally from the synchrotron burn-off effect.
It was pointed out in that work that the reduced dimensionality of the simulation reflects in an artificial compression of the magnetic field around the polar axis, such that in order to reproduce the high-energy morphology of the Crab nebula, a smaller value of the average magnetic filed is found (
${\sim}50~\unicode[STIX]{x03BC}\text{G}$
), which is approximately
$3/4$
times smaller than that estimated from the integrated emission spectrum (de Jager & Harding Reference de Jager and Harding1992). With this low average field, electrons are affected by weaker synchrotron losses and a steeper spectral index (of approximately
$0.8$
) is required in order to fit the high-energy component of the spectrum. At the same time, with a lower field, a larger number of particles is needed in order to reproduce the synchrotron emission and this reflects into an overestimate of the IC component by the corresponding factor. The conclusion of that work was that, in order to reproduce the whole spectrum self-consistently, a larger and more diffuse magnetic field is mandatory, and 3-D models are indeed needed.
Two-dimensional MHD models have also dealt with the question of time variability at different wavelengths, showing that wisps can be reproduced in terms of time scales, morphological patterns and outward velocities (Camus et al. Reference Camus, Komissarov, Bucciantini and Hughes2009; Olmi et al. Reference Olmi, Del Zanna, Amato, Bandiera and Bucciantini2014). In these models, wisps originate very close to the TS shock front, in the region where eddies of high velocity develop.
3.1 Investigating the origin of radio emitting particles
For many years, MHD numerical models have mainly been concerned with the properties of high-energy emission, while ignoring emission at lower energies. As previously mentioned, clarifying the origin of radio emitting particles is fundamental to obtain a correct estimate of the value of the pulsar pair multiplicity
$\unicode[STIX]{x1D705}$
and, as a direct consequence, constraints on the particle acceleration process.
In a recent study Olmi et al. (Reference Olmi, Del Zanna, Amato, Bandiera and Bucciantini2014) have considered different possibilities for the origin of radio particles, testing whether it is possible or not to discriminate their origin based on the resulting emission morphology. The question is whether radio particles are part of the pulsar outflow, and accelerated at the TS together with high-energy emitting particles, or if they are of different origin and accelerated somewhere else in the nebula.
Three different hypotheses were then considered: radio particles are injected in the nebula as part of the pulsar outflow (case A); radio particles are uniformly distributed in the nebula, produced somewhere else, for instance in the thermal filaments, and continuously accelerated by the interaction with local turbulence (case B); radio particles are a relic population, born in a primordial outburst of the pulsar during its early life, and still present thanks to their very long lifetimes against synchrotron losses (case C).
Each scenario was tested by simulating the entire evolution of the nebula, with the suitable numerical tracers following the energy evolution of radio particles. Emission properties were then computed on top of the numerical data corresponding to an age of the Crab nebula of approximately 1000 yrs. Results obtained in the three different scenarios are shown in figure 1. When particles are injected only for a short time during the initial stages of the nebular evolution (case C), with the injection limited to
$t\leqslant 100$
yr, emission maps are clearly different from what observations show, as can be seen from the right-most panel of figure 1. Here radio emission mainly comes from the outer shell of the nebula, while the inner region results to be much fainter than observed. This is the result of the advection of radio particles during the early stages of evolution. Particles are advected by large-scale flow vortexes to the outer region of the nebula, and the absence of freshly injected particles at later times results into a lack of emitting particles in the central region, since all of them are confined to the outer shell. Kinetic diffusion, which is not included in the model, is however expected to have negligible effects on these low-energy particles, while turbulent flow diffusion might be more relevant. The inner nebula can be somehow enriched in primordial electrons by the effect of the Rayleigh–Taylor (RT) instabilities that set in at its outer boundary (Porth, Komissarov & Keppens Reference Porth, Komissarov and Keppens2014a
). In the discussed simulations the resolution is not sufficiently high to allow this mechanism to produce significant mixing, and effectively enrich the inner nebula with primordial electrons. In any case we do not expect the effects of turbulent RT mixing to be sufficient to reconcile our case C with radio observations (see e.g. Bandiera, Neri & Cesaroni Reference Bandiera, Neri and Cesaroni2002), which require that the electrons fully concentrated in the outer layers should become spread more or less uniformly in the whole nebula, down to the very inner regions, where the filaments do not even reach (they are seen to penetrate inside to
$1/3$
of the radius).
Figure 1. Radio intensity maps at
$\unicode[STIX]{x1D708}=1.4$
GHz, with intensity given by the colour bar in units of
$\text{mJy}~\text{arcsec}^{-2}$
. Emission is computed at
${\sim}1000$
yr and is plotted in linear scale and normalized to its maximum value. (a) Radio emitting particles are injected as part of the pulsar outflow into the nebula (case A). (b) Radio particles are considered as uniformly distributed in the nebula (case B). In all cases the maps are treated as in Bandiera et al. (Reference Bandiera, Neri and Cesaroni2002), subtracting small-scale structures. (c) Radio particles are injected in the nebula at early times and advected with the flow, without diffusion or any further acceleration (hypothesis C).
Therefore the hypothesis of an early burst is not a viable description of the origin of the radio emission from the Crab nebula. A caveat to this conclusion is that it is obtained by means of 2-D simulations. While in the inner region of the nebula we expect the 2-D description to be a good approximation of the dynamics, given that the kink instabilities that are suppressed by the reduced dimensionality are unlikely to have the time to develop, on the larger scales the dynamics is shown to be rather different in three dimensions, as already mentioned. Notice that the axisymmetric model is robust within a few TS radii, while moving further from the symmetry axis, the differences with the full 3-D model become more and more remarkable. A more detailed discussion on this can be found later on in § 4. The possibility exists that in 3-D, case C will look similar to our current case B, thanks to more effective mixing and wandering of magnetic field lines. It is clear that an analogous study will need to be done in three dimensions before the final word against a primordial origin of the radio emitting particles can be spelled. On the other hand, cases A and B lead to almost undistinguishable emission maps, and comparison with observations does not allow one to discriminate between the two cases. Moreover, the analysis of the wisp time variability in the inner region of the nebula is also found to be very similar in the two cases (figure 2). Wisps appear with similar features and outward velocities whether radio particles are considered as uniformly distributed in the nebula or injected and accelerated at the TS as the high-energy ones.
Figure 2. Appearance of wisps at radio frequency in case A (a) and case B (b). Wisps maps are obtained by subtracting images of the emission separated by a time interval of approximately 2 months. Outward motion of the wisps can be observed as consecutive darker and lighter arcs.
This proves that the existence of wisps does not imply that emitting particles are injected and accelerated at the shock. In fact wisps arise in the simulated maps as a consequence of the intrinsic properties of the MHD flow: bright moving features are associated with local enhancements of the magnetic field and Doppler boosting of material moving toward the observer.
3.2 Investigating the particle acceleration mechanism at the shock
Within a MHD description wisps arise thus as a natural consequence of the underlying flow properties. Therefore the fact that wisps are observed to be not coincident at the different wavelengths suggests a difference in the acceleration sites of the particles responsible of such emission.
As previously mentioned, different acceleration mechanisms require very different physical conditions to be effective, and thus identifying the location at which particles have been accelerated can put strong constraints on the mechanism at work.
In Olmi et al. (Reference Olmi, Del Zanna, Amato and Bucciantini2015) different possibilities were considered for acceleration of radio, optical and X-ray particles. In particular two distribution functions were defined for accounting for radio and X-ray particles. Since the distribution functions were assumed in such a way as to fit the complete spectrum of the Crab nebula (i.e. cutoff energies and normalization constants are chosen in order to obtain the best fit of the spectrum), the optical component is automatically determined as a superposition of the radio and X-ray contributions.
Three different hypotheses were considered for injection of particles: in case I particles are injected uniformly at the shock surface; in case II particles are injected in a wide equatorial sector or in the complementary narrow polar one; in case III, emitting particles are injected in a narrow equatorial region, or in the complementary wide polar one.
The entire evolution of the Crab nebula was then simulated by introducing as many numerical tracers as necessary in order to follow the evolution of particles injected at the different locations. In particular, each family could be injected in five different zones: wide or narrow polar caps, wide or narrow equatorial belt or along the entire TS surface. The entire system was evolved up to the actual age of the Crab (
${\sim}1000$
years), and wisps’ profiles were extracted from synthetic intensity maps during a period of 10 years at the end of the evolution, with monthly frequency.
For a direct comparison with wisps’ observations published by Hester et al. (Reference Hester, Mori, Burrows, Gallagher, Graham, Halverson, Kader, Michel and Scowen2002), Bietenholz et al. (Reference Bietenholz, Hester, Frail and Bartel2004) and Schweizer et al. (Reference Schweizer, Bucciantini, Idec, Nilsson, Tennant, Weisskopf and Zanin2013), emission maps were convolved with the appropriate instrumental point spread functions (PSFs). The adopted full width at half maximum (FWHM) was:
$0.5^{\prime \prime }$
for X-rays (from Chandra),
$1.8^{\prime \prime }$
for radio (Very Large Array) and
$0.75^{\prime \prime }$
for optical emission (North Optical Telescope).
Following the analysis by Schweizer et al. (Reference Schweizer, Bucciantini, Idec, Nilsson, Tennant, Weisskopf and Zanin2013) intensity peaks were extracted from a
$3^{\prime \prime }$
wide slice around the polar axis in the upper hemisphere of each map, where wisps are more prominent. For ease of comparison only peaks with intensity greater than
$1/3$
of the intensity maximum were taken into account. Results are shown as plots of the radial distance (from the pulsar, in arc seconds) of the local intensity maxima as a function of time (in months).
Figure 3. Radial positions of the local intensity maxima (in arcseconds) as a function of time (in months) with orange diamonds identifying radio wisps (
$\unicode[STIX]{x1D708}_{r}=5$
GHz), green crosses optical ones (
$\unicode[STIX]{x1D708}_{o}=3.75\times 10^{14}$
Hz) and light blue circles for X-rays (1 keV). On (a) case I is shown. On (b) case III is shown, with X-ray particles injected in the equatorial zone and radio ones injected in the complementary sector.
In figure 3 wisps obtained in the case of uniform injection or injection in a wide (narrow) equatorial (polar) region are compared. As expected, when particles of different families are injected at the same location, i.e. across the entire shock surface as in case I, wisps appear to be coincident at the different wavelengths. But as soon as they are injected at distinct locations (case II and III) wisps at radio, optical and X-ray frequencies are no more coincident. Here we only show wisps obtained under case III hypothesis (with X-ray particles injected in the equatorial region and radio ones in the opposite polar zone), while for a complete discussion of the wisp properties obtained in the different scenarios we refer to Olmi et al. (Reference Olmi, Del Zanna, Amato and Bucciantini2015). The main result of the analysis presented in the cited paper is that, in order to reproduce the absence of X-ray wisps in the region within
${\sim}6^{\prime \prime }$
from the pulsar, as observed in Schweizer et al. (Reference Schweizer, Bucciantini, Idec, Nilsson, Tennant, Weisskopf and Zanin2013), X-ray particles must be injected in a narrow equatorial sector, approximately coincident with the striped region of the wind where dissipation of magnetic field is expected to be most efficient.
This may indicate that X-ray particles are produced via Fermi I acceleration in the striped zone, where magnetization can be lowered enough by dissipation so as to allow this mechanism to be effective.
In Schweizer et al. (Reference Schweizer, Bucciantini, Idec, Nilsson, Tennant, Weisskopf and Zanin2013) it was also found that optical wisps are characterized by narrower angular profiles than X-ray wisps. Already in that article, two possible explanations were suggested for this finding: either a stronger Doppler boosting of the optical emitting region, or a higher degree of anisotropy of the optical emitting particles. A relevant investigation in this sense was carried out by Yuan & Blandford (Reference Yuan and Blandford2015), who showed that very different emission patterns can result from different assumptions on the emitting particle pitch angle distributions. In the MHD modelling by Olmi et al. (Reference Olmi, Del Zanna, Amato and Bucciantini2015) the particle distribution function was assumed to be isotropic at all energies and the model was found to correctly reproduce the angular size of the X-ray wisps, but overestimated the extension of the optical ones. This favours the idea that if the wisps are to be explained as a result of the properties of the MHD flow alone, then the X-ray emitting particles must be more isotropic than optical emitting ones, a fact that can be used in principle to put constraints on the acceleration mechanism, though this goes beyond the scopes of the present work.
Moving to lower frequencies, strong constraints on radio emission are again difficult to derive: the only case that can be easily excluded is the one in which radio particles are injected in a narrow polar cone, since no wisps appear to be produced. Since radio wisps must not be coincident with X-ray ones, the remaining possibilities are that acceleration of low-energy particles happens in the complementary polar region or in a wider equatorial stripe.
A possible interpretation of this finding is that radio particles could result from driven magnetic reconnection primarily occurring at moderately high latitudes, where conditions for driven magnetic reconnection to operate as an acceleration mechanism might be locally satisfied.
Going beyond the qualitative conclusions this study has provided requires an effort to compute the particle spectrum that the study by Sironi & Spitkovsky (Reference Sironi and Spitkovsky2011) would predict as a function of latitude along the shock surface for a given value of the inclination between magnetic and rotation axis. In the meantime, however, recent developments in pulsar and pulsar wind theory (Philippov, Spitkovsky & Cerutti Reference Philippov, Spitkovsky and Cerutti2015; Tchekhovskoy, Philippov & Spitkovsky Reference Tchekhovskoy, Philippov and Spitkovsky2016) are somewhat modifying the picture of the wind that was assumed within the modelling summarized above, showing that the energy flux from the pulsar has a stronger latitude dependence than previously thought. Early results from 3-D modelling have shown that also our assumptions on the level of magnetization at the shock need to be modified, allowing for much more magnetized winds than 2-D models (see the results and the discussion in the next section). All this wealth of information will have to be incorporated in the models to obtain a reliable assessment of the particle acceleration process/es at work in the different energy ranges.
4 A step forward: 3-D models
The 2-D modelling of PWNe, while very successful at reproducing the morphology of the inner region of these nebulae, is clearly not able to provide at the same time a correct description of the nebular magnetic field on the large scales, as demonstrated by the lack of agreement between the simulated emission spectrum and multiwavelength observations. In 2-D models, the lack of the third degree of freedom suppresses the kink mode, and the magnetic field is allowed to artificially pile up around the symmetry axis. As a consequence, the value of the wind magnetization parameter cannot be raised above
$\unicode[STIX]{x1D70E}\simeq 3\times 10^{-2}$
, otherwise the effect of the compression affects also the evolution of the TS, reducing its size to much smaller values than observed. The small magnetization of the wind that one is forced to assume also leads to an average value of the field in the nebular volume which is lower than that inferred from spectral modelling (
${\sim}200~\unicode[STIX]{x03BC}\text{G}$
) by a factor of
$3\div 4$
. As already discussed, this strongly affects the properties of the nebula at large scales, such as the total integrated spectrum and the emission maps of the entire nebula. The first 3-D models of the Crab nebula via relativistic MHD simulations were presented by Porth, Komissarov & Keppens (Reference Porth, Komissarov and Keppens2013), Porth et al. (Reference Porth, Komissarov and Keppens2014b
). Since these simulations are computationally very demanding, they could only evolve the system for
${\sim}70$
years, so that the self-similar expansion phase was not completely reached. Nevertheless, many interesting results were found. Working in three dimensions allowed the authors to increase the magnetization to values greater than unity, since now the kink instability is free to set in and leads to efficient dissipation of the field inside the nebula, as already noticed in previous simulations of jets (Mignone et al.
Reference Mignone, Striani, Tavani and Ferrari2013). Notice that in these simulations the dissipation is only numerical. They also found that the innermost structure is perfectly analogous to that obtained with 2-D axisymmetric models (of lower magnetization), indicating that MHD models are robust in describing the jet–torus morphology of PWNe. Jets are retrieved in the downstream region of the TS as the result of flow collimation by magnetic hoop stresses. The loss of axisymmetry reflects in a complex structure of the field in the simulated PWN: the highly ordered field no longer survives in the nebula, but it becomes more or less randomized. The toroidal component of the field is still dominant near the TS, while around the polar jets and close to the PWN boundaries (the contact discontinuity with the supernova ejected material) the poloidal component becomes more important.
Despite the important results introduced by this work, some problems are still unsolved. In particular, the value of the magnetic field extrapolated to the actual age of the system appears to be well below what deduced observationally (Porth, private communication). This unfortunately leaves open the question of whether the comparison of emission maps and spectra with real data will be satisfactory. However, it is clear that these kind of 3-D simulations are the primary tool to reach a comprehensive description of PWNe in all their aspects.
We have recently performed our own 3-D simulation of the Crab nebula, with the aim of reaching at least the self-similar expansion phase, if not to follow the complete evolution of the system (
${\sim}1000$
years). Since the simulation by Porth et al. (Reference Porth, Komissarov and Keppens2014b
) shows a faint emission in the torus region, and a complete absence of emission from the polar jets, we have decided to use a slightly different model for the injected pulsar wind, with the main goal of reproducing the observed emission from the Crab nebula, as we discuss below. For ease of comparison with Porth’s results, we maintain some initialization choices, namely at
$t=0$
wind conditions are defined for
$r\leqslant r_{wind}=1$
ly, the slowly expanding supernova ejected material extends up to
$r_{ej}=5$
ly, whereas static interstellar medium (ISM) conditions apply at larger radii, up to the simulation box boundaries. This set-up corresponds to an initial age of the system of
$t_{0}\simeq 250$
years, thus, in order to follow the complete evolution of the Crab nebula, approximately 700 years of simulated evolution would be needed. In the following we will refer to two different times:
$t$
as the simulated time and
$t_{eq}=t+t_{0}$
as the equivalent evolution time of the system.
Our simulation is performed with the shock-capturing, finite-volume numerical code Pluto (Mignone et al.
Reference Mignone, Bodo, Massaglia, Matsakos, Tesileanu, Zanni and Ferrari2007; Mignone et al.
Reference Mignone, Zanni, Tzeferacos, van Straalen, Colella and Bodo2012), which supports adaptive mesh refinement (AMR), and it has been run on 2048 processors at the Tier0 CINECA (Bologna, Italy) FERMI cluster. We define a cubic numerical grid of
$[-10\,\text{ly},+10\,\text{ly}]^{3}$
, that can contain the whole nebula up to the end of its evolution, with the pulsar wind continuously injected in a sphere with radius
$0.01$
ly (which allows the TS to detach from the inner boundary of the wind). Seven AMR levels are employed to guarantee the required accuracy close to the centre of the domain, where the wind is injected. The resolution is then decreased dynamically with distance from the pulsar, preserving a higher resolution around the wind TS and in the inner part of the nebula, more or less in the region in which the jet–torus morphology develops. Since the TS is not static but its position changes with time, we cannot use a static decomposition of the domain, and AMR is a necessary requirement.
The solenoidal condition of the magnetic field (
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}=0$
) is ensured by the divergence-cleaning method, which augments the relativistic MHD system with an extra equation for a generalized Lagrange multiplier, allowing the propagation and progressive damping of the divergence errors (Dedner et al.
Reference Dedner, Kemm, Kröner, Munz, Schnitzer and Wesenberg2002). The high Lorentz factor flow is managed by the selective use of the Harten–Lax–van Leer (HLL) Riemann solver in the most relativistic regions, or the HLLD (where ‘D’ stands for discontinuities) solver (Mignone, Ugliano & Bodo Reference Mignone, Ugliano and Bodo2009) for the smoother part of the system. Additional smoothing in the proximity of the strong shock is achieved by a shock flattening option.
As previously mentioned, the fundamental difference between our model and the one by Porth et al. (Reference Porth, Komissarov and Keppens2014b
) is in the wind parametrization and angular dependence. Following our previous works in two dimensions (see e.g. Olmi et al.
Reference Olmi, Del Zanna, Amato, Bandiera and Bucciantini2014, Reference Olmi, Del Zanna, Amato and Bucciantini2015), we assume an energy flux in the wind which is still axisymmetric and depends on radius
$r$
and polar angle
$\unicode[STIX]{x1D703}$
(measured from the
$z$
symmetry axis, which is assumed to be the rotation axis of the central pulsar) approximately as predicted by the split-monopole model (Michel Reference Michel1973). The observed jet–torus morphology and the expected oblate shape of the TS are accounted for by considering an anisotropic distribution of the energy flux in the wind. The level of anisotropy is governed by an appropriate free parameter of the model (
$\unicode[STIX]{x1D6FC}$
) which represents the ratio between the polar and equatorial energy fluxes in the wind, and which must be chosen
$\unicode[STIX]{x1D6FC}\ll 1$
(Del Zanna et al.
Reference Del Zanna, Amato and Bucciantini2004).
Since the numerical cost of 3-D simulations is very large, we did not include additional numerical tracers in the present simulations, with the exception of the tracer for the particle maximum energy (the
$\unicode[STIX]{x1D716}_{\infty }$
defined in Del Zanna et al. (Reference Del Zanna, Volpi, Amato and Bucciantini2006)), which is absolutely mandatory to account for the synchrotron emission. As a consequence we could only inject particles uniformly along the shock front.
The striped morphology of the magnetic field is ensured by choosing
$B(r,\unicode[STIX]{x1D703})\propto f(b,\unicode[STIX]{x1D703})\sin \unicode[STIX]{x1D703}$
, where
$f(b,\unicode[STIX]{x1D703})$
is an appropriate function of latitude.
$f$
also depends on the free parameter
$b$
, which represents the width of the striped zone of the wind. In particular, for large values of
$b$
, the pure split-monopole configuration is recovered, while in the case of
$b\sim 1$
, corresponding to a large striped region (of approximately
$60^{\circ }$
), due to a highly tilted pulsar magnetosphere, dissipation before the TS is supposed to have shaped the field strength at basically all latitudes (see figure 4 and discussion in Del Zanna et al. (Reference Del Zanna, Volpi, Amato and Bucciantini2006)). The initial magnitude of the magnetic field is governed by the initial magnetization
$\unicode[STIX]{x1D70E}_{0}$
, which is contained in the proportionality constant.
Figure 4. Shape of the wind magnetization
$\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D703})$
for the value of
$\unicode[STIX]{x1D70E}_{0}=1$
employed in our 3-D simulation and for three values of
$b$
. The upper (blue) curve is the case with a narrow striped wind region used here, with
$\unicode[STIX]{x1D70E}_{eff}=1.5$
, while the bottom (red) curve refers to a case that approximately reproduces the magnetization function as in the B3D run of Porth et al. (Reference Porth, Komissarov and Keppens2014b
).
The set of these three free parameters
$(\unicode[STIX]{x1D6FC},b,\unicode[STIX]{x1D70E}_{0})$
is chosen based on the model that appears to best fit the Crab observations in 2-D MHD simulations. In particular we use
$b=10$
and
$\unicode[STIX]{x1D6FC}=0.1$
, and, based on the findings of the previous 3-D study, we raise the magnetization to unity (
$\unicode[STIX]{x1D70E}_{0}=1$
). The initial value of the wind Lorentz factor is chosen to be
$\unicode[STIX]{x1D6FE}_{0}=10$
. While the actual wind Lorentz factor is certainly much larger, say up to
${\sim}10^{6}$
(Kennel & Coroniti Reference Kennel and Coroniti1984a
), the adopted value is high enough to guarantee that the flow is highly relativistic, in which case the post-shock dynamics is known to be independent on the exact value of
$\unicode[STIX]{x1D6FE}_{0}$
.
The chosen parametrization corresponds to an effective magnetization of
$\unicode[STIX]{x1D70E}_{eff}=1.5$
, that can be obtained by integrating over
$\unicode[STIX]{x1D703}$
the quantity
$\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D703})$
defined in equation (5) of Olmi et al. (Reference Olmi, Del Zanna, Amato, Bandiera and Bucciantini2014). Different possible choices of
$\unicode[STIX]{x1D70E}_{eff}$
are shown in figure 4, with our current one being blue in colour. This quantity can be directly compared with the quantity
$\bar{\unicode[STIX]{x1D70E}}$
in table 1 of Porth et al. (Reference Porth, Komissarov and Keppens2014b
), noticing that our parametrization roughly corresponds to their case D3D, in which the pulsar obliquity is of order
$10^{\circ }$
.
In figure 5, the time evolution of the nebula and its TS are shown at the simulation time of 250 yr, corresponding to the physical time
$t_{eq}\simeq 500$
yr. In particular, in the left-hand side of the figure, the ratio between the polar radius of the TS, defined as the maximum extent of the TS in the polar direction, and the radius of the nebula is shown. Assuming a polytropic equation of state with index
$4/3$
and a radial power law for the evolution of the nebula radius, the hydrodynamical evolution can be obtained from the momentum conservation at the TS (see equation (44) from Porth et al. (Reference Porth, Komissarov and Keppens2014b
)). This trend is also drawn as a dashed line in the plots, as a basis for comparison. In the right-hand panel the same plot is produced for the ratio between the equatorial radius of the TS, defined as the maximum extent of the TS in the equatorial plane, and the radius of the nebula. The fits to both ratios with a suitable function is also shown (as a solid line) to help the reader to see the long-term trend within the highly dynamical evolution.
Figure 5. Evolution with time of the ratio between the equatorial/polar radius of the TS and the radius of the nebula. Full circles represent data extracted from the simulation, solid lines our fit to these data and dashed ones the expected behaviour in the hydrodynamical case, as a comparison. In (a,b), the ratio between the polar (equatorial) radius of the TS and the outer radius of the nebula is shown. As can be noticed from these plots, the self-similar phase of evolution appears to be fully reached.
Figure 6. Evolution with time of the magnetic field strength, averaged over the entire nebula volume. Filled circles represents data extracted from simulations, while the line is the fit to the data.
As can be easily understood from these plots, the system seems to have fully reached the self-similar phase of expansion.
In figure 6 the evolution with time of the magnetic field strength averaged over the nebula volume is also shown. The magnetic field appears to decrease down to values of the order of
${\sim}100~\unicode[STIX]{x03BC}\text{G}$
in the first 100 years of simulation, which means that the average field in the nebula is approximately a factor of 2 lower than expected. Anyway the trend seems to indicate that the strong magnetic dissipation of the first period of evolution slows down after around 100 years: the field strength decreases down to a minimum value and afterwards, for the following 100 years, it stays almost constant, or it even increases somewhat. This might indicate that magnetic dissipation progressively becomes less important, leading the magnetic field to saturate to the equipartition value at some point of the evolution. For times up to those investigated by Porth et al. (Reference Porth, Komissarov and Keppens2014b
), the trend we find for the magnetic field time evolution perfectly agrees with their results for the time evolution of the magnetic energy in the nebula.
As expected, after some time, the magnetic topology becomes more complex, with the development of a poloidal component. In figure 7 the structure of magnetic field lines at
$t=250$
yr is shown, with the colour bar indicating the ratio of the toroidal to the total field magnitude, measured by the quantity
$\unicode[STIX]{x1D6FC}_{B}=B_{tor}/B$
. Preponderance of different components of the field are expressed as red/yellow colours and green/blue ones, indicating, respectively, the toroidal and poloidal components. As expected, the field is mainly toroidal in the inner part of the nebula, as only this component is injected in the wind and amplified post-shock, while the poloidal component is mostly important outside along the polar jets and in the external regions, where the field structure is more efficiently modified by polar high-speed flow or by the turbulent motions, respectively. In the same figure, we also show isosurfaces of constant velocity, dark blue indicating the wind region and light blue to highlight the jets, where we have chosen a reference value of
$v=c/2$
.
Figure 7. Magnetic field lines drawn with seed points laying on a sphere of radius 2.5 ly, at the time
$t=250$
yr. The contrast between toroidal and poloidal components (
$\unicode[STIX]{x1D6FC}_{B}=B_{tor}/B$
) is highlighted by the colour scale, with red colour meaning that the field is completely toroidal and the blue one that it is predominantly poloidal. The velocity field is also represented as a 3-D contour plot, with colour levels corresponding to
$0.95c$
and
$0.5c$
. Dimensions of the nebula are indicated by the box axes, in units of ly.
The relative importance of the poloidal and toroidal component of the field can also be seen in a 2-D slice of the data cube, shown in figure 8(a). Here the variable
$\unicode[STIX]{x1D6FD}_{B}\equiv B_{pol}/B$
is represented at
$t=250$
yr in the plane
$(0,y,z)$
, and red colour indicates complete dominance of the poloidal component of the magnetic field (i.e.
$\unicode[STIX]{x1D6FD}_{B}\sim 1$
). As expected, the poloidal field is mostly concentrated near the polar jets. The variable
$\unicode[STIX]{x1D6FD}_{B}$
is also a good marker for estimating the region where results obtained with 2-D axisymmetric models are robust. In fact, where the ratio of the polar component to the total magnetic field is small (blue and green regions, where the polar component is below 50 % of the total) the axisymmetric description (where the field is completely toroidal) is reliable. In yellow and red regions the two descriptions are substantially different. This means that the 2-D description provides a good approximation within approximately
$1/5$
of the nebula radius. Looking at figure 8(b), it is evident that the magnetic field magnitude is stronger in the inner nebula, where the field is mostly toroidal, and where the average value is close to what is expected from spectral modelling
$\left<B\right>\simeq 200\,\unicode[STIX]{x03BC}\text{G}$
.
Figure 8. (a) The ratio of the polar to the total magnetic field magnitude in a 2-D slice of the data cube at
$t=250$
yr. Data are also cropped within a circle of radius equal to the nebular radius (2.8 ly) for ease of interpretation. As expected, the polar component is mostly important near the polar jets. (b) Image of the 2-D distribution of the total magnetic field magnitude in the
$0,y,z$
plane, shown in units of Gauss at the time 250 yr. The magnetic field appears to be stronger close to the TS, where the mean value is around 1 mG.
Figure 9. (a) 2-D map of the particle pressure distribution, with intensity varying according to the colour (logarithmic) scale. (b) Magnitude of the velocity field in a 2-D slice at
$t=250$
yr. (c) Blow-up of the inner nebula in a plot of the velocity magnitude with arrows indicating the velocity field. Effects of the hoop stresses, guiding the flow in the polar direction to form jets, are clearly visible.
Figure 10. (a) Map of the radio surface brightness at 5 GHz (
$t=250$
yr, corresponding to the evolution time
$t_{eq}\simeq 500$
yr). Intensity is indicated by the colour bar, in linear scale. (b) Linear map of the optical surface brightness at
$10^{15}$
Hz (
$t\simeq =250$
yr,
$t_{eq}\simeq 500$
yr). (c) Map of the emission at X-rays (1 keV). In all the maps the intensity is expressed in
$\text{mJy}~\text{arcsec}^{-2}$
units and distance from the pulsar is in ly.
In figure 9 2-D plots of the velocity magnitude and the magnetic pressure at
$t=250$
yr are shown, again in the
$(0,y,z)$
plane. High-velocity fluxes in the polar direction are present, with velocity of the order of
${\gtrsim}0.7c$
. These are also regions of high magnetic pressure. In this regions jets are launched thanks to the magnetic hoop stresses, which collimate the high-velocity flow escaping from the polar sector of the TS in the polar jets. This behaviour can also be easily recognized in the velocity map, shown in (c) of the same figure.
The resulting emission morphology computed in radio, optical and X-ray bands at
$t=250$
yr is shown in figure 10. As can be seen, the synchrotron burn-off effect is clearly apparent at increasing photon energies, and emission maps are encouraging. Polar jets are well formed but the torus is quite under-luminous. Most of the emission comes indeed from the jets bases, where the particle density is higher, while in the torus zone, although the magnetic field is still high, the number of particles is considerably lower (compare with figure 9
a).
The inner nebula also shows the expected time variability, at least at radio frequencies (5 GHz). Since the outputs of our 3-D simulation are very large, and their damping is rather time consuming, we currently only save snapshots every 5 years of evolution. This time interval is too large to address the X-ray variability, which occurs on time scales of weeks to months, but is adequate to discuss variability in the radio band, where times scales are found to be of order a few years (Bietenholz et al. Reference Bietenholz, Hester, Frail and Bartel2004).
In figure 11(b) wisps are shown as a sequence of consecutive light and dark arcs. They form near the pulsar, the position of which is indicated by the black cross, and then move outward in the bulk of the nebula. As expected, wisps appear to be more prominent in the upper hemisphere of the nebula (with respect to the pulsar equatorial plane) than in the lower one. Properties of wisps are here obtained with the same technique as in the 2-D case, by subtracting emission maps separated by 5 years (in radio the expected periodicity of wisps is of the order of years).
Figure 11. Map of the wisp-like structures at radio frequencies (5 GHz) in the inner nebula (a) and in the whole nebula (b). The images are obtained as the difference between two surface brightness maps separated by a time interval of 5 years around
$t=250$
years. The position of the pulsar is identified by the black crosses and intensity is normalized to its maximum.
The results reported in figure 11, when compared with those reported in figure 2 (from 2-D MHD simulations by Olmi et al. (Reference Olmi, Del Zanna, Amato, Bandiera and Bucciantini2014)), immediately highlight the much greater richness of variable structures that is found in 3-D simulations. A number of non-axisymmetric structures appear on the large scales, and the map shows a much closer resemblance to the time variability data reported in figure 2 of Bietenholz et al. (Reference Bietenholz, Hester, Frail and Bartel2004). However, in the inner nebula the 2-D and 3-D results appear to be rather similar, lending support to the validity of the conclusions that were obtained for that region within the 2-D description.
5 The ideal tearing model for the Crab nebula gamma-ray flares
As discussed in § 2.4, the sudden gamma-ray flares observed in the Crab nebula are most likely due to stochastic reconnection events occurring around the TS region or at the jet base, where the strongly magnetized plasma is supposed to convert its energy into heat and particle acceleration. Unless peculiar magnetic configurations leading to forced (driven) reconnection are considered (Lyutikov et al. Reference Lyutikov, Sironi, Komissarov and Porth2016), the most natural process leading to spontaneous reconnection of current sheets is the tearing instability. However, even within relativistic resistive MHD or force-free electrodynamics, as appropriate for the hot, highly magnetized plasma in the Crab nebula, the tearing instability is usually regarded as a slow process (Lyubarsky Reference Lyubarsky2005; Komissarov, Barkov & Lyutikov Reference Komissarov, Barkov and Lyutikov2007), precisely as in classical MHD, with a growth time scale which is the geometric mean of the ideal Alfvénic one and the much longer diffusion time.
This is actually true if one refers to the current sheet width
$a$
as the characteristic length. However, for very thin current sheets, the natural size of the system is rather their length
$L$
, so that, by rescaling the quantities with the (inverse) aspect ratio
$\unicode[STIX]{x1D716}=a/L\ll 1$
, we find a maximum growth rate of the tearing instability modified as:
$$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{max}\,a/c_{A}\sim S_{a}^{-1/2}\Rightarrow \unicode[STIX]{x1D6FE}_{max}\,L/c_{A}\sim S_{L}^{-1/2}\unicode[STIX]{x1D716}^{-3/2},\end{eqnarray}$$
where
$c_{A}$
is the background Alfvén velocity and
$S_{a}$
and
$S_{L}$
are the Lundquist numbers (
$S_{l}=l/c_{A}\unicode[STIX]{x1D702}$
, with
$l$
a characteristic length scale and
$\unicode[STIX]{x1D702}$
the magnetic diffusivity), referred to
$a$
or
$L$
respectively (
$S_{a}=\unicode[STIX]{x1D716}\,S_{L}$
). Suppose now that
$\unicode[STIX]{x1D716}\sim S_{L}^{-\unicode[STIX]{x1D6FC}}$
. For instance, in the steady-state, incompressible, non-relativistic Sweet–Parker model we have
$\unicode[STIX]{x1D6FC}=1/2$
, and the evolution is then expected to be of explosive type. Provided
$S_{L}$
is large enough, the current sheet is seen to be rapidly disrupted by the so-called plasmoid instability (Loureiro, Schekochihin & Cowley Reference Loureiro, Schekochihin and Cowley2007; Bhattacharjee et al.
Reference Bhattacharjee, Huang, Yang and Rogers2009; Samtaney et al.
Reference Samtaney, Loureiro, Uzdensky, Schekochihin and Cowley2009; Uzdensky, Loureiro & Schekochihin Reference Uzdensky, Loureiro and Schekochihin2010; Sironi et al.
Reference Sironi, Giannios and Petropoulou2016). In any realistic environment, where current sheets may form as a result of convective motions, occurring on ideal or Alfvèn time scales, it is natural to expect their progressive thinning. As a result, the system becomes tearing and plasmoid unstable before the Sweet–Parker configuration is realized. The maximum growth rate
$\unicode[STIX]{x1D6FE}_{max}$
in equation (5.1) becomes independent on
$S_{L}$
, and we can speak of an ideal tearing instability (Pucci & Velli Reference Pucci and Velli2014; Landi et al.
Reference Landi, Del Zanna, Papini, Pucci and Velli2015; Tenerani et al.
Reference Tenerani, Velli, Rappazzo and Pucci2015; Del Zanna et al.
Reference Del Zanna, Landi, Papini, Pucci and Velli2016a
).
The theory of the tearing instability has been recently revisited, for the first time, in the case of resistive relativistic MHD, and both analytical and numerical investigation of the classical and ideal tearing have been carried out (Del Zanna et al.
Reference Del Zanna, Papini, Landi, Bugli and Bucciantini2016b
). In particular, when current sheets with
$\unicode[STIX]{x1D716}\sim S_{L}^{-1/3}$
are considered, the expected instability dispersion relation, independent on the value of
$S_{L}$
, is retrieved, and both the linear and nonlinear stages are seen to follow a universal law, if quantities are normalized against the relativistic Alfvén speed (here
$4\unicode[STIX]{x03C0}\rightarrow 1$
)
$$\begin{eqnarray}c_{A}=\frac{B}{\sqrt{e+p+B^{2}}}\,c,\end{eqnarray}$$
where
$B$
is the magnetic strength,
$e$
the energy density and
$p$
the thermal pressure, for the unperturbed fluid (
$c_{A}\rightarrow c$
for magnetically dominated plasmas).
Figure 12. (a) The growth in time of root-mean-square magnetic fluctuations from the initial seed value for different runs, symbols and colours refer to different values of the background Alfvén velocity
$c_{A}$
. (b) Colour map of the flow velocity and magnetic field lines for the case
$c_{A}=0.5c$
at time
$t=20\,\unicode[STIX]{x1D70F}_{c}\equiv 10\,\unicode[STIX]{x1D70F}_{A}$
.
In figure 12(a) the averaged magnetic fluctuations across the current sheet (initially set to a seed small value) are shown, as a function of
$t/\unicode[STIX]{x1D70F}_{A}\equiv t\,c_{A}/L$
, for five different numerical simulations at increasing
$c_{A}/c$
(from 0.50 up to 0.98): notice the
$e$
-folding time of order
$\unicode[STIX]{x1D70F}_{A}$
(thus an ideal time scale), and a very similar trend even for the nonlinear stage when the plasmoid instability develops, with the merging of magnetic islands and secondary reconnection events. The final stage of the evolution is, due to the periodic boundary conditions assumed in the sheet direction, a single X-point with Petschek-type funnels containing super-(magneto)sonic jets feeding a large magnetic island. The maximum velocity achieved in jets corresponds to the external Alfvén speed,
$c_{A}=0.5\,c$
in the case shown in figure 12(b), and we find a quite high reconnection rate which depends only mildly on the Lundquist number,
${\mathcal{R}}\sim \ln (S_{L})^{-1}$
, as predicted in the relativistic case for magnetically dominated plasmas (Lyubarsky Reference Lyubarsky2005).
The result that the ideal tearing instability develops within a few light-crossing times of the macroscopic length
$L$
, independently on microphysics, is of course very important in our attempt to model the gamma-ray flares of the Crab nebula. The measured time scale
$\unicode[STIX]{x1D70F}$
for growth in the light curve is of approximately one day, so we need current sheets to develop on a scale
$L\sim \unicode[STIX]{x1D70F}c\sim 10^{15}$
cm. This length is
${\sim}1\,\%$
of the TS radius, compatible with the size of the small-scale variable features observed in the inner regions by HST and Chandra, like the knot or the anvil, but also with the base of the polar jets, whose kinks are likely to produce sheared magnetic regions. We cannot certainly determine the precise location for reconnection in the Crab nebula from this highly idealized model, but we have demonstrated that the reconnection process is able to occur on ideal time scales, independently on
$S_{L}$
and without invoking kinetic mechanisms, provided that during the current sheet dynamical evolution (thinning) the critical threshold of
$\unicode[STIX]{x1D716}\sim (S_{L})^{-1/3}$
can be reached.
6 Summary and conclusions
In the last decade, MHD numerical modelling has been widely recognized as a powerful tool for the investigation of the physics of the relativistic magnetized plasma in PWNe.
While many open questions still remain, the advent of relativistic MHD simulations has allowed us to answer a number of questions, and has opened a promising way to answer more. The initially puzzling jet–torus structures observed in the X-ray emission of a number of nebulae have found a straightforward interpretation as the result of the pulsar wind structure and in particular of the anisotropy of its energy flux. The striking agreement between the simulated and observed high-energy morphology of the Crab nebula, which is the PWNe prototype, suggests indeed that the flow structure in the inner regions must be very close to what 2-D MHD simulations predict (Del Zanna et al. Reference Del Zanna, Volpi, Amato and Bucciantini2006).
The open problems concern the exact origin of the low-energy emitting particles, that can be either continuously injected in the nebula as part of the pulsar wind or a relic population, or finally the result of acceleration somewhere else.
This problem was approached by comparing observations with simulations done under different hypotheses for the radio particles origin. The main conclusion was that the global radio emission is basically insensitive to the spatial distribution of the particles, as long as they are not pure relics. The latter is the only case that can be easily excluded, not only the surface brightness maps are almost identical, but also the variability of the inner regions (radio wisps) is very similar too. This is a result of the fact that within the MHD framework wisps naturally arise from the properties of the flow. This result should be validated in the future also in three dimensions, since the stronger mixing of the magnetic field lines could in principle allow for the relic scenario to be reconsidered.
It is a matter of fact that wisps are observed at radio, optical and X-ray frequencies, and that they are not coincident at the different wavelengths, with differences measured both in terms of spatial locations and of outward velocities (Bietenholz et al. Reference Bietenholz, Hester, Frail and Bartel2004; Schweizer et al. Reference Schweizer, Bucciantini, Idec, Nilsson, Tennant, Weisskopf and Zanin2013). As already discussed, within the framework of MHD, the only way to account for the different properties of the wisps at different frequencies is to assume some differences in the spatial distribution of the particles that are responsible for the emission at the different wavelengths. An obvious difference among particles of different energies (and hence emitting in different frequency bands) is the role of energy losses. However, losses are not very important in the innermost regions of the nebula from which wisps are observed. Any other difference in the spatial distribution of the particles of different energies suggests differences in the acceleration sites. In other words, if wisps are different at radio and X-ray frequencies, the particles responsible for radio emission must have had a different history than the particles responsible for X-ray emission, and the simplest explanation is that they were accelerated in different locations.
This possibility has been studied by performing 2-D MHD simulations which assume a non-uniform particle acceleration at the TS front. The TS was divided in complementary regions (an equatorial band and a polar cone, with varying angular extent) and particles of different energies were injected in non-coincident ones. The conclusion is that X-ray wisps are best reproduced if injection in a narrow belt around the pulsar rotational equator is considered. The equatorial region is presumably where the flow magnetization is lower, if effective magnetic dissipation takes place in the striped wind. Then the mechanism responsible for the acceleration of those highest-energy particles could be Fermi I, which has been proven to be very effective at relativistic shocks of low enough magnetization (Spitkovsky Reference Spitkovsky2008; Sironi & Spitkovsky Reference Sironi and Spitkovsky2011), and which naturally provides a particle spectral index close to that inferred from X-ray observations.
Radio emission properties are again much more complicated to constrain: the only scenario that radio observations directly exclude is one in which low-energy particles are injected in a narrow polar cone, since almost no wisps are found. But emission properties in the case of uniform injection, or injection in a wide polar or equatorial band, are basically indistinguishable.
Since 2010, when the Crab gamma-ray flares were revealed for the firs time (Tavani et al. Reference Tavani2011), many attempt to explain this explosive events were done. Up to now the most plausible explanation seems to be violent stochastic reconnection events occurring in the TS region or at the jet base.
This can be obtained thanks to either collapsing magnetic flux bundles or other peculiar configurations, leading to explosive reconnection and particle acceleration (Lyutikov et al.
Reference Lyutikov, Sironi, Komissarov and Porth2016), or via the spontaneous reconnection in thinning current sheets, the so-called ideal tearing instability, which is also proved to develop on extremely rapid time scales, basically the light-crossing time of the sheet length. Del Zanna et al. (Reference Del Zanna, Papini, Landi, Bugli and Bucciantini2016b
) have shown that this mechanism can account for the correct time variation scaling observed in the light curves of gamma-ray flares by considering current sheets which develop on scales of
${\sim}1\,\%$
of the TS radius, which are compatible with the scales of the bright variable features observed in the inner nebula and at the base of the polar jets.
Although axisymmetric 2-D models have allowed us to investigate many of the problems connected to the physics at work in the inner nebula, they were not able to reproduce the bulk properties of the nebula at larger scales. In fact they point to values of the wind magnetization of the order of
$10^{-2}$
, which are incompatible with the integrated emission spectrum, since they lead to a nebular magnetic field which is too low to account for the whole emission. As already highlighted by Volpi et al. (Reference Volpi, Del Zanna, Amato and Bucciantini2008), reproducing the synchrotron part of the spectrum requires conversion of the flow energy flux into accelerated particles with efficiencies larger than one, over-predicting the IC emission. Moreover the suppression of the third dimension also produces an artificial pile up of the magnetic field morphology, which prevents to compute reliable emission properties of the outer nebula, and also to reproduce the kinking jets.
The first 3-D modelling seems to provide a solution to this problem (Porth et al. Reference Porth, Komissarov and Keppens2013, Reference Porth, Komissarov and Keppens2014b ). These models show that kink instabilities, which are artificially suppressed in two dimensions, efficiently reduce the hoop stresses that force the magnetization to have such low values, allowing reaching values of the order of unity. Recent works about current driven instabilities in three dimensions (Mizuno et al. Reference Mizuno, Lyubarsky, Nishikawa and Hardee2011; Mignone et al. Reference Mignone, Striani, Tavani and Ferrari2013) confirm in fact the formation of kinks in jets, even if the flow arising from the termination shock and the hoop stresses in the nebula itself are not taken into account.
These first 3-D simulations seem to provide a reasonable solution to the long-standing
$\unicode[STIX]{x1D70E}$
-paradox and to account for the formation of a PWN with features similar to those observed in the Crab nebula, although the jets are very weak and almost under-luminous. Moreover they only cover a very short time (
${\sim}70$
yr) after the SN explosion, and at the end of the simulation the self-similar expansion phase has not been reached yet. Here we have presented our preliminary results of the first long-term 3-D simulations of the Crab nebula. We show that, in just
${\sim}200$
yrs the nebula reaches the self-similar phase, and thus the inferred information should be directly compared with observed properties.
As also shown by Porth et al. (Reference Porth, Komissarov and Keppens2014b
), the onset of the kink instability leads the initial toroidal field to be strongly modified during the evolution of the system, with the poloidal component which suddenly develops and becomes as important as the toroidal one near jets. As soon as the self-similar phase has been reached the strong magnetic dissipation, which lowers the initial field to average values of approximately
$100~\unicode[STIX]{x03BC}\text{G}$
, slows down. Extrapolation of this trend is promising to lead the field to reach the equipartition value at the final stage of the evolution. Nevertheless the volume averaged field is still lower than expected (by a factor of
${\sim}2$
), which may indicate that the dissipation is still too high in the first stages of the evolution.
Our simulations show that the dynamics of the inner nebula is in good agreement with what was found in 2-D simulations, suggesting that 2-D models are absolutely robust in describing the physics of the inner regions. This is obviously not true anymore when one moves from the pulsar vicinity to the outer nebula, where the effects of the poloidal component of the field become important.
We compute emission properties at radio, optical and X-ray frequencies. The synchrotron burn-off effect is seen through the decrease of the size of the emitting area with increasing observation frequency. However the brightness contrast of the torus is rather low when compared with observations, while the brightest features appears to lie at the bases of the jets, where the particle pressure is higher.
We also compute maps of the time variability at radio frequencies for both the wisps region (the inner nebula) and for the entire nebula, obtaining very encouraging results. The simulated radio wisps confirm again the similarities between results from 2-D and 3-D simulations in reproducing the properties of the inner nebula. In the outer regions, on the other hand, the 3-D results resemble the observations much more closely than those in two dimensions, with the appearance of moving structures near the polar jets and a number of additional arc-like features in the body of the nebula.
As soon as the simulation will have reached an evolved enough stage, we plan to address the issue of the integrated synchrotron and IC spectrum, and gather insight on the actual number of particles in the wind. The goal of constraining the pulsar wind physics, by finding a suitable set of parameters that allows to reproduce both the morphology and multi-wavelength spectrum of the Crab Nebula, finally seems within reach.
Acknowledgements
B.O. and L.D.Z. acknowledge support from the INFN – TEONGRAV initiative (local PI: L.D.Z.) and from the University of Florence grant ‘Fisica dei plasmi relativistici: teoria e applicazioni moderne’. N.B. acknowledges support from the NSMAG (EU FP7 - CIG) project (PI: N.B.).
