Relaxation of Pulsar Wind Nebula via Current-Driven Kink Instability

1
Relaxation of Pulsar Wind Nebula via
Current-Driven Kink Instability

• Yosuke Mizuno (?? ??)
• Institute of Astronomy
• National Tsing-Hua University
• Collaborators
• Y. Lyubarsky (Ben-Gurion Univ), K.-I. Nishikawa
  (NSSTC/UAH), P. E. Hardee (Univ. of Alabama,
  Tuscaloosa)

Mizuno et al., 2011, ApJ, 728, 90
2
Pulsar Wind Nebulae
Pulsar magnetosphere
Termination Shock
Pulsar wind
Pulsar wind nebula
electromagnetic fields
Synchrotron IC radiation

• Pulsar wind nebulae (PWNe) are considered as
  relativistically hot bubbles continuously pumped
  by e-e- plasma and magnetic field emanating from
  pulsar
• Pulsar loses rotation energy by generating
  highly magnetized ultra-relativistic wind
• Pulsar wind terminates at a strong reverse shock
  (termination shock) and shocked plasma inflates a
  bubble with in external medium
• From shocked plasma Synchrotron and
  Inverse-Compton radiation are observed from radio
  to gamma-ray band (e.g., Gaensler Slane 2006)

3
Pulsar Wind Nebulae (obs.)
Vela (Pavlov et al. 2001)
3C58 (Slane et al. 2004)
G54.10.3 (Lu et al. 2002)
G320.4-1.2 (Gaensler et al. 2002)
4
Simple Spherical Model of PWNe

• Close to pulsar, energy is carried mostly by
  electromagnetic fields as Poynting flux
• Common belief at termination shock, wind must
  already be very weakly magnetized
• Magnetization parameter s (ratio of Poynting to
  kinetic energy flux) needs to be as small as
  0.001-0.01 at termination shock (e.g., Rees
  Gunn 1974, Kennel Coroniti 1984)
• Such low value of s is puzzling because it is not
  easy to invent a realistic energy conversion
  mechanism to reduce s to required level (s
  problem) (reviews by Arons 2007 Kirk et al.
  2009)

5
Dependence on s to shock downstream structure
Kennel Coroniti 1984
Postshock speed
At shock downstream
c/3
sgtgt1 effectively weak (magnetic energy
dominated) sltlt1 significant fraction of total
energy in upstream converted to thermal energy in
downstream
sgtgt1 almost constant with relativistic
speed sltlt1 velocity just after shock becomes c/3
limit, then decreasing From radio observation of
Crab nebula, expanding velocity is 2000km/s at
2pc (s0.003)
6
Axisymmetric RMHD Simulations of PWNe
Del Zanna et al.( 2004)
Synchrotron emission map
Flow magnitude

• Extensive axisymmetric RMHD simulations of PWNe
  show that the morphology of PWNe including
  jet-torus structure with s0.01(e.g., Komissarov
  Lyubarsky 2003, 2004, Del Zanna et al. 2004,
  2006)
• If magnetization were larger, then the nebula
  would be elongated by magnetic pinch effect
  beyond observational limits

7
Termination Shock structure
Axisymmetric RMHD simulations of PWNe Komissarov
Lyubarsky 2003, 2004 Del Zanna et al. 2004,
2006 Bogovalov et al. 2005
Del Zanna et al. 2004
Flow magnitude
F?sin2(?) ??sin2(?) B??sin(?)G(?)
A ultrarelativistic Pulsar wind B subsonic
equatorial outflow C supersonic equatorial
funnel D bright arch a termination shock
front b rim shock c FMS surface
8
Constraining ? in PWNe
Smaller s, jet does not formed
?0.03
?0.003
Larger s, PWNe elongates
?gt0.01 required for Jet formation (a factor of
10 larger than within 1D spherical MHD models)
?0.01
(Del Zanna et al. 2004)
9
Dependence on Field Structure
?0.03
b100
b10
B(?)
(Del Zanna et al. 2004)
10
Synchrotron Emission maps
X-rays
optical
?0.025, b10
(Weisskopf et al 00)
(Hester et al 95)
?0.1, b1
Emax is evolved with the flow f(E)?E-?, EltEmax
(Del Zanna et al. 2006)
(Pavlov et al 01)
11
Obliquely rotating Pulsar magnetosphere

• In pulsar wind, most of energy transferred by
  waves, which an obliquely rotating magnetosphere
  excites near the light cylinder
• In equatorial belt of wind, the sign of magnetic
  field alternates with pulsar period, forming
  stripes of opposite magnetic polarity (striped
  wind Michel 1971, Bogovalov 1999)
• Theoretical Modeling of pulsar wind suggest
  that most of wind energy is transported in
  equatorial belt (Bogovalov 1999 Spitkovsky 2006)
• In the equatorial belt, magnetic dissipation of
  the striped wind would be a main energy
  conversion mechanism

Spitkovsky (2006)
12
Obliquely rotating Pulsar magnetosphere
(Kirk Lyubarsky 01)

• In pulsar wind, most of energy transferred by
  waves, which an obliquely rotating magnetosphere
  excites near the light cylinder
• In equatorial belt of wind, the sign of magnetic
  field alternates with pulsar period, forming
  stripes of opposite magnetic polarity (striped
  wind Michel 1971, Bogovalov 1999)

13
Obliquely rotating Pulsar magnetosphere (cont.)
Spitkovsky (2006)

• Theoretical Modeling of pulsar wind suggest that
  most of wind energy is transported in equatorial
  belt (Bogovalov 1999 Spitkovsky 2006)
• In the equatorial belt, magnetic dissipation of
  the striped wind is main energy conversion
  mechanism

14
Dissipation of Alternating Fields

• For simple wave decay, due to relativistic time
  dilation, complete dissipation could occur only
  on a scale comparable to or larger than radius of
  termination shock (Lyubarsky Kirk 2001 Kirk
  Skjaeraasen 2003)
• But, alternating fields can annihilate at
  termination shock by strong deceleration of wind
  via magnetic reconnection (Petri Lyuabrsky
  2007)
• After waves decay via magnetic reconnection s
  lt 1 (0.1)
• At quantitative level, s problem is partially
  solved if Poynting flux is converted into plasma
  energy via dissipation of oscillating part of
  field

1D RPIC simulation with s 45, G 20
(dissipation occurs)
Petri Lyubarsky 2007
15
Dissipation of Alternating Fields in Far Zone of
Wind
Waves decay
wind accelerates
dissipation rate ?
proper wavelength
proper time
The wave dissipation scale is about or larger
than the termination shock radius (Lyubarsky
Kirk 2001 Kirk Skjæraasen 2003)
The flow sharply decelerates at the shock
dissipation of alternating fields at the
termination shock
16
Magnetic Reconnection at Termination Shock of
striped pulsar wind
(l wavelength of striped wind, g1 Lorentz
factor at upstream)
Full dissipation at
dissipation with s 45, g1 20
Initial condition
1D RPIC Simulation (Pétri Lyubarsky, AA, 2007)
17
Energy conversion at high latitudes

• At high latitudes, magnetic field does not change
  sign (no reconnection occurs)
• Fast magnetosonic waves may transport significant
  amount of energy
• These waves can decay relatively easily
  (Lyubarsky 2003) but can release only a fraction
  of the Poynting flux into plasma (because at high
  latitudes, most of energy is carried by mean
  magnetic field)
• Even though this fraction is still not known,
  this fraction is less than ½ because angular
  distribution of Poynting flux in pulsar wind is
  maximum at rotational equator, where mean field
  is zero

18
Another Possibility CD Kink Instability in PWNe

• At quantitative level, s problem is partially
  solved if Poynting flux is converted into plasma
  energy via dissipation of oscillating part of
  field (Petri Lyubarsky 2007)
• But, from residual magnetic field, s still
  cannot be as small as required (0.11).
• Question still remains how the residual mean
  field s could become extremely small
  (0.0010.01) need another mechanism
• Begelman (1998) proposed that problem can be
  solved if current-driven kink instability
  destroys concentric field structure in pulsar
  wind nebula
• As first step, we perform 3D evolution of simple
  cylindrical model of PWNe (Begelman Li 1992)
  with growing CD kink instability using 3D RMHD
  simulation code

19
CD Kink Instability

• Well-known instability in laboratory plasma
  (TOKAMAK), astrophysical plasma (Sun, jet, pulsar
  etc).
• In configurations with strong toroidal magnetic
  fields, current-driven (CD) kink mode (m1) is
  unstable.
• This instability excites large-scale helical
  motions that can be strongly distort or even
  disrupt the system
• For static cylindrical force-free equilibria,
  well known Kurskal-Shafranov (KS) criterion
• Unstable wavelengths
• l gt Bp/Bf 2pR
• However, rotation and shear motion could
  significant affect the instability criterion

Schematic picture of CD kink instability
3D RMHD simulation of CD kink instability in
helical force-free field (Mizuno et al. 2009)
20
Purpose of Study

• Begelman (1998) proposed that s problem can be
  solved if current-driven kink instability
  destroys concentric field structure in pulsar
  wind nebula
• As first step, we perform 3D evolution of simple
  cylindrical model of PWNe (Begelman Li 1992)
  with growing CD kink instability using 3D RMHD
  simulation code RAISHIN

21
4D General Relativistic MHD Equation

• General relativistic equation of conservation
  laws and Maxwell equations
• ?n ( r U n ) 0
  (conservation law of particle-number)
• ?n T mn 0 (conservation
  law of energy-momentum)
• ?mFnl ?nFlm ?lF mn 0
• ?mF mn - J n
• Ideal MHD condition FnmUn 0
• metric ds2-a2 dt2gij (dxib i dt)(dx jb j dt)
• Equation of state p(G-1) u

(Maxwell equations)
r rest-mass density. p proper gas pressure.
u internal energy. c speed of light. h
specific enthalpy, h 1 u p / r. G specific
heat ratio. Umu velocity four vector. Jmu
current density four vector. ?mn covariant
derivative. gmn 4-metric. a lapse function,
bi shift vector, gij 3-metric Tmn energy
momentum tensor, Tmn pgmn r h Um UnFmsFns
-gmnFlkFlk/4. Fmn field-strength tensor,
22
Conservative Form of GRMHD Equations (31 Form)
(Particle number conservation)
(Momentum conservation)
(Energy conservation)
(Induction equation)
U (conserved variables)
Fi (numerical flux)
S (source term)
v-g determinant of 4-metric vg determinant of
3-metric
Detail of derivation of GRMHD equations Anton et
al. (2005) etc.
23
3D GRMHD code RAISHIN
Mizuno et al. 2006a, astro-ph/0609004 Mizuno et
al. 2011, ApJ

• RAISHIN dode utilizes conservative,
  high-resolution shock capturing schemes
  (Godunov-type scheme) to solve the 3D General
  Relativistic MHD equations (metric is static)
• Reconstruction PLM (Minmod MC
  slope-limiter), CENO, PPM, WENO, MP, MPWENO,
  WENO-Z, WENO-M, Lim03
• Riemann solver HLL, HLLC, HLLD approximate
  Riemann solver
• Constrained Transport Flux CT, Fixed Flux-CT,
  Upwind Flux-CT
• Time evolution Multi-step TVD Runge-Kutta
  method (2nd 3rd-order)
• Recovery step Noble 2 variable method,
  Mignore-McKinney 1 variable method
• Equation of states constant G-law EoS,
  variable EoS for ideal gas

Numerical Schemes
24
Ability of RAISHIN code

• Multi-dimension (1D, 2D, 3D)
• Special and General relativity (static metric)
• Different coordinates (RMHD Cartesian,
  Cylindrical, Spherical and GRMHD Boyer-Lindquist
  of non-rotating or rotating BH)
• Different spatial reconstruction algorithms (10)
• Different approximate Riemann solver (3)
• Different constrained transport schemes (3)
• Different time advance algorithms (2)
• Different recovery schemes (2)
• Using constant G-law and variable Equation of
  State (Synge-type)
• Parallelized by OpenMP (shared memory) and MPI
  (distributed memory)

25
Relativistic Regime

• Kinetic energy gtgt rest-mass energy
• Fluid velocity light speed
• Lorentz factor gtgt 1
• Relativistic jets/ejecta/wind/blast waves
  (shocks) in AGNs, GRBs, Pulsars, etc
• Thermal energy gtgt rest-mass energy
• Plasma temperature gtgt ion rest mass energy
• p/r c2 kBT/mc2 gtgt 1
• GRBs, magnetar flare?, Pulsar wind nebulae
• Magnetic energy gtgt rest-mass energy
• Magnetization parameter sgtgt 1
• s Poyniting to kinetic energy ratio B2/4pr
  c2g2
• Pulsars magnetosphere, Magnetars

26
Cylindrical Model of PWNe

• This model (Begelman Li 1992) quasi-static
  cylindrical configuration with purely toroidal
  magnetic field
• The plasma within cylinder is relativistically
  hot and hoop stress is balanced by thermal
  pressure
• Cylinder is confined on outside by non-magnetized
  plasma
• Linear analysis shows that such configuration is
  unstable with respect to CD kink instability
  (Begelman 1998)

27
Initial Condition for Simulations
Radial profile
pressure
Toroidal field

• We solve 3D RMHD equations in Cartesian
  coordinates
• We consider hydrostatic hot plasma column
  containing a pure toroidal magnetic field with
  radius R and height Lz (magnetic hoop stress is
  balanced by gas pressure)
• At Rgt1, hot plasma column is surrounded by a hot
  static unmagnetized medium with constant gas
  pressure
• p0105 r0c2 (relativistically hot, rc2 ltlt pg),
  G4/3 (adiabatic index)
• Put small radial velocity perturbation
• Computational domain Cartesian box of size 6R x
  6R x Lz (Lz1R) with grid resolution of N/R,L60
• Boundary periodic in axis direction, reflecting
  boundary in x, y direction

N total number of modes, fk random phase,
akx,y, random direction
28
Cylindrical Model of PWNe
Based on cylindrical model of PWNe (Begelman Li
1992), radial gas pressure and toroidal magnetic
field profiles in hot plasma column are given by
Where xr/R, h is found for any x from equation
In this solution, magnetic hoop stress is
balanced by gas pressure
At xgt1, hot plasma column is surrounded by hot
static unmagnetized medium with constant gas
pressure
where h0 is solution of eqs at x1
29
Results (2D gas prssure)
Case A perturbation N2, fk0, n1 mode in
x-direction, n2 mode in y-direction
Gas pressure

• Initial small velocity perturbation excites CD
  kink instability n1 mode in x-direction and n2
  mode in y-direction
• radial velocity increases with time in linear
  growth phase
• At about t6R/c, CD kink instability shifts to
  nonlinear phase
• In nonlinear phase, two modes interact and lead
  to turbulence in hot plasma column
• Gas pressure within column, which was initially
  high to balance magnetic hoop stress, decreases
  because hoop stress weakens

30
Results (2D magnetic field)
Case A perturbation N2, fk0, n1 mode in
x-direction, n2 mode in y-direction
As a result of CD kink instability, magnetic
loops come apart and release magnetic stress
31
Time Evolution of Volume Averaged Quantities
Eprhg2-p, EmB2/2, EtEpEm

• Initial slow evolution in linear growth phase
  lasts up to t6R/c, and is followed by a more
  rapid evolution in nonlinear growth phase
• In nonlinear phase, rapid decrease of magnetic
  energy ceases about t11R/c
• While magnetic energy declines, plasma energy
  increases because growth of CD kink instability
  leads to radial velocity increases which
  contributes kinetic energy

magnetic energy
Plasma energy
Total energy

• At about t11R/c, increase in plasma energy
  nearly ceases and hot plasma column is almost
  relaxed
• Multiple-mode (dashed lines) lead to more gradual
  interaction, slower development of turbulent
  structure, and later relaxation of hot plasma
  column

32
Time Evolution of s
Volume-averaged magnetization parameter s in hot
plasma column (Rlt1)
sB2/rh (for hot plasma definition)

• Initially, volume-averaged magnetization s 0.3
  in hot plasma column
• In linear growth phase, s gradually decreases
• After transition to nonlinear phase, s rapidly
  decreases because the magnetic field strongly
  dissipates by the turbulent motion
• When CD kink instability saturates, s0.01

33
Radial Profile
Case A
Radial profile of toroidal- and axial- averaged
quantities for case A
Radial field
Toroidal field

• In linear phase, Br Bz grow, while Bf pg
  decline gradually beginning from near the axis
• In nonlinear phase, Bf pg decrease rapidly,
  and Br Bz increase throughout hot plasma column
• At end of nonlinear phase (t11R/c), all
  magnetic field components become comparable and
  field totally chaotic
• In saturation phase, magnetized column begins
  slow radial expansion (relaxation)

Gas pressure
Axial field

• For different initial perturbation profiles,
  evolutionary timescale is different but physical
  behavior is similar (not shown here)

34
Discussion Elongation of PWNe

• Our simulation confirm scenario envisaged by
  Begelman (1998)
• Toroidal magnetic loops come apart, hoop stress
  declines, and pressure difference across the
  nebula is washed out in nonlinear phase of CD
  kink instability
• For this reason, elongation of PWNe cannot be
  correctly estimated by axisymmetric models
• Because axisymmetric models retain a concentric
  toroidal magnetic field geometry
• To understand the morphology of PWNe correctly,
  we should perform 3D RMHD simulations

35
Discussion Radiation

• Radiation from Crab nebula is highly polarized
  along axis of nebula (e.g., Michel et al. 1991,
  Fesen et al. 1992)
• It is indicated that the existence of ordered
  toroidal magnetic field in PWNe
• From our simulation results, we see that even
  though instability eventually destroys toroidal
  magnetic field structure, magnetic field becomes
  completely chaotic only at the end of nonlinear
  stage of development
• Therefore toroidal magnetic field should
  dominate in central part of nebula that are
  filled by newly injected plasma

36
Summery

• We have investigated development of CD kink
  instability of a hydrostatic hot plasma column
  containing toroidal magnetic field as a model of
  PWNe
• CD kink instability is excited by a small
  initial velocity perturbation and turbulent
  structure develops inside the hot plasma column
• At end of nonlinear phase, hot plasma column
  relaxes with a slow radial expansion
• Magnetization s decreases from initial valule
  0.3 to 0.01
• For different initial perturbation profiles,
  timescale is a bit different but physical
  behavior is same
• Therefore relaxation of a hot plasma column is
  independent of initial perturbation profile
• Our simulation confirm the scenario envisaged by
  Begelman (1998)

37
Crab Nebula