Transcript PPT

3D GRMHD code RAISHIN and
Relativistic Jet Simulations
Yosuke Mizuno
Center for Space Plasma and Aeronomic Research
University of Alabama in Huntsville
CFD-MHD seminar, ASIAA, Taiwan, 3/31/10
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
Relativistic Jets in Universe
Mirabel & Rodoriguez 1998
Requirement of Relativistic MHD
• Astrophysical jets seen AGNs show the relativistic
speed (~0.99c)
• The central object of AGNs is supper-massive black
hole (~105-1010 solar mass)
• The jet is formed near black hole
Require relativistic treatment (special or general)
• In order to understand the time evolution of jet
formation, propagation and other time dependent
phenomena, we need to perform relativistic
magnetohydrodynamic simulations
Applicability of MHD Approximation
• MHD describe macroscopic behavior of plasmas
if
– Spatial scale >> ion Larmor radius
– Time scale >> ion Larmor period
• But MHD can not treat
– Particle acceleration
– Origin of resistivity
– Electromagnetic waves
1. Development of 3D GRMHD
Code “RAISHIN”
Mizuno et al. 2006a, preprint, Astro-ph/0609004
Mizuno et al. 2006, proceedings of science, MQW6, 045
Numerical Approach to Relativistic MHD
• RHD: reviews Marti & Muller (2003) and Fonts (2003)
• SRMHD: many groups developed their own code
– Application: relativistic Riemann problems, relativistic jet propagation, jet
stability, pulsar wind nebule, relativistic shock/blast wave etc.
• GRMHD
– Fixed spacetime (Koide, Shibata & Kudoh 1998; De Villiers &
Hawley 2003; Gammie, McKinney & Toth 2003; Komissarov 2004; Anton
et al. 2005, 2010; Annios, Fragile & Salmonson 2005; Del Zanna et al.
2007, Nagataki 2009…)
– Application: The structure of accretion flows onto black hole and/or
formation of jets, BZ process near rotating black hole, the formation of
GRB jets in collapsars etc.
– Dynamical spacetime (Duez et al. 2005; Shibata & Sekiguchi 2005;
Anderson et al. 2006; Giacomazzo & Rezzolla 2007 )
Propose to Make a New GRMHD Code
(statement at 2006)
• The Koide’s GRMHD Code (Koide, Shibata & Kudoh 1999;
Koide 2003) has been applied to many high-energy
astrophysical phenomena and showed pioneering results.
• However, the code can not perform calculation in highly
relativistic (g>5) or highly magnetized regimes.
• The critical problem of the Koide’s GRMHD code is the
schemes can not guarantee to maintain divergence free
magnetic field.
• In order to improve these numerical difficulties, we have
developed a new 3D GRMHD code RAISHIN (RelAtIviStic
magnetoHydrodynamc sImulatioN, RAISHIN is the Japanese
ancient god of lightning).
Finite Difference Method
Linear wave equation
Hydrodynamic equations
are a set of wave equations
Finite difference in each term
df/dx as a function of fj+1,n fj,n fj-1,n
Forward derivative
Backward derivative
Central derivative
Finite Difference Method
Conservative form of wave equation
flux
Finite difference
FTCS scheme
Upwind scheme
Lax-Wendroff scheme
4D General Relativistic MHD Equation
•
General relativistic equation of conservation laws and Maxwell equations:
∇n ( r U n ) = 0
∇n T mn = 0
(conservation law of particle-number)
(conservation law of energy-momentum)
∂mFnl  nFlm  lF mn = 0
∇mF
mn
=-J
n
(Maxwell equations)
FnmUn = 0
•
Ideal MHD condition:
•
metric： ds2=-a2 dt2+gij (dxi+b i dt)(dx j+b j dt)
•
Equation of state : p=(G-1) u
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.
mu
U : velocity four vector. Jmu : current density four vector.
mn
∇ : covariant derivative. gmn : 4-metric. a : lapse function, bi : shift vector, gij : 3-metric
mn
mn
mn
m n
ms n
lk
T : energy momentum tensor, T = pg + r h U U +F F s -gmnF Flk/4.
Fmn : field-strength tensor,
Conservative Form of GRMHD
Equations (3+1 Form)
(Particle number conservation)
(Momentum conservation)
(Energy
conservation)
(Induction equation)
U (conserved variables)
Fi (numerical flux) S (source term)
√-g : determinant of 4-metric
√g : determinant of 3-metric
Detail of derivation of
GRMHD equations
Anton et al. (2005) etc.
Detailed Features of the Numerical
Schemes
Mizuno et al. 2006a, astro-ph/0609004
and progress
• RAISHIN utilizes conservative, high-resolution shock capturing
schemes (Godunov-type scheme) to solve the 3D GRMHD
equations (metric is static)
* Reconstruction: PLM (Minmod & MC slope-limiter), convex ENO,
PPM, Weighted ENO5, Monotonicity Preserving5, MPWENO5
* Riemann solver: HLL, HLLC approximate Riemann solver
* Constrained Transport: Flux CT, Fixed Flux-CT, Upwind Flux-CT
* Time evolution: Multi-step Runge-Kutta method (2nd & 3rd-order)
* Recovery step: Koide 2 variable method, Noble 2 variable method,
Mignore-McKinney 1 variable method
* Equation of states: constant G-law EoS, variable EoS for ideal gas
Reconstruction
Cell-centered variables (Pi)
→ right and left side of Cell-interface
variables(PLi+1/2, PRi+1/2)
Piecewise linear interpolation
PLi+1/2
PRi+1/2
Pni-1
Pn i
Pni+1
• Minmod and MC Slope-limited
Piecewise linear Method
• 2nd order at smooth region
• Convex CENO (Liu & Osher 1998)
• 3rd order at smooth region
• Piecewise Parabolic Method (Marti
& Muller 1996)
• 4th order at smooth region
• Weighted ENO5 (Jiang & Shu
1996)
• 5th order at smooth region
• Monotonicity Preserving (Suresh &
Huynh 1997)
• 5th order at smooth region
• MPWENO5 (Balsara & Shu 2000)
HLL Approximate Riemann Solver
• Calculate numerical flux at cell-inteface from reconstructed
cell-interface variables based on Riemann problem
• We use HLL approximate Riemann solver
• Need only the maximum left- and right- going wave speeds
(in MHD case, fast magnetosonic mode)
HLL flux
FR=F(PR), FL=F(PL); UR=U(PR), UL=U(PL)
SR=max(0,c+R, c+L); SL=max(0,c-R,c-L)
If SL >0
FHLL=FL
SL < 0 < SR , FHLL=FM
SR < 0
FHLL=FR
HLLC Approximate Riemann Solver
Mignore & Bodo (2006)
Honkkila & Janhunen (2007)
• HLL Approximate Riemann solver: single state in Riemann fan
• HLLC Approximate Riemann solver: two-state in Riemann fan
• (HLLD Approximate Riemann solver: six-state in Riemann fan)
(Mignone et al. 2009)
HLL
HLLC
Constrained Transport
Differential Equations
- The
evolution
equation can keep
divergence free
magnetic field
• 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 method
• Evans & Hawley’s Constrained Transport (Komissarov (1999,2002,2004),
de Villiers & Hawley (2003), Del Zanna et al.(2003), Anton et al.(2005))
• Toth’s constrained transport (flux-CT) (Gammie et al.(2003), Duez et
al.(2005))
• Fixed Flux-CT, Upwind Flux-CT (Gardiner & Stone 2005, 2007)
• Diffusive cleaning (Annios et al.(2005) etc) (better method for AMR or RRMHD)
Flux interpolated Constrained Transport
２D case
Toth (2000)
Use the “modified flux” f that is such a linear
combination of normal fluxes at neighbouring
k+1/2 interfaces that the “corner-centred” numerical
representation of divB is kept invariant during
integration.
k-1/2
j-1/2
j+1/2
Time evolution
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 (a=2, b=1)
RK3: second and third step (a=4, b=3)
Recovery step
• The GRMHD code require a calculation of primitive variables
from conservative variables.
• The forward transformation (primitive → conserved) has a closeform solution, but the inverse transformation (conserved →
primitive) requires the solution of a set of five nonlinear equations
Method
• Koide’s 2D method (Koide, Shibata & Kudoh 1999)
• Noble’s 2D method (Noble et al. 2005)
• Mignone & McKinney’s method (Mignone & McKinney 2007)
Noble’s 2D method
quantities(D,S,t,B) → primitive variables (r,p,v,B)
• Solve two-algebraic equations for two independent variables
W≡hg2 and v2 by using 2-variable Newton-Raphson iteration
method
• Conserved
W and v2 →primitive variables r p, and v
• Mignone
& McKinney (2007): Implemented from Noble’s method
for variable EoS
Variable EoS
Mignone & McKinney 2007
• In the theory of relativistic perfect 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:
• G: constant specific heat ratio
• Taub’s fundamental inequality
(Taub 1948)
Q → 0, Geq → 5/3, Q → ∞, Geq → 4/3
• TM EoS
(Mathew 1971, Mignone et al. 2005)
Ability of RAISHIN code
(current status)
• 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 (7)
• Different approximate Riemann solver (2)
• Different constrained transport schemes (3)
• Different time advance algorithms (2)
• Different recovery schemes (3)
• Using constant G-law and variable Equation of State (Synge-type)
• Parallelized by OpenMP
Relativistic MHD Shock-Tube Tests
Exact solution: Giacomazzo & Rezzolla (2006)
Relativistic MHD Shock-Tube Tests
Mizuno et al. 2006
Balsara Test1 (Balsara 2001)
FR
SR
SS
FR
CD
Black: exact solution, Blue: MC-limiter,
Light blue: minmod-limiter, Orange: CENO,
red: PPM
400 computational zones
• 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 slopelimiter and PPM reconstructions.
• Although MC slope limiter and
PPM reconstructions can resolve
the discontinuities sharply, some
small oscillations are seen at the
discontinuities.
2. Highlights of Jet Simulations
Hardee, Mizuno & Nishikawa (2007)
Spine-Sheath Relativistic Jets
Non-rotating BH
Fast-rotating BH
(GRMHD Simulations)
Color: density
• Two component (Spine-Sheath)
jet structure is seen in recent
GRMHD simulations of jet
formation in black hole-accretion
disk system (e.g., Hawley &
Krolik 2006, McKinney 2006,
Hardee et al. 2007)
Color: total velocity
• jet spine: Formed by twisted
magnetic field by framedragging effect of rotating
black hole
• broad sheath wind:
Formed by twisted magnetic
field by rotation of accretion
disk
Disk Jet/Wind
BH Jet
Disk Jet/Wind
2D GRMHD Simulation of jet formation
Radiation Images of Black HoleDisk System
• Calculation of thermal free-free
emission and thermal synchrotron
emission by ray-tracing method
considered GR radiation transfer
from a relativistic flows in black
hole systems (2D GRMHD
simulation, rotating BH cases).
• The radiation image shows the
front side of the accretion disk and
the other side of the disk at the top
and bottom regions because the GR
effects.
• We can see the formation of twocomponent jet based on
synchrotron emission and the
strong thermal radiation from hot
dense gas near the BHs.
Wu, Fuerst, Mizuno et al. (2008)
Schematic picture of Ray-tracing method
Radiation image seen from Radiation image seen from
q=85 (optically thick)
q=85 (optically thin)
Stability of Magnetized Spine-Sheath
Relativistic Jets
Mizuno, Hardee, & Nishikawa (2007)
• We investigate the stability of magnetized two-component (spine-sheath)
relativistic jets against Kelvin-Helmholtz (KH) instability by using 3D relativistic
MHD simulations.
T=0
• Cylindrical super-Alfvenic jet
established across the
computational domain with a
parallel magnetic field
• Put precession perturbation
from jet inlet to break symmetry
T=60 (Weakly magnetized, static external medium case)
vj
• The jet is disrupted by the
growing KH instability
Effect of magnetic field and sheath wind
Mizuno, Hardee & Nishikawa (2007)
1D radial velocity profile along jet
ue=0.0
ue=0.5c
ue=0.0
ue=0.5c
•The sheath flow reduces the growth rate of KH modes and slightly increases the
wave speed and wavelength as predicted from linear stability analysis.
•The magnetized sheath reduces growth rate relative to the weakly magnetized case
•The magnetized sheath flow damped growth of KH modes.
Criterion for damped KH modes:
(linear stability analysis)
Current-Driven Instability:
Static Plasma
•We studied the development of
current-driven (CD) kink
instability of a static force-free
helical magnetic field
configuration by using 3D
RMHD simulations.
•We found the initial
configuration is strongly
distorted but not disrupted.
•The linear growth and nonlinear evolution depends on the
radial density profile and
strongly depends on the
magnetic pitch profile.
Mizuno et al. (2009b)
Decrease density
with Constant pitch
case: CD kink
instability leads to a
helically twisted
density and
magnetic filament
Increase pitch
Decrease pitch
CD kink instability of Sub-Alfvenic Jets:
Spatial Properties
Initial Condition
• Cylindrical sub-Alfvenic jet established
across the computational domain with a
helical force-free magnetic field (stable
against KH instabilities)
– Vj=0.2c, Rj=1.0
• Radial profile: Decreasing density with
constant magnetic pitch
• Jet spine precessed to break the symmetry
Preliminary Result
• Precession perturbation from jet inlet
produces the growth of CD kink
instability with helical density distortion.
• Helical structure propagates along the jet
with continuous growth of kink amplitude
in non-linear phase.
Mizuno et al. 2010, in preparation
3D density with magnetic field lines
Magnetic Field Amplification by Relativistic
Shocks in Turbulent Medium
Initial condition
• Density: mean + small inhomogenity
with 2D Kolmogorov-like power-law
spectrum
• Relativistic flow in whole region with
constant magnetic field (parallel to
shock propagation direction)
• a rigid reflecting boundary at x=xmax
to create the shock. (shock propagates
in –x direction)
Preliminary result
• Density inhomogenity induces turbulent
motion in shocked region
• Turbulence motion stretch and deform the
magnetic field lines and create filamentary
structure with strong field amplification.
Mizuno et al. 2010 in prep
Time evolution
Future Implementation of RAISHIN
• Resistivity (extension to non-ideal MHD; e.g., Watanabe & Yokoyama
2007; Komissarov 2007; etc)
• 2 fluid MHD with resistivity (Zenitani et al. 2008)
• Couple with radiation transfer (link to observation: collaborative
works with Dr. Wu)
• Kerr-Schild coordinates (to avoid singularity at BH radius in
GRMHD simulations)
• Improve the realistic EOS
• Include Effect of radiation and neutrino emisson (cooling, heating)
• Include Nucleosysthesis (post processing)
• Couple with Einstein equation (dynamical spacetime)
• Adaptive mesh refinement (AMR)
• Parallerization by MPI for PC cluster type supercomputer
• Apply to astrophysical phenomena in which relativistic
outflows/shocks and/or GR essential (AGNs, microquasars,
neutron stars, and GRBs etc.)