太陽フレア・黒点の数値モデリング研究
金子岳史 (新潟大学)
2023年 天文・天体物理夏の学校

太陽・恒星セッション 2023.08.02

自己紹介
金子 岳史 （新潟大学 教育学部 講師）
学位: 博士（理学）（東京大学 2017年）
研究: 数値シミュレーションによる太陽物理学研究
学歴・職歴:
2004 - 明治学園高等学校卒 (福岡県)2007 - 代々木ゼミナール 小倉校 (福岡県)2008 - 東京大学 理科一類
2010 - 東京大学 理学部 地球惑星物理学科
2012 - 東京大学大学院 理学系研究科 地球惑星科学専攻
2017 - 名古屋大学 宇宙地球環境研究所 研究員
2020 - 名古屋大学 宇宙地球環境研究所 特任助教

夏の学校には2012, 2013年度に参加。写真は2012年
度。

2022 - ロッキード・マーティン太陽天体物理学研究所 ポスドク (アメリカ)
2023 - 新潟大学 教育学部 講師

これまでにやったシミュレーション
太陽プロミネンス形成(KY2015)

太陽プロミネンス形成(KY2017)

プロミネンス RTI (KY2018)

浮上磁場/黒点形成(Kaneko+, 2022)
プロミネンス形成・噴出(未出版)

プラズマ噴出/データ駆動(Kaneko+, 2021)

太陽フレアSDO/AIA 131Å

GOES X-ray

太陽表面磁場SDO/HMImagnetogram

拡大

強い磁場が密集太陽黒点
SDO/HMI
continuum

SDO/HMImagnetogram

太陽は、表面の様子を高空間分解能で観測できる唯一の恒星

観測 理論シミュレーション : 例1
観測 (Yohkoh X-ray)
フレアのカスプ構造

磁気流体シミュレーション
Yokoyama et al. (1998)

Tsuneta et al. (1992)

観測：プラズモイド噴出
Ohyama & Shibata (1998) Yohkoh

Glesener et al. (2013) SDO/AIA 131

Masuta et al. (1994)

リコネクションモデル

見える時もある。
カスプ構造は説明。
プラズモイド噴出は見える？

なぜいつも見えない？

観測 理論シミュレーション : 例2
観測：フレア発生時のスペクトル(RHESSI, Grigis & Benz, 2004)

様々な物理過程を含めた磁気流体シミュレーション(Cheung et al. 2018)

多温度プラズマの熱的スペクトルを重ねると
冪スペクトルになる（現在も論争中）

磁気流体シミュレーション

磁気流体方程式の解析解は分からないので、数値的に解く。

磁気流体シミュレーション
観測(SDO/AIA 304Å)

シミュレーション(MHD+熱伝導+輻射冷却)

適当な初期条件と境界条件を与えて磁気流体方程式を数値的に解くと、太陽で観測される現象をある程度は再現できる。もちろん完全ではない。

磁気流体シミュレーション
完全には再現できない理由• 厳密な物理定数が分からない（磁気拡散係数、熱伝導係数など）• 何らかの物理項が足りない（熱伝導, 輻射, 電離非平衡, 非熱的粒子など）• 空間・時間解像度が足りない• 初期条件、境界条件が簡略化されている

観測磁場からコロナの3次元磁場を推定する観測された太陽表面磁場→2次元マップ

観測:太陽表面の磁場(青が負極、赤が正極)

仮定: ローレンツ力=0!×# = 0 ⇔ ∇×# ×# = (3次元磁場構造を推定
フォースフリー磁場
表面からコロナへかけて伸びる3次元磁場(数値モデリング)

Inoue et al. (2018, Nat. Comm.)

データ拘束シミュレーション
RESEARCH LETTERa

Amari et al. (2014)b

Figure 4 | Evolution and eruptionof the twisted flux rope. Selected
field lines of some configurations
obtained after an evolution driven by
photospheric changes, with theconfiguration of December 12,
20:30 UT being taken as the initialcondition. a, In accordance with
Fig. 3b, at around 20:30 UT the
configuration is almost out of
equilibrium and the rope rapidly
rises. b, Sometime later, it erupts
through the arcades that were able toconfine it several hours before.
c, Observed typical shape exhibitedby eruptions on the solar limb (image
taken at a wavelength of 304 Å by
SDO/AIA on 2012 October 14, at
2:32 UT; courtesy of NASA/SDO andthe AIA, EVE and HMI scienceteams). d, Modelled appearance ofthe major eruption of December 13,
projected onto the solar limb,
showing agreement with the typicalshape shown in c (for convenience,
the arcades have been removed).

磁気流体シミュレーションの初期条件に非線形フォースフリー磁場を利用。

より現実的な、複雑な磁場構造から計算をスタートできる。

cd

直接観測できないコロナの3次元磁場の時間発展を推定できる。

課題

初期は観測と一致しているが、時間発展
した後の表面磁場は観測値からズレる。

データ駆動シミュレーション
Pomoell et al. (2019)

下部境界条件に観測された太陽表面磁場の時系列データを利用

時間発展した後も、シミュレーションの下部境界付近の磁場は、観測と一致。

データ駆動シミュレーションの課題
Figure 2. Top view of magnetic ﬁelds for GT, PF, and the four coronal ﬁeld reconstruction models at t/τ0=240. The bottom boundary shows the vertical GT
magnetic ﬁeld (i.e., Bz/B0 at z/H0 = 0), saturating at −0.1 (black) and 0.1 (white). The tubes indicate the ﬁeld lines, where the tubes with reddish (bluish) colors areintegrated from the seeds placed in the positive (negative) polarity. The seeds are identical for all six cases.

異なるデータ駆動手法の比較（Toriumi et al., 2020）

磁場駆動 (Jiang)
正解

電場駆動

(Cheung & DeRosa)

磁場駆動(Yang)

電場駆動(Hayashi)Figure 3. Same as Figure 2 but for the side views. The animated version of this ﬁgure presents the side views from different angles.

これまで、いくつかのデータ駆動手法が提案されている。
手法によって、得られる数値解が異なることが分かっている。
データ駆動手法はまだ発展途上。

データ駆動手法の概要
1. 磁場駆動
観測磁場データ(2次元)を直接、下部境界へインプットする。
2. 電場駆動
観測磁場時系列データから電場を計算し、下部境界へインプットする。
*!

*!"#

観測
時間幅
)

誘導方程式の時間微分の項を、
観測時系列データの差分で置き換え、電場を求める。

!!"# − !!
= −∇×'$#

新たなデータ駆動手法の提案'×!-driven method (Kaneko et al. 2021)
電場(と磁場)から速度場を導出し、下部境界へ与える方法。
1. ゲージ変換により、電場を磁場に直交させる(Fisher et al. 2010)+ = +$ − ∇- → + / # = 0
2. +×#ドリフトの式で、磁場に直交する速度成分を導出
0$ =

+×#
*%

求めた1 $ を下部境界へインプット

ちなみに、この方法は電離圏の研究をしている人の発表を聞いてる時に思いついた（思い出した）。
自分の分野以外の話を聞くのも大事…。

基礎方程式

ptions. The source region consists of two main opposite-polarity magnetic
patches that, in@B= r a⇥ 5(vh⇥interval
B ⌘J ),
ng motion as well as a decrease in the magnetic flux of both polarities
over@t
wo eruptions occurred. Moreover, as shown in panels (d) and (f) of Figure 2, the1 strengthsJ=r ⇥ B,
nts Bx and By are found to increase after the first eruption.
4⇡

where t, ⇢, v, METHOD
B, J , ⌘, and I denote time, mass density, velocity fields, magnet
3. NUMERICAL
ガス圧=0と仮定した磁気流体方程式
異常抵抗モデルand unit vector, respectively. We used theanomalous resistivity in the following f(コロナでは磁気圧>>ガス圧)
he zero-beta MHD equations as follows.

⌘ = 0, (J < Jc ),@⇢+ r · (⇢v) = 0,(1)@t
⌘ = ⌘0 (J/Jc 1)2 , (J Jc ),✓◆2@ (⇢v)
B
BB
11
2
where
Jc +
= 10 I9 G/cm, ⌘0==0,1011 cm2 /s, and we restrict ⌘  ⌘max+r
· ⇢vv(2)= 10 cm /s.
@t8⇡velocity
4⇡ fields which reproduce
We inverted
the observational photospheric mag
この形式で抵抗を与えるとリコネクショ
ンのタイムスケールが現実的な値になるこ@Binverted velocity fields were computed by
the following three steps:= r ⇥ (v ⇥ B ⌘J ),(3)とが、経験的に知られている。
@t
1. Inversion of the induction equation
1J=r ⇥ B,
We solved
an inverse problem of the induction equation as (4)2& や3' の値はテキトーに決める。
4⇡(本当はミクロスケールのプラズマ物理でn+1nd I denote time, mass density, velocity fields, magnetic fields,
currentBdensity,
resistivity,
B
I
obs
obs
決まる。)
=
r
⇥
E
,
y. We used the anomalous resistivity in the following form:⌧

再現を試みたプラズマ噴出イベント
上段: Hα blue wing (-0.5Å)
下段: Hα line center(SMART/SDDI)

再現を試みたプラズマ噴出イベント

再現を試みたプラズマ噴出イベント

再現を試みたプラズマ噴出イベント

シミュレーション結果

観測との比較: 鉛直磁場
2 hr

2 hr

観測
SDO/HMI 鉛直磁場 Bz

シミュレーション

観測との比較: 水平磁場
2 hr

2 hr

観測
SDO/HMI 鉛直磁場 Bx

シミュレーション

観測との比較: 水平磁場
2 hr

2 hr

観測
SDO/HMI 鉛直磁場 By

シミュレーション

エネルギーの時間変化
実線:磁気自由エネルギー

破線:運動エネルギー

磁気自由エネルギー:(熱・運動エネルギーへ変換可能な
磁気エネルギー)$( $ − (%&'%!"## = ',8+.%&' : ポテンシャル磁場
電流が0の磁場。電流が流れていな
いので、ジュール熱も発生しない
し、ローレンツ力も働かない。

エネルギーの時間変化
実線:磁気自由エネルギー

破線:運動エネルギー

赤破線：
観測でプラズマ噴出が発生した時刻
青破線
シミュレーションでプラズマ噴出が発生した時刻シミュレーションは噴出時刻が数十分ずれている

今後の課題

今後の課題
完全には再現できない理由• 厳密な物理定数が分からない（磁気拡散係数、熱伝導係数など）• 何らかの物理項が足りない（熱伝導, 輻射, 電離非平衡, 非熱的粒子など）• 空間・時間解像度が足りない• 初期条件、境界条件が簡略化されている

黒点

フレアを発生させる黒点
bipolarcomplexMaximum X-ray ﬂux

unipolar32

β-type (AR11241)

ﬂare productive

δ-type (AR11492)

Area(Sammis et al., 2000; Toriumi & Wang, 2020)

Simulations of sunspots
Rempel & Cheung (2014)(a) B=10.6 kG

The Astrophysical Journal, 785:90 (17pp), 2014 April 20(b) B=21.2 kG

(a)

(b)

(c)

(d)

Cross section along height (color:|B|)

The Astrophysical Journal, 785:90 (17pp), 2014 April 20

Rempel & Cheung(c),(d) B=21.2 kG but diﬀerent bottom boundary
conditions to the velocity

Figure 12. Comparison of results from four different flux emergence simulations. The top panels compare two simulations with different initial field strength, but thesame overall flux: (a) 10.6 kG; (b) 21.2 kG. In both cases the bottom boundary condition after the flux emergence is an open boundary as described above. In panels(c) and (d) the initial field strength is 21.2 kG, but we used different bottom boundary conditions. In panel (c) we allow for horizontal flows, but no vertical flows inregions with more than 2.5 kG field strength; in panel (d) we continued the inflow in regions with more than 2.5 kG field strength. Except for panel (d), where we did
not observe spot formation, results are comparable.

The bottom boundary was
at -15 Mm.
underlying dynamics differ substantially from the deeper layers

likely not affected, this is not necessarily true during the decay
phase. We find in the presented simulation a rather rapid decay
driven by subsurface flows fragmenting the magnetic field on a33

due to the diminishing pressure scale heights near the surface.
An AR scale flux bundle moves at most a distance comparable toRempel &Inclusion of deep convection zone
RMHD simulation of ﬂux emergence(Toriumi & Hotta, 2019)
Code: R2D2 (RSST and Radiation for Deep
Dynamics, Hotta & Iijima, 2019)•

The bottom boundary was at -200 Mm•

δ-spot was reproduced.

approach

downﬂow

34

Another example using R2D2
Hotta and Iijima (2020)

β-type magnetic distribution

35

What drives the collision of opposite polarity?

7pp), 2019 November 20

36

Vertical velocity at a certain height (red:downﬂow,
blue:upﬂow)
Toriumi & Hotta
Rising flux tube and formation of sunspots

Toriumi &Hotta (2019)

Hotta &
Iijima (2020)

δ-spotβ-spot2525−1/3

velocity scales with ρ0
in the mixing length theory. This setting
approximately maintains the Mach number. Fig. 1 shows the RSST
factor ξ (panel a) and the reduced speed of sound cs /ξ (panel b).
Hotta et al. (2012) show that if the Mach number estimated with the
reduced speed of sound is smaller than 0.7, the RSST cause no side
effect on the result.
We limit the Alfvén velocity to 40 km s−1 to deal with the lowβ region above the photosphere (Rempel, Schüssler & Knölker,
2009b). This does not affect what is occurring at the photosphere.

Figure 2. Initial setup of the magnetic flux tube at z = −35 Mm. The colour
What is the essential diﬀerence?
contour shows the vertical velocity v z , and the contour lines show the axial
magnetic
field Bx . Each contourvelocity
line shows 0, 2.5,ﬁeld?
5, and 7.5 kG.
How can we extract the essence from the complicatedis located about 35 Mm below the photosphere and Btb = 104 G

Downloaded from https://academic.oup.

Figure 1. (Panel a) The RSST factor ξ and (panel b) the reduced speed ofsound cs /ξ are shown.

Aim & Method
Aim
To reveal the impact of convective ﬂows in the convection zone on the formation of
δ-spot magnetic ﬁelds and their energy build-up
MethodParameter surveys (93 cases) by RMHD simulations and a statistical analysis

37

Simulation settings
9 6 Mm

1. Reproduction of solar convection
without magnetic ﬁeld (achieving
statistical equilibrium)
2. Introduce a ﬂux tube
ﬂux transport• Numrcial code: R2D2
(Hotta & Iijima, 2019)• Grids:
horizontal: 768 x 768 (uniform)
depth:
512 (nonuniform)

9 6 Mmﬂux tube
2 0 0 Mm

Basic equations

39

RSST parameter

Solving Radiative Transfer Eq.
(only in z-direction)

Radiative heating
Equa4on of State (OPAL, Rogers et al., 1996)

Background stratiﬁcation(Model S, Christensen-Dalsgaard, 1996)(Hotta & Iijima, 2019)

Supercomputer Fugaku (富岳)Supercomputer Fugaku (@RIKEN, Japan)

40

Typical case of δ-spot formation
Vertical velocity（Red: downﬂow, Blue: upﬂow）

Magnetic ﬁeld strength41

Parameter survey
Parameters of the initial ﬂux tube
Field strength: ﬁxed
Twist: ﬁxed
Horizontal positions: 31 cases(changed uniformly by 3 Mm)• Depth: 3 cases (-22, -26, -30Mm)•••

total 93 cases

42

Vertical velocity at a certain heightResults of parameter survey*( @ ) = 1

Various kinds of magnetic distributions were created only by thediﬀerence in the convective ﬂows surrounding the ﬂux tubes.

43Nonpotential ﬁeld (as a proxy of free energy))/ @# = 1

Nonpotential ﬁeld*)* = # − #*+,

44

Distribution trend of !!" :)01

2 3456

Distribution trend of *)*Temporal meanof *)* in each case
*)*

45

*)*

-

sum all cases& normalized- ./01

Seeking for relationship with convective ﬂows
Mean 1( at diﬀerent heights

*!2 -

vs
correlation?

./01

at the photosphereSeeking for relationship with convective ﬂows
Vertical velocity:≪ 1( ≫Horizontal divergence:≪ ∇3 / 0 ≫47

Vertical vorticity:≪ ∇×0 4 ≫Correlations:)01

2 3456

vs +/ , ∇7 - ., ∇×. /48

Correlation Coeﬃcients (CCs) at diﬀerent depth• Black:*)*- ./01vs1( - ./01→ max |CC| ~ 0.7: good correlation• Red:*)*- ./01vs∇3 / 0 - ./01

→ max |CC| ~ 0.3: weak correlation• Blue:*)*- ./01vs∇×0 ( - ./01

→ max |CC| < 0.1: no correlation

Convective ﬂows (lowpass ﬁlter applied)
Vertical velocity: 1(

Horizontal divergence: ∇3 / 0

Vertical vorticity: ∇×0 (49

Correlations:)01

2 3456vs+/ 2 3456 (lowpass ﬁlter)(a) CCs at 7.5,+66 = 1.44×107% Mm7#
|CC| > 0.8 in z > -50 Mm (max |CC| ~ 0.9)
Very good correlation(b) CCs as a function of < and 7.5,+66
CC is the maximum at 7.5,+66 = 1.44×107% Mm7#
~70 Mm50

Correlations:)01

2 3456vs∇7 - . 2 3456 (lowpass ﬁlter)(a) CCs at 7.5,+66 = 1.44×107% Mm7#
|CC| > 0.8 in z > -30 Mm (max |CC| ~ 0.9)
Very good correlation(b) CCs as a function of < and 7.5,+66
CC is the maximum at 7.5,+66 = 1.44×107% Mm7#
~70 Mm51

Correlations:)01

2 3456vs∇×. / 2 3456 (lowpass ﬁlter)(a) CCs at 7.5,+66 = 1.44×107% Mm7#
|CC| does not exceed 0.4.
(b) CCs as a function of < and 7.5,+66|CC| does not exceed 0.4.
Not good correlation

52

Interpretation

Energybuild-upcoronaphotosphere

RelatedFlux transport:
Gathering to downﬂow region
by converging ﬂow

strong downﬂow
large-scale convection
Within one supergranular cell?

Prediction of high free energy region?
CCs of*)*- ./01

vs 1( , ∇3 / 0 at @ = 0• At @ = 0, convection ﬂows are notaﬀected by magnetic ﬁelds• 1( , ∇3 / 0 at 7.5,+66 = 1.44×107% Mm7#
have still high CC values (~0.7)• The times when the ﬂux emerged in
the photosphere were @~20 − 30 hr• We can constrain high free energy
region from the downﬂow plumeeven before ﬂux emergence(leading time ~ 1 day?)

54

将来研究
境界条件: d001

55

境界条件: d021

R2D2は上空のコロナを含んでいない。
R2D2の結果を下部境界としてデータ駆動手法を適用し、計算領域をコロナへ拡張。
フレアが発生するか検証する。

まとめ• 太陽研究では、観測とシミュレーションの組み合わせが重要- 高時空間分解能の観測データが大量に蓄積されている。
- 計算機の性能上昇により、シミュレーションも高解像度・高精度化。
観測データとシミュレーションの融合（まだ発展途上）• 大規模シミュレーション- スーパーコンピュータにより大規模計算が可能になっているが、出てくる結果も複雑になる。
- 結果をどのように解析し、解釈するかを考える必要。

これからシミュレーションをやってみたい方
CANS+ (シミュレーションコード＋解析コード+解説)

https://www.astro.phys.s.chiba-u.ac.jp/cans/doc/index.html
CANS (シミュレーションコード+解説, ちょっと古い)
hEps://www.astro.phys.s.chibau.ac.jp/netlab/astro/index2.html

これからシミュレーションをやってみたい方
スーパーコンピュータ利用環境
国立天文台 CfCA: Cray XC50 「アテルイII」
修士・博士学生向けカテゴリ → XC-MD (1つの計算あたり400コアまで利用可)研究計画が採択されれば、無料で使える（指導教員の推薦書、成果報告は必要）これからシミュレーションをやってみたい方
流体学校