シミュレーション天文学最前線２００２ （於：三鷹） ２００２．１２．２５

セッション４：磁気流体力学（MHD)−天体MHDシミュレーションの現状と今後ー
柴田一成
京大理花山天文台

天体MHD simulation 論文数の変遷• Astrophysical Journal のみ対象
Simulation および MHD Simulationという用語を Abstract/Keywords に含む論文の数•simulation• 1982-86
273• 1987-91
438• 1992-96 1056
• 1997-01 1451

MHD simulation all papers
236416577444127993923211737

ちなみに K.Shibata の論文数は？
• １９８２−２００１ 共著論文数
MHD simulation in ApJ
39
MHD simulation in other Journal 32(世界の天体MHDシミュレーション論文の１割くらいを、私の周辺で生産）観測の論文
31
理論・流体simulation・ほか 18

Contents(Korea meeting 2001 + α）

• はじめに（基本問題、問題の位置づけ）– 対象となる天体・天体現象– 基本問題の例• 基礎方程式系
– 定式化
– 近似

• 数値計算法– 基本アルゴリズム– 他手法との比較– 計算の難しさと工夫• 最前線の天体シミュレーション• 今後の課題

１．はじめに（基本問題、問題の位置づけ）• 天体ＭＨＤ現象– 惑星磁気圏現象（サブストームなど）– 太陽活動現象（フレア、コロナ、太陽風、コロナ質量放出など）– 恒星活動現象（フレア、コロナ、恒星風、、、）– 降着円盤、宇宙ジェット（ＡＧＮ、連星系、原始星）– 星形成– 星間物質– 銀河団プラズマ太陽コロナ・フレア・ジェット（ようこう軟Ｘ線）

太陽フレ
アから放出される
電磁波
Gamma ray burst

コロナ質量放出（ＣＭＥ）（ＳＯＨＯ／ＬＡＳＣＯ，可視光／人工日食）

かに星雲中の波（ＨＳＴ）(Chandra)astrophysical jets

Protostellar jet (HH 1/2)
AGN jet (Cyg A)

Close binary system (SS433)

宇宙ジェットと降着円盤（原始星ジェット）天体MHD現象共通の問題• 磁場生成ダイナモ （ＭＨＤ乱流）• 構造形成ーループ、磁気圏、星間雲
ＭＨＤ不安定性• プラズマ加熱

磁気リコネクション
ＭＨＤ乱流／shock／non-MHD（非定常ーフレア、定常ーコロナ）• プラズマ加速
磁気遠心力／磁気圧、ガス圧、輻射圧（非定常ーコロナ質量放出／ジェット、定常ー天体風／ジェット）
粒子加速（non-MHD)• 非熱的粒子

天体MHD基礎過程

天体ＭＨＤ現象の共通点（磁場を散逸させるのは、きわめて困難）
磁気拡散時間（電流散逸時間）t D = L / η ≈ 10 L9 T6214

η = ηSpitzer ≈ 10 T64Ａｌｆｖｅｎ時間
磁気レイノルズ数

3/ 2−3 / 2s

cm / s2t flare = 10 − 10 sec
2フレア発生時間

23t A = L / VA = 10 sec
Rm = t D / t A ≈ 10 >> 113MHD近似• 流体近似特徴的長さ＞＞平均自由行程、またはイオンのラーモア半径• ゆっくりした現象（変位電流無視＝非相対論）
特徴的時間＞＞衝突時間、またはイオンのラーモア周期• 準中性
粒子数密度＞＞Goldreich-Julian 密度（n = div (v x B)/e)太陽コロナ・プラズマの特徴的長さ• ラーモア半径−11/ 2mi vc
 B   T ≈ 10cmrLi =  6 eB 100G   10 K • 平均自由行程

22n1  kT  T  8lmfp =  2  ≈ 10 cm 6   9 −3 n e 
 10 K   10 cm 

• フレアのサイズ

rflare ≈ 109 cm−1

フレア・プラズマ（１０＾７Ｋ）における熱伝導の重要性
23nkLn L   T 
≈ 1.4  9 −3  9   7 
tcond =
5/ 2
κ 0T
 10 cm  10 cm   10 K 
2−13kTn  T 4
≈ 4 ×10  9 −3   7 
t rad =nQ(T ) 10 cm   10 K 
1/ 2Ln  B 
 L 
≈ 10  9  9 −3  tA =VA 10 cm  10 cm   10G −5 / 2sec

3/ 2sec−1 / 2sec• フレア・プラズマ（T〜10^７K)では熱伝導の時間スケールは、
Alfven 時間より短い＝＞implicit treatment の必要性

2．基礎方程式
電磁流体（ＭＨＤ）方程式
未知数８個：密度（ρ）、速度ベクトル（ｖ）、磁場ベクトル（Ｂ）、圧力（ｐ）
方程式８個：

MHDの難しさの根源＝電磁流体波（Alfven, fast, slow)
ω − k VA cos θ = 02222

ω − (Cs + VA )k ω + 4Cs VA k 4 cos 2 θ = 0422Group velocity diagram

22v g = ∂ω / ∂k22

MHD wave 特性曲面

２次元： fast ＋−, slow ＋−, 流線＝５本の射影特性曲線
２．５次元／３次元：
fast ＋−, Alfven ＋−, slow ＋−, 流線＝７本の射影特性曲線

Difficulties in Numerical
Astrophysical MHD••

Why difficult ?
because there is gravity(self-gravity, external gravity)
1. Large dynamic range in density, gaspressure, and Alfven speed
2. Wave amplification through verticalpropagation
3. Various instabilities driven by gravitationalenergy– Boundary condition is most difficult

1. Large dynamic range in density,gas pressure, and Alfven speed• For example, in the solar atmosphere, thedensity decreases from the photosphere to thecorona by more than 7 ~ 9 orders of magnitude.
Alfven speed

VA = B / 4πρ

increases rapidly with height.
• Hence we need small grid size (and thus many
grid points in vertical directions) and short time
step due to CFL condition∆t < ∆z /(VA + Cs + V 2 )1/ 222

2. Wave amplification throughvertical propagation• The amplitude of MHD wave rapidly increases
with height when it propagates upward. For
example, the amplitude of slow mode MHD
wave propagating along vertical flux tube
increases asV|| ∝ ρ−1 / 2A−1 / 2• Hence even small amplitude wave at the bottom
quickly become large amplitude wave to formshock and influence upper layers significantly.

Alfven wave model of spicules:
numerical simulation (Kudoh-Shibata 1999)

3. Various instabilities driven bygravitational energy• Various instabilities occur in a
gravitationally stratified gas layer, e.g.,– convective instability, Rayleigh-Taylor
instability, Magnetic buoyancy instability,
Parker instability, etc• These instabilities generate MHD waves
and electric currents, which are the sourceof energy to heat the corona and flares.

Solar emerging flux due toParker instability(Shibata et al. 1989, 1990)

3. Numerical method• Main Newtonian MHD codes in our group– modified Lax-Wendroff code
(originally developed in 1978 by myself, and refinedand extended by Dr. Matsumoto, Dr. Yokoyama)
simple, fast, applicable to many MHD
problems, but often unstable and suffersfrom various numerical errors– CIP-MOCCT code
(developed in 1995 by Dr. Kudoh)• Note: Higher order Godunov scheme is becoming
popular even in MHD code
(Balsara,Ryu,Hanawa,Sano&Inutsuka,,,,)

CIP scheme• A universal solver for hyperbolic equations byCIP(Cubic Interpolated
Pseudoparticle/Propagation)• Developed by
Prof. T. Yabe (Tokyo Institute of Technology)– Yabe and Aoki (1991) Computer Physics
Communications, 66, 219-232• Good at handling contact discontinuity• Can treat all kind of matter (gas, fluid, solid)
simultaneously• Can also be applied to Vlasov equation
See Ogata’s poster

Basic concept of CIP scheme• We want to solve advection equation∂f∂f+u=0∂t∂x• Exact solution is f(x,t)=f(x-ut,0) ifu=constant.

Cubic interpolation• CIP method uses cubic interpolation todetermine the profile between i-1 and i′

Fi ( x) = ai ( x − xi ) + bi ( x − xi ) + fi ( x − xi ) + fi
32nn

Advection of gradient• CIP scheme utilize advection of gradient∂f ′∂f ′+u=0∂t∂xto determine ai and bi

図３

図４

Example 1
Simulation of a log slamming on the water surface. Moving
body is captured with fixed grid system (by Xiao)
150 x 150 x 150

Example 2(simulation of milk crown on
100 x 100 x 35 grids: Yabe et al. )Example 3
Comet Shoemaker-Levy 9 on entry into
Jovian atmosphere (Yabe et al. 1994)

CIP-MOCCT scheme• Developed by Kudoh, Matsumoto, Shibata(1999) Computational Fluid Dynamics
Journal 8, 56-68• For astrophysical MHD problems• Fluid part => CIP scheme• Magnetic field part => MOCCT scheme(Stone and Norman 1992,Evans and Hawley 1988)

MOCCT scheme• CT (Constrained-Transport) scheme(Evans-Hawley 1988)– satisfy divergence free condition (div B =0)• MOC (Method of Characteristics) scheme(Stone-Norman 1992)– stable for Alfven wave

MOCーCT： （１）CTスキーム(Evans & Hawley 1988)

e= internal energy perunit massepsilon = electric fieldCT スキーム : div B=0 を満たす
ように induction eq. を解く

（２） MOC スキーム
(Stone & Norman 1992) :Alfven wave を精度良く解く

ε = − (v B − v B )
*x

*y*
y*x

CIPMOCCT
scheme(Kudoh,Matsumoto,Shibata,1999,
CFDJ 8, 56)

(CIP)

MHD
shock
tube
(Brioand Wu)

Application of CIP-MOC code toMagnetically driven jet from accretion disk
(Kudoh, Matsumoto, Shibata 1998)
Jet

Central object

Accretion disk

Jet velocity ~ Kepler velocity

MHD turbulence driven bymagnetorotational instability(Kudoh, Matsumoto, Shibata 2002)

Summary of Merits of CIP-MOCCT
scheme• To solve contact discontinuity with high accuracy• Applicable to either non-conservative orconservative equations• Applicable to any fluid with complicated physics,
such as gas+fluid+solid, etc.
• Applicable to very low beta plasma(gas pressure << magnetic pressure)• Simple and fast
See Ogata’s poster for astrophysicalapplication

最新研究のレビュー• 磁気リコネクション• 太陽フレア• 宇宙ジェット
• 降着円盤磁気リコネクション

磁気リコネクションの基本問題：
Reconnection rate は何が決めているか？
• Magnetic flux reconnected per unit timeVinηJdΦ / dt
E
MA ====VA B0
VA B0 VA B0 VAt rec = L / V in = t A M A−1

磁気リコネクションの基本問題をめぐる激しい論争•Petschek

spontaneous

driven （境界条件）

Ugai-Tsuda(1977)
Scholer(1989)

Sato-Hayashi(1979)
Priest-Forbes(1986)

YokoyamaShibata(1994)
Sweet-Parker

Biskamp(1986)

SatoHayashi(1979)
Strong inflow is assumed at
the external boundary,which drives reconnection.
＝＞driven reconnection

Tanuma et al. (2001) ApJ

Reconnetion rateFractal Reconnection(Shibata and Tanuma 2001)

太陽フレア2D MHD Simulation of Reconnection withHeat Conduction(Yokoyama and Shibata 1998, 2001)

観測的可視化（理論的に予想されるフレアX線像）

フレア温度の scaling law の発見（Yokoyama and Shibata 1998）∝ B −2T∝B

L

6/7 2/7

フレアの温度は何で決まっているか？
• 磁気リコネクション加熱＝熱伝導冷却のバランスで決まる（Yokoyama and Shibata1998, 2001)

B VA / 4π = κT
2

T∝B

7/2

L/ 2L

6/7 2/7

宇宙ジェット宇宙ジェットの速度は何が決めているか？
• Many MHD simulations– Shibata and Uchida 1986, 1987, 1990, Stoneand Norman 1994, Matsumoto et al 1996,
Hirose et al. 1997, Kudoh et al. 1998,
Ustyugova et al. 1996, …
confirmed thatJet velocity ~Kepler velocity
Why ?

Mass flux and Jet velocity(Kudoh et al. 1998 for thick disk case)
Jetvelocity／Ｖｋ
Mass flux

max2

 VA 2Emg =   ∝ B p
 VK Kato, S. X et al. 2002
(for thin disk case)Summary of Characteristics of non-relativistic
mag. Accretion
centrifugal Disksｖｓ mag. pressure
MHD Jets fromacceleration
Poloidal field

centrifugal
strongmagnetic pressure
weakMag. Field config.

straighthighly twisted

Mass flux

ρC s rρC s r

Jet velocity V
V k 
 C sV k 

Range of application

0.01 < Emg < 122AEmg = mag energy/gra v energy

1/ 32

BpBϕ

VA 

V k 
Cs 

1/ 3Emg < 0.01(Kudoh & Shibata)

Tomisaka (2002) ApJ ５７５，３０６• After adiabatic core isformed, outflow start to beejected• In the case of strong
magnetic field (beta ~ 1),
outflow is accelerated bycentrifugal force plus
magnetic pressure force,but in the weak field case,outflow becomes moreturbulent, and isaccelerated by turbulent
magnetic pressure

Intermittent Accretion and Ejection(K. Sato: poster) (see also Kuwabara)
Mass Ejection Rate

Mass Accretion Ratetime

Magnetic Line of Forceacross Ergospehere
Kerr black holeErgospheret = 7τS

Propagation of Alfven wave:
Electromagnetic energy
transportation

Magnetic field lines

Koide et al. (2002)
Sicence降着円盤Magneto-rotational instability• Balbus-Hawley (1991)• [Chandrasekhar (1961), Velikhov (1959)]• Explains α viscosity of accretion disks(Hawley-Balbus 1991, Hawley et al. 1995,
Brandenburg et al. 1995, Matsumoto-Tajima 1995)

Time variability of 3D magnetized accretion
disk (Kawaguchi et al. 1999, Machida et al. 2001
Hawley et al. 2001, Stone et al. 2001)

Sano and Inutsuka (2001) ApJ 561, L179•

Local simulation ofaccretion disk
• Saturation ofmagnetorotational
Instability=> Reconnection

今後の課題• 物理・天体物理– 磁気リコネクションの基本問題
3D tearing instability の非線形発展＝＞ フラクタル？ 乱流？
– MHD＋additional physics(e.g., heat conduction, cosmic ray,,,)– ねじれた磁束管の浮上にともなう３次元リコネクション(太陽フレアのモデル）– ３次元MHDジェット（安定性、コリメーション、リコネクション）– 降着円盤（MRI saturation mechanism,ダイナモ）

今後の課題（続）• 数値MHD– CIP-MOCCT code をあらゆるMHD問題に適用してみる （e.g., 一般相対論的MHDコード）– 多階層統合コードーAdaptive mesh refinement– MHD＋プラズマ粒子／Vlasov統合コード