極高エネルギー宇宙線空気シャワーのジオメトリ再構成のためのシャワー構造の研究(A study of the structures of ultra-high energy cosmic ray air
showers for their geometry reconstruction)

理学研究科
数物系専攻
山根涼

目次
第 1 章 極高エネルギー宇宙線
1.1 宇宙線 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.1 宇宙線のエネルギースペクトル . . . . . . . . . . . . . . . . . . . . . . . .
1.1.2 宇宙線の化学組成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.3 宇宙線の加速機構 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 極高エネルギー宇宙線 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.1 GZK カットオフ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.2 磁場による偏向 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 空気シャワーと宇宙線観測 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.1 空気シャワー現象 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.2 空気シャワーの縦方向発達 . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.3 空気シャワーの横方向分布 . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.4 空気シャワー粒子の到来時刻分布 . . . . . . . . . . . . . . . . . . . . . . .
1.3.5 空気シャワーの観測方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4 観測結果の現状 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4.1 エネルギースペクトル . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4.2 質量組成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4.3 到来方向分布 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
1
2
4
8
8
10
12
12
13
15
16
17
19
20
23
23

第 2 章 宇宙線望遠鏡実験
27
2.1 宇宙線望遠鏡実験の装置 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.2 地表検出器 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.2.1 地表検出器基本構成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.2.2 地表検出器電子回路 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.2.3 通信塔基本構成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.2.4 無線 LAN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.2.5 トリガー DAQ サイクル概観 . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.2.6 各検出器におけるデータバッファリングとレベル 1 トリガー . . . . . . . . 33
2.2.7 通信塔におけるレベル 2 トリガー . . . . . . . . . . . . . . . . . . . . . . . 35
2.2.8 境界トリガー . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.2.9 SD アレイのトリガー効率 . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.2.10 モニター情報とモニタリングシステム . . . . . . . . . . . . . . . . . . . . 36
2.2.11 バックグラウンドミューオン信号に対応するエネルギー損失 . . . . . . . . 41
2.3 大気蛍光望遠鏡 (BRM ステーション/LR ステーション) . . . . . . . . . . . . . . . 43
2.3.1 光学系 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
2.3.2 PMT カメラ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
2.3.3 大気蛍光望遠鏡用データ収集システム . . . . . . . . . . . . . . . . . . . . 46

i

2.4
2.5

2.3.4 ハイブリッドトリガーシステム . . . . . . . . . . . . . . . . . . . . . . . .
2.3.5 大気蛍光望遠鏡のキャリブレーション . . . . . . . . . . . . . . . . . . . .
SD アレイと FD による同時検出事象 . . . . . . . . . . . . . . . . . . . . . . . . .
Electron Light Source(ELS) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47
49
55
56

第 3 章 モンテカルロシミュレーション
59
3.1 CORSIKA シミュレーション . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.1.1 シニング . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.1.2 デシニング . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.2 GEANT4 シミュレーション . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
第4章
4.1
4.2
4.3

空気シャワー粒子の到来時刻の遅れ及びその標準偏差の従う関数作成
65
空気シャワー構造 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
シャワー粒子の到来時刻の遅れの関数化 . . . . . . . . . . . . . . . . . . . . . . . 66
シャワー粒子の到来時刻の標準偏差の関数化 . . . . . . . . . . . . . . . . . . . . . 67

第 5 章 ジオメトリ再構成プログラム86
5.1 検出器レベルでの解析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
5.1.1 粒子の到来時刻の決定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
5.1.2 粒子数の決定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
5.2 検出器レベルでのクオリティーカット . . . . . . . . . . . . . . . . . . . . . . . . . 87
5.3 モデル関数によるジオメトリ再構成 . . . . . . . . . . . . . . . . . . . . . . . . . . 88
5.4 イベントレベルでのクオリティーカット . . . . . . . . . . . . . . . . . . . . . . . 91
第6章
6.1
6.2
6.3

6.4
6.5

モンテカルロシミュレーションによる到来方向及びジオメトリ決定精度の評価
93
イベントセット . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
FD の観測結果を用いたエネルギー補正 . . . . . . . . . . . . . . . . . . . . . . . . 93
ジオメトリ再構成プログラム内の種々のしきい値決定 . . . . . . . . . . . . . . . . 93
6.3.1 χ2/dof のしきい値決定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
6.3.2 Ncluster 及び NSD のしきい値決定 . . . . . . . . . . . . . . . . . . . . . . . 95
クオリティーカットにおける χ2 /dof 及び pointing error のしきい値決定 . . . . . . 97
到来方向及びジオメトリ決定精度の評価 . . . . . . . . . . . . . . . . . . . . . . . 99

第 7 章 結論

103

ii

概論

1019 eV 以上のエネルギーを持つ宇宙線 (極高エネルギー宇宙線) は，爆発的な天体現象や現在の宇宙に残された初期宇宙の痕跡粒子を起源に持つと推定され，この起源を解明することは宇宙物理学の大きな課題のひとつである．
これらを解決するために，我々は日米韓露共同で宇宙線望遠鏡実験 (Telescope Array Experiment，TA 実験) を米国ユタ州に建設し，2008 年 3 月より定常運転を続けている．これまで TA 実験では，約 700 km2 の領域に展開された 507 台の地表検出器からなる空気シャワーアレイが常時 97%以上の稼働率 (全 507 台に対する正常安定稼動検出器の割合) で運用されてきた．さらにその周囲の 3 つのステーションに設置された合計 38 台の大気蛍光望遠鏡による観測を加え，到来する宇宙線をこれら 2 種類の検出方法で同時観測している．そして地表検出器と大気蛍光望遠鏡とのハイ
ブリッド観測によって 15% 以下のエネルギー決定誤差を達成することを目標としている．
極高エネルギー宇宙線の粒子種が陽子であり GZK カットオフが存在するならば，一次エネル
ギーが 1020 eV 以上の宇宙線の発生源は 50 Mpc 以内に存在する．さらに，極高エネルギー宇宙線
は磁場による偏向をあまり受けない．すなわち，極高エネルギー宇宙線の発生源は近傍に存在し，宇宙線は到来方向の情報を失わない．しかし，現在極高エネルギー領域の宇宙線と候補天体現象との間に，有意な相関は得られていない．
極高エネルギー宇宙線の起源を特定するためには，到来する宇宙線のエネルギー，到来方向を高精度で決定する必要がある．地表検出器アレイは高精度で宇宙線の到来方向を測定することができるので，極高エネルギー宇宙線の到来方向と予想される発生源との相関を確認することができると期待される．
地表検出器アレイでは地表に到来する空気シャワーの粒子数から一次宇宙線のエネルギーを，到来時刻差から到来方向を決定している．これらの決定精度を決めるものの一つはシャワー構造の取り扱いである．
空気シャワーフロントは平面ではなく，曲率を持って地表に到来する．そのため，粒子の到来時刻差から到来方向を決定するためには，このシャワーフロント曲率の到来時刻への影響を考慮
しなければならない．また，その曲率は一次宇宙線のエネルギーや到来天頂角に依存する．シャワー面の曲率は平面からの時刻の遅れとして Linsley により実験的な式が作られ，さらに AGASA
グループにより改良が加えられている．しかし AGASA グループの求めた到来時刻の遅れの関数
は，空気シャワー中心から比較的近傍かつ，天頂角 45◦ までのシャワーのみを対象とするため，この関係式を使用するためには TA 実験に適した形に改良を加える必要がある．
本研究では，地表検出器アレイに対するモンテカルロシミュレーションと空気シャワーシミュ
レーションから，空気シャワー構造の詳細を調べ，粒子の到来時刻の遅れの関数と，その標準偏差が従う関数を作成した．これらの関数は天頂角 60◦ までのシャワーを対象としているため，統計量の増加が期待される．さらに，視野が広がることにより，宇宙のより広い領域で宇宙線の到来方向と候補天体現象との相関の探索が可能になる．また，その関数を用いて空気シャワーの到来方向を求める解析プログラムを作成し，さらに，シミュレーションイベントの解析から到来方
向決定精度を求めた．その結果，各エネルギーでそれぞれ 1.96◦ (1018.5 eVの時)，1.60◦ (1019.0 eV)，
1.13◦ (1019.5 eV)，0.86◦ (1020.0 eV) の精度で到来方向を再構成することが可能となった．

第1章

1.1

極高エネルギー宇宙線

宇宙線

宇宙線とは宇宙空間に存在する高エネルギー放射線であり，1912 年にオーストリアの物理学者
V.F.Hess によって発見された．その成分は約 90% が陽子で，約 9% が α 粒子，残りはそれより重い原子核である．これまでに観測された宇宙線のエネルギーは 108 eV から 1020 eV 以上までの広範囲にわたっており，その到来頻度はエネルギーのほぼ 3 乗に比例して小さくなることが知られ
ている．

1.1.1

宇宙線のエネルギースペクトル

現在までに観測されている一次宇宙線のエネルギースペクトルを図 1.1 に示す．

図 1.1: 一次宇宙線のエネルギースペクトル．右下の丸で囲まれた範囲は極高エネルギー宇宙線と呼ばれる領域である．

太陽活動の影響を受けない 1010 eV を越えるエネルギー領域では，宇宙線のエネルギースペク
トルはエネルギーの冪乗 F (E) = K × E −γ で近似される．この冪は 1015 eV 付近までは γ ∼ 2.7
であり，そこから γ ∼ 3.0 と変化し，1019 eV 付近からまた γ ∼ 2.7 となる．この二ヶ所のスペク
トルの折れ曲がりをそれぞれ “knee”，“ankle” と呼んでいる．このように，宇宙線の到来頻度は

1

エネルギーが大きくなるにつれて冪乗で小さくなるため観測量を増やすことが難しい．特に，極高エネルギー宇宙線と呼ばれる図 1.1 の丸で囲んだ 1020 eV 以上のエネルギー領域の宇宙線は，到来頻度が年間 100 km2 あたり 1 事象と非常に少なく，観測量が不足しているためにエネルギース
ペクトルがよく決まっていない．
knee より低いエネルギー領域の宇宙線は，超新星爆発の衝撃波で加速されると考えられている．
また，1015 eV 以上の宇宙線はこれ以下のエネルギー領域とは別の機構で宇宙線が生成，加速されているためスペクトルの冪が変化すると考えられている．
knee でのスペクトルの折れ曲がりの原因には全く別の考え方もある．それは，このエネルギー以上では荷電粒子のラーモア半径 rL が銀河円盤の厚さよりも大きくなるため，粒子が銀河円盤の外に逃げ出していく効果が現れる，というものである．銀河磁場は約 3 µG なので，1015 eV の陽子では rL ∼ 0.3 pc となる．この値は銀河円盤の厚さより小さいが，拡散運動の効果を合わせると，これ以上のエネルギーを持つ宇宙線に対しては銀河系内への閉じ込め効果が効きにくくなるため，エネルギースペクトルの冪が大きくなるとされている．
ankle 領域のエネルギーの宇宙線についても複数の理論がある．一つは knee と同様に宇宙線の曲率半径が銀河系の厚みより大きくなるため銀河内に閉じ込めることができず，主要な宇宙線の起源が銀河系内から銀河系外へと遷移するという考え方である．この場合は ankle 周辺での銀河系内起源の宇宙線の到来方向に異方性が検出されると期待される．また，化学組成は宇宙線の起源の変化に伴って鉄から陽子へと変化していると期待される．
一方，別の考え方では ankle 構造が p γ → p e+ e− の対生成によるエネルギー損失によって生じるという解釈もある．これが正しければ ankle 領域よりも低エネルギーで既に宇宙線の起源が銀河系外とならなければならない．そうすると化学組成は ankle 領域をまたいで陽子のままであると予想される．この場合，1017.5 eV 付近の “2nd knee” と呼ばれるエネルギー領域で宇宙線の起源は銀河系内から銀河系外に変化していると考えられ，それに伴って化学組成は鉄から陽子へ変化していると予測される．

1.1.2

宇宙線の化学組成

宇宙線の化学組成は，宇宙の化学組成と比べると，重い元素が H，He に比べて相対的に多い．
これは H，He が加速領域に入る際に比較的イオン化されにくいことと関係があると考えられている．しかし，大まかに見ると宇宙線の化学組成と，宇宙の化学組成はよく一致している．原子は星の中で核融合によって合成され，星風や爆発によって星間物質と星とを何代にもわたって経て宇宙の化学組成を形成している．つまり，宇宙線の化学組成が宇宙の化学組成から大きくずれていないのは，星の中で合成された原子核が何らかの形で加速されて宇宙線となっていることを示している．宇宙線の化学組成には，宇宙線の原子核がどこで合成され，どこで加速され，どのよ
うにして伝播してきたのかが反映されているはずである．従って，宇宙線中の放射性同位体比や重元素比を知ることで，宇宙線の起源について知ることができる．
knee より低いエネルギー領域の宇宙線の観測は人工衛星や気球，惑星探査機などの飛翔体技術の発達と測定技術の向上により Fe グループまでの元素に対しては精度よく行われている．その構
成比は数 10 MeV∼ 数 10 GeV まではほとんど変わらない．図 1.2 に He から Ni までの領域での，宇宙線の組成と太陽系の組成との比較を示す [1]．ここから，Li，Be，B と Sc，Ti，V，Cr，Mnの二つの元素群が，太陽系物質に比べて宇宙線には圧倒的に多いことがわかる．これらの元素は星の核合成の最終生成物としては存在しないが，炭素や酸素核の破砕 (Li，Be，B) や鉄の破砕 (Sc，Ti，Cr，Mn) によって生成されるためである．

2

図 1.2: Si を 100 としたときの，宇宙線の組成と太陽系の組成．• は低エネルギーの宇宙線，◦ は高エネルギーの宇宙線，⋄ は太陽系の組成を表している [1]．

3

1.1.3

宇宙線の加速機構

宇宙線の加速機構についても多くの議論がなされている．宇宙線の加速機構を考える上で，先に述べたような冪型のエネルギースペクトルを導き出せることが必須条件となる．現在，この条件を満たす加速機構として最も有力とされているのが “フェルミ加速機構” である．
フェルミ加速フェルミ加速は 1949 年に Fermi によって提唱された荷電粒子加速理論で，荷電粒子が磁場を持った星間雲との衝突を繰り返すことによって運動エネルギーを得るという統計的な加速モデルである．これによれば，星間雲と荷電粒子の衝突方向は統計的にランダムであるが，粒子の 1 回の衝突あたりの平均的なエネルギー増加量 ∆E は正であり，∆E = αE と粒子のエネルギーに比例することが示される．星間雲との衝突による加速は，この α が星間雲の移動速度の 2 乗に比例することから，二次のフェルミ加速と呼ばれる．実際の宇宙線加速モデルとしては，この二次フェルミ加速よりもプラズマ衝撃波による加速が有力であると考えられている．この場合，加速領域の形状の考察から α が衝撃波の移動速度の 1 乗に比例することが示され，一次のフェルミ加速と呼ばれる．荷電粒子が衝撃波との 1 回の衝突で ∆E = αE のエネルギーを獲得するとすると，n回の衝突後の粒子のエネルギー En は以下のようになる．

En = E0 (1 + α)n(1.1)

ここで，E0 は粒子の初期エネルギーである．よって，エネルギーが E になるのに必要な衝突回
数は，ln (E/E0 )(1.2)
n=ln (1 + α)となる．ここで，1 回の衝突で加速領域から逃げ出す確率を Pesc とおけば，n 回衝突後に加速領域に留まる確率は (1 − Pesc )n であるので，E 以上のエネルギーに加速される粒子の割合は，

N (≥ E) ∝

∞
X

m=n

(1 − Pesc )m =(1 − Pesc )n
Pesc

(1.3)

である．式 (1.2) を式 (1.3) に代入すると，N (> E) ∝となる．ここで，1 E −γPesc E0ln (1 − Pesc )−1
Pesc
γ=
≈
ln(1 + α)
α

(1.4)

(1.5)

である．このようにしてフェルミ加速では冪型のエネルギースペクトルが導き出される．
超新星残骸での衝撃波加速
超新星爆発残骸による加速での最大加速エネルギーを計算してみよう．
衝撃波の拡散係数を D，速度を u とし，拡散の特徴的な長さ D/u は衝撃波の曲率半径よりも十分小さく，平面と近似する．そのため，衝撃波による加速が行われる時間は放出する物質の総量から見積もることができる．例えば，質量 10 M⊙ ，平均速度が 5×108 cm/sec，放出する物質の平4均密度が 1 proton/cm3 の場合，拡散の特徴的な時間は TA ∼1000 yrs となり，広がりの大きさは
5 pc 程度となる．
このように，超新星爆発からの強い衝撃波が有限な寿命 TA を持つことによって，この加速機構
による 1 粒子あたりの加速エネルギーが制限される．単位時間あたりのエネルギーの増加率は，

dE
αE
=
dt
Tcycle(1.6)

ここで，Tcycle は 1 回の衝突にかかる時間，α は衝突ごとに得られるエネルギーの獲得率である．
また，α は以下の式で表される．
4 u1 − u2
α=
(1.7)
3
c
ここで，u1 は上流領域での流体の速度，u2 は下流領域での流体の速度である．
式 (1.6) を積分することで最大加速エネルギー Emax を見積る．
Z Tcycle
dE
· dt
(1.8)
Emax =
dt0
ここではまず，Tcycle を見積もる．衝撃波の上流での粒子数の流れは，

J = −D∇N + uN

(1.9)

より与えられる．u1 は衝撃波面に対して相対的に負の方向になるため，dN= −u1 Ndz

(1.10)

N (z) = ρcr exp (−zu1 /D1 )(1.11)

D1となる．それゆえ上流領域では，となる．ここで，D1 は衝撃波の拡散係数，ρcr は衝撃波での宇宙線の数密度である．上流領域の単位面積あたりの総粒子数は，ρcr D1 /u1 となる．またショック面と宇宙線との単位面積あたりの衝突率は ρcr c/4 である．そのため上流領域での滞在時間は(ρcr D1 /u1 ) (ρcr c/4)−1 = 4D1 / (u1 c)(1.12)となる．下流領域では滞在時間の計算に粒子が脱出するかどうかを考慮にいれる必要があるので，
簡単ではない．ここでは単純に上流領域と同じであるとして，滞在時間を 4D2 / (c u2 ) とすると [2]，
衝突にかかる時間は，
4 D1 D2
Tcycle =+(1.13)c u1
u2となる．ここで，D2 は下流領域での衝撃波の拡散係数である．
次に拡散係数の見積もりが必要である．拡散の平均自由行程 λD は，ラーモア半径よりは小さくできない．ラーモア半径は荷電粒子の磁場中でのらせん運動の半径であり，次の式 (1.14) で与えられる．
pc(1.14)rL =ZeB
ここで，e，Z ，p，B はそれぞれ，電気素量，宇宙線の原子番号，磁場に垂直な方向の運動量，磁場の強さである．従って拡散係数の最小値 Dmin はこのラーモア半径を使って，Dmin =

rL c
1 Ec
∼
3
3 ZeB
5

(1.15)

となる．ここでは，Tcycle ≥ 20E/(3u1 ZeB) となるような強いショック u2 = u1 /4 を考えている．
D1 = D2 = Dmin を代入すると最大加速エネルギー Emax はEmax ≃3 u1
ZeB (u1 TA )20 c

(1.16)となる．質量 10 M⊙ の天体の超新星爆発の放出速度 5×108 cm/sec, TA ∼ 1000 yrs，銀河磁場を
3 µG とすると，Emax ≃ Z × 3 × 1013 [eV](1.17)となる．
ボトムアップモデル
フェルミ加速あるいはその他のメカニズムによって，低いエネルギーの宇宙線をより高いエネ
ルギーへと加速し，極高エネルギー宇宙線を放射している放射源モデルを総称してボトムアップ
モデルと呼ぶ．
荷電粒子の衝撃波加速を考えると，荷電粒子の進行方向を曲げるのは磁場であるから，加速領域への閉じ込め条件から加速されるエネルギーの限界が決まり，

Emax ≤ γeZBR(1.18)となる (前節と異なり衝撃波の寿命は十分長いものとする)．ここで，γ は衝撃波のローレンツファクター，Z は宇宙線原子核の原子番号，B は加速領域の磁場の強さ，R は加速領域の大きさである．図 1.3 に，宇宙線源候補天体の大きさと磁場強度，そして最大加速エネルギーの関係を示し
た [3]．この図から，銀河系そのものを含むほとんどの銀河系内天体は極高エネルギー宇宙線の加速源候補となりえない (ガンマ線バースト (GRB)，中性子星が例外) ことがわかる．銀河系外の候
補天体は活動銀河核，ガンマ線バースト，電波銀河，銀河団などが挙げられる．
以下，いくつかの候補天体での加速と観測される特徴について簡単にまとめる．

• 中性子星：表面の磁場が 1013 G を越えるような中性子星は，鉄原子核を 1020 eV 以上まで加速することが可能である [4]．しかしこのような非常に強い磁場が存在すると，荷電粒子はシンクロトロン放射によってエネルギーを失い，実際にはそのような高エネルギーまでの加速は不可能であると考えられている．ただし，シンクロトロン放射によるエネルギー損失率
は中性子星表面からの距離に強く依存し，表面から 109 cm 以上離れた領域で加速された粒子は大きなエネルギー損失を伴わず外部へ脱出できる，という指摘もある [5]．
このモデルが正しければ，極高エネルギー宇宙線の到来方向分布は銀河面に集中し，鉄原子核が主成分となるはずである．

• AGN：活動銀河核 (Active Galactic Nuclei，AGN) では，中心にあるブラックホールに降着する物質の重力エネルギーを粒子の運動エネルギーとして解放している [6]．粒子加速は回転する磁場によって生ずる電場によって起こる．ただし，AGN での加速では，その周囲の放射場との衝突による荷電粒子のエネルギー損失が顕著であり，最大加速エネルギーは1019 eV までと考えられている [7]．
• 電波ローブ：AGN の電波ローブと呼ばれる領域では，中心核から放出されるジェットと銀河間空間ガスの衝突によって強い衝撃波が生じており，100 kpc 以上にも広がっている [6]．
6Magnetic Field Strength
12

10 G

Neutron Stars

6

10 G
LHC
TEVATRON
SppS

GRBs
20 e

10

V
at
β

at

1G

20 eV

10

White Dwarfs
AGNs
Sunspots

=

β

1/

=

30

1

12

10

0

eV
Radio Galaxy
SNRs

-6

10 G

Interplanetary space
Galactic Cluster
Galactic Halo

1km

6

10 km
1AU

1pc

1kpc

1Mpc 1Gpc

Size

図 1.3: 磁場の強さと天体の大きさと最大加速エネルギーの関係 (Hillas プロット)[3]．

その中で特に強い電磁波が放射されている領域 (radio hot-spot) では，光子，物質ともに密度が低いためエネルギー損失が少なく，かつ粒子を封じ込めるだけの十分な強さの乱流磁場が存在していることが電波観測から予想されている [8]．しかし，hot-spot を伴う電波ローブを持つ AGN はまれで，地球からの距離も遠く，しかもこれまでに観測された極高エネル
ギー宇宙線の到来方向分布と合わない．仮に M87(地球からの距離：約 18 Mpc) がそのような放射源であれば，M87 の方向に観測事象が集中するはずである．また，次に近い既知の類似な天体は NGC315 で，これは 100 Mpc と GZK 限界 (1.2.1 節詳述) を越えており遠すぎる．

• GRB：ガンマ線バースト (Gamma Ray Bursts，GRB) は，既知の天体現象の中でも最も激烈なものの一つであり，極高エネルギー宇宙線の到来方向と同様に等方的に分布している．
GRB によるエネルギー放出率は，極高エネルギー宇宙線源が一様に分布していると考えた
場合の 1019 eV 以上の宇宙線としてのエネルギー放出率にほぼ等しく，非常に有力な候補天体である [9]．しかし，GZK 機構のために地球近傍 50 Mpc 以内に存在しなければならない．
このような GRB は 100 年に 1 回程度しか発生しないと予想されているので，観測される宇宙線の到来方向分布は GRB の方向に集中することが容易に予想できる．
• 銀河団：数 µG オーダーの磁場が 500 kpc 程度の広がりを持って観測されており，これは最大加速の条件 (式 1.18) を満たしている [10]．しかし，1020 eV 以上まで宇宙線が加速されたとしても，銀河団内部での GZK 機構によるエネルギー損失が大きく，銀河団の外には放射されないと考えられている．

7

トップダウンモデル
前述のように，これまでに観測された極高エネルギー宇宙線を既知の天体現象や素粒子物理，宇宙論の標準的な枠組の中で説明するのではなく，素粒子物理の標準的な理論を越えた新しい物理に基づくモデルや，未知の超高エネルギー天体現象に基づくモデルが提案されている．こうしたモデルは，従来のような低エネルギー粒子が加速されて極高エネルギーに達するというボトムアップモデルに対して，トップダウンモデルと呼ばれる．
以下はそのいくつかの例である．

• 宇宙紐や超重粒子の崩壊・対消滅：宇宙の初期，すなわちインフレーションの最後または直後に生成された宇宙紐などの位相欠陥 (Topological Defect)，または超重粒子 (Super Heavy
Relic Particles) が，銀河近傍で崩壊して極高エネルギー宇宙線を生成しているとするモデルである [11]．位相欠陥は互いの衝突や対消滅によって，超重ゲージボソンやヒッグスボソンに崩壊すると考えられている．これらはパイオンのハドロンジェットを作り最終的にガンマ線，電子，ニュートリノに崩壊し，ハドロンのうち数 % だけが核子となる．GZK 機構を考慮すると，このシナリオでは，GZK カットオフより低いエネルギーでは陽子が，高いエネ
ルギーではガンマ線が主成分になると予想される．超重粒子も崩壊によってハドロンジェッ
トを作るが，超重粒子は cold dark matter として銀河ハロー内にも存在すると考えられる
ため，GZK 機構の影響を受けず，粒子種がエネルギーとともに変化することはないと予想されている．
• Z バースト模型：極高エネルギーのニュートリノが，重力によって銀河近傍に集積している
宇宙背景ニュートリノと相互作用して Z0 粒子を作り，その崩壊生成物が極高エネルギー宇
宙線として観測される，とするモデルである [12]．この反応で生成されるのはほとんど光子である．また，宇宙背景ニュートリノは銀河近傍に大きな半径を持つクラスターを作るので，到来方向分布は等方的になると予想される．
もしトップダウンモデルが正しければ，標準的な素粒子物理と天体物理の理論を越える新たな物理への糸口となる可能性がある．

1.2

極高エネルギー宇宙線

宇宙線物理学という研究分野の中で今，最も興味深いトピックスの一つが「極高エネルギー宇
宙線」である．極高エネルギー宇宙線とは，図 1.1 の円で示された 1020 eV 付近，あるいはそれ以上のエネルギーを持つ宇宙線である．この節では極高エネルギー宇宙線の伝播において重要ない
くつかの項目を説明する．

1.2.1

GZK カットオフ

1963 年に宇宙背景放射 (Cosmic Microwave Background Radiation，CMBR) の存在が発見された直後，Greisen，Zatsepin と Kuzmin は，極高エネルギー宇宙線は伝播中にこの CMBR と衝突して大きなエネルギー損失を起こす，と指摘した [13][14]．2.7 K の CMBR(∼ 10−3 eV) は，極高エネルギー宇宙線陽子の静止系では 100 MeV を越える高エネルギーのガンマ線と同等であり，

8

陽子は高い確率で核子共鳴状態に励起される．共鳴状態は以下のような反応で π 中間子を放出して崩壊するので，大きなエネルギー損失が起こる．

γ(∼ 150MeV) + p(静止) → ∆(1232) → p + π 0 , n + π + , · · ·(1.19)

また，この共鳴が起こるエネルギーよりも低いエネルギーで，電子対生成 (γ p → e+ e− p) によるエネルギー損失も起こる．

図 1.4: 極高エネルギー宇宙線の減衰長．細い実線が光子，太い実線が陽子，点線が鉄の減衰距離である．

図 1.4 に銀河間空間における光子数分布モデルにもとづいて推定された宇宙線の減衰長を示す．
これによれば，1020 eV 以上の陽子の減衰長は 50 Mpc 以下であることがわかる．従って，地球で観測された 1020 eV 以上の宇宙線は，その起源が地球から 50 Mpc 以内になければならない．逆にそれより遠方の源から放射された宇宙線は，どんなに高い初期エネルギーを持っていても，伝播
中にエネルギーを失い，地球では 1020 eV 以下の宇宙線として観測されることになる．
宇宙線が陽子よりも重い原子核である場合は，それ以外のエネルギー損失も起きる．2×1019 eV以上の宇宙線原子核は CMBR 光子との衝突による photo-disintegration によって 1 Mpc あたり
3∼4 個の核子を失う．従って，観測された極高エネルギー宇宙線が陽子より重い原子核である場
合は，その源は 20 Mpc 以内に存在していなければならない．
以上のような反応により極高エネルギー宇宙線は宇宙に広く分布する CMBR によって遮蔽されるため，地球で観測されるエネルギースペクトルには 1020 eV 付近に鋭いカットオフが現れると考えられている．特に一次陽子の場合 GZK カットオフと呼び，このカットオフが現れるエネルギー
限界を GZK 限界と呼ぶ．図 1.5 に様々な宇宙線源の分布モデルから予想されるエネルギースペク
トルを示した．例えば，源が宇宙空間に一様に分布する場合，5×1019 eV に急激なカットオフが9

現れると予想されている．したがって，観測される極高エネルギー宇宙線のスペクトルが実際にどのような構造を持つかを明らかにすることは，宇宙線源を研究する上で重要な鍵となる．

図 1.5: 宇宙線のエネルギースペクトルに予想される GZK カットオフ [13][14]．一つの宇宙線源がある red shift(z = 0.004 ∼ 1.0) の位置にあると仮定した場合に予想されるスペクトル．右は宇宙線源が空間分布しているとしてその分布がいくつかの異なる evolution model に従うと仮定した場合に予想されるスペクトル．

1.2.2

磁場による偏向

宇宙線は銀河系磁場や銀河系外磁場による偏向を受けることで，到来方向の情報を失ってしまう．すなわち，発生源を同じとする宇宙線であっても，磁場による大きな偏向を受けると，我々には異なる方向から到来したように観測される [15]．磁場による偏向を理解することは，宇宙線と既知の天体現象との相関を調べる上で非常に重要である．
銀河磁場
銀河磁場の大きさは，ファラデー効果により測定されている．

RM ∝

Z d
0

ne Bk · dl(rad/m2 )(1.20)

ここで，RM は偏向面の回転の大きさ，d はソースまでの距離，ne は電子数密度である．
また，銀河磁場の大きさは，銀河中の物質密度に比例し，銀河中心から離れるほど減少する．銀
河磁場と銀河中心からの距離との関係は以下の式のようになる．

B ∝ eRGC (kpc)/8.5(1.21)

ここで，RGC は銀河中心からの距離である．しかし，銀河中心近傍の磁場の研究はまだあまり進んでおらず，この関係式を使用できるのは RGC ≥ 3 kpc の領域のみである．これより，太陽系の磁場の大きさは約 2 µG である．

10

これらの通常磁場に加えてランダム磁場も存在し，太陽系では通常磁場の 2 倍の約 4 µG である．しかし，ランダム磁場の大きさのスケールは 50 ∼ 100 pc 程度であり，この磁場による極高エネルギー宇宙線への偏向の影響は重要ではない．
また，非常に重要であるが未解決の問題として，銀河ハローの磁場の取扱いがある．近年，銀河ハローの磁場が極高エネルギー宇宙線の散乱に大きな影響を与えるとする研究も報告されている [16]．
これらの磁場の影響で，1019 eV の宇宙線は銀河中で 2◦ ∼ 4◦ 程度の偏向を受けると見積もられている [15]．
銀河磁場下の極高エネルギー宇宙線陽子の軌跡を図 1.6 に示す [17]．
銀河間磁場
銀河間磁場の大きさは宇宙空間中の銀河の密度より，上限が与えられる [18]．µG 程度の磁場を持つ銀河が，1 Mpc の距離で分布していると仮定すると，平均銀河間磁場の上限は 1 nG である．
ランダムウォークによる偏向角 θ は次式で表される．
1/2−1
B−9 d100 LC
θ = 2.5◦ E20(1.22)

ここで，E20 は 1020 eV を基準としたエネルギー，B−9 は 1 nG を基準とした磁場の大きさ，d100
は 100 Mpc を基準としたソースまでの距離，LC は銀河間の平均距離 (1 Mpc) である．
また，実際の宇宙線の伝播距離 ∆d はソースまでの距離よりも大きく，以下の式で表される．
−2 2∆d = 0.047E20
B−9 d2100 LC [Mpc](1.23)

銀河間磁場下の極高エネルギー宇宙線陽子の軌跡を図 1.7 に示す [19]．
図 1.6，図 1.7 からわかるように，1020 eV を超えるエネルギーの宇宙線は銀河内外ともにほとんど偏向されない (銀河間磁場下で 1020 eV の宇宙線は 50 Mpc 伝播する間に約 2◦ しか偏向されない [20])．
以上より，1020 eV を超えるエネルギーの宇宙線の源は比較的近傍に存在し，宇宙線は到来方向の情報を失わない．しかし，現在このエネルギー領域の宇宙線と候補天体現象との間に，有意な相関は得られていない．

図 1.6: 銀河磁場下の極高エネルギー宇宙線陽子の軌道を表す [17]．

11図 1.7: 銀河間磁場下の極高エネルギー宇宙線陽子の軌道を表す [19]．

1.3

空気シャワーと宇宙線観測

1.1 節で述べたように，宇宙線のフラックスはほぼ E −3 で減少するので，109 eV∼1020 eV のエ
ネルギー領域ではフラックスの変化は 1030 にもなる．フラックスの大きい 1014 eV 以下の低エネ
ルギー宇宙線は，気球や人工衛星などの飛翔体に搭載された観測装置によって直接観測される．しかし 1014 eV 以上の高エネルギー宇宙線は到来頻度が少ないため，観測には大きな検出面積と長い露光時間が必要であり，そのため飛翔体搭載検出器によって直接観測をすることは難しい．そこで，高エネルギー宇宙線に対しては一次宇宙線が大気と相互作用してできる二次宇宙線 (空気シャワー) を検出するという間接測定が行われている．

1.3.1

空気シャワー現象

宇宙線が大気に入射すると，大気中の原子核と相互作用して二次粒子を生成し，さらにそれらの二次粒子も粒子生成を行なう．この過程をくり返し，一次宇宙線が結果的に大量の二次粒子群となる現象を空気シャワー現象と呼ぶ．通常空気シャワーは高エネルギーハドロンによる核カス
ケードと，高エネルギーガンマ線から生成される電磁カスケードからなる．
宇宙線の主成分である陽子などのハドロンが大気原子核と相互作用した場合には，π ，K などのハドロンの多重生成 (核カスケード) が起きる．生成された二次粒子の大部分は結果的に π となる．そのうち π 0 は短い寿命 (8.4 × 10−17 s) で 2 個のガンマ線に崩壊し，電磁カスケードを形成す
る．π ± は τ = 2.60 × 10−8 s で，π ± → µ± + νµ(1.24)のように崩壊する．µ± は τ = 2.20 × 10−6 s で，µ+ → e+ + ν¯µ + νeµ− → e− + νµ + ν¯e

12(1.25)(1.26)

のように崩壊する．大気の厚みは核相互作用の平均自由行程の約 10 倍であるため，相互作用が繰り返されることで多数の粒子が生成される．これを核カスケードと呼ぶ．図 1.8 はこの様子を模式的に表したものである．
高エネルギーのガンマ線は，電子対生成により電子と陽電子を作る．これらの電子，陽電子は制動放射によってガンマ線を放射する．エネルギー E の電子による制動放射の断面積 σbrems. は，次の Bethe-Heitler の式で表される．
 


 1
4Z 2 re2 dv
22− 13
σbrems. (E, v)dv =+ (1 − v)(1.27)
1 + (1 − v) − (1 − v) ln 184Z
137 v
3
9

ここで，v は hν/E ，ν は制動放射される光子の振動数, Z はターゲットとなる原子核の原子番号，reは古典電子半径である．従って，制動放射による単位厚さ当たりのエネルギー損失 (dE/dX)brems.
は以下のように表される．


E
dE
≃−
(1.28)
dX brems.
X0

ここで，X0 は電子のエネルギーが 1/e になる厚さで輻射長と呼ばれ，以下のように定義される．

1
1
4Z 2 re2 N (1.29)=
ln 184Z − 3
X0
137 A

ここで，N はアボガドロ数，A はターゲット原子核の質量数である．X0 は大気ではおよそ 38 g/cm2
である．一方，電子対生成の断面積 σpp は，

 

 1
2
4Z 2 re222− 13
du u + (1 − v) − (1 − v) ln 184Z
(1.30)σpp (hν, u)du =+ (1 − u)137
39となる．ここで，u = E/hν, E は生成される電子のエネルギーをあらわす．この断面積から電子対生成の interaction length が求められ，結果として，


dE
7 E
=−
(1.31)
dX pair
9 X0となる．従って，電子対生成と制動放射の interaction length は同程度であることがわかる．この二つの過程を繰り返して，多数の電子，陽電子，ガンマ線がつくられる現象が電磁カスケードである．粒子数が増大し 1 粒子あたりのエネルギーが減少すると，大気中での電離損失が優勢とな
り，電子成分が大気に吸収されて粒子総数が減少する．この境目のエネルギーを臨界エネルギーと呼び，空気中では 84 MeV である．

1.3.2

空気シャワーの縦方向発達

空気シャワー中に含まれる粒子の総数が増大するにつれて，個々の粒子の持つエネルギーは減少する．空気シャワー中の電子成分の場合，臨界エネルギー Ec (空気中では Ec =84 MeV) に近付
くと大気中の原子，分子による電離損失が大きくなり，急速にエネルギーを失って大気に吸収されてしまう．そのため，空気シャワー中の粒子の総数はある段階で減少に転じる．このような空気シャワーの通過した物質量に対する粒子数の変化は縦方向発達 (longitudinal development) と呼ばれる．
エネルギー E0 の 1 個のガンマ線が大気に入射したときに生じる空気シャワーの縦方向発達については，次の近似式が用いられる．
 

0.31
3
Ne (t) ∼ √ exp t 1 − ln s
(1.32)y213図 1.8: 空気シャワーの模式図．

y = ln
s=



E0
Ec



(1.33)

3t
t + 2y(1.34)

ここで，t は電子の大気中での輻射長 38 g/cm2 を単位とするシャワーが通過した大気厚さである．
s はシャワーの発達段階を表すパラメータでエイジパラメータと呼ばれる．このシャワーエイジはシャワーの発達とともにその値は増加し，Ne (全電子・陽電子数) が最大になるとき，s ∼ 1 となる [21]．
宇宙線陽子がつくる空気シャワーについては，Gaisser と Hillas の提案に基づく以下の式が用いられることが多い．

Ne (X) = Nmax

X − X0
Xmax − X0 Xmax −X0
λ

exp



Xmax − Xλ(1.35)

ここで，Nmax は最大発達時の粒子数，Xmax はシャワー最大発達深さ [g/cm2 ]，X0 は一次宇宙線の最初の相互作用点の深さ [g/cm2 ]，λ は縦方向発達の attenuation length で 70 g/cm2 である．
一次宇宙線のエネルギーの違いは Nmax の違いとなって現れ，およそ Nmax ∼ 2 × E0 × 10−9 (E0は一次宇宙線のエネルギー，単位：eV) の関係がある．また，同一の一次エネルギーである場合，一次宇宙線粒子種の違いは X0 ，Xmax に現れ，重い原子核ほど X0 , Xmax の値が小さくなる傾向がある．これは，重い原子核ほど大気との相互作用断面積が大きいため，より上空で発達が始まること，さらに核子当たりのエネルギーが小さいため，エネルギー細分化が早く進むことによる．
簡単なモデルでは，Xmax ∝ ln(E0 /A)(A は一次宇宙線の質量数) である．陽子，鉄原子核による空気シャワーをシミュレーションして求めた縦方向発達曲線を図 1.9 に示す．

14

図 1.9: 空気シャワーの縦方向発達．

1.3.3

空気シャワーの横方向分布

電磁カスケードは大気中で電子の多重クーロン散乱により横方向に広がっていく．電磁カスケードの横方向分布は解析的に求められ [22]，Greisen により整理された NKG function が一般的に使われる [23]．
空気シャワー中心からの距離 R の位置での電子数密度 ρe (R) は次式で表される．

ρe (R) = C

Ne
2
RM



R
RM

s−2.0 
R s−4.5
1+RM(1.36)

ここで，Ne は総荷電粒子数，C は規格化定数，RM はモリエールユニットである．モリエールユ
ニットは，電子が X0 進む間にクーロン散乱により最初の方向から垂直方向にずれた厚さであり，
RM ∼ 9.3 g/cm3 である．また s はエイジパラメーターであり，電磁カスケードの発達の度合いを示すものであり，
3ts=(1.37)
t + 2ln (E0 /ǫ0 )である．ここで，E0 はカスケードの親粒子のエネルギー，t は X0 を単位とした大気の厚さである．
ハドロンにより引き起こされる空気シャワーは純粋な電磁カスケードとは少し様相が異なる．ハ
ドロンにより引き起こされる空気シャワーの電磁成分は，核相互作用によって生成された多世代の π0 の崩壊からの電磁カスケードの重ね合わせである．シャワーの中心部は核カスケードが発達する過程で次々に生成される π0 の崩壊による発達途上の電磁カスケードが主成分であるが，シャワー軸からはなれたところでは核カスケードの初期段階で生成された π0 の崩壊からの粒子の増殖が終わった電磁カスケードが主成分となる．このため，シャワー中心から遠く離れた場所での横方向分布は，NKG function の外挿により得られる分布よりフラットになることが予想されてお15

り，このシャワー中心から遠方での横方向分布を考慮した式が Linsley により与えられている [24]．
シャワー中心からの距離 R における電子数密度 ρe (R) は

 

R α
R −(η−α)ρe (R) ∝1+(1.38)
RM
RMと表される．AGASA グループの結果によると，電子成分の横方向分布は 1 km 以遠で粒子数密度が減少しているので，上記の Linsley による 1.38 ではその構造を十分に表すことができず，横方
向分布に以下のように補正項を加える [24]．

2 !−δR −(η−1.2)
R
R −1.21+
1.0 +
ρe (R) = C(1.39)
RM
RM
1[km]
次に高エネルギーミューオンは，上空でのハドロン相互作用の情報を持っているが，その多くはコア近傍 (∼ 数 10 m) に集中している．しかし一般に観測有効面積を大きくとる必要のある大規
模な空気シャワーアレイによる観測では高エネルギーミューオンの観測は難しく，遠方まで広がる低エネルギーミューオン (∼ 1 GeV 領域) が観測対象となる．その横方向分布は Greisen による実験式で表される [25]．
  Nµ Cµ R −αR −δρµ (≥ Eµ , R) = g(Eµ , R)
(1.40)
1.0 +
R0
R0
R02

Cµ =

Γ(β)
2πΓ(2 − α)Γ(α + β − 2)

(1.41)
0.17g (Eµ , R) = [51/(50 + Eµ )][3/(2 + Eµ )]0.14R(1.42)

α = 0.75, β = 2.5(1.43)

ここで，Eµ はミューオンのエネルギー閾値，Nµ は総ミューオン数，R0 はミューオンの平均的な広がりを示す距離である．
いくつかのエイジにおける横方向分布を図 1.10 に示した．

1.3.4

空気シャワー粒子の到来時刻分布

図 1.11 に空気シャワーの概略図を示す．地表へのシャワー粒子の到来時刻はシャワー軸からの距離に応じてシャワー平面 (シャワー軸に垂直でシャワー軸上のシャワー粒子を含む仮想平面) に対して時間差が生じる．また，シャワーフロントは厚みを持っているため，その厚みも考慮する必要がある．粒子の平均到来時刻分布 hσt i は Linsley により実験的に与えられており，シャワー
中心からの距離を R，空気シャワーの天頂角を θ，一次粒子のエネルギーを E として次式で表される [26]．
R bhσt i = σt0 1 +[sec](1.44)
Rt

ここで，σt0 =1.6 ns，Rt =30 m，b = (2.08 ± 0.08) − (0.40 ± 0.06)secθ + (0 ± 0.06)log(E/1017 eV)である．
この関数は AGASA グループにより改良されており，シャワー平面に対するシャワー粒子の遅れの平均 TD と，その標準偏差 TS は
1.5R
TD (ρ, R) = 2.6 1 +ρ−0.5 × 10−9 [sec](1.45)
3000[cm]16

図 1.10: 空気シャワーの横方向分布．


TS (ρ, R) = 2.6 1 +R
3000[cm]

1.5

ρ−0.3 × 10−9 [sec](1.46)

で表される [27]．ここで，ρ は粒子数密度 [m−1 ] である．

1.3.5

空気シャワーの観測方法

現在の空気シャワーの観測方法は大きく分けて，地表粒子検出器による観測法と大気蛍光法の2 種類がある．
地表粒子検出器による観測法
多数の粒子検出器を地表に配置して空気シャワー粒子を検出し，宇宙線を間接観測する装置は空気シャワーアレイと呼ばれ，高エネルギー宇宙線の観測装置として伝統的に用いられてきた．
空気シャワー粒子は，約 90% の電磁成分 (電子，陽電子，光子) と約 10% のミューオン成分，さらに少量のハドロン成分からなる．このうち電磁成分とミューオン成分を観測する検出器として最も一般的に用いられているのが，プラスチックシンチレータを用いたシンチレーション検出器である．また，特にミューオン成分のみを選択的に検出するために，シールドされたシンチレー
ション検出器や比例計数管でミューオン検出器を構成し，別に配置する場合もある．さらに，空気シャワー軸に近いハドロン成分を観測するために，ハドロンカロリメーターを備えている空気
シャワーアレイもある．
空気シャワーアレイの場合，一次宇宙線の到来方向は比較的簡単に求められる．空気シャワーは図 1.11 のようにシャワー面は粒子群 (シャワー平面) を形成して地表に到来する．この厚みを考17慮しつつ各検出器への到来時刻差から空気シャワーの到来方向を求めることができる．一方，一次宇宙線のエネルギーは，各検出器で検出された粒子数密度から空気シャワー全体での横方向分
布と全粒子数を推定することで得られる．しかし，横方向分布からエネルギーを推定するためには，ハドロン相互作用を考慮した大規模なモンテカルロ計算を必要とする．その際のハドロン相
互作用の不確定性の影響を受けることは避けられない．また，空気シャワー発達のゆらぎの影響も無視できない．

図 1.11: 地表付近の空気シャワー断面の概念図．

大気蛍光法
大気蛍光法は，空気シャワー粒子の電離損失の一部が光に変換されることを利用するものである．高エネルギーの荷電粒子が大気中を通過すると，その軌跡の近傍にある窒素などの大気分子は励起され，これらの励起された大気分子から蛍光が発せられる．すなわち，大気は宇宙線に対してシンチレーターとして働く．80 MeV の電子が 1 気圧の大気中を 1 m 通過した場合に発せられる大気蛍光光子数はおよそ 4 個であるが，膨大なシャワー粒子数のためにシャワー全体での光子数は 2 × 1015 個 (一次エネルギーが 1020 eV の場合) にも達し，その発光の継続時間は 30 µs 程度になる．これを光学望遠鏡で集光し高感度で高速の光センサーカメラで撮像すれば，およそ 30 km
遠方からでも観測することができる．さらに大気蛍光の放射方向は等方的であるため一次宇宙線の到来方向によらず観測を行うことができるという利点がある．このような大気蛍光法は，新たな空気シャワー観測の方法として 1960 年代初めに菅，Greisen，Chudakov などによって提案された．
大気蛍光法による空気シャワー観測結果からの一次宇宙線の各種情報の推定は，大きく分けて

1. 空気シャワー軸の再構成．
2. 空気シャワーの縦方向発達の再構成．
の 2 段階を経て行われる．宇宙線の到来方向 (=シャワー軸) は蛍光の軌跡から簡単に求めることができる．特に，遠く離れた 2 台の望遠鏡でシャワーをステレオ観測すると，それぞれの望遠鏡

18

が決めたシャワー平面の交差した線がシャワー軸となる (図 1.12)．空気シャワー軸を決定した後，さらに次のような手順で空気シャワーの縦方向発達の再構成をする．

i. 空気シャワー中の電子 1 個から光電子増倍管に到達する光子数の期待値を，その電子の空間的な位置，望遠鏡までの距離から大気による散乱・吸収，装置の検出効率などを考慮して計算する．この光子数には大気蛍光だけでなく，チェレンコフ光も考慮する．
ii. 実際に光電子増倍管で検出された光電子数と i．で計算される値を比較し，空気シャワーの各発達点での電子数を見積もり，式 (1.35) と比較することでシャワー発達関数 Ne (X) を求める．
空気シャワーが地表に到達するまでに失うエネルギーの 90% 以上は，シャワー中の電子成分の電離損失によるものであるので，観測されたシャワー発達関数 Ne (X) から式 (1.47) により一次宇宙線のエネルギー E0 が決まる．
Z
ǫ0
E0 =
Ne (X) dX(1.47)Xrad
ここで，Xrad は大気中での輻射長で 38.0 g/cm2 ，ǫ0 は大気に対する臨界エネルギーで 84 MeV である．また，シャワーの縦方向発達の様子 (特に最大発達深さ Xmax ) から，一次粒子の質量組成の測定が可能である．
このように，大気蛍光法による観測では，シャワーの縦方向発達が推定され，一次エネルギーはモンテカルロ計算に依存することなく決定される．ただし，この方法でエネルギーを測定するには，大気中の大気蛍光発光効率，大気中での紫外光の散乱・吸収率，鏡の反射率，フィルターの吸収率，光電子増倍管の量子効率などをすべて積み上げて掛け合わせていく必要があり，系統誤差を正しく管理することは容易ではない．
大気中での発光過程と散乱・吸収において考慮する必要のある系統誤差として，チェレンコフ
光，レイリー散乱，ミー散乱がある．チェレンコフ光は大気蛍光のバックグラウンド光となり系
統誤差の原因となるため，データ解析の際には事象ごとに混入するチェレンコフ光を見積もる必要がある．レイリー散乱は大気分子による散乱で，波長より十分小さな直径の誘電体球による光の散乱である．それに対し，波長に対し無視できない大きさの誘電体球による散乱をミー散乱と呼び，主に大気中のダスト，エアロゾル (霧，もや，塵，煙など) によって起こる．レイリー散乱
による大気の透過率は大気分子の密度から簡単な計算で求められるが，ミー散乱による透過率は
大気中のダストやエアロゾル分布の様子に依存するため，大気モニターによる実測から測定する必要がある．

1.4

観測結果の現状

極高エネルギー宇宙線は宇宙空間での爆発的な天体現象や未知の物理と結びついていると考えられ，観測結果が注目されている．この極高エネルギー宇宙線を，広大な領域の地表検出器アレ
イと大気蛍光望遠鏡で同時観測する実験として，Pierre Auger 観測所と宇宙線望遠鏡実験が現在稼働している．両者の実験結果にはある程度の一致が見られるものの，到来頻度の少なさと両者の系統的なズレから，極高エネルギー宇宙線の起源についての統一的見解はまだ得られていない．

19

図 1.12: 大気蛍光による空気シャワーのステレオ観測の概念図．

Pierre Auger 観測所
Pierre Auger 観測所 (以下 Auger) は南米アルゼンチンのパンパに展開されている．大型水タン
ク (10 m2 ×1.2 m) を一辺 1.5 km の三角形の網目状に 1600 台並べ，3000 km2 をカバーし，その周り四ヶ所に望遠鏡を設置し大気蛍光を観測する．大型水タンクでは粒子の通過距離が長いため，シャワー中の電磁成分よりもミューオンに対しての感度が高い．また 1.2 m の高さを持つためニュートリノから生じる水平方向からのシャワーに対する感度がプラスチックシンチレーター検出器に比べ高い．大気蛍光望遠鏡は視野角 29◦ を解像度 1.5◦ でカバーしている．
地表検出器でのエネルギー決定にはシャワー中心から距離 1000 m でのミューオン粒子数を用い，これを大気蛍光望遠鏡と同時観測した事象を使って較正している．エネルギー決定精度に約 22%の系統誤差を持つ [28]．
宇宙線望遠鏡実験
宇宙線望遠鏡実験 (以下 TA) は米国ユタ州の砂漠に展開されている．3 m2 ，2 層のプラスチック
シンチレーター 507 台を格子状に並べ，約 700 km2 をカバーし，その周り三ヶ所から仰角 3 ∼ 33◦のアレイ上空を解像度 1.0◦ で観測している．また，40 MeV 電子線形加速器を観測サイトに設置
し，大気中に射出される電子ビームを基準光源として，望遠鏡の一括較正を行う試みが進行中である．
地表検出器でのエネルギー決定にはシャワー中心から距離 800 m での粒子数密度を用い，これを大気蛍光望遠鏡と同時観測した事象を使って較正する．エネルギー決定精度に約 21% の系統誤差を持つ [28]．

1.4.1

エネルギースペクトル

図 1.13 に各実験で報告されているエネルギースペクトルを示す．また，表 1.1 に各実験で観測されたエネルギースペクトルの指数と，その前後でのエネルギースペクトルの折れ曲がりのエネ

20ルギーを示す．

図 1.13: 各実験で観測されたエネルギースペクトル．それぞれ Auger(Hybrid/SD)[29]，TA SD[30]，
Yakutsk SD[31]，HiRes I[32]，HiRes II[32]．
．

Yakutsk
HiRes
Auger
TA

γ1
3.29±0.17
3.25±0.01
3.27±0.02
3.33±0.04

γ22.74±0.20
2.81±0.03
2.68±0.01
2.68±0.04

γ2
5.1±0.7
4.2±0.1
4.2±0.7

log10 EA
19.01±0.01
18.65±0.05
18.61±0.01
18.69±0.03

log10 ES19.75±0.04
19.41±0.02
19.68±0.09

表 1.1: 各実験で観測されたエネルギースペクトルの指数と折れ曲がりのエネルギー．それぞれYakutsk SD[31]，HiRes I & HiRes II[32]，Auger(Hybrid/SD)[29]，TA SD[30] による報告．

ここで，EA は ankle，ES はスペクトルの suppression が始まる点である．また，γi はスペクト
ルの羃指数である．それぞれ γ1 は ankle 以下，γ2 は ankle と suppression の間，γ3 は suppression
より上のエネルギー領域での羃指数である．
図 1.14 に図 1.13 の各実験のエネルギースケールを調整してプロットした結果を示す．それぞれについて，Auger は +10.2%，TA は −9.4%，Yakutsk は −43.9%，HiRes I は −8.9%，HiRes IIは −9.7% の調整が行われている．図 1.15 に図 1.14 の縦軸の頻度にエネルギーの 3 乗を掛けてス
ペクトルの構造をわかりやすくしたプロットを示す．
このようにエネルギースケールの調整を行うと，各実験のスペクトルはよく一致する．各実験
スペクトルには約 1018.7 eV で羃指数が −3.3 から −2.7 に変わる曲がり (ankle) が，約 1019.6 eV から上では宇宙線の到来が急激に減る suppression が見られる．
このように異なった観測方法を持つ三つの実験 (大気蛍光望遠鏡 [HiRes]，水タンクアレイ [Auger]，
シンチレーターアレイ [TA]) で結果が一致することから ankle と suppression の存在は実験的に確21

図 1.14: 図 1.13 を各実験のエネルギースケールを調整したプロット．
．

図 1.15: 図 1.14 の縦軸の頻度にエネルギーの 3 乗を掛けたプロット．
．

22立したといって良い．
しかし，ankle と supression の原因に関わる一次宇宙線粒子種の統一的見解はまだ得られていな
い．これは一次粒子種の観測結果が各実験で異なるからである．後述するが (1.4.2 節)，Auger の測定では約 1018.6 eV から上の領域で徐々に陽子から重い原子核へと変化しているのに対し，TAの測定では 1018.2 ∼ 1019.6 eV の領域で陽子であり変化がない．一次粒子種が陽子であれば，ankle
は CMB との衝突による電子・陽電子対生成により起こり，suppression は GZK カットオフにより起こる．一次粒子種が原子核，例えば鉄の場合，suppression は CMB との衝突による原子核の光分解から起こるが，ankle の発生に関しては別の原因を考える必要がある．

1.4.2

質量組成

knee 領域以上の超高エネルギーの宇宙線の場合，その到来頻度の少なさから，直接観測が難しいためその質量組成を調べることは大変困難である．そこで，空気シャワーの観測を通じて間接的にその組成を調べる幾つかの方法が試みられている．
ミューオンの存在比
ミューオンは電子に比べ重く，ほぼ電離以外のエネルギー損失をしないため，空気シャワーの最大発達時の粒子数を保ちながら地上に到達する．その一方で電子成分は最大発達後の減衰が早い．また，原子番号 A の原子核はエネルギー E0 /A を持つ A 個の核子の組み合わせと考えられる
ため，一次宇宙線が重核のとき，空気シャワーの発達は早く減衰も早い [15]．一方，重核のシャ
ワーでは総ミューオン数 Nµ は大きい値となるので，宇宙線質量数が大きくなると，シャワー中の総ミューオン数 Nµ と総電子数 Ne の比は大きくなる [15]．
さらに Nµ /Ne からハドロンによる空気シャワーと γ 線による空気シャワーを区別することもできる．γ 線の場合，発生するミューオンは電磁カスケード中の γ 線と大気中の原子核との光パイオ
ン生成によって生じる π を親粒子としている．しかし，光パイオン生成の断面積は原子核のパイオン生成断面積と比較して約 1/1000 であるため，生成されるミューオン数は少ない．これより，
γ 線起源のシャワーの選別が可能である．
空気シャワーの最大発達深さ

Xmax は一次粒子種に依存するパラメーターであり [15]，極高エネルギーの粒子種を判別するのに重要である．観測によって求められた平均 Xmax を核種ごとのシミュレーションと比較すること
で質量組成解析を行うことができる．Auger と TA がそれぞれ大気蛍光法で得られたデータを用いて行った質量組成解析の結果を図 1.16 に示す [33][34]．Auger の測定では約 1018.6 eV から上の領域で徐々に陽子から重い原子核へと変化しているのに対し，TA の測定では 1018.2 ∼ 1019.6 eVの領域で陽子であり変化がない．

1.4.3

到来方向分布Auger は 2007 年に極高エネルギー宇宙線の到来方向分布は，VCV カタログに載っている 75 Mpc
以内の AGN との相関を示すと発表した [35]．しかしその後統計量の増加により，AGN との相関は弱まっている．図 1.17 は，Auger で観測された極高エネルギー宇宙線の到来方向分布を銀河座23図 1.16: 上図は Auger による平均 Xmax 測定結果 [33]．下図は TA のステレオ解析による平均 Xmax
測定結果 [34]．

24標で示したものである．青丸は 75 Mpc 以内の AGN の方向を表し，大きさは半径 3.1◦ の範囲を示す．この円の大きさは Auger の方向決定精度と銀河磁場の影響を考えた値である．
TA は現在までの観測から得られる統計では，VCV カタログに載っている AGN と到来方向との相関，および等方的に宇宙線が到来していると仮定したときの到来方向の分布との相関の両方に対して矛盾しないと発表している [20]．図 1.18 は，TA の 3 年間 (2008 年 5 月 11 日 ∼2011 年 9月 15 日) で得られた 57 EeV 以上のエネルギーの宇宙線の到来方向を示している．青丸と赤丸は到来方向を表し，赤丸は到来方向に AGN が存在している場合を示している．

図 1.17: Auger で観測された 55 EeV 以上の極高エネルギー宇宙線の到来方向を銀河座標系で示している．図の中の黒点が到来方向を表し，青丸は VCV カタログに載っている 75 Mpc 以内に存在する AGN を中心に 3.1 度に相当する大きさの丸で示されている [35]．

25図 1.18: TA で観測された 57 EeV 以上の極高エネルギー宇宙線の到来方向を銀河座標系で示している．図の中の青丸と赤丸が到来方向を表し，その内赤丸は到来方向に AGN が存在している場
合を示している．黒点は VCV カタログに載っている 75 Mpc 以内に存在する AGN を示す [20]．

26

第2章

宇宙線望遠鏡実験

極高エネルギー宇宙線の到来頻度は極めて低く，AGASA グループの観測によれば 100 km2 の地表に年に 1 度飛来する程度である．AGASA グループは 100 km2 の範囲に 111 台のシンチレー
ション検出器を配置し，その範囲に降り注ぐ空気シャワーを観測した．これにより，検出感度を
当時の世界レベルから一挙に約 5 倍広げることに成功し，大きな成果を導いた．しかし，建設以来十数年を経て AGASA は感度，装置性能ともに限界に達し，これ以上に起源解明に迫ることはできず，2004 年に観測を終了した．
そこで日米韓露共同の新たな観測プロジェクトである Telescope Array(TA) 実験が立案され，建設された．この TA 実験は，エネルギー決定精度にシミュレーション依存性の少ない大気蛍光望
遠鏡 (Fluorescence Detector，FD) と，稼働時間が長く大きな統計量の蓄積が見込める地表検出
器 (Surface Detector，SD) による空気シャワーアレイのお互いの利点を生かし，極高エネルギー
宇宙線の空気シャワーを観測するものである．
また，FD は HiRes 実験で，プラスチックシンチレーターを利用した SD アレイは AGASA 実験
で，それぞれ用いられた技術である．さらに，三ヶ所の望遠鏡ステーションのうち一ヶ所は，HiRes
実験で使われた望遠鏡 14 台がそのまま移設され稼働しており，AGASA 実験と HiRes との齟齬を理解する上で重要な検出器となっている．
また，TA 実験の主目的は以下の通りである．

• AGASA 実験の約 7 倍の感度を有する SD アレイによる定常観測によって，最高エネルギー
領域での観測事象数を飛躍的に増やす．
• 最高エネルギー宇宙線の粒子種を同定し，これを手がかりにその起源を解明する．
• エネルギー決定の系統誤差を約 13% に向上させ，最高エネルギー領域での宇宙線エネルギー
スペクトルを確定する．
• 約 2◦ の角度分解能と十分な統計量によって極高エネルギー宇宙線の点源天体を探索する．

2.1

宇宙線望遠鏡実験の装置

SD アレイの設置場所は設置，運用，修理作業に適した地形，植生，立地が望まれる．さらに太陽光電池パネルの有効な稼働のため高い晴天率，良好で安定した検出器通信のため平坦な地形が望ましい．また，FD の建設場所は観測時間を長く取るため高い晴天率，測定のノイズ源となる
人工光が少ない，大気の透明度が高い，周囲に視界を遮るものが無い地形が望ましい．これらの条件を考慮して TA 実験の実験場として米国ユタ州ミラード郡の西部砂漠地帯 (平均標高 1400 m，
大気深さ 860 g/cm2 ) が選ばれた．この地域は年間降雨量 250 mm，晴天率 60% と理想的な気候であり，人口光も極めて少ない．HiRes グループの大気透明度測定から，FD の建設地域は米国の標準的な砂漠よりも大気透明度が高いということがわかっている [36]．

27SD アレイは 1.2 km 間隔で並べられた 507 台のプラスチックシンチレーション検出器からなり，検出面積は約 700 km2 である．FD ステーションはアレイ周辺の三ヶ所に約 35 km 間隔で設置されている．一つのステーションに 12 台または 14 台の望遠鏡が設置され，全部で 38 台である．地
表近くで発生する霧などの影響を避けるために，ステーションは周囲から 100 ∼ 200 m 高い所に置かれている．TA の検出器配置を図 2.1 に示す．SD アレイ上空を FD の視野で覆い，到来する
宇宙線を二つの異なる検出器で観測することが可能となっている．

図 2.1: TA 実験における装置の配置図．赤丸が大気蛍光望遠鏡ステーション，青丸が SD のデー
タ収集に使われる通信塔を表す．左上は SD，左下は大気蛍光望遠鏡ステーションの写真．

2.2

地表検出器

TA の SD アレイは間隔 1.2 km の碁盤目上にならべられた 507 台の SD からなり，入射角 45◦ 以
下で検出立体角 1300 km2 sr となっている．この SD アレイは，図 2.9 のように三つの領域に分けられ，各領域ごとに一つの通信塔と無線 LAN 通信を通じて接続されている．従って，TA の SD
アレイは 3 領域の部分アレイから構成されている．各通信塔はそれぞれ Smelter Knoll(SK) 通信
塔，Black Rock Mesa(BRM) 通信塔，Long Ridge(LR) 通信塔と呼ばれる．
運用当初は，三つの部分アレイはそれぞれ独立に稼働していたが，後述する「境界トリガー」の導入により，2008 年 11 月から 3 基の通信塔は互いに無線 LAN を通じて接続され，507 台の SD
は一つの SD アレイとして運用されている．

28

2.2.1

地表検出器基本構成

SD には，1.2 cm 厚，2.0 m×1.5 m のプラスチックシンチレータ (CI 工業社製 CIMS-G2) を 2 層重ねたものが用いられている．2 枚のシンチレータからの信号のコインシデンスをとることにより，環境ガンマ線からの雑音信号と空気シャワー荷電粒子からの信号を効果的に弁別できる．なお，1 層のシンチレータは 8 枚の部分シンチレータで構成されている．図 2.2 に，SD シンチレー
タの構成図を示す．
シンチレータで発生した光は，各層ごとに 104 本の波長変換ファイバー (Kuraray 社製 Y-11)
を通じて光電子増倍管 (Electron-tubes 社製 9124SA，PhotoMultiplier Tube，以下 PMT) に導かれる．
上記のシンチレータ，ファイバー，PMT はシートで遮光した上で，厚さ 1.2 mm のステンレス製の箱に収納されている．
PMT の出力はフラッシュADC(AD9235BRU-65) により，上下層のチャンネルごとに，12 bit，
50 MHz のサンプリング周波数でデジタル化して記録されている．検出器のタイミングは，各検出器に搭載された GPS 受信機 (Motorola M12+) によって，約 14 ns の相対精度で同期されている．
電力源としては 120 W の太陽電池パネル (京セラ社製 KC120J) と蓄電池 (DYNASTY 社製 DCS100L) を用い，365 日 24 時間の運用が可能となっている．上記のように個々の SD は電力線，通信線などを用いずに完全に自立して運用されるので，環境負荷は非常に小さい．設置の際も環境負荷を最小にするべくヘリコプターで運搬される．
図 2.3 にソーラーパネルの裏側を示す．ここにはバッテリー，クーラーボックス，エレクトロニクスが置かれ，これらは砂塵，雨滴を防ぐためステンレスボックスに収納されている．バッテ
リーは温度に依存して充放電性能が変化するので，クーラーボックス内に収納し温度変化を抑えている．また，バッテリー側面にはフィルムヒーターを貼付してあり，必要に応じてバッテリーを温める．

図 2.2: TASD．左は現地に設置した SD の写真．右はシンチレータ等の設置図．

2.2.2

地表検出器電子回路

SD の電子回路は，主に以下の項目からなる．
• PMT 波形のデジタル化のための FADC．12 bit，50 MHz サンプリングのものが上下層のシンチレータ用に 2 チャンネルあるため，情報量は 1.2 Gbps となる．
29

図 2.3: ソーラーパネルの下の箱にエレクトロニクスとバッテリーが収納されている．

• 高速信号処理のための FPGA(Xilinx 社製 XC3S1000)．ゲート数は 100 万である．主に PMT
波形の処理と保存，PMT 印加電圧コントロール，シリアルインターフェースコントロール，
無線 LAN 通信の初段処理，GPS 通信の初段処理に用いられる．
• 全体機能統括用 CPLD(Xilinx 社製 XC2C256)．主に CPU のブートコントロール，FPGA
ブートコントロール，電源コントロールに用いられる．
• 複雑な処理を行うための CPU(ルネサステクノロジ社製 SH7750, 通称 SH4)．動作周波数は150 MHz である．大半の情報処理はこれを用いて行われる．
• 位置情報，絶対時間情報取得のための GPS．
• 無線 LAN 通信モデム．
図 2.4 に SD のエレクトロニクスを示す．

図 2.4: (左) 収納箱に収められた SD 用エレクトロニクス．各 SD と通信塔に設置されている．無線 LAN モジュールは，メインボードの下段にある．(右)SD 用エレクトロニクスのブロック図．

旧来の ADC+TDC ではなく FADC を用いている理由は，ゲート幅を可変にすることでシャワー
フロントから遅れてくる中性子の影響を減少させることができるためである．中性子とシンチレー

30

タとの反応確率は荷電粒子やガンマ線と比べて小さいが，反応した際のエネルギー損失は他の粒子よりも大きい．これによりエネルギー再構成の精度が落ちてしまう．また，空気シャワーが生成する中性子のフラックスは高エネルギー核反応モデルの依存性が大きく系統誤差の要因になっ
てしまう．

2.2.3

通信塔基本構成

図 2.5 に通信塔の外観を示す．通信塔は，主に SD と通信するための無指向性アンテナと電子機
器，トリガー用電子機器，DAQ 用の産業用 PC，通信塔同士が通信するための指向性アンテナと電子機器 (Motorola 社製 Canopy PTP100)，発電用太陽電池と蓄電池からなる．通信塔では，部
分アレイのレベル 2 トリガー判定を行っている．これに加え，SK 通信塔では部分アレイ境界のト
リガー判定とトリガー配信を行っている．

図 2.5: LR 通信塔の外観．太陽電池パネルにより電源供給している．

2.2.4

無線 LAN

通信塔と検出器との通信
無線 LAN モデムには IEEE802.11g 規格の機器を使用している．ただし，ノートパソコンなどで通常用いられている TCP/IP では通信オーバーヘッドをが大きくなるため，ここでは使われていない．1 台の通信塔は 1 秒間に 200 台程度の SD と通信せねばならず，許容されるオーバーヘッ
ドはわずか 5 ms である．CPU の動作周波数は 150 MHz であり，TCP/IP を用いて，5 ms 以下の31

オーバーヘッドを達成するのは極めて困難であると考え，シリアル転送装置として無線 LAN 機器を使用することとした．
この機器が持つ通信速度はストリーミングモードを除くと，1.0 Mbps，2.0 Mbps，5.5 Mbps である．通信速度が大きい方が DAQ システム構築に有利であり，かつトリガースレッショルドを下げることができるのだが，パケットロスの確率が大きければ実効的な通信速度は小さくなる. 図
2.6 に 10 km の距離で測定した通信速度を示す．送受信のパケットサイズは同じにしてある．試験の結果，オーバーヘッドは 0.9 ms，実効通信速度は機器の持つ最大通信速度の 75% 程度であることが分かった．また，5.5 Mbps モードではパケットロスとなる確率が高いことも分かった．この試験では 5.5 Mbps モードでのパケットロスは許容範囲であったが，パケットロスの確率は気象条
件その他の外的要因で変化するので，安全を見て，トリガー用の通信では 2.0 Mbps，DAQ では
1.0 Mbps で運用することとした．DAQ では 1 度の通信におけるパケットサイズが大きいために，
パケットロスが起こった際の時間のロスが大きいためである．
一般に用いられる無線 LAN 機器では自動的に通信モードを切り替える機能がついている．我々も同様の機能を実装して試験したのだが，切り替えに必要な時間がハードウェア側で 50 ms 以上かかっていたため，個別の検出器ごとに通信モードを切り替えることは行わなかった.

図 2.6: 左の図は，横軸が送受信パケットサイズ，縦軸が所要時間である．従って，直線でフィットした際の傾きの逆数が通信速度となり切片がオーバーヘッドとなる．右の図は，横軸が送受信
パケットサイズ，縦軸が通信に失敗した試行の割合である．

通信塔同士の通信
通信塔と検出器とは，2.4 GHz 帯で通信しているので，干渉を避けるため通信塔同士は異なる帯域で通信することが望ましい．そこで，通信塔同士は 5.7 GHz 帯で通信することとした．通信塔には産業用 PC(Technologic System 社製 TS-5700) が設置されているので市販の無線 LAN 機器を使用している．

2.2.5

トリガー DAQ サイクル概観

SD のトリガー DAQ サイクルを時系列にそって書き下すと，図 2.7 のようになる．上のブロッ
クはタワー PC での処理，下のブロックは各 SD における処理を表す．1∼5 はトリガーサイクル，
6∼9 は DAQ サイクルを意味する．SD は，PMT の波形信号の処理 (1)，レベル 1 トリガーの判定
32

(2) を行い，レベル 1 トリガーが発行されるとトリガーをタワー PC へ送信する (3)．レベル 1 ト
リガーを受けたタワー PC はレベル 2 トリガーの判定を行い (4)，レベル 2 トリガーが発行されると SD へトリガーを送信する (5)．SD はタワーからのレベル 2 トリガーを受けると，トリガーに対応する波形を探索し (6)，タワー PC へ探索結果を送信する (7)．探索結果を受けたタワー PC
は SD に波形送信命令を出し (8)，SD は波形をタワー PC に対して送出する (9)．ただし，DAQ
が 1 秒で終わらなかった場合は，8，9 のサイクルのみが継続される．また，オーバーヘッドを減らすために，異なる情報の送受信をまとめて同時に行っている．

図 2.7: SD トリガー DAQ パイプライン概念図．上のブロックはタワー PC での処理，下のブロックは各 SD における処理を表す．1∼9 の番号で表されている手続きについては 2.2.5 節に記述した．

2.2.6

各検出器におけるデータバッファリングとレベル 1 トリガー

物質中での粒子のエネルギー損失量には下限が存在し，このときの粒子を最小電離粒子 (Minimum
Ionizing Particle:MIP) と呼ぶ。図 2.8 に物質中での粒子のエネルギー損失を示す [37]．シンチレー
タ中での粒子のエネルギー損失が最小となるエネルギーは粒子によって異なるが，地表に到来す
る宇宙線は主にミューオンであり，様々なエネルギーを持つので，これらを平均した形での 1 MIPを測定することで SD を較正することが出来る．SD250 台で試験した結果，1 MIP の粒子の入射の結果 PMT の光電面から放出される光電子数は 24.6±0.4 であった．また，1 個のミューオンによるエネルギー損失 (=1 MIP) は一定であるため，各検出器の 1 MIP に対する応答を常にモニター
し求めることで，全 SD の応答を規格化しそろえることができる．
空気シャワーに含まれる荷電粒子の大半は放射性物質起源の放射線に比べずっとエネルギーが大きいため，上下層のシンチレータを貫通する．従って，シャワー粒子に対するトリガーは上下層のコインシデンスでかければ良い．この空気シャワー粒子入射を判定するトリガーがレベル 0
トリガーである．レベル 0 トリガーのための波形積分時間は 8 タイムスライス (=160 ns) である．
スレッショルドは 0.3 MIP に相当する 15 ADC カウントである．トリガー条件を満たした波形の記録長は 2.56 µs で，継続時間の長い信号は単に連続した複数の波形として保存される．これは無線 LAN を通じたデータ転送の際のパケットサイズを大きくしすぎないためである．
以下に，データバッファリングの手順について示す．なお，以下の手順はすべて FPGA に実装された論理回路によって行われる．

33

図 2.8: さまざまな高エネルギー粒子のさまざまな物質中でのエネルギー損失 [37]．

1. 8 タイムスライス分 (幅 160 ns) の FADC 積算値の計算を，8 タイムスライス (幅 160 ns) ず
つ移動させながら行い，GPS の 1 PPS 信号からの 1 秒間における FADC 積算値の最頻値を計算する．この値を次の 1 秒間のペデスタルレベルとする．
2. 8 タイムスライス分 (幅 160 ns) の FADC 積算値を計算し，1 タイムスライス (幅 20 ns) ずつ移動させながらこの計算を繰り返し (移動和)，移動和の値からペデスタルレベルを差し引いた値が上下層共にスレッショルド (0.3 MIP 相当) よりも大きければ，レベル 0 トリガーを発行する．トリガーから 1280 ns 以内に起こった次のトリガーは無効化される．
3. レベル 0 トリガーが発行されるとその波形を波形記録専用の SDRAM に転送する．波形はトリガーの 640 ns 以前から 1920 ns 以後まで保存される．トリガーの 1280 ns 以後から 2560 ns
以後の間に二つ目のトリガーがかかった場合，波形の重複を防ぐため二つ目のトリガータ
イミングは，一つ目のトリガータイミングの 2560 ns 後まで遅らされる．波形記録専用の
SDRAM への転送と同時に，波形の積分値，波形の最大値，1 PPS 間のクロック数，波形番号を保存する．これらの情報は 10 ms ごとに CPU と共有の SDRAM に DMA 転送される．
また，これらの情報をレベル 0 トリガーテーブルと呼ぶ．

34

レベル 1 トリガー
レベル 1 トリガーのスレッショルドは 150 カウント (=3 MIP)，積分ゲート幅は波形の記録幅と同じ 128 タイムスライス (=2560 ns) である．このスレッショルドは物理上の要請ではなく，通信バンド幅で決まっている．
レベル 0 トリガーテーブルを元にして，通信塔に送信するレベル 1 トリガーのリストを作る手
順について示す．なお，以下の手順はすべて CPU で行われる．

1. GPS の 1 PPS と同時に CPU に割り込みがかかる．割り込み処理と同時に 1 PPS 間のクロッ
ク数，トリガーレート，ペデスタルレベル，波形番号の最初と最後を FPGA から読み込む．
2. DMA 転送中であればそれが終了するまで待った後，レベル 0 トリガーの処理を開始する．
処理の内訳は，クロック数から µs への変換，DMA 転送が正しく行われたかどうかの確認，
モニターデータの生成である．クロック数をそのまま使うのではなく µs とする理由は，デー
タサイズを小さくするため，また，クロック周波数が少しずつ異なる SD 同士で同期をとるためである．同時に，波形の積分値からペデスタルレベルを差し引いた値がスレッショルド
よりも大きければレベル 1 トリガーを発行する．
3. レベル 1 トリガーがあると，µs(20 bit) と波形の積分値の上下層の平均 (4 bit) を組み合わせ，
3 byte とのデータとして記録する．また，このリストをレベル 1 トリガーテーブルと呼ぶ．

2.2.7

通信塔におけるレベル 2 トリガー

通信塔では各 SD からレベル 1 トリガーテーブルを収集し，レベル 1 トリガーが 8 µs の時間内
で 3 台以上隣接する SD から見つかると，レベル 2 トリガーを発行する．レベル 2 トリガーがあるとその時刻 (レベル 2 トリガーの判定に使用された SD の内，最も早くレベル 0 トリガーを発行した SD のトリガー時刻) は DAQ リストに登録される．即座にトリガーを各 SD に配信しない理由
は通信オーバーヘッドを減らすためである．
以下に，レベル 2 トリガー判定の手順を示す．

1. レベル 1 トリガーテーブルを SD から取得する．これと同時にモニターデータ (後述) も取得
する．DAQ リストが空でなければ最も古い DAQ リストの送信も行う．
2. 境界トリガーのため，BRM 通信塔と LR 通信塔が，部分アレイの境界に接している SD の
レベル 1 トリガーテーブル (境界トリガー情報) を SK 通信塔に送出する．
3. GPS の 1 PPS と同時に CPU に割り込みがかかる．この後, レベル 1 トリガーテーブルのµs(20 bit) と SD 座標番号 (12 bit) を組み合わせ，さらに時刻順に並び替えてアレイ全体のトリガーテーブルを生成する．1 PPS と同期させて処理を行う理由は，SD のトリガーサイ
クルと同期させることで無駄な待ち時間を発生させないためである．
4. このトリガーテーブルを 8 µs のゲート幅で走査し，レベル 1 トリガーが三つ以上ある場合
は隣接条件を満たすかどうかの確認を行う．3 台以上の隣接条件を満たしていればレベル 2
トリガーを発行し，DAQ リストへの登録を行う．また，レベル 2 トリガーは 2.2.8 節で述べる，境界トリガーシステムを通じて他の部分アレイに送信される．
レベル 2 トリガーが発行されると，各通信等はレベル 2 トリガー時刻から ±64 µs の範囲のレベ
ル 0 トリガーテーブルを SD から収集する．

35

2.2.8

境界トリガー

三つの部分アレイを独立に運用すると，主に以下の二つの問題がある．

• 部分アレイの境界ではトリガー効率が低下する．
• 境界付近に落ちた空気シャワーイベントで一つの部分アレイでのみトリガーが発行された場
合には情報が欠損する．
従って，アレイ全域でトリガー効率が一定となるように部分アレイの境界に特化したトリガーを構築する必要がある．また，他の部分アレイにトリガーを配信する機能が必要である．
SK 通信塔が三つの部分アレイから境界トリガー情報とレベル 2 トリガーを収集し，境界トリガーの発行とトリガーの配信を行う．
部分アレイ間の境界領域の設定を図 2.10 に示す [38]．
以下に，境界トリガー判定の手順を示す．

1. レベル 1 トリガーテーブルとレベル 2 トリガーを各部分アレイから収集する．
2. レベル 1 トリガーテーブルを用いて，レベル 2 トリガーと同じ方式で境界トリガーの判定(図 2.10 の水色領域) を行う．境界トリガーが見つかればレベル 2 トリガーと同時に各部分
アレイに送信する．
3. さらに, 上記のトリガー判定にかからなかったが境界をまたがって隣接 2 台の SD がトリガーをかけている場合は (図 2.11)，ライントリガーを発行し対応する二つの部分アレイに他のトリガーと同時に送信する．これは，境界に接していない SD(図 2.10 の黄色領域) を含めるとレベル 2 トリガー判定にかかる可能性があるからである．
4. 各部分アレイはライントリガーを受け取ると，境界から 2 列目 (図 2.10 の黄色領域) のトリ
ガーテーブルを走査し，レベル 2 トリガーと同じ方式で境界トリガー判定 (図 2.10 の水色領
域+黄色領域) を行う．

2.2.9

SD アレイのトリガー効率

図 2.12 に SD アレイのトリガー効率のエネルギー依存性を示す [39]．トリガー効率は E >1018.7 eV の領域で 100% となる．

2.2.10

モニター情報とモニタリングシステム

各検出器が定期的に収集 · 更新している自らの動作状況 · 状態を表すモニター情報を，任意のタイミングですべての SD から中央へ転送することは，通信帯域の狭さから不可能である．そこで，
レベル 1 トリガーテーブルの転送の際，各 SD はヘッダー情報と共にモニター情報を分割して送信している．モニター情報は 9 kByte 程度であり，トリガー DAQ に負荷をかけないよう 600 分割
して毎秒 16 Byte ずつ送信される．
以下に，モニター情報として取得されている項目を示す．収集時間間隔が，1 秒のもの，1 分のもの，10 分のものがある．ただし，収集時間間隔が 1 秒の情報は毎秒取得されるため，送信時に分割されない．つまり, 上記の 9 kByte には含まれていない．

36

図 2.9: 各通信塔 (白丸) が担当する部分アレイ領域．領域名はそれぞれ通信塔の名前からつけられ，SK 領域 (赤線)，BRM 領域 (黒線)，LR 領域 (青線) と呼ばれる．

37

図 2.10: BRM，LR，SK 各部分アレイ境界設定図．水色 SD のトリガーテーブルを SK 通信塔にある PC で転送し境界トリガー判定を行う．黄色 SD のトリガーテーブルはライントリガーパター
ンを参照するためタワーエレクトロニクスで記録しておく [38]．

図 2.11: ライントリガーパターンの確認の例．青線は境界を，赤 SD はレベル 1 トリガーのかかった SD を表す．(左)(中) の場合，赤 SD は境界をまたがって隣接しているので黄色で丸のかかれた
SD にレベル 1 トリガーがあるかを確認する．(右) の場合，赤 SD は境界をまたがっていないため，
ライントリガーパターンではない．また，ライントリガーの判定前に部分アレイにおけるレベル
2 トリガー判定は終わっており，隣接したレベル 1 トリガーがこれ以上なければどのトリガーも発行されない [38]．

38図 2.12: トリガー効率のエネルギー依存性．E > 1018.7 eV でトリガー効率 100% となる [39]．

• 毎秒取得される情報– 1 PPS 間のクロック数．GPS 及びクロックの安定性の確認に用いる．
– GPS タイムスタンプ．1 PPS が継続して出力されつづけているか，つまり SD 間で同期がとれているかどうか確認するために用いる．
– レベル 1 トリガーテーブルの数．ホットチャンネル，デッドチャンネルの識別に用いる．
• 収集時間間隔 1 分の情報
– レベル 0 トリガー数．トリガーサイクルの安定性，ゲインが適切に設定されているかの確認に用いる．
– バッテリー電圧．解析には用いられないが，SD の安定稼働のために重要な項目である．
– 充電電流．太陽電池による充電が正常かを確認する．
– 太陽電池の出力電圧．
– 温湿度．
• 収集時間間隔 10 分の情報
– バックグラウンドミューオン信号のスペクトル．SD のゲイン較正に用いられる．約
50 ADC カウントが最頻となるようにゲインが調整されている．
– ペデスタルのヒストグラム．SD の較正に用いられる．
– PMT 出力信号のピーク波高のヒストグラム．これを用いてダイナミックレンジを決める．
– GPS の状態を表すフラグ．
– 通信できている GPS 衛星数．
これらのモニター情報を確認し，不具合が生じていればメンテナンスする．モニター用アプリ
ケーションの出力グラフの一部を図 2.13 に示す．

39

図 2.13: SD0318 のモニター結果．上から，通信できている GPS 衛星数，通信状態，バッテリー
電圧，気温，ペデスタル，ミューオンピーク，トリガーレートを表す．

40

ペデスタルのヒストグラム
ペデスタルレベルを得るため，FPGA 内部では 8 タイムスライスごとに FADC カウントの 8 タ
イムスライス分の積分値を求め，それらを 1 秒分ためてヒストグラムとしている．CPU はこのヒストグラムを 1 秒ごとに読み取り，さらに 10 分間積算する．これが「ペデスタルのヒストグラム」と呼ばれるモニター情報で，10 分ごとに中央へ転送される．分解能は 1 カウント，レンジは 0 から 256 カウントまでである．
バックグラウンドミューオン信号のスペクトルレベル 0 トリガーごとに FPGA 内部では波形の積分値を求め，それらを 1 秒分ためてヒストグ
ラムとしている．積分ゲート幅は 240 ns，波高の最大値の前 80 ns，後 160 ns である．CPU はこのヒストグラムを 1 秒ごとに読み取り，10 分間積算する．レベル 0 トリガーを満たす信号の原因の大半は，シャワーと無関係なバックグラウンドミューオンによるものと考えられるので，これが，「バックグラウンドミューオン信号のヒストグラム」と呼ばれ，10 分ごとに中央に送られる．
これは，バックグラウンドミューオンのヒストグラムがモニター情報として取得されるのが 10 分
ごとであるためである．分解能は 1 カウント，レンジは 0 から 512 カウントまでである．
図 2.14 にモニター情報から得たバックグラウンドミューオンによる信号の大きさとペデスタルのヒストグラムを示す．

PMT 出力信号のピーク波高のヒストグラム
FPGA から DMA 転送された信号波形のうちピーク波高の値を CPU 側でヒストグラム化する．
分解能は 32 カウント，レンジは 0 から 4096 までである．

2.2.11

バックグラウンドミューオン信号に対応するエネルギー損失

SD ではシンチレーターの発光量から粒子数を決定するため，発光量を詳細に調べる必要がある．
シンチレーターの発光量は入射粒子によるエネルギー損失量とシンチレーターの発光効率で決まる．そこでシミュレーションを用いて検出器における入射粒子のエネルギー損失特性を調べた．
物質中でのエネルギー損失は Bethe-Bloch の式で表される．

dE
1 2me c2 β 2 γ 2 Tmax
δ
2Z 1
2
= Kz
ln
−β −
(2.1)
dx
A β2 2I2
2
ここで，ze は入射粒子の電荷量，K/A = 4πNA re2 me c2 /A = 0.307[MeV/(g/cm2 )]，Tmax は 1 回の衝突で自由電子に与えることができる最大の運動エネルギー，I は平均励起エネルギー，δ は電離損失に対する密度効果の補正である．−dE/dx の単位は [MeV/(g/cm2 )] であり，この式より物質中でのエネルギー損失量は (粒子が通過する物質の密度)×(通過距離) に比例する．
粒子が物質中を通過する距離は空気シャワーの天頂角を θ として，cosθ に反比例し，エネルギー
損失量が通過距離に比例する．従って物質中における粒子のエネルギー損失 dE について以下の式が成り立つと考えられる．

dE ∝

1
cosθ

41

(2.2)

図 2.14: SD のモニター情報から得たバックグラウンドミューオン信号のヒストグラムとペデスタ
ルのヒストグラム．緑色と赤色のヒストグラムはそれぞれ上層と下層のペデスタルのヒストグラ
ムである．また，青色と紫色のヒストグラムはそれぞれ上層と下層のバックグラウンドミューオ
ン信号のヒストグラムである．それぞれのヒストグラムの右側のエントリー (250 ADC カウントの後ろの立ち上がり) は over flow を示す．積分幅は 240 ns．バックグラウンドミューオン信号のヒストグラムではペデスタル分はすでに引いてある．

GEANT4 を用いて，TA 実験の SD にミューオンが入射した場合のエネルギー損失を計算した．
その結果の一つとして，図 2.15 に θ=0◦ ，60◦ で 1 GeV のミューオンを入射した場合のシンチレー
ター上下層でのエネルギー損失分布を示す．ランダウ分布で近似した場合のピーク値、標準偏差はそれぞれ• θ=0◦ :上層 2.00 MeV，σ=0.10 MeV，下層 2.01 MeV，σ=0.10 MeV• θ=60◦ :上層 4.15 MeV，σ=0.22 MeV，下層 4.15 MeV，σ=0.20 MeVであり，ピーク値を cosθ で補正すると• θ = 60◦ (補正後) :上層 2.08 MeV，下層 2.08 MeVとなる [40]．
また 1 GeV のミューオンは物質中でほぼ電離損失しか起こさない [41]．制動放射の反応断面積σ は以下の式に示すように粒子の質量に反比例するので，ミューオン (m=105.7 MeV) の反応断面積は電子 (m=0.511 MeV) の約 1/40000 と小さい．

σ∝



e2
mc2

2

(2.3)

そのため，物質中で制動放射を起こす確率は非常に小さく，二次粒子放出によるシンチレーター
上下層でのエネルギー損失の違いは見られない [40]．
以降では 1 MIP=2.0 MeV という値を用いる．

42

図 2.15: µ− を入射した場合のエネルギー損失分布．E=1 GeV，左が θ=0◦ ，右が θ= 60◦ である
[40]．

2.3

大気蛍光望遠鏡 (BRM ステーション/LR ステーション)

TA 実験の FD は空気シャワーによる大気蛍光を球面鏡で集光し，それを PMT カメラで撮像して観測する装置である．
TA 実験では，FD は三ヶ所のステーションに集中して設置されているが，それらのうち MD ス
テーションにある 14 基は，HiRes 実験で使用されていた望遠鏡を移設したものである．一方でBRM ステーションと LR ステーションにはこの実験のために新たに設計，製作された口径 3 mの球面鏡を持つ望遠鏡が各々12 台ずつ設置されていて，ステーション全体で方位角 108◦ ，仰角3◦ ∼33◦ を観測する (図 2.16)．
以下，本節では BRM ステーション，LR ステーションに設置されている FD について紹介する．
MD ステーションに設置されている FD については米田泰久の修士論文を参照 [42]．

2.3.1

光学系

望遠鏡の主鏡は合成球面鏡である．望遠鏡 1 基の視野は，方位角方向に 18◦ ，仰角方向に 15.5◦である．この望遠鏡を 2 段に重ねることによって，仰角は上段の望遠鏡が 3◦ ∼ 18.5◦ ，下段が18.5◦ ∼ 33◦ となり，全体では 3◦ ∼ 33◦ がカバーされる．各ステーションには，合計 12 基の望遠鏡が扇形に配置され，方位角方向の視野の広がりは 108◦ である．望遠鏡に使われる球面鏡は，六角形のセグメントミラー 18 枚からなる口径 3.3 m の合成球面鏡である (図 2.17)．各セグメン
トミラーの PMT カメラに対する取付角度は調節可能である．また，合成球面鏡の中心にはセグ
メントミラーは配置されておらず，PMT カメラ内の全 PMT のゲインを相対的にそろえるために用いる標準光源装置 (Xe フラッシャー) が取り付けられている [43]．セグメントミラーのデザインは対辺間隔が 660 mm の六角形で，厚みは 11 mm，曲率半径は 6067 ± 40 mm である．素材はTempax glass で，表面には Anodization によって Al2 O3 がコーティングされている．反射率が波
長域 330 nm ∼ 450 nm で 80% 以上が維持されるように定期的に (およそ年 1 回) 洗浄されている．
焦点面，すなわち PMT カメラの光電面におけるスポットサイズは，二次元ガウス分布の標準偏差を換算して約 12 mm 以下になるように各セグメント鏡の方向が調整されており，スポットサイ
ズは PMT の大きさ (六角形の対辺間の距離 61 mm) に比べ十分小さい．

43

図 2.16: BRM ステーションの外観 (左上)，BRM ステーション内の望遠鏡 (右上)．MD ステーションの外観と望遠鏡 (下)．

44

図 2.17: FD(左) とセグメントミラー (右) 構造と寸法 (単位は mm)．

45

2.3.2

PMT カメラ

PMT(光電子増倍管) カメラは，16 × 16 のアレイ状に並べられた 256 本の PMT から構成されており，反射鏡の焦点面に設置される (図 2.18)．PMT カメラ筐体は PMT とカメラ内のエレクトロニクスを砂や埃から保護するために密閉容器となっている．PMT は浜松ホトニクス製 R9508 で，1 本で約 1◦ × 1◦ の視野を持ち，これが撮像の分解能に相当する．PMT の出力信号は後部に取り付けられたプリアンプにより増幅され，ボックス背面にあるパッチパネルを経由してデータ収集
系に送られる．また，観測対象である大気蛍光の波長域は主に 300 nm ∼ 400 nm であるので，その他の波長域のバックグラウンド光をカットするために，個々の PMT の光電面上にはバンドパ
スフィルターとして SCHOTT 製 BG3 が取り付けられている．さらに，PMT カメラの前面窓には 300 nm 以上の波長域に対して透過率が 90% 以上のクラレ製パラグラス-UV00 を用いている．

図 2.18: PMT カメラ (左) と BG3 が取り付けられた PMT(R9508)(右)．

2.3.3

大気蛍光望遠鏡用データ収集システム

大気蛍光望遠鏡システムのためのデータ収集システムの概略図を図 2.19 に示す [44]．PMT の電源電圧は各 PMT 独立に印加されて，ゲインが 8 × 104 となるように調整されている．PMT
出力のプリアンプのゲインは 5 倍で，PMT を含めたトータルでのゲインが 4 × 105 となっている．PMT の出力信号はプリアンプで増幅され，カメラ内部に設けられたパッチパネルを経由してSignal Digitizer and Finder(SDF) モジュールに送られる．イベントトリガーの生成と信号記録には，SDF，Track Finder(TF)，Central Trigger Distributor(CTD) の 3 種類のモジュールが使用される (図 2.20)．

SDF SDF はメインアンプ，FADC からなる信号記録系と，データ収集制御のためのトリガー生成の最初の段階である信号認識のための論理回路系が実装された VME9U1 幅のモジュールである．1 モジュールで 16 入力を持ち，1 カメラあたり 16 枚の SDF を有する．メインアンプ
出力信号は逐次 12 bit，40 MHz サンプリングの FADC によってデジタル化された後，連続する 4 サンプル (100 ns) ごとに足しあわされてメインバッファメモリーに記録される．バッ
ファリングされたデータに対してフレーム (12.8 µs) ごとに信号認識判定が行われる．S/N
比が逐次計算され，これが外部から任意に設定されるしきい値を越えるか否かで大気蛍光信

46

号の有無 (ファーストレベルトリガー) が判断され，その結果が TF へ転送される．イベン
トトリガーが成立すると，時間的に対応する波形データがバッファメモリーからオンライン計算機へと 1 フレーム (12.8 µs) を 1 単位として 4 単位転送され保存される．

TF TF は 1 モジュールで 1 基の望遠鏡 (=1 カメラ=256 ch PMT) を担当する VME9U1 幅のモジュールであり，12.8 µs ごとにまとめられた画像の中に空気シャワー事象のような空間的広がりを持ち，時間的には連続した発光現象が無いかどうか判定する「トラック判定」の論
理回路を持つ．
カメラごとに全 PMT のファーストレベルトリガーの判定結果がヒットパターンとして TFに集められトラック判定される．TF はファーストレベルトリガーのパターンが次の何れかの条件を満たした場合に「空気シャワートラック」であると判断し，セカンドレベルトリ
ガーを生成する．

1. ファーストレベルトリガーが生じている PMT が 5 本以上隣接して連結している場合．
2. カメラの視野の端に，ファーストレベルトリガーが生じている PMT が 3 本以上隣接して連結している場合．
3. ファーストレベルトリガーよりもずっと大きなシグナルを記録した PMT があった場
合 (この閾値は外部から設定/変更可能である)．
TF によるセカンドレベルトリガー判定結果は逐次 CTD に転送される．
CTD CTD は VME9U1 幅のモジュールである．1 モジュールに 12 台の TF が接続される．すな
わち一つの FD ステーションに 1 台の CTD モジュールがあり，トリガー信号の分配とトリガー時刻の記録を行っている．ステーション内の全ての TF からのセカンドレベルトリガー
情報を集め，データ取得の最終判断を行い，データ収集開始のためのトリガー信号 (ファイナルトリガー) をステーション内の全信号記録系に向けて送出する．さらに CTD はカメラ
間にまたがったトラックを判定するための論理回路を持ち，またトリガー生成だけでなく全エレクトロニクスを同期させるための共通クロックの供給，初期化信号の送出などの働きも担う．また，GPS モジュールを持っておりイベント発生の絶対時刻も CTD が記録する．

2.3.4

ハイブリッドトリガーシステム

ハイブリッドトリガーは FD のトリガーシグナルを用いた SD トリガーである．ハイブリッド
トリガー導入以前，TA 実験の FD と SD アレイはそれぞれ独立のトリガーシステムにより別々にデータ収集を行っていた．しかし，それぞれの検出器の検出可能エネルギー領域は異なり，1017 eVの宇宙線は FD ではトリガーされるが，SD アレイではほとんどトリガーされない．また，このようなイベントは FD1 ステーションのみで観測される場合がほとんどであり，FD 単眼の情報から
ではシャワーのジオメトリを十分な精度で再構成できない．そこで，従来の FD と SD アレイの完
全に独立なトリガーアルゴリズムに加え，FD のトリガーと同期して SD1 台以上の波形を収集するのがハイブリッドトリガーである．ここで，波形収集の対象となるのは各 FD ステーションに対応した SD アレイ領域である．BRM ステーションは BRM 領域，LR ステーションは LR 領域，
MD ステーションは SK 領域に対応する．図 2.21 にハイブリッドトリガーの概念図を示す．ハイ
ブリッドトリガーを導入することで，低エネルギー側でのイベント再構成精度が向上し，ankle 付
近でのエネルギースペクトルの形状や，組成の解析を高精度で行うことが可能になる．ハイブリッ
ド解析は 1017.5 eV 以上のイベントに対して，コア位置 50 m，到来方向 1◦ の決定精度を持つ [45]．

47図 2.19: FD エレクトロニクスのブロック図．

(a) SDF ボード(b) TF ボード(c) CTD ボード

図 2.20: FD のための VME モジュール．(a) の SDF は望遠鏡 1 台あたり 16 枚．(b) の TF は望遠
鏡 1 台あたり 1 枚．(c) の CTD はステーション 1ヶ所あたり 1 枚．

48

図 2.21: ハイブリッドトリガーの概念図．FD によりトリガーが発行されると，FD ステーションに対応した SD アレイ領域の通信塔に波形収集命令が送られ，通信塔は領域内の SD から波形を収集する．

SD の波形検出効率
ハイブリッドトリガーを用いてレベル 0 トリガーを発行した SD の波形収集を行うと，1017.5 eVの低エネルギーのイベントについても，97.9% の確率で 1 台以上の SD の波形情報を得ることが可能である (図 2.22)[45]．

SD の波形探索時間
FD と SD アレイでは同一イベントの検出時刻が異なる．そこで FD はこの差を考慮した時刻でハイブリッドトリガーを生成しなければならない．この時刻の決定には FD と SD アレイの両方で観測されたイベントの時刻差が使用されている．図 2.23 は同一のイベントに対する FD の観測
時刻と SD アレイの観測時刻の差である．この分布から，FD の観測時刻から −35 ± 64 µs がハイ
ブリッドトリガーで指定する SD アレイの波形探索時間幅とて採用されている．

2.3.5

大気蛍光望遠鏡のキャリブレーション

大気蛍光法による空気シャワー観測から一次宇宙線のエネルギーを精度よく推定するためには，
FD の感度較正や大気状態のモニターが重要である．SDF によって記録された信号から PMT に入射した光子の個数，さらには大気中の各深さでの空気シャワーの粒子数を算出しなければならない．このためには，エレクトロニクスの応答，各部の増幅率，量子効率，望遠鏡光学系の透過率，反射率などが精度良く求められていなければならない．また，空気シャワーが発する光の発光量とそれらの光の大気中の物質による吸収，散乱を推定するためにさまざまな方法で地表から10 km 以上の上空までの大気の状態を観測中常時測定し，知っておく必要がある．そこで TA 実
験では以下のような望遠鏡光学系のキャリブレーションや大気モニターを行っている．

49

図 2.22: COSMOS を用いたシミュレーションで，1017.5 eV の陽子によるイベントでレベル 0 ト
リガーを発行した SD の台数．97.9% で 1 台以上の SD がレベル 0 トリガーを発行している [45]．

図 2.23: 同一イベントに対する FD の観測時刻と SD アレイの観測時刻の差 [45]．

50

PMT ゲインの較正
CRAYS 1 台の望遠鏡あたり 3 本の PMT については，光量の絶対値が較正されている光源を用いて，ゲインの絶対値とその印加電圧依存性が実験室において測定されている．この測定のために開発された光源が CRAYS(Calibration using RAYleigh Scattering) である (図 2.24)．
これは窒素ガスを充填したチェンバーと N2 パルスレーザー (Thermo Laser Science 製，VSL337ND-S) からなる．レーザーの発光波長は 337.1 nm で光量は 300 nJ/pulse，パルス幅は4 µs である．レーザーの出力エネルギーは放射パルスごとにエネルギーメーター (Laser Probe
Inc. Rjp-435) でモニターされており，パルスごとのエネルギーの精度は 3% である．また，我々の使用しているレーザーでは発光量のパルスごとのばらつきは ±5% であることがわかっている．PMT の視野の中心軸がレーザーの光路と正確に垂直になるように PMT がセットされ，レーザー光のレイリー散乱光のみを検出する．PMT への入射光量を正確に見積もるため，レーザー光量以外にも実験装置の位置関係，ガスの温度，圧力など正確に記録される．
結果として本測定の系統誤差は ±8% と見積もられている [46]．

図 2.24: CRAYS の外観．左図の右側からレーザーが照射され，レーザーの出力を較正するための光学系を通り，円で囲まれたチェンバーにレーザー光が入射する．左図の円で囲まれたチェン
バーの内側を右図に示す．右図の矢印はレーザー光の光路を表し，PMT は図の右側の絞りの奥に設置される．絞りは構造物による多重散乱などを遮蔽し，ビーム軸に対して 90◦ 方向の散乱光のみが PMT に入射するように設置されている．

YAP YAP は基準 PMT の光電面の中心前面に貼付されている絶対光量光源である (図 2.25)．観
測中の絶対ゲインをモニターするために用いられる．TA 実験で使われているこの光源は
YAlO3 :Ce シンチレーターを窓とする直径 4 mm，高さ 1 mm の円筒容器内に 241 Amα 線源を封入したものである．ピーク波長は 370 nm，パルス幅約 20 ns，PMT 光電面で発生する
光電子数は約 450 個で，パルスごとのばらつきは 5% 以下である．CRAYS を使った PMTのゲインの測定の際に，同時にその PMT に貼付されている YAP の光量の絶対値も CRAYSの光量との比較によって測定されている [43]．
Xe フラッシャー CRAYS による較正で増幅率と印加高電圧の関係が明らかになっている基準 PMTを基準にして，同じカメラ内の残り全ての PMT のゲインは同じ 8 × 104 にそろえられている．具体的には，基準 PMT のゲインが 8 × 104 になるように高電圧を印加したのち，カメラ面上のどの点でも同じ光子フラックスになるように調整された一様パルス光源でカメラ全体を照らし，同一のパルス光に対してすべての PMT が同じレスポンスをするように基準
PMT 以外の PMT の印加電圧を調整する．
51

図 2.25: YAP．PMT 前面中央に設置され，絶対光量光源として使用される．

この相対ゲインキャリブレーションに用いられる光源，Xe フラッシャー，は Xe ランプ (浜松ホトニクス製 L4646)，ソケット (同社製 E2418)，電子回路 (同社製 C3684) とテフロン拡散板からなる．この光源は各望遠鏡の合成鏡面中心に設置されており，主鏡から約 3 m 離れた位置からカメラ感面を照らす．図 2.26 に Xe フラッシャーの照射イメージを示す．
この光源は外部トリガーで制御されており，パルス幅 2 µs，光源から照射される光量は 2 ×
104 pe/pulse に相当し，繰り返し周波数 25 Hz で発射される．カメラ感面における光量の一様からのばらつきは 2% であることが測定されている．
相対ゲインの調整 (=印加電圧の変更) は年 3 回程度の頻度で行われている．通常はこの光源を用いて全 PMT のゲインがモニターされている (頻度は観測中の 1 時間に 1 回)[43]．

図 2.26: Xe フラッシャーの照射イメージ．

XY-スキャナー 合成鏡のスポットサイズは 12 mm 以下になるように調整されていて，これはひ
とつの PMT 感面の大きさ 61 mm より小さい．したがって，一様光源を照射して PMT 単
位でゲイン合わせをするだけでなく，PMT 感面上のゲイン分布，すなわち応答の二次元不52図 2.27: XY-スキャナー．

均一性を，スポットサイズより高い位置分解能であらかじめ知っておく必要がある．そこでUVLED(日亜化学，NSHU590B) 光源を搭載し，カメラに装着することができる大型 XY-ス
キャナーが開発された (図 2.27)．XY-スキャナーの位置決め精度は ±2 mm(実測取り付け精度の系統誤差)，±0.1 mm(位置決め精度の統計誤差) である．

この装置に搭載されている UVLED は波長 365 ± 10 nm の単色光を発し，そのパルス幅は400 ns，スポットサイズは 3.83 mm(1.65 σ ，90%)，光量は約 1000 pe/pulse である．不均一性測定は UVLED 光源をカメラ有効感度面上を 4 mm 間隔で 60 パルスずつ照射することで行われ，各測定ごとの相対ゲインの決定精度は約 2.6% である．この較正で求められたカメ
ラ感面の二次元不均一性は，観測から一次元宇宙線の情報を推定する際に利用される．具体的には，モンテカルロ法によって空気シャワーイメージと各 PMT の信号をシミュレートする際に二次元不均一性が考慮され，これらの較正結果が利用される．

大気状態のモニター
空気シャワーによって生成された大気蛍光は大気中を伝播する際に，大気中の原子，微粒子によって散乱・吸収される．この散乱・吸収量を知るために，我々はいくつかの装置を大気蛍光望
遠鏡の視野中に設置し，運用している．以下ではそれらの装置の概要を紹介する．

LIDAR TA 実験では，BRM ステーションに LIDAR 装置を設置し大気モニターに利用している．
LIDAR は LIght Detection And Ranging の略で，強いパルスレーザー光を大気に向けて発射し，その最後方散乱光強度の時間プロファイルを測定することで，ビーム射線に沿った大気の状態を推定する大気モニター法である．我々の LIDAR 装置はコンピューター制御経緯台式架台に搭載された Nd:YAG レーザー (esi 製 Orion，波長 355 nm，最大出力 5 mJ，パルス
幅 5 ns，繰り返し周波数 1 Hz)，エネルギーメーター (Coherent 社製，J50LP-1A，波長範囲
300 nm ∼12 µm)，集光鏡 (Meade 社製，LX200GPS-30)，PMT(浜松ホトニクス製，R3479)，
PMT 出力波形記録のためのデジタルオシロスコープ (Lecroy 社製，WaveRunner6039) からなる (図 2.28)．LIDAR 装置は BRM ステーションから 100 m 離れた場所にある．現在のところ LIDAR による大気モニター観測は，FD 観測の開始前と終了時に行われている [47][48]．
53図 2.28: LIDAR ドーム (右) とドーム内部の装置 (中)．ドームの中にある望遠鏡にはレーザーとPMT が取り付けられている (右)．

IR Camera 上空が低い雲でおおわれているときは，大気蛍光が雲で隠れてしまう．そのため
BRM ステーションではステーションの視野全体を赤外線カメラ (Infra-Red Camera，IR
Camera) で 1 時間ごとに撮影し監視している．雲が存在する部分は温度が高くなるので，IR
Camera を使うことで視野内の雲を撮像することができる．撮像する方向は FD 視野と同じ方向を 12 枚と水平，垂直の計 14 方向である．
図 2.29 に IR Camera とこれが撮像した写真を示す．右図は望遠鏡の視野順に 12 枚を並べた写真であり，同じ日の異なる時間に撮像したものである．上段の 12 枚の写真では雲が全体に存在するが，中段でしだいに減少し，下段の写真では雲がなくなっている様子がわかる．

図 2.29: IR Camera(左) と IR Camera で撮像された写真 (右)．右の赤外線写真では青色に比べ緑色のほうが温度が高く，雲の存在を表している．上段，中段，下段の中の各 12 枚は同時刻に撮像された写真である．上段から下段に進むにつれて時間が経過し，時間とともに雲が少なくなっていく様子がわかる．12 枚の撮影方向はそれぞれ FD の各カメラの視野の方向である．

CLF LIDAR と同様の紫外線レーザーを使用した CLF(Central Laser Facility) は三ヶ所の FD ス
テーションからほぼ等距離となる地点に設置されている (図 2.30)．CLF から垂直に射出さ54れたレーザーからそれぞれの FD ステーションの方向に散乱された側方散乱光をそれぞれのFD ステーションで観測する．この側方散乱を観測した結果を解析することでミー散乱の消散係数を求めることができる．レイリー散乱による散乱量は計算によって求めることができるので，計算で求めたレイリー散乱量と実際に観測された光量との差がミー散乱による散乱量を表すことになる．大気の状態が一次元的 (大気状態が高さにのみ依存する) と考えた場
合には LIDAR によって得られたミー散乱の消散係数と比較することが可能である．
現在 CLF は毎 FD 観測中，30 分ごとに運用され，FD で観測されている．

図 2.30: TA サイト中央に設置された CLF(左) と CLF 内部 (右)．

2.4

SD アレイと FD による同時検出事象

SD アレイはほぼ 100% の効率で稼動しているので，FD による観測が行われている時間は，SD
アレイと FD で同時観測が行われているとみなせる．しかも，エネルギー 1019 eV 以上の空気シャワーについて言えば，FD の有効検出範囲は SD アレイの展開領域を超えている．そのため，このエネルギー以上の宇宙線は同時検出されることになる．図 2.31 に SD アレイと FD による空気シャワーの同時検出の概念図を示す．
SD アレイによる観測では，計測された粒子数と粒子到来時刻を基に一次宇宙線の到来方向が決定され，シミュレーションとの比較から一次宇宙線のエネルギーが決定される．FD 観測では大気中を通過する空気シャワーの軌跡を撮像し，信号の時間差から一次宇宙線の到来方向を決定し，軌跡の輝度 (光量) とシミュレーションを比較することで一次宇宙線のエネルギーが推定される．このように SD アレイと FD は測定原理が異なり，観測された事象を解析する上で相補的な役割を果たす．SD アレイのデータ解析でのエネルギーはシミュレーションと一次組成への依存性が大きくなり，系統誤差が大きくなる．一方で，FD のデータ解析でのエネルギー推定はシミュレーショ
ンや一次宇宙線の組成への依存性は小さく，これに加えて ELS を用いたエネルギーの絶対較正をすることができるため，エネルギー決定の系統誤差は SD アレイに比べて小さい．一方で，FD 単
独での一次宇宙線の到来方向の決定は難しく (1 ステーションでのみ観測された事象の場合．二つ以上の FD ステーションで同時検出される「ステレオ事象」の場合の方向決定精度は SD アレイ
よりも良い．)，その系統誤差は ∼ 5◦ 程度であるのに対して，SD 到来方向の決定法は原理が単純
55

図 2.31: SD アレイと FD による空気シャワーの同時検出．

で，1019 eV 以上の宇宙線に対しては ∼ 1.4◦ の到来方向決定精度を持つ．SD アレイと FD による同時検出事象では，それぞれの長所を組み合わせた観測により高精度なエネルギー推定が行える．

2.5

Electron Light Source(ELS)

FD のエネルギースケールの較正には，Electron Light Source(ELS) と呼ばれる電子線形加速器を用いる．ELS は高エネルギー加速器研究機構との共同開発で製作され，2009 年夏に TA 実験サ
イトの BRM ステーション前方 100 m の位置に設置された．この ELS は，1% 以下の誤差で既知のエネルギーの電子ビームを射出できるように設計されている．
TA 実験以前の大気蛍光観測実験では，実験室などの環境で測定された大気蛍光発光スペクト
ルが，観測地の高度や気象条件でも成り立つと仮定していた．さらに大気蛍光発光効率に始まり，鏡の反射率，フィルター透過率，量子効率，増幅率など様々な係数をそれぞれ個別に測定によって求め，それらを掛け合わせて一つの係数とせざるを得なかった．しかし，我々は世界で初めて空中射出型電子線型加速器=ELS を空気シャワー観測装置に隣接して設置し，これによる電子ビー
ム，すなわちエネルギーと粒子数が実験者の指定した値にあらかじめ調整されている「擬似空気シャワー」を使って望遠鏡の各要素の係数の掛け算をひとまとめの係数として知ることを可能にし，結果としてエネルギー決定の系統誤差を小さくすることを可能にした．
ELS は BRM ステーションの正面から 100 m 離れた場所に設置されている．FD からの距離が近いため，ELS による疑似空気シャワーイベントに対する大気による蛍光の散乱，減衰は無視することができる．ELS の仕様を表 2.1 に示す．ELS から鉛直上方に射出された電子ビームは，大気分子によって散乱されながら進み，ほとんどの電子が望遠鏡の視野内で電離損失によりエネル
ギーを失って止まる．パルス中の全電子のエネルギー損失は 4×1016 eV のエネルギー損失に相当し，このとき FD で観測される大気蛍光量は 10 km 先での 4×1020 eV のエネルギー損失と等価である．空気シャワーの電磁成分はハドロン成分に比べてシミュレーションでよく再現できるため，
実際に観測された光量とシミュレーションによって得られた光量とを比較することで，大気蛍光
発光効率と FD の装置定数を一括して絶対較正できる．
56

ビーム中の電子のエネルギー

40 MeV

電荷量109 e− /pulse

ビーム強度

6.4 mJ/pulse≤ 1 Hz≤ 2 µs(typical 1 µs)

繰り返しレート
パルス幅

表 2.1: Electron Light Source の仕様．

ELS の外観とシステムの概略図を図 2.32 と図 2.33 に示す．
ELS によって生成され BRM が観測したイベントディスプレイを図 2.34 に示す．

図 2.32: TA サイトに設置された ELS．

57

図 2.33: ELS のシステム概略図．

図 2.34: ELS によって生成され BRM が観測したイベントディスプレイ．

58

第3章

モンテカルロシミュレーション

SD による極高エネルギー宇宙線観測では各 SD で検出された空気シャワーの粒子数，及び時間差から，空気シャワーのジオメトリや一次宇宙線のエネルギーを決定するので，SD の応答について常にモニターし詳しく調べ，理解しておく必要がある．しかし，入射粒子に対する SD の詳細な応答を完全に実験的に調べることは困難なので，モンテカルロシミュレーションを用いて行う．
空気シャワーシミュレーションには CORSIKA(ver.6.960)[49] を用いた．原子核相互作用については高エネルギー領域では QGSJET-II-03[50]，低エネルギー領域では FLUKA2008.3C[51] を使用し，電磁成分のシミュレーションについては EGS4[52] を使用した．また，SD 応答については
GEANT4[53] を用いた．

3.1

CORSIKA シミュレーション

CORSIKA は様々な種類の宇宙線により発生する二次粒子の振る舞いを計算できる空気シャワー
シミュレーションプログラムである [49]．プログラムは FORTRAN で書かれ，原子核，ハドロン，
ミューオン，電子，及び光子の大気中での相互作用と崩壊をシミュレートできる．

3.1.1

シニング

TA 実験で対象としている宇宙線のエネルギーは 1018 eV を超えるものであり，単一の 1018 eVの陽子による空気シャワーの二次粒子を全てトラッキングするとなると，3 Gflops の CPU を用いて 80 日程度必要である [36]．そこで一般に空気シャワー現象のシミュレーションではシニングと呼ばれる手法が使われる [54][55]．この手法は多量の二次粒子が生成されるとき，その生成数を抑制するかわりに，生成された粒子には生成されなかった粒子数に相当する重みをもたせるというものである．
しかし，この近似は粒子数や粒子の到来時刻に大きな fluctuation を生みだし，シニングしていないシャワーとは異なる数密度を与える [56]．この問題を図示したのが図 3.1(a) である．ここでは，実線で表された粒子が入射した検出器のみが粒子を検出し，点線で表された粒子 (抑制されている) の射線上にある検出器からでは，粒子は検出されない．

3.1.2

デシニング

シニング情報を持った粒子から地表での粒子情報を復元する方法をデシニングと呼ぶ．今節では TASD 用に開発されたデシニング手法について述べる [56]．はじめに，デシニングにより復元可能な情報と，CPU の処理速度との兼ね合いで，二次粒子をどれだけシニングするか (シニングレベル) を決定する．ここでは，最適なシニングレベルとして 10−6 に設定されている [56]．重みw を持つ粒子の軌道上の，任意の点を基準にしたガウス分布に従う広がりを持たせる．このガウス分布にしたがい w-1 の粒子を復元させ射出する．ここで復元された粒子のエネルギーは，重み59図 3.1: (a) シニングにより発生する問題．点線はシニングにより破棄された粒子 (実際のシャワー
では存在する)，実線は重みを持たせた粒子であり，線が太いものほど重みが大きい．黒四角は SDを表す．(b) デシニングのアルゴリズム [56]．

を持った基の粒子のエネルギー周りのガウス分布に従う．また，それぞれの粒子の速度は光速度で近似する．ここで，粒子の軌道に持たせるガウス分布の幅によって，復元された粒子のもつエネルギー分布や粒子の軌道に広がりを持たせる点の地上からの距離の最大値が決まるが，このガウス分布の幅はシニングレベル (10−6 ) に合わせて決定されている [56]．
図 3.2(a) にシニングされたシャワーと，ノンシニングシャワーの粒子数の横方向分布の比較を，
図 3.2(b) にデシニングされたシャワーと，ノンシニングシャワーの粒子数の横方向分布の比較を示す [57]．ここで一次粒子は 1019 eV の陽子，シャワーの天頂角は 45◦ である．この図より，デシニングによって粒子数の横方向分布は平均，RMS ともに良く再現されていることがわかる．

3.2

GEANT4 シミュレーション

GEANT4 は様々な素粒子の物質中での反応や振る舞いをシミュレートすることのできるソフト
ウェアである [53]．シミュレーションを C++ で記述し，物理プロセスや検出器のジオメトリー，
データ生成やその保存，視覚化等の多くのツールキットからなる．そのため，検出器の素材や形状を設定することができ，各ユーザーの目的に沿ったシミュレーションが可能である．CORSIKA
により生成された空気シャワー粒子の検出器中でのエネルギー損失の計算はこの GEANT4 によるシミュレーションで行った．今節では TASD 用に開発された GEANT4 シミュレーションについて述べる [58]．
図 3.3 にシミュレーションで設定した SD の形状と，SD のステンレスボックス内の構成を簡潔に説明するために，各構成要素をずらして，見やすく配置したものを示す [40]．各構成要素は以下のように色で種類分けされている．
• 紫色:屋根 (鉄)(1.4 mm)，屋根の支持具 (鉄)．
• 灰色:ステンレスボックス上蓋 (1.5 mm) と下蓋 (1.2 mm)，上下層分離ステンレス板 (1.0 mm)．
60

図 3.2: 1019 eV，天頂角 45◦ の陽子による空気シャワーの粒子数の横方向分布．(a) シニングレベ
ル 10−6 のシニングシャワーとノンシニングシャワーとの比較 (b) デシニングシャワーとノンシニ
ングシャワーとの比較 [57]．

• 黒色:ブラックシート (0.8 mm)• 赤色:タイベックシート 4 層 (各層 0.4 mm)．
• 青色:プラスチックシンチレーター 2 層 (各層 12 mm 厚，1500 mm×2000 mm) (1 層 8 枚 (各12 mm 厚，1500 mm×250 mm))．
• 水色:発泡スチロール (5.0 mm)
実際の SD では，ブラックシートは 0.2 mm を 4 枚重ねたもの，タイベックシートは 0.2 mm を 2枚重ねたものだが，シミュレーションではこれらを 1 層として扱っている．
さまざまな粒子種，エネルギー，天頂角の粒子を SD に入射することで，シンチレーター上下層でのエネルギー損失を計算する．図 3.4 は SD の応答を GEANT4 を用いて計算する際に，シミュレーションプログラム内に SD を記述するセットアップの概略図である [57]．ここで地面の密度は2.0 g/cm3 ，厚みは 2.3 m である．以下に入射粒子，及びその二次粒子のエネルギー損失を計算する手順を示す．

1. SD を中心とする 6 m×6 m の範囲にランダムに粒子の到来点を決定する．
2. 粒子の到来方向にランダムな方位角を与える．
3. 粒子を入射し，SD の各層でのエネルギー損失を求め記録する．
4. これらの作業を 1.2×106 回行う．
ここで粒子の到来点を SD の面積 (1.5 m×2 m) よりも 12 倍大きい 6 m×6 m としたのは，検出器の縁での粒子の相互作用や，SD 付近の地表に落ちた粒子により生成される二次粒子の影響を考慮するためである．

61

図 3.3: シミュレーションプログラム内に記述した SD の形状 (上)．SD のステンレスボックス部
分の構成 (下)[40]．

62図 3.4: SD の応答を GEANT4 を用いて計算する際に，シミュレーションプログラム内に SD を記
述するセットアップの概略図．ある粒子種，天頂角，エネルギーを与えた粒子を 6 m×6 m の範囲にランダムに降らし，方位角をランダムに与える．矢印は入射粒子のジオメトリ例を表す [57]．

このシミュレーションによって得られたエネルギー損失の二次元ヒストグラムの例を図 3.5 に示す [57]．横軸は上層でのエネルギー損失，縦軸は下層でのエネルギー損失である．このヒストグ
ラムは，両層ともにエネルギー損失がゼロである場合は含んでいないことに注意する．図 3.5(a)は 1 GeV のミューオンを天頂角 0◦ で入射した時のヒストグラムである．上下層ともに約 2 MeVに強いピークをもつ．図 3.5(b) は 1 GeV のミューオンを天頂角 60◦ で入射した時のヒストグラム
である．上下層ともに約 4 MeV にピークをもつ．片方の層のエネルギー損失が約 4 MeV であり，
もう一方の層のエネルギー損失が幅をもっているのは，粒子が片方の層を通過し，もう一方の層をかすめる，もしくは全く通過しないケースが存在するためである．図 3.5(c) は 1 GeV の γ 線を
天頂角 0◦ で入射した時のヒストグラムである．上下層ともに約 4 MeV にピークをもつ．エネル
ギー損失がミューオンに比べて 2 倍となるのは，γ 線は SD の屋根や内部で対生成を起こすためである．上層のエネルギー損失が幅をもっているのは，上層で対生成された e+ や e− が下層を通過するケースによるものである．図 3.5(d) は 1 GeV の γ 線を天頂角 60◦ で入射した時のヒストグラ
ムである．この時，SD 付近の地表に落ちた γ 線により生成される二次粒子の影響は 0◦ の場合に比べ大きい．また，図 3.5(c) と図 3.5(d) の上下層ともに 0.3 MeV 以下の領域での頻度上昇 (左下)は，検出器の縁での相互作用や，SD 付近の地表に落ちた γ 線により生成された二次粒子の影響である．

63

図 3.5: SD 内でのエネルギー損失の二次元ヒストグラムの例．横軸は上層でのエネルギー損失，縦軸は下層でのエネルギー損失を，色は頻度を示す．(a) エネルギー 1 GeV，天頂角 0◦ の µ+ ．(b)エネルギー 1 GeV，天頂角 60◦ の µ+ ．(c) エネルギー 1 GeV，天頂角 0◦ の γ 線．(d) エネルギー
1 GeV，天頂角 60◦ の γ 線 [57]．

64

第4章

空気シャワー粒子の到来時刻の遅れ及びその標準偏差の従う関数作成

今章では，前章までに述べたシミュレーションを用いた空気シャワー粒子の到来時刻分布関数と，その分散関数の導出について述べる．

4.1

空気シャワー構造

図 4.1 にシャワーの模式図を示す．ここで T0 はシャワーコアの地表到来時刻，TD はシャワーの曲率による粒子のシャワー平面からの到来時刻の遅れ，TPlane は T0 の時点でのシャワー平面と SD間の時間に換算した距離，R は SD のシャワー軸からの距離，ρ は粒子数密度 (検出粒子数/3[m2 ])，ti は SD にシャワー粒子が入射した時刻．またここで幾何学的に TD = ti − T0 − TPlane である．
ここでシミュレーションで生成した空気シャワーを分析することにより，シャワー粒子の到来
時刻分布 (平均 TD ) を SD のシャワー中心からの距離 R，粒子数密度 ρ，シャワーの天頂角 θ，空気シャワー軸を地面に射影した半直線を基線とした際の SD の方位角 φ(図 4.2) をパラメータとして関数化する．またシャワー粒子の到来時刻の標準偏差分布 (TD の標準偏差) をシャワー曲面の前方と後方それぞれについて，同様のパラメータで関数化する．

図 4.1: 空気シャワーの模式図．T0 はシャワーコアの地表到来時刻，TD はシャワーの曲率による粒子のシャワー平面に対する到来時刻の遅れ，TPlane は T0 の時点でのシャワー平面と SD 間の時間に換算した距離，R は SD のシャワー軸からの距離，ρ は粒子数密度 (検出粒子数/3[m2 ])，ti はSD にシャワー粒子が入射した時刻．

図 4.3 に，1019.0 eV 陽子のシャワーについて，シミュレーションから得られた TD とその標準偏差を示す．天頂角が大きい，または SD の方位角 φ が大きいことはシャワーの通過してきた物質

65

図 4.2: 座標系の模式図．太線は空気シャワー軸，点 O は空気シャワーのコア位置，点 A は SD の位置を示す．点線は空気シャワー軸を地面に射影したもの．

が多いこと (シャワーが年老いていること) を意味する．図 4.3 から，シャワーエイジが大きいほど TD とその標準偏差は小さくなることがわかる．これは，シャワーエイジが大きいほど，検出されるシャワー粒子の発生点が相対的に遠くなることが原因であると考えられる．また，シャワー
エイジが大きくなることでシャワー中の電磁成分の比率が減少することも原因の一つであると考えられる．これは，シャワーエイジが大きくなると，シャワーフロントの持つ厚みの後方部分に多く分布している電磁成分が大気に吸収されて減少してしまうためである．
同様に 1019.5 eV と 1020.0 eV の陽子のシャワーについても図 4.4 と図 4.5 に示す．これらの図の比較から，一次粒子のエネルギーが大きくなると TD とその標準偏差が小さくなっていることがわかる．これは，粒子数密度が大きくなることでシャワーフロント前方の粒子が検出される確率が高くなるためである．

4.2

シャワー粒子の到来時刻の遅れの関数化

TD を関数化するにあたって参考にしたものは，以下の AGASA グループが実験的に求めた経験式である [27]．
TD = 2.6 1 +R
3000[cm]aρb × 10−9 [sec](4.1)

ここで，a=1.5，b=−0.5．
a，b を θ，φ の関数として，1019.0 eV の場合の TD (図 4.3) をフィットした．

TD = 2.6 1 +R
3000[cm]

1.78−0.25secθ−0.00065φ

ρ−1.1+0.5secθ−0.0005φ × 10−9 [sec](4.2)

また，この関数が使用できる範囲はシャワー軸からの距離が 500 ∼ 2500 m，シャワーの天頂角が0◦ ∼ 60◦ である．フィット時の粒子数密度はシミュレーションから求めた各領域 (R，θ，φ) での平均値を用いた．図 4.6 に 1019.0 eV 陽子のシャワーについて，シミュレーションから得られた粒子数密度とその標準偏差を示す．

66

4.3

シャワー粒子の到来時刻の標準偏差の関数化

TD の標準偏差を関数化するにあたって参考にしたものも，以下の AGASA グループが実験的に求めた経験式である [27]．
aRTS = 2.6 1 +ρb × 10−9 [sec](4.3)
3000(cm)
ここで，a=1.5，b=−0.3．
図 4.3 に示すように，すなわち図 4.3 の誤差棒の長さがマーカーの上下で異なることからわかる
ように，TD の標準偏差はシャワー曲面に対して前後非対称であり，それぞれに対してこの標準偏差の関数を導出する必要がある．ここでシャワー曲面の前方と後方に対する TD の標準偏差の関数をそれぞれ TSFront ，TSBehind とする．a，b を θ，φ の関数として，1019.0 eV の場合の TSFront ，TSBehind (図 4.3 の誤差棒) をフィットした．

TSFront = 2.6 1 +TSBehind = 2.6 1 +R
3000[cm]

R
3000[cm]

1.49−0.19secθ−0.0008φρ−0.45+0.002(−1+secθ)φ × 10−9 [sec]1.38+(0.00079−0.0011secθ)φ(4.4)

ρ−0.3−0.2secθ+0.003(−1+secθ)φ × 10−9 [sec](4.5)
また，これらの関数が使用できる範囲はシャワー軸からの距離が 500 ∼ 2500 m，シャワーの天頂
角が 0◦ ∼ 60◦ である．フィット時の粒子数密度は TD をフィットした時と同様に，シミュレーショ
ンから求めた各領域 (R，θ，φ) での平均値を用いた (図 4.6)．
図 4.7 と図 4.8 に，それぞれ 1019.0 eV 陽子と 1020.0 eV 陽子のシャワーについて式 (4.2)，式(4.4)，式 (4.5) より得られた TD ，TSFront ，TSBehind を，シミュレーションから得られた TD ，TSFront ，TSBehind と比較したものを示す．これらより，一次エネルギーが異なる場合でも，作成した関数がシミュレーションにより得られた値をフィットできていることがわかる．これは R と ρ によって一次エネルギー依存を適切に表現できているためである．

67

68

図 4.3: エネルギー 1019.0 eV 陽子による TD とその標準偏差 (誤差棒で表示)．それぞれ 0◦ ≤ φ <30◦ (赤丸)，30◦ ≤ φ < 60◦ (黄緑四角)，60◦ ≤ φ < 90◦ (青三角)，90◦ ≤ φ < 120◦ (桃丸)，120◦ ≤φ < 150◦ (水色四角)，150◦ ≤ φ ≤ 180◦ (緑三角)．

69

70

図 4.4: エネルギー 1019.5 eV 陽子による TD とその標準偏差 (誤差棒で表示)．それぞれ 0◦ ≤ φ <30◦ (赤丸)，30◦ ≤ φ < 60◦ (黄緑四角)，60◦ ≤ φ < 90◦ (青三角)，90◦ ≤ φ < 120◦ (桃丸)，120◦ ≤φ < 150◦ (水色四角)，150◦ ≤ φ ≤ 180◦ (緑三角)．

71

72

図 4.5: エネルギー 1020.0 eV 陽子による TD とその標準偏差 (誤差棒で表示)．それぞれ 0◦ ≤ φ <30◦ (赤丸)，30◦ ≤ φ < 60◦ (黄緑四角)，60◦ ≤ φ < 90◦ (青三角)，90◦ ≤ φ < 120◦ (桃丸)，120◦ ≤φ < 150◦ (水色四角)，150◦ ≤ φ ≤ 180◦ (緑三角)．

73

74

図 4.6: エネルギー 1019.0 eV 陽子によるシャワーの粒子数密度とその標準偏差 (誤差棒で表示)．それ
ぞれ 0◦ ≤ φ < 30◦ (赤丸)，30◦ ≤ φ < 60◦ (黄緑四角)，60◦ ≤ φ < 90◦ (青三角)，90◦ ≤ φ < 120◦ (桃丸)，120◦ ≤ φ < 150◦ (水色四角)，150◦ ≤ φ ≤ 180◦ (緑三角)．

75

76

77

78

79

図 4.7: 1019.0 eV 陽子のシャワーについて，作成した関数により得られた TD ，TSFront ，TSBehindと，シミュレーションから得られた TD ，TSFront ，TSBehind との比較を表す．

80

81

82

83

84

図 4.8: 1020.0 eV 陽子のシャワーについて，作成した関数により得られた TD ，TSFront ，TSBehindと，シミュレーションから得られた TD ，TSFront ，TSBehind との比較を表す．

85

第5章

ジオメトリ再構成プログラム

地表検出器 (SD) から得られる情報は主に，粒子の到来時刻，シンチレーター内部でのエネル
ギー損失，SD の座標である．これらの情報を用いて，宇宙線の到来方向や空気シャワーコアの落
ちた位置などのジオメトリを再構成する．このジオメトリ再構成は以下の手順で行われる．

5.1

検出器レベルでの解析

ここでは，SD に記録された波形の FADC 値の情報から粒子の到来時刻と粒子数を決定する手
順について述べる．

5.1.1

粒子の到来時刻の決定

SD に記録された波形の FADC 値の情報から粒子の到来時刻を決定する手順について述べる．
各 SD はレベル 0 トリガーがかかってから，−640 ns∼1920 ns までの範囲の FADC 値を記録しており，記録された波形の先頭時刻は粒子の到来時刻よりも早くなるため，波形の先頭時刻をそのまま粒子の到来時刻としてはいけない．そのため，信号の立ち上がり時刻を粒子の到来時刻とし，信号の立ち上がりから立ち下がりまでを一つの波形とする．以下にその手順を示す．
1. 波形の先頭から 8 タイムスライス分 (幅 160 ns) の FADC 積算値 (以下，本節では「和」と呼ぶ) を計算する．
2. 1 タイムスライス (幅 20 ns) ずつ移動させながら和を計算し，和の値から (ペデスタルレベル)×8 を差し引いた値が 7 よりも大きければ最初のビンを波形の立ち上がり時刻とする．ここで値を 7 とした理由は，ペデスタルレベルの標準偏差よりも十分大きい値としたためである (ペデスタルレベルの標準偏差の典型値 0.5×8(タイムスライス)×1.5=6 より大きい)．
3. ここで決定した到来時刻の 1 タイムスライス前の FADC 値が (ペデスタルレベル)+(ぺデスタルレベルの標準偏差)×1.5 以上であれば，波形の立ち上がり時刻を 1 タイムスライス戻す．
4. 手順 3．を条件を満たす限り繰り返す．
5. 手順 3．と同様の作業を時刻を進める方向についても行い，波形の立ち下がり時刻を決める．
6. これらの作業を波形の記録幅 (2560 ns) が終了するまで行い，上記の条件を満たす波形が複数存在する場合はそれぞれの波形について波形の立ち上がり時刻と，立ち下がり時刻を決定する．
この作業は上下層に対し，それぞれ行う．
図 5.1 に上記の手順で選択した波形の範囲の 1 例を示す．

86

図 5.1: SD で記録された波形の 1 例．緑色矢印は信号と認識された範囲，赤色矢印は波形の立ち上がりを表す．

5.1.2

粒子数の決定1. 上記の方法で決定した波形と認識された範囲，すなわち信号の立ち上がりから立ち下がりまでの間の FADC 値の積算値を求め，さらにペデスタル分を差し引く．
2. この値を検出器ごとのゲイン (較正係数 [FADC カウント/MeV]) を用いてエネルギー表記(MeV) に直す．
3. エネルギー表記から粒子数に換算する (1 MIP=2.0 MeV)．
4. 複数の波形が存在する場合はそれぞれの波形について 1．∼3．の計算を行う．
この作業は上下層に対し，それぞれ行う．
ただし以下では，SD の記録した波形の立ち上がり時刻と粒子数とに上下層の平均値を用いる．

5.2

検出器レベルでのクオリティーカット

検出器レベルでのクオリティーカットは以下の項目について行った．

• 不調と判断される SD を除去する．SD が正常であるかどうかの判断基準は以下の通りである．
– SD のゲイン (FADC カウント/MeV) が 6 以上 54 以下 (典型値 17)– SD の発光量 (光電子数/MeV) が 1 以上 (典型値 10)– FADC カウントのペデスタルが 1 以上 20 以下 (典型値 5)– レベル 0 トリガーレートが 500 Hz 以上 1000 Hz 以下 (典型値 740Hz)– 1 PPS 間のクロック数のばらつきによる時刻の不定性が 100 ns 以下• 通信エラーが発生していた SD を除去する．
表 5.1 に 2010 年の不調 SD の出現率を示す．年間の不調 SD の出現率は 1.9% であり，出現率の最も高い 7 月でも 95% の SD は常時安定稼働している．

87

月

不調 SD 出現率

月

不調 SD 出現率

月

不調 SD 出現率

1月
2月
3月
4月

1.4%
1.4%
1.6%
1.4%

5月
6月
7月
8月

0.9%
1.3%
5.0%
2.0%

9月
10 月
11 月
12 月

4.3%
0.7%
1.1%
1.9%

表 5.1: 2010 年の不調 SD の出現率を表す．年間を通じて 1.9% の出現率であり，多くは通信エラーによるものである．7 月，9 月の出現率が他の月に比べて高いのは，雷雨に伴う SD のリブートのためである．

5.3

モデル関数によるジオメトリ再構成

前章で求めた空気シャワー関数と，前節までで得た粒子の到来時刻，粒子数を用いて空気シャ
ワーのジオメトリを再構成する．
まず，初めにコア付近の SD に関して再構成に使用する波形を決定するために以下の作業を行う．

a. 粒子数 1.4 以上の SD を全て選択する．
b. 3×3 台の領域で，a．の条件を満たす SD が最も多く存在する領域を選択する．この際，3×3台の領域の中心は a．の条件を満たす SD である必要がある (図 5.2(左))．ここで，a．の条件を満たす SD の台数が同じ領域が二つ以上存在する場合，SD で記録された粒子数の合計が最も大きい領域を選択する．
c. 手順 b．で選択した領域内の SD それぞれについて，隣接している SD との波形の立ち上がり
時刻の差が，SD 間の距離を R，光速度を c として R/c 以下であるかを確認する．もし R/c以上の波形の立ち上がり時刻の差がある SD が存在すればその台数を記録する．もし領域内の SD 全てに関してこれが 1 台以下であるならば，手順 d．は行わず以下の解析へ進む (図5.2(右))．
d. 手順 c．で最も多くの SD と波形の立ち上がり時刻の差が R/c 以上だった SD に対して，複数の波形が存在するかを確認する．存在していた場合は，現在使用している波形の直後に存在する波形を使用して，手順 c．をもう 1 度行う．この際，波形の立ち上がり時刻の差がR/c 以上となる隣接 SD が 1 台以下となる波形が存在する場合は，以下の解析でその波形を使用する (図 5.3(中))．他の波形が存在しない，または立ち上がり時刻の差が R/c 以上とな
る隣接 SD が 1 台以下になる波形が存在しない場合は，その SD を除去して手順 c．をもう 1度行う (図 5.3(右))．ここで除去された SD は後の解析でも使用しない．この作業は手順 c．
で示した作業の終了条件を満たすまで繰り返し行う．
ジオメトリ再構成は以下の χ2 関数を用いて，χ2 /dof を最小化するパラメーターを探すことに相当する．この関数におけるフリーパラメーターはシャワーコアの地表到来時刻 T0 ，シャワー軸の天頂角 θ，方位角 φ の 3 つである．

χ2 =

X

ξi2 (T0 , θ, φ)

(5.1)

i

ξi2 =

(ti − T0 − TPlane − TD )2TS2 + (20 × 10−9 )2
88

(5.2)

図 5.2: (左) 粒子数 1.4 以上の SD が最も多く存在する領域選択の例．赤四角は粒子数 1.4 以上のSD を表す．緑四角で囲まれた領域には，緑丸の SD を含む 6 台の粒子数 1.4 以上の SD が存在す
る．一方，青四角で囲まれた領域には，粒子数 1.4 以上の SD が 7 台存在するが，その中心である
青色丸の SD が粒子数 1.4 を検出していないので，このようなケースでは緑四角の領域を選択する．(右)SD 間の波形の立ち上がり時刻差が R/c 以上となっている例．数字は波形の立ち上がり時
刻の差が R/c 以上の隣接 SD の台数を表す．この絵は中心の SD と右上の SD の波形の立ち上がり
時刻の差が R/c 以上であり，その他の SD 同士では R/c 以下であるケースである．この場合，全ての SD に関して数字が 1 以下なので，再構成へと進む．

図 5.3: (左) 波形の立ち上がり時刻の差が R/c 以上である隣接 SD が 2 台以上の SD が存在する
ケース．数字は波形の立ち上がり時刻の差が R/c 以上の隣接 SD の台数を表す．この場合青丸のSD が隣接 4 台全ての SD と波形の立ち上がり時刻の差が R/c 以上である．このようなケースでは青丸の SD に複数の波形が存在しているかを確かめる．(中) 左図の青丸の SD に複数の波形が存在
し，隣接 SD との波形の立ち上がり時刻の差が R/c 以下となる波形が存在していたケース．(右)左図の青丸の SD に他の波形が存在しない，または立ち上がり時刻の差が，R/c 以上となる隣接SD が 1 台以下になる波形が存在しなかったケース．この場合，青丸の SD は以下の解析から除去
される．

89

ここで分母の 20 × 10−9 [sec]) は SD の時間分解能である．また，TS は ti − T0 − TPlane − TD が負の時は TSFront (式 (4.4))，正の時は TSBehind (式 (4.5)) を使用する．
初めに，パラメーターの初期値を決定する．T0 の初期値は最大粒子数を検出した SD の波形の立ち上がり時刻とする．θ，φ の初期値は，検出粒子数の多い SD3 台を用いて決定する．図 5.4 の→ようにフロントが平面の空気シャワーを仮定して，以下の方程式から法線ベクトル −n を求め，ジオメトリの初期値とする．

図 5.4: シャワーフロントを平面とする空気シャワーの模式図．シャワー軸の方向ベクトルの初期値を計算する際に，このように仮定する．

−→→ + c(t − t ) · −→→) = 0n · (−
x
n −−
x
1
1
0
0

−
→
→ + c(t − t ) · −→→) = 0n · (−
x
n −−
x
2
2
1
1
−
→
→
n ·−
n =1

(5.3)
(5.4)
(5.5)

→
ここで −
xi ，ti はそれぞれ各 SD の位置ベクトル，波形の立ち上がり時刻である．
ここで得られた T0 ，θ，φ を初期値とし，以下の手順でジオメトリ再構成を行う．またフィッティ
ングツールには MINUIT を使用した [59]．
P√
P √
ρi )，
1. シャワーのコア位置 X ，Y を粒子数の重みから次のように求める．X = ( xi ρi )/(
P√
P √
ρi )．
Y = ( yi ρi )/(2. シャワー中心から 2500 m の範囲で図 5.5(左) に示すように 3×3 台の SD で中心を含む 6 台
以上がレベル 0 トリガー条件を満たしている領域があるかを確かめ，これを満たすイベントのみを選択する．このとき，その領域内でレベル 0 トリガーを満たした SD の台数を，以下では Ncluster とする．ここで Ncluster のしきい値を 6 台にした理由は次章で述べる．
3. シャワー中心から 500 m 以遠かつ 2500 m 以内で，検出粒子数 1 以上の SD を到来方向の再構成に用いる．以下では，このような SD の台数を NSD とする．この選択は，以下で複数回行うフィッティングの度に行う．またこれが 6 台以上であるイベントのみを選択する．NSDのしきい値を 6 台にした理由は次章で述べる．

90

4. χ2 /dof を最小にするパラメーターを探す．この際に動かすパラメーターは θ，φ の二つである．
5. χ2 /dof を最小にするパラメーターを探す．この際に動かすパラメーターは T0 ，φ の二つである．
6. χ2 /dof を最小にするパラメーターを探す．この際に動かすパラメーターは T0 ，θ の二つである．
7. 手順 4．で χ2 /dof が収束しなかった場合のみ，手順 5．6．を再び行う．χ2 /dof の収束条件
は EDM(expected distance to minimum)≤ 10−5 である [59][60]．
8. 手順 5．で χ2 /dof が収束しなかった場合のみ，手順 4．6．を再び行う．
9. 手順 6．で χ2 /dof が収束しなかった場合のみ，手順 4．5．を再び行う．
10. 手順 4．を再び行う．
11. χ2 /dof が 4.5 以上かつ ξi2 が最も大きい SD の ξi2 が 4.5 以上のとき，その SD に再構成に使用した波形より後ろに別の波形が存在しないかを確認する．もし波形が存在すれば，現在使用している波形の一つ後ろの波形を用いて手順 3．以降の再構成をもう 1 度行う．存在しな
い場合はその SD を除去して手順 3．以降の再構成をもう 1 度行う．この作業は 10 回を上限に，χ2 /dof が 4.5 以上かつ ξi2 が最も大きい SD の ξi2 が 4.5 以上である限り行う．その場合10 回のジオメトリ再構成の中で，最も χ2 /dof が小さくなった場合の結果を選択する．また，
χ2 /dof のしきい値を 4.5 とした理由は次章で述べる．

図 5.5: イベント再構成のための条件．赤四角はレベル 0 トリガー条件を満たした SD を示す．左の場合は再構成を行う．右の場合は再構成を行わない．

5.4

イベントレベルでのクオリティーカット

以下の項目を満たす事象のみを，再構成に成功した事象とする．

• 空気シャワーのコア位置がアレイ境界の内部にあり，かつコア位置が境界から 1.2 km 以上離れている．
• 天頂角が 60◦ 以下である．
91

• χ2 /dof が 4.5 以下である．
• pointing error が 2.5◦ 以下である．
q
ここで pointing error は到来方向の fiting error σθ2 + sin2 θσφ2 である．また，χ2 /dof 及び pointing error のしきい値をそれぞれ 4.5 と 2.5◦ に決定した理由は次章で述べる．

92

第6章

モンテカルロシミュレーションによる到来
方向及びジオメトリ決定精度の評価

前章までに作成したジオメトリ再構成プログラムを使用して，シミュレーションにより生成されたイベントの到来方向及びコア位置を含むジオメトリを再構成し，その精度を定量的に評価した．

6.1

イベントセット

ここでは，再構成精度評価のために用いたシミュレーションのイベントセットについて述べる．
一次粒子は陽子を仮定した．一次エネルギーは 1018.5 eV，1019.0 eV，1019.5 eV ，1020.0 eV それぞれに対して log スケールで ±0.1 の範囲でランダムとした．天頂角は図 6.1 上図に示すように，
0◦ ≤ φ ≤ 60◦ の範囲で，sinθcosθ 分布に従う乱数で生成した．方位角は図 6.1 中図に示すように，
0◦ ≤ φ ≤ 360◦ の範囲でランダムに分布するよう決めた．コア位置は図 6.1 下図に星で示す CLFの位置を中心とする半径 25 km の範囲にランダムに分布させた．

6.2

FD の観測結果を用いたエネルギー補正

大気蛍光望遠鏡 (FD) による観測から決定される一次エネルギーには，モンテカルロ計算に必然的に含まれるハドロン相互作用の系統的不定性が含まれない．一方，SD 観測からの一次エネル
ギー決定はハドロン相互作用による不定性を持つ．TA 実験では，ハイブリッド観測を行うこと
で，SD アレイで測定されたエネルギーを FD 観測から補正することができる．TA 実験のこれまでの観測から，TA 実験の SD アレイは FD に比べエネルギーを 27% 大きく見積もることが明らかになっている [57]．SD アレイ観測による一次エネルギーにはこの補正を加える．図 6.2 にこの補正を加えた後の，SD アレイと FD により測定されたエネルギーの比較を示した．

6.3

ジオメトリ再構成プログラム内の種々のしきい値決定

5.3 節の再構成手順 2．及び 3．で示した SD の台数条件 Ncluster ，NSD のしきい値と，手順 11．
での χ2 /dof のしきい値を解析結果から決定する．
それらのしきい値のためのシミュレーションでは，一次粒子種は陽子とし，エネルギーは到来頻度と検出効率からイベント数が多くなる 1019 eV とした．
ここで注意がある．一般に空気シャワーアレイによる観測から空気シャワーのコア位置を決定するためには，シャワーサイズと粒子数の横方向分布を用いる必要がある．そのため，後述する
単純な粒子数の重心をコア位置とする方法では，200 m 程度の精度しか得られない．しかし，本研究では到来方向再構成の精度を向上させることに重点を置くため，Ncluster と NSD に課す条件，
及び χ2 /dof のしきい値は，シミュレーションでのコア位置を再構成ソフトウェアにもそのまま与えた場合，すなわちコア位置が厳密にわかっているとした場合の解析から決定した．
93

図 6.1: シミュレーションイベントを生成するために用いたシャワーパラメーターの分布．上図は
天頂角分布，中図は方位角分布，下図はコア位置分布を示す．

94

図 6.2: 27% の補正を加えた後の SD と FD により測定 · 解析された実イベントのエネルギーの比較．左は散布図，右はエネルギー比の分布．

6.3.1

χ2/dof のしきい値決定

初めに，前章の再構成プログラムの手順 11．で用いる χ2 /dof のしきい値を決定する．ここでは，
手順 2. 及び 3. における Ncluster ，NSD によるカット及び，手順 11．は行わない．図 6.3 に 1018.5 eV，1019.0 eV，1019.5 eV，1020.0 eV の一次陽子の場合の χ2 /dof のヒストグラムを示す．1019.0 eV の場
合の χ2 /dof 分布から，総数の 68% が含まれる χ2 /dof=4.5 をしきい値とした．

6.3.2

Ncluster 及び NSD のしきい値決定

次に，Ncluster と NSD の条件を決定するために，Ncluster ，NSD のセレクション条件を変えながら，手順 2．3．のカットに対する通過率と到来方向決定精度を求めた．この際，手順 11．も含めた再構成を行っている．図 6.4 に，Ncluster と NSD によるカットがない場合の到来方向決定精度のヒストグラムを示す．また，全イベントの 68% が含まれる到来方向決定精度の上限 (以下，到来方向決定精度) と，カットの通過率を表 6.1 に示す．

Ncluster ≥ 5
≥6
≥7
≥8
=9

NSD ≥ 41.80◦ (97.1%)
1.78◦ (93.9%)
1.74◦ (77.1%)
1.68◦ (41.4%)
1.69◦ (9.9%)

≥5

≥6

≥7

≥8

1.73◦ (91.1%)

1.62◦ (76.8%)

1.49◦ (53.6%)

1.34◦ (29.9%)

1.72◦ (89.1%)
1.69◦ (74.6%)
1.65◦ (40.7%)
1.68◦ (9.8%)

1.62◦ (75.9%)
1.61◦ (66.5%)
1.60◦ (38.2%)
1.65◦ (9.4%)

1.48◦ (53.4%)
1.48◦ (49.2%)
1.51◦ (31.2%)
1.58◦ (8.4%)

1.34◦ (29.8%)
1.34◦ (28.6%)
1.37◦ (20.3%)
1.47◦ (6.2%)

表 6.1: Ncluster ，NSD に課すセレクション条件ごとの 1019.0 eV 陽子のシャワーの到来方向決定精度の 68% 上限レベル．括弧内はこれらのセレクションを通過したイベントの割合を表す．

最終的に，到来方向決定精度と，セレクションを通過するイベントの割合から，NSD =6 とNcluster =6 を再構成のセレクション条件とした．

95

図 6.3: 陽子の到来方向再構成の χ2 /dof のヒストグラム．それぞれ (左上)1018.5 eV，(右上)1019.0 eV，(左下)1019.5 eV，(右下)1020.0 eV．破線は全イベントの 68% が含まれる χ2 /dof の値である．

図 6.4: 1019.0 eV 陽子のシャワー再構成の到来方向決定精度のヒストグラム．

96

6.4

クオリティーカットにおける χ2 /dof 及び pointing error のしきい
値決定

5.4 節のクオリティーカットにおける χ2 /dof 及び pointing error のしきい値を解析結果から決定する．ここでもしきい値はシミュレーションでのコア位置を再構成ソフトウェアにも与えた場合の解析から決定した．図 6.5 に 1018.5 eV，1019.0 eV，1019.5 eV，1020.0 eV の陽子に対する χ2 /dofのヒストグラムを示す．NSD と Ncluster によるセレクションを行っているため，図 6.3 と比べてイベント数が減少している．さらに，χ2 /dof≥4.5 の場合に手順 11．(ξ 2 が大きい SD が他に波形を持っていればその波形を使用して，持っていなければその SD を除去して到来方向再構成を再び行う) を行っているため，図 6.3 で χ2 /dof ≥ 4.5 に分布していたイベントの多くが χ2 /dof < 4.5となり，χ2 /dof のヒストグラムの形が変形している．χ2 /dof のしきい値ごとの到来方向決定精度と，セレクションを通過したイベントの割合を表 6.2 に示す．

図 6.5: 陽子の到来方向再構成の χ2 /dof のヒストグラム．それぞれ (左上)1018.5 eV，(右上)1019.0 eV，(左下)1019.5 eV，(右下)1020.0 eV．

ここでは，到来方向決定精度と，セレクションを通過するイベントの割合から，χ2 /dof=4.5 をクオリティカットのしきい値とする．
図 6.6 に 1018.5 eV，1019.0 eV，1019.5 eV，1020.0 eV の陽子に対する到来方向再構成の pointing
error ヒストグラムを示す．また，pointing error のしきい値ごとの到来方向決定精度と，セレク
ションを通過したイベントの割合を表 6.3 に示す．
図 6.6 と表 6.3 より，極端に大きい pointing error のイベントのみを除去するために，pointing
error=2.5◦ をクオリティカットのしきい値とする．

97

χ2 /dof < 3.0
1.83◦ (83.6%)
χ2 /dof < 3.0
1.51◦ (81.6%)
χ2 /dof < 3.0
1.05◦ (73.5%)
χ2 /dof < 3.0
0.79◦ (66.0%)

1018.5 eV
< 3.5
< 4.0
◦
1.88 (89.4%) 1.93◦ (94.4%)
1019.0 eV
< 3.5
< 4.0
◦
1.54 (88.5%) 1.57◦ (94.3%)
1019.5 eV
< 3.5
< 4.0
1.08◦ (83.4%) 1.11◦ (92.1%)
1020.0 eV
< 3.5
< 4.0
◦
0.82 (78.1%) 0.84◦ (89.5%)< 4.5
1.96◦ (98.8%)

no cut1.98◦< 4.5
1.60◦ (99.2%)

no cut1.62◦< 4.5
1.13◦ (99.3%)

no cut1.13◦< 4.5
0.86◦ (98.8%)

no cut0.86◦

表 6.2: クオリティーカットにおける χ2 /dof のしきい値ごとの到来方向決定精度．括弧内はセレク
ションを通過したイベントの割合を表す．エネルギーごとに上から 1018.5 eV，1019.0 eV，1019.5 eV，1020.0 eV の陽子の解析結果を表す．

図 6.6: 陽子の到来方向再構成の pointing error のヒストグラム．それぞれ (左上)1018.5 eV，(右上)1019.0 eV，(左下)1019.5 eV，(右下)1020.0 eV．

981018.5 eVpointing error< 1.5◦
1.99◦ (1.6%)
pointing error< 1.5◦
1.19◦ (44.1%)
pointing error< 1.5◦
1.11◦ (97.6%)
pointing error< 1.5◦
0.86◦ (99.9%)< 2.0◦< 2.5◦◦
1.71 (75.4%)
1.97◦ (98.8%)
1019.0 eV< 2.0◦< 2.5◦
1.59◦ (97.9%)
1.62◦ (≃ 100.0%)
1019.5 eV< 2.0◦< 2.5◦
1.13◦ (99.9%)
1.13◦ (≃ 100.0%)
1020.0 eV< 2.0◦< 2.5◦
0.86◦ (≃ 100.0%) 0.86◦ (≃ 100.0%)

no cut1.98◦ (≃ 100.0%)
no cut1.62◦ (≃ 100.0%)
no cut1.13◦ (≃ 100.0%)
no cut0.86◦ (≃ 100.0%)

表 6.3: pointing error のしきい値ごとの到来方向決定精度．括弧内はセレクションを通過したイ
ベントの割合を表す．それぞれ上から 1018.5 eV，1019.0 eV，1019.5 eV，1020.0 eV の陽子の解析結
果を表す．

6.5

到来方向及びジオメトリ決定精度の評価

到来方向決定精度のヒストグラムを図 6.7 に示す．ここで，コア位置はシミュレーションでのコア位置を再構成ソフトウェアにも与えている．ここで一次粒子は陽子である．到来方向決定
精度は各一次エネルギーでそれぞれ，1.96◦ (1018.5 eVの時)，1.60◦ (1019.0 eV)，1.13◦ (1019.5 eV)，0.86◦ (1020.0 eV) となった (表 6.4)．
次に，コア位置を粒子数の重心とする手順 1．を行った場合の到来方向決定精度を図 6.8 に示す．
ここで一次粒子は陽子である．到来方向決定精度は各エネルギーでそれぞれ，2.35◦ (1018.5 eVの時)，1.91◦ (1019.0 eV)，1.40◦ (1019.5 eV)，1.00◦ (1020.0 eV) となった (表 6.4)．
同様にコア位置決定精度のヒストグラムを図 6.9 に示す．到来方向決定精度は各エネルギーでそれぞれ，209 m(1018.5 eVの時)，193 m(1019.0 eV)，186 m(1019.5 eV)，185 m(1020.0 eV) となった
(表 6.4)．これよりコア位置決定精度が約 200 m の時に，到来方向決定精度は約 20% 悪化することがわかった．実際に，TA 実験でシャワーサイズと粒子数の横方向分布からコア位置を求める方法
では，コア位置は 50 m 程度の精度で決定することができる．そのため，コア位置をシャワーサイ
ズと粒子数の横方向分布から決定する場合，到来方向決定精度の悪化が単純に 5% 程度になると仮定すると，到来方向決定精度は各エネルギーでそれぞれ 2.06◦ (1018.5 eVの時)，1.68◦ (1019.0 eV)，
1.19◦ (1019.5 eV)，0.90◦ (1020.0 eV) 程度になると予想される．

99

図 6.7: 陽子のシャワー再構成の到来方向決定精度のヒストグラム．コア位置はシミュレーショ
ンでのコア位置を与えている．それぞれ (左上)1018.5 eV，(右上)1019.0 eV，(左下)1019.5 eV，(右下)1020.0 eV．破線は全イベントの 68% が含まれる到来方向決定精度の値である．

コア位置 = シミュレーションでのコア位置

1018.5 eV
1.96◦

1019.0 eV 1019.5 eV 1020.0 eV
1.60◦
1.13◦
0.86◦
コア位置 = 粒子数の重心
18.510
eV 1019.0 eV 1019.5 eV 1020.0 eV
2.35◦
1.91◦
1.40◦
1.00◦
208 m
193 m
186 m
185 m
表 6.4: 到来方向決定精度，及びジオメトリ決定精度を表す．一次粒子は陽子である．上はシミュレーションでのコア位置を再構成ソフトウェアにもそのまま与えた場合の到来方向決定精度を表す．下はコア位置を粒子数の重心とした時の，到来方向決定精度及びコア位置決定精度を表す．

100

図 6.8: 陽子のシャワー再構成の到来方向決定精度のヒストグラム．コア位置は粒子数の重心である．それぞれ (左上)1018.5 eV，(右上)1019.0 eV，(左下)1019.5 eV，(右下)1020.0 eV．破線は全イベ
ントの 68% が含まれる到来方向決定精度の値である．

101

図 6.9: 陽子のシャワーコア位置決定精度のヒストグラム．それぞれ (左上)1018.5 eV，(右上)1019.0 eV，(左下)1019.5 eV，(右下)1020.0 eV．破線は全イベントの 68% が含まれるコア位置決
定精度の値である．

102

第7章

結論

TA 実験における空気シャワー粒子の到来時刻の遅れの関数，及びその標準偏差の関数を，SDを用いたシミュレーションイベント解析 (一次粒子は陽子) より作成した．関数は SD のシャワー
軸からの距離，粒子数，シャワーの天頂角，空気シャワー軸を地面に射影した半直線を基線とした際の SD の方位角の 4 変数で作成した．

1.78−0.25secθ−0.00065φ
R
TD = 2.6 1 +ρ−1.1+0.5secθ−0.0005φ ×10−9 [sec](4.2)
3000[cm]
TSFront = 2.6 1 +R
3000[cm]TSBehind = 2.6 1 +1.49−0.19secθ−0.0008φ

R
3000[cm]ρ−0.45+0.002(−1+secθ)φ ×10−9 [sec]1.38+(0.00079−0.0011secθ)φ(4.4)

ρ−0.3−0.2secθ+0.003(−1+secθ)φ ×10−9 [sec] (4.5)

この関数は一次エネルギー依存性を SD のシャワー軸からの距離と粒子数で表現することで，一次エネルギーが未知の状態でも使用することができる．また，これらの関数は天頂角 60◦ までの
空気シャワーを対象としており，統計量の増加が期待される．さらに，視野が広がることにより，宇宙のより広い領域で宇宙線の到来方向と候補天体現象との相関の探索が可能になる．
作成した関数を使用した空気シャワージオメトリの再構成プログラムを作成し，その再構成精度をシミュレーションイベント解析 (一次粒子は陽子) から評価した．
シミュレーションでのコア位置を再構成ソフトウェアにも与えた場合の到来方向決定精度は，各エネルギーでそれぞれ 1.96◦ (1018.5 eVの時)，1.60◦ (1019.0 eV)，1.13◦ (1019.5 eV)，0.86◦ (1020.0 eV)となった．
また，コア位置を粒子数の重心とした場合のコア位置の決定精度は，各エネルギーでそれぞれ
209 m(1018.5 eVの時)，193 m(1019.0 eV)，186 m(1019.5 eV)，185 m(1020.0 eV) となった．このコア位置を与えた場合の到来方向決定精度は各エネルギーでそれぞれ，2.35◦ (1018.5 eVの時)，1.91◦ (1019.0 eV)，
1.40◦ (1019.5 eV)，1.00◦ (1020.0 eV) となった．これよりコア位置の決定精度が約 200 m の時に，到来方向決定精度は約 20% 悪化することがわかった．実際に，TA 実験でシャワーサイズと粒子数の横方向分布からコア位置を求める方法では，コア位置は 50 m 程度の精度で決定することができる．そのため，コア位置をシャワーサイズと粒子数の横方向分布から決定する場合，到来方向決定精度の悪化が単純に 5% 程度になると仮定すると，到来方向決定精度は各エネルギーでそれぞ
れ 2.06◦ (1018.5 eVの時)，1.68◦ (1019.0 eV)，1.19◦ (1019.5 eV)，0.90◦ (1020.0 eV) 程度になると予想さ
れる．

103

謝辞
本研究に携わる機会を与えて頂き，研究を進めるにあたっては，全てにおいて数多くのご助言，ご指導いただきました大阪市立大学の荻尾彰一准教授に深く感謝致します．また, 幾度となく的確なご助言，ご指導をいただきました大阪市立大学の林嘉夫教授，小島浩司客員教授 (愛知工業大学) に深く感謝致します．
大阪市立大学の櫻井信之特任助教には本研究における全てについて御指導、御助言して頂きましたことを深く感謝致します．
大阪市立大学の藤井俊博博士 (現 東京大学宇宙線研究所)，山崎勝也氏には毎日の研究生活の中で研究に関することのみならず多くの相談に応じて頂き，的確なご助言を頂きましたことを深く感謝致します．
毎日の研究生活の中で支えとなっていただいた大阪市立大学の倉本和幸氏，松宮大輔氏，後藤昂司氏，小林翔悟氏，小林隼矢氏，寺農篤氏，山口泰弘氏に感謝致します．また，同大学の同期として常に良い刺激を与えてくれ，時には相談にのって頂いた米田泰久氏，南平兵衛氏に深く感謝致します．
東京大学宇宙線研究所の野中敏幸助教には SD に関する知識全般や，アメリカにおける研究生
活でお世話になりましたことを感謝致します．
東京大学宇宙線研究所の川田和正特任助教には空気シャワー構造解析や SD の製作、明野観測所での生活でお世話になりましたことを感謝致します。
立命館大学の奥田剛司特任助教には空気シャワー再構成プログラムについてご助言を頂きまし
たことを感謝致します．
東京大学宇宙線研究所の福島正己教授，大木薫氏, 山本邦之氏，芝田達伸特任助教，多米田裕一
郎博士，池田大輔博士，理化学研究所の冨田孝幸博士，京都大学の長滝重博准教授，東京工業大学の北村星爾氏，山梨大学の向井啓兒郎氏，ユタ大学の Tiffany Wong 氏，漢陽大学の Bok Kyun
Shin 氏，CRC 管理人の Dave Barr 氏，Cindy Barr 氏にはユタ現地での生活，並びに観測作業でお世話になりましたことを感謝致します．
東京大学宇宙線研究所の佐川宏行准教授，竹田成宏助教，木戸英治博士，東京工業大学の常定
芳基助教，得能久生特任助教，神奈川大学の有働慈治特別助手には多くのご指導，ご助言を頂きましたことを感謝致します．
大阪市立大学大学運営本部学務企画科の川口幸子氏，東京大学宇宙線研究所の木次敦子氏，岡村宜子氏には数多くの事務処理を行って頂きましたことを感謝致します．
本研究は宇宙線望遠鏡実験なしには成り立ちませんでした．宇宙線望遠鏡実験コラボレーター，
関係研究機関の皆様に深く感謝します．
最後に，研究生活を支えてくれた家族，友人に深く感謝致します．
なお，本研究は日本学術振興会科学研究費補助金 (特別推進研究) H21 ∼ H25 「最高エネルギー
宇宙線で探る極高宇宙現象」，日本学生支援機構 H23 留学生交流支援制度 (ショートステイ・ショー
トビジット) 「数学研究所をハブとする理学の大学院留学生国際研究交流ネットワークプログラ
ム」の支援を受けて行いました．これらの関係機関の皆様に感謝致します．

104

参考文献
[1] J. A. Simpson, Ann. Revs. Nucl & Part. Sci., 33,323(1983).
[2] L. Drury, Rep. Prog. Phys, 46,973(1983).
[3] A. M. Hillas, Astron. Astrophys., 22,425(1984).
[4] C. Thompson, R. Duncan, Monthly Notices of the Royal Astron. Society, 275,255(1995).
[5] P. Blasi, R. Epstein, et al., Astrophys. J. Lett., 533,L123(2000).
[6] A. Hillas, Annual Review of Astronomy and Astrophysics, 22,425(1984).
[7] E. Blodt, P. Ghosh, Mothly Notices of the Royal Astronomical Society, 307,491(1999).
[8] J. Rachen, P. Biermann, Astron. Astrophys., 272,161(1993).
[9] M.J.R. Rees, R. Mészáros, Astrophys. J., 93,430(1994).
[10] H. Kang, D. Ryu, T.W. Jones, Astrophys. J., 456,422(1996).
[11] V. Kuzmin, I. Tkachev, Soviet Journal of Experimental and Theoretical Phys. Lett.,
68,271(1998).
[12] T.J. Weiler, Astropart. Phys., 3,303(1999).
[13] K. Greisen, Phys. Rev. Lett, 16,748(1966).
[14] T. Zatsepin, V.A. Kuzmin, JETP Lett, 4,178(1966).
[15] Antoine Letessier-Selvon, Todor Stanev, Rev. Mod. Phys., 83,907(2011).
[16] Y.Y. Jiang, L. Hou, J. Han, X. Sun, W. Wang, Astrophys. J., 719,459(2010).
[17] The TA Collaboration, The Telescope Array Project Design Report (2000)[18] P.P Kronberg, Rept. Prog. Phys., 57,325(1994).
[19] L.W. Cronin, Nucl. Phys. B., 138,465(2005).
[20] T. Abu-Zayyad, et al., arXiv: 1205.5984v1(2012).
[21] T.K. Gaisser, Cosmic Rays and Particle Physics (1990).
[22] K. Kamata, J. Nishimura, Suppl. Prog. Theor. Phys., 6,93(1958).
[23] K. Greisen, Progress in Cosmic Ray Physics III.(ed.by J.G. Wilson), 27(1956).
105

[24] J. Linsley et al., J. Phys. Soc. Japan., 17 Suppl A-III 91(1962).
[25] K. Greisen, Ann. Rev. Nucl. Sci., 10,63(1960).
[26] J. Linsley, Phys. G: Nucl. Phys., 12,51(1986).
[27] T. Hara, Proc. of 18th ICRC., 11,27(1983).
[28] 一般社団法人日本物理学会, 日本物理学会誌, 第 67 巻(2012).
[29] F. Salamida, et al., Proc. 32nd Int. Conf. Cosmic Rays., 2,145(2011).
[30] D. Ivanov, B.T. Stokes, G.B. Thomson, et al., Proc. 32nd Int. Conf. Cosmic Rays.,2,258(2011).
[31] A.A. Ivanov, et al., Moscow University Physics Bulletin., 65(2010).
[32] R. Abbasi, et al., Phys. Rev. Lett., 100,101101(2008).
[33] U. Micheal, et al., UHECR2012. (2012).
[34] Y. Tameda, et al., Proc. 32nd Int. Conf. Cosmic Rays., HE1.3,1268(2011).
[35] The Pierre Auger Collaboration, Science, 9,938(2007).
[36] 武多昭道, 博士論文 東京大学 (2011).
[37] W.M. Yao, et al., J.Phys.G. and Nucl.Phys., 33(2006).
[38] 山川雄一, 修士論文 東京大学 (2009).
[39] N. Sakurai, et al., Proc. of 30nd ICRC., vol.5,1159(2007).
[40] 宮田孝司, 修士論文 東京理科大学 (2008).
[41] P.K.F. Grieder, Extensive Air Showers, vol.1(2010).
[42] 米田泰久, 修士論文 大阪市立大学 (2013).
[43] H. Tokuno, et al., Nucl. Instr. and Meth. A, 601,364(2009).
[44] Y. Tameda, et al., Nucl. Instr. and Meth. A, 609,227(2009).
[45] 石森理愛, 修士論文 東京工業大学 (2011).
[46] S. Kawana, et al., Nucl. Instr. and Meth. A, 681,68(2012).
[47] T. Tomida, et al., Nucl. Instr. and Meth. A, 654,653(2011).
[48] T. Tomida, et al., International Journal of Nanomanufacturing, Vol.8 No.3,247(2012).
[49] D. Heck, G. Schatz, T. Thouw, J. Knapp, J.N. Capdevielle, CORSIKA: A Monte Carlocode to simulate extensive air showers. Technical Report6019, FZKA (1998).
[50] S. Ostapchenko, QGSJET-II: Towards reliable description of very high energy hadronicinteractions. Nucl. Phys. Proc. Suppl., 151,143(2006).
106

[51] A. Ferrari, P.R. Sala, A. Fasso, J. Ranft, FLUKA: A multi-particle transport code (Program
version 2005). Technical Report 2005-010, CERN (2005).
[52] W.R. Nelson, H. Hirayama, D.W.O. Rogers, The EGS4 code system. Technical Report
0265, SLAC (1985).
[53] S. Agostinelli, et al.,A506,250(2003).

GEANT4:[54] A.M. Hillas, Shower simulation:
52B,29(1997).

A simulation toolkit. Nucl. Instrum. Meth.,

Lessons from MOCCA. Nucl. Phys. Proc. Suppl.,[55] M. Kobal, A thinning method using weight limitation for air-shower simulations. Astropart.Phys., 15,259(2001).
[56] B.T. Stokes, D. Ivanov, et al., Dethinning extensive air shower simulations in CORSIKA.
2011. Submitted to Astroparticle Physics.
[57] D. Ivanov, doctor thesis (2012).
[58] J. Allison, et al., Geant4 developments and applications. IEEE Trans. Nucl. Sci.,
53,270(2006).
[59] Minuit Home Page, http://seal.web.cern.ch/seal/work-packages/mathlibs/minuit/home.html(2013/01 現在)[60] F. James, 1972 CERN computing and data processing school, Pertisau, Austria, (1972)107