後期特別理論演習ロゴ

作成時期・97年10月

画像作成・GATEWAY2000 P5-133

使用言語・VC++




赤ランプ   赤ランプ   赤ランプ



  わが東京大学理学部物理学科では、学部4年次に「特別実験」「特別理論演習」という授業があります。 各教授の元について、実験や研究の訓練をする授業です。大学院進学後の専門授業と同じ内容を取っても、 全く無関係な内容を選択しても、いっこうにかまいません。私は現在、数値計算をメインに 扱う研究室にいて、毎週一回の結果報告の時にレポートを提出しています。 レポート群の前半はグラフィカルで、HP化する意義があると思ったので、ここに転載しました。 後半は数式のみの内容なので、省略しました。

  他の項目と違い、専門知識を要するテーマです。美術や芸術という要素は少ないですが、 学校ではこんなことをやっている、という自己紹介のつもりのページです。


§1 井戸型ポテンシャルと電子間相互作用

10月8日報告分



  量子ホール効果を考えるにあたって、電子間相互作用を実感するため、 簡単な模型でその効果を見てみよう。
  1次元井戸型ポテンシャル中に2個の電子がある時の定常解を求める。 井戸の幅はπで、下図のような形とする。また、スピンはそろっているものとする。

1次元井戸型ポテンシャル

  2個の電子の座標をそれぞれ x 、y とすると、シュレディンガー方程式は定数を省略して

シュレディンガー方程式

である。ただし a はクーロン斥力のパラメーターである。
  a = 0 の時は2次元正方形井戸型ポテンシャルの一体問題の 定常解に帰着し、その解は

Ψ=Σsin(nx)sin(my) ただし n,m=1,2,3,…

であるから、反対称化した解

Ψ=sin x sin 2y - sin 2x sin y

がその解で、図示すると下図のようになる。井戸の中で 二つの電子が離れている状態の確率が大きくなっているのがわかる。

等高線

  実際、観測した時の電子の位置

実際の電子の位置

(注)上の式は誤りです(98.2/3 記)をプロットすると下図のようになる。クーロン斥力が働かなくとも、 電子がフェルミオンであるがゆえの影響が斥力となって表れているのがわかる。

フェルミオンの影響


  さて、次に a≠0 の時を考える。この時は2次元正方形井戸型ポテンシャルで 1体問題を考える時に下図のようなポテンシャルを与えればよい。

ポテンシャル

このポテンシャルに対して、下から6つまでの固有値と固有ベクトルを、 71×71のメッシュ上で数値計算すると

数値計算結果

となる( a = 5.0 )。図において、明るい程 |Ψ| が大きいことを表し、符号はΨの符号である。 図下の数値は固有値である。

  今、問題としている解は左から2番目の解で、これをもとに電子の位置を プロットすると下図のようになる。

クーロン斥力

a≠0 では明らかにクーロン斥力によって2個の電子が反発しているのがわかる。
  今は1次元で考えたが、2次元に拡張し、磁場をかけると、魔法数などの概念が 生まれると思われる。


§2 電子の古典平衡配置と運動の温度依存性

10月22日報告分



  回転放物面ポテンシャル中の電子の、古典平衡配置を求めてみる。 計算方法はいたってオーソドックスな「最急降下法」を使い、電子数3個から10個までの時の、 基底状態といくつかの準安定状態を求めた。

  最急降下法を用いる時には、局所最適解に陥る危険性を十分に考慮しなければならない。 そのためにはパラメーターの初期値(ここでは電子の配置)をうまく取る必要がある。 そこで、(配置を求めたい電子数?1)個の平衡配置に電子を一つ付け加える、 という方針で、3個の時から4個、5個、と順に求めた。新しく加える電子の位置は、 大した根拠はないが、平衡配置のVoronoi図の境界線同士の交点(図の赤丸)、とした。

  こうして得られた結果が下図である。反復回数は各々1万1千回である。 緑字の 10-3 という数字は電子数10個の第3励起状態、ということを意味する。 黒数字は無次元化したエネルギーを電子数で割った値である。

2次元ボロノイ構造

  6-1 における配置のひずみはVoronoi図から明らかである。 また、電子が10個になると平衡配置が突然増える。10-0と10-1のエネルギーの差は 驚くほど小さく、10-2に見られる「わずかなひずみ」や、10-4のような非対称配置など、 電子10個の時は特徴が多彩である。

  電子4つの時の「3角形の中心に1個の電子」という配置は鞍点で、 計算過程ですぐに4-0配置になったので、ここには挙げていない。この種の鞍点の配置は 他に膨大にあると思われる。例えば、一列に並んだ配置など。

  今回私が取った方法は電子を一つ一つ付け加えていくので、 電子数が多い時は法外な時間がかかる。さらに電子数が何10個、何100個となると、 星の数ほど準安定状態が得られるであろう。全てのパターンを調べつくすのは現実的に 不可能である。また、今回の結果が真の答えという保証もない。しかしながら、 上図を見てわかるように、配置の微妙な差というのはそのままエネルギー差が小さいことを意味し、 このことは量子論ではむしろエネルギーの不確定性によって無視されるかもしれない。 電子の古典配置の問題は「どこまで追求すればよいか」ということも議論の対称となるであろう。

  なお、Voronoi図の作図を、コンピュータープログラムで自動作画させることは 極めてやっかいであることを付け加えておく。


  次に、電子の古典的運動に着目してみた。先の図でいう 6-1 配置を初期出発配置(t=0)として、 一つの電子の動径方向の初期速度をいろいろ変えてみて、その後の運動がどのようになるかを見てみた。 初期速度は系の温度に対応する。すなわち、温度を変えるとどうなるか、ということを見てみた。

  電子が6個の時の平衡配置は、先に見たように 6-0 配置と 6-1 配置がある。これらはどちらも 準安定状態なので、6-1 配置から出発した場合、温度が低い時は 6-1 配置から変化しないが、 6-0 配置と 6-1 配置の間にあるであろう壁をこえるエネルギーが温度として与えられると、 よりエネルギーの低い 6-0 配置に移動すると思われる。

  これらの予想をたてて、実際にシミュレーションしたのが下図である。

分子動力学法と対称性

  図は縦軸に回転対称性の目安を与える量を取り、横軸には温度を取り、電子の存在位置を角度で フーリエ変換した量をさらに時間で積分して、その値を全領域でプロットした。 明るいほど積分量が大きいことを表している。すると、矢印Aまでの低温領域では対称性が6の時、 すなわち 6-1 配置にじっとしているが、矢印Aを越えるとしばらくのあいだ対称性5、すなわち 6-0 配置に落ちてあまり動かないのがわかる。さらに温度を上げると 6-0 配置と 6-1 配置を 行ったり来たりするようになり、矢印B付近の温度を越えると0以外の対称性が突如として消える。 これは固体→液体の相転移に対応するものと思われる。

  下側の図は軸はそのままで、積分量を大きく(暗部を明るく)した図である。 対称性が1.5、2.5、…、の時に明るくなっているのはフーリエ変換特有のゴミであると思われるが、 温度が極端に高い時に対称性2.5が若干明るくなっているのが気になる(矢印C)。 単なるゴミか、意味があるのか。

(授業終了後の捕捉)

矢印Aから矢印Bの領域で対称性が5と6の間を行ったり来たりしているのは物理的に誤りとのことである。この原因は電子への初期速度の与え方に問題があるそうだ。

また、上記の対称性について議論をした。小数の対称性は離散的フーリエ変換をしている限り物理的には全く意味がない(これは勘違いによる初等的なミス)。


§3 電子の安定配置の3次元への拡張と最急降下法の考察

10月29日報告分



3次元安定配置

   閉じ込めポテンシャルが3次元調和ポテンシャルの時の、電子の3次元安定配置を求めた。今回も最急降下法を用いたが、 時間の都合から、局所最適解か真の解かの考察がなされていない。従って、数値的な事柄は 記入していない。

   今回は3次元なので、2次元のボロノイ図に相当する図は作図できない。 電子間の「垂直二等分面」同士の交線や交点を作図すれば可能かもしれないが、 極めて複雑なプログラムになると思われるので、代わりに電子の位置を頂点とする多面体を作図した。 しかし、配置のひずみや対称性はこれだけからは明らかではないので、別に幾何学的な考察が 必要であろう。たとえば、13電子の図(最下段)は見た目は正20面体であるが、 ひずんでいるかもしれない(今の段階では不明)。

   3次元では一つの殻に入る電子数が2次元のそれに比べて大きいようである。 今回の結果を見る限り、電子数が12個になって初めて殻の内側に電子が現われた。 ここには示していないが、このあと19個までずっと外側の殻に電子が入っていき、 20個になって初めて内側に2個の電子が入る。電子の数がもっと大きくなれば、 外側の球面上で3角格子の構造を取ることが期待される。

   図の見方について。左の列の青い多面体は、電子の位置を頂点とする多面体を表している。 電子が内外に別れる時は外側の電子についてのみ考えている。 右の二つの図はワイヤーフレーム表示である。ワイヤーは多面体の辺に相当する部分について作図した。 これは交差法で立体視できるようになっている。ぜひ試して頂きたい。

   テクニカルなこと。コンピューターが計算ではじき出す値は、電子の位置座標のみである。 これらは順序がばらばらで、互いの位置関係という概念が含まれていないので、 そこから多面体をコンピューターに自動作成、作画させるのは、ボロノイ図の作画と同様に面倒で やっかいな作業である。なお、多面体の色は、各面の法線と光線ベクトルとのなす角度から 算出している。ちなみに、この小論は紙なので、これ限りだが、コンピューターが稼動する環境では これらの多面体がなめらかにグリグリ上下左右に回転し、見ていて気持ちがよい。 (マセマティカのアニメーションなど比ではない。)

◆   ◆   ◆

   電子の安定配置を求めるにあたって、私は最急降下法を使ったが、 収束速度において最急降下法はニュートン法に比べ段違いに劣る。しかし、 最急降下法はニュートン法に比べてかなり直感的な方法で、しばしば初期配置から 安定配置が推測できる。もし、初期配置と安定配置との間に何かしら明確な相関があれば (結果的には無かったが)、より厳密な解を得られるようになると思い、 初期配置と安定配置との関係を計算してみた。

初期配置と等高線

   右図の等高線はそれぞれ、前のプリントで言う所の4?0配置と8?0配置について、 他の電子が感じるポテンシャルの等高線である (配置された電子のごく近傍は値がオーバーフローしていて、等高線は書かれていない)。

   さて、この配置に新しい電子を付け加えて、最急降下法で収束させることを考える。 4?0配置に電子を付け加える時は、その電子の位置によって5?0配置と5?1配置の どちらかに行き着く。そこで、5?1配置になった時の付け加えた電子の初期位置を赤点で プロットした。同様に8?0配置に対しては9?0配置と9?1配置が考えられるが、 9?0配置に行き着く電子の初期位置を赤点でプロットした。そしてそれを等高線に重ねた。 (なお、赤点のプロットは円形内部で行ったので、図の四隅では計算していない)

   図において、赤点でまだらになっている領域があるが、これは本質的ではなく、誤差に相当する。 最急降下法では勾配ベクトルに対してどの程度電子を進めるかというパラメーターを設定する必要が あるが、このパラメーターは今回0.05とした。これをもっと小さくとれば、 まだらはなくなると思われる。というのは、まだらになる原因は、付け加えた電子と 配置電子との距離が極端に小さいと、勾配ベクトルが非常に大きくなって、 計算の一回のステップで、電子が互いに飛んで離れてしまい、最急降下にならないからである。 このことは、一見数値計算において「悪さ」をしているように思えるが、 逆に利用することを考えて、この領域を集中的に初期配置として採用すれば、 簡単に多くの安定状態を得ることができる可能性がある。

   しかし正直な所、実際には収束速度の点で、ニュートン法に切り替えた方がよいかもしれない。



§4 3次元ボロノイ構造の作図と磁場中の電子の運動

11月5日報告分



3次元ボロノイ構造   3次元ボロノイ構造の作図について。前回3次元ボロノイ構造の作図は極めて複雑と書いたが、 複雑と難解は別物と解釈して、気合のままにプログラムを作成した。電子が3次元中に配置されているとき、 そのボロノイ構造はいくつかの平面で囲まれる。そこで電子同士の垂直2等分面同士の交点などを求めて、 うまく作図したのが右図である。普通に作画すると平面はポリゴンにならず無限大の大きさになるので、 電子を包み込む有限の大きさの球面で全体をカットしてある。

  右下のジャガイモのような図はその左隣の図において、内側の2個の電子の周りに生成される多面体を表している。 眺めている角度は隣の図と少し異なっている。また、光のあたる角度がこれだけ異なって見えるのは ポリゴンの裏表にかかわる問題で、バグではない。しかしながら、3次元ボロノイ構造生成プログラムはとても長く冗長なので、 別のバグが少し残っている。

  数値について。各図左上の数字はエネルギーを電子数で割った値である。小数点以下第15桁が安定するまで反復計算を行い、 14桁目を四捨五入して13桁目まで表示した。右上の2つの数字のうち左の数字は電子の数で、 右の数字はエネルギーの小さい順に並べた状態の種類である。

  今回、電子の配置とエネルギーを計算し直した。方法はあいかわらず最急降下法で、初期配置はランダム配置とした。 この一連の計算をある電子数につき80回行い、その中から基底状態、準安定状態を拾った。 ちなみに、最急降下方法は電子数が多い時収束が悪く、電子数22個を越えると極端に収束が遅くなったので、 今回も20個までとした。次なる課題はやはりニュートン法のマスターであろう。

  これは紙へのプリントアウトなので、これきりだが、実際は上下左右に回転させることができる。 紙面ではわかりにくい対称性なども、回転して様々な角度から眺めてみると直感的に理解しやすい。 もっとも、ボロノイ構造を見ることの目的からして、紙面へのプリントは最初から意味をなさないかもしれない。

◆   ◆   ◆

  放物面ポテンシャルによって2次元平面に閉じ込められた10個の電子の運動を考える。 初期配置として内側に2個外側に8個という基底状態を採用し、全電子にランダムな初期速度を与える。 さらに平面に垂直な方向に磁場をかける。すると電子は電子間相互作用とローレンツ力により、 平面内を複雑に動き始める。この状況において、磁場の強さと初期速度をいろいろ変えてみて、 電子の座標の分散の合計をプロットした。

磁場中の運動と温度   右図において、横軸はある程度時間がたった後の平均運動エネルギー、縦軸は磁場の強さである。 色は分散の大きさに対応し、青い丸は分散が小さく、黒い丸は分散が大きい。図から、 磁場が強いと分散が小さくなることがわかる。このことは磁場が強くなると系は融けにくくなることを意味する。

上図の等高線   上のデータを粗いメッシュごとにまとめて平均を取り、マセマティカで等高線を書くと右図のようになる。 分散が急激に変わる境界線が欲しかったのだが、それは得られなかった。

  ちなみに、磁場中の電子の集団運動は街灯に集まる夜の虫に似ていた。


§5 磁場中2体問題の厳密解の考察

11月19日報告分


問題としている状況は、磁場中の電子の運動方程式で、2体問題に話を限定してみる。 古典平衡配置からのズレを変数に取り、その振幅が小さいとして調和近似を行うと、運動方程式は

運動方程式

ただし

ダイナミカルマトリクスの定義

これを解けばよい。磁場が無い時は対角化された A を用いて 微分方程式 という2階の微分方程式にかけて、ノーマルモードが簡単に求まる。磁場がある時はAとB を同時に対角化しようとすると混乱するので以下のように新しい変数を作って、 変形できる。

変形後の方程式

ただし

yの定義

閉じ込めポテンシャルとして回転放物面ポテンシャルを採用し、ポテンシャルの微分を行って 古典平衡配置のデータ座標データを代入すると

解くべき方程式

と書き下せる。これの固有値を求めると

固有値

となる。対応する固有ベクトルの厳密解は煩雑極まりないので、 固有値 のみについて考えると、任意定数をα、βとして

固有ベクトル

となる。この式において b→0 とすると、ゼロ磁場のノーマルモードの一つ 角速度 に帰着する。磁場が生じると、徐々に円運動に近づき、磁場が強いときは完全な円運動になる。 この意味で、求めた8個の固有値のうち、虚数と符号を除いた値に対応する状態をノーマルモードと 呼ぶことにすると、磁場があっても、2個の粒子には4個のノーマルモード、4個の角速度が 対応していることがわかり、角速度は明確な磁場依存性を持つ。その様子は

磁場と固有値

となる。角速度が磁場に比例するモードとそうでないモードがあることがわかり、 これは今日読んだPHYSICAに出てくる19体の場合の数値計算の結果と似ている。 強磁場の極限については先の8×8の行列においてb のみ計算に効くと考えて、他は全て0として固有値を求め、 そこから角速度を求めると0とbだけになる。 このことは行列の一般性から粒子数によらないので、強磁場極限ではかならず角速度の 分布は2個(0とb)に分かれることになる。



さて、先の8個の固有値に対応する固有ベクトルの具体的な表式はマセマティカで見たところ、かなり煩雑である。 しかし、そのグラフを見ることは容易である。実際にbを変数として固有ベクトルのグラフを書かせたところ、いくつかのことがわかった。ゼロ磁場の時の4個のノーマルモードは別の考察より以下のようになる。

固有ベクトル
(この図の矢印の長さは意味がない。また、実際にはω=0を除いて、矢印の方向に微少振動している。)

ところが、先の8×8の行列を対角化して得られる「ノーマルモード」は、ゼロ磁場の時、

固有ベクトル

となっている。角速度ω=√2の「ノーマルモード」は縮退しているので、このようにも書けるということが自然に出てきた。これらを足したりひいたりすると、その上の図の本来の意味でのノーマルモードになる。

では、磁場をかけると、どうなるか。
?ω=0は系の移動を表すので、変化を受けない。
?ω=√6は先に考察したとおり、徐々に円運動になっていく。
?ω=√2は右回りと左回りとで、角速度が速くなるのと、次第に0になるのとに分かれる。磁場が紙面下からかかっているなら、左回りが磁場に比例して速くなり、右回りは遅くなる。

このようにして、強磁場の極限では、ω=√2の左回りとω=√6の「ノーマルモード」は下図のように発展する。

ノーマルモード

では、角速度が0に発展する残り二つの「ノーマルモード」はどうなるのか?数式の上では以下のような結果が得られる。

固有ベクトル?ただし、Cは任意

これは何かというと、平衡配置から少しならいかようにズレていてもOKということである。しかし、ズレ方が古典配置から遠ざかるズレ方は巨視的に見るとき意味が無く、意味のあるズレとは、結局、系の回転のズレである。

こうして、強磁場極限の「ノーマルモード」とは何か、と言われたら、上の二つの「ノーマルモード」に加えて、一度出てきた下図のような「ノーマルモード」を形式的に付け加えるのがよかろう。

ノーマルモード



§6 3次元多体系でのニュートン法の式と遊び

11月26日報告分


3人のうち、誰が3次元系をやるのか知らないが、その人がニュートン法を使うとき、 役に立つ式を書いておく。

ポテンシャルの式を

ハミルトニアン

とすると、xによる一階微分は

一階微分

となる。他の座標についても同様の形となる。次に2階微分を求めたい。 2階微分の「式の形」は以下の4パターンのみである。

同じ座標で2回微分する時:

    添え字が同じ時: 2階微分

    nとmが等しくないの時: 2階微分

異なる座標で2回微分する時:

    添え字が同じ時: 2階微分

    nとmが等しくないの時: 2階微分

これらを丁寧にコーディングしてdynamical matrixを作ればよい。分母は皆共通なので、 それなりの最適化ができる。2次元については既にプログラムしている。


◆   ◆   ◆


 授業とはまったく関係ない所で、私は趣味としてコンピューターシミュレーションを やることがある。浅瀬の「波」のシミュレーションがしたくて、液体の振る舞いを まったく勝手にモデル立てて、計算してみた。これはこの授業とは無関係なため、 ここには挙げないつもりだったが、どこか電子の液体と似ていると思うので、参考として、 ここにその途中結果を載せる。ただし、昨日の夜に始めたことなので、 十分な考察はされていない。

 粒子数は1万個を考える。働く力は、重力と、空気抵抗と、ごく近くの粒子との クーロン型の斥力である。初期条件として、空中に初期速度0の四角状の塊を考える。この後はニュートンの 運動方程式を解いていく。

 60ステップごとの計算結果を次に載せる。最初はランダムに配置していた塊が、 60ステップ後にはほぼ等間隔に配置されているのに気付く。また、底にたまるにつれて 粒子の密度の差が見られる。

 波のような現象は見られず、まだ液体とは言えない。ほとんど砂だ。

まるで砂




赤ランプ   赤ランプ   赤ランプ



   リンク:[ メニュー ] [ 入口 ] [ 作者 ]
   メール:Morikawa_Hiroshi@yahoo.co.jp