プラズマ粒子シミュレーションにおけるさまざまな速度分布関数の乱数生成方法

Loading various velocity distributions in particle simulation
銭谷誠司 (Seiji Zenitani),1, 2 中野慎也 (Shin’ya Nakano),3, 4 and 岡光夫 (MitsuoOka)5, 61)

オーストリア科学アカデミー 宇宙科学研究所

2)Space Research Institute, Austrian Academy of Sciences, Schmiedlstraße 6,8042 Graz, AUSTRIAa)
3)〒 190-8562 東京都立川市緑町 10-3 統計数理研究所

4)

The Institute of Statistical Mathematics, 10-3 Midori-cho, Tachikawa,Tokyo 190-8562, JAPAN
5)

カリフォルニア大学バークレー校 宇宙科学研究所

6)Space Sciences Laboratory, University of California, Berkeley, California 94720,

USA(Dated: 12 March 2024)

宇宙空間物理学では、さまざまな形状の速度分布関数を取り扱う。本稿ではその中で、カッパ分布・フラットトップ分布・ロスコーン分布（subtracted Maxwellian）・
相対論的マクスウェル分布・相対論的カッパ分布の５つの分布の性質とプラズマ粒子
シミュレーションでの乱数生成方法を紹介する。

※本稿を頻繁に更新しています。随時、以下の URL で最新版を確認してください。

https://sci.nao.ac.jp/MEMBER/zenitani/research-j.html#PIC
Keywords: カッパ分布, フラットトップ分布, ロスコーン分布, 相対論的マクスウェ
ル分布, 相対論的カッパ分布, プラズマ粒子シミュレーション, 乱数

a)Electronic mail: seiji.zenitani@oeaw.ac.at1

CONTENTS
I. カッパ分布2

II. フラットトップ分布

10III. ロスコーン分布

16IV. 相対論的マクスウェル分布

20

A. Swisdak 法

21

B. 修正 Canﬁeld 法

25

C. 比較

28

V. 相対論的カッパ分布29

VI. 謝辞32

References

I.

32

カッパ分布
宇宙空間物理学では、高エネルギー部分に「べき乗」成分をもつカッパ分布（kappa dis-tribution）がよく利用される。カッパ分布は 1960 年代に Vasyliunus, Olbert らのグループによって宇宙空間のプラズマ研究に導入された。18,27 その後、太陽圏のプラズマ観測が進むにつれて、Maxwell 分布に収まらない高エネルギー成分が普遍的に観測されるようになり、
熱的コア成分と高エネルギーの「べき乗」成分をまとめて扱うことができるカッパ分布の有用性も認知されてきた。現在では、太陽・太陽風・太陽圏深部などさまざまな環境のプラズ
マ分布関数の議論でカッパ分布が使われている。13カッパ分布の表式は以下の通りである。

NΓ(κ + 1)fκ (v)d v =2
3/2(πκθ ) Γ(κ − 1/2)
3



v2
1+ 2
κθ

−(κ+1)d3 v(1)

ここで N は数密度、κ はカッパ分布のべき指数、θ は最頻速度、Γ はガンマ関数である。特に、微分エネルギーフラックス（F(E)dE = N (E)vdE = fκ (v)v 2 vd3 v ）が高エネルギー側でF(E) ∝ E −κ となる。また、κ → ∞ の極限では、Maxwell–Boltzmann 分布に漸近する。
2

最初に規格化定数を確認する。規格化定数を A とおいて、v → v に変換する。
−(κ+1)−(κ+1)
ZZZ 
Z ∞
v2
v2
3
A
1+ 2
d v = 4πA
1+ 2v 2 dv = Nκθκθ0
積分式の中で x ≡ v 2 /κθ2 とおいて、ベータ関数 B(x, y) の公式を使うと


Z ∞
1
x1/2
1
3
1
π 1/2
Γ(κ − 1/2)
2 3/2
2 3/2
(κθ )
dx
=
(κθ
)
B
,
κ
−
=
(κθ2 )3/2
κ+12(1 + x)22
2
4
Γ(κ + 1)
0

(2)

(3)

となり、これを (2) 式に代入・整理して以下を得る。

Γ(κ + 1)N2
3/2(πκθ ) Γ(κ − 1/2)

(4)

Fκ (v) ≡ fκ (v)4πv 2(5)A=
等方カッパ分布の radial 密度分布はとなるので、(∂/∂v)Fκ (v) = 0 より最頻値が v = θ であることがわかる。
最後にエネルギー密度を計算する。同様にベータ関数の公式を利用して
−(κ+1) −(κ+1)
Z ∞
Z ∞
v2
v2
1 2 2
E = 4πA
1+ 2
mv v dv = 2πmA
1+ 2
v 4 dv
κθ
2
κθ
0
0


3
κ
5
3
,κ −
N mθ2
=
= πmA(κθ2 )3/2 B
22
4 κ − 3/2

(6)(7)

となる。途中の式展開の都合上、κ > 3/2 でなければならない。自由度 3、比熱比 γ = 5/3のプラズマ系のエネルギーは、温度（ボルツマン定数を含む）を T としてE=

1
3
NT = NT
γ−1
2

(8)

κ
mθ2
2κ − 3

(9)

であるから

T =
を得る。

カッパ分布の例を図 1 に示す。フォーマットやべき指数 κ の値は、当該分野で有名な

Summers & Thorne 24 論文の図１に合わせたものである。κ の値が大きくなるにつれて、分布が Maxwell 分布に近づいていくようすがわかる。
図 2 では、κ = 3.5, ∞ のカッパ分布（黒線）に加えて、κ = 3.5 分布に内接する２種類の Maxwell 分布を青線で示した。密度 NM 、熱速度 vM の Maxwell 分布で


 32

v2
1
3
exp − 2 d3 v
fM (v)d v = NM
2
πvM
vM

3

(10)

Normalized Distributions

101
10 1
10 3
10 5
10 7
10 9

= 2.0
= 6.0
= 25.0Maxwellian ( = )

4
2 0
2
4
Normalized Velocity (v/ )FIG. 1. カッパ分布の関数形のサンプル。カッパ分布は κ = ∞ の極限で Maxwell 分布（点線）に収Normalized Distributions

束する。

10 1
10 2
10 3
10 4

= 3.5
Maxwellian ( = )Inscribed curve (v0 = )Inscribed curve (v0 = 0)

2
1 0
1
2
Normalized Velocity (v/ )FIG. 2. カッパ分布（κ = 3.5, ∞; 黒）および内接 Maxwell 分布（v0 = θ, 0; 青）。

4以下のようにパラメーターをとると

vM = θ,−(κ+1)
Γ(κ + 1) − 3
1
NM ≡ eκ 2 1+NΓ(κ − 1/2)κ(11)

最頻速度 v0 = θ でカッパ分布に内接する Maxwell 分布を得ることができる。15 e はオイ
ラー数 e = 2.718 . . . である。また
2
vM
=

κ 2
θ ,
κ+1


NM =

1
κ+1

3/2

Γ(κ + 1)NΓ(κ − 1/2)(12)とすると、v0 = 0 でカッパ分布に内接する分布が得られる。これらの内接分布の導出方法
は Zenitani & Nakano 29 の Appendix を参照されたい。

1.0

Nnt/N
nt/

Fraction

0.8
0.6
0.4
0.2
0.0 2

4

6

8

10

FIG. 3. べき指数 κ のカッパ分布における高エネルギー成分の数密度（Nnt /Nκ ）およびエネルギー（Ent /Eκ ）の比率。曲線は理論値（(14) 式および (15) 式）で、黒点は数値検証結果。

Oka et al. 15 は、カッパ分布から内接 Maxwell 分布を引いた高エネルギー側の成分
Z ∞
Nnt =


Fκ (v) − FM (v) dv

(13)

v0

をカッパ分布固有の非熱的高エネルギー成分と見ることを提案している。ここでは、内接点として v0 = θ を採用する。図 3 は、この高エネルギー成分が数密度（Nnt /Nκ ）およびエネ5ルギー（Ent /Eκ ）で全体に占める割合を、 κ の関数としてプロットしたものである。15 後日、導出を追記する予定であるが、密度比は
R∞
R∞


Fκ (v)dv − θ FM (v)dv
Nnt
1 3
e κκ−1/2 Γ(κ + 1) Γ̄(3/2, 1)
θκR∞== I 1+κ
κ− ,
−
Nκ
2 2
(κ + 1)κ+1 Γ(κ − 1/2)Γ(3/2)
Fκ (v)dv
0(14)であり、エネルギー比は

Ent
=
Eκ

R∞
θ

κ
= I 1+κFκ (v)



R∞
1
2
mv
dv
−
FM (v)2R
θ∞
1Fκ (v) 2 mv 2 dv
0

κ−3/2

1 5
κ− ,
2 2

−

1
mv 2
2



dv

eκ
Γ(κ + 1) Γ̄(5/2, 1)
κ+1(κ + 1) Γ(κ − 3/2)Γ(5/2)(15)となる。ここで、Γ̄(α, x) は第二種不完全ガンマ関数、Ix (α, β) は正規化された不完全ベータ関数である。図からは、べき指数 κ が小さいほど高エネルギー成分の比率が大きくなることが明らかである。エネルギー面の寄与はさらに大きく、κ = 10 でおよそ２割、κ ≲ 4 で半分以上、κ → 3/2 でほぼ全てを高エネルギー成分が担っている。図 2 からも想像できるが、内接点として v0 = 0 を採用すると、密度比は

3/2
NM
1
Γ(κ + 1)
Nnt
=1−
=1−
Nκ
Nκ
κ+1
Γ(κ − 1/2)

(16)エネルギー比はEnt
NM T M=1−=1−
Eκ
Nκ T κ



1
κ+1

5/2

Γ(κ + 1)Γ(κ − 3/2)

(17)となる。このとき、最頻速度 v = θ 付近でも多くの粒子がカウントされるため、高エネル
ギー成分の数密度比・エネルギー比はさらに大きくなる。16 従って、内接点 v0 は問題の性質にあわせて注意深く選ぶ必要があるだろう。
プラズマ粒子（Particle-in-cell; PIC）シミュレーションやモンテカルロシミュレーション
でカッパ分布プラズマを扱うためには、カッパ分布に従う速度分布を乱数生成する必要がある。カッパ分布は、統計学の基本分布の１つである t 分布（「スチューデントの t 分布」ともいう）の多変量版に対応する。変数が１つの単変量の t 分布を乱数生成する方法は、例えばBailey 4 法などが知られている。そして、多変量 t 分布を乱数生成する方法も統計学のソフ
トウェアパッケージで使われている。10 本稿では、３変量の t 分布に対応させてカッパ分布
を乱数生成する方法 1 を少しアレンジして紹介する。
準備のため、Maxwell–Boltzmann 分布の生成方法を考える。
 m  23
mv 2
3
d3 v
exp −
f (v)d v = N
2πT
2T

6

(18)

ここで、熱速度を

r
vth ≡

2T
m

(19)

と定義する。Maxwell–Boltzmann 分布に従う乱数を生成する標準的な方法は Box–Muller 法であろう。5 Box–Muller 法は２つの独立した一様乱数 u1 , w1 ∼ U (0, 1) を使って
pr1 = −2 ln u1 sin(2πw1 )pr2 = −2 ln u1 cos(2πw1 )

(20)
(21)

のように２つの正規乱数 r1 ∼ N (0, 1), r2 ∼ N (0, 1) を生成する。このとき、r1 , r2 はそれ
2
ぞれ独立した分布である。Maxwell–Boltzmann 分布は分散 σ 2 = (1/2)vthの正規分布に対応

するから、(20), (21) 式を σ 倍すれば速度分布を得ることができる。実用上は、これらの手
順を２回繰り返して、４つの一様乱数（u1 , u2 , w1 , w2 ）から速度ベクトル３成分（r1 , r2 , r3 ）と余りの１つ（r4 ）を生成することが多いだろう。また、多くの言語では、一様乱数 u の定義域は 0 ≤ u < 1 である。粒子数が多い場合、u = 0 が現れて ln u = − inf となるため、例
えば u = 1 − u とおいて u = 0 を避ける工夫も必要である。
本稿では、べき指数 κ が 2 以上の整数あるいは半整数のカッパ分布を考える。ν = 2κ−1、

σ 2 = κθ2 /ν とおくと、(1) 式を多変量（p 変量）の t 分布に対応させることができる。カッ
パ分布の場合は p = 3 である。

NΓ[(ν + p)/2]
f (v)d3 v = p/2 p/2 pν π σΓ(ν/2)



1 v2
1+
ν σ2

− (ν+p)
2

d3 v

(22)

一般に、p 次元の多変量 t 分布に従う乱数は、p 次元の正規分布（(r1 , r2 , · · · rp ) ∼ Np (0, 1)）と自由度 ν のカイ二乗分布 r5 ∼ χ2ν の組み合わせで生成することができる。10,12 そこで、分
散 σ 2 = 1 の３変量 t 分布に従う速度ベクトル v̄ の３成分を

r1v̄x = p,r5 /νr2v̄y = p,r5 /νr3v̄z = p
r5 /ν

(23)

で生成する。自由度 ν のカイ二乗分布 r5 は、独立した正規乱数 ni ∼ N (0, 1) を ν 個足し合わせたものである。

r5 =

ν
X

n2i

(24)

i=1

ここで Box–Muller 法の表式（(20) 式・(21) 式）を利用して (24) 式をさらに簡略化する。

ν が偶数で ν = 2k となる場合は、Ui 、Wi を一様乱数として
k 
k

X
Y

r5 =
− 2 ln Ui sin2 (2πWi ) + cos2 (2πWi ) = −2 ln
Ui
i=1

i=1

7

となるから、Wi は不要で Ui のみが必要である。この一様乱数 Ui についても Ui = 0 を避ける工夫が必要である。また、ν = 2k + 1 が奇数の場合は

r5 = −2 lnkY

Ui + n̄2

(25)

i=1

となり、もう１つ正規乱数 n̄ ∼ N (0, 1) が必要になる。
最後に、分散 σ 2 を考慮して目的の v ベクトルを得る。

vx = σv̄x ,vy = σv̄y ,vz = σv̄z(26)

温度が非等方であったり、重心速度がゼロでない場合のカッパ分布の一般形は以下のようになる。磁場方向と垂直方向で温度が異なるバイカッパ分布（bi-kappa distribution）もここに含まれる。

Γ[(ν + p)/2]
f (v)d v = p/2 p/2 1/2
ν π |Σ|Γ(ν/2)
3N

1
1 + (v − µ)T Σ−1 (v − µ)ν

− (ν+p)
2

d3 v

(27)

ここでも p = 3 で、µ は平均速度ベクトル、Σ は分散共分散行列、|Σ| はその行列式である。
そして、コレスキー分解 Σ = AAT で得られる行列 A を使って、(26) 式の代わりに

v = µ + Av̄

(28)とする。この行列 A を Maxwell 分布の段階でかけて、事前に非等方 Maxwell 分布を作ってから、カイ二乗分布を使ってカッパ分布を生成しても良い。
この手順を表 V の擬似コードにまとめた。一様乱数のいくつか（ui , Ui ）については、定義域が 0 を含まないよう注意する。最初に Box–Muller 5 の手続きを繰り返して、３成分の速度ベクトルを得る。次に、κ が半整数（5/2, 7/2, · · · ）の場合は κ − 1/2 個の一様乱数を使って自由度 ν のカイ二乗分布 χ2ν を計算する。κ が整数の場合は、κ − 1 個の一様乱数と１個の正規乱数を使って自由度 ν のカイ二乗分布 χ2ν を計算する。１個の正規乱数には、Maxwell–Boltzmann 分布を作る際に余った r4 を流用すれば良い。最後に、カイ二乗分布と分散を考慮して速度３成分を得る。一連の手続きで κ + 3 あるいは κ + 7/2 個の一様乱数を使っている. (26) 式の代わりに (28) 式を使う場合は、適宜、定数 A を行列 A に書き直せば良い。
このアルゴリズムを使って κ = 3.5 の等方カッパ分布を乱数生成した結果を図 4、図 5 に示す。速度の関数 F (|v|)、あるいはエネルギーの関数 F (E) の形で、数値結果（青色のヒストグラム）と解析解（黒の実線）を比較している。また、点線は前述の内接 Maxwell 分布8Normalized density

100

v 7.0

10 2
10 4
6
10 10
1

100
|v|/

101

FIG. 4. κ = 3.5 のカッパ分布に従って乱数生成した粒子速度分布（ヒストグラム）と解析解（黒のNormalized density線）。黒の点線は内接 Maxwell 分布である。

100
10 2

E 4.0

10 4
6
10 10
1

100

E/m 2

101

102

FIG. 5. κ = 3.5 のカッパ分布に従って乱数生成した粒子エネルギー分布（ヒストグラム）と解析解（黒の線）。黒の点線は内接 Maxwell 分布。

9

TABLE I. カッパ分布の乱数生成法
ν ← 2κ − 1generate uniform random u1 , u2 , w1 , w2 ∈ (0, 1]√v̄x ← −2 ln u1 sin(2πw1 )√v̄y ← −2 ln u1 cos(2πw1 )√v̄z ← −2 ln u2 sin(2πw2 )√n ← −2 ln u2 cos(2πw2 )if κ is an integer thengenerate uniform random U1 , U2 , · · · , Uκ−1 ∈ (0, 1]
Qκ−1 
2
χ2ν ← −2 ln
i=1 Ui + nelseif κ is a half integer thengenerate uniform random U1 , U2 , · · · , Uκ−1/2 ∈ (0, 1]Qκ−1/2 Uiχ2ν ← −2 ln
i=1
endif
s
s
p
κθ2 /νκθ2(2κ − 3)T
A← p
=
=
2
2
χ
mχ2ν
χν /ν
ν
vx ← µx + Av̄x
vy ← µy + Av̄yvz ← µz + Av̄z
return vx , vy , vz

（v0 = θ）である。両図からは数値結果と解析解がよく一致することがわかる。解析解はそれぞれ ∼ v −2κ , ∼ E −κ−1/2 に漸近する。図 4 では、カッパ分布の最頻速度が |v| = θ である
ことを確認できる。

II.

フラットトップ分布
フラットトップ分布は宇宙空間で観測される、低エネルギー側でフラットな一定値をとる速度分布関数である。2,8,14,17,26 高エネルギー側は（急峻な）べき乗型になることが多い。
この独特の性質を表現するためによく用いられるのが以下の関数形である。規格化定数のガ10ンマ関数の使い方によって、多少、数式表現の細部が異なることがある。


1

Γ 1+
 v 2κ −(κ+1)/κ
3Nft
κ
3

 
 1+
fft (v)d v =
d3 v
3
1
4πθ3
θ
Γ 1+
Γ 1−
2κ
2κ

(29)

θ はフラットトップ分布の特徴速度、そして κ は高エネルギー側のべき指数をコントロール
する係数である。分布関数の形状は図 6 を参照していただきたい。異方性を考慮して、以下のようなバイ・フラットトップ分布を定義することもできる。

1!κ !−(κ+1)/κ
Γ 1+
2
2
v
3N
v
κ
∥
ft
 
 1+

fbft (v)d3 v =
+ ⊥
d3 v
2
2
3
1
4π(θ∥ θ⊥)θ∥2 θ⊥
Γ 1+
Γ 1−
2κ
2κ
等方分布に戻って極座標を考えると


1

Γ 1+
 v 2κ −(κ+1)/κ
3Nft
κ
2
 
 1+
Fft (v) ≡ fft (v)4πv = 3 
v2
3
1
θ
θ
Γ 1+
Γ 1−
2κ
2κ

(30)



(31)

となる。
べき指数 κ は、κ > 3/2 を満たす必要がある。カッパ分布でも κ > 3/2 だったことを思い出そう。(1) 式と (29) 式は、どちらも高エネルギーのべき乗成分 f (v) ∼ v −2(κ+1) を持っR ∞ −2κ+2ている。この成分のエネルギー ∝v
dv が κ = 3/2 で無限大に発散するため、双方とも κ > 3/2 でなければならないわけである。逆の極限 κ → ∞ では、(29) 式右辺の右半部はステップ関数のように振る舞う。ガンマ関数部も Γ(1) = 1 に漸近するから分布関数はfft (v) →

3Nft
4πθ3

fft (v) → 0(|v| < θ)(32)(θ < |v|)(33)となり、プラズマが半径 θ の球内に集中することがわかる。
次に、一般化ベータプライム分布を導入する。一般化ベータプライム分布はベータ関数

B(α, β) を用いて
 αp−1  p −(α+β)
x
x
1+
q
q

 p −(α+β)αp Γ(α + β)
x
= αp
1+
xαp−1 .
q Γ(1 + α)Γ(β)
qpB (x; α, β, p, q) =qB(α, β)′Z ∞B′ (x) dx = 1

(34)
(35)

0

11

100
10 1
10 2
10 3
10 4 0

Flattop dist.
Sampling dist.
f1(v)
f2(v)

1

v/

2

3

FIG. 6. フラットトップ分布（κ = 3, θ = 1）の位相空間密度の一次元カット（太線）。サンプリング
分布（細線）と２つの基線 f1 (v), f2 (v) については本文で記述。

と表される。これを用いると、フラットトップ分布のプラズマ密度および粒子エネルギー密
度はZ ∞Fft (v)dv

Z ∞ 
3
1
′
= Nft
B v; , 1 − , 2κ, θ dv
2κ
2κ
0

N=

0

= Nft


Z ∞
1 2
E=
Fft (v)mv dv
20 Z ∞ 
 
5
3
3
1
5
2 B 2κ , 1 − 2κ′
= N mθ
B v; , 1 − , 2κ, θ dv
3
1
2
2κ
2κ
B 2κ
, 1 − 2κ01
Γ(5/(2κ))Γ(1 − 3/(2κ))= Nft mθ22
Γ(3/(2κ))Γ(1 − 1/(2κ))

1 2 3 Γ(1 + 5/(2κ))Γ(1 − 3/(2κ))
= Nft
mθ
,
2
5 Γ(1 + 3/(2κ))Γ(1 − 1/(2κ))

(36)

(37)

と表すことができる。途中、Γ(1 + α) = αΓ(α) の関係を使うと便利である。こうして検証す
ると、(29) 式の係数は Oka et al. 17 のものと一致するが、最初に関数形を導入した Thomsen
31
et al. 26 とは一致しないことがわかる。(36) 式の 2κ
のところに 2κを代入すれば Thomsen et12al. 26 の係数が得られるので、おそらく計算ミスだろう。また、右辺の積分の中のベータプ
ライム分布の引数は正でなければならないが、既に述べたように κ > 3/2 なので問題ない。

1.0

Nnt/N
nt/

Fraction

0.8
0.6
0.4
0.2
0.0 2

4

6

8

10

FIG. 7. べき指数 κ のフラットトップ分布における高エネルギー成分の数密度（Ntail /N ）およびエネルギー（Etail /E ）の比率。曲線は理論曲線（(38) 式および (39) 式）で、黒点は数値検証結果。

続けて、フラットトップ分布の内部の熱的成分・非熱的成分を評価したい。ここでは Oka

et al. 17 に倣って、E < 21 mθ2 以下を熱的エネルギー、E > 21 mθ2 以上を非熱的エネルギーとする。前述のようにフラットトップ分布はベータプライム分布で表すことができ、ベータ
プライム分布の累積密度関数は正規化された不完全ベータ関数 Ix (α, β) = Bx (α, β)/B(α, β)x (α, β) と書くことができる。今回は v = θ を境に熱的・非熱的成分を区切るを用いて、I 1+xので、この関数に x = 1 を代入したものが熱的成分の比率を表す。非熱的成分の比率は、密
度では

R∞





Fft (v)dv
13
1 3= 1 − I1
,1 −= I1 1 − ,
2
2
N
2κ
2κ
2κ 2κ
R∞




1
2
F
(v)(
mv
)dv
3
5
5
3
ft
2
θ
= 1 − I1
,1 −= I1 1 − ,
2
2
E
2κ
2κ
2κ 2κ
θ

(38)
(39)

となる。途中、Ix (α, β) = 1 − I1−x (β, α) の関係を用いた。これらの κ に対する依存性を図7（Oka et al. 17 の Fig. 1(d) に対応）に示す。カッパ分布（図 3）と同様の傾向を確認する
ことができる。

13

次に、粒子シミュレーションで、フラットトップ分布に従う速度分布関数を乱数生成することを考える。少なくとも２通りの方法があり、そのひとつは、フラットトップ分布がベー
タプライム分布であることを利用してベータプライム乱数を生成する方法である。31 この方法は別のガンマ乱数生成ルーチンを必要とするため、本稿では単体で動作する方法を紹介する。
図 6 の実線（細線）で示すように、フラットな分布 f1 (v) と power-law 分布 f2 (v) を組
み合わせて、「サンプリング分布」を用意する。θ ̸= 1 の場合は後でスケールし直せば良いので、θ = 1 とする。A は計算を簡単にするための定数である。

1
A(|v| < 1)4π
1f2 (v) =A v −2(κ+1) (1 ≤ |v|)
4π


1
Γ 1+
3Nft
κ
 

A≡ 3 
3
1
θ
Γ 1+
Γ 1−
2κ
2κ
f1 (v) =

(40)
(41)

(42)

極座標に直すとF1 (v) = 4πv 2 f1 (v) = Av 2

(v < 1)(43)

F2 (v) = 4πv 2 f2 (v) = Av −2κ

(1 ≤ v)(44)である。これに対応して、0 ≤ r1 (v) ≤ 1 あるいは 0 ≤ r2 (v) ≤ 1 を満たす棄却関数を準備しておく。

−(κ+1)/κ
Fft (v)= 1 + v 2κF1 (v)
−(κ+1)/κ
Fft (v)
1r2 (v) ==+1F2 (v)v 2κ
r1 (v) =(v < 1)(45)(1 ≤ v)(46)(31) 式の v の生成手順に移る。まず、サンプリング関数 F1 (x), F2 (x) が作る面積Z 1
1
F1 (v)dv = A
(47)
3
0
∞

Z ∞
1
1
1−2κ
=
F2 (v)dv = Av
A(48)
1 − 2κ
2κ − 1
1
1の比 (2κ − 1) : 3 に応じて、一様乱数 X1 を用いて F1 (v) あるいは F2 (v) に処理を割り振る。F1 (v), F2 (v) は v のべき乗分布であるから、これに従う乱数を容易に生成することができる。次に、もう１つの一様乱数 X2 と棄却関数を用いて、X2 < r1 (v) あるいは X2 < r2 (v)14TABLE II. フラットトップ分布の乱数生成法
2κ − 13
p1 ←, p2 ←
2κ + 22κ + 2
repeatgenerate X1 , X2 ∼ U (0, 1)
if X1 ≤ p1 thenx ← (X1 /p1 )1/3
if X2 < (1 + x2κ )−(κ+1)/κ , breakelsex ← ((1 − X1 )/p2 )1/(1−2κ)if X2 < (x−2κ + 1)−(κ+1)/κ , break
endif
end repeat
generate X3 , X4 ∼ U (0, 1)pv⊥1 ← 2θ⊥ x X3 (1 − X3 ) cos(2πX4 )pv⊥2 ← 2θ⊥ x X3 (1 − X3 ) sin(2πX4 )v∥ ← θ∥ x (2X3 − 1)
return v⊥1 , v⊥2 , v∥

を満たしていれば、その粒子を採用して次のステップに進む。満たしていなければ、その粒子を捨てて処理をやり直す。最後にこれをカルテシアンの３成分に戻せば良い。これらのアルゴリズムをまとめたものが 表 II である。最初の条件分岐に使った一様乱数 X1 を次の処理にも流用して、乱数の数を節約している。アルゴリズムの最後の部分で v⊥1 , v⊥2 , v∥ 成分に θ⊥ , θ⊥ , θ∥ を掛けて、(30) 式のバイ・フラットトップ分布を生成する。
棄却法の採択効率は



 

3
1
R∞
Γ 1+
Γ 2−
Fft (v)dv
2κ
2κ
0


=
eﬀ ft (κ) = R 1
R∞1
F1 (v)dv + 1 F2 (v)dv
0
Γ 2+
κ

(49)

である。これを図 8 に示す。ちょうど eﬀ ft (3/2) = 3/5 であり、κ > 3/2 での採択効率は> 60% である。
表 II のアルゴリズムで生成したフラットトップ分布のエネルギーヒストグラムを図 9 に示す。各パラメーターは κ = 3, θ∥ = θ⊥ = θ = 2、粒子数は 106 個である。ヒストグラムの

15

1.0
0.8
0.6
0.4
0.2
0.00

Acceptance efficiency

1

2

3
4
index

5

6

FIG. 8. フラットトップ分布生成法の棄却法の採択効率。曲線は理論曲線（(49) 式）で、黒点は数値検証結果。

数値解と、実線の解析解が一致していることがわかる。低エネルギー側のビンが荒いためわ
かりにくいが、エネルギー空間の分布関数は低エネルギーのフラットトップ側では ∝ E 1/2となり、高エネルギー側では ∝ E −(κ+1/2) となる。

III.

ロスコーン分布

内部磁気圏のプラズマ速度分布関数は、磁力線に沿った方向の成分が消えた「ロスコーン」を備えていることが多い。このようなロスコーン分布を扱うためにさまざまな分布関数モデルが使われているが、その中で最もよく使われているのが Ashour-Abdalla & Kennel 3が提案した subtracted Maxwell 分布である。関数形は以下の通り。
!

 2

 2


2
v∥2
N0
1
v⊥
v⊥
1−∆
v⊥
f (v∥ , v ⊥ ) = 1/2 exp − 2 × 2 ∆ exp − 2 +exp − 2 − exp − 2
π θ∥
θ∥
πθ⊥
θ⊥
1−β
θ⊥
βθ⊥

(50)
ここで、β ∈ [0, 1] と ∆ ∈ [0, 1] はロスコーンの形状をコントロールするパラメーターである。β = 0 のときは分布は bi-Maxwellian になり、β = 1 に近づくにつれてロスコーンが顕著になる。また、∆ = 1 のときもロスコーンが埋まって bi-Maxwellian になり、∆ = 0 の

16

Normalized density

101

E 0.5

10 1

E 3.5

10 3
5
10 10
2

10 1

100
E/m 2

101

102

FIG. 9. フラットトップ分布のエネルギーヒストグラムと解析解（黒線）。

ときにロスコーンが顕著になる。本稿ではまず、∆ = 0 の場合を議論する。θ∥ , θ⊥ は磁力線√
水平／垂直方向の熱速度である。熱速度の定義が 2 倍違う場合は適宜補正する。v∥ 成分と v⊥ 成分は独立で、v∥ 成分は Maxwell 分布である。v⊥ 成分では、Maxwell 分布から別のMaxwell 分布を引いてロスコーンの形を模倣している。∆ = 0 のときの v⊥ 成分を円柱座標系で考えると2πv⊥
gV (v⊥ ) = 2πθ⊥ (1 − β)


 2


2
v⊥
v⊥
exp − 2 − exp − 2
θ⊥
βθ⊥

(51)

2
2
となる。x ≡ v⊥
/θ⊥とおいて変数変換すると、以下の表式を得る。

1gX (x) =
1−β





x
exp (−x) − exp −
β


(52)

ここで subtracted Maxwellian の圧力を見てみよう。v∥ 方向の圧力 P∥ は Maxwell 分布のものと同じだから、P∥ = N
17

1 2
mθ
2 ∥


(53)

となる。v⊥ 方向の圧力 P⊥ はZZZno NZ ∞
n
o
1
3
2
2
2
P⊥ =
d v f (v∥ , v ⊥ )m(v⊥1 + v⊥2 ) =dv⊥ gV (v⊥ )mv⊥
22 0
 −x/β 
xe
(
)
−x
2


Z ∞ nZ ∞xe − βo
1 2
1 2
β2
dx xgX (x) = N
dx
=N
mθ
mθ
2 ⊥
2 ⊥
1−β
0
0


1 2
mθ (1 + β)
(54)
=N
2 ⊥
となる。β = 0 のときは Maxwell 分布と圧力と一致する。そして、β が大きくなると、v⊥ の小さな粒子が減って v⊥ の大きな粒子が増えるため、垂直方向の実効温度が増大する。圧力の非等方性をパラメーター表示すると

A=

P⊥
θ2
− 1 = ⊥2 (1 + β) − 1P∥θ∥となる。また、エネルギー密度は以下のようになる。
ZZZ


1
1
2
2
d3 v f (v∥ , v ⊥ ) · m(v∥2 + v⊥
) = N m θ∥2 + 2θ⊥(1 + β)
2
4

(55)

(56)

次に、subtracted Maxwellian を乱数生成する方法を議論する。v∥ 成分は既存の方法で生成できるから、我々は v⊥ 成分のみを考えればよい。以下のように、２つの一様乱数 U1 , U2を使って新しい変数 x を作ってみよう。30x ← − log U1 − β log U2(57)s ≡ − log U1(58)t ≡ −β log U2(59)s ∼ Gs (s) = exp (−s) , s ≥ 01tt ∼ Gt (t) = exp −
, t≥0
β
β

(60)

右辺の各項をとおくと、これらは(61)の指数分布に従って分布する。両者を足した x = s + t の分布は (60), (61) 式の畳み込み



Z x
x
1
exp (−x) − exp −(62)
Gs (s) · Gt (x − s) ds =G(x) =1−ββ018TABLE III. ロスコーン分布（Subtracted Maxwellian）の生成法 30

generate U1 , U2 , U3 ∼ U (0, 1)
generate N ∼ N (0, 1)
x ← − log U1 − β log U2√v⊥1 ← θ⊥ x cos(2πU3 )√v⊥2 ← θ⊥ x sin(2πU3 )pv∥ ← θ∥ 1/2 N
return v⊥1 , v⊥2 , v∥となるから、(57) 式の手順で (52) 式の x を乱数生成できることがわかる。そして x を生成したのち、もう１つの一様乱数 U3 を使って方位角を決定し、v⊥1 と v⊥2 を生成することは容易である。v∥ 成分も含めた速度３成分の生成手順を表 III に示す。
表 III のアルゴリズムで、106 個の粒子を使って β = 0.5 の subtracted Maxwellian を乱数生成した結果を図 10(a) に示す。位相空間密度を v⊥ -v∥ 平面に表示している。分布関数の内側の特徴を拡大して見るため、熱速度を θ∥ = θ⊥ = θ = 1.5 とおいた。v∥ 方向に穴が空いた、ロスコーン型の分布関数を生成できていることがわかる。
最後に ∆ ̸= 0 のケースを考える。(50) 式の v⊥ 成分は、∆ の重みの Maxwell 分布と、1−∆の重みの通常の subtracted Maxwellian が重なっている。そこで、∆ の確率で bi-Maxwellianを生成し、残り 1 − ∆ の確率で subtracted Maxwellian を生成してやれば良い。そのようなロジックは、新しい一様乱数 U4 を使った条件分岐で実装することができる。しかし本稿では、技巧的だが一様乱数 U4 を使わずに済む方法を紹介する。(59) 式を



U2
t ≡ −β log min
,1
1−∆

(63)

とおくと、∆ の確率で t = 0 となり、1 − ∆ の確率で t は指数分布に従って分布する。言い換えると、∆ の確率で β = 0 相当である bi-Maxwellian を生成し、残り 1 − ∆ の確率でsubtracted Maxwellian を生成できるわけである。つまり、(57) 式を



U2
,1
x ← − log U1 − β log min
1−∆

(64)

と置き換えれば、(50) 式に従う分布を生成することができる。この修正法を使って、β =0.5, ∆ = 0.1 の subtracted Maxwellian を生成した結果を図 10(b) に示す。図 10(a)（∆ = 0）と比べて、ロスコーン内に 10%ほどプラズマが埋まっていることがわかる。

19

Loss-cone(AK) distribution
6

Loss-cone(AK) distribution
6

10 2

2
00

2

v

4

6

10 3
10 4
10 5
10 6

4

v

v

4

2
00

2

v

4

6

10 2
10 3
10 4
10 5
10 6

FIG. 10. (a) subtracted Maxwellian（β = 0.5, ∆ = 0.0）のモンテカルロ計算結果。v⊥ -v∥ 平面での位相空間密度。(b) subtracted Maxwellian（β = 0.5, ∆ = 0.1）のモンテカルロ計算結果。

IV.

相対論的マクスウェル分布

プラズマの温度 T が静止質量エネルギー mc2 と比べて無視できない場合（電子の場合、512 keV が静止質量エネルギーに相当する）には、相対論的なプラズマ分布を考慮する必要がある。その中でも最も基本的な分布は、マクスウェル分布を拡張した相対論的マクスウェ
ル分布である。この分布を Maxwell–Jüttner 分布 11 あるいは Jüttner–Synge 分布と呼ぶこともある。相対論的マクスウェル分布の関数形は



N
γmc2
f (p)d p =
exp −
d3 p
2
2
4πm cT K2 (mc /T )
T
3

(65)

である。p = mγv は粒子の運動量、γ = [1 − (v/c)2 ]1/2 はローレンツ因子、Kn (x) は第２種
修正 Bessel 関数である。

(65) 式を極座標 (p, θ, φ) で書き直すとN
exp −f (p)dp = 2
m cT K2 (mc2 /T )!
p
1 + (p/mc)2
p2 dpTとなる。以降、t = T /mc2 とおいたうえで、簡単のため m = c = 1, N = 1 とする。
!
p γ1 + p2
1
1
f (p)dp =exp −p2 dp =exp −p2 dp
tK2 (1/t)ttK2 (1/t)
t

(66)

(67)

ここで、文献 25 を参考に相対論的マクスウェル分布の基本量を確認しよう。Digital Li-

brary of Mathematical Functions 19 の公式 10.32.8 より、√Z ∞π(x/2)ν
exp (−x cosh(s)) sinh2ν (s)ds
Kν (x) =Γ(ν + 1/2) 0
20

(68)

この (68) 式の両辺を xν で割って x で微分する。さらに Bessel 関数の基本関係式
d Kn (x)Kn+1 (x)
=−
n
dx
x
xn
を使うと、以下の有用な関係式を得ることができる。25√Z ∞π(x/2)ν
Kν+1 (x) =exp (−x cosh(s)) sinh2ν (s) cosh(s)ds
Γ(ν + 1/2) 0

(69)

(70)

次に p = sinh(s) とおくと、容易に dp = cosh(s)ds, γ = cosh(s) が得られる。(68) 式と (70)
式から

√

Z
 γ
dp
π(1/2t)ν ∞Kν (1/t) =
exp −
p2ν
Γ(ν + 1/2) 0
t
γ
√
Z


∞
γ 2ν
π(1/2t)ν
exp −Kν+1 (1/t) =p dp
Γ(ν + 1/2) 0
t

(71)
(72)

の 2 式が得られる。(72) 式に ν = 1 を代入すると
Z
 γ
1 ∞
exp −
p2 dp
(73)
K2 (1/t) =t 0tR∞となるから、(67) 式の規格化定数の妥当性 0 f (p)dp = 1 を確認することができる。次に、相対論的マクスウェル分布の圧力を考える。気体の圧力は特定方向へのモーメントフラック
ス、言い換えると 3 方向へのモーメントフラックスの総和の 1/3 に相当するから、(71) 式に

ν = 2 を代入したものを用いて
Z
Z
1 ∞
1 ∞
2 dp
P =
f (p) pf (p) · (pv)dp = t
=
3 0
γ
3 0

(74)

となる。これは規格化前の単位系の P = N T に相当する。
プラズマ粒子シミュレーションやモンテカルロシミュレーションで相対論的マクスウェ
ル分布を乱数生成するために、これまでさまざまな数値解法が提案されてきた。6,21–23,28,29 これら全てが棄却法を使ったアルゴリズムである。本稿では Swisdak 法 23 と修正 Canﬁeld 法29の 2 つのアルゴリズムを紹介する。

A.

Swisdak 法
Swisdak 23 は、対数凸な分布に対して適用できる Devroye 7 の棄却法を相対論的マクスウェ

ル分布に応用した。この手法を図 11 を用いて解説する。温度 t = 1 の相対論マクスウェル
分布の理論曲線を黒線で、棄却法のサンプリング関数を赤線で示している。図の縦軸を log
で取っていることに注意する。

21

<latexit sha1_base64="STuflwa6ATi0xbdS/9cSbAfAZUM=">AAADY3ichVJNb9NAEH2JC7Qp0M8DUoVkiCoVCUUbt03aWwWXHosgbaU0imxnG6z6S7YTVEX+A733CueKP8If4MAP4NBzD71w4O3GOVQoZC3bM2/em9ndGSf2vTQT4lepbMw9evxkfqGy+PTZ86XlldXjNBokrmy5kR8lp46dSt8LZSvzMl+exom0A8eXJ87FexU/Gcok9aLwU3YZy05g90Pv3HPtjFB7K+4Gb83zbvDG7C5XRW1/r2HtWqaoCdG0thvKsJo71rZZJ6JWFcU6ilZKr3CGHiK4GCCARIiMtg8bKZ826hCIiXUwIpbQ8nRcIkeF2gFZkgyb6AW/fXrtAg3pq5ypVrus4vNNqDSxKX6KG3Enfojv4lb8mZprpHOovVzy74y1Mu4uXb34eE/VdJ1Dfo5N8h1duUcrYDTD55nVJrz/6RN9Qokv+mSBjoVkjRgbEnepVvfX0cgkx+Qcyhuhynhe3OS0XH1Wsmnn2lN36uj+5DN30WOvhv/oTHqu7nk2o3JScB7qz/QMjPuq9BVO3WS0zOnGsVWrN2q7H3aqB++K+ZvHBl5ji3fQxAEOcYSW3sU1vuJb+bexaKwa62NquVRo1vBgGS//AkF6svc=</latexit>(pm , fm ) (xR , fm )<latexit sha1_base64="6N6jP7oxZtFHzvtBAADR7k1oo2Y=">AAADY3ichVLLTttQED2Jyys8ApRFpQrJJUICCUU3VYEtgk2XtBBAClFkO5dg4ZdsJ4Ai/wD7bmGN+JH+AAs+oIuuWbDpguMbZ0GjkGvZnsc5M3NnxgwcO4qFeMrltQ9j4xOTU4Xpmdm54vzC4lHkt0NLVi3f8cMT04ikY3uyGtuxI0+CUBqu6chj82Iv9R93ZBjZvncYXwey7hotzz6zLSOmqbZ21fi5oZ813HW9MV8SZaGOPihUMqGE7Oz7C7kvOEUTPiy04ULCQ0zZgYGITw0VCAS01dGlLaRkK79EggK5baIkEQatF/y2qNUyq0c9jRkptsUsDt+QTB2r4lHci2fxWzyIv+Lf0FhdFSOt5Zp/s8eVQaN48+nghazhPJP4BKvEmypzk5JLb4zzkdn6uPf4obqhxKW6mat8HlFd+jq0W2Sn/asrSz9G/x6p1kWJ/iTr5LBYLWYyKCdKS3tqqvkkI6tocladAZ5OzVIzj0dkDjPMW/6p2oHeXFN+gVtX+X/HBoWjr+XKVnnzx7fSzm62f5P4jBWssQfb2MF37KOqqviFW9zl/2jT2qK21IPmcxnnI94cbfkVQ9qyoA==</latexit><latexit sha1_base64="njXtUlrm7MlupaKO7HKcv0JXBzg=">AAADY3ichVLLTttQED2JgUJ4hEcXlSok0wgJpCq6QTy2iG666IIKAkghimznEiz8ku0EUOQfYM8W1ogf6Q900Q9gwbqLblhwfOMsIEpzLdvzOGdm7syYgWNHsRB/cnltbHziw+RUYXpmdq44v7B4FPnt0JJVy3f88MQ0IunYnqzGduzIkyCUhms68ti8+Jb6jzsyjGzfO4yvA1l3jZZnn9mWEdNUW7tq/PiqnzXcdb0xXxJloY4+KFQyoYTs7PsLuRWcogkfFtpwIeEhpuzAQMSnhgoEAtrq6NIWUrKVXyJBgdw2UZIIg9YLflvUapnVo57GjBTbYhaHb0imjlXxWzyIv+KXeBTP4mVorK6KkdZyzb/Z48qgUbz5dPCPrOE8k/gEq8SbKnOTkktvjPOR2fq4//FDdUOJS3UzV/k8orr0dWi3yE77V1eWfoz+PVKtixL9SdbJYbFazGRQTpSW9tRU80lGVtHkrDoDPJ2apWYej8gcZpi3/FO1A725pvwCt67yfscGhaONcmW7vPVzs7S7l+3fJD7jC9bYgx3s4jv2UVVV3OIO9/knbVpb1D72oPlcxlnCm6MtvwIx8rKa</latexit>(xL , fm )<latexit sha1_base64="mxsGFmkIg2bi4d0meS8vf0ccynU=">AAADZXichVLLTttQED2JaUvTUqBUCIlF3UZIIKH0OmoD2SHYdMkrgARRZDs3qYVftZ1UKPIP9Ae6bZeoP9IfYMEXINZdsGHB8Y2zQFWaa9mexzkzc2fGCl0nToS4LhS1qSdPn00/L714OfNqdm7+9VEc9CJbNuzADaITy4yl6/iykTiJK0/CSJqe5cpj63wn8x/3ZRQ7gX+YXISy6Zld3+k4tpnQ1FwNW/vreqflfZBremuuLCqCp1bTM8HYFAaFen2zWq3rhnIJUUZ+doP5wjucoY0ANnrwIOEjoezCRMznFAYEQtqaGNAWUXKUXyJFidweUZIIk9ZzfrvUTnOrTz2LGSu2zSwu34hMHSviSlyKv+KP+C1uxf3YWAMVI6vlgn9ryJVha/b70sEdWeN5FvEpVoi3VOY2JY/eBF8mZhvh/seP1A0lvqmbecrnEzWgr0+7TXbWv6ayjGKM7pFpA5TpT/NOjovVZSaTcqq0rKeWmk86sYo2Z9X/h6dTs9XMkwmZoxzzmH+mdmA414xf4taNVksfLxxVK0at8mnvY3lrO9+/aSzjPVbZgw1s4TN20WAVX/EDP/GreKPNaG+0xSG0WMg5C3h0tLcP5JezfQ==</latexit>(pR , fm /e)<latexit sha1_base64="AMiQyHR/E+Eb5WGlYknI+9WkhwY=">AAADZXichVLLTttQED2JoaVpKa8KIbGoS4REpSpcRzSQHaIbFl1QQQAJosh2LsHCr9pOKhT5B/oD3bZLxI/wAyz6BVXXLLph0eMbZ4FQyLVsz+OcmbkzY4WuEydC/C4UtYnJZ8+nXpRevpp+PTM7N38YB93Ilg07cIPo2DJj6Tq+bCRO4srjMJKmZ7nyyLr4lPmPejKKncA/SC5D2fTMju+cObaZ0NRcC1ufP+hnLW9dvtdbs2VRETy1mp4JxpYwKNTrW9VqXTeUS4gy8rMXzBXe4RRtBLDRhQcJHwllFyZiPicwIBDS1kSftoiSo/wSKUrkdomSRJi0XvDboXaSW33qWcxYsW1mcflGZOpYFbfiStyJG3Et/or7kbH6KkZWyyX/1oArw9bM96X9f2SN5lnEp1gl3lKZ25Q8ehOcj802xD3Fj9QNJb6pm3nK5xPVp69Hu0121r+msgxjDO+RaX2U6U/zTo6K1WEmk3KqtKynlppPOraKNmfVe8TTqdlq5smYzFGOecg/VTswmGvGL3HrhquljxYOqxWjVvn4ZaO8vZPv3xSWsYI19mAT29jFHhqs4it+4Cd+Ff9o09obbXEALRZyzgIeHO3tf9Kjs3c=</latexit>(pL , fm /e)FIG. 11. 相対論的マクスウェル分布（t = 1; 黒線）に対する、Swisdak 法 23 のサンプリング関数（赤線; 本文参照）とモンテカルロ計算結果。

簡単のため、相対論的マクスウェル分布の (67) 式の比例係数を落として!
p
1 + p2p2
f (p) = exp −t(75)

とおく。この分布の最頻値（ピーク位置）pm は


√
p2m = 2t t + 1 + t2

(76)

fm = f (pm )(77)

であり、分布の最大値をとおく。さらに、図 11 のように、分布関数の左側と右側のスロープ上に２つの点 pL , pR を考える。それぞれの点で (log f ) の勾配を得ておく。

ppt 1 + p2
f (p)pλ(p) = ′=f (p)2t 1 + p2 − p2
22

(78)

ここで、以下のような 3 つのサンプリング分布を考える。

p − pLhL (p) = f (pL ) exp(−∞ < p < xL )λL

hC (p) = fmhR (p) = f (pR ) exp −p − pRλR

(79)(xL < p < xR )(80)

(xR < p < ∞)(81)

これらを図 11 に赤線で示す。この図では、2 つの指数分布は p = pL , pR で相対論的マクス
ウェル分布に接する直線になる。対数空間では相対論的マクスウェル分布は上に凸であるため、接点以外で接することはない。xL , xR は各分布の交点の座標で

fm
f (pL )
fm
xR = pR − λR log
f (pR )
xL = pL + λL logである。各区間の面積とその総和は
Z xL
hL (p)dp = fm λL
SL =
−∞
Z xR
SC =
hC (p)dp = fm (xR − xL )

(82)
(83)

(84)
(85)

x
Z L+∞

hR (p)dp = fm λR

SR =
X

(86)

xR

S = SL + SC + SR = fm (xR − xL + λL + λR )(87)であるから、以下の確率で 3 つのサンプリング分布に処理を分配してやれば良い。

qL =

λL
,
x R − x L + λL + λ R

qR =

λR
,
x R − x L + λ L + λR

(88)

qC = 1 − qL − qR

(89)

Devroye 7 によると、2 つの接点がf (pL ) = f (pR ) =

fm
e

(90)

P
を満たすとき、 S を最小に抑えることができる。このとき、(82), (83) 式および (88) 式も
x L = p L + λL ,
qL =

λL
,
pR − pL23

x R = p R − λR ,

(91)

λR
,
pR − pL

(92)

qR =



p2m ← 2t t +√TABLE IV. Swisdak 法
1 + t2compute pL and pR for Eq. (90)f (pR )f (pL )λL ← ′, λR ← − ′f (pL )f (pR )x L ← pL + λ L , x R ← pR − λ R ,λL
λRS ← pR − p L , q L ←, qR ←, qC ← 1 − qL − qR
S
S
repeatgenerate X1 , X2 ∼ U (0, 1)
if X1 ≤ qC then// Central box

p ← xL + (xR − xL )(X1 /qC )if X2 < f (p)/f (pm ) breakelse if X1 ≤ qC + qL then// Left tailU ← (X1 − qC )/qLp ← xL + λL log U
if 0 < p and U X2 < f (p)/f (pm ) breakelse// Right tailU ← (X1 − qC − qL )/qRp ← xR − λR log U
if U X2 < f (p)/f (pm ) break
endif
end repeat
generate X3 , X4 ∼ U (0, 1)
px ← p (2X3 − 1)ppy ← 2p X3 (1 − X3 ) cos(2πX4 )ppz ← 2p X3 (1 − X3 ) sin(2πX4 )return px , py , pzのようにシンプルになる。
指数分布 f (x) = exp (−x/λ) に従う乱数は一様乱数 U を用いて

x ← −λ log U
24

(93)

のようにすれば生成できる。従って、少しの変更で (79), (81) 式のサンプリング分布を生成することができる。あとは棄却法を使って、サンプリング分布から目的の分布（(75) 式）を切り出してやれば良い。棄却法の採択率は

eﬀ sw (t) =1
1=
SL + S C + S R
fm (pR − pL )

(94)

である。
一連の乱数生成アルゴリズムを表 IV にまとめた。理解しやすいよう、原論文から表記を変更している。このアルゴリズムの注意点は次の 3 点である。1 つめは、区間を分配するときに使った一様乱数 X1 を、例えば (X1 /qC ) のような形で再利用していることである。これは、乱数の利用個数を減らして計算コストを削減する技巧である。2 つめは、左側のブロック（left tail）で p < 0 が生成される可能性があることである。プログラミング上の扱いが厄介なので、ここでは明示的に p < 0 を省いている。3 つめは、何らかの手段で (90) 式を満たす pL , pR を求めなければならないことである。Swisdak 23 は root ﬁnder で (90) 式の数値解を求めており、この点が Swisdak 法の隠れたコストだと言える。

Swisdak 法で、106 個の粒子を使って相対論的マクスウェル分布を乱数生成した結果を、
図 11 の青色のヒストグラムに示す。数値解が黒線の厳密解とよく一致していることがわかる。このときの棄却法の採択率は 89.9% であった。

B.

修正 Canﬁeld 法
本項では著者ら 29 が最近提案した修正 Canﬁeld 法を紹介する。静止質量で規格化された粒子の運動エネルギーを x とおく。

x≡

Ekin
= γ − 1,
mc2

(95)

このとき (67) 式は以下のようになる。

f (x)dx =pe−1/t −x/te
(1 + x) x(x + 2) dx
tK2 (1/t)
25(96)

次に、重み係数 wi (t)、重み総和 S(t)、規格化された重み係数 πi (t)、棄却関数 R(x) を定義する。a = 0.56、b = 0.35 は定数である。Canﬁeld 論文 6 に倣って添字を i = 3, 4, 5, 6 とする。
√
√
√
3b πt
w3 (t) = π, w4 (t) = a 2t, w5 (t) =, w6 (t) = (2t)3/2 ,
(97)
2
√
6
X
√
√
3b πt
S(t) =wi (t) = π + a 2t ++ (2t)3/2
(98)
2
i=3

wi (t)
πi (t) =S(t)6X(i = 3, 4, 5, 6),πi (t) = 1,

(99)

i=3

√
(1 + x) x + 2√
R(x) ≡ √2 + ax1/2 + b 2x + x3/2

(100)

さらに形状ファクター k 、スケールファクター λ のガンマ分布 Ga(x; k, λ) を考える。

Ga(x; k, λ) =(97)–(101) 式を用いて (96) 式を√ −1/tteS(t)f (x)dx = √2K2 (1/t)1
xk−1 e−x/λk
λ Γ(k)!
i 
R(x) dxπi (t) Ga x; , t
2
i=3

6
X

(101)



と書き換える。さらに R(x) を取り除いたサンプリング分布 F (x) を定義する。
!
√ −1/t6 i teS(t) XF (x)dx = √πi (t) Ga x; , tdx
2
2K2 (1/t) i=3

(102)

(103)

P6

i=3 πi (t) = 1 となるよう πi (t) を定義したため、(102),(103) 式の右辺の大括弧部は、４種類のガンマ分布が確率 πi (t) で重ね合わさっている状態に相当する。そこで、粒子１つ１つについて、一様乱数 X1 ∼ U (0, 1) と確率 πi (t) を使って、ガンマ分布を１つを選び、そのガンマ分布に従う乱数を生成すれば、(103) 式の分布を得ることができる。ガンマ分布の乱数生成法は、数値計算の教科書や文献 29 の Appendix A などを参照されたい。この確率 πi (t)は、温度 t に応じて図 12 のように変化する。t ≪ 1 または t ≫ 1 の極限では、粒子分布は

k = 3/2 あるいは k = 3 のガンマ分布に漸近する。
次に、得られた粒子分布に対して棄却法を適用する。(100) 式の棄却関数は 0.957 ≲ R(x) <1 の範囲で変化する。そこで一様乱数 X2 ∼ U (0, 1) を用いてX2 < R(x)(104)

を満たすかどうかを判定する。技巧的だが、0.957 ≲ R(x) であることがわかっているので

X2 < 0.95 or X2 < R(x)
26

(105)

Probability

1.0
0.8
0.6
0.4
0.2
0.0
10 3

4

6

3

5

10 1
101
TM/mc2

103

FIG. 12. 確率密度 π3 , π4 , π5 , π6 の温度依存性という条件を使っても良い。こうすると、95% の場合は前半の条件だけで判定できて、残り5% の場合だけ棄却関数 R(x) を呼び出せば良いので、計算を高速化することができる。この技法を「圧搾法」と呼び、地味だが有効な技法である。これらの条件を満たしたときに、x の値を採用する。条件を満たさなかった場合は、粒子を棄却して X1 の生成からやり直す。
これらの手順を繰り返すことで (102) 式の分布を得ることができる。さらにp
p = x(x + 2)

(106)

を経て運動量 p の分布（(66) 式）を得ることができる。最後に一様乱数 X3 , X4 ∼ U (0, 1) を使って p を球面上に展開し、x, y, z 成分を得る。

p = p (2X3 − 1) xppy = 2p X3 (1 − X3 ) cos(2πX4 )p pz = 2p X3 (1 − X3 ) sin(2πX4 )(107)

ここまでの流れを疑似コードにまとめたものが表 V である。

(102) 式の分布を作るために (103) 式のサンプリング分布を生成したわけだから、棄却法の採択効率はR∞√ 1/t
f (x)dx2e K2 (1/t)
0.
eﬀ MC (t) = R ∞= · · · = √ √√√
F (x)dxt
π + a · 2t + b · 3 πt + (2t)3/2
0
2

27

(108)

TABLE V. 修正 Canﬁeld 法 29
compute π3 , π4 , π5 for given t using Eqs. (97)–(99)repeatgenerate X1 , X2 ∼ U (0, 1)
X1 < π 3then i ← 3

elseif X1 < π3 + π4

then i ← 4

ifelseif X1 < π3 + π4 + π5 then i ← 5else i ← 6
endifgenerate x ∼ Ga(i/2, t)
until X2 < 0.95 or X2 < R(x)// (105) 式（圧搾法）

generate X3 , X4 ∼ U (0, 1)pp ← x (x + 2)
px ← p (2X3 − 1)ppy ← 2p X3 (1 − X3 ) cos(2πX4 )ppz ← 2p X3 (1 − X3 ) sin(2πX4 )return px , py , pz

である。

C.

比較
棄却法の採択率の相対論的温度 t = T /mc2 に対する依存性を図 13 に示す。Swisdak 法（(94) 式）
、修正 Canﬁeld 法（(108) 式）や、他のいくつかの方法（Sobol 法 22 、Canﬁeld 法 6 ,Schnittman–Krolik 法 21 ）の採択効率を比較している。Swisdak 法は非相対論側で 88%〜相対論側で 90%、修正 Canﬁeld 法は全パラメーター領域で 95% 以上の良い効率が出ていることがわかる。その他の手法は、温度パラメーターに応じて得意・不得意が大きいうえ、生成効率もこの 2 手法に及ばない。
実際の利用状況では、採択率に加えて、さまざまな要素によってコードのパフォーマン
スが決まるはずである。例えば、比較的計算負荷が大きな、乱数の生成個数も大事な要素である。Sobol 法と修正 Canﬁeld 法は乱数を 6 つ使うのに対して、Swisdak 法は乱数を 4 つし28

Efficiency

1.0
0.8
0.6
SZ & Nakano (2022)
0.4
Swisdak (2013)
Schnittman+ (2013)
0.2
Canfield+ (1987)
Sobol (1976)
0.0
10 3 10 2 10 1 100 101 102 103
TM/mc2

FIG. 13. 相対論的マクスウェル分布の乱数生成効率。Swisdak 法 23（紫; (94) 式）、修正 Canﬁeld 法（赤; (108) 式）および他の手法（Sobol 法 22,28 , Canﬁeld 法 6 , Schnittman–Krolik 法 21 ）の採択効率を示す。

か使わない。一方で Swisdak 法は、前述のように root ﬁnder を使う必要があり、その実装と計算負荷が未知数である。これらの数値解法の性能を、さまざまな利用状況と実行環境で検証すべきであろう。

V.

相対論的カッパ分布
カッパ分布を拡張した相対論的カッパ分布の表式は 9


1
N
Γ
κ
+
κ
2


fRK (p)d3 p ≡κTκ
3/2(2πmκTκ ) (κ + 1) Γ(κ − 2) 2 F1 − 23 , 52 ; κ + 12 ; 1 − 2mc2−(κ+1)(γ − 1)mc2
× 1+
d3 p
κTκ

(109)

とされている。ここで Nκ はプラズマ密度、Tκ は特徴温度、2 F1 は超幾何関数、κ は高エネ
ルギーのべき指数を決めるパラメーターである。以降、本稿では簡単のため、m = c = 1、t ≡ Tκ /mc2 とおく。κ パラメーターは、非相対論の場合は κ > 3/2 であったのに対して相対論では κ > 3 であることに注意する。これは、相対論的高エネルギー領域では p ≈ γ と29

TABLE VI. 相対論的カッパ分布の乱数生成法 29
compute π3 , π4 , π5 for given κ, t using Eqs. (112)–(114)
repeatgenerate X1 , X2 ∼ U (0, 1)
X1 < π 3then i ← 3

elseif X1 < π3 + π4

then i ← 4

ifelseif X1 < π3 + π4 + π5 then i ← 5else i ← 6
endifgenerate X3 ∼ Ga(i/2, 1), X4 ∼ Ga(κ + 1 − i/2, 1)
X3
x ← κt ×
X4
until X2 < 0.95 or X2 < R(x; a, b)
generate X5 , X6 ∼ U (0, 1)pp ← x (x + 2)
px ← p (2X5 − 1)ppy ← 2p X5 (1 − X5 ) cos(2πX6 )ppz ← 2p X5 (1 − X5 ) sin(2πX6 )なるため、エネルギー密度R∞0γfRK (p)dp ∝R∞0p3−κ dp が κ = 3 で発散するからである。

本節では、我々が最近、修正 Canﬁeld 法を応用して開発した手法を紹介する。29 本稿執筆時点では、これは相対論的カッパ分布向けに作られた唯一の解法である。(109) 式を極座標で書き表すと、以下のようになる。

γ−1
fRK (p)dp = A(κ, t) 1 +κt−(κ+1)4πp2 dp(110)

紙面節約のため、規格化定数を A(κ, t) とおいた。さらに前節と同様に運動エネルギー

x≡

E
=γ−1
mc2

(111)

を導入し、重み関数 wi (κ, t) とその総和 S(κ, t)、確率関数 πi (κ, t) を定義する。多くが κ お
よび t の関数であることに注意する。定数 a = 0.56、b = 0.35 および棄却関数 R(x)（(115)30

式← (100) 式）は、前節のものと同一である。
√
2π 
1
w3 (κ, t) =
Γ κ−
,
2√
2

3b 2π
3
w5 (κ, t) =(κt) Γ κ −
,
4
2
6
X
S(κ, t) ≡wi (κ, t),√
w4 (κ, t) = a κt Γ(κ − 1),
w6 (κ, t) = 2(κt)3/2 Γ(κ − 2),

(112)
(113)

i=3

πi (κ, t) =wi (κ, t)S(κ, t)(i = 3, 4, 5, 6),

6X√(1 + x) x + 2√.
R(x) ≡ √2 + ax1/2 + b 2x + x3/2

πi (κ, t) = 1,

(114)

i=3

(115)

(111)–(115) 式と、前の節（(34) 式）で触れた一般化ベータプライム分布 B ′ (x) を使うと、(110) 式を以下のように書き換えることができる。同様に、R(x) を含まないサンプリング分
布 FRK (x) を定義する。· · · の部分には複雑な係数が現れるが、今回のアルゴリズムでは使用しないので省略した。

fRK (x)dx =···FRK (x)dx =

···

!

i
i
πi (κ, t) B′ x; , κ + 1 − , 1, κt
R(x) dx
2
2
i=3
!
6
 i

X
i
πi (κ, t) B′ x; , κ + 1 − , 1, κt
dx.
2
2
i=3


6
X

(116)
(117)

FRK (x) は πi の確率（重み）で４種類のベータプライム分布を重ね合わせた分布であるから、これを生成するのは容易である。また、ベータプライム分布 B′ (x; 2i , κ + 1 − 2i , 1, κt) に従う乱数 x は、ガンマ分布 Gamma(x; k, λ) に従うガンマ乱数 XGa(k,λ) を２つ使って

x = κt ×

XGa(i/2,1)
XGa(κ+1−i/2,1)

.

(118)

とすれば求められる。29 そして前節と同様に、棄却関数 R(x) を使って fRK (x)（(116) 式）を得たのち、x を p に変換して、p を空間３成分に分配すれば良い。ここまでのロジックを表VI の疑似コードにまとめた。乱数生成効率は、前節で解説した Maxwell–Jüttner 分布の乱数生成効率と同様、95 ∼ 100 % である。これは性質のよい棄却関数 R(x)（(100)・(115) 式）を利用したためである。生成効率の理論式はかなり複雑なので本稿では省略する。
この方法で相対論的カッパ分布を乱数生成した結果を図 14 に示す。相対論的カッパ温度

t = 1.0, カッパ指数 κ = 3.5 のカッパ分布を 106 個の粒子を用いて乱数生成した。粒子の運動量 |p| のヒストグラム（青）が、黒線の厳密解（(110) 式）とよく一致する。また、右側の高エネルギー部分に f (p)dp ∝ p−(κ−1) dp（この場合は p−2.5 ）の「べき」乗成分が綺麗に見え31

Normalized density

10 1

p 2.5

10 3
5
10 10
1

100

101
p/mc

102

103

FIG. 14. 相対論的カッパ分布の乱数生成結果。運動量ヒストグラム（青）と厳密解（黒; (110) 式）。

ている。低エネルギーの熱的コア成分から高エネルギーの「べき」乗成分に変化する、カッパ分布の特徴がよく現れている。

VI.

謝辞

本研究は、科学研究費助成事業 若手研究 (B) 25871054 および基盤研究 (C) 17K05673,(C) 21K03627, (S) 17H06140 の助成のもとで行われた。

REFERENCES
1

R. F. Abdul and R. L. Mace, Phys. Plasmas 22, 102107 (2015).

2

Y. Asano et al., J. Geophys. Res. 113, A01207, doi:10.1029/2007JA012461 (2008).

3

M. Ashour-Abdalla, C. F. Kennel J. Geophys. Res. 83, 1531 (1978).

4

R. W. Bailey, Math. Comp. 62, 779 (1994).

5

G. E. P. Box and M. E. Muller, Ann. Math. Stat. 29, 610 (1958).

6

E. Canﬁeld, W. M. Howard, and E. P. Liang, Astrophys. J. 323, 565 (1987).
32

7

L. Devroye, Non-Uniform Random Variate Generation, Springer-Verlag (1986), available
at http://luc.devroye.org/rnbookindex.html.

8

W. C. Feldman et al., Phys. Rev. Lett. 49, 199 (1982).

9

L. Han-Thanh, K. Scherer, and H. Fichtner, Phys. Plasmas 29, 022901 (2022).

10

M. Hofert, R Journal 5, 129 (2013)11

F. Jüttner, Ann. Phys. 339, 856 (1911).

12

D. P. Kroese, T. Taimre, and Z. I. Botev, Handbook of Monte Carlo methods, John Wiley& Sons. (2011), （邦訳：D. P. クローゼ、T. タイマー、Z. I. ボテフ著／伏見正則・逆瀬川浩孝監訳、モンテカルロ法ハンドブック、朝倉書店 2014）13

G. Livadiotis (ed.), Kappa Distributions: Theory and Applications in Plasmas, Elsevier,
Amsterdam (2017).

14

M. D. Montgomery J. R. Asbridge, and S. J. Bame, J. Geophys. Res. 75, 1217 (1970).

15

M. Oka, S. Ishikawa, P. Saint-Hilaire, S. Krucker, and R. P. Lin, Astrophys. J. 764, 4(2013).

16

M. Oka, S. Krucker, H. S. Hudson, and P. Saint-Hilaire, Astrophys. J. 799, 129 (2015).

17

M. Oka et al., Phys. Plasmas 29, 052904 (2022).

18

S. Olbert, in Physics of the Magnetosphere, eds. R. L. Carovillano, J. F. McClay, & H. R.
Radoski (Astrophysics and Space Science Library, Vol. 10; Dordrecht: Reidel), 641 (1968).

19

F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C.
W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds., NIST
Digital Library of Mathematical Functions. https://dlmf.nist.gov/, Release 1.1.12 of2023-12-15.

20

L. A. Pozdnyakov, I. M. Sobol, and R. A. Sunyaev, Astrophys. Space Phys. Rev. 2, 189(Translation) (1983).

21

J. D. Schnittman, and J. H. Krolik, Astrophys. J. 777, 11 (2013).

22

I. M. Sobol, in Monte Carlo Methods in Computational Mathematics and Mathematical
Physics (Novosibirsk, 1976), pp. 24–29 [in Russian] (1976).

23

M. Swisdak, Phys. Plasmas 20, 062110 (2013).

24

D. Summers and R. M. Thorne, Phys. Fluids B 3, 1835, doi:10.1063/1.859653 (1991).

25

A. Suzuki, Phys. Plasmas 15, 072107 (2008).

26

M. F. Thomsen, H. C. Barr, S. P. Gary, W. C. Feldman, and T. E. Cole, J. Geophys. Res.
88, 3035 (1983).
33

27

V. M. Vasyliunas, J. Geophys. Res. 73, 2839 (1968).

28

S. Zenitani, Phys. Plasmas 22, 042116 (2015).

29

S. Zenitani and S. Nakano, Phys. Plasmas 29, 113904 (2022).

30

S. Zenitani and S. Nakano, J. Geophys. Res. 128, e2023JA031983 (2023).

31

S. Zenitani, Res. Notes AAS 8, 30 (2024).

32

四辻 哲章、計算機シミュレーションのための確率分布乱数生成法、プレアデス出版,

ISBN:4903814351 (2010).

34