ビームシミュレーション法について

 

 粒子加速器は、各種磁石(双極磁石、四極磁石、六極磁石、 ソレノイドなど) ・高周波空胴・挿入光源等々、様々な線形および非線形要素から 構成されており、その中で運動する荷電粒子の軌道を正確に計算することは極めて 難しい問題である。例えば、リング状の加速器に蓄積されたビームが長時間安定で いられるかどうかを評価しようと思ったら、そのリングを天文学的回数周回した後 の各粒子の位置座標を知らねばならない。(実際、粒子がほぼ光速度で走っている 場合、仮にリングの直径が1キロメートルあったとしても、その粒子はたった1時間 の間にリングを数億回周回する。)粒子ビームの物理的挙動をモデルを使って理論的 に解析することはよくあるが、しかしながら加速器の最終的なデザインの決定は通常 高速のコンピュータを用いたシミュレーション計算の結果に基づいて行われる。課題 演習A4では、ビームシミュレーション法の基礎の基礎についても学習する。

 

I. 調和振動子

 

 まず、調和振動子の運動を考えてみよう。運動方程式は、ωを振動子の角周波数として、

 

(1)


 のように書ける。この問題の解は既にわかっており、

 

(2)


 で与えられる。ここでc1およびc2は積分定数である。

 さて、運動方程式(1)を二つの一階の微分方程式に分離してみよう:

 

(3)

 pは明らかにx方向の運動量成分に対応している。t=0における座標x(0)と運動量p(0)を用いて解(2)を書き直すと、

 

(4)

が得られる。よって、(3)の第一式を使えば、

 

(5)

であることもわかる。このように、初期座標x(0)と初期運動量p(0)さえ知ることができれば、その後の任意の時間における粒子の座標および運動量を予言することが可能である。

 

問題1 微分方程式(1)を解け。

問題2 式(4)が実際に方程式(1)を満足していることを確かめよ。

 

II. 調和振動子運動のシミュレーション

 

 加速器中の粒子運動は、最も単純には、調和振動子によって表現することができる。前節で、調和振動子を扱ったのはそのためである。解析解は既に式(4)および(5)で与えられているが、今度はコンピュータを用いて数値的に運動方程式を積分してみよう。この場合も、式(1)よりもむしろ式(3)を出発点に採る方が都合がよい。

 まず、高校の教科書に書いてあるように、

 

              (6)

であることに注意する。Δtが十分に小さいならば、上式は単に

 

(7)

のように差分化することができる。ここで、xnはある時刻における粒子の座標で、またxn+1はその時刻から微小時間Δtだけ経過した後の座標値を表している。Δtが無限小となる極限において式(7)は厳密であるが、コンピュータ上では有限の数値しか取り扱えないので、シミュレーションの結果には必ず幾分かの誤差が含まれている。Δtを小さくとればとるほど誤差は減るが、その分計算時間が増大してしまう。

 式(7)のような差分化を式(3)に適用すると、

 

(8)

となる。さらに書き直すと、

 

(9)

を得る。したがって、座標と運動量の初期値x0およびp0が与えられれば、式(9)より、その後の座標xnと運動量pn (n=1,2,3,....) がΔtの時間間隔で近似的に計算できる。

 

問題3 差分方程式(9)をプログラム化せよ。

 

III. ハミルトン形式

 

 以上で取り扱ってきた振動子の運動は、ハミルトン関数

 

   (10)

から導かれる。ハミルトン関数Hに支配された粒子の軌道は、いわゆる正準運動方程式

 

(11)

 
から求めることができる。このように、ハミルトン関数は保存系における粒子運動を解析するための基礎であると言える。事実、加速器中の荷電粒子運動は(いわゆる放射減衰のようなある種の摩擦力を無視すれば)通常ハミルトン関数によって記述することが可能である。

 今、式(10)を多少一般化してみる。即ち、V(x)を座標xの任意関数として

 

(12)

のようなハミルトン関数を考える。(V(x)はポテンシャル関数と呼ばれる。)ここまで考えてきた調和振動子は、ポテンシャル関数V(x)がxの二次関数で表される場合に対応している。調和振動子よりもやや複雑な例として

 

(13)

を採ると、今度は振り子の運動を表すことができる。もし、xが常に小さい値しか採らないならば、のように近似できるので、振り子の運動は調和振動子の運動に帰着する。

 物理現象を記述する運動方程式は一般に二階の微分方程式と なるため、その解には必ず二つの未定定数が含まれている。したがって、二つの独立 な物理量の初期値をこれらの未定定数と関係づけてやると都合がよい。 第I節や第II節において、座標だけでなく 運動量も同時に考慮してきたのはそのためである。ハミルトン関数は座標と運動量の 関数であり、粒子運動の時間発展を座標と運動量によって記述することは全く道理 にかなっている。加速器物理におけるシミュレーションでは、縦軸を運動量、横軸 を位置座標にとった空間上で各粒子の運動を追跡する場合が多い。このような空間は “位相空間”と呼ばれており、位相空間上で粒子運動を観察すると、様々な物理的 知見が得られる。

 

問題4 ハミルトン関数(10)による運動方程式が式(3)で与えられることを確認せよ。

問題5 |x|<<1の場合、振り子の運動が調和振動子の運動に帰着することを確かめよ。

問題6 振り子の運動と調和振動子の運動の違いについて論ぜよ。

問題7 V(x)=xn/n (nは整数)の場合の粒子の軌道を図示せよ。

 

IV. シンプレクティック数値積分

 

 各粒子がハミルトン関数に基づいて運動していることが あらかじめわかっている場合、ハミルトン関数に支配された運動に特有の物理的性質 が数値積分の過程で損なわれないよう配慮してやるべきである。前述したように、 微分方程式を式(9)のように差分化した時点で、シミュレーション の結果には不可避的に誤差が混入してしまう。この数値積分の誤差は場合によっては どんどん蓄積し、物理的に全く意味のない結果を与えることもあり得るのである。 例えば、先の調和振動子の例では、ハミルトン関数(10)は時間に 陽には依存していないので、当然ハミルトン関数自体が運動の定数でなければ ならない。事実、厳密解(4)と(5)を式 (10)の右辺に代入すると、

 

(14)

となることがわかる。時刻tが任意にとれる ことに留意すると、式(14)はHの初期値が、 その後のあらゆる時間において保存されることを意味している。つまり、 Hは運動の定数である。この事実を念頭に置き、 式(9)より得られたxおよびp の数値積分値を使って各時刻におけるHを計算してみるとよい。 積分間隔Δtを大きく設定すればするほどHの値が急速に発散することを見るで あろう。あるいは、位相空間上で粒子の軌道を追って みるのもよい。 Hは定数なので、式(10)は

 

(15)

と書くこともでき、粒子が位相空間上で楕円軌道を描くことが わかる。ところが、式(9)に基づく数値積分により得られた粒子 の軌道は決して閉じることはない。これらの事実は、保存すべきエネルギーが保存 せず、無制限に増加することを意味しており、物理的に不合理である。このような 物理的不合理を最小限に抑える方法として“シンプレクティック数値積分”という 便利な方法が考案されており、加速器物理の分野では頻繁に利用されている。 一般的には、ルンゲ・クッタ法と呼ばれる数値積分のアルゴリズムがよく知られて いるが、ルンゲ・クッタ法はシンプレクティックではない。したがって、 ハミルトン系の粒子運動にルンゲ・クッタ法を適用することはできる限り避ける のが賢明である。さて、式(3)の数値積分をシンプレクティック に行うことは極めて簡単である。即ち、式(9)のかわりに

 

(16)

のような差分式を考えてやるだけでよい。つまり、まず旧座標 xnおよび旧運動量pn から新座標xn+1を求めておき、次に旧運動量 pnと新座標xn+1から 新運動量pn+1を計算するのである。 ハミルトン関数が式(12)で与えられる、より一般的な場合も 同様に

 

(17)

 
を用いることによって、シンプレクティックな数値積分を実行することが可能である。
 

問題8 差分方程式(9)によるシミュレーション結果が含む誤差について検討せよ。

問題9 差分方程式(9)によるシミュレーション結果を式(16)による結果と比較せよ。

問題10 式(13)においてω2が定数ではなく、xの関数である場合について考察せよ。