3D GRMHD code RAISHIN and Relativistic Jet Simulations

1. 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

2. 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

3. Relativistic Jets in Universe Mirabel & Rodoriguez 1998

4. 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

5. 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

6. 1. Development of 3D GRMHD Code “RAISHIN” Mizuno et al. 2006a, preprint, Astro-ph/0609004 Mizuno et al. 2006, proceedings of science, MQW6, 045

7. 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 )

8. 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).

9. 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

10. Finite Difference Method Conservative form of wave equation flux Finite difference FTCS scheme Upwind scheme Lax-Wendroff scheme

11. 4D General Relativistic MHD Equation • General relativistic equation of conservation laws and Maxwell equations: 　　　　　　　∇n (r Un) = 0(conservation law of particle-number) 　　　　　　　∇nTmn= 0(conservation law of energy-momentum) 　　　　　　　∂mFnl + ∂nFlm + ∂lF mn= 0 　　　　　　　∇mFmn= - Jn • Ideal MHD condition:FnmUn= 0 • metric：ds2=-a2 dt2+gij (dxi+b i dt)(dx j+b 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 Un+FmsFns -gmnFlkFlk/4. Fmn : field-strength tensor,

12. 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.

13. 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

14. Reconstruction Cell-centered variables (Pi) → right and left side of Cell-interface variables(PLi+1/2, PRi+1/2) • Minmod and MCSlope-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) Piecewise linear interpolation PLi+1/2 PRi+1/2 Pni-1 Pni Pni+1

15. 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

16. 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

17. 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)

18. Flux interpolated Constrained Transport ２D 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. k+1/2 k-1/2 j-1/2 j+1/2

19. 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)

20. 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 • Koide’s 2D method (Koide, Shibata & Kudoh 1999) • Noble’s 2D method (Noble et al. 2005) • Mignone & McKinney’s method (Mignone & McKinney 2007)

21. Noble’s 2D method • Conserved quantities(D,S,t,B) → primitive variables (r,p,v,B) • Solve two-algebraic equations for two independent variablesW≡hg2and v2 by using 2-variable Newton-Raphson iteration method W and v2 →primitive variables r p, and v • Mignone & McKinney (2007): Implemented from Noble’s method for variable EoS

22. 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)

23. 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

24. Relativistic MHD Shock-Tube Tests Exact solution: Giacomazzo & Rezzolla (2006)

25. Relativistic MHD Shock-Tube Tests 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

26. 2. Highlights of Jet Simulations

27. Hardee, Mizuno & Nishikawa (2007) Spine-Sheath Relativistic Jets (GRMHD Simulations) Non-rotating BH Fast-rotating BH 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) • jet spine: Formed by twisted magnetic field by frame-dragging effect of rotating black hole • broad sheath wind: Formed by twisted magnetic field by rotation of accretion disk Color: total velocity Disk Jet/Wind BH Jet Disk Jet/Wind 2D GRMHD Simulation of jet formation

28. Radiation Images of Black Hole-Disk System Wu, Fuerst, Mizuno et al. (2008) • Calculation of thermal free-free emissionandthermal 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 two-component jet based on synchrotron emission and the strong thermal radiation from hot dense gas near the BHs. Schematic picture of Ray-tracing method Radiation image seen from q=85 (optically thick) Radiation image seen from q=85 (optically thin)

29. 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

30. 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)

31. Mizuno et al. (2009b) Current-Driven Instability: Static Plasma Decrease density with Constant pitch case: CD kink instability leads to a helically twisted density and magnetic filament • 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 non-linear evolution depends on the radial density profile and strongly depends on the magnetic pitch profile. Decrease pitch Increase pitch

32. CD kink instability of Sub-Alfvenic Jets:Spatial Properties Mizuno et al. 2010, in preparation 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 3D density with magnetic field lines • 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.

33. Magnetic Field Amplification by Relativistic Shocks in Turbulent Medium Mizuno et al. 2010 in prep Time evolution 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.

34. 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.)