The Role of The Equation of State in Resistive Relativistic Magnetohydrodynamics

1
The Role of The Equation of State in Resistive
Relativistic Magnetohydrodynamics

• Yosuke Mizuno
• Institute of Astronomy
• National Tsing-Hua University

Mizuno 2013, ApJS, 205, 7
ASIAA CompAS Seminar, March 19, 2013
2
Contents

• Introduction Relativistic Objects, Magnetic
  reconnection
• Difference between Ideal RMHD and resistive RMHD
• How to solve RRMHD equations numerically
• General Equations of States
• Test simulation results (code capability in
  RRMHD, effect of EoS)
• Summery

3
Relativistic Regime

• Kinetic energy gtgt rest-mass energy
• Fluid velocity light speed
• Lorentz factor ggtgt 1
• Relativistic jets/ejecta/wind/blast waves
  (shocks) in AGNs, GRBs, Pulsars
• 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 s gtgt 1
• s Poyniting to kinetic energy ratio B2/4pr
  c2g2
• Pulsars magnetosphere, Magnetars
• Gravitational energy gtgt rest-mass energy
• GMm/rmc2 rg/r gt 1
• Black hole, Neutron star
• Radiation energy gtgt rest-mass energy
• Er /rc2 gtgt1
• Supercritical accretion flow

4
Relativistic Jets
Radio observation of M87 jet

• Relativistic jets outflow of highly collimated
  plasma
• Microquasars, Active Galactic Nuclei, Gamma-Ray
  Bursts, Jet velocity c
• Generic systems Compact object(White Dwarf,
  Neutron Star, Black Hole) Accretion Disk
• Key Issues of Relativistic Jets
• Acceleration Collimation
• Propagation Stability
• Modeling for Jet Production
• Magnetohydrodynamics (MHD)
• Relativity (SR or GR)
• Modeling of Jet Emission
• Particle Acceleration
• Radiation mechanism

5
Relativistic Jets in Universe
Mirabel Rodoriguez 1998
6
Ultra-Fast TeV Flare in Blazars

• Ultra-Fast TeV flares are observed in some
  Blazars.
• Vary on timescale as sort as
• tv3min ltlt Rs/c 3M9 hour
• For the TeV emission to escape pair creation
  Gemgt50 is required (Begelman, Fabian Rees 2008)
• But PKS 2155-304, Mrk 501 show moderately
  superluminal ejections (vapp several c)
• Emitter must be compact and extremely fast
• Model for the Fast TeV flaring
• Internal Magnetic Reconnection inside jet
  (Giannios et al. 2009)
• External Recollimation shock (Bromberg
  Levinson 2009)

PKS2155-304 (Aharonian et al. 2007) See also
Mrk501, PKS122221
Giannios et al.(2009)
7
Magnetic Reconnection in Relativistic
Astrophysical Objects

• Pulsar Magnetosphere Striped pulsar wind
• obliquely rotating magnetosphere forms stripes
  of opposite magnetic polarity in equatorial belt
• magnetic dissipation via magnetic reconnection
  would be main energy conversion mechanism

• Magnetar Flares
• May be triggered by magnetic reconnection at
  equatorial current sheet

Spitkovsky (2006)
8
Purpose of Study

• Quite often numerical simulations using ideal
  RMHD exhibit violent magnetic reconnection.
• The magnetic reconnection observed in ideal RMHD
  simulations is due to purely numerical
  resistivity, occurs as a result of truncation
  errors
• Fully depends on the numerical scheme and
  resolution.
• Therefore, to allow the control of magnetic
  reconnection according to a physical model of
  resistivity, numerical codes solving the
  resistive RMHD (RRMHD) equations are highly
  desirable.
• We have newly developed RRMHD code and
  investigated the role of the equation of state in
  RRMHD regime.

9
Applicability of Hydrodynamic Approximation

• To apply hydrodynamic approximation, we need the
  condition
• Spatial scale gtgt mean free path
• Time scale gtgt collision time
• These are not necessarily satisfied in many
  astrophysical plasmas
• E.g., solar corona, galactic halo, cluster of
  galaxies etc.
• But in plasmas with magnetic field, the effective
  mean free path is given by the ion Larmor radius.
• Hence if the size of phenomenon is much larger
  than the ion Larmor radius, hydrodynamic
  approximation can be used.

10
Applicability of MHD Approximation

• MHD describe macroscopic behavior of plasmas if
• Spatial scale gtgt ion Larmor radius
• Time scale gtgt ion Larmor period
• But MHD can not treat
• Particle acceleration
• Origin of resistivity
• Electromagnetic waves

11
Ideal / Resistive RMHD Eqs
Ideal RMHD
Resistive RMHD
Solve 11 equations (8 in ideal MHD) Need a
closure relation between J and E gt Ohms law
12
Ohms law

• Relativistic Ohms law (Blackman Field 1993
  etc.)

isotropic diffusion in comoving frame (most
simple one)
Lorentz transformation in lab frame
Relativistic Ohms law with istoropic diffusion

• ideal MHD limit (s gt infinity)

Charge current disappear in the Ohms
law (degeneracy of equations, EM wave is decupled)
13
Difference in Ideal Resistive RMHD

• Evolution of EM fields

ideal
resistive

• unnecessary to solve Amperes law in ideal MHD
• E field can be determined directly from Ohms law

• with Ohms law
• 2 additional equations should be solved

• Disadvantage of RRMHD
• Courant condition
• ?E4pq should be satisfied as well as ?B0

14
Basic Equations for Ideal / Resistive RMHD
Resistive RMHD
Ideal RMHD
Hyperbolic equations
Source term
Stiff term
15
Numerical Integration
Resistive RMHD
Constraint
Hyperbolic equations
Solve Relativistic Resistive MHD equations by
taking care of 1. stiff equations appeared in
Amperes law 2. constraints ( no monopole,
Gausss law) 3. Courant conditions (the largest
characteristic wave speed is always light speed.)
Source term
Stiff term
16
For Numerical Simulations
Physical quantities
Basic Equations for RRMHD
Primitive Variables
Conserved Variables
Flux
Source term
Operator splitting (Strangs method) to divide
for stiff term
17
Basics of Numerical RMHD Code
Non-conservative form (De Villier Hawley
(2003), Anninos et al.(2005))
UU(P) - conserved variables, P primitive
variables F- numerical flux of U
where

• Merit
• they solve the internal energy equation rather
  than energy equation. ? advantage in regions
  where the internal energy small compared to total
  energy (such as supersonic flow)
• Recover of primitive variables are fairly
  straightforward
• Demerit
• It can not applied high resolution
  shock-capturing method and artificial viscosity
  must be used for handling discontinuities

18
Basics of Numerical RMHD Code
Conservative form
System of Conservation Equations
UU(P) - conserved variables, P primitive
variables F- numerical flux of U, S - source
of U

• Merit
• Numerically well maintain conserved variables
• High resolution shock-capturing method (Godonuv
  scheme) can be applied to RMHD equations
• Demerit
• These schemes must recover primitive variables P
  by numerically solving the system of equations
  after each step (because the schemes evolve
  conservative variables U)

19
Finite Difference (Volume) Method
Conservative form of wave equation
flux
Finite difference
FTCS scheme
Upwind scheme
Lax-Wendroff scheme
20
Difficulty of Handling Shock Wave
Numerical oscillation (overshoot)
Diffuse shock surface

• Time evolution of wave equation with
  discontinuity using Lax-Wendroff scheme (2nd
  order)

initial

• In numerical hydrodynamic simulations, we need
• sharp shock structure (less diffusivity around
  discontinuity)
• no numerical oscillation around discontinuity
• higher-order resolution at smooth region
• handling extreme case (strong shock, strong
  magnetic field, high Lorentz factor)
• Divergence-free magnetic field (MHD)

21
Flow Chart for Calculation

• Reconstruction
• (Pn cell-center to cell-surface)
• 2. Calculation of Flux at cell-surface
• 3. Integrate hyperbolic equations gt Un1
• 4. Integrate stiff term (E field)
• 5. Convert from Un1 to Pn1

Primitive Variables
Conserved Variables
Flux
Fi-1/2
Ui
Fi1/2
Pi1
Pi-1
Pi
22
Reconstruction

• Minmod MC Slope-limited Piecewise linear
  Method
• 2nd order at smooth region
• Convex ENO (Liu Osher 1998)
• 3rd order at smooth region
• Piecewise Parabolic Method (Marti Muller 1996)
• 4th order at smooth region
• Weighted ENO, WENO-Z, WENO-M (Jiang Shu 1996
  Borges et al. 2008)
• 5th order at smooth region
• Monotonicity Preserving (Suresh Huynh 1997)
• 5th order at smooth region
• MPWENO5 (Balsara Shu 2000)
• Logarithmic 3rd order limiter (Cada Torrilhon
  2009)

Cell-centered variables (Pi) ? right and left
side of Cell-interface variables(PLi1/2, PRi1/2)
Piecewise linear interpolation
PLi1/2
PRi1/2
Pni-1
Pni
Pni1
23
Approximate Riemann Solver
lR, lL fastest characteristic speed
Primitive Variables
Conserved Variables
Flux
HLL flux
Hyperbolic equations
t
lL
lR
M
Ui
Fi-1/2
Fi1/2
R
L
x
Pi1
Pi-1
Pi
If lL gt0 FHLLFL lL lt 0 lt lR ,
FHLLFM lR lt 0 FHLLFR
24
Approximate Riemann Solver

• HLL Approximate Riemann solver single state in
  Riemann fan
• HLLC Approximate Riemann solver two-state in
  Riemann fan (Mignone Bodo 2006, Honkkila
  Janhunen 2007)
• HLLD Approximate Riemann solver six-state in
  Riemann fan (Mignone et al. 2009)
• Roe-type full wave decomposition Riemann solver
  (Anton et al. 2010)

25
Wave speed (for ideal RMHD)

• To calculate numerical flux at each cell-boundary
  via Riemann solver, we need to know wave speed in
  each directions

lEvi, entropy wave
Alfven waves
Magneto-acoustic waves are found from the quartic
equation

• Some simple estimation for fast magnetosonic wave
• gt Leismann et al. (2005), no numerical iteration

26
Constrained Transport (for Ideal MHD)

• The evolution equation can keep divergence free
  magnetic field

Differential Equations

• If treat the induction equation as all other
  conservation laws, it can not maintain divergence
  free magnetic field
• ? We need spatial treatment for magnetic field
  evolution
• Constrained transport scheme
• Evans Hawleys Constrained Transport (need
  staggered mesh)
• Toths constrained transport (flux-CT) (Toth
  2000)
• Fixed Flux-CT, Upwind Flux-CT (Gardiner
  Stone 2005, 2007)
• Other method
• Diffusive cleaning (GLM formulation) (better
  method for AMR or RRMHD)

27
Flux interpolated Constrained Transport
2D case
Toth (2000)
Use the modified flux f that is such a linear
combination of normal fluxes at neighbouring
interfaces that the corner-centred numerical
representation of divB is kept invariant during
integration.
k1/2
k-1/2
j-1/2
j1/2
28
Difficulty of RRMHD
1. Constraint
should be satisfied both constraint numerically
2. Amperes law
Equation becomes stiff at high conductivity
29
Constraints
Approaching Divergence cleaning method (Dedner et
al. 2002, Komissarov 2007)
Introduce additional field F Y (for numerical
noise) advect decay in time
30
Stiff Equation
Komissarov (2007)
Problem comes from difference between dynamical
time scale and diffusive time scale gt analytical
solution
Amperes law
diffusion (stiff) term
Operator splitting method
Hyperbolic source term Solve by HLL method
Analytical solution
source term (stiff part) Solve (ordinary
differential) eqaution
31
Time Evolution (ideal RMHD)
System of Conservation Equations
We use multistep TVD Runge-Kutta method for time
advance of conservation equations (RK2
2nd-order, RK3 3rd-order in time)
RK2, RK3 first step
RK2 second step (a2, b1)
RK3 second and third step (a4, b3)
32
Flow Chart for Calculation (RRMHD)
Strang Splitting Method
Step1 integrate diffusion term in half-time step
Step2 integrate advection term in half-time step
Un(En1/2, Bn)
Step3 integrate advection term in full-time step
Step4 integrate diffusion term in full-time step
(En1, Bn1)Un1
33
Recovery step

• The GRMHD code require a calculation of primitive
  variables from conservative variables.
• The forward transformation (primitive ?
  conserved) has a close-form solution, but the
  inverse transformation (conserved ? primitive)
  requires the solution of a set of five nonlinear
  equations
• Method
• Nobles 2D method (Noble et al. 2005)
• Mignone McKinneys method (Mignone McKinney
  2007)

34
Nobles 2D method

• Conserved quantities(D,S,t,B) ? primitive
  variables (r,p,v,B)
• Solve two-algebraic equations for two
  independent variables Whg2 and v2 by using
  2-variable Newton-Raphson iteration method

W and v2 ?primitive variables r p, and v

• Mignone McKinney (2007) Implemented from
  Nobles method for variable EoS

35
General (Approximate) EoS
Mignone McKinney 2007

• In the theory of relativistic perfect single
  gases, specific enthalpy is a function of
  temperature alone (Synge 1957)

Q temperature p/r K2, K3 the order 2 and 3 of
modified Bessel functions

• Constant G-law EoS (ideal EoS)
• G constant specific heat ratio

• Taubs fundamental inequality(Taub 1948)

Q ? 0, Geq ? 5/3, Q ? 8, Geq ? 4/3
Solid Synge EoS Dotted ideal G5/3 Dashed
ideal G4/3 Dash-dotted TM EoS

• TM EoS (approximate Synges EoS)
• (Mignone et al. 2005)

c/sqrt(3)
36
2. 1. Numerical Testsideal RMHD
37
1D Relativistic MHD Shock-Tube
Exact solution Giacomazzo Rezzolla (2006)
38
1D Relativistic MHD Shock-Tube
Mizuno et al. 2006
Balsara Test1 (Balsara 2001)

• The results show good agreement of the exact
  solution calculated by Giacommazo Rezzolla
  (2006).
• Minmod slope-limiter and CENO reconstructions
  are more diffusive than the MC slope-limiter and
  PPM reconstructions.
• Although MC slope limiter and PPM
  reconstructions can resolve the discontinuities
  sharply, some small oscillations are seen at the
  discontinuities.

FR
SR
SS
FR
CD
Black exact solution, Blue MC-limiter, Light
blue minmod-limiter, Orange CENO, red PPM
400 computational zones
39
Advection of Magnetic Field Loop
2D
No advection

• Advection of a weak magnetic field loop in an
  uniform velocity field
• 2D (vx, vy)(0.6,0.3)
• 3D (vx,vy,vz)(0.3,0.3,0.6)
• Periodic boundary in all direction
• Run until return to initial position in advection
  case

B2
Advection
Volume-averaged magnetic energy density (2D)
3D
B2
Nx512
256
128
Advection
No advection
40
Cylindrical Explosion

• Cylindrical explosion in magnetized medium
• (propagation of strong shock waves)
• Pc1.0, rc0.01 (Rlt1.0)
• Pe5x10-4, re10-4
• Bx0.1 (uniform)

41
Numerical Tests
42
1D CP Alfven wave propagation test

• Aim Recover of ideal RMHD regime in high
  conductivity
• Propagation of large amplitude circular-polarized
  Alfven wave along uniform magnetic field
• Exact solution Del Zanna et al.(2007) in ideal
  RMHD limit

BxB0, vx0, k wave number, zA amplitude of wave

• p1, B00.46188 gt vA0.25c, ideal EoS with G2

Using high conductivity s105
43
1D CP Alfven wave propagation test
Numerical results at t4 (one Alfven crossing
time)
Solid exact solution Dotted Nx50 Dashed
Nx100 Dash-dotted Nx200
New RRMHD code reproduces ideal RMHD solution
when conductivity is high
L1 norm errors of magnetic field By almost 2nd
order accuracy
44
1D Self-Similar Current Sheet Test

• Assumption Magnetic pressure ltlt gas pressure
• Magnetic field configuration B0, By(x,t),0,
  By(x,t) changes sign within a thin current layer
  (thickness Dl)
• The evolution of thin current layer is a slow
  diffusive expansion due to resistivity and
  described by diffusion equations
• As the thickness of the layer becomes much larger
  than Dl, the expansion becomes self-similar
• ct/x2, erf error function
• Test simulation
• Initial solution at t1 with p50, r1, Ev0,
  s100 with G2

45
1D Self-Similar Current Sheet Test
dotted dashed analytical solution at t1
10 Solid numerical solution at t10
Numerical Simulation shows good agreement with
exact solutions with moderate conductivity regime
46
1D Shock-Tube Test (Brio Wu)

• Aim Check the effect of resistivity
  (conductivity)
• Simple MHD version of Brio Wu test
• (rL, pL, ByL) (1, 1, 0.5), (rR, pR,
  ByR)(0.125, 0.1, -0.5)
• Ideal EoS with G2

Orange solid s0 Green dash-two-dotted s10 Red
dash-dotted s102 Purple dashed s103 Blue
dotted s105 Black solid exact solution in
ideal RMHD
Smooth change from a wave-like solution (s0) to
ideal-MHD solution (s105)
47
1D Shock-Tube Test (Balsara 2)

• Aim check the effect of choosing EoS in RRMHD
• Balsara Test 2
• Using ideal EoS (G5/3) approximate TM EoS
• Changing conductivity from s0 to 103
• Mildly relativistic blast wave propagates with
  1.3 lt g lt 1.4

48
1D Shock-Tube Test (Balsara 2)
SS FS
CD
FR
SR
Purple dash-two-dotted s0 Green dash-dotted
s10 Red dashed s102 Blue dotted s103 Black
solid exact solution in ideal RMHD
The solutions Fast Rarefaction, Slow
Rarefaction, Contact Discontinuity, Slow Shock
and Fast Shock.
49
1D Shock-Tube Test (Balsara 2)
SS FS
FR
CD
SR

• The solutions are same but quantitatively
  different.
• rarefaction waves and shocks propagate with
  smaller velocities
• lt lower sound speed in TM EoSs relatives to
  overestimated sound speed in ideal EoS
• these properties are consistent with in ideal
  RMHD case

Purple dash-two-dotted s0 Green dash-dotted
s10 Red dashed s102 Blue dotted s103 Black
solid exact solution in ideal RMHD
50
2D Cylindrical Explosion

• Cylindrical blast wave expanding into uniform
  magnetic field
• Standard test for ideal RMHD code
• No exact solution in multi-dimensional tests
• Initial Condition
• Density, pressure
• Rlt0.8, p1, r0.01
• 0.8ltRlt1.0 decrease exponentially to ambient gas
• Rgt1.0 pr0.001
• Magnetic field uniform in x-direction with
  B0.05
• Using different conductivity s0-105
• Using ideal EoS with G4/3 and TM EoS

51
Global Structure in Cylindrical Explosion
s105

• qualitatively similar to ideal RMHD results
• No different between ideal EoS and TM EoS

52
1D cut of Cylindrical Explosion
Purple dash-two-dotted s0 Green dash-dotted
s10 Red dashed s102 Blue dotted s103 Black
solid s105

• At high conductivity (s gt 103) no difference
  (recovers ideal MHD solution)
• conductivity ? maximum gas and mag pressure ?
• No mag pressure increase for s0

53
2D Kelvin-Helmholtz Instability

• Linear and nonlinear growth of 2D
  Kelvin-Helmholtz instability (KHI) magnetic
  field amplification via KHI
• Initial condition
• Shear velocity profile
• Uniform gas pressure p1.0
• Density r1.0 in the region vsh0.5, r10-2 in
  the region vsh-0.5
• Magnetic Field
• Single mode perturbation
• Simulation box
• -0.5 lt x lt 0.5, -1 lt y lt 1

a0.01, characteristic thickness of shear
layer vsh0.5 gt relative g2.29
mp0.5, mt1.0
A00.1, a0.1
54
Growth Rate of KHI
Amplitude of perturbation
Volume-averaged Poloidal field

• Initial linear growth with almost same growth
  rate
• Maximum amplitude transition from linear to
  nonlinear
• Poloidal field amplification via stretching due
  to main vortex developed by KHI
• Larger poloidal field amplification occurs for
  TM EoS than for ideal EoS

Purple dash-two-dotted s0 Green dash-dotted
s10 Red dashed s102 Blue dotted s103 Black
solid s105
55
2D KHI Global Structure (ideal EoS)

• Formation of main vortex by growth of KHI in
  linear growth phase
• secondary vortex?
• main vortex is distorted and stretched in
  nonlinear phase
• B-field amplified by shear in vortex in linear
  and stretching in nonlinear

56
2D KHI Global Structure (TM EoS)

• Formation of main vortex by growth of KHI in
  linear growth phase
• no secondary vortex
• main vortex is distorted and stretched in
  nonlinear phase
• vortex becomes strongly elongated in nonlinear
  phase
• Created structure is very different in ideal and
  TM EoSs

57
Field Amplification in KHI

• Field amplification structure for different
  conductivities
• Conductivity low, magnetic field amplification
  is weaker
• Field amplification is a result of fluid motion
  in the vortex
• B-field follows fluid motion, like ideal MHD,
  strongly twisted in high conductivity
• conductivity decline, B-field is no longer
  strongly coupled to the fluid motion
• Therefore B-field is not strongly twisted

58
2D Relativistic Magnetic Reconnection

• Consider Pestchek-type reconnection
• Initial condition Harris-like model

Uniform density gas pressure outside current
sheet, rbpb0.1
Density gas pressure
Magnetic field
Current
Resistivity (anomalous resistivity in rltrh)
hb1/sb10-3, h01.0, rh0.8
Electric field
59
Global Structure of Relativistic MR
Plasmoid
t100
Slow shock
Strong current flow
60
Time Evolution of Relativistic MR

• Outflow gradually accelerates and saturates t60
  with vx0.8c
• TM EoS case slightly faster than ideal EoS case
• Magnetic energy converted to thermal and kinetic
  energies (acceleration of outflow)
• TM EoS case has larger reconnection rate than
  ideal EoS.
• Different EoSs lead to a quantitative difference
  in relativistic magnetic reconnection

Reconnection outflow speed
Solid ideal EoS Dashed TM EoS
Magnetic energy
Reconnection rate
time
61
Summary

• In 1D tests, new RRMHD code is stable and
  reproduces ideal RMHD solutions when the
  conductivity is high.
• 1D shock tube tests show results obtained from
  approximate EoS are considerably different from
  ideal EoS.
• In KHI tests, growth rate of KHI is independent
  of the conductivity
• But magnetic field amplification via stretching
  of the main vortex and nonlinear behavior
  strongly depends on the conductivity and choice
  of EoSs
• In reconnection test, approximate EoS case
  resulted in a faster reconnection outflow and
  larger reconnection rate than ideal EoS case