1. No category

宇宙電磁力学
2D Polar Cap(II)
柴田 晋平 (しばたしんぺい)
1
平成 26 年 6 月 26 日
1
山形大学理学部物理 〒 990 山形市小白川町 1-4-12 Email: [email protected] URL: http://astrwww.kj.yamagata-u.ac.jp tel 0236-28-4552 fax 0236-28-4567
155
表 13.1: 無次元化の表 (添字 u を単位量とする)
13.3.3
物理量
単位
無次元量
t˜
意味
t
r
tu
Lu
∇
E, B
q
r˜
˜
∇
長さ
L−1
u
Bu
qu
˜ B
˜
E,
q˜
m
n
mu
nu
m
˜
n
˜
質量
ρe
j
p
q u nu
q u nu c
mu
電荷密度
f
mu c/tu
ρ˜e
j˜
˜
p
f˜
時間
空間微分
電磁場
電荷
数密度
電流密度
運動量
力
Non-dimensional Form
基礎方程式は、サイクロトロン振動数とプラズマ振動数が等しく 1 になるように単位時間を定め、そし
て、速度単位が光速になるように規格化されている:
√
qu Bu
4πqu2 nu
−1
ωcu ≡
= tu =
, ≡ ωpu
mu c
mu
ctu = Lu
(13.55)
この時点で、単位をどのようにとるかは上の式を満たすようにさえすればどのようにとっても良い (表 13.1)。
初期のプログラムでは磁場の単位を極磁場 B∗ ∼ 1012 G にとっていたがこれは非常に大きな値であった。
このために、サイクロトロン振動数と等しくなるプラズマ振動数は非常に高く、よって密度単位は GJ 密度
に較べると非常に大きな値になっていた。(GJ 密度の無次元表現 ρ˜c は非常に小さな値になる。) また、長
さの単位である Lamore 半径は非常に小さいので無次元表現の polar cap 半径がとてつもなく大きな値に
なってしまう。いっぽう、モデルは基本的に一次元静電問題なので、磁場は粒子の運動を一次元に拘束する
以上に意味が無いので、実際は 1012 G でなくても、103 G でもなんでも良かったのであった。実際、一次元
コード (EM1) では磁場は方程式系に現れず、磁場強度の単位は不要である。今回の円柱モデルでは磁場は
現れるので、磁場の単位は必要ではあるが、それは電磁波の成分としての磁場が重要であって、一次元方向
に粒子を拘束する磁場の強度は、別に B∗ である必要はない1 。
Polar Cap model で重要なのは、グリッドサイズ ∆z 、粒子の慣性長 λi = c/ωp,GJ 、polar cap size Rpc =
(Ω∗ R∗ /c)1/2 R∗ の間に
(13.56)
∆z ≪ λi ≪ Rpc
の関係が成立していることであろう。ここで、ωp,GJ は GJ 密度に対するプラズマ振動数、R∗ は星の半径、
Ω∗ は星の自転角速度である。とくに、慣性長と polar cap size の比は現実に近いものにして計算したい。
また、グリッドサイズが慣性長よりも十分に小さくないとノイズの多い計算になってしまう。
今回の改訂では、以上を考慮して、GJ 密度の M 倍を単位電荷密度にすることにする。M はペア生成の
増幅率 (multiplicity) を意識している。離散化のときに、∆z の出現を避けたいこともあり、単位長さのグ
1 但し、B
z ez
× E のドリフト運動の評価のときのみに注意する必要は有る。
156
リッド幅にすることにする。以上の方針で規格化する。
密度は以下で規格化する：
q u nu = M
対応するプラズマ振動数は、
√
ωpu =
2M ·
Ω∗ B∗
2πc
(13.57)
qu B∗
1/2
· Ω∗ = 4.486 × 1010 sec−1 M1/2 B12 P −1/2 m
¯ −1/2
q
mu c
(13.58)
ここで、
m
¯q =
mu /me
qu /e
(13.59)
単位粒子の質量電荷比 (mass-to-charge ratio) の電子に対する比であり、B12 ,P はそれぞれ、極磁場と周期
を 1012 G および sec で表した値である。また、長さの単位は、
c
−1/2
= 2.019 cm M−1/2 B12 P 1/2 m
¯ 1/2
q
ωpu
Lu =
(13.60)
共回転速度は速度単位 c で規格化されるので u
˜c = r × Ω∗ /c であるが位置ベクトルが Lu というミクロな
長さで定義されるため、無次元形では
u˜c =
Lu
r˜ × Ω∗ ez = r˜ ×
c
と言う形になる。ここで、
(
Lu Ω∗
c
)
˜
ez ≡ r˜ × Ω
˜ = Lu Ω∗ = Ω∗ = √ ϵ
˜ = |Ω|
Ω
c
ωpu
2M
(13.61)
(13.62)
である。また、小さなパラメーター
√
ϵ=
Ω∗
−1/2
¯ 1/2
= 5.979 × 10−10 B12 P −1/2 m
q
qu B∗ /mu c
(13.63)
はマクロな振動数とミクロな振動数 (非相対論的) の比の平方根である。規格化は ωcu = ωpu の下で行われ
たのでこれから直ちに
˜ ∗ = B∗ = √ 1 1
B
Bu
2M ϵ
(13.64)
であることも分かる。
規格化された星の半径と Polar Cap の半径を調べてみると、
ω∗ =
Ω∗ R∗
= 2.094 × 10−4 R6 P −1
c
(13.65)
を導入して、
˜ ∗ ≡ R∗ = Ω∗ R∗ /c = ω∗
R
˜
Lu
Ω∗ Lu /c
Ω
(13.66)
また、
3/2
˜ pc
R
となる。
3/2
ω∗
(2M)1/2 ω∗
1/2 ˜
= ω∗ R
=
∗ =
˜
ϵ
Ω
3
1/2 1/2 3/2 −1 −1/2
= 7.169 × 10 M B12 R6 P m
¯q
(13.67)
(13.68)
157
˜ ∗ 、Polar Cap の半径 R
˜ pc 、光半径 c/Ω∗ Lu に対
ここで、(13.61) を考えると、u
˜ の大きさは、星の半径 R
3/2
して、それぞれ、ω∗ 、ω∗ 、1 となる。
電荷密度の表式は、
˜ c = −∇
˜ = −2Ω
˜ φ ≈ −2Ω
˜ ·E
˜ · (˜
˜B
˜z + Ω
˜ ϖ[
˜ B]
˜B
˜z
ρ˜c = ∇
u × B)
˜ ∇×
(13.69)
˜z = B
˜∗ とすると、
となる。ここで、B
1
M
で予定通り、GJ 電荷密度は単位電荷密度の M−1 である。
ρ˜c = −
(13.70)
質量 m、電荷 q を持つプラズマが GJ 電荷密度を持つときの慣性長は、 無次元形で
√
c/ωp,GJ
ωpu
4πqu2 Mρc /mu
mM
˜
˜
λi =
=
=
=
2
Lu
ωp,GJ
4πq ρc /mu
q˜
(13.71)
である。
次に、電磁場の大きさについて考える。まず、共回転電場は、
˜ c = −˜
˜ = −Ω
˜
˜ ϖe
E
uc × B
˜ φ×B
(13.72)
であるので、
2/3
˜ c| ∼ Ω
˜R
˜ pc B
˜z =
|E
ω∗
(2M)1/2 ϵ
(13.73)
程度の大きさであることが分かる。また、Polar Cap に真空ギャップができたときの典型的な電場の大きさ
′
˜ = ρ˜′ から見積もると、磁場に沿った規格化された長さを z˜ とすると
˜ ·E
は∇
˜
dE
= −˜
ρc
d˜
z
(13.74)
˜gap = −˜
˜B
˜∗ R
˜ pc
E
ρc z˜ = 2Ω
(13.75)
より
˜ c | と同程度である。また、Goldreich-Julian 電流が流れたときに発生するトロイダル磁場の大きさは、
で |E
H
˜ = j˜ より、 Bφ dℓ˜ = 2πϖBφ = π ϖ
˜ ×B
∇
˜ 2 ρ˜c から、
3/2
ω∗
˜φ = 1 ρ˜c ϖ
B
˜ =
2
(2M)1/2 ϵ
(13.76)
である。
パルサーによる加速電圧の上限は、共回転電場と Polar cap 半径の積で見積もることができる：
(
)
Ω∗ Rpc B∗
B∗ R∗3 Ω2∗
µB Ω2∗
ϕmax =
Rpc =
=
(13.77)
2c
2c2
c2
˜ c |R
˜ pc /2 =
ここで、µB = B∗ R∗3 /2 は中性子星の磁気モーメント。同じ結果は、無次元の式で表して、|E
˜ R
˜ pc = Ω
˜R
˜ 2 /2 としても得られる。単位電位は Bu Lu なので、
|˜
u × B|
pc
ϕmax
1 ω∗3
ϕ˜max =
=
Bu Lu
2 ϵ2
となる。
(13.78)
158
無次元形の運動方程式は、
d˜
p
˜ +v
˜
˜ × B)
= q˜(E
dt˜
(13.79)
˜ = p/mu c = mγ
˜ とした。これから、
である。ここで、p
˜ v
dmγ
˜
˜
˜·E
= q˜v
dt˜
(13.80)
∫
となり、積分すると
˜ · d˜
q˜E
s = q˜∆ϕ˜
[mγ]
˜ =
(13.81)
になる。こうして、
γmax = q˜ϕ˜max /m
˜
を得る。その値は、
(
γmax = 1.284 × 10
7
q˜
m
˜
)
B12 R63 P −2 m
¯ −1
q
(13.82)
(13.83)
慣性長と polar cap size を比較すると、
˜ pc
R
1/2
= 2γmax
λ˜i
(13.84)
Multiplicity parameter
(13.85)
maximum avairable voltage (no dimension)
(13.86)
を得る。
以上をまとめると次のようになる：
キーパラメータは
M
ϕ˜max
である。このとき、GJ density は
1
M
対応する慣性長とプラズマ振動数は単位質量、単位電荷密度の粒子に対して
ρ˜c = −
˜i =
λ
√
M
(13.88)
√
Mϕ˜max
(13.89)
1
ω
˜ p,GJ
になる。ポーラーキャップのサイズは
˜ pc = 2
R
である。典型的な電磁場の大きさは
(13.87)
=
√
˜0 =
E
ϕ˜max
M
(13.90)
である。最大電位差に相当するローレンツ因子は単位質量、単位電荷密度の粒子に対して γmax = ϕ˜max で
˜∗ = 1 で代用する。
ある。計算上は便宜上 B
159
13.3.4
ペア形成過程
主なガンマ線放射プロセスは curvature radiation である。質量 m、電荷 e、ローレンツ因子 γ の粒子が
放射する curvature photon の典型的なエネルギー hνc は
Eph =
hνc
3 ¯h/mc 3
=
γ =
2mc2
4 Rc
(
γ
γc
)3
(13.91)
1/6
であり (ここで、γc = 1.337 × 106 R6 P 1/6 )、放射パワーは
2 e2 4
γ
Pc =
3 c3
(
c2
Rc
)2
=
2 e2 c 4
γ
3 Rc
(13.92)
で与えられる。Rc は粒子の軌道の曲率半径。これより、粒子が ℓ 進む間に放射される光子数は
Np =
4 ℓ
α γ
9 Rc
(13.93)
と表せる。ここで、α = e2 /¯
hc = 1/137 は微細構造定数である。
これらの量を数値計算のなかで取り扱うために、数値計算のなかで定義されている polar cap size Rpc と、
単位質量、単位電荷にたいする最大ローレンツ因子 γmax を用いて規格化して、現実世界と simulation の対
応を取ることにする。
放射される曲率光子の典型的エネルギー (13.91) は
(
Eph = Eph,max
Eph,max =
γ
)3
γmax
3¯
h/mc 3
15/2
3
γ
= 4.233 × 106 R6 P −11/2 B12
m
¯ −3
q
4 Rc max
(13.94)
(13.95)
と表せる。
curvature photon がペアを作る mean free path は
ℓp =
[ ]
4.4 ¯h Bq
4
exp
α mc B⊥
3χ
(13.96)
hνc B⊥
2mc2 Bq
(13.97)
で与えられる。ここで、
χ=
であり、Bq = m2 c3 /e¯
h = 4.41 × 1013 G、光子の進行方向と磁場のなす角を θ とすると、B⊥ = B∗ sin θ で
ある。ここで、sin θ = ℓp /Rc と近似して ℓp を求めることにする。(13.96) を、simulation で使いやすいよ
うに ℓp を Ppc で規格化すると、
X≡
3χ
4
3 B∗ Rpc
ℓp
Eph
4 Bq R c
Rpc
ℓp
= Xmax
= Xmax Y
Rpc
=
(13.98)
(13.99)
ここで、Y = ℓp /Rpc はポーラーキャップサイズを単位とした pair mean free path である。また、Xmax は
種光子のエネルギーによってきまる係数で
Xmax = X0 Eph
(13.100)
160
と定義され、
X0 =
3 B∗ Rpc
= 3.561 × 10−6 R6 P −1
4 Bq R c
(13.101)
−1/2
1/2
である。ここで、Rpc ≈ R∗ ω∗ 、Rc ≈ ω∗
で評価した。(13.96) は、
(
)2
(
)2
[ ]
[ ]
ℓp
4.4 ¯
h/mc Bq Rc
1
1
2
Y =
=
exp
= A exp
Rpc
α Rc B∗ Rpc
X
X
(13.102)
ここで、
−5/2
A = 3.385 × 10−7 R6
P 3/2
(13.103)
である。結局、
Y2
=
A exp(1/X)
(13.104)
Y
=
X/Xmax
(13.105)
を連立して解くことなる。(13.104) は、X の増加とともに急激に 0 になる関数であり、Y < 1 のみペアを
polar cap 内に作れることから、X の最小値 Xmin = 1/(− ln A) が存在する。典型的には Xmin = 10−2 で
ある。
一方、(13.105) は、原点から右上がりで、X = Xmax で Y = 1 になる。Xmax は光子のエネルギー Eph
の関数になっていて、たとえば、Eph ≈ 10 のへなへなの光子については Xmax = 10−6 で大変小さい値で、
Xmax < Xmin になってしまい、Y < 1 の解は存在しない。(Rpc より内側ではペアを作らない。) この境
界線は、Xmin = Xmax より、Eph,min = Xmin /X0 であり、このエネルギーより低いエネルギーの光子は
polar cap 内部 (Y < 1) でペアを作らない。
図 13.4 な二つの関数の交点が解である。
根は、xmin < x < xmin にあるとして、二分法あるいは
do i=1,n
ell = x/x_max
x = 1/dlog(ell*ell/ell0)
enddo
のような反復により解くことができる。
二分法のプログラムも作ってから重大なことに気がついた：Eph を与えて ℓ(あるいは Y ) を求めるのは確
かに数値的にしないと解けないので面倒であるけれど、その逆、Y を与えて、Eph を求めることは簡単に
できる。従って、Y (Eph という関数を求めておいて、表または fitting 関数を使って、Eph から Y を求めれ
ば非常に計算時間を短くできる。
(13.104), (13.105), (13.100) より
Eph =
1
X0 Y (2 ln Y − ∈ A)
(13.106)
を得る。ここで、Y の上限は１と考えることにする。これに対応して最小の Eph が存在する。この光子の
エネルギーより低い光子は polar cap 内部でペアを作らない。逆に,Y の最小値はグリッド幅くらいでこれ
を Ymin とすると、対応して Eph のしきい値が存在し、これよりエネルギーの光子はすぐにペアを作ること
になる。ここで、
Eph,min
=
1
X0 (− ln A)
(13.107)
Eph,max
=
1
X0 Ymin (2 ln Ymin − ln A)
(13.108)
161
図 13.4: mean free path ℓ
図 13.5: Eph の関数としての mean free path
である。
上記の典型的値では、図 13.5 のように fit できて、
log Y = 4.99976 − 1.37842 log Eph + 0.0499655(log Eph )2
(13.109)
で十分良い近似になっていると考えられる。
以上まとめると、Eph > Eph,min であれば上の式で mean free path Y = ℓp /Rpc で評価し、Eph < Eph,min
であればガンマ線として出てゆくとすることにしよう。
13.3.5
ふろく：曲率光子の生成
プログラム内で曲率光子の生成ルーチンを設ける。
162
図 13.6: (a) は高田君提供の f (x)/x のグラフ。(b) は N (x)、(c) は N (x) から生成した光子 (200 万光子)。
(e) は (c) の対数スケールプロット
曲率光子のスペクトルを与える関数を高田君より頂く。
dN
= f (x)/x
dx
(13.110)
を与えるサブルーチン FK をファイル curv.f に格納。中央部分は数表なので hokan.f にある多項式による補
間法と併用する。ここで、
∫
∞
f (x) = 2x
K5/3 (ξ)dξ
(13.111)
dN ′
dx
dx′
(13.112)
2x
である。
積分して
∫
x
N (x) =
0
を求めておく。(実際は数表の形で保存 curvN.data し、補間法によって、N → x の変換が出来るようにし
ておく。) 指数関数的に減少する関数なので x < xmax = 7 程度で十分である。一様乱数を [0, N (xmax ) で
振ってそれを N として、N → x の変換してそれを光子のエネルギーとすると、dN/dx が大きいところほ
どそのあたりの x に落ちる確率が高くなる。
以上から、数表 N = N (x) (ファイル curvN.f) を読み込む乱数を振る乱数を数表から補間法で光子のエネ
ルギー x に変換する。というサブルーチン (curfgen) を作成する。今後、これを一回コールするごとに光子
を一つ得ることが出来る。グラフィカルな結果を図??に示す。
normalization ついては、200 万光子の実験で < x >= 0.3059 を得ているのでこれを規格化因子として
使う。