340 likes | 353 Vues
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.
E N D
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 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,
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) • 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
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 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. k+1/2 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 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)
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
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) • 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
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
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)
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)
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
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.
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.
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.)