「物理学」カテゴリーアーカイブ

1次元時間依存シュレーディンガー方程式を解かねばならぬ

1次元時間依存シュレーディンガー方程式を確実に解くことを目的とします。

本稿の方針は
計算時間が掛かってもいいので、高精度に解く
ことを主眼においています。対象とする問題は、
ポテンシャルに空間的、時間的に不連続な点があっても高精度に解く
事を目的とします。

ポテンシャルが滑らかで、変なことが起こらないのならば、空間や時間を離散化して解く
クランク=ニコルソン法による時間発展が簡単で優秀な方法です。

本稿では導出する過程では\(N\)次元の問題として扱っていますが、プログラムは1次元の問題のみです。

方針


空間に関する積分を適応型数値積分、時間に関する積分を適応刻み幅ルンゲ=クッタ法で行います。
要は空間も時間も離散化しないで解く、ということです。

時間依存しないシュレーディンガー方程式は適応型数値積分であるQUADPACKを使うことで求められます。
詳しくはこちらへ。

解法


原子単位系における時間依存シュレーディンガー方程式は

として書けます。
波動関数を

のように展開します。ここで、基底関数\(\varphi_n(\mathbf{r})\)はハミルトニアン\(H_0\)の固有関数で、以下の固有値問題

の固有値\(E_n\)に属する固有関数です。
この固有関数は以下の通り規格化されているとします。

関数で展開するのはハミルトニアンに含まれる空間に対する二階微分を解析的に行い、消す為です。
式(1)に代入すると、

となり、左から\(\varphi^*_m(\mathbf{r})\)を掛けて全空間で積分すると

という式が得られます。ここで、\(U_{m,n}(t)\)を

と書きました。

この\(U_{m,n}(t)\)は

の性質を満たします。すなわち、ポテンシャル\(V(\mathbf{r},t)\)が実数値関数であれば、\(U_{m,n}(t)\)はエルミート行列となります。行列形式で書くとすれば、

と表せられます(補遺1)。

ここでは行列形式で書いても見やすくなるだけで、今回は使わないので、式(6)を高精度に解くことに集中します。

高精度に解くために

この問題を解くために高精度に求めなければならない重要な部分は

1. 空間に対する積分
\(\displaystyle
\int \varphi^*_m(\mathbf{r})f(\mathbf{r})\varphi_n(\mathbf{r})d\mathbf{r}
\)


2. 時間に対する積分(連立微分方程式)
\(\displaystyle
\frac{d c_m(t;c_1,~\cdots,c_{N})}{dt}=F_m(t;c_1,~\cdots,c_{N})
\)

という箇所を、不連続箇所があっても高精度に求める手法を使えばよい訳です。

空間に関してはQUADPACKによる適応型の数値積分, 時間に関しては刻み幅制御のルンゲ=クッタ法を用いることで対処します。
各々の詳細は以下のページをご覧ください。
最速・高精度の数値積分
ルンゲクッタ法の説明と刻み幅制御

ここまで確実に解けると仮定すると、精度の不安は基底関数の数だけに依存します。
基底関数の数を積めば積むほど精度が高い計算が出来ますが、計算時間が莫大になります。
ここは計算機の性能と折り合いをつけるしかありません。

ここから考える数値計算上の基底関数\(\varphi_n(x)\)は計算区間の端でゼロになるsin基底関数で、区間[x_a,x_b]で定義され、
\(\displaystyle
\varphi_n(x)=\sqrt{\frac{2}{x_b-x_a}}\sin\left(n\pi\frac{x-x_a}{x_b-x_a}\right)
\)

です。これは、固有値問題
\(\displaystyle
-\frac{1}{2}\frac{d^2}{dx^2}\varphi_n(x)=E_n\varphi_n(x),~~E_n = \frac{1}{2}\left(\frac{n\pi}{x_b-x_a}\right)^2
\)

境界条件
\(\varphi_n(x_a)=\varphi_n(x_b)=0\)
の元での解となっています。

さて、定式化は終わりましたので、あとはプログラムして解くだけです。

莫大な計算時間


問題は計算時間が膨大にかかりすぎる事です。

計算時間を評価してみます。

・次元数を\(D\)
・1次元辺りの基底関数の数は\(N\)
・合計\(N^D\)元1次連立微分方程式を解く
・1つの連立微分方程式の右辺を計算するのに\(N^D\)個の積分
・適応刻み幅ルンゲクッタ法を\(N_{t}\)回繰り返す
・適応刻み幅ルンゲクッタ法は時間1ステップ辺り\(5\)回評価
・積分は空間の分点数\(N_{\mathbf{r}}^D\)回に比例
・ポテンシャルが対称行列であることを利用すれば\(1/2\)倍
・CPU個数を\(N_{\mathrm{cpu}}\)とし、完璧に並列計算が行われれば\(1/N_{\mathrm{cpu}}\)倍
・関数を1回呼び出すのに必要な時間を\(t_f\)

であるので、全計算時間(関数が評価される回数に比例する、と考える)は
\(\frac{5}{2N_{\mathrm{cpu}}}\times N^D_{\mathbf{r}}N_{t}N^{2D}t_f\)と求められます。

現実的な量で計算時間を見積もってみます。
10[原子単位]秒の時間発展を区間20[原子単位]メートルで計算することを考えます。

\(D=1,~ N_{\mathbf{r}}=1000,~N_{t}=2000,~N=50,~N_{\mathrm{cpu}}=1\)とすると、\(1.25\times 10^{10}\)回関数を呼び出さなければなりません。

手持ちのノートパソコン(1.9GHz、1コア使用)で、\(2.5\times 10^8\)回の関数\(f(x)=\sin(x)\)の呼び出しに掛かるCPU時間を測ると約7[cpu秒]かかりました。
ということは、関数を1回呼び出すのにかかる時間\(t_f\)は約\(3\times 10^{-8}\)秒です。

よって、\(1.25\times 10^{10}\)回の関数の呼び出しだけで約400秒、すなわち8分かかります。

現実には他の処理も入るので、8分は絶対にかかる、という事です。
仮にCPUが10個あれば40秒、CPU速度が倍程度の物を使えば20秒位になるでしょうか。
まぁまぁ受け入れられる時間になります。

2次元を計算しようと思ったら基底関数が二乗、かつ空間の分点も二乗になるので、CPUが10個,4GHz近くのCPUを使えても\(20\times 1000\times 50=1000000\)秒程度、すなわち11日です。

11日の計算は現実的ではありません。
二乗で効いてくる基底関数の数を減らすか、空間の積分の精度を犠牲にする、CPUを増やす等の対策をしなければなりません。

また、ここで示した時間は最低かかる時間ですので、実際に動かすのであればこれ以上は確実にかかるでしょう。関数の中の関数を呼び出す等の操作があれば倍々に増えていくのです。

計算量を減らす工夫


ポテンシャルが変わらない場合、\(U_{m,n}(t)\)の積分は、前の時刻と値が変わらないのでいちいち計算するのは無駄です。

数値的に”ポテンシャルが変わっていない”という情報を取り出すことが出来れば、この無駄を排除することが出来ます。

”ポテンシャル関数が変わっていない”というのはどうすれば評価できるのでしょうか?
私が考えたのは、Quadpackの自動積分を行う時に評価された引数の点における関数の値が全て一致すれば関数が変わっていない、と判断するという評価方法です。

実際にこれはうまく行きまして、ポテンシャルが一度だけ切り替わる、という問題に対して計算時間が
1085[CPU sec] → 25[CPU sec]
にまで減少しました。約1/50の計算時間になりました。

先ほどの11日で~と言っていたのはこの工夫をし無ければ、ですので、ポテンシャルの変化が殆どない問題に限れば1時間程度で計算が終わる、ということです。

もちろん、ポテンシャルが変化し続ける場合も試しましたが、変化はほぼありません。
83[CPU sec] → 84[CPU sec]
程度になりました。
この場合は2次元の例で出した通り、11日間掛かる推測です。

初期状態の準備


式(6)を解こうとしても初期条件が与えられなければ数値的に解くことは出来ません。
初期条件で良く使われるのは、1)解析的な解を与える場合と 2)系の固有状態を与える場合の2種類があります。

  1. 解析的な解を与える場合

    この場合、時間発展を考える場合、上記の時間発展方法ではsin関数基底に射影して係数\(c_n(t)\)を求めなければなりません。
    時刻\(t=t_a\)で初期状態\(\psi(x,t_a)\)が解析的に与えられた場合、係数\(c_n(t)\)は、

    として与えられます。

  2. 系の固有状態を与える場合

    この場合は時刻\(t_a\)のポテンシャル中の時間依存しないシュレーディンガー方程式
    \(\displaystyle
    \left[H_0 + V(\mathbf{r},t_c)\right]\phi_n(\mathbf{r})=E’_n\phi_n(\mathbf{r})
    \)

    を実際に解かなければなりません。
    そのためには、時間発展に使う基底と同じ基底\(\varphi_n(\mathbf{r})\)で展開し、対角化するのが賢い方法です。
    すなわち、

    として展開し、ハミルトニアンを行列表示にし、対角化して固有値、固有ベクトルを得ればいいのです。

計算の解析


計算結果を解析するにあたり、ある時刻の固有関数系の基底にどの位存在確率があるかという情報が良く使われます。

初期状態で利用した固有状態の存在確率を調べましょう。
今、求めたいのは時刻\(t_c\)のポテンシャル中の固有関数で展開した時の係数\(a_m(t)\)を求める事です。

ここで、\(\phi_n(\mathbf{r})\)は時刻\(t_c\)のポテンシャル中の固有関数であり、

を満たしています。また、\(\varphi_n(\mathbf{r})\)は数値計算を行う際の基底関数であり、

を満たしています。
任意の時刻\(t\)における波動関数を

として表したいので、\(a_n(t)\)について求めていけば、

という関係式を導くことが出来ます。

数値計算プログラム


プログラムは以下のリンクにあります。
http://slpr.sakura.ne.jp/qp/supplement_data/highprecision_1Dtdsetise.tar.gz
そのままのプログラムでは、以下の設定となっています。
プログラムの利用に関してはhttp://slpr.sakura.ne.jp/qp/about/をお読みください。

解凍すると、4つのファイル(main.f90, integrate.f90,anime.plt,anime_project.plt)が入っています。計算パラメータの設定は、main.f90に入っているので各自調整してください。

繰り返しますが、全て原子単位系です。デフォルトでは、
————————————————–
計算時間の範囲\([ta,tb]\) … \([-5, 20]\)秒
出力時間枚数\(Nt\) … \(100\)枚
1ステップ当りの時間積分の精度\(\mathrm{RKeps}\) … \(10^{-5}\)
1ステップ当りの時間の最小刻み幅\(\mathrm{hmin}\) … \(10^{-4}\)
※理想は最小刻み幅は倍精度ならば\(10^{-14}\)程度が良いと思いますが、時間変化するポテンシャルを計算しようとする時、いつまでも最小刻み幅で進む時間がたまにあります。こういう時間はさっさと抜け出たいので、それを考慮して設定できるようにしています。

位置空間の範囲\([xa,xb]\) … \([-10, 10]\)
出力時の分点数\(Nx\) … 200点
時間積分の精度\(\mathrm{QPeps}\) … \(10^{-6}\)

初期状態を時間依存しないシュレーディンガー方程式を解くか、解析解で与えるかiname … “tise”
inameが”tise”の場合、どの時刻\(t_{ini}\)の固有状態を使うか … \(t_a\)

基底関数の数\(N\) … \(30\)個
計算に用いるスレッド数\(N{\mathrm{cpu}}\) … 1つ
空間積分の分点数を強制的に増やすか否かQPforce … 0

初期状態波動関数… 時間依存しないシュレーディンガー方程式の基底状態
ポテンシャル\(
\begin{eqnarray}
V(x,t)=\left\{
\begin{aligned}
0.5 x^2 (t\lt 0)\\
0.1 x^2 (t\ge 0)
\end{aligned}
\right.
\end{eqnarray}\)
————————————————–
で解くプログラムとなっています。

lapackと、必要があればopemMPを使います。

ifort -mkl -openmp integrate.f90 main.f90
./a.out

もしくは、gfortranでlapackにリンクしてintegrate.f90 main.f90を一緒にコンパイル、実行してください。
実行後、以下のファイルを生成します。

  • wf_***.d
    波動関数の各時刻におけるデータです。時間順にwf_0.d(初期状態), wf_1.d, wf_2.d,…という風に指定された等間隔の時刻で出力します。アウトプットのファイル数はNtで指定した数です。
  • project_ratio.d
    各時刻でのsin関数基底への射影した際の式(2)の係数\(c_n(t)\)の絶対値二乗を出力します。最も高い状態が振幅を持ってしまう場合、基底数が足りていないことを示しています。\( |c_n(t)|^2\)の事です。
  • unitarity.d
    ユニタリー性を確認します。要は数値計算の発散とかで全計算区間内の存在確率が変化していないか見るためのものです。\(\sum_{n=1}^N |c_n(t)|^2\)の事です。
  • timerange.d
    時間グリッドを確認します。物理的な意味は無く、どの時間で計算負荷が大きいかを見るためのものです。
  • timeindex_correspondence.d
    時間インデックス(wf.***.dの***の値と、原子単位系の時刻の変換)を出力します。

もしも初期状態を時間依存しないシュレーディンガー方程式を解いて準備していた場合、以下のファイルを追加で出力します。

  • initial_eigenvalue.d
    初期状態のエネルギー固有値を出力します。式(14)の\(E’_n\)です。
  • project_tdse_to_tisebasis.d
    各時刻における初期状態の固有状態に射影した時の係数の絶対値の二乗を出力します。式(13)の係数\(a_n(t)\)の絶対値二乗である\(|a_n(t)|^2\)を出力します。

波動関数をgnuplot上でアニメーションとして見たい場合、同封してあるファイル”anime.plt”をgnuplot上でloadすれば良いです。

gnuplot> load "anime.plt"

また、gifアニメが欲しい場合、

gnuplot> load "anime_project.plt"

とすれば、波動関数と時間依存しないシュレーディンガー方程式を解いた時の基底関数への係数の絶対値二乗を共に出力します。

実行例


デフォルトの設定で計算を行います。
実行すると

$ ./a.out
0 / 100        1.0000000000
10 / 100        1.0000001064
20 / 100        1.0000001082
30 / 100        1.0000012466
40 / 100        1.0000079960
50 / 100        1.0000082264
60 / 100        1.0000114047
70 / 100        1.0000154242
80 / 100        1.0000154993
90 / 100        1.0000193449
     1.948[CPU sec]
$

という文が出力されます。
左の***/100とは、今計算が100枚中***枚まで終わった、ということを示しており、枚数は変数Ntに対応しています。
右の1に近い数字は、全空間の存在確率を表しており、計算が発散しないかどうかを確認できます。これはルンゲクッタ法の要求精度\(RKeps\)に大きく依存します。
最後のCPUsecは計算に要したCPU時間であり、1コアではおおよそ実時間と同一です。なので、この計算は約2秒で終わる、ということです。ちなみにCPUは1.9GHzの物です。

anime_project.pltによってgifアニメを生成すると、以下のgifアニメが得られます。

その他の出力されるデータは以下の通り。

計算例


ポテンシャルによる反射


※この上図の縦軸は数値計算上の基底である\(\varphi_n(x)\)への射影した時の係数の絶対値二乗です。

ラビ周波数


基底状態と第一励起状態の間の周波数でポテンシャルを揺らしています。

鬼畜なポテンシャル


こんな変なポテンシャルでも計算できます、というデモです。

補遺1

もし仮に、\(U_{m,n}(t)\)が時間によって変化しないのであれば、
\(\displaystyle
\frac{\partial \mathbf{x}}{\partial t}=\mathbf{A}\mathbf{x}
\)

の形の微分方程式であるので、指数関数の解を仮定して\(\mathbf{A}\)を対角化、その固有値固有ベクトルを求めて、線形結合で表現すればokです(線形代数II/連立線形微分方程式等にあるのでそちらを参照してください)。

二体衝突後に物体が持ち得る速度

2体が弾性衝突(エネルギー保存)を起こす時を考えます。
目的は、衝突後に2体が持つ速度ベクトルを知ることです。

2体を
物体1… 重さ\(m_1\)
物体2… 重さ\(m_2\)
衝突前の速度ベクトルを\(\mathbf{v}^{(i)}_1, \mathbf{v}^{(i)}_2\)、
衝突後の速度ベクトルを\(\mathbf{v}^{(f)}_1, \mathbf{v}^{(f)}_2\)
と書くと、衝突後に考えられる速度ベクトルは、
\(
\begin{align}
\mathbf{v}^{(f)}_1&=R(\theta)(\mathbf{v}^{(i)}_1-\mathbf{V})+\mathbf{V}\\
\mathbf{v}^{(f)}_2&=-R(\theta)(\mathbf{v}^{(i)}_2-\mathbf{V})+\mathbf{V} \\
\end{align}
\)

と表せられます。ここで、\(\mathbf{V}\)は重心の速度ベクトルで
\(
\begin{align}
\mathbf{V}=\frac{m_1\mathbf{v}^{(i)}_1+m_2\mathbf{v}^{(i)}_2}{m_1+m_2}
\end{align}
\)

を表し、\(R(\theta)\)はユニタリーな回転行列を表します。衝突が起こる際の角度はこの式からは求められません、あくまで、衝突によって起こり得る角度です。
ここで示しているのは、いろんな方向に衝突が起こることを考えた時の表記で、衝突によって変化することが出来る速度ベクトルを表しています。

特に二次元上での衝突を考える場合、回転行列は
\(
\begin{equation}
R(\theta)=
\left(
\begin{array}{rr}
\cos\theta& -\sin\theta\\
\sin\theta& \cos\theta\\
\end{array}
\right)
\end{equation}
\)

と書くことが出来ます。

二体の衝突


ある2つの物体が衝突することを考えます。
下添え字を物体を区別する番号、上添え字を衝突前\(i\) (initial), 衝突後\(f\) (final)として記述すると、
物体1… 重さ\(m_1\), 位置ベクトル\(\mathbf{r}_1\), 速度ベクトル\(\mathbf{v}_1\)
物体2… 重さ\(m_2\), 位置ベクトル\(\mathbf{r}_2\), 速度ベクトル\(\mathbf{v}_2\)
と書けます。但し、ここでは衝突前、後を表す上添え字\(i,f\)は省いています。

また、重心の位置,速度ベクトル\(\mathbf{R},\mathbf{V}\)は
\(
\begin{align}
\mathbf{R}&=\frac{m_1\mathbf{r}_1+m_2\mathbf{r}_2}{m_1+m_2} \\
\mathbf{V}&=\frac{m_1\mathbf{v}_1+m_2\mathbf{v}_2}{m_1+m_2}
\end{align}
\)

と書けます。

続いて重心から見た物体1,2の位置、速度ベクトルを\(‘\)を付けて表しますと、

物体1(重心系)… 重さ\(m_1\), 位置ベクトル\(\mathbf{r}’_1\), 速度ベクトル\(\mathbf{v}’_1\)
物体2(重心系)… 重さ\(m_2\), 位置ベクトル\(\mathbf{r}’_2\), 速度ベクトル\(\mathbf{v}’_2\)
と書くことが出来ます。
ここで、重心系における速度ベクトルは明示的に
\(
\begin{align}
\mathbf{v}’_1&=\mathbf{v}_1-\mathbf{V} \\
\mathbf{v}’_2&=\mathbf{v}_2-\mathbf{V}
\end{align}
\)

と表すことが出来ます。

重心系から見ると重心の位置、速度ベクトルはもちろんゼロとなります。
\(
\begin{align}
\mathbf{R}’&=\frac{m_1\mathbf{r}’_1+m_2\mathbf{r}’_2}{m_1+m_2}=\mathbf{0}\\
\mathbf{V}’&=\frac{m_1\mathbf{v}’_1+m_2\mathbf{v}’_2}{m_1+m_2}=\mathbf{0}
\end{align}
\)

重心系から見た系の全運動量\(\mathbf{p}’\)もゼロにならなければなりません。
このことから、
\(
\mathbf{p}’=m_1\mathbf{v}_1+m_2\mathbf{v}_2=\mathbf{0}
\)

という関係式が成り立ちます。
この式が言っているのは、重心系を考える限り、物体1,2の速度ベクトルは\(\mathbf{v}’_1,\mathbf{v}’_2\)は衝突が起こる起こらないにかかわらず、上式が満たさなければならない、ということです。
繰り返しますが、衝突が起こる起こらないにかかわらず、です。

衝突が起こり、どの角度に飛んでいくか?を与えられた情報からでは決めることは出来ません。
実験結果や、物体の形状を表す別の条件式が必要となります。

ここで出来ることは、起こりうる角度を知ることだけです。
衝突後に起こりうる角度を求めてみましょう。
重心系では系の合計の運動量\(\mathbf{p}’\)がゼロベクトルにならなければなりません。
このことは衝突前後でノルムが変わる衝突は許されないことを示しています。
これは、重心系においてユニタリーな変換が許されることを意味しています。

衝突後に重心系で取りうる角度は重心系で考えられる角度方向全てなので、
ユニタリーな回転行列\(R(\theta)\)を用いて、衝突後に重心系で速度ベクトルが
\(
\begin{align}
\mathbf{p}’_1&\to R(\theta) \mathbf{p}’_1\\
\mathbf{p}’_2&\to -R(\theta) \mathbf{p}’_2
\end{align}
\)

を満たす変化のみ許される、ということです。マイナスの符号は、\(\mathbf{p}’\)がゼロでなければならないという条件からきているものです

実験系での表現に戻せば、運動量ベクトルは
\(
\begin{align}
\mathbf{p}^{f}_1&=m_1R(\theta) (\mathbf{v}^{(i)}_1-\mathbf{V})+m_1\mathbf{V}\\
\mathbf{p}^{f}_2&=-m_2R(\theta) (\mathbf{v}^{(i)}_2-\mathbf{V})+m_2\mathbf{V}
\end{align}
\)

であり、速度ベクトルは
\(
\begin{align}
\mathbf{v}^{f}_1&=R(\theta) (\mathbf{v}^{(i)}_1-\mathbf{V})+\mathbf{V}\\
\mathbf{v}^{f}_2&=-R(\theta) (\mathbf{v}^{(i)}_2-\mathbf{V})+\mathbf{V}
\end{align}
\)

と書き表すことが出来ます。
特に二次元の場合、デカルト座標において回転行列は
\(
\begin{equation}
R(\theta)=
\left(
\begin{array}{rr}
\cos\theta& -\sin\theta\\
\sin\theta& \cos\theta\\
\end{array}
\right)
\end{equation}
\)

と書くことが出来るので、
\(
\begin{equation}
\left(
\begin{array}{r}
v^{(f)}_x\\
v^{(f)}_y
\end{array}
\right)
=
\left(
\begin{array}{rr}
\cos\theta& -\sin\theta\\
\sin\theta& \cos\theta\\
\end{array}
\right)
\left(
\begin{array}{r}
v^{(i)}_x-V_x\\
v^{(i)}_y-V_y
\end{array}
\right)
+
\left(
\begin{array}{r}
V_x\\
V_y
\end{array}
\right)
\end{equation}
\)

と表現することが出来ます。

衝突角度について


衝突角度を簡単に推測してみましょう。
衝突に際し重要なのは、物体と物体が接する点における法線ベクトル\(\mathbf{n}\)だと推測できます。
恐らく、実験系において法線ベクトルを軸として、反転させた入射ベクトル\(-\mathbf{v}^{(i)}_{1,2}\)を\(\pi\)回転させる角度になるでしょう、と予測できます。
これは私の直感なので間違っているかもしれません。

ともあれ、必要な関係式は\(\mathbf{v}^{(i)}_1,\mathbf{v}^{(i)’}_1 \)を軸とした時の、実験系の角度\(\theta_1\)と重心系の角度\(\theta’_1\)との間の関係式です。

衝突後の速度ベクトルは
\(
\begin{align}
\mathbf{v}^{(f)}_1&=R(\theta) (\mathbf{v}^{(i)}_1-\mathbf{V})+\mathbf{V}\\
\mathbf{v}^{(f)}_2&=-R(\theta) (\mathbf{v}^{(i)}_2-\mathbf{V})+\mathbf{V}
\end{align}
\)

で表されていますので、実験系での角度、すなわち\(\mathbf{v}^{(i)}_1\)と\(\mathbf{v}^{(f)}_1\)のなす角度\(\theta_1\)は
\(
\begin{align}
\theta=\text{acos}
\left(
\frac{\mathbf{v}^{(i)}_1\cdot \mathbf{v}^{(f)}_1}{|\mathbf{v}^{(i)}_1||\mathbf{v}^{(f)}_1|}
\right)
\end{align}
\)

となります。
今、弾性散乱を考えているので、必要であれば\(|\mathbf{v}^{(i)}_1|=|\mathbf{v}^{(f)}_1|\)が成り立っています。
なので、
\(
\begin{align}
\theta=\text{acos}
\left(
\frac{\mathbf{v}^{(i)}\cdot [R(\theta) (\mathbf{v}^{(i)}_1-\mathbf{V})+\mathbf{V}]}{|\mathbf{v}^{(i)}_1||\mathbf{v}^{(f)}_1|}
\right)
\end{align}
\)

です。
変形すれば、
\(
|\mathbf{v}^{(i)}_1||\mathbf{v}^{(f)}_1|\cos\theta_1-\mathbf{v}^{(i)}_1\mathbf{V}=\mathbf{v}^{(i)}_1R(\theta’)(\mathbf{v}_1^{i}-\mathbf{V})
\)
という関係式が導けます。
この式は\(\theta\)と\(\theta’\)を結びつける式になっていますので、実験系での角度が分かれば、重心系での角度に焼直せるため、衝突した時の方向も決められるのです(簡単には求められなさそうですが・・・)。

量子コンピュータを使用する素因数分解アルゴリズム

量子コンピュータを使うと素因数分解が高速にできる。

が、実際はどうやって素因数分解をするんでしょう?
そして具体的に何が高速化するんでしょう?これを知るのが本稿の目的です。

  1. 素因数分解アルゴリズム(古典的)
  2. 実際のアルゴリズム(量子的)
    1. 初期状態の準備
    2. レジスタ1の状態に量子フーリエ変換を施す
    3. レジスタ2に、レジスタ1を元にした関数を作用させる
    4. レジスタ2を観測する
    5. レジスタ1を量子フーリエ変換する
    6. レジスタ1の観測
    7. 連分数展開によって周期\(r\)を見つける
  3. 最大公約数を求めるアルゴリズム(古典)

素因数分解アルゴリズム


量子的に高速に出来るのであれば、そのアイデア自体は古典的に考えられるはずです。

素因数分解する古典的アルゴリズムは以下のように記述されます。
\(m\)を素因数分解したい数とし、\(x\)をmと互いに素な整数(\(\text{gcd}(x,m)=1\)を満たす数)と表記します。

ここで、\(\text{gcd}(a,b)\)は、\(a\)と\(b\)の最大公約数(greatest common divisor)を出力する関数です。

  1. \(a(=0,1,\cdots )\)を引数として、関数\(f_m(a)=x^a~ \text{mod}~ m\)を計算
  2. \(f_m(a)\)の周期\(r\)を見つける。実は\(f_m(a)\)は周期\(r\)の周期関数(周期\(r\)は実際に計算してみるまで分からない。数式で表せば、\(f_n(a)=f_n(a+r)\))。
  3. 周期\(r\)が判明したら、\(p=\text{gcd}(x^{r/2}+1,m),~q=\text{gcd}(x^{r/2}-1,m)\)を計算
  4. \(p,q\)は\(m\)の素因数となっている

という風に素因数分解を行うことが出来ます。

具体的に\(m=21\)を素因数分解してみましょう。\(m\)と互いに素な整数として、\(x=11\)を選びます。
実際に横軸に\(a\)をとり、関数\(f_m(a)\)を計算すると以下のようになります。

この計算結果を見ますと、周期\(r\)が6だと分かります。

周期が\(6\)が分かったので、\(p=\text{gcd}(x^{r/2}+1,m),~q=\text{gcd}(x^{r/2}-1,m)\)を計算します。すると、
\(
\begin{align}
p&=\text{gcd}(x^{r/2}+1,m)=\text{gcd}(11^{6/2}+1,21)=3 \\
q&=\text{gcd}(x^{r/2}-1,m)=\text{gcd}(11^{6/2}-1,21)=7
\end{align}
\)

と分かります。\(21\)は\(3\times 7\)と書けるので、素因数分解を行うことが出来ました。

さて、量子的に素因数分解を行う場合、上のアルゴリズムはShorのアルゴリズムと呼ばれます。
上の計算は、

  • 最大公約数を求めるアルゴリズム(Euclidの互除法, 計算量\(O(\text{log} n)\))
  • 周期\(r\)を見つけるアルゴリズム(
    古典的:\(O(\exp[(\text{log}n)^{1/3}(\text{ln}\text{ln}n)^{2/3}])\),
    量子的:\(O((\text{log}n)^{2}(\text{log}\text{log}n)(\text{log}\text{log}\text{log}n))\)

です。最も時間が掛かる部分は、周期\(r\)を見つける部分で、これを見つけるために量子コンピュータを用いるのです。Shorのアルゴリズムは、この部分だけ量子コンピュータを使うのです。

量子的に高速に計算する、という部分は、
\(f_m(a)\)の周期\(r\)を見つける
という部分です。

実際のアルゴリズム


ここで紹介するモデルは、量子ビット(キュービット)を出さないで説明するものです。

キュービットの考えはShorのアルゴリズムの本質ではなく、実現できるか?の部分で重要です。
すなわち、ある状態を指定するために
・1つの系で、非常に高い量子番号に属する状態まで制御する(\(1[つの系]\times N[量子状態]\))
ではなく、
・複数の系で、基底状態と1つの励起状態を制御する(\(Q[つの系]\times 2[量子状態]\))
の方が現実で実現しやすいという考えです。

1. 初期状態の準備



2つの系(2つの記憶素子、レジスタ)を用意し、初期状態として各々の基底状態を表します。
この時点で、2つの系は独立に存在するもので、系同士に相互作用や、もつれ等は存在しません。

それぞれの系を添え字\(1,2\)で表すと、系全体の状態\(|\psi\rangle\)は(直積の記号を省略して)
\(
|\psi\rangle=|0\rangle_1|0\rangle_2
\)

と書き表せます。
ここで、\(|a\rangle\)は量子数\(a\)に属する固有状態を表しています。

図示すると、上の図のようになり、左はレジスタ1, 右はレジスタ2の量子状態を表していて、それぞれ、基底状態であることを表しています。

2. レジスタ1の状態に量子フーリエ変換を施す


レジスタ1に量子フーリエ変換を施すことで、レジスタ1の全状態の確率振幅を等確率にします。
量子フーリエ変換は量子状態に対して施す変換で、

\(
\begin{align}
\mathcal{F}|a\rangle&=\frac{1}{\sqrt{N}}\sum_{b=0}^{N-1}e^{-2\pi i a b/N}|b\rangle \\
\mathcal{F}^{-1}|a\rangle&=\frac{1}{\sqrt{N}}\sum_{a=0}^{N-1}e^{2\pi i a b/N}|a\rangle
\end{align}
\)

で定義します。基底状態\(|0\rangle\)に対して量子フーリエ変換を施すと、考えたい全励起状態を等確率にすることが出来ます。
\(
\begin{align}
\mathcal{F}|0\rangle&=\frac{1}{\sqrt{N}}\sum_{b=0}^{N-1}e^{-2\pi i 0 b/N}|b\rangle \\
&=\frac{1}{\sqrt{N}}\sum_{b=0}^{N-1}|b\rangle
\end{align}
\)

量子フーリエ変換を作用させると、上図のようになります。

3. レジスタ2に、レジスタ1を元にした関数を作用させる


レジスタ2に、レジスタ1を参照して計算した結果をレジスタ2に出力する、という演算を行います。
この操作を行うことで、レジスタ1と2をもつれさせます(もつれ、という言葉が正しいかは分かりません。ですが、言いたいことはレジスタ1と2を結びつける、という操作を行います)。

この演算を行う事が出来る装置があるかはわかりませんが、あると仮定します。この ”あるかどうか分からない演算を行う装置” はオラクル(神託装置)といいます。

このオラクルを用いて、レジスタ1の量子数\(a\)の値を引数として、レジスタ2の量子状態を\(|f_m(a)\rangle\)にしていきます。

この作用を\(\hat{U}_{1\to2}\)という演算子で表すと、レジスタ1の結果をレジスタ2に格納するので、

\(
\begin{align}
\hat{U}_{1\to2} \hat{\mathcal{F}}_1|0\rangle_1|0\rangle_2
&=\frac{1}{\sqrt{N}}\sum_{a=0}^{N-1}\hat{U}_{1\to2}|a\rangle_1|0\rangle_2 \\
&=\frac{1}{\sqrt{N}}\sum_{a=0}^{N-1}|a\rangle_1|x^a~ \text{mod}~ m\rangle_2
\end{align}
\)

4. レジスタ2を観測する


レジスタ2を観測し、レジスタ2の量子数を知ります。
この観測は以下の2つの影響を系に与えます。

・レジスタ2の量子状態を確定させる
・レジスタ1が取り得る量子状態に制限を加える

1回の”観測”の操作はただ1点だけを与えますが、数式上扱いづらいです。
なので、数式では何度も観測を行った時に得られる期待値を数式で表すと分かりやすくなります。
レジスタ2の状態を量子数\(k\)の状態に見出す確率を知るために、左から\(\langle k|_2\)を作用させます。
すると、以下のように式変形をすることが出来ます。

\(
\begin{align}
\langle k|\hat{U}_{1\to 2} \hat{\mathcal{F}}_1|0\rangle_1|0\rangle_2
&=\frac{1}{\sqrt{N}}\sum_{a=0}^{N-1}|a\rangle_1|x^a ~\text{mod}~ m\rangle_2 \\
&\approx\frac{1}{\sqrt{N}}\sum_{a=0}^{N-1}|a\rangle_1\langle k|x^a ~\text{mod}~ m\rangle_2 \\
&=\frac{1}{\sqrt{N}}\sum_{a=0}^{N-1}|a\rangle_1\left(\frac{1}{\sqrt{J}}\sum_{j=0}^J \delta_{a,a_j^{(k)}}\right)
\end{align}
\)
ここで、\(a_j^{(k)}\)は、
\(k=x^{a_{j}^{(k)}}~ \text{mod} ~m,~~(j=0,1,\cdots,J-1)\)を満たす\(J\)個の\(a\)です。

この観測により、レジスタ1の状態を特定の状態を\(|a_j^{(k)}\rangle,~~(j=0,1,\cdots, J-1)\)のみ存在させることになります。レジスタ1のとれる状態は周期的で周期\(r\)をもち、
\(
\displaystyle a_j^{(k)}=a_0^{(k)}+j\cdot r,~~(j=0,1,\cdots, J-1)
\)

という形で表すことが出来ます。

レジスタ1の状態を図示するとこのようになっているはずです。

係数をまとめると分かりやすくなります。
単なる式変形ですが、状態に適当な係数が掛かった状態の重ね合わせで表現されるので、なじみ深い形…だと私は思います。
\(
\begin{align}
\langle k|\hat{U}_{1\to 2} \hat{\mathcal{F}}_1|0\rangle_1|0\rangle_2
\approx=\sum_{a=0}^{N-1}\left(\frac{1}{\sqrt{NJ}}\frac{1}{\sqrt{J}}\sum_{j=0}^{J-1} \delta_{a,a_j^{(k)}}\right)|a\rangle_1
\end{align}
\)

5. レジスタ1を量子フーリエ変換する



レジスタ1の状態に対して量子フーリエ変換を行うことで、レジスタ1の持つ周期を知ることが出来ます。

量子フーリエ変換を行う前は量子の属する状態そのものに周期\(r\)の情報は含まれていませんが、量子フーリエ変換を行った後は量子数そのものに周期\(r\)に関係する値が含まれるようになります。

実際に観測を行うと、

\(
\begin{align}
\hat{\mathcal{F}}^{-1}_1\langle k|\hat{U}_{1\to 2} \hat{\mathcal{F}}_1|0\rangle_1|0\rangle_2
&\approx\sum_{a=0}^{N-1}\left(\frac{1}{\sqrt{NJ}}\frac{1}{\sqrt{J}}\sum_{j=0}^J \delta_{a,a_j^{(k)}}\right)\hat{\mathcal{F}}^{-1}_1|a\rangle_1 \\
&=\sum_{b=0}^{N-1}\left(\frac{1}{N\sqrt{J}}\sum_{j=0}^{J-1} e^{2\pi i a_j^{(k)}b/N}\right)|b\rangle_1 \\
&=\sum_{b=0}^{N-1}\left[\frac{1}{N\sqrt{J}}e^{2\pi i a_0^{(k)}b/N}\sum_{j=0}^{J-1} e^{2\pi i jrb/N}\right]|b\rangle_1
\end{align}
\)

レジスタ1の量子数\(b’\)に見出す確率は、左から\(\langle b’|_1\)を作用させて絶対値の2乗をとると分かります。

まず、左から\(\langle b’|_1\)を作用させると、
\(
\begin{align}
C'(b’)=\langle b’|_1\hat{\mathcal{F}}^{-1}_1\langle k|\hat{U}_{1\to 2} \hat{\mathcal{F}}_1|0\rangle_1|0\rangle_2
&\approx\frac{1}{N\sqrt{J}}e^{2\pi i a_0^{(k)}b’/N}\sum_{j=0}^{J-1} e^{2\pi i jrb’/N}
\end{align}
\)

となります。確率は上記\(C'(b)\)の値の絶対値の2乗で与えられるので(表記を\(b’\to b\)に変更)、
\(
\begin{align}
|C'(b)|^2&=\left|\frac{1}{N\sqrt{J}}e^{2\pi i a_0^{(k)}b/N}\sum_{j=0}^{J-1} e^{2\pi i jrb/N}\right|^2 \\
&=\frac{1}{N^2 J}\left|\sum_{j=0}^{J-1} e^{2\pi i jrb/N}\right|^2
\end{align}
\)
となります。確率密度は
\(
\displaystyle \left|\sum_{j=0}^{J-1} e^{i(2\pi jrb/N)}\right|
\)
に依存します。この形は離散フーリエ変換の場合に似ています。
つまり、この項が大きくなる時の状態がレジスタ1の量子状態として観測されることが期待され、
観測される状態とは、\(j\)が変化しても位相が一致しているとき、すなわち

\(2\pi rb/N=2\pi \times (整数)\)

を満たすような\(r\)の時に存在確率密度が大きくなる、ということになります。
整数をsと表記すると、

\(rb/N= s~\to~ b/N= s/r\)

と表されます。

6. レジスタ1の観測



確率密度分布は\(b=a_0^{(k)}+j\cdot r\)の量子状態を中心に分布しますが、あくまで、確率分布ですので、観測を行った場合は若干ずれた所に見出したりします。

実際にレジスタ1の観測を行った際に得られた量子数を\(A\)とします。すると、関係式

\(\frac{A}{N}\approx \frac{s}{r}\)

が近似的に成立していると考えることが出来ます。
今、左辺は完全に分かっている量ですが、右辺の分子・分母は整数を持つ、ということくらいしかわかりません。これを求めるために連分数展開を用います。

7. 連分数展開によって周期\(r\)を見つける


レジスタ1の観測によって得られた状態の量子数を\(A\)とします。今、
\(\displaystyle
\frac{A}{N}=\frac{s}{r}
\)

ここで、\(A\)は観測された量子数、
\(N\)はレジスタの全量子状態数、
\(s (\lt N)\)は任意の整数、
\(r\)は求めたい周期です。

これから行いたい操作は、\(\frac{A}{N}\)の分子分母を出来るだけ小さい数にする、という作業を行いたいのです。
これは単純に割るだけで済む話ではありません。なぜなら、\(\frac{A}{N}\)は綺麗に割り切れる数ではない可能性があるからです。

例えば、\(k=184, N=243\)という値であった場合、直感的に\(\frac{184}{243}\approx \frac{3}{4}\)だということが分かります。

これを機械的に行うためには連分数展開を利用します。
連分数展開とは、ある実数\(x\)を整数\(a_n\)を用いて以下のように展開することです。
\(\displaystyle
x=a_0+\frac{1}{a_1+\frac{1}{a_2+\frac{1}{\cdots}}}
\)

ここで、\(a_n~(n=0,1\cdots)\)は、
\(
\begin{align}
&a_0=\lfloor x \rfloor,~~b_0=\frac{1}{x-a_0} \\
&a_n=\lfloor b_{n-1} \rfloor,~~b_n=\frac{1}{b_{n-1}-a_n}
\end{align}
\)
と得られます。ここで\(\lfloor x \rfloor\)はガウスの記号で、実数の整数部分を表しています。
\(x\)が有理数であれば、連分数展開は有限の項数で終わります。

連分数展開後、\(x\)の分数による近似は
\(\displaystyle
x\approx\frac{d_n}{r_n}
\)

ここで、
\(
\begin{align}
& d_0=a_0,~d_1=1+a_0a_1,~d_n=a_n d_{n-1}+d_{n-2},\\
& r_0=1,~r_1=a_1,~r_n=a_n r_{n-1}+r_{n-2},
\end{align}
\)

として順々に求める事が出来ます。
また、連分数の打切りに関しては定理

「\(\frac{q}{p}\)を\(\left|\frac{q}{p}-x\right|\lt \frac{1}{2q^2}\)を満たす任意の有理数とする時、\(\frac{q}{p}\)は\(x\)の近似分数である。さらに、その近似分数は\(p,q\)の最大公約数は1である。」

を用いて、条件式\(\left|\frac{q}{p}-x\right|\lt \frac{1}{2q^2}\)から決めます。

連分数展開によって分子・分母を出力するプログラムはこんな感じで実装できます。

最大公約数を求めるアルゴリズム(古典)


最大公約数を求めるアルゴリズムはEuclidのアルゴリズムと呼ばれ、他のアルゴリズムと比べても早く終わります。そのプログラムは以下で実装できます。

function gcd(a,b)
  implicit none
  integer,intent(in)::a,b
  integer::gcd
  !compute the greatest common divisor
  integer::x,r
 
  x=a
  gcd=b
  r=mod(x,gcd)
  do while( r .gt. 0)
     x=gcd
     gcd=r
     r=mod(x,gcd)
  enddo
 
  return
end function gcd

参考文献

Elisa Bäumer, Jan-Grimo Sobez, Stefan Tessarini, Shor’s Algorithm
https://qudev.phys.ethz.ch/content/QSIT15/Shors%20Algorithm.pdf

C. P. ウィリアムズ、S. H. クリアウォータ著、西野哲郎、新井隆、渡邉昇訳「量子コンピューティング」(springer, 2000)

水素様原子の波動関数のコード

非相対論的な水素様原子の電子の波動関数を求めるfortranコードを書きます。

解く問題

電荷\(+Ze\)の原子核の周りに存在する、電荷\(-e\)を持つ電子が満たす時間依存しないシュレーディンガー方程式は、
\(
\displaystyle \left[-\frac{1}{2}\Delta-\frac{Z}{r}\right]\psi({\bf r})=E\psi({\bf r})
\)

です。
\(
\left. \psi({\bf r})\right|_{|{\bf r}|\to \infty}\to 0
\)

の境界条件を満たす正規直交化された解は、
\(
\begin{align}
\psi({\bf r})&=R_{n,l}(r)\cdot Y_{l,m}(\theta,\varphi)\\
&=\sqrt{\left(\frac{2}{n}\right)^3\frac{(n-l-1)!}{2n[(n+l)!]}} e^{-Zr/n} \cdot \left(2\frac{Zr}{n}\right)^l\cdot L_{n-l-1}^{2l+1}\left(2\frac{Zr}{n}\right)\cdot Y_{l,m}(\theta,\varphi)
\end{align}
\)

と表されます。ここで、

  • \(n,l,m\)  主量子数,方位量子数,磁気量子数
  • \(r,\theta,\varphi\)  動径方向,z軸方向の角度,xy平面の角度
  • \(L_n^{(\alpha)}(x)\)  ラゲール陪多項式[1]
  • \(Y_{l,m}(\theta,\varphi)\)  球面調和関数[2]

を意味します。

動径関数\(R_{n,l}(r)\)は微分方程式
\(
\displaystyle \left[-\frac{1}{2}\left(\frac{d^2}{dr^2}+\frac{2}{r}\frac{d}{dr}\right)+\frac{l(l+1)}{2r^2}-\frac{Z}{r}\right]R_{n,l}(r)=E_n R_{n,l}(r)
\)

を満たす正規直交化された関数で、固有値は
\(
\displaystyle E_n=-\frac{1}{2n^2}
\)
です。

波動関数を出力するFortranコードはこれ。

計算結果


3次元の波動関数を表示するには3つの座標と波動関数の値の計4つが必要です。
しかし、4次元空間を分かりやすく描画する方法はありません。

そこで、ここでは\(\mathrm{Re}{\psi(x,0,z,t)}\)のように、\(y=0\)にしてグラフを描画します。
\(z\)軸周りの構造がうまく表現できませんが、良しとします。

\(n=1\)

\(n=2\)

\(n=3\)

”color”と書いてあるのは、波動関数の形がはっきり映るようにカラーバーの範囲を調整したものです。
なので、色が飽和しているところはそのカラーバーの範囲で表現することはできない領域です。

\(n=20\)まで同じような画像を作ったのですが、画像の容量が大きくなりすぎました。
なので、見たい方は以下のリンクから元ファイルをダウンロードしてください。

\(n=1\)  Hatomc_n1.png(48KB)

\(n=2\)  Hatomc_n2.png(223KB)

\(n=3\)  Hatomc_n3.png(612KB)

\(n=4\)  Hatomc_n4.png(1MB)

\(n=5\)  Hatomc_n5.png(2MB)

\(n=6\)  Hatomc_n6.png(3MB)

\(n=7\)  Hatomc_n7.png(4MB)

\(n=8\)  Hatomc_n8.png(3MB)

\(n=9\)  Hatomc_n9.png(6MB)

\(n=10\)  Hatomc_n10.png(6MB)

\(n=11\)  Hatomc_n11.png(9MB)

\(n=12\)  Hatomc_n12.png(9MB)

\(n=13\)  Hatomc_n13.png(10MB)

\(n=14\)  Hatomc_n14.png(14MB)

\(n=15\)  Hatomc_n15.png(15MB)

\(n=16\)  Hatomc_n16.png(16MB)

\(n=17\)  Hatomc_n17.png(18MB)

\(n=18\)  Hatomc_n18.png(19MB)

\(n=19\)  Hatomc_n19.png(21MB)

\(n=20\)  Hatomc_n20.png(21MB)

n=1~n=20までのファイル hatom.zip(164MB)

例えば、n=13, n=20はこんな感じです(荒くしてます)。

特殊関数に関するプログラミングは特殊関数へ。

——————————————-

How to plot in gnuplot?
Just for note for myself.

You need 2 files which named : multi.plt and atommulti.plt

Before ploting, you must make directory “dat” and move every data files to “dat”.

and in gnuplot, type

call "multi.plt" 2

then, you can get figure.

multi.plt

atommulti.plt

参考


[1]ラゲール陪多項式\(L_m^{(\alpha)}(x)\)

  • 定義域

    \(\;\;\;\;
    0 \lt x \lt \infty
    \)

  • 微分方程式

    \(\;\;\;\;
    \displaystyle
    \left[x\frac{d^2}{dx^2}+(\alpha+1-x)\frac{d}{dx}+n\right]L_m^(\alpha)(x)=0
    \)

  • 具体例

    \(\;\;\;\;
    \begin{align}
    & L^{(0)}_0(x)=1 \\
    & L^{(0)}_1(x)=-x+1 \\
    & L^{(0)}_2(x)=\frac{1}{2}(x^2-4x+2) \\
    & L^{(0)}_3(x)=\frac{1}{6}(-x^3+9x^2-18x+6) \\
    & L^{(1)}_0(x)=1 \\
    & L^{(1)}_1(x)=-x+2 \\
    & L^{(1)}_2(x)=\frac{1}{2}\left(x^2-6x+6\right) \\
    & L^{(1)}_3(x)=\frac{1}{6}\left(-x^3+12x^2-36x+24\right)
    \end{align}
    \)

  • 漸化式

    \(\;\;\;\;
    \begin{align}
    L_n^{(\alpha)}(x)&=\binom{n+\alpha}{n}a_0(x) \\
    &\hspace{1em}a_{m-1}(x)=1-\frac{n-m+1}{m(\alpha+m)}\cdot x\cdot a_m(x) \\
    &\hspace{5em}(m=n,n-1,n-2,\cdots,1,\;\;a_n(x)=1)
    \end{align}
    \)
    [3]
    1階微分
    \(\;\;\;\;
    \begin{align}
    \frac{d}{dx}L_n^{(\alpha)}(x)=-L_{n-1}^{(\alpha+1)}(x)
    \end{align}
    \)
    2階微分
    \(\;\;\;\;
    \begin{align}
    \frac{d^2}{dx^2}L_n^{(\alpha)}(x)=L_{n-2}^{(\alpha+2)}(x)
    \end{align}
    \)

  • 直交性

    \(\;\;\;\;
    \displaystyle
    \int_{0}^{\infty} e^{-x}x^{\alpha}L^{(\alpha)}_n(x)L^{(\alpha)}_{n’}(x) dx =\frac{\Gamma(\alpha+n+1)}{n!}\delta_{n,n’},\;\;(\alpha>-1)
    \)

[2]球面調和関数\(Y_{l,m}(\theta,\varphi)\)
ルジャンドル多項式(\(P_l^m(x)\),[3])を用いて
\(
Y_{l,m}(\theta,\varphi)=(-1)^{(m+|m|)/2}\sqrt{\frac{2l+1}{4\pi}\frac{(l-|m|)!}{(l+|m|)!}}P_l^{|m|}(\cos\theta)e^{im\phi}
\)

[3]ルジャンドル陪多項式\(P_l^m(x)\)

  • 定義域

    \(\;\;\;\;
    1 \lt x \lt -1
    \)

  • 微分方程式

    \(\;\;\;\;
    \displaystyle (1-x^2)\frac{d^2}{dx^2}P_l^m(x)-2x\frac{d}{dx}P_l^m(x)+\left[l(l+1)-\frac{m^2}{1-x^2}P_l^m(x)\right]=0
    \)

  • 具体例

    \(\;\;\;\;
    \begin{align}
    & P_0^0(x)=1 \\
    & P_{1}^{-1}(x)=\frac{1}{2}(1-x^2)^{1/2} \\
    & P_{1}^{0}(x)=x \\
    & P_{1}^{1}(x)=-(1-x^2)^{1/2} \\
    & P_{2}^{-2}(x)=\frac{1}{8}(1-x^2) \\
    & P_{2}^{-1}(x)=\frac{1}{2}x(1-x^2)^{1/2} \\
    & P_{2}^{0}(x)=\frac{1}{2}(3x^2-1) \\
    & P_{2}^{1}(x)=-3(1-x^2)^{1/2} \\
    & P_{2}^{2}(x)=3(1-x^2)
    \end{align}
    \)

  • 漸化式

    用いている式
    \(\;\;\;\;
    \begin{align}
    & P_{|m|}^{|m|}(x)=(-1)^{|m|}(2|m|-1)!!(1-x^2)^{\frac{|m|}{2}} \\
    & P_{l}^{-m}(x)=(-1)^m\frac{(l-m)!}{(l+m)!}P_l^m(x) \\
    & P_{|m|}^{|m|+1}(x)=(2|m|+1)xP_{|m|}^{|m|}(x) \\
    & P_{|m|}^{|m|+q}(x)=\left(\frac{2|m|-1}{q}+2\right)xP_{|m|}^{|m|+q-1}(x)-\left(\frac{2|m|-1}{q}+1\right)P_{|m|}^{|m|+q-2}(x)
    \end{align}
    \)[1]
    1階微分
    \(\;\;\;\;
    \begin{align}
    m=\pm 1&\; \\
    &:\hspace{1em} \frac{d}{dx}P_l^m(x)=\frac{lx}{(x^2-1)}P_l^m(x)-\frac{l+m}{(x^2-1)} P_{l-1}^m(x) \\
    m\ne\pm 1\; \\
    &:\hspace{1em} \frac{d}{dx}P_l^m(x)=c_{2}P^{m+2}_{l-1}(x)+c_0 P^{m}_{l-1}(x)+c_{-2} P^{m-2}_{l-1}(x) \\
    &\hspace{1em} c_{2}=\frac{1}{4(m+1)}\\
    &\hspace{1em} c_{0}=\frac{l+m}{2}\left(1+\frac{l}{1-m^2}\right)\\
    &\hspace{1em} c_{-2}=-\frac{(l+m)(l+m-1)(l+m-2)(l-m+1)}{4(m-1)}
    \end{align}
    \)
         ※\(m=\pm 1\)の時、\(x=\pm 1\)で発散します。
    2階微分
    \(\;\;\;\;
    \begin{align}
    m=\pm 1,\pm 3&\; \\
    &:\hspace{1em} \frac{d^2}{dx^2}P_l^m(x)=\frac{1}{1-x^2}\left[2x\frac{dP_l^m(x)}{dx}+\frac{(l+1)(l+m)}{2l+1}\frac{dP_{l-1}^m(x)}{dx}-\frac{l(l-m+1)}{2l+1}\frac{dP_{l+1}^m(x)}{dx}\right] \\
    m\ne\pm 1,\pm 3&\; \\
    &:\hspace{1em} \frac{d^2}{dx^2}P_l^m(x)=c_{2}\frac{dP^{m+2}_{l-1}(x)}{dx}+c_0 \frac{dP^{m}_{l-1}(x)}{dx}+c_{-2} \frac{dP^{m-2}_{l-1}(x)}{dx}
    \end{align}
    \)
         ※\(m=\pm 1,\pm 3\)の時、\(x=\pm 1\)で発散します。

  • 直交性

    \(l\)に対する直交性
    \(\;\;\;\;
    \displaystyle
    \int_{-1}^{1} P_{m}^{l}(x)P_{m}^{l’}(x) dx =\frac{2}{2l+1}\cdot \frac{(l+m)!}{(l-m)!}\delta_{l,l’}
    \)

[3] Abramowitz & Stegun (1965), ISBN 978-0486612720 (Numerical Method, 22.18, Page 789)
[4] Abramowitz & Stegun (1965), ISBN 978-0486612720 (Numerical Method, 22.2, Page 774,775)

時間依存しないシュレーディンガー方程式と変分原理 2/2 (計算例)

計算例


具体的に適当なポテンシャルの形を考えて解いてみましょう。

プログラムは
時間依存しないシュレーディンガー方程式と変分原理 1/2
にあります。

おさらいですが、間接法による変分原理は、

  • 固有値の上限を与える(変分原理によるエネルギー固有値は必ず厳密なエネルギー固有値よりも大きくなる)
  • 基底関数の数を増やすほど精度は上がる(完全系の場合)

ということを想定しています。
そして、時間依存しないシュレーディンガー方程式と変分原理 1/2 (計算手法の説明)
では、これから束縛状態のみを考えるので、sin関数系の場合の変分原理について計算を行っていきました。
上に想定した変分原理の想定は、

正しい境界条件の下で誤差なく行列要素が求められ、誤差なく対角化が行われた

場合にのみ成り立ちます。sin関数系の場合の境界条件は固定端条件、つまり区間[a,b]で
\(
\psi(x=a)=\psi(x=b)=0
\)

ということです。

ここでは、解析的な解が存在する

  1. 井戸型ポテンシャル
  2. 三角量子井戸
  3. 調和ポテンシャル
  4. F.モースによるポテンシャル[3]
  5. \(1/(\cosh^2)\)型ポテンシャル[3]
  6. クーロンポテンシャル(動径方向, \(l=0\))

について比較していきます。

シュレーディンガー方程式は全て原子単位系(\(m=1, \hbar=1\))で考え、
\(
\displaystyle \left[-\frac{1}{2}\frac{d^2}{dx^2}+V(x)\right]\psi(x)=E\psi(x)
\)

の形を考えます。

無限に深い井戸型ポテンシャル


ポテンシャル

\(
V(x)=\left\{
\begin{eqnarray}
0 \hspace{2em}&(& 0\lt x \lt L,)\\
\infty \hspace{2em}&(& x\le 0 \ or\ L\le x)
\end{eqnarray}
\right.
\)

※井戸の区間が\([a,b]\)であっても変数変換\(y=x-a,\;\;L=b-a\)によって上記ポテンシャルの形に変換されます。

固有値

\(
\displaystyle E_n=\frac{(n+1)^2\pi^2}{2 L^2},\;\;(n=0,1,2,\cdots)
\)

固有関数

\(
\psi_n(x)=\left\{
\begin{eqnarray}
\sqrt{\frac{2}{L}}\sin\left(\frac{(n+1)\pi}{L}(x+\frac{L}{2})\right) &(& 0\lt x \lt L,)\\
0 \hspace{6em}&(& x\le 0 \ or\ L\le x)
\end{eqnarray}
\right.
\)

解法

区間\(0\lt x \lt L,\)ではシュレーディンガー方程式の\(V(x)=0\)の解、sin,cosの形をしてなければなりません。
また、\(\psi(x=0)=0,\;\psi(x=L)=0\)では境界条件より、波動関数の値はゼロです。
よって、解はsin,cosの形で、そのゼロ点がちょうど\(x=0,x=L\)にくるような関数となります。

数値計算解と解析解との比較

数値計算の設定:
Gauss-Lobatto求積法の次数(Ngl):12
計算区間(a,b):[0,\(\pi\)]
基底関数の数(N):200

エネルギー固有値

量子数\(n\) 数値計算解 解析解
0 0.5000000000000000 0.5000000000000000
1 2.0000000000000000 2.0000000000000000
2 4.5000000000000000 4.5000000000000000
3 8.0000000000000000 8.0000000000000000
4 12.5000000000000000 12.5000000000000000
10 60.500000000000000 60.500000000000000
20 220.50000000000000 220.50000000000000
30 480.50000000000000 480.50000000000000
40 840.50000000000000 840.50000000000000
50 1300.5000000000000 1300.5000000000000
考察

行列は元々対角です。なので、対角化する前から対角行列です。
誤差は入りようがありません。

三角量子井戸


ポテンシャル

\(
V(x)=\left\{
\begin{eqnarray}
\infty \hspace{2em}&(& x\lt 0)\\
\alpha x \hspace{2em}&(& 0 \lt x)
\end{eqnarray}
\right.
\)

固有値

\(
\displaystyle E_n=-\left(\frac{\alpha^2}{2}\right)^{1/3}a_n
\)

ここで、\(a_n\)はAiry関数[4]のゼロ点です。

固有関数

\(
\displaystyle \psi_n(x)=C\cdot \text{Ai}\left(\left(\frac{2}{\alpha^2}\right)^{1/3}(\alpha x-E_n)\right)
\)

※積分区間がエアリー関数のゼロ点から無限までであり、これは解析的に実行できません。なので規格化定数をあらわに書くことはできません。

解法

区間\(0\lt x \lt \infty,\)ではシュレーディンガー方程式の\(V(x)=\alpha x\)の解、Airy関数の形をしてなければなりません。
また、\(\psi(x=0)=0, \psi(x\to \infty)=0\)の境界条件を満たします。
よって、解はAiry関数の形で、\(x\)が漸近で減衰していく方の解(\(Ai(x),Bi(x)\)のAiの方)で、そのゼロ点がちょうど\(x=0\)にくるような関数となります。

数値計算解と解析解との比較

数値計算の設定:
Gauss-Lobatto求積法の次数(Ngl):12
計算区間(a,b):[0,60]
基底関数の数(N):200
\(\alpha\):1

エネルギー固有値

量子数\(n\) 数値計算解 解析解
0 1.8557572021449626 1.8557570814892386
1 3.2446077578572123 3.2446076240031596
2 4.3816713853293896 4.3816712392861303
3 5.3866139387958034 5.3866137807905003
4 6.3052631766282223 6.3052630065857747
10 10.8669422991860358 10.8669420487522821
20 16.8461591354108187 16.8461586901918032
30 21.8969186983512394 21.8969179157570224
40 26.4182318860569758 26.4182304793452509
50 30.5803380635310234 30.5803354172938597

考察

上記計算では基底状態から励起状態までほとんど精度が変わっていません。
また、変分原理が想定したとおり、全ての数値計算解は、厳密な解よりも高い結果を与えています。
厳密な解よりも高い結果を与えていることから、
行列要素は高精度で計算されているが、基底関数が足らない、もしくは計算区間が足らない
事が考えられます。

精度がなぜ良くないのかを考えましょう。
\(x=0\)の境界条件は\(\sin\)関数で良く表現されます。が、\(x\to\infty\)側は指数関数で減少していきます。
(一次の)指数関数で減衰していくことを\(\sin\)関数系ではうまく表現できないのではないか、と思います。
実際、Airy関数の指数関数での減衰は漸近で
\(
\displaystyle \frac{e^{-\frac{2}{3}x^{3/2}}}{x^{1/4}}
\)

と表されます。
この後、調和ポテンシャルとクーロンポテンシャルの場合も考えますが、
調和ポテンシャル(漸近で\(\exp(-x^2)\)、基底状態の精度14,15桁)

三角量子井戸(漸近でおおよそ\(\exp(-x^{2/3})\)、基底状態の精度7,8桁)

クーロンポテンシャル(漸近で\(\exp(-x)\)、基底状態の精度2,3桁)
の順に精度が悪くなっていきます。
指数関数的な減衰が早い≒無限に深い量子井戸
を意味するはずなので、恐らくこうなのではと思います。

調和ポテンシャル



\(
\displaystyle V(x)=\frac{1}{2}x^2
\)

固有値

\(
\displaystyle E_n=n+\frac{1}{2},\;\;(n=0,1,2,\cdots)
\)

固有関数

\(
\displaystyle \psi_n(x)=\frac{\pi^{-1/4}}{\sqrt{2^n n!}}e^{-\frac{x^2}{2}} H_n\left(x\right)
\)
ここで\(H_n(x)\)はエルミート多項式を表します。

解法

1次元調和振動子の直接解法

数値計算解と解析解との比較

数値計算の設定:
Gauss-Lobatto求積法の次数(Ngl):12
計算区間(a,b):[-30,30]
基底関数の数(N):200

量子数\(n\) 数値計算解 解析解
0 0.5000000000000038 0.5
1 1.5000000000000340 1.5
2 2.4999999999998268 2.5
3 3.4999999999996776 3.5
4 4.4999999999999165
10 10.4999999999999680 10.5
20 20.4999999999997691 20.5
30 30.5000000000065015 30.5
40 40.5000064794500432 40.5
50 50.5426649449273313 50.5
考察

変分原理の典型的な傾向を見せています。というのは高い励起状態になればなるほど精度が悪くなっているためです。
高励起状態になるにつれて精度が悪くなっていくのは基底関数が足らない為でしょう。計算区間に関してはポテンシャルは境界で高いまま(\(V(a)=V(b)=450\))であり、求められたエネルギー固有値の値\(E_{50}\approx 50\)と比べても十分大きいです。なので、計算区間は足りています。

また、変分原理の固有値の上限を与えるに関しては必ずしも満たしているわけでは無いことに気が付きます。これは行列要素を求める際の誤差や丸め誤差に起因するものであると考えられます。なので、原理が破られているのではなく、数値計算上の問題に起因しています。\(n=30,40,50\)では厳密解よりも大きな値を与えています。これは、数値計算上の誤差よりも変分原理の原理的な値の方が大きくなっているためです。

F.モースによるポテンシャル[3]



\(
\displaystyle V(x)=A\left(e^{-2\alpha x}-2e^{-\alpha x}\right)
\)

固有値

\(
\displaystyle E_n=-A\left[1-\frac{\alpha}{\sqrt{2A}}\left(n+\frac{1}{2}\right)\right]^2
\)

ここで、\(n\)は正の整数で、ゼロから始まり、
\(
\displaystyle\frac{\sqrt{2A}}{\alpha}\gt n+\frac{1}{2}
\)

を満足する最大値\(n_{\text{max}}\)まで。

固有関数

\(
\begin{align}
\psi(x)&=e^{-\xi/2} \xi^s w(\xi) \\
& w(\xi)=F(-n,2s+1,\xi) \\
& \xi=\frac{2\sqrt{2A}}{\alpha}e^{-\alpha x} \\
& s=\frac{\sqrt{-2E}}{\alpha}
\end{align}
\)
ここで\(F(a,b,x)\)は合流型超幾何関数です。

この問題では離散スペクトルは有限個だけ存在します。もしも
\(
\displaystyle\frac{\sqrt{2A}}{\alpha}\gt n+\frac{1}{2}
\)

であれば、離散スペクトルは一般に存在しません。

解法

[3]を参照してください

数値計算解と解析解との比較

数値計算の設定:
Gauss-Lobatto求積法の次数(NGL):12
計算区間(a,b):[-3,80]
基底関数の数(N):200
\(A\):5
\(\alpha\):1

エネルギー固有値

量子数\(n\) 数値計算解 解析解
0 3.5437173567132851 -3.5438611699158100
1 1.3814252908344455 -1.3815835097474309
2 0.2185797147047042 -0.2193058495790518
3 -0.0049226123802256
4 0.0081675032299395
10 0.0527783252864115
20 0.2666761210195411
30 0.6404251854103465
40 1.1703593933291645
50 1.8543552927405405
考察

厳密解は3つしか存在しません。この3つの状態に対しては数値計算解の方が厳密解よりも大きくなり、想定通りです。
しかし精度があまり良くありません。これはポテンシャルが\(V(x\to\infty)=0\)であり、波動関数が速やかにゼロに向かわないためだろうと考えられます。
また、基底状態でもあまり良い結果が得られていません。上記ポテンシャルの閉じ込めが強くなく、多くの基底関数が必要となり、数値計算で用いた200個では十分ではないのでしょう。

4つ目の解に注目しましょう。厳密解は存在しないのに数値計算ではあたかも解が存在するように見えます(\(E\le 0\))。また、変分原理は上限を与えるため、4つ目の解が存在すると考えてしまうかもしれません。しかし、これは違います。明らかにするために4つ目の解の波動関数の形(青色の線)を拡大してみますと、

となり、物理的に意味を成していないことが分かります。なぜなら、波動関数は\(x\to\infty\)に向かうにつれて徐々に小さくなっていかなければなりません。もしもこれが正しければ、\(x=80\)で節を持っている解ということになり、これは束縛状態の解としては適しません。

エネルギー固有値だけを考えては求められないものだったということですね。

さらに、束縛状態以上ではすべて連続状態です。それにもかかわらず\(E\gt 0\)でも離散的に固有値が得られているのは基底関数に固定端の境界条件を課しているためです。これらに物理的な意味は無いので評価する際には気を付けましょう。

\(1/(\cosh^2)\)型ポテンシャル[3]



\(
\displaystyle V(x)=-\frac{V_0}{\cosh^2 \alpha x}
\)

固有値

\(
\displaystyle E_n=-\frac{\alpha^2}{8}\left[-(1+2n)+\sqrt{1+\frac{8V_0}{\alpha^2}}\right]^2
\)

で、\(n\)は
\(
\displaystyle n\lt s = \frac{1}{2}\left(-1+\sqrt{1+\frac{8V_0}{\alpha^2}}\right)
\)

から有限個の準位が決まります。

固有関数

\(
\begin{align}
\psi(x)&=(1-\xi^2)^{\xi/2}F\left[\varepsilon-s, \varepsilon+s+1, \varepsilon+1, \frac{1-\xi}{2}\right]\\
& \varepsilon=\frac{\sqrt{-2E}}{\alpha}\\
& \xi=\tanh \alpha x
\end{align}
\)

ここでの\(F(a,b,c,x)\)は超幾何関数です。

解法

[3]を参照してください

数値計算解と解析解との比較

数値計算の設定:
Gauss-Lobatto求積法の次数(NGL):12
計算区間(a,b):[-30,30]
基底関数の数(N):200

エネルギー固有値

量子数\(n\) 数値計算解 解析解
0 3.8286053565505553 -3.885009803959261
1 1.9750293441823026 -1.9750294118777854
2 0.6459145040771800 -0.7050490197963089
3 0.0750649726143545 -0.0750686277148326
4 0.0076704714372625
10 0.1139873400292487
20 0.4849841979196881
30 1.1689850012563472
40 2.1388129318397651
50 3.3888397726738235
考察

エネルギー固有値の上限を確かに与えています。そのほかはおおよそF.モースによるポテンシャルと同じような結果を与えています。

精度に関して、エネルギー固有値の精度は\(n=0,2\)の時、1~2桁程度の精度ですが、\(n=1,3\)の時、は6~7桁一致と著しく精度が上がっています。計算区間を変えたりしたのですが、この傾向は変わりません。
なぜなのか…詳しくは分かりません。偶関数、奇関数の問題だと思います。すなわち、奇関数の場合はなにか\(x\gt 0\)の領域に誤差があっても同じ量の誤差で\(x\lt 0\)で打ち消し合ってくれるのですが、偶関数の場合は打ち消し合いが起きないのではないか、ということです。
ただ、これが行列要素の精度からくるものなのかは詳しく見ていないので分かりません…

クーロンポテンシャル(動径方向, \(l=0\))


ポテンシャル

\(
\displaystyle V(x)=-\frac{1}{x}
\)

※ここでは、一階微分の存在しない形にした時のシュレーディンガー方程式を考えています。

固有値

\(
\displaystyle E_n=-\frac{1}{2n^2},\;\;(n=1,2,\cdots,)
\)

です。

固有関数

\(
\begin{align}
\psi(x)=C \cdot e^{-x/n} x^l L_{n-1}^{(1)}(x)
\end{align}
\)

ここで,\(C\)は規格化定数, \(L_n^k(x)\)はラゲール倍多項式[5]で
\(
\begin{align}
L_0^k(x)&=1 \\
L_1^k(x)&=-x+k+1 \\
L_2^k(x)&=\frac{1}{2}[x^2-2(k+2)x+(k+1)(k+2)] \\
L_3^k(x)&=\frac{1}{6}[-x^3+3(k+3)x^2-3(k+2)(k+3)x+(k+1)(k+2)(k+3)]
\end{align}
\)
です。もちろん規格化定数をあらわに決めることが出来ますが、通常、この問題は3次元のシュレーディンガー方程式を解いた時の動径方向として出てくるので、その時の規格化は
\(
\displaystyle \int_0^\infty x^2 \psi(x)^2 dx =1
\)

で規格化されます。しかし、本稿の規格化は\(\displaystyle \int_0^\infty \psi(x)^2 dx =1\)
で規格化しているので、エネルギー固有値のみの比較を行います。

解法

\(x=0, x\to infty\)の漸近形を考えて、本当の解を
(\(x=0\)の漸近形)×(\(x=\infty\)の漸近形)×(未知関数)
と仮定します。これをシュレーディンガー方程式に代入すると未知関数がラゲール陪多項式だと分かります。

数値計算解と解析解との比較

数値計算の設定:
Gauss-Lobatto求積法の次数(NGL):12
計算区間(a,b):[0,80]
基底関数の数(N):200

エネルギー固有値

量子数\(n\) 数値計算解 解析解
1 0.4970775696080274 -0.5000000000000000
2 -0.1246297940063756 -0.1250000000000000
3 0.0554455545682885 -0.0555555555555556
4 0.0312035408553899 -0.0312500000000000
5 0.0199704293764349 -0.0200000000000000
6 0.0134298208377274 -0.0138888888888889
7 0.0068151982214300 -0.0102040816326531
8 0.0018616612074302 -0.0078125000000000
10 0.0251247065596211 -0.0050000000000000
20 0.2439421094228725 -0.0012500000000000
考察

クーロンポテンシャルの場合、\(x=\infty\)では波動関数は指数関数で減衰していきます。
このため、減衰がゆっくりになってしまい、sin関数でうまく表現できないために精度が落ちているのだと推測できます。
また、計算区間を十分に大きく取ると励起状態の表現はうまくいきますが、反対に基底状態の表現ができなくなります。なぜなら、多くの区間で波動関数はゼロですが、\(x=0\)の近傍だけに大きなピークを持つため、sin関数ではうまく表現できません。その兆候は上の数値計算でも見えてます。基底状態は2桁程度ですが、量子数の増加と共に1桁くらい精度が上がっています。
クーロンの場合、ラゲール関数を基底関数としてとる方が良いでしょう。

参考文献


[1]P. J. Davis, P. Rabinowitz著, 森正武訳, 「計算機による数値積分法」, Gauss-Lobatto求積法 p.88
[2]P. J. Davis, P. Rabinowitz著, 森正武訳, 「計算機による数値積分法」, ゼロ点に挟まれる区間の積分 p.131
[3]ランダウ=リフシッツ,「量子力学1」第3刷, p81-84
[4]Abramowitz and Stegun, HANDBOOK OF MATHEMATICAL FUNCTIONS, p.446,478 http://people.math.sfu.ca/~cbm/aands/page_446.htm

Airy関数のゼロ点は
エアリー関数 Ai, Bi(ゼロ点)-Keisan
より22桁の精度で求めたものを用いています。
ちなみにそれらの値(Airy関数のゼロ点)は、

0 -2.338107410459767038489
1 -4.087949444130970616637
2 -5.52055982809555105913
3 -6.78670809007175899878
4 -7.944133587120853123138

10 -13.6914890352107179283
20 -21.2248299436420969552
30 -27.58838780988244481195
40 -33.28488468190140187962
50 -38.52880830509424882263

です。ゼロ点を0,1,2,…と、今回の計算の都合上0からインデックスを付けています。
もしかすると最後の1~2桁はあってないかもしれませんが、本稿では16桁まであっていればいいので、問題ありません。

[5]
Associated Laguerre Polynomial -wolfram mathworld式(22)-(25)

時間依存しないシュレーディンガー方程式と変分原理 1/2 (計算手法の説明)

時間依存しないシュレーディンガー方程式を数値的に解くための、変分原理について述べていきます。
変分原理はエネルギー固有値の上限を与えます(厳密な解よりも変分原理で求めた解は必ず大きくなる)。
基底関数を適切に選べば高精度の結果が与えられます。ですが、基底関数をどのように選べばよいか?の指標は無く、経験則に頼らざるを得ない状況が出てきます。

ここでは、変分原理によってなぜエネルギー固有値が求まるのかを述べて行きます。

もしも変分原理より正しい解を得たければ、shooting methodと呼ばれる方法も推奨します。
恐らくこちらの方が色々なポテンシャルに対して有効ですが、
・計算速度が遅い
・発散の可能性がある
・変分原理よりも経験値が必要
です。

まとめ


区間[a,b]で定義された時間依存しないシュレーディンガー方程式
\(
\displaystyle \left[-\frac{1}{2}\frac{d^2}{dx^2}+V(x)\right]\psi(x)=E\psi(x)
\)

の解である固有関数を
\(
\displaystyle \psi(x)=\sum_n c_n \varphi_n(x) \ \ …(6)
\)

の形で固定端条件(\(\psi(a)=\psi(b)=0\))の下で探します。
ここで、\(\varphi_n(x)\)は区間\([a,b]\)で正規直交化されたsin関数系
\(
\displaystyle \varphi_n(x)=\sqrt{\frac{2}{b-a}}\sin\left(\frac{n\pi}{b-a}(x-a)\right),\;\;\;(n=1,2,\cdots,N)
\)

です。この時、行列要素\(H_{ij}=\langle\varphi_i|-\frac{1}{2}\frac{d^2}{dx^2}+V(x)|\varphi_j\rangle\)は
\(
\begin{align}
H_{ij}&=K_{ij}+V_{ij} \\
&\hspace{1em}K_{ij}=\frac{1}{2}\left(\frac{j\pi}{b-a}\right)^2\delta_{ij} \\
&\hspace{1em}V_{ij}=\frac{1}{j}\sum_{m=2}^{M-1}\omega_m^{(GL)} \cos\left(\frac{\pi}{2}x_m^{(GL)}\right)
\sum_{n=0}^{j-1}(-1)^n \sin\left(\frac{i\pi}{2j}(x_m^{GL}+2n+1)\right)V(X_n^{(m)}) \\
&\hspace{2em}X_{n}^{(m)}=\frac{h}{2}(x_m^{(GL)}+2n+1)+a,\;\;h=\frac{b-a}{j}
\end{align}
\)
と求められます(\(j\ge i\)を推奨)。ここで、\(\omega_m^{(GL)},\ x_m^{(GL)}\)はM点のGauss-Lobatto求積法の区間[-1,1]における重みと求積点で、
\(
\displaystyle \omega_m=\frac{2}{M(M-1)[P_{M-1}(x_m)]^2},\ \ x_m\ne \pm 1
\)

であり、求積点は関数
\(
\displaystyle \frac{P_{M-1}(x)+P_{M}(x)}{x-1}
\)

のゼロ点です。\(P_{n}(x)\)は\(n\)次Legendre多項式を表します。

  1. 解説
  2. 時間依存から時間依存しないシュレーディンガー方程式へ
  3. 間接法
    1. 行列要素の計算
    2. K_{nm}の計算
    3. V_{nm}の計算
    4. Gauss-Lobatto求積法による積分の計算
    5. ゼロ点に挟まれる区間の積分
  4. 数値計算プログラム
  5. 正規直交系を用いない場合
  6. 研究分野での数値計算法について
  7. 参考文献

計算の具体例は
時間依存しないシュレーディンガー方程式と変分原理 2/2 (計算例)をご覧ください。


スポンサーリンク



解説


対象とする問題を導出します。何はともあれ、時間依存シュレーディンガー方程式から出発します。

時間依存から時間依存しないシュレーディンガー方程式へ


時間依存シュレーディンガー方程式は
\(
\displaystyle i\frac{\partial \Psi(x,t)}{\partial t}=\hat{H}(x)\Psi(x,t)
\)

と記述されます。ここで、簡単ために原子単位系(\(\hbar=1, m=1\))を考えます。
変数分離の形で解を仮定します。
\(
\displaystyle \Psi(x,t)=\psi(x)U(t)
\)

代入して移行すると以下の式を得ます。

\(
\begin{align}
i\psi(x)\frac{\partial U(t)}{\partial t}=U(t)\hat{H}\psi(x) \\
i\frac{1}{U(t)}\frac{\partial U(t)}{\partial t}=\frac{1}{\psi(x)}\hat{H}\psi(x) = E
\end{align}
\)

2番目の式の右辺\(E\)は何らかの定数です。
左辺は時間のみを含む式、真ん中の式は位置のみを含む式なので、いつどんな時でもこの関係式が成り立っていないといけないという条件からこれらはなんらかの定数でなければなりません。その定数を\(E\)と置きました。

すなわち、2つの式
\(
\begin{align}
i\frac{\partial U(t)}{\partial t}=E U(t)\;\;\; \cdots(1)\\
\hat{H}\psi(x) = E\psi(x)\;\;\; \cdots(2)
\end{align}
\)
が導けます。

\(t\)に関する式は\(\hat{H}(x)\)に依らず解けてしまい、
\(U(t)=e^{-i E t}\)(1)
です。位置に関する式は
\(
\hat{H}\psi(x) = E\psi(x)
\)
(2)
を満たさなければなりません。この(2)が時間依存しないシュレーディンガー方程式と呼ばれているものです。
これを数値計算で解いていくことを考えます。

左から\(\psi^*\)を掛けて全空間で積分しますと、
\(
\begin{align}
\int \psi^* \hat{H} \psi dx – E \int \psi^* \psi dx = 0\;\;\; \cdots(3)
\end{align}
\)
という式を得ます。この問題は
条件\(\int \psi^* \psi dx=1\)の下で、\(\int \psi^* \hat{H} \psi dx\)の極値を求めよ、
という問いになります。
上記問題を考えた時、(3)と同一の式になります。

すなわち、原理的にはどんな\(\psi(x)\)を与えても\(E(=E[\psi])\)の値を計算できるわけです。
その無数に考えられる\(\psi\)の中で、特に(3)の解が停留する\(\psi\)、すなわち
\(
\displaystyle \delta \int \psi^* (\hat{H}-E) \psi dx = 0
\)
となる\(\psi\)を探せ、という問題になります。

↑がもしも停留する\(\psi\)が見付かれば、それはシュレーディンガー方程式と同一になるということは以下の通りに示せます。

\(
\displaystyle \delta \int \psi^* (\hat{H}-E) \psi dx =0 \;\;\;\cdots (4)
\)
厳密な\(\psi\)の変分に対して停留、すなわち\(\psi\to \psi+\delta \psi\)、もしくは複素平面で\(\frac{\pi}{2}\)ずれた方向\(\psi\to \psi+i\delta \psi\)に対して0でなければなりません。
(4)に対して2つの方向に対して変分を取ると、
\(
\begin{align}
\int (\delta \psi)^* (\hat{H}-E)\psi dx &+ \int \psi^* (\hat{H}-E)(\delta\psi) dx =0 \\
-i\int (\delta \psi)^* (\hat{H}-E)\psi dx &+ i\int \psi^* (\hat{H}-E)(\delta\psi) dx =0
\end{align}
\)
と2式が導けます。2式は変分をとった方向にかかわらず同時に満たされていなければなりません。同時に成立するためには
\(
(\hat{H}-E)\psi=0,\;\; \text{or}\;\; \psi^*(\hat{H}^{\dagger}-E^*)=0
\)
を満たしている必要があります。これは時間依存しないシュレーディンガー方程式と一致します。

以上をまとめますと、

厳密な\(\psi\)の変分はシュレーディンガー方程式に一致する

近似的に\(\psi\)を変分するとシュレーディンガー方程式の近似的な解に帰着されるだろう

近似的に\(\psi\)を構成し、変分問題を解く

代表的な方法として、
手法1)直接法…\(\psi\)の具体的な形を変分パラメータと共に与える
手法2)間接法…変分パラメータとしてある正規直交基底関数系を採用する
の2つがありますが、数値計算で使われる間接法を説明します。

間接法


間接法による変分原理は、

  • 固有値の上限を与える(変分原理によるエネルギー固有値は必ず厳密なエネルギー固有値よりも大きくなる)
  • 基底関数の数を増やすほど精度は上がる(完全系の場合)

という特性があります。

さて、今解きたい問題は
条件
\(
\displaystyle \int \psi^* \psi dx =1
\)

の下で
\(
\displaystyle \int \psi^* \hat{H} \psi dx
\)

の極値を与える\(\psi\)を探すこと、です。

すなわち
\(
\displaystyle \int \psi^* \hat{H} \psi dx – E \int \psi^* \psi dx =0 \ \ …(5)
\)

を考えて変分します。求めたい波動関数は
\(
\displaystyle \psi(x)=\sum_n c_n \varphi_n(x) \ \ …(6)
\)

と展開できる、と考えます。ここで、\(c_n\)は変分パラメータ、\(\varphi_n(x)\)は区間[a,b]で定義されている正規直交関数系です。正規直交関数系\(\{\varphi_n(x)\}\)は以下の性質を持ちます。
\(
\displaystyle \int_a^b \varphi_n(x) \varphi_m(x)dx =\delta_{nm}
\)

この制約が無くても変分は可能ですが、以降の式が簡略になるので、正規直交関数だけを考えます。

(6)を(5)に代入すると、
\(
\begin{align}
\sum_{n}\sum_{m} c_n^* c_m \left[H_{nm}-E \delta_{nm}\right]=0 \\
\end{align}\ \ …(7)
\)

なります。ここで、
\(
\begin{align}
H_{nm}&=\int_a^b \varphi_n(x) \hat{H} \varphi_m(x)dx
\end{align}
\)

とおきました。(7)式を満たす係数\(\{c_n\}\)の内、自明な解(\(c_n=0\ (for\ all\ n)\))はすぐ見つかります。しかし、この解は\(\psi(x)=0\)であり、条件\(
\displaystyle \int \psi^* \psi dx =1\)を満たしません。なので解としては不適切です。

では非自明な解を見つけましょう。(7)式の解を探すため、極値を探します。今、係数の組を変分パラメータとみなしているので、(変分パラメータで微分したもの)=0を計算します。すなわち
\(
\displaystyle \frac{\partial}{\partial c_n^*} (for\ all\ n)
\)

を(7)式に作用させます。すると、
\(
\begin{align}
\sum_{n}\sum_{m} \left[H_{nm}-E\delta_{nm}\right]c_m=0
\end{align}
\)

となります。\(n,m\)を行列の行、列ととらえれば、

\(
\begin{align}
\left(
\begin{array}{cccc}
H_{11} & H_{12} & \cdots & H_{1N}\\
H_{21} & H_{22} & \cdots & H_{2N}\\
\ddots & \ddots & \ldots & \ddots \\
H_{N1} & H_{N2} & \cdots & H_{NN}
\end{array}
\right)
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
\ddots \\
c_{N}
\end{array}
\right)
=
E
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
\ddots \\
c_{N}
\end{array}
\right)
\end{align}
\)

となるわけです。ということは、左辺の正方行列を対角化すれば右辺のE,すなわち固有値が求まる、という事になります。

行列要素の計算


行列を対角化するためには行列要素\(H_{nm}\)が分かっていないとなりません。
これからこの行列要素を計算していきます。

ハミルトニアンが
\(
\displaystyle \hat{H}=-\frac{1}{2}\frac{d^2 }{d x^2} + V(x)
\)

として表現されている場合、行列要素は
\(
\begin{align}
H_{nm}&=K_{nm}+V_{nm} \\
&\hspace{3em} K_{nm}=-\frac{1}{2}\int_a^b \varphi_n(x)\frac{d^2}{dx^2} \varphi_m(x)dx \\
&\hspace{3em} V_{nm}=\int_a^b \varphi_n(x)V(x) \varphi_m(x)dx
\end{align}
\)
と書き表され、これらを数値計算で求めます。

\(K_{nm}\)の計算


まず、微分を含んだ
\(
\displaystyle K_{nm}=-\frac{1}{2}\int_a^b \varphi_n(x)\frac{d^2}{dx^2} \varphi_m(x)dx
\)

の項を考えます。微分は数値計算で行いません。
というのも正規直交関数系はこちらで用意する知っている関数系なので、求めることが出来てしまうからです。
多くの正規直交関数系には漸化式と微分方程式が存在します。一番簡単な基底としてsin関数系を考えてみましょう。
区間\([a,b]\)で正規直交化されたsin関数系は
\(
\displaystyle \varphi_n(x)=\sqrt{\frac{2}{b-a}}\sin\left(\frac{n\pi}{b-a}(x-a)\right),\;\;\;(n=1,2,\cdots,N)
\)

と書き表せます。sin関数なので、2階微分はもちろん
\(
\displaystyle \frac{d^2}{dx^2}\varphi_n(x)=-\left(\frac{n\pi}{b-a}\right)^2\varphi_n(x)
\)

を満たします。これを積分の中に代入して、
\(
\begin{align}
K_{nm}&=-\frac{1}{2}\int_a^b \varphi_n(x)\frac{d^2}{dx^2} \varphi_m(x)dx \\
&=\frac{1}{2}\left(\frac{n\pi}{b-a}\right)^2\int_a^b \varphi_n(x) \varphi_m(x)dx \\
&=\frac{1}{2}\left(\frac{n\pi}{b-a}\right)^2\delta_nm
\end{align}
\)

となります。積分を計算することなく求める事が出来ました。

より一般には一回微分の存在しない形、もしくは漸化式によって手計算で微分を含む積分を片付けてしまったり、部分積分を用いて何とか計算します。
正規直交系であれば、うまく整理することによって誤差が全く入らず、かつ高速に計算できます。

\(V_{nm}\)の計算


さて、残る問題は
\(
\displaystyle V_{nm}=\int_a^b \varphi_n(x)V(x) \varphi_m(x)dx
\)

の計算です。基底関数の次数が低い場合、この計算は何も問題がありません。
調和ポテンシャル\(V(x)=x^2\)の時、\(n=1,m=2\)の時の被積分関数は

と非常になだらかですが、\(n=5,m=99\)の時の被積分関数は

となります。なので細かい分点をとることになり、非常に積分が難しいです。

離散フーリエ変換によって求められそうな気もしますが、周期境界条件が満たされていないことは明らかなのであきらめましょう。
有限区間で早く振動する関数の積分は難しいです。二重指数関数型数値積分も振動する関数は苦手です。

こういった場合の計算方法として、
ゼロ点の位置が分かっている場合、そこで区間を分割して積分する
という直観的な方法があります[2]。
ゼロ点の位置が分かっていれば、その間の被積分関数の振る舞いはシンプルな形になるので滑らかな関数として考えられます。

また、端点での値がゼロだと分かっている場合、Gauss-Lobatto求積法が非常に有用です。

Gauss-Lobatto求積法による積分の計算


Gauss-Lobatto求積法はGauss-Legendre求積法と同じく被積分関数を多項式で近似します[1]。
Gauss-Legendre求積法と異なる点は端点で分点を持つということです。端点を明示的に書けば、積分は
\(
\displaystyle \int_{-1}^1 f(x)dx\sim \frac{2}{N(N-1)}(f(-1)+f(1))+\sum_{i=2}^{N-1}\omega_i f(x_i)
\)

と書けます。ここで、\(\omega_i,x_i\)はGauss-Lobatto求積法における重み、求積点です。
Gauss-Lobatto求積法の重みは
\(
\displaystyle \omega_i=\frac{2}{N(N-1)[P_{N-1}(x_i)]^2},\ \ x_i\ne \pm 1
\)

であり、求積点は
\(
\displaystyle \frac{P_{N-1}(x)+P_{N}(x)}{x-1}
\)

のゼロ点です。\(P_{n}(x)\)は\(n\)次Legendre多項式を表します。

Gauss-Lobatto求積法は
\(f(\pm 1)\)が既知の値(特に0)を持つ場合
に便利なことが分かります。
関数\(f(\pm 1)\)がゼロであれば、その点を計算する必要性がなくなりますし、見かけ上の特異点であってもこちらで値を指定すれば参照されることはなく、計算できないという状況に陥りません。

ゼロ点に挟まれる区間の積分


ゼロ点に挟まれる区間の積分は、結局以下の式になります[2]。
\(
\displaystyle \varphi_j(x)=\sqrt{\frac{2}{b-a}}\sin\left(\frac{j\pi}{b-a}\right)
\)

とすると、積分は
\(
\begin{align}
&\int_a^b \varphi_i(x) V(x) \varphi_j(x) dx \\
&=\sum_{n=0}^{j-1} \int_{x_n}^{x_{n+1}} \varphi_i(x) V(x) \varphi_j(x)dx \\
&\approx \sum_{n=0}^{j-1}\sum_{m=2}^{M-1}\frac{h}{2} \omega_m^{(GL)}\varphi_i(X_n^{(m)}) V(X_n^{(m)}) \varphi_j(X_n^{(m)}) \\
&=\frac{h}{2}\sum_{m=2}^{M-1}\omega_m^{(GL)}\varphi_j(X_0^{(m)})S_j^{(m)} \\
\end{align}
\)
ここで、
\(
\begin{align}
&x_n=nh+a,\ h=\frac{b-a}{j},\ (n=0,1,\cdots,j) \\
&X_n^{(m)}=\frac{h}{2}x_m^{(GL)}+\frac{x_{n+1}+x_n}{2} \\
&X_{n+1}^{(m)}=X_n^{(m)}+h \\
&S_j^{(m)}=\sum_{n=0}^{j-1}(-1)^n \varphi_i(X_n^{(m)})V(X_n^{(m)})
\end{align}
\)
と置きました。
\(x_m^{(GL)},\ \omega_m^{(GL)}\)は区間\([-1,1]\)でのM点Gauss-Lobatto求積法の求積点と重みです。
計算を可能な限り減らすため、\(\varphi_j(x)\)の周期性を用いています。

数値計算プログラム


行列要素を組み立てた後は適当なルーチンを用いて対角化します。
今回の場合、基底関数をsin関数に選んだため、強制的に固定端条件が課されます。
なので、この基底関数系を用いるといかなる問題も端点\(a,b\)に無限に高いポテンシャルの壁がある、と読み替えることが出来ます。
よって行列は対称行列、エルミートです。故に固有値は必ず実数です。

この問題を解く場合、実対称行列を対角化するルーチン、lapackのdsyevdを用いるのが良いでしょう。
対角化のプログラムに関する詳細は対角化へどうぞ。

以下のプログラムは行列要素を上記手法で作り、対角化し、得られた波動関数を規格化し、出力するものです。

スポンサーリンク

コンパイルと実行は

$ ifort -mkl main.f90
$ ./a.out

などとしてください。lapackを用いますので、それと一緒にコンパイルしてください。

ラグランジュの未定乗数について


ラグランジュの未定乗数は特に物理的な意味を持たない魔法の変数、として置きますが、今回の場合は未定乗数がエネルギー固有値というちゃんとした意味を持つのです。

正規直交系を用いない場合


\(
\begin{align}
\sum_{n}\sum_{m} c_n^* c_m \left[H_{nm}-E I_{nm}\right]=0 \\
\sum_{n}\sum_{m} c_n^* c_m \left[J_{nm} H_{nm}-E \right]=0
\end{align}\ \ …(7)
\)

となります。ここで、
\(
\begin{align}
H_{nm}&=\int_a^b \varphi_n(x) \hat{H} \varphi_m(x)dx\\
J_{nm}&=I_{nm}^{-1}\\
I_{nm}&=\int_a^b \varphi_n(x) \varphi_m(x)dx
\end{align}
\)
と置きました。行列の形では、
\(
\begin{align}
\left(
\begin{array}{cccc}
J_{11}H_{11} & J_{12}H_{12} & \cdots & J_{1N}H_{1N}\\
J_{21}H_{21} & J_{22}H_{22} & \cdots & J_{2N}H_{2N}\\
\ddots & \ddots & \ldots & \ddots \\
J_{N1}H_{N1} & J_{N2}H_{N2} & \cdots & J_{NN}H_{NN}
\end{array}
\right)
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
\ddots \\
c_{N}
\end{array}
\right)
=
E
\left(
\begin{array}{cccc}
1 & 0 & \cdots & 0\\
0 & 1 & \cdots & 0\\
\ddots & \ddots & \ldots & \ddots \\
0 & 0 & \cdots & 1
\end{array}
\right)
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
\ddots \\
c_{N}
\end{array}
\right)
\end{align}
\)
と書けます。

\(
\begin{align}
\sum_{n}\sum_{m} c_n^* c_m \left[H_{nm}-E I_{nm}\right]=0 \\
\sum_{n}\sum_{m} c_n^* c_m \left[I_{nm}^{-1} H_{nm}-E \right]=0
\end{align}
\)

研究分野での数値計算法について


実際の研究の場では対角化はしますが、上記方法によって行列要素を作ることはしません。遅いからです。
特別な基底関数を用いると行列要素\(H_{nm}\)が1~2項、しかも非常に簡単な式に置き換えることが出来ます。精度も速度も上記方法以上です。
この計算方法に関してはこれ以上は述べることはしませんが、こういった方法も存在する、とだけ述べておきます。

参考文献


[1]P. J. Davis, P. Rabinowitz著, 森正武訳, 「計算機による数値積分法」, Gauss-Lobatto求積法 p.88
[2]P. J. Davis, P. Rabinowitz著, 森正武訳, 「計算機による数値積分法」, ゼロ点に挟まれる区間の積分 p.131
[3]ランダウ=リフシッツ,「量子力学1」第3刷, p81-84

水素原子の解析解1/3 -相対座標と変数分離

ボーアモデルにおける水素原子の解析解、すなわち電子1個(電荷-e)と陽子1個(電荷+e)の系でスピンを無視する場合を考えます。

戦略1)
この問題をそのまま解こうとすると電子で3次元、陽子で3次元なので、3+3=6次元の問題になります。
私達は今、クーロンポテンシャルが電子陽子間でしか依存しないことを知っています。
なので、古典力学で知られている通り、相対座標と重心座標を使い分離できそうだ、と考えられます。
やってみましょう。

2体問題における変数分離


2体系を考えます。(参考文献[1])
2つの物体(質量\(m_1, m_2\))には2体間の相対距離だけで決まる、時間に依存しないポテンシャル\(V({\bf r_1}-{\bf r_2})\)で相互作用しています。
よって、この系のハミルトニアン\(H\)は
\(
\displaystyle H=\frac{{\bf p_1}^2}{2m_1}+\frac{{\bf p_2}^2}{2m_2}+V({\bf r_1}-{\bf r_2})
\)

となるわけです。それに対応するシュレディンガー方程式は、
\(
\displaystyle i\hbar\frac{\partial}{\partial t}\Psi({\bf r_1},{\bf r_2},t)=\left[ -\frac{\hbar^2}{2m_1}\nabla^2_1-\frac{\hbar^2}{2m_2}\nabla^2_2 +V({\bf r_1}-{\bf r_2})\right]\Psi({\bf r_1},{\bf r_2},t)
\)

となります。ここで相対座標系
\(
\displaystyle {\bf r}={\bf r_1}-{\bf r_2},\;\;\; {\bf R}=\frac{m_1{\bf r_1}+m_2{\bf r_2}}{m_1+m_2}
\)

をとります。この座標系ではシュレディンガー方程式は以下のように書き直されます。
\(
\displaystyle i\hbar\frac{\partial}{\partial t}\Psi({\bf R},{\bf r},t)=\left[ -\frac{\hbar^2}{2M}\nabla^2_{\bf R}-\frac{\hbar^2}{2\mu}\nabla^2_{\bf r} +V({\bf r})\right]\Psi({\bf R},{\bf r},t)
\)

ここで\(M=m_1+m_2\)は系の全質量、\(\mu=\frac{m_1m_2}{m_1+m_2}\)は換算質量です。
今、ハミルトニアンが\({\bf R}と{\bf r}\)のそれぞれの和で書けています。
このことは変数分離が可能であることを示しています。なので波動関数を
\(
\Psi({\bf R},{\bf r},t)=\Phi({\bf R})e^{-iE_{CM}t/\hbar}\cdot\psi({\bf r})e^{-iEt/\hbar}
\)

として分離し、シュレディンガー方程式に代入、整理すれば\(\Phi({\bf R})\)と\(\psi({\bf r})\)に関する方程式
\(
\displaystyle -\frac{\hbar^2}{2M}\nabla^2_{\bf R}\Phi({\bf R})=E_{CM}\Phi({\bf R}) \;\;\;\; \ldots(A)
\)

\(
\displaystyle \left[-\frac{\hbar^2}{2\mu}\nabla^2_{\bf r} +V({\bf r})\right]\psi({\bf r})=E\psi({\bf r})\;\;\;\; \ldots(B)
\)

が導けます。
式\((A)\)はエネルギー\(E_{CM}\)の質量\(M\)の自由粒子に対するシュレディンガー方程式とみなせ、
式\((B)\)はポテンシャル\(V({\bf r})\)中を質量\(\mu\)を持った粒子に対するシュレディンガー方程式であるとみなせます。
系の全エネルギー\(E_{tot}\)は
\(
E_{tot}=E_{CM}+E
\)

です。
適切な座標変換により、2体問題2つの1体問題に置き換えた。そういう事を今やったのです。

——
水素原子に適応しましょう[2]。
クーロンポテンシャル\(V(r)\)は
\(
\displaystyle V(r)=-\frac{e^2}{4\pi\varepsilon_0 r}
\)

と書けます。ここで、\(r\)は電子陽子間距離を表します。
電子の質量を\(m\), 陽子の質量を\(M\)と記述します。

重心と共に動く座標系で考えましょう。
今、重心は止まっています(重心の運動量\({\bf P}=0\))。
なので、解くべき問題は1体の時間依存しないシュレディンガー方程式だけで、
\(
\displaystyle \left[-\frac{\hbar^2}{2\mu}\nabla^2 -\frac{e^2}{4\pi\varepsilon_0 r}\right]\psi({\bf r})=E\psi({\bf r}) \;\;\;\;\; \ldots(1)
\)

を解けばよい、ということになりました。\(\mu=\frac{mM}{m+M}\)はこの系の換算質量を表します。

戦略2)
6次元の問題から重心の自由度を3つ減らして3次元の問題に落とせたわけですが、3次元の問題を解くのはまだまだ大変です。相互作用ポテンシャルが中心力ポテンシャルであることから、3次元極座標で考えればrに関する項は変数分離できて、少なくともrに関する次元が更に1つ落とせそうです。あわよくば、残りの2次元も変数分離出来たら1次元の問題が3つになって嬉しいのですが、果たしてそううまくいくのでしょうか・・・

3次元極座標を用いた変数分離


ここでは3次元極座標、特に球面座標を使って変数分離を試みましょう(実は球面座標以外でも放物座標系と呼ばれる座標を使っても変数分離できます)

球面座標は、空間を\(r,\theta,\phi\)を用いて指定する座標です。デカルト座標との間には、

\(
\begin{eqnarray}
\left\{
\begin{aligned}
r &=\sqrt{x^2+y^2+z^2} \\
\theta&= \arccos{\frac{z}{r}}\\
\phi&=\tan^{-1}\frac{y}{x}
\end{aligned}
\right.
\;\;\;\;\;\;\;\;
\left\{
\begin{aligned}
x &= r\sin{\theta}\cos{\phi} \\
y &= r\sin{\theta}\sin{\phi} \\
z &= r\cos{\theta}
\end{aligned}
\right.
\end{eqnarray}
\)

の関係があり、体積要素\(d{\bf V}\)は
\(
d{\bf V}=r^2\sin{\theta}dr d\theta d\phi
\)

と表されます。

波動関数\(\psi({\bf r})\)を動径方向について分離します。
\(
\psi({\bf r})=R(r)Y(\theta,\phi)
\)

この座標系でのラプラシアン\(\nabla^2\)は
\(
\displaystyle \nabla^2=\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial}{\partial r}\right)+\frac{1}{r^2\sin{\theta}}\frac{\partial}{\partial \theta}\left(\sin{\theta}\frac{\partial}{\partial \theta}\right)+\frac{1}{r^2\sin^2{\theta}}\frac{\partial^2}{\partial \phi^2}
\)

であるので、これらを式(1)に代入すると、
\(
\begin{align}
\displaystyle \left[-\frac{\hbar^2}{2\mu}\left\{\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial}{\partial r}\right)+\frac{1}{r^2\sin{\theta}}\frac{\partial}{\partial \theta}\left(\sin{\theta}\frac{\partial}{\partial \theta}\right)+\frac{1}{r^2\sin^2{\theta}}\frac{\partial^2}{\partial \phi^2}\right\} -\frac{e^2}{4\pi\varepsilon_0 r}\right]R(r)Y(\theta,\phi)\\
=ER(r)Y(\theta,\phi)
\end{align}
\)

を得ます。


調べた限りではここからの解法は2つに分かれます。
①数学的に直接解く
②物理的意味を持つ角運動量演算子を考えながら解く
です。
どちらでも最後は同じになるので構わないと思います。
もしも深く知りたい場合は①を解いてから②を解くことをお勧めします。
長くなるのでページを移します。

①直接解く場合は後ほど

②物理的意味を踏まえながら解く場合は後ほど

参考文献

[1]Physics of Atoms and Molecules 2nd Edition pp.109~111 by B.H. Bransden C.J. Joachain(2003)
[2]Physics of Atoms and Molecules 2nd Edition pp.147~148 by B.H. Bransden C.J. Joachain(2003)
※参考した本はこれです。

exp(ikx)が進行波であることの証明

\(e^{i{\bf k\cdot r}}\)が\({\bf k}\)の方向に進行すること、\(e^{-i{\bf k\cdot r}}\)が\(-{\bf k}\)の方向に進行することの所以を証明します。

まとめ

波は、実部と虚部を並べた時、虚部がある方向に進行していきます。

証明


一次元自由粒子の時間発展を考えます。
自由粒子の時間依存シュレーディンガー方程式は
\(
\displaystyle i\hbar\frac{\partial}{\partial t}\Psi(x,t)=\left(-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}\right)\Psi(x,t)
\)

と書けます。時間依存しないシュレーディンガー方程式を求めるために、
\(\Psi(x,t)=\psi(x)e^{-i\frac{E}{\hbar}t}\)
と置いて代入すると
\(
\displaystyle -\frac{\hbar^2}{2m}\frac{d^2}{dx^2}\psi(x)=E\psi(x)
\)

となります。
上の微分方程式の解は、
\(
\psi(x)=e^{\pm i{\sqrt{\frac{2mE}{\hbar^2}}x}}
\)

です。通常、波数\(k=\sqrt{\frac{2mE}{\hbar^2}}\)を定義して、
\(
\psi(x)=e^{\pm ikx}
\)

と書きます。連続状態を考えているので、エネルギー\(E\)は、
\(
\displaystyle E=\frac{\hbar^2k^2}{2m}
\)

です。

時間依存しないシュレーディンガー方程式が解けたので、時間依存するシュレーディンガー方程式の解\(\Psi(x,t)\)は、角周波数\(\omega=\frac{E}{\hbar}\)を用いて、
\(
\begin{align}
\Psi(x,t)&=e^{\pm ikx}e^{-i\frac{E}{\hbar}t} \\
&=e^{i(\pm kx-\omega t)}
\end{align}
\)

と記述されます。

\(e^{ikx}\)の解について


時間依存しない解が
\(e^{ikx}\)
に対する時間発展を考えます。
この時、時間依存する解は\(e^{i(kx-\omega t)}\)ですが、変形して、
\(
\displaystyle e^{ik(x-\frac{\omega}{k} t)}
\)

と書きます。
時刻\(t=0\)に対して、波動関数は\(e^{ikx}\)
時刻\(t=\Delta t\)に対して、波動関数は\(e^{ik(x-\frac{\omega}{k}\Delta t)}\)
です。よって、時刻\(t=0\to \Delta t\)の間に、関数が\(f(x)\to f(x-a)\)になった、つまり位置\(+a\)だけ平行移動された関数になった、とみることができます。
なので進行する波であると言えるのです。
また、その時の進行速度を考えますと、時間\(\Delta t\)の間に距離\(\frac{\omega}{k}\Delta t\)だけ進んだわけですから、速度\(v\)は、
\(
\displaystyle v=\frac{\frac{\omega}{k}\Delta t}{\Delta t}=\frac{\omega}{k}
\)

です。
グラフで時間経過を表せば、以下のようになります。
時刻t=0で\(e^{ikx}\)で表される波形(実線、実部:紫線、虚部:緑線)が時間経過して\(e^{ik(x-\frac{\omega}{k}\Delta t)}\)となると破線になります。
expikx_c
↑画像のgnuplotスクリプト

\(e^{-ikx}\)の解について


次に時間依存しない解が
\(e^{-ikx}\)
に対する時間発展を考えます。
この時、時間依存する解は\(e^{-i(kx+\omega t)}\)ですが、変形して、
\(
\displaystyle e^{-ik(x+\frac{\omega}{k} t)}
\)

と書きます。
時刻\(t=0\)に対して、波動関数は\(e^{-ikx}\)
時刻\(t=\Delta t\)に対して、波動関数は\(e^{-ik(x+\frac{\omega}{k}\Delta t)}\)
です。よって、時刻\(t=0\to \Delta t\)の間に、関数が\(f(x)\to f(x+a)\)になった、つまり位置\(-a\)だけ平行移動された関数になった、とみることができます。
なので後進する波であると言えるのです。
また、その時の進行速度を考えますと、時間\(\Delta t\)の間に距離\(-\frac{\omega}{k}\Delta t\)だけ進んだわけですから、速度\(v\)は、
\(
\displaystyle v=\frac{-\frac{\omega}{k}\Delta t}{\Delta t}=-\frac{\omega}{k}
\)

です。
グラフで時間経過を表せば、以下のようになります。
時刻t=0で\(e^{-ikx}\)で表される波形(実線、実部:紫線、虚部:緑線)が時間経過して\(e^{-ik(x+\frac{\omega}{k}\Delta t)}\)となると破線になります。
exp-ikx_c
↑画像のgnuplotスクリプト


波動関数は、
進行する波は虚部が先に来て後から実部(実部の位相が\(\pi/2\)遅れる)
後進する波は実部が先に来て後から虚部(虚部の位相が\(\pi/2\)遅れる)
のように見えます。
言い換えれば、波は虚部が先行している方向に進む、ということです。

3次元の場合


3次元でも同じことが言えます。
3次元自由粒子のハミルトニアンは
\(
\displaystyle \hat{H}=-\frac{\hbar^2}{2m}\left(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2}\right)
\)

で変数分離可能なので、3次元の1つの方程式ではなく、1次元の3つの方程式に変換可能であり、x,y,z方向は独立に決められます。
なので、
\(
\displaystyle \Psi(x,y,z,t)=e^{ik_x(x-\frac{\omega_x}{k_x} t)+ik_y(y-\frac{\omega_y}{k_y} t)+ik_z(z-\frac{\omega_z}{k_z} t)}
\)

と表されます。
ここで\(\omega_x, \omega_y, \omega_z\)は、1次元の3つの方程式に分けた時のそれぞれのエネルギー固有値\(E_x,E_y,E_z\)に対応する角周波数を示しています(\(\omega_i=E_i/\hbar,\;\; (i=x,y,z)\))。(エネルギーは\(E=E_x+E_y+E_z\)です。)

上の場合は\({\bf k}=(k_x,k_y,k_z)\)方向に進行する波を表します。
通常は波数ベクトルとして書いて、
\(
\begin{align}
& \Psi({\bf r},t)= e^{i{\bf k \cdot r}-\omega t} \\
&\;\;\;\;\; {\small (\omega=\omega_x+\omega_y+\omega_z)}
\end{align}
\)

と書かれます。

1次元調和振動子の直接的解法

1次元調和振動子の良く知られた解法は3種類あります。

  1. 昇降演算子を利用し\(\hat{a}\psi(x)=0\)から求める方法
  2. 微分方程式を境界条件の下、直接解く方法
  3. ハイゼンベルグ方程式から解く方法

です。ここでは直接的に解く2番目の方法を紹介します。

対象にする微分方程式は一次元調和振動子の原子単位系での表現
\(
\displaystyle \left[-\frac{1}{2}\frac{d^2}{dx^2}+Ax^2\right]\psi(x)=E\psi(x)
\)

です。ここでAは正の実数です。
この微分方程式を境界条件\(\left. \psi(x)\right|_{x\to \pm\infty}=0\)の下で解きます。

まとめますと、エネルギー\(E_n\)と固有関数\(\psi_n(x)\)は、
\(
\begin{align}
E_n&=\sqrt{2A}\left(n+\frac{1}{2}\right)\\
\psi_n(x)&=\frac{\left(2A\right)^{1/8}}{\pi^{1/4}\sqrt{2^n n!}}H_n( {\scriptsize (2A)^{1/4}}x)\exp\left(-\sqrt{\frac{A}{2}} x^2\right)\;\;\;\; (n=0,1,2,\cdots)
\end{align}
\)

です。


導出


両辺を\(-2\)を掛けて式変形して
\(
\displaystyle \left[\frac{d^2}{dx^2}-2Ax^2+2E\right]\psi(x)=0
\)

常套手段ですが、この微分方程式を満たすべき関数\(\psi(x)\)の漸近形を考えます。
何とか解ける部分だけを解いていこうという方針です。

漸近での解


もしも\(x\to \infty\)だったならば、上式の\(2E\)の項は定数であるため、\(-2Ax^2\)と比べて無視できるようになるはずです。
ただし、微分の項\(\frac{d^2}{dx^2}\)は比較することはできないため無視はしないで置いておきます。
つまり、関数\(\psi(x)\)は\(x\to \infty\)の漸近で
\(
\displaystyle \left[\frac{d^2}{dx^2}-2Ax^2\right]\psi(x)=0
\)

の形の微分方程式を満足しなければなりません。
式変形して
\(
\displaystyle \frac{d^2}{dx^2}\psi(x)=2Ax^2\psi(x)
\)

関数を2回微分したら\(x^2\times\)(自分自身)が出てくる、という性質を持つ関数は
\(
\displaystyle \psi(x)= C \exp(Bx^2)
\)

という形で掛けそうです。実際にこの形の2階微分は
\(
\begin{align}
\frac{d\psi}{dx}&=C\cdot 2Bx \exp(Bx^2) \\
\frac{d^2\psi}{dx^2}&=C\cdot 2B(1+2Bx^2)\exp(Bx^2) \\
&\sim C\cdot 4B^2x^2\exp(Bx^2)
\end{align}
\)
です。今、2階微分で括弧内に定数”1”が出てきていますが、\(x\to \infty\)の漸近を考えているのでこの定数は\(2Bx^2\)と比較して無視されるはずです。よって最後の式変形をしています。

漸近での微分方程式に代入すると、
\(
\begin{align}
\frac{d^2}{dx^2}\psi(x)&=2Ax^2\psi(x) \\
C 4B^2x^2\exp(Bx^2)&=2Ax^2 C \exp(Bx^2) \\
2C(2B^2-A)x^2\exp(Bx^2)&=0
\end{align}
\)
上記の方程式は\(x\)の値に依らず常に満たされなければなりません。
よって、\(2C=0\)もしくは\(2B^2-A=0\)である必要があります。
\(C=0\)の解は波動関数が至る所でゼロである自明な解です。なので物理的には適しません。

※漸近でゼロであるだけで漸近じゃないところではokなんじゃないの?という疑問が出てきますが、\(C=0\)の場合、漸近で\(\psi(x)=0\)であり、その微分\(\frac{d\psi}{dx}=C 2Bx \exp(Bx^2)\)もゼロです。波動関数は全領域に渡って滑らかに続いていなければならないため、いつまでたっても波動関数はゼロです。ということは、存在確率がゼロ、有限にはならないつまらない解なのです。通常、\(x\to \infty\)の漸近で\(\psi(x)=0\)という条件は入れますが、その微分\(\frac{d\psi}{dx}\)もゼロという条件は入れません。

なので物理的に意味のある解は
\(
\begin{align}
2B^2-A=0 \\
\to B=\pm \sqrt{\frac{A}{2}}
\end{align}
\)
というBでなければなりません。よって波動関数\(\psi(x)\)は漸近で2つの解の線形結合として書かれ、
\(
\displaystyle \psi(x)= C_1 \exp\left(-\sqrt{\frac{A}{2}} x^2\right)+C_2 \exp\left(\sqrt{\frac{A}{2}}x^2\right)
\)

という形で記述できます。
しかし、境界条件\(\left. \psi(x)\right|_{x\to \pm\infty}=0\)を満たさなければならない場合、係数\(C_2=0\)である必要があります。よって波動関数の漸近形は
\(
\displaystyle \psi(x)= C\exp\left(-\sqrt{\frac{A}{2}} x^2\right)
\)

と求められます。

漸近ではないところの解


漸近では波動関数は
\(
\displaystyle \psi(x)= C\exp\left(-\sqrt{\frac{A}{2}} x^2\right)
\)

という形を満たさなければならないことが分かりました。残りは漸近ではない場所の解です。
定数変化法を使います。定数をxに依存する関数として全領域に渡る波動関数を
\(
\displaystyle \psi(x)= f(x)\exp\left(-\sqrt{\frac{A}{2}} x^2\right)
\)

と仮定します。この操作は\(\psi(x)\)に関する方程式から\(f(x)\)に関する微分方程式に焼直す操作です。
全領域に渡る波動関数の満たすべき微分方程式は
\(
\displaystyle \left[\frac{d^2}{dx^2}-2Ax^2+2E\right]\psi(x)=0
\)

なので、これに代入します。
計算過程で出てくる文字を減らすために、\(B=-\sqrt{\frac{A}{2}}\)として、\(B\)で記述していきます。
\(
\begin{align}
&\left[\frac{d^2}{dx^2}-2Ax^2+2E\right]\psi(x)= \\
&\;\;\;\;\;\;\;\left\{\frac{d^2f}{dx^2}+4Bx\frac{df}{dx}+(4B^2-2A)x^2+(2B+2E)f\right\}e^{Bx^2}=0
\end{align}
\)
と求める事ができ、\((4B^2-2A)=0\)を利用すると、\(f(x)\)の満たすべき微分方程式は、
\(
\displaystyle \frac{d^2f}{dx^2}+4Bx\frac{df}{dx}+(2B+2E)f=0
\)

と求まります。
この微分方程式の解はエルミートの微分方程式と呼ばれる形をしています。
エルミートの微分方程式は
\(
\displaystyle \frac{d^2 H(x)}{dx^2}-2x\frac{dH(x)}{dx}+2nH(x)=0
\)

の微分方程式を満たす直交多項式です[1]。
そのままエルミートの微分方程式を今回の問題に適応することはできません。なぜならば1階微分の係数が違っているからです。
なので変数変換により、\(x=\alpha y\)とおいて、一階微分の係数\(4B\)をちょうど\(-2\)にする定数\(\alpha\)を探しましょう。
\(
\begin{align}
\frac{d}{dx}&=\frac{dy}{dx}\frac{d}{dy}=\frac{1}{\alpha}\frac{d}{dy}\\
\frac{d^2}{dx^2}&=\frac{1}{\alpha^2}\frac{d^2}{dy^2}
\end{align}
\)
を代入して、
\(
\displaystyle \frac{1}{\alpha^2}\frac{d^2 f(x)}{dy^2}+4By\frac{df(x)}{dy}+(2B+2E)f(x)=0
\)

xとyが混在しているので気持ち悪いです。\(f(x)\to f(y)\)とおきましょう。そして最後に求められたyを用いたf(y)をxの式f(x)に焼直します。
両辺に\(\alpha^2\)を掛けると
\(
\displaystyle \frac{d^2 f(y)}{dy^2}+4B\alpha^2 y\frac{df(y)}{dy}+(2B+2E)\alpha^2 f(y)=0
\)

です。今、一階微分の係数を\(-2\)にする定数\(\alpha\)を探していました。
なので、
\(
\begin{align}
4B\alpha^2&=-2 \\
\alpha&=\pm\sqrt{-\frac{1}{2B}}=\pm\left(\frac{1}{2A}\right)^{1/4}
\end{align}
\)
となります。ここの\(\alpha\)のプラスマイナスはどちらでもいいです。
プラスに取ることにすれば、\(f(y)\)についての微分方程式は、
\(
\displaystyle \frac{d^2 f(y)}{dy^2}-2 y\frac{df(y)}{dy}-\left(\frac{B+E}{B}\right) f(y)=0
\)

であるため、エルミートの微分方程式の形になりました。
エルミートの微分方程式は\(f(y)\)に掛かる項、すなわち\(-\left(\frac{B+E}{B}\right)\)の部分が\(2n\), (\(n=0,1,2,\cdots\))である時だけ解を満たします。
よって、
\(
\begin{align}
-\frac{B+E}{B}&=2n\\
E=E_n&=-2B(n+\frac{1}{2})\\
&=\sqrt{2A}(n+\frac{1}{2})\;\;\; (n=0,1,2,\cdots)
\end{align}
\)
を満たすエネルギー\(E=E_n(n=0,1,2,\cdots)\)でなければ解は存在しません。
その時の解\(f(y)\)も\(n\)で番号付けされて、
\(
f(y)=H_n(y)\;\;\; (n=0,1,2,\cdots)
\)

と記述され、右辺はエルミート多項式と呼ばれます。
エルミート多項式\(H_n(y)\)は
\(
\begin{align}
H_0(y)&= 1\\
H_1(y)&= 2y\\
H_2(y)&= 4y^2-2\\
H_3(y)&= 8y^3-12y\\
H_4(y)&= 16y^4-48y^2+12\\
H_5(y)&= 32y^5-160y^3+120y\\
H_6(y)&= 64y^6-480y^4+720y^2-120\\
\ldots
\end{align}
\)
と言う関数です。三項漸化式
\(
H_n(y)=2yH_{n-1}(y)-2(n-1)H_{n-2}(y)
\)

を用いて\(n=2\)以上のエルミート多項式は記述されます。

以上より微分方程式を満たす解とエネルギーは、整数(\(n=0,1,2,\cdots\))で番号付けされた解を持つことが分かりました。
\(y\)を変数変換に従い\(x\)に戻します。
\(
\begin{align}
x&=\alpha y \\
&\alpha=\left(\frac{1}{2A}\right)^{1/4}
\end{align}
\)
なので、
\(H_n(y)=H_n(x/\alpha)\)
です。
漸近形での波動関数と合わせて、本当に求めたかった解\(\psi(x)\)は、
\(
\begin{align}
\displaystyle \psi(x) &= f(x)\exp\left(-\sqrt{\frac{A}{2}} x^2\right) \\
&=C H_n(x/\alpha)\exp\left(-\sqrt{\frac{A}{2}} x^2\right)
\end{align}
\)

であることが分かりました。定数\(C\)の自由度がありますが、これは物理的に合った条件を用いてなされます。

おさらい


解きたかった問題は
\(
\displaystyle \left[-\frac{1}{2}\frac{d^2}{dx^2}+Ax^2\right]\psi(x)=E\psi(x)
\)

を境界条件\(\left. \psi(x)\right|_{x\to \pm\infty}=0\)の下で解くことでした(A>0)。

この微分方程式を何とか解くために、まず\(x\to \pm\infty\)で解の満たすべき形を考えました。
その後、全領域で境界条件を満たす解を定数変化法によって求めました。
計算の結果、微分方程式の解は与えられた適当なエネルギー\(E\)に対していつも存在するわけではなく、特別な\(E=E_n\)で番号付けされた値でなければ解が存在しないことが分かりました。
\(n\)番目のエネルギー\(E_n\)と関数\(\psi_n(x)\)は、
\(
\begin{align}
E_n&=\sqrt{2A}\left(n+\frac{1}{2}\right) \\
&\psi_n(x)= C_n H_n(x/\alpha)\exp\left(-\sqrt{\frac{A}{2}} x^2\right) \\
&\;\;\;\;\alpha=\left(2A\right)^{-1/4}
\end{align}
\)

です。与えられた条件下では規格化までは追求していないのでここでお仕舞いですが、規格化を考えてみましょう。

規格化


波動関数\(\psi_n(x)\)を規格直交化(正規直交化とも言う)します。
規格化は、全空間での存在確率
\(
\displaystyle \int_{-\infty}^{\infty}|\psi(x)|^2dx
\)

を1にすること。
直交化は異なる状態間をかけ合わせて積分した時ゼロになること・・・すなわち内積をゼロにします。
\(
\displaystyle \int_{-\infty}^{\infty}\psi_m^*(x)\psi_n(x)dx=C\cdot\delta_{mn}
\)

ここで\(C\)は任意の定数で、\(\delta_{mn}\)はクロネッカーのデルタを表します。
両方を考慮すると、規格直交化とは、
\(
\displaystyle \int_{-\infty}^{\infty}\psi_m^*(x)\psi_n(x)dx=\delta_{mn}
\)

を満たすように係数を決定します。

実際に波動関数を代入すると、
\(
\begin{align}
\int_{-\infty}^{\infty}\psi_m^*(x)\psi_n(x)dx&= C_m^*C_n \int_{-\infty}^{\infty} H_m(x/\alpha)H_n(x/\alpha)\exp\left(-\sqrt{2A} x^2\right)dx
\end{align}
\)

変数変換を行います。
\(
\begin{align}
x&=\alpha y \\
&\alpha=\left(\frac{1}{2A}\right)^{1/4}
\end{align}
\)
より、
\(
\begin{align}
& C_m^*C_n \int_{-\infty}^{\infty} H_m(x/\alpha)H_n(x/\alpha)\exp\left(-\sqrt{2A} x^2\right)dx \\
&= C_m^*C_n \int_{-\infty}^{\infty} H_m(y)H_n(y)\exp\left(-\sqrt{2A}\alpha^2y^2\right)\cdot \alpha dy \\
&= C_m^*C_n\alpha \int_{-\infty}^{\infty} H_m(y)H_n(y)\exp\left(-y^2\right) dy \\
\end{align}
\)
ここで、エルミート多項式の特性を使います。エルミート多項式の関係式の中に、
\(
\displaystyle \int_{-\infty}^{\infty} H_m(x)H_n(x)\exp\left(-y^2\right) dx=2^n n!\sqrt{\pi}\delta_{mn}
\)

という関係式があります。
なので、これを使うと、
\(
\displaystyle \int_{-\infty}^{\infty}\psi_m^*(x)\psi_n(x)dx=C_m^*C_n\alpha 2^n n!\sqrt{\pi}\delta_{mn}
\)

と求められます。今、\(({\mbox 右辺})=\delta_{mn}\)にしたいので、係数について、
\(
\begin{align}
C_m^*C_n\alpha 2^n n!\sqrt{\pi}\delta_{mn}&=\delta_{mn} \\
(C_m^*C_n\alpha 2^n n!\sqrt{\pi}-1)\delta_{mn}&=0
\end{align}
\)
にするように\(C_n\)を決定すればよいことになります。
クロネッカーのデルタより、\(m=n\)以外の時は係数に依らずゼロになるため、自明です。
なので、\(m=n\)の時,
\(
\begin{align}
|C_n|^2\alpha 2^n n!\sqrt{\pi}&=1\\
\to |C_n|&=\sqrt{\frac{1}{\sqrt{\pi}2^n n!\alpha}} \\
C_n&=\frac{1}{\sqrt{\alpha}\pi^{1/4}\sqrt{2^n n!}}\cdot \exp(i\theta)
\end{align}
\)
と規格化定数が求められます。位相\(\exp(i\theta)\)の任意性がありますが、この位相に特に意味は無く、自由に選べます。よって、式が簡単になる\(\theta=0\)を取ることにします。
以上の結果から、規格直交化された波動関数は\(\psi_n(x)\)は、
\(
\begin{align}
\psi_n(x)&=\frac{1}{\sqrt{\alpha}\pi^{1/4}\sqrt{2^n n!}}H_n(x/\alpha)\exp\left(-\sqrt{\frac{A}{2}} x^2\right)\\
\alpha&=\left(2A\right)^{-1/4}
\end{align}
\)
となります。\(\alpha\)を使わず、\(A\)だけで表せば、
\(
\displaystyle \psi_n(x)=\frac{\left(2A\right)^{1/8}}{\pi^{1/4}\sqrt{2^n n!}}H_n( {\scriptsize (2A)^{1/4}}x)\exp\left(-\sqrt{\frac{A}{2}} x^2\right)
\)

です。

参考文献

[1]小野寺嘉孝著「物理のための応用数学」裳華房 p.71

バレル内部でのBB弾の方程式

バレル内部でのBB弾の運動方程式です。

目的は、

BB弾は何秒間バレル内部に存在しているのか?
バレルが長いと減速になりうるのか?

を知ることです。詳細は弾道計算本をご覧ください。

弾道計算(BB弾)の理論
弾道計算(BB弾)の結果
弾道計算の結果2, 比較と詳細データ
弾道計算(BB弾)のコード(fortran90)
バレル内部でのBB弾の方程式←今ここ


電動ガンの場合です。
バレル内部1
空気抵抗は初速0[m/s]から90[m/s]前後にまで加速されるわけですから、空気抵抗の粘性抵抗、慣性抵抗どちらの項も無視することはできないでしょう。

空気の漏れも考えます。

この条件下では、バレル内部の運動方程式は以下のように導くことができます。

内部圧力変化と力


ピストン-BB弾間の圧力変化による力\(F_V(t)\)は、
\(
F_V(t)=S\cdot P(t)
\)

と書けます。ここで\(P(t)\)はピストン-BB弾間の内部圧力です。
\(P(t)\)は、断熱変化を仮定すると、ピストン-BB弾間体積\(V(t)\)を用いて
\(
P(t)V^{7/5}(t)=\text{const}
\)

の関係があります。体積\(V(t)\)は
\(
V(t)=x(t)S_b-x_0(t)S_s
\)

と書けます。ここで、ピストンの位置による時間変化を\(x_0(t)\)、BB弾の位置を\(x(t)\)、バレルの断面積を\(S_b\)、シリンダーの断面積を\(S_s\)としました。
この式は時刻に依存せずに決まるので、時刻\(t=0\)の時に大気圧\(P_0\)で、体積\(V_0\)であれば、任意の時刻での圧力\(P(t)\)は
\(
\displaystyle P(t)=P_0 \left(\frac{V_0}{V(t)}\right)^{7/5}
\)

と書けます。

空気の漏れについて


十分短い時間の間、時刻\(t\)において内部圧力\(P(t)\), 空気の密度\(\rho\)とすると、漏れ出る空気の速度\(v(t)\)はベルヌーイの定理から
\(
\begin{align}
P(t)&=\frac{1}{2}\rho v(t)^2+P_0 \nonumber \\
&\to v(t)=\pm \sqrt{\frac{2}{\rho}(P(t)-P_0)}
\end{align}
\)
が成立するとします。厳密には、ピストンの速度は漏れ出る空気の速度に対して無視できるほど小さい、という仮定の下で成立します。漏れ出る流量\(Q\)は、バレル-BB弾の隙間の断面積\(A\)、実験とのズレを調節する無次元の流量係数\(c’\)を用いて\(Q=c’ A v(t)\)と書けます。
また、断熱過程を圧力-体積の関係を用いたいので、漏れ出る空気はBB弾-ピストン間の体積の増加として扱います。

フルシリンダーの場合の運動方程式


その他、空気抵抗による力を入れると、BB弾の位置\(x(t)\), ピストンの位置\(x_0(t)\), 体積\(V(t)\)の運動方程式は
\(
\begin{align}
m\frac{d^2}{dt^2}x(t)&=\left[P(t)-P_0\right]S_b-\left(6\pi \eta R |v(t)|+\frac{1}{2}C_d \rho \pi R^2 |v(t)|^2\right)\frac{v(t)}{|v(t)|} \label{bbin1}\\
m_s\frac{d^2}{dt^2}x_0(t)&=-k\left[x_0(t)-x_B-l\right]-[P(t)-P_0]S_s -F_f\label{bbin2} \\
\frac{d}{dt}V(t)&=v(t)S_b-v_0(t)S_s+c'(S_b-\pi R^2)\mathrm{sgn}(P(t)-P_0)\sqrt{\frac{2}{\rho}|P(t)-P_0|} \label{bbin3}
\end{align}
\)
と導くことが出来ます。

ここで、\(P_0\)は大気圧、\(\eta\)は粘性率、\(R\)はBB弾の半径、\(C_d\)は抗力係数、\(\rho\)は空気の密度, \(S_b\)はバレルの断面積、\(S_b\)はシリンダーの断面積、\(k\)はピストンのばね定数、\(l\)はばねの自然長、\(F_f\)はピストン-シリンダーの摩擦、\(v(t),v_0(t)\)はそれぞれBB弾の速度、ピストンの速度を意味します。

空気抵抗に関する詳細は球体の空気抵抗と係数をご覧ください。

ここで1つ言えることは、外部と内部の圧力が一定になる最適なバレル長というものが存在する、ということです。

さて、計算結果はこのページには載せません。弾道計算本として計算結果をまとめた本を出しているので、そちらをご参照ください。本の詳細については弾道計算本の自家通販をご覧ください。

補遺


補遺1
理想気体として扱える条件は、
低い気圧(分子の数が少なく、衝突等が無視できる)かつ高い温度(分子間力が分子の運動エネルギーに比べて無視できる)
です。
どうやら実在気体では10気圧以下、室温以上でこの条件は良く満たされ、理想気体とのずれは1%以内のようです[1]。
[1]を引用すると、

一般に,沸点の低い酸素・窒素・水素・ヘリウム等は,室温またはそれ以上の温度で10atm以下の圧力の場合,理想気体の値の1%以内で理想気体に近い性質を示す。

とあります。ピストンで空気が圧縮されたとき、BB弾とピストン間の圧力が10気圧以上にならなければ良い近似だと言えるでしょう。

参考文献


[1] 実在気体 -第2節 気体の状態方程式
[2]PERFECT HIT -TOKYO MARUI

『集弾性アップへの道』 BB弾とバレル内部に隙間があることが写真で確認できます。