Development of General Relativistic Magnetohydrodynamic Code and its Application to Central Engine of Long Gamma-Ray Bursts
aa r X i v : . [ a s t r o - ph . H E ] F e b Draft version October 29, 2018
Preprint typeset using L A TEX style emulateapj v. 04/20/08
DEVELOPMENT OF GENERAL RELATIVISTIC MAGNETOHYDRODYNAMIC CODE AND ITSAPPLICATION TO CENTRAL ENGINE OF LONG GAMMA-RAY BURSTS
Shigehiro Nagataki Draft version October 29, 2018
ABSTRACTIn order to investigate formation of relativistic jets at the center of a progenitor of a long gamma-ray burst (GRB), we develop a two-dimensional general relativistic magnetohydrodynamic (GRMHD)code. We show the code passes many, well-known test calculations, by which the reliability of thecode is confirmed. Then we perform a numerical simulation of a collapsar using a realistic progenitormodel. It is shown that a jet is launched from the center of the progenitor. We also find that themass accretion rate after the launch of the jet shows rapid time variability that resembles to a typicaltime profile of a GRB. The structure of the jet is similar to the previous study: a poynting flux jetis surrounded by a funnel-wall jet. Even at the final stage of the simulation, bulk Lorentz factor ofthe jet is still low, and total energy of the jet is still as small as 10 erg. However, we find that theenergy flux per unit rest-mass flux is as high as 10 at the bottom of the jet. Thus we conclude thatthe bulk Lorentz factor of the jet can be potentially high when it propagates outward. It is shownthat the outgoing poynting flux exists at the horizon around the polar region, which proves that theBlandford-Znajek mechanism is working. However, we conclude that the jet is launched mainly bythe magnetic field amplified by the gravitational collapse and differential rotation around the blackhole, rather than the Blandford-Znajek mechanism. Subject headings: gamma rays: bursts — relativity — black hole physics — accretion, accretion discs— supernovae: general INTRODUCTION
Gamma-Ray Bursts (GRBs; in this study, we consideronly long GRBs, so we refer to long GRBs as GRBshereafter) have been mysterious phenomena since theirdiscovery in 1969 (Klebesadel et al. 1973). Last decade,observational evidence for supernovae (SNe) and GRBsassociation has been reported (e.g. Woosley and Bloom2006, and references therein).Some of the SNe that associate with GRBs were veryenergetic and blight. The estimated explosion energy wasof the order of 10 ergs, and produced nickel mass was ∼ . M ⊙ . Thus they are categorized as a new type ofSNe (sometimes called as hypernovae). The largeness ofthe explosion energy is very important, because it can notbe explained by the standard core-collapse SN scenario,and other mechanism should be working at the center ofthe progenitors.The promising scenarios are the collapsar sce-nario (Woosley 1993) and the magnetar sce-nario (Usov 1992). In the collapsar scenario, arapidly rotating black hole (BH) is formed at thecenter, while a rapidly rotating neutron star with strongmagnetic fields ( ∼ G) is formed in the magnetarscenario. Many numerical simulations have been donefor the collapsar scenario (MacFadyen & Woosley 1999;Proga et al. 2003; Proga & Begelman 2003;Mizuno et al. 2004a; Mizuno et al. 2004b;Proga 2005; Fujimoto et al. 2006; Shibata et al. 2006;Nagataki et al. 2007; Sekiguchi & Shibata 2007;Suwa et al. 2007; Barkov & Komissarov 2008a) Yukawa Institute for Theoretical Physics, Kyoto University,Oiwake-cho Kitashirakawa Sakyo-ku, Kyoto 606-8502, Japan,[email protected] and the magnetar scenario (Takiwaki et al. 2004;Komissarov & Barkov 2007; Burrows et al. 2007;Bucciantini et al. 2008; Dessart et al. 2008;Takiwaki et al. 2008; Bucciantini et al. 2009). Inthis study, we investigate the collapsar scenario.In the collapsar scenario, a BH is formed as a resultof gravitational collapse. Also, rotation of the pro-genitor plays an essential role. Due to the rotation,an accretion disk is formed around the equatorialplane. On the other hand, the matter around therotation axis falls into the BH almost freely. It ispointed out that the jet-induced explosion along therotation axis may occur due to the heating throughpair annihilation of neutrinos and anti-neutrinos thatare emitted from the accretion disk (Woosley 1993;MacFadyen & Woosley 1999; Fryer & M´esz´aros 2000).Effect of extraction of rotation energy from the accretiondisk by magnetic field lines that leave the disk surface(Blandford-Payne effect (Blandford & Payne 1982)) isalso investigated by several authors (Proga et al. 2003;Proga & Begelman 2003; Mizuno et al. 2004a;Mizuno et al. 2004b; Proga 2005; Fujimoto et al. 2006;Nagataki et al. 2007; Suwa et al. 2007). Recently,the effect of extraction of rotation energy from theBH through outgoing poynting flux (Blandford-Znajek effect (Blandford & Znajek 1977)) is inves-tigated (Barkov & Komissarov 2008a). In order toinvestigate the collapsar scenario completely, a high-quality numerical code including effects of a lot ofmicrophysics (neutrino physics, nuclear physics, andequation of state for dense matter) and macrophysics(magneto-hydrodynamics, general relativity) has to bedeveloped. Although many numerical studies have beenreported, such a numerical code has not been developed S. Nagatakiyet. Thus we have to develop our numerical code stepby step.In this study, we investigate the dynamics of collapsarstaking into account the general relativistic effects. Ex-traction of rotation energy from a rotating BH is one ofthem. Also, even when the rotation energy is extractedfrom the accretion disk, the properties of the accretiondisk should depend on the properties of the BH: if theBH is rotating, the inner region of the accretion diskshould be enforced to co-rotates with the BH. We inves-tigate how a jet is launched at the center of a progenitor,and how the property of the jet is. Effects of rotationof the BH on the formation of GRB jet have not beeninvestigated so much. Barkov and Komissarov (2008)is a pioneering study. However, only one case is inves-tigated in their study, and the initial progenitor modelthey used is a simplified one-dimensional model withoutrotation and magnetic fields (Bethe 1990). Since thereshould be many initial conditions of progenitors (pro-genitor mass, metallicity, angular momentum, magneticfields), it should be important to investigate the generalrelativistic effects using a different initial condition fromthe previous study. In this study, we use a realistic initialcondition for the progenitor model that is developed byWoosley and Heger (2006), in which rotation and mag-netic fields are taken into account.When we investigate the general relativistic effects,one has to develop a General Relativistic Magneto-Hydro Dynamic (GRMHD) code. So far, there aremany studies on GRMHD code for fixed backgroundspace times using high-order conservative schemes basedon either approximate or full wave-decomposition Rie-mann solvers (Gammie et al. 2003; Komissarov 2005;Anninos et al. 2005; Ant´on et al. 2006;Del Zanna et al. 2007; Tchekhovskoy et al. 2007) ornon-conservative schemes (De Villiers & Hawley 2003;Anninos et al. 2005). Since the accreted mass onto theBH is still less than the initial BH mass in this study,we take the GRMHD code for the fixed background.Especially, we develop our code using the conservativescheme of Gammie et al. (2003) with the method ofNoble et al. (2006) for transforming conserved variablesto primitive variables.The plan of the paper is as follows. In section 2, wepresent the formulation of the GRMHD code. In sec-tion 3, we show results of many, well-known test calcula-tions to confirm the reliability of the code. After we showthe reliability, we present results of numerical simulationsof collapsars in section 4. Summary and discussion arepresented in section 5. DEVELOPMENT OF GRMHD CODE
We have developed a two-dimensional GRMHD codefollowing Gammie et al. (2003) and Noble et al.(2006). We have adopted a conservative, shock-capturingscheme with Harten, Lax, and van Leer (HLL) fluxterm (Harten et al. 1983) with flux-interpolated con-strained transport technique (T´oth 2000). We use athird-order Total Variation Diminishing (TVD) Runge-Kutta method for evolution in time, while monotonizedcentral slope-limited linear interpolation method is usedfor second-order accuracy in space (van Leer 1977). 2Dscheme (2-dimensional Newton-Raphson method) is usu-ally adopted for transforming conserved variables to primitive variables (Noble et al. 2006).When we perform simulations of GRMHD, ModifiedKerr-Schild coordinate is basically adopted with massof the BH ( M ) fixed where the Kerr-Schild radius r isreplaced by the logarithmic radial coordinate x = ln r .When we show the result, the coordinates are sometimestransfered from Modified Kerr-Schild coordinate to Kerr-Schild one for convenience. In the following, we use G = M = c = 1 unit. G is the gravitational constant, c is thespeed of light, and M is the gravitational mass of theBH at the center. Throughout this paper we follow thestandard notation (Misner et al. 1970). Formalism
Number of variables that appear in the equations ofGRMHD is 13: rest-mass density ( ρ ), internal energydensity ( u ), pressure ( p ), four-velocity of fluid ( u µ ),and Faraday tensor ( F µν ). Note that Faraday tensorhas only 6 independent components due to the relation F µν = − F νµ . We can reduce the number of independentvariables to 8 using the MHD condition ( u µ F µν = 0),equation of state ( p = ( γ − u : γ -law gas is assumed),and the unit length of the four velocity ( u µ u µ = − ρ , u , u i , B i ) as the 8 in-dependent variables where u i is the space component ofthe four velocity. B i can be written as ℜ i / α where α isthe lapse function ( α = p − /g tt ) and ℜ i is the magneticfield measured by the Fiducial observer (FIDO) whosefour velocity is n µ = ( − α, , , ∂ t ( √− gρu t ) = − ∂ i ( √− gρu i ) (1) ∂ t ( √− gT tν ) = − ∂ i ( √− gT iν ) + √− gT κλ Γ λνκ (2) ∂ t ( √− gB i ) = − ∂ j (cid:2) √− g ( b i u j − b j u i ) (cid:3) , (3)where T µν is the stress energy tensor that is composed ofthe sum of the matter part ( T µν Matter = ( ρ + u + p ) u µ u ν + pg µν ) and electromagnetic part ( T µν EM = F µα F να − g µν F αβ F αβ / √ π is absorbed into thedefinition of the Faraday tensor (Gammie et al. 2003). b µ is introduced so that Eq.(3) looks simple, andit is defined as b µ = ǫ µνκλ u ν F λκ where ǫ µνκλ =( − / √− g ) [ µνλκ ]. [ µνλκ ] is the completely antisymmet-ric symbol. In the fluid-rest frame, b µ becomes (0, B i ).In this study, we adopt the conservative scheme forintegration of the GRMHD equations. In this case, theleft terms of Eq.(1)-(3) are considered to be fundamentalvariables and called as the conserved variables. The rightterms of Eq.(1)-(3) are flux terms with a source term (thesecond right term of Eq.(2)).Since we have to estimate pressure of the fluid, wehave to estimate the primitive variables from the con-served variables at each time step. The problem is thatevelopment of GRMHD Code and Application to GRBs 3the primitive variables can not be expressed analyticallyby the conserved variables. Thus we have to use theNewton-Raphson method to obtain the primitive vari-ables from the conserved ones (Noble et al. 2006).Basically, we adopt the 2D scheme introduced by No-ble et al. (2006) to calculate the primitive variables.However, it sometimes happens that the 2D scheme failsto converge well, and the primitive variables can not beobtained precisely. In such a case, we first adopt the1D W scheme introduced by Noble et al. (2006) and seewhether the 1D W scheme converges. If it converges well,we adopt the primitive variables obtained by the 1D W scheme for the next time step. Otherwise, we adopt thesecond choice explained in the following subsection. Supplemental Method to Calculate PrimitiveVariables
Following Noble et al. (2006), we introduce convenientvariables v , W , Q µ , and ¯ Q µ . These variables, of course,depend on the primitive variables. The definition of thesevariables are: v = v i v i , W = ω Γ , Q µ = αT tµ , and ¯ Q µ = j µλ Q λ , where v i is the fluid velocity relative to FIDO, ω = ρ + u + p , Γ = 1 / √ − v , and j µλ = g µλ + n µ n λ . Itis apparent that Q µ and ¯ Q µ can be written analyticallyby the conserved variables. On the other hand, v and W can not be expressed analytically by the conservedvariables. Thus, we have to solve v and W numericallyin order to determine the proper, corresponding primitivevariables.Here we show that an upper limit and a lower limitfor W can be obtained before searching for a solution of W and v numerically. Thanks to this fact, all we haveto do is to seek the solution with the condition W min ≦ W ≦ W max . From Eq.(28) and Eq.(29) in Noble et al.(2006), W and v satisfy the following equations: v = ¯ Q W + ( Q µ ℜ µ ) ( ℜ + 2 W )( ℜ + W ) W (4) v = 2 ℜ (cid:20) ( Q µ ℜ µ ) W − W + p − ( Q µ n µ ) (cid:21) − . (5)From these equations and the relation 0 ≦ v < v and W satisfy the following relations: f ( W ) = W + 2 ℜ W + ( ℜ − ¯ Q ) W − Q µ ℜ µ ) W − ℜ ( Q µ ℜ µ ) ≧ g ( W ) = W + { ℜ + ( Q µ n µ ) − p } W −
12 ( Q µ ℜ µ ) ≦ h ( W ) = W + {ℜ + ( Q µ n µ ) − p } W −
12 ( Q µ ℜ µ ) ≧ . (8)Since f (0) ≦ f ′ (0) ≦
0, and at least one of the solutionfor f ′′ ( W ) = 0 is less than 0, there is only one positivesolution W a that satisfies f ( W a ) = 0. Thus, from Eq.(6), W has to be greater than W a .We can understand the behavior of g ( W ) from its firstderivative for W : g ′ ( W ) = W (cid:20) W + 2 { ℜ + ( Q µ n µ ) − p } (cid:21) . (9) -1 0 100.20.40.60.81 -1 0 105001000-1 0 1-0.200.20.40.60.81 -1 0 10.511.522.53 Fig. 1.—
Simulation of 1D shock tube test (Komissarov 1999).The state at t = 1 . − ≤ x ≤
2. The upperleft panel shows density, the upper right panel shows pressure, thelower left panel shows the velocity in the x-direction, and the lowerright panel shows the bulk Lorentz factor.
It is apparent that W = 0 is a solution for g ′ ( W ) = 0.As for the other solution(s), it is not so obvious becausethe pressure p depends on W and v . However, it willbe natural to consider that the monotonic relation holdsbetween W and p . It means that the pressure rises when W becomes larger. If this assumption is adopted, aslong as g ′ ( W ) = 0 has another solution, it is a positiveone W = W α ≧
0. This is because when W = 0, p should be also 0 and (cid:2) W + 2 { / ℜ + ( Q µ n µ ) − p } (cid:3) isa positive value. Thus, g ( W ) = 0 has only one positivesolution W b . This holds even if g ′ ( W ) = 0 has only onesolution at W = 0. Also, same conclusion can be derivedfor h ( W ): there is only one positive solution W c thatsatisfies h ( W c ) = 0.Since h ( W ) ≧ g ( W ), the relation W c ≦ W b holds.In conclusion, W has to be in the range W min =Max( W a , W c ) ≦ W ≦ W b = W max . Thus all we have todo is to find a solution of W that satisfies v = v in this range. This procedure is more expensive than the2D scheme and the 1D W scheme, but the solution for W and v is more likely to be found because the range forthe solution of W is determined apriori. Thus we use thismethod as a supplementary one to obtain the primitivevariables. TEST CALCULATIONS
Using the GRMHD code that is developed in thisstudy, we check whether it can pass many, well-knowntest calculations. The first three tests are special rel-ativistic hydrodynamic (SRHD) or special relativisticmagnetohydrodynamic (SRMHD) calculations, while therest of three tests are GRMHD ones.
Shock Tube Problems
S. Nagataki -1 0 10102030 -1 0 1-1000100200300-1 0 1-1-0.500.51 -1 0 1123456 -1 0 1-1-0.500.51 -1 0 1-20020
Fig. 2.—
Simulation of 1D collision test (Komissarov 1999). Thestate at t = 1 . − ≤ x ≤
2. The left panelsshow density, velocity in the x-direction, and bulk Lorentz factor(from top to bottom), while right panels show pressure, velocityin x-direction, and y-component of magnetic field (from top tobottom).
1D shock tube tests are the most basic test problemsfor SRHD/SRMHD. We have carried out a number ofthe test simulations introduced in Komissarov (1999)and Balsara (2001). Here we describe only two of them.One is the shock tube test1 and the other is the collisiontest (Komissarov 1999; Mizuno et al. 2006). The initialleft and right states are summarized in Table.1. Numberof grid points is 600 for both simulations. The resultsare shown in Fig.1 and Fig.2, which show that the testcalculations are well solved as in the previous studies.
Double Shock Problems
Here 2D shock tube problem is done to con-firm whether the shock dynamics in the multidimen-sional flow can be solved safely. This problem in-cludes the interactions of shocks, rarefactions, con-tact discontinuities. Initially a square computationaldomain is prepared in x-y plane and divided intofour quarter boxes. Initial condition in each boxis summarized in Table.2. This condition is samewith previous study (Del Zanna & Bucciantini 2001;Zhang & MacFadyen 2006; Mizuta et al. 2006). We use400 ×
400 uniform grid points in a square computationalbox. Boundary condition is open ones. Density contourat the final stage of the simulation is shown in Fig.3,which shows that our code reproduces the previous stud-ies very well.
Cylindrical Explosion Test
Here we go to a SRMHD test. A famous, cylindri-cal blast explosion test is done (Del Zanna et al. 2003;Leismann et al. 2005). We use the [0 , × [0 ,
1] Carte-sian grid with a resolution of N x = N y = 250 grid Fig. 3.—
Simulation of 2D shock tube problem. Density contourat t = 0 . × ≤ x ≤ ≤ y ≤ points. We define an initially static background with ρ = 1 . , p = 0 .
01, and B x = 4 .
0. The relativistic flowcomes out by setting a much higher pressure, p = 10 within a circle of radius r = 0 .
08 placed at the center ofthe domain. γ for the equation of state is set to be 4/3.Final time is set to be 0.4. The result is shown in Fig.4.The upper left panel shows the density contour in loga-rithmic scale. The upper right panel shows the pressurecontour in logarithmic scale. The lower left panel showscontour of the bulk Lorentz factor. The lower right panelshows the divergence of the magnetic fields in logarithmicscale with magnetic field lines. These results are consis-tent with the previous studies. Especially, the divergenceof the magnetic fields is kept as small as 10 − . Gammie’s Flow
Next we consider a GRMHD test. A steady,magnetized inflow solution on the equatorial planearound a Kerr BH is considered (Takahashi et al. 1990;Gammie 1999). Initially, the steady inflow solution forthe Kerr parameter a = 0 . F θφ = 0 . ≦ r ≦ . − − ≦ θ/π ≦ . − ]. The model is run for t = 1 .
5. Thephysical values at boundaries are fixed throughout thesimulation. Results are shown in Fig.5: density, ra-dial component of the 4-velocity, the φ component ofthe 4-velocity, and ℜ φ at the final stage of the simu-lation. When the initial state is written in the samefigure, we can see that the final state coincides withthe initial state. To show it more quantitatively, weintroduce the norms of the errors for these values asa function of the number ( N ) of grid points in the ra-dial coordinate. The definition of the norm of the errorevelopment of GRMHD Code and Application to GRBs 5 Fig. 4.—
A RMHD 2D cylindrical explosion test with a pressurejump as high as 10 . The state at t = 0 . × ≤ x ≤ ≤ y ≤
1. The upper left panel shows thedensity contour in logarithmic scale. The upper right panel showsthe pressure contour in logarithmic scale. The lower left panelshows contour of the bulk Lorentz factor. The lower right panelshows divergence of the magnetic fields in logarithmic scale withmagnetic field lines. is Σ i = Ni =1 | a (final) − a (initial) | / Σ i = Ni =1 | a (initial) | . In Fig.6,the norms of errors are shown. We can see that these val-ues converge roughly proportional to N − , as expected. Blandford-Znajek Monopole Solution
Further we continue to test the GRMHD code.We consider the Blandford-Znajek monopole solu-tion (Blandford & Znajek 1977). This analytic solu-tion has been investigated numerically by previousstudies (Komissarov 2004b; McKinney & Gammie 2004;Tanabe & Nagataki 2008).The computational domain is axisymmetric, with agrid that extends from r in = 0 . r + to r out = 230 andfrom θ = 0 to θ = π where r + = 1 + √ − a is the outerevent horizon. The numerical resolution is 300 × ℜ µ = − n ∗ ν F µν = (0 , α sin θ/ √− g, ,
0) in theKerr-Schild coordinate where ∗ F µν and g are the dualfield tensor and determinant of the Kerr-Schild met-ric. The plasma velocity relative to the FIDO is setto zero initially, and its pressure and density are setto small value ( P = ρ = ℜ / . B ≥ Γ ρ + ( γ Γ − ( γ − u is satisfied, densityand internal energy are increased by the same factor sothat the critical condition holds (Komissarov 2004b). γ is set to be 4/3. We have performed numerical simula-tions with the Kerr parameters 0, 0.01, 0.1, 0.2, 0.3, 0.4,0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99, and 0.995 until time t =200.The total energy flux, which is the integrated outgoing Fig. 5.—
Gammie’s equatorial inflow solution in the Kerr metricwith a = 0 . F θφ = 0 .
5. Numberof grid point is 1024. The state at t = 1 . φ component of the 4-velocity, and ℜ φ at the final stage of thesimulation. Boyer-Lindquist coordinate is used for the simulation. Fig. 6.—
Convergence results for the Gammie’s equatorial in-flow solution in the Kerr metric with a = 0 . F θφ = 0 .
5. Norms of the error for ρ , u r , u φ , and ℜ φ at the final stage of the simulation are shown in the figure. Thestraight line represents the slope expected for second-order conver-gence. The definition of the norm of the error is written in thetext. S. Nagataki
Fig. 7.—
Outgoing poynting fluxes as a function of the zenith an-gle for the Blandford-Znajek monopole solution with Kerr param-eter a = 0 . , . , . , . , . , . , . , . , . , . , .
9, and 0 . r = 20 and t = 200. Numbers of gridpoint are 300 × poynting flux over the zenith angle, can be written as˙ E = 2 π Z π dθ √− g ( − T rt ) = 2 π Z dx (cid:18) dθdx (cid:19) √− g ( − T rt )= 2 π Z dx F E , (10)where x = θ/π is introduced as a convenient vari-able (Gammie et al. 2003).In Fig.7, the outgoing poynting fluxes ( F E ) as a func-tion of zenith angle are shown. The fluxes are measuredat r = 20 at the final stage of the simulations. We wouldlike to note that the outgoing poynting flux hardly de-pends on the radius where it is evaluated. This meansthat the conservation of the outgoing poynting flux isconfirmed numerically.In Fig.8(a), we plot the total energy flux ( ˙ E ) at the fi-nal stage for small Kerr parameters (0 ≤ a ≤ .
2) by rect-angular points. Dashed line is just the interpolation ofthe calculated values. For comparison, the second-orderanalytical solution is shown by dotted line and the forth-order analytical solution is shown by solid line. Fromthis comparison, we can see that all of them coincidewith each other. Thus the results of the numerical simu-lations by the GRMHD code are confirmed by analyticalsolutions.The situation becomes different for large Kerr parame-ters. In Fig.8(b), we plot the same values with Fig. 8(a),but for wide range of the Kerr parameters (0 ≤ a ≤ Fig. 8.—
Upper panel (a): Comparison of the derived, conserved,total energy flux. Dashed line with rectangular points is numericalresult for small Kerr parameter (0 ≤ a ≤ . ≤ a ≤ t = 200. more reliable than the analytical estimation (see Tanabeand Nagataki (2008) for detailed discussion). Fishbone and Moncrief ’s Test
Here we present a final test of the GRMHD code. Asteady and stationary torus (Fishbone & Moncrief 1976;Abramowicz et al. 1978) around a Kerr BH that is sup-ported by both centrifugal force and pressure is solvednumerically. Of course, it should be solved as a steadyand stationary state.We have integrated a Fishbone-Moncrief solutionaround a Kerr BH with a = 0 .
9. We set u t u φ = 4 .
45 and R in = 6 .
0. The grid extends radially from r in = 1 .
40 to r out = 100. The same floors with Gammie et al. (2003)are used for ρ and u . The numerical resolution is N × N and the solution is integrated for t = 10. The resultingnorm of the error, which converges roughly proportionalto N − , is shown in Fig.9.Next we follow the time evolution of the Fishbone-Moncrief solution with magnetic fields. The vector po-tential, A φ ∝ max( ρ/ρ max − . ,
0) where ρ max is the peakdensity in the torus, is introduced (Gammie et al. 2003).The field is normalized so that the minimum value of p gas /p mag becomes 10 . The time integration extends for t = 2000. The number of grid points is 256 × r in = 1 .
40 to r out = 300while it extends in the zenith angle from θ = 0 to θ = π .The density contours in logarithmic scale (from 10 − to 10 ) are shown in Fig.10. These are projected on the(r sin θ ,r cos θ )-plane. The upper left panel shows theinitial state. The upper right panel shows the final stateof the simulation without magnetic fields. The lower leftevelopment of GRMHD Code and Application to GRBs 7 Fig. 9.—
Convergence results for the Fishbone-Moncrief problemfor a Kerr BH with a = 0 .
9. The straight line represents the slopeexpected for second-order convergence. Kerr-Schild coordinate isused for the simulation. panel shows the final state of the simulation with mag-netic fields. The lower right panel is same with the lowerleft one, but for a wide region. Due to the presence ofthe magnetic fields, the angular momentum in the torusis conveyed outward and the the torus starts to accrete,and the jet is launched from the BH around the polarregion. This result is consistent with the previous stud-ies (Gammie et al. 2003; McKinney & Gammie 2004;McKinney 2006a; Mckinney 2006b). SIMULATION OF COLLAPSAR
Since our code has passed the many test calculationsshown in the previous section, we now simulate the dy-namics of a collapsar using the code. However, we haveto say beforehand that no microphysics is included inthe code such as nuclear reactions, neutrino processes,and equation of state for dense matter. So this is theFIRST STEP of our project to simulate the dynamics ofa collapsar and formation of a relativistic jet of a GRB.
Method of Calculation
We have done a 2D GRMHD simulation of a col-lapsar using the Modified Kerr-Schild coordinate and G = c = M = 1 units. When we show results, thecoordinate is transfered from the Modified Kerr-Schildcoordinate to the Kerr-Schild one, and the units are fre-quently transfered to cgs units. The calculated regioncorresponds to a quarter of the meridian plane under theassumption of axisymmetry and equatorial symmetry.The spherical mesh with 256( r ) × θ ) grid points isused for all the computations. The calculated region cov-ers from r =1.8 to 3 × (that corresponds to 5.3 × cmand 8.9 × cm in cgs units, as explained below) withuniform grids in the Modified Kerr-Schild space.We adopt the model 12TJ in Woosley and Heger(2006). This model corresponds to a star that has 12 M ⊙ Fig. 10.—
Density contour in logarithmic scale (from 10 − to10 ) for the Fishbone-Moncrief problem for the Kerr parameter a = 0 .
9. Number of grid points is 256 × t = 2000. The upper left panel shows the initialstate. The upper right panel shows the final state of the simulationwithout magnetic fields. The lower left panel shows the result withmagnetic fields. The lower right panel is same with the lower leftone, but for a wide region. These results are projected on the (rsin θ ,r cos θ )-plane. initially with 1% of solar metallicity, and rotates rapidlyand does not lose its angular momentum so much byadopting small mass loss rate. As a result, this starhas a relatively large iron core of 1 . M ⊙ , and rotatesrapidly (the estimated Kerr parameter that a BH formingof mass and angular momentum of the inner 3 M ⊙ wouldformally have is 0.57) at the final stage. Of course, whatkind of stars are appropriate for progenitors of GRBs isstill under debate (Yoon et al. 2006). Thus we chosedthe model 12TJ as a first example of our study becausethe iron core is large and rotating rapidly, which seemsto form a rapidly-rotating BH, among the models listedin Woosley and Heger (2006). We assume that the cen-tral part of the star with 2 M ⊙ has collapsed and formeda BH at the center with the Kerr parameter a = 0 . M = 2 M ⊙ , r = 1 corresponds to 2.95 × cm, as explained above.Also, the inner boundary r = 1 . r + = 1 + √ − a = 1 . u θ = 0 initially. As for u φ , we extrapolate its value suchas u φ ( r, θ ) = u φ ( r, π/ × sin θ .Effects of magnetic fields are taken into account inthe model 12TJ. However, again, since 1-D calculation isdone, we do not know the configuration of the magneticfields. It is difficult to extrapolate magnetic fields thatsatisfy the condition div B = 0 everywhere. Also, thereare much uncertainty on the amplitude of the magnetic S. Nagataki Fig. 11.—
Contours of rest mass density at the central regionin logarithmic scale, in which cgs units are used assuming thatthe gravitational mass of the BH is 2 M ⊙ . The length unit in thevertical/horizontal axes corresponds to 2.95 × cm. Upper panel(a) shows the state at t = 110000 (that corresponds to 1.0835 sec),while lower panel (b) shows the one at t = 180000 (that correspondsto 1.773 sec). These results are projected on the (r sin θ ,r cos θ )-plane. fields in a progenitor. Thus we do not use the informa-tion on magnetic fields of the model 12TJ. Rather, weadopt the same treatment in section 3.6. That is, thevector potential A φ ∝ max( ρ/ρ max − . ,
0) sin θ where ρ max is the peak density in the progenitor (after extract-ing the central part of the progenitor that has collapsedand formed a BH). The field is normalized so that theminimum value of p gas /p mag becomes 10 . The defini-tion of p mag is p mag = b /
2. The reason why we adoptthe strong dependence on the zenith angle for A φ is soas not to suffer from discontinuity of magnetic fields atthe polar axis. The resulting biggest amplitude of themagnetic fields is 7.4 × G at r = 950 (2.8 × cm).We use a simple equation of state p gas = ( γ − u where we set γ =4/3 so that the equation of state roughlyrepresents radiation gas.As for the boundary condition in the radial direction, Fig. 12.—
Same with Fig.11(b), but for a wider region. we adopt the outflow boundary condition for the innerand outer boundaries (Gammie et al. 2003). As for theboundary condition in the zenith angle direction, axisof symmetry condition is adopted for the rotation axis,while the reflecting boundary condition is adopted for theequatorial plane. As for the magnetic fields, the equato-rial symmetry boundary condition, in which the normalcomponent is continuous and the tangential componentis reflected, is adopted.
Results
In Fig.11, color contours of rest mass density at thecentral region are shown. Colors represent the densityin units of g cm − in logarithmic scale. These resultsare projected on the (r sin θ ,r cos θ )-plane. The length r = 200 corresponds to 5.9 × cm. The time unitcorresponds to 9.85 × − sec. Upper panel (a) repre-sents the contours of rest mass density at t = 110000(that corresponds to 1.0835 sec), while lower panel showsthe contours at t = 180000 (that corresponds to 1.773sec). Fig.12 is the same figure with Fig.11(b), but for awider region. A jet is clearly seen along the rotation axis.In Fig.13, mass accretion rate history on the horizon isshown. The definition of the mass accretion rate is˙ M = 2 × π Z θ dθ √− gρu r . (11)It takes about 0.15 sec for the inner edge of the matterto reach the horizon. When the matter reaches there,there is an initial spike of the mass accretion rate. Afterthat, there is a quasi-steady state like Fig.11(a) is real-ized. Then, the jet is launched at ∼ p gas /p mag )in logarithmic scale at t = 180000 in Fig.14. As ex-pected, the plasma beta is low in the jet region while itis high in the accretion disk region. We show color con-tours of bulk Lorentz factor around the central regionat t = 180000 in Fig.15(a) (upper panel, in logarithmicscale). Color contours of the energy flux per unit rest-evelopment of GRMHD Code and Application to GRBs 9 Fig. 13.—
Mass accretion rate history on the horizon. The unit M ⊙ sec − is used assuming that the gravitational mass of the BHis 2 M ⊙ throughout the calculation. mass flux ( E = − T rt / ( ρu r )), which is conserved for aninviscid fluid flow of magnetized plasma, are also shownin Fig.15(b) (lower panel, in logarithmic scale). Thisvalue represents the bulk Lorentz factor (Γ ∞ ) of the in-vischid fluid element when all of the internal and mag-netic energy are converted into kinetic energy at largedistances (McKinney 2006a). We can see that the bulkLorentz factor of the jet is still low (Fig.15(a)), but it canbe potentially as high as 10 at large radius (Fig.15(b)).At t = 180000, the strength of the magnetic field ( √ πb )at the bottom of the jet is found to be ∼ G, and u φ / u t is ∼ . r ms on the equatorial plane. Here r ms is the marginally stable orbit. For the Kerr BH with a = 0 . r ms is 4.23. As stated in section 4.1, the ini-tial biggest amplitude of the magnetic fields is 7.4 × Gat r = 950 where the initial density is ∼ g cm − ,the expected amplification factor of the magnetic fieldsdue to the gravitational collapse and differential rota-tion around the BH is ( ρ/ρ ) / × ( d Φ /dt/ π ) ∗ ∆ t ∼ × × ∼ × . Thus the initial mag-netic field can be amplified as large as several times of10 G, which is roughly consistent with the amplitudeof the magnetic fields at the bottom of the jet. At latephase, the magneto-rotational instability (MRI) may bealso working, which is discussed in the next section.In Fig.16, contours of the φ component of the vec-tor potential ( A φ ) at t = 180000 are shown. Levelsurfaces coincide with poloidal magnetic field lines, andfield line density corresponds to poloidal field strength.As expected, the magnetic fields are strong at the jetregion, which makes the plasma beta very low. FromFig.16, the opening angle of the jet is estimated as5 ◦ − ◦ . From Fig.14, 15(b), and 16, this jet should corre-spond to the poynting flux jet (Hawley & Krolik 2006).This jet is surrounded by the funnel-wall jet re-gion (Hawley & Krolik 2006), which is shown in Fig.17 Fig. 14.—
Contour of the plasma beta ( p gas /p mag ) at t = 180000in logarithmic scale. below.In Table.3 and 4, the integrated energies of matter andelectromagnetic field at t = 180000 are shown. The in-tegrated region is between the horizon and r = 200 (forTable.3) or r = 40 (for Table.4), and within the zenithangle measured from the polar axis. As for the mattercomponent, the contribution of the rest mass energy issubtracted. That is, E Matter = 2 × π Z r =200 or 40 r + dr Z θ dθ √− g ( T , Matter − ρu u ) . (12)Factor 2 is coming from the symmetry of the system withrespect to the equatorial plane. The field part can bewritten as E EM = 4 π Z r =200 or 40 r + dr Z θ dθ √− gT , EM . (13)It can be seen that the energy in electromagnetic fielddominates that in matter within r ≤
40, while they be-come comparable within r ≤ ∼ × erg; Frail et al. 2001).Finally, we show the rest-mass density, outgoing massflux, and outgoing poynting flux in Fig.17. The toppanel (a) shows the rest-mass density (g cm − ) as afunction of the zenith angle at r = 10 r ms = 42 . t = 180000. It is seen that the low-density regionis realized around the polar axis, which correspondsto the jet region (0.1 radian corresponds to 5.7 ◦ ).The middle panel (b) shows the outgoing mass flux ρu r (g cm − s − ) at r = 10 r ms and t = 180000. Itis seen that the outgoing mass flux exists around0 . ≤ θ ≤ .
14, which corresponds to the funnel-wall jet (De Villiers et al. 2003; Hirose et al. 2004;McKinney & Gammie 2004; Kato et al. 2004;De Villiers et al. 2005; Hawley & Krolik 2006). Thebottom panel (c) shows the outgoing poynting flux (inunits of 10 erg s − rad − ) at r = r + and t = 180000.Definition of the outgoing poynting flux is F BZ = − × π √− gT tr, EM . (14)0 S. Nagataki Fig. 15.—
Upper panel (a): contours of bulk Lorentz factoraround the central region at t = 180000. Lower panel (b): con-tours of the energy flux per unit rest-mass flux at t = 180000 thatrepresent the bulk Lorentz factor (Γ ∞ ) of the invischid fluid ele-ment when all of the internal and magnetic energy are convertedinto kinetic energy at large distances. The contours are written inlogarithmic scale. Factor 2 is coming from the assumption of the symmetryof the system with respect to the equatorial plane. Wecan see that the positive outgoing poynting flux exists atthe jet region (0 ≤ θ ≤ .
23 in radian). The integratedenergy of the outgoing poynting flux at r = r + and t =180000 is 4 . × erg s − . Since the duration of thejet in this study is ∼ . SUMMARY AND DISCUSSION
In order to investigate the formation of relativistic jetsat the center of a progenitor of a GRB, we have devel-oped a two-dimensional GRMHD code. In order to con-firm the reliability of the code, we have shown that thecode passes many, well-known test calculations. Thenwe have performed a numerical simulation of a collap-sar using a realistic progenitor model. We have followed
Fig. 16.—
Contours of the φ component of the vector potential( A φ ) at t = 180000. Level surfaces coincide with poloidal mag-netic field lines, and field line density corresponds to poloidal fieldstrength. The biggest amplitude of the magnetic fields is ∼ G. the time evolution of the system for 1.773 s, and it wasshown that a jet is launched from the center of the pro-genitor. We also found that the mass accretion rate isin quasi-stable state before the launch of the jet, while itshows rapid time variability that resembles to a typicaltime profile of a GRB after the launch. The structureof the jet is similar to the previous study: a poyntingflux jet is surrounded by a funnel-wall jet. Even at thefinal stage of the simulation, the bulk Lorentz factor ofthe jet is still low, and the total energy of the jet is stillas small as 10 erg. However, we found that the en-ergy flux per unit rest-mass flux ( E = − T rt / ( ρu r )) is ashigh as 10 at the at the bottom of the jet. Thus weconclude that the bulk Lorentz factor of the jet can bepotentially high when it propagates outward. Also, aslong as the duration of the activity of the central engineis long enough, the total energy of the jet can be largeenough to explain the typical explosion energy of a GRB( ∼ × erg). It is shown that the outgoing poyntingflux exists at the horizon around the polar region, whichproves that the Blandford-Znajek mechanism is working.However, we conclude that the jet is launched mainly bythe magnetic field amplified by the gravitational collapseand differential rotation around the BH, rather than theBlandford-Znajek mechanism in this study.When we apply the Blandford-Znajek for-mula (Barkov & Komissarov 2008b), the integratedoutgoing poynting flux is˙ E = 3 . × f ( a )Ψ M − erg s − (15)where M = M BH / M ⊙ = 1, Ψ = Ψ / G cm ,and f ( a ) = a / (1 + √ − a ) = 0 . . × B erg s − for the jet with openingangle θ = 5 ◦ , which is comparable to our numerical resultevelopment of GRMHD Code and Application to GRBs 11 Fig. 17.—
The top panel (a): the rest-mass density (g cm − ) as afunction of the zenith angle at r = 10 r ms = 42 .
3. The middle panel(b): the outgoing mass flux ρu r (g cm − s − ) at r = 10 r ms . Thebottom panel (c): the outgoing poynting flux (in units of 10 ergs − rad − ) at r = r + . These values are evaluated at the final stageof the simulation ( t = 180000). (Fig.17(c)).As for the efficiency of converting the released gravi-tational energy to the jet’s energy, it can be estimatedas follows: the mass accretion rate is ∼ . M ⊙ s − (Fig.13), the total energy of the jet at the final stageis ∼ erg (Table.3), and the duration of the jet is ∼ . ∼ − . When we use the outgoing poynting flux atthe horizon (4.6 × erg s − ), the efficiency is as low as ∼ × − . These values seem to be very small com-pared with the previous study (De Villiers et al. 2005;Mckinney & Narayan 2007a). One of the reason will bebecause they used an almost steady disk model. On theother hand, we used a realistic progenitor model thatcollapses gravitationally. Thus the resulting mass accre-tion rate is pretty high. Second reason may be becausethe efficiency is still low even at the final stage of thesimulation. If we perform the simulation further, the ef-ficiency may become higher with time: mass accretionrate will become smaller, and the jet energy might belarger due to the amplification of the magnetic fields dueto winding-up (and MRI) effects. Also, when the initialamplitude of the magnetic field is set to be larger (as inBarkov and Komissarov 2008), the efficiency may be en-hanced. Further, we should investigate the dependenceof the dynamics on progenitor models as well as the Kerrparameter of the BH. We are planning to investigate thispoint systematically in the next paper.It is well known that the system is unsta-ble against MRI when there is a strong nega-tive shear profile ( d Ω /d ln r ) (Balbus & Hawley 1991;Balbus & Hawley 1994), where Ω is the angular velocity.The saturation toroidal magnetic field strength is roughly expected to be B φ ∼ (4 πρ ) / r Ω (Akiyama et al. 2003;Akiyama & Wheeler 2005), which is confirmed by semi-global simulations (Obergaulinger et al. 2008). The sat-uration poloidal magnetic field strength is roughly anorder of magnitude smaller (Obergaulinger et al. 2006).Thus B φ may be amplified by MRI as strong as 1 . × G ρ r ms Ω where Ω is Ω / rad s − . The char-acteristic timescale for saturating the MRI is the Alfv´encrossing time: t A ∼ . R ρ / B − where R , ρ ,and B are the radius in units of 10 cm, the density inunits of 10 g cm − , and the amplitude of magnetic fieldsin units of 10 G, respectively. Thus this characteristictimescale can be shorter than the winding-up timescalefor strong magnetic fields. However, the length scale ofthe mode with the largest MRI growth rate is approx-imately λ MRI ∼
50 cm P . B ρ − / where P . is therotation period in units of 0.5 ms, which is too shortto be resolved numerically. At least, it is not resolvedin the beginning of the simulations. After the magneticfield is amplified to a certain value due to gravitationalcollapse and winding-up effect, MRI may be working inthis study (Obergaulinger et al. 2006b; Ott et al. 2006;Burrows et al. 2007; Dessart et al. 2008). It will be nec-essary to develop a sophisticated code that takes intoaccount the MRI effectively with help of semi-global sim-ulations (Obergaulinger et al. 2008) in order to evaluatethe influence of MRI on the dynamics of a collapsar.It is well known that it becomes difficult toobtain the matter part of the primitive vari-ables ( ρ , u , u i ) precisely by the Newton-Raphsonmethod (Noble et al. 2006) due to numerical trun-cation errors (Komissarov 2002; Komissarov 2004a;Komissarov 2004b; McKinney & Gammie 2004;Komissarov 2005; Mckinney 2006b) when the elec-tromagnetic part of the stress energy tensor ( T µν EM )greatly exceeds the matter part ( T µν Matter ). The problemis that the time integration of the electromagnetic partdoes not become so reliable, either. This is because theMHD condition ( u µ F µν = 0) is implicitly assumed inthe basic equations, and the resulting basic equationof electromagnetic part depends on the velocity offluid (Eq.3). Such pathological conditions may berealized at the bottom of the jet in our study. Inorder to confirm the validity of our results in thisstudy, we are planning to develop a general relativisticforce-free code that is coupled with the GRMHD codesophisticatedly (Mckinney 2006b).It is very important to evaluate the terminal bulkLorentz factor, because GRBs are considered to beemissions from relativistic flows with their bulk Lorentzfactors greater than 10 (e.g. (Lithwick & Sari 2001)).Although an ad-hoc thermal (and kinetic) energy depo-sition into the polar region seems to lead to relativisticjets with bulk Lorentz factors ∼
100 (Aloy et al. 2000;Aloy et al. 2002; Zhang et al. 2003; Zhang et al. 2004;Cannizzo et al. 2004; Mizuta et al. 2006;Mizuta & Aloy 2008; Morsony et al. 2007;Wang et al. 2008), it is still controversial whethersuch ad-hoc energy deposition is justified bynumerical simulations with proper neutrinophysics (Nagataki et al. 2007). On the otherhand, numerical study on the acceleration of elec-tromagnetically powered jet requires quite high-2 S. Nagatakiresolution (Komissarov et al. 2007; Narayan et al. 2007;Tchekhovskoy et al. 2008; Komissarov et al. 2008).Due to the reason, a simplified jet model withan idealized boundary condition is used atpresent in order to investigate whether the ini-tial poynting flux can be effectively convertedinto kinetic energy (Tchekhovskoy et al. 2008;Komissarov et al. 2008). According to their results, aslong as confinement of the jet is realized, accelerationoperates over several decades in radius and considerablefraction of the poynting flux can be converted intothe kinetic energy. Thus, from the high-ratio of thepoynting flux relative to rest-mass flux seen in our study(Fig.15(b)), a relativistic jet with high bulk Lorentzfactor may be realized at large radius.It is true that the two-dimensional restriction canbe a significant limitation. First, anti-dynamo theo-rem (Moffat 1978) prevents the indefinite maintenanceof the poloidal magnetic field in the face of dissipation.Second, axisymmetric simulations tend to overem-phasize the channel mode (Hawley & Balbus 1992),which produces coherent internal magnetized flowsrather than the more generic MHD turbulence. Hy-drodynamic instability in the azimuthal direction maybe also very important (Nagakura & Yamada 2008;Nagakura & Yamada 2009). Thus we are plan-ning to develop a three-dimensional GRMHDcode (De Villiers et al. 2003; Hirose et al. 2004;De Villiers et al. 2005; Hawley & Krolik 2006;Beckwith et al. 2008; Shafee et al. 2008;Mckinney & Blandford 2008) and investigate thedifference between two-dimensional simulations ofcollapsars with three-dimensional ones.In this study, we assumed that the central regionof the progenitor has collapsed and a BH is formedat the center with surrounding envelope unchanged.Thus we solved the GRMHD equation on a fixedbackground. But the final goal of our project is tostudy how a GRB is formed from the gravitationalcollapse of a massive star. Thus we are planning todevelop a GRMHD code on a dynamical background,which makes the study on the gravitational collapseand BH formation at the center of a massive starpossible (Shibata 2003; Sekiguchi & Shibata 2004;Sekiguchi & Shibata 2005; Baiotti et al. 2005;Duez et al. 2006; Sekiguchi & Shibata 2007).In this study, photo-disintegration of nuclei andneutrino processes are not taken into account. Photo-disintegration absorb considerable amount of thermalenergy, and cooling/heating due to neutrino processesmay have great influence on the dynamics of a collap-sar (Di Matteo et al. 2002; Kohri & Mineshige 2002;Nagataki et al. 2003a; Surman & McLaughlin 2004;Lee et al. 2005; Gu et al. 2006; Nagataki et al. 2007;Kawanaka & Mineshige 2007; Kawabata et al. 2008;Rossi et al. 2008; Zhang & Dai 2009;Cannizzo & Gehrels 2009). Especially, pair-annihilation of electron-type neutrinos may be akey process to drive a GRB jet (Woosley 1993;MacFadyen & Woosley 1999; Asano & Fukuyama 2000; Asano & Fukuyama 2001; Miller et al. 2003;Surman & McLaughlin 2005; Kneller et al. 2006;Shibata et al. 2007; Birkl et al. 2007). Thus we areplanning to include such microphysics in our code, andperform more realistic simulations of collapsars.The SNe associated with GRBs often show pe-culiar properties. Some are very energetic andblight (Galama et al. 1998; Iwamoto et al. 1998;Hjorth et al. 2003; Malesani et al. 2004), but oth-ers prohibit such blight SNe from being accom-panied (Fynbo et al. 2006; Della Valle et al. 2006;Gal-Yam et al. 2006). Since the brightnessof SNe depends on the mass of produced Ni (Woosley et al. 1999; Nakamura et al. 2001),it is suggested that there is a huge variety of the amountof Ni in a SN that associates with a GRB. It is stillcontroversial where and when Ni is produced in aSN accompanied by a GRB (Nagataki et al. 2003b;Nagataki et al. 2006). It may be produced in aGRB jet (Maeda et al. 2002; Maeda & Nomoto 2003;Tanaka et al. 2007; Maeda et al. 2008; Tominaga 2009;Maeda & Tominaga 2009; Bucciantini et al. 2009), orit may be produced in (or outflow from) the accretiondisk around the BH (MacFadyen & Woosley 1999;Pruet et al. 2004; Fujimoto et al. 2004;Surman et al. 2006; Hu & Peng 2008), or itmay be synthesized around a proto-neutronstar (Uzdensky & MacFadyen 2007). At present, itis impossible to investigate the explosive nucleosynthesisin a collapsar in our code because nuclear reactions arenot taken into account. We are planning to includethis effect, and study the site of Ni production.Also, study of a GRB as a possible site where veryheavy elements and light elements are synthesizedis very important (Lemoine 2002; Beloborodov 2003;Suzuki & Nagataki 2005).This research was supported by Grant-in-Aid for Sci-entific Research on Priority Areas No. 19047004 by Min-istry of Education, Culture, Sports, Science and Technol-ogy (MEXT), Grant-in-Aid for Scientific Research (S)No. 19104006 by Japan Society for the Promotion ofScience (JSPS), and Grant-in-Aid for young Scientists(B) No.19740139 by Japan Society for the Promotionof Science (JSPS). The computation was carried out onNEC SX-8 at Yukawa Institute for Theoretical Physics(YITP) in Kyoto University and Cray XT4 at Centerfor Computational Astrophysics (CfCA) in National As-tronomical Observatory of Japan (NAOJ). A portion ofthis work was carried out while the author was at theKavli Institute for Particle Astrophysics and Cosmology(KIPAC). I am grateful to Roger Blandford, the direc-tor of KIPAC, for continuous encouragement, useful dis-cussions, and warm hospitality at KIPAC. I also thankWeiqun Zhang, Jonathan C. McKinney, Dmitri A. Uz-densky, Ruben Krasnopolsky for useful discussions. Iwish to thank all of the colleagues in YITP and KIPACfor useful discussions and kind support. Finally, I wouldlike to thank my family for continuous, warm support.
REFERENCESAbramowicz, M., Jaroszinski, M., Sikora, M. 1978, A&A, 63, 221 Aloy, M.A., M¨uller, E., Ib´a˜nez, J.M., Mrt´i, J.M., MacFadyen, A.2000, ApJ, 531, L119 evelopment of GRMHD Code and Application to GRBs 13
Aloy, M.A., Ib´a˜nez, J.M., Miralles, J.-A., Urpin, V. 2002, A&A,396, 693Anninos, P., Fragile, P.C., Salmonson, J.D. 2005, ApJ, 635, 723Ant´on, L., Zanotti, O., Miralles, J.A., et al. 2006, ApJ, 637, 296Akiyama, S., Wheeler, J.C., Meier, D.L., Lichtenstadt, I. 2003,ApJ, 584, 954Akiyama, S., Wheeler, J.C. 2005, ApJ, 629, 414Asano, K., Fukuyama, T. 2000, ApJ, 531, 949Asano, K., Fukuyama, T. 2001, ApJ, 546, 1019Baiotti, L. et al. 2005, Phys. Rev. D, 71, 024035Balbus, S.A., Hawley, J.F. 1991, ApJ, 376, 214Balbus, S.A., Hawley, J.F. 1994, MNRAS, 266, 769Balsara, D. 2001, ApJS, 132, 83Barkov, M., Komissarov, S.S. 2008, MNRAS, 385, L28Barkov, M., Komissarov, S.S. 2008, arXiv:0809.1402Beckwith, K., Hawley, J.F., Krolik, J.H. 2008, ApJ, 678, 1180Beloborodov, A.M. 2003, ApJ, 588, 931Bethe, H.A. 1990, Rev. Mod. Phys., 62, 801Birkl, R., Aloy, M.A., Janka, H.-Th., M¨uller, E. 2007, A&A, 463,51Blandford, R.D., Znajek, R.L. 1977, MNRAS, 179, 433Blandford, R.D., Payne, D.G. 1982, MNRAS, 199, 883Bucciantini, N., Quataert, E., Arons, J., Metzger, B.D.,Thompson, T.A. 2008, MNRAS, 383, 25Bucciantini, N., Quataert, E., Metzger, B.D., Thompson, T.A.,Arons, J., Del Zanna, L. 2009, arXiv:0901.3801Burrows, A., Dessart, L., Livne, E., Ott, C.D., Murphy, J. ApJ,664, 416Cannizzo, J.K., Gehrels, N., Vishniac, E.T. 2004, ApJ, 601, 380Cannizzo, J.K., Gehrels, N. 2009, arXiv:0901.3564De Villiers, J.-P., Hawley, J.F. 2003, ApJ, 589, 458De Villiers, J.-P., Hawley, J.F., Krolik, J.H. 2003, ApJ, 599, 1238De Villiers, J.-P., Hawley, J.F., Krolik, J.H., Hirose, S. 2005, ApJ,620, 878Del Zanna, L., Bucciantini, N. 2002, A&A, 390, 1177Del Zanna, L., Bucciantini, N., Londrillo, P. 2003, A&A, 400, 397Del Zanna, L., Zanotti, O., Bucciantini, N., Londrillo, P. 2007,A&A, 473, 11Della Valle, M., et al. 2006, Nature, 444, 1050Dessart, L., Burrows, A., Livne, E., Ott, C.D. 2008, ApJ, 673, L43Di Matteo, T., Perna, R., Narayan, R. 2002, ApJ, 579, 706Duez, M.D., Liu, Y.T., Shapiro, S.L., Shibata, M., Stephens, B.C.2006 Phys. Rev. Lett., 96, 031101Fishbone, L.G., Moncrief, V. 1976, ApJ, 207, 962Frail, D.A., et al. 2001, ApJ, 562, L55Fryer, C., M´esz´aros, P. 2000, ApJ, 588, L25Fujimoto, S., Hashimoto, M., Arai, K., Matsuba, R. 2004, ApJ,614, 847Fujimoto, S., et al. 2006, ApJ, 644, 1040Fynbo, J.P.U., et al. 2006, Nature, 444, 1047Galama, T.J., et al., 1998, Nature, 395, 670Gal-Yam, A., et al. 2006, Nature, 444, 1053Gammie, C.F. 1999, ApJ, 522, L57Gammie, C.F., McKinney, J.C., T´oth, G. 2003, ApJ, 589, 444Gu, W.-M., Liu, T., Lu, J.-F. 2006, ApJ, 643, L87Harten, A., Lax, P.D., van Lerr, B. 1983, SIAM Rev. 25, 35Hawley, J.F., Balbus, S.A. 1992, ApJ, 400, 595Hawley, J.F., Krolik, J.H. 2006, ApJ, 641, 103Hirose, S., Krolik, J.H., De Villiers, J.-P., Hawley, J.F. 2004, ApJ,606, 1083Hjorth, J., et al. 2003, Nature, 423, 847Hu, T., Peng, Q. 2008, ApJ, 681, 96Iwamoto, K., et al., 1998, Nature, 395, 672Kato, Y., Mineshige, S., Shibata, K. 2004, ApJ, 605, 307Kawabata, R., Mineshige, S., Kawanaka, N. 2008, ApJ, 675, 596Kawanaka, N., Mineshige, S. 2007, ApJ, 662, 1156Klebesadel, R.W., Strong, I.B., Olson, R.A. 1973, ApJ, 182, L85Kneller, J.P., McLaughlin, G.C., Surman, R.A. 2006, J. Phys. G.,32, 443Kohri, K., Mineshige, S. 2002, ApJ, 577, 311Komissarov, S.S. 1999, MNRAS, 303, 343Komissarov, S.S. 2002, MNRAS, 336, 759Komissarov, S.S. 2004b, MNRAS, 350, 427Komisarrov, S.S. 2004a, MNRAS, 350, 1431Komissarov, S.S. 2005, MNRAS, 359, 801Komissarov, S.S., Barkov, M.V., Vlahakis, N., K¨onigl, A. 2007,MNRAS, 380, 51 Komissarov, S.S., Barkov, M.V. 2007, MNRAS, 382, 1029Komissarov, S.S., Vlahakis, N., K¨onigl, A., Barkov, M.V. 2008,arXiv:0811.1467Lee, W.H., Rmirez-Ruiz, E., Page, D. 2005, ApJ, 632, 421Leismann, T., Ant´on , L., Aloy, M.A., M¨ u ller, E., Mart´i, J.M.,Miralles, J.A., Ib´a˜nez, J.M. 2005, A&A, 436, 503lemoine, M. 2002, A&A, 390, L31Lithwick, Y., Sari, R. 2001, ApJ, 555, 540MacFadyen, A.I., Woosley, S.E., ApJ, 1999, 524, 262Maeda, K., et al., 2002, ApJ, 565, 405Maeda, K., Nomoto, K. 2003, ApJ, 598, 1163Maeda, K., et al. 2008, Science, 319, 1220Maeda, K., Tominaga, N. 2009, arXiv:0901.0410Malesani, D., et al. 2004, ApJ, 609, L5McKinney, J.C., Gammie, C.F. 2004, ApJ, 611, 977McKinney, J.C. 2006a, MNRAS, 368, 1561McKinney, J.C. 2006b, MNRAS, 367, 1797McKinney, J.C., Narayan, R. 2007a, MNRAS, 375, 513McKinney, J.C., Narayan, R. 2007b, MNRAS, 375, 531McKinney, J.C., Blandford, R. 2008, arXiv:0812.1060Melesani, D., et al. 2004, ApJ, 609, L5Miller, W., George, N.D., Kheyfets, A., McGhee, J.M. 2003, ApJ,583, 833Misner, J.C., Thorne, K., Wheeler, J.A., Gravitation (W.H.Freeman and Co., San Francisco, 1970)Mizuno, Y., Yamada, S., Koide, S., Shibata, K. 2004a, ApJ, 606,395Mizuno, Y., Yamada, S., Koide, S., Shibata, K. 2004b, ApJ, 615,389Mizuno, Y., Nishikawa, K., Koide, S., Hardee, P., Fishman, G.J.2006, astro-ph/0609004Mizuta, A., Yamasaki, T., Nagataki, S., Mineshige, S. 2006, ApJ,651, 960Mizuta, A., Aloy, M.A. 2008, arXiv:0812.4813Moffat, K. 1978, Magnetic Field Generation in electricallyConducting Fluids (Cambridge: Cambridge Univ. Press)Morsony, B.J., Lazzati, D., Begelman, M.C. 2007, ApJ, 665, 569Nagakura, H., Yamada, S. 2008, ApJ, 689, 391Nagakura, H., Yamada, S. 2009, arXiv:0901.4053Nagataki, S., Kohri, K., Ando, S., Sato, K. 2003a, APh, 18, 551Nagataki, S., Mizuta, A., Yamada, S. Takabe, H., Sato, K. 2003b,ApJ, 596, 401Nagataki, S., Mizuta, A., Sato, K. 2006, ApJ, 647, 1255Nagataki, S., Takahashi, R., Mizuta, A., Takiwaki, T. 2007, ApJ,659, 512Nakamura, T., et al. 2001, ApJ, 555, 880Narayan. R., McKinney, J.C., Farmer, A.J. 2007, MNRAS, 375,548Noble, S.C., Gammie, C.F., McKinney, J.C., Del Zanna, L. 2006,ApJ, 641, 626Obergaulinger, M., Aloy, M.A., M¨uller, E. 2006, A&A, 450, 1107Obergaulinger, M., Aloy, M.A., Dimmelmeier, H., M¨uller, E.2006, A&A, 457, 209Obergaulinger, M., Cerd´a-Dur´an, P., M¨uller, E., Aloy, M.A. 2008,arXiv:0811.1652Ott, C.D., Burrows, A., Thompson, T.A., Livne, E., Walder, R.2006, ApJS, 164, 130Proga, D., MacFadyen, A.I., Armitage, P.J., Begelman, M.C.2003, ApJ, 599, L5Proga, D., Begelman, M.C. 2003, ApJ, 582, 69Proga, D. 2005, ApJ, 629, 397Pruet, J., Thompson, T.A., Hoffman, R.D. 2004, ApJ, 606, 1006Rossi, E.M., Armitage, P.J., Menou, K. 2008, MNRAS, 391, 922Skiguchi, Y., Shibata, M. 2004, Phys. Rev. D, 70, 084005Skiguchi, Y., Shibata, M. 2005, Phys. Rev. D, 71, 084013Skiguchi, Y., Shibata, M. 2007, PThPh, 117, 1029Shafee, R., et al. 2008, ApJ, 687, L25Shibata, M. 2003, Phys. Rev. D, 67, 024033Shibata, M., Liu, Y.T., Shapiro, S.L., Stephens, B.C. 2006,Phys. Rev. D, 74, 104026Shibata, M., Sekiguchi, Y., Takahashi, R. 2007, PThPh, 118, 257Surman, R., McLaughlin, G.C. 2004, ApJ, 603, 611Surman, R., McLaughlin, G.C. 2005, ApJ, 618, 397Surman, R., McLaughlin, G.C., Hix, W.R. 2006, ApJ, 643, 1057Suwa, Y., Takiwaki, T., Kotake, K., Sato, K. 2007, PASJ, 59, 771Suzuki, T.K., Nagataki, S. 2005, ApJ, 628, 914 TABLE 1Initial Conditions for Shock Tube Test1 and Collision Test
Test Type ρ p v x v y v z B x B y B z Shock Tube Test1 left state 1 1000 0 0 0 1 0 0right state 0.1 1 0 0 0 1 0 0Collision Test left state 1 1 5/ √
26 0 0 10 10 0right state 1 1 -5/ √
26 0 0 10 -10 0
Note . — γ for the equation of state is set to be 4/3. Final time is set to be 1.0 for Shock Tube Test 1, and 1.2 for Collision Test. TABLE 2Initial Conditions for 2D Shock Tube Problem
Region x y ρ p v x v y A 0 ≦ x ≦ ≦ y ≦ ≦ x ≦ ≦ y ≦ ≦ x ≦ ≦ y ≦ ≦ x ≦ ≦ y ≦ Note . — γ for the equation of state is set to be 5/3. Final time is set to be 0.4. TABLE 3Integrated Energies of Matter and Field (r ≤ θ . ◦ . ◦ . ◦ . ◦ . ◦ Matter 1.44E+46 7.09E+46 1.70E+47 3.14E+47 5.10E+47Field 2.96E+46 1.11E+47 2.39E+47 4.12E+47 6.30E+47 θ . ◦ . ◦ . ◦ . ◦ . ◦ Matter 7.66E+47 1.10E+48 1.52E+48 2.05E+48 2.69E+48Field 8.91E+47 1.19E+48 1.52E+48 1.88E+48 2.26E+48
Note . — The energy is written in units of erg. As for the matter component, the contribution of the rest mass energy is subtracted.
TABLE 4Integrated Energies of Matter and Field (r ≤ θ . ◦ . ◦ . ◦ . ◦ . ◦ Matter 6.89E+43 3.15E+44 8.96E+44 2.03E+45 3.95E+45Field 8.60E+45 3.44E+46 7.73E+46 1.37E+47 2.14E+47 θ . ◦ . ◦ . ◦ . ◦ . ◦ Matter 6.76E+45 1.04E+46 1.47E+46 1.96E+46 2.50E+46Field 3.08E+47 4.19E+47 5.46E+47 6.91E+47 8.52E+47
Note . — Same with Table.3, but integration is done for R ≤≤