スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。
スポンサー広告 | --/--/-- --:--

一次元クーロン場の固有状態の数値計算 (4)

これまでの計算では、基底準位が -0.67 a.u. となっていましたが、
水素原子を模擬するには、 -0.5 a.u. ( = 13.6 eV) となってほしいところです。
たとえば、光を当てた時などに、光子エネルギー hν との関係が重要だから。

soft-core potential\[
V(x) = -\frac{1}{\sqrt{x^2 + a^2}}
\tag{1}
\]の a を調整すればよさそうと思っていたのですが、
論文 [3] で、a = √2 (a.u.) にすれば、E = -0.5 (a.u.) になる
と書かれていました。

さっそく、計算してみた。
Coulomb1DEValuesAsqrt2Xinf100Log.png
基底準位の計算値は・・・
E = -0.499999999995472

見事に、0.5 a.u. になってます!

こんなぴったりなるというのはきっと偶然じゃないだろうから、
実は解析的に解けるんでしょうか???

ポテンシャルを描画してみると、
1DCoulombSoftCorePotential01.png

青が a = 1 で、赤が a = √2 、緑はクーロン。
ポテンシャルの底が上がった分、基底準位が持ち上がったんですね。
それにしても、ちょうど√2 でうまくいくというのは不思議!


[1] J. H. Eberly et al. "Numerical Experiments in Strong and Super-Strong Fields"
in M. Gavrila Ed. "Atoms in Intense Laser Fields"
[2] J. Javanainen, J. H. Eberly, Q. Su, "Numerical simulations of multiphoton ionization and above-threshold electron spectra," Phys. Rev. A 38, 3430 (1988)
[3] A. Gordon, R. Santra, F. Z. Kartner, "Role of the Coulomb singularity in high-order harmonic generation," Phys. Rev. A 72, 063411 (2005)
スポンサーサイト
数値計算>量子力学 | コメント(0) | 2015/02/06 20:23

一次元クーロン場の固有状態の数値計算 (3)

計算結果の正確な数値が載っている論文 [2] を見つけたので、
これをプロットして、僕の計算結果と比較してみた。
Coulomb1DEValuesEVXinf100n2scaledLogPaper.png

n = 10 までは合ってるようですね!

上の方の励起状態が直線からはずれているのは、
無限遠の設定が悪いから。

今度は、無限遠を10倍の距離 x = 1000 に設定して、
-0.1 ≦ E ≦ 0 の領域で、エネルギー間隔を 1/10 の ΔE = 0.0001 にして、
固有値を求めてみた。

Coulomb1DEValuesEVXinf1000Log.png

今度は、 E = -0.001 付近の n = 40 あたりまで、直線に乗っている。
$1/\sqrt{|E|}$ のプロットも分かりやすいので。
Coulomb1DEValuesEVXinf1000n2scaled.png

n = 45 ぐらいまでの励起状態は正しく求められてそうですね。


参考文献
[1] J. H. Eberly et al. "Numerical Experiments in Strong and Super-Strong Fields"
in M. Gavrila Ed. "Atoms in Intense Laser Fields"
[2] J. Javanainen, J. H. Eberly, Q. Su, "Numerical simulations of multiphoton ionization and above-threshold electron spectra," Phys. Rev. A 38, 3430 (1988)
数値計算>量子力学 | コメント(0) | 2015/02/04 18:37

一次元クーロン場の固有状態の数値計算 (2)

Soft-core pontential で a = 1 (ボーア半径)の場合を計算してみる。\[
V(x) = -\frac{1}{\sqrt{x^2 +1}}
\tag{1}
\]まずは、ポテンシャルの形を実際のクーロンポテンシャル\[
V(x) = -\frac{1}{|x|}
\tag{2}
\]と比較すると、こんな感じ。
Coulomb1D01.png
結構違ってるけど、まあ、遠く離れると、一致してます・・・

調和振動子の時と同様に、無限遠における零点を直線補間で求めてみる。

まずは、-0.9 ≦ E ≦ 0 の範囲を間隔 ΔE = 0.001 刻みで。
無限遠点は x = 100 に設定。

下から順番に固有エネルギーをプロット。
Coulomb1DEValuesXinf100.png
分かりやすいように、1 a.u. = 27.2114 eV を使って、
eV 単位に変換したものはこちら。
Coulomb1DEValuesEVXinf100.png
文献[1] によると、基底状態は、E = -0.67 a.u. と書かれていて、
こちらの計算結果は E = -0.669761 a.u. となっているので、
基底状態については、計算結果あってそうです。

eV で言うと、 約 18 eV でネオンの基底準位に近いらしいですね。
水素は 13.6 eV ( = 0.5 a.u.) なので、a の値を調整すれば近づけられるんだろうか・・・

文献 [1] によれば、3次元の場合と同様に、
エネルギー固有値は、$1/n^2$ のスケールに乗るらしい。

[1] に載っているように、$1/\sqrt{|E_n|}$ をプロットしてみた。
Coulomb1DEValuesEVXinf100n2scaled.png
確かに乗っている。上の方は、E = -0.01 に近づいてきて、
無限遠の設定値 x = 100 が古典的許容領域に入ってきてしまっているからでしょう。

上の方の励起状態は、無限遠をもっと大きい値に設定して、
エネルギー間隔も狭くなってくるので、さらに細かく刻んで計算する必要いがありそう。

両対数プロットしても、直線上に乗るはず。
Coulomb1DEValuesEVXinf100n2scaledLog.png
基底状態 (n=1) のところで少し曲がってるようにも見えますね・・・
この付近は、クーロンポテンシャルからはずれているので、
直線に乗らなくても問題はないのですが。

次回は、上の方の励起状態を細かく計算してみたいと思います。

参考文献
[1] J. H. Eberly et al. "Numerical Experiments in Strong and Super-Strong Fields"
in M. Gavrila Ed. "Atoms in Intense Laser Fields"
数値計算>量子力学 | コメント(0) | 2015/02/02 12:07

一次元クーロン場の固有状態の数値計算 (1)

解析的に解けない例として(たぶん解けないと思う)、計算してみたかったのが
一次元クーロン場の固有状態。

一次元クーロンポテンシャルはこんな感じのもの。\[
V(x) = -\frac{e^2}{4\pi\varepsilon_0 x}
\tag{1}
\]
こんなものは実際には存在しませんが、
3次元クーロンポテンシャルを疑似的に扱うのに使えるようです。

ただし、原点に特異性があると扱いにくいので、\[
V(x) = -\frac{e^2}{4\pi\varepsilon_0 \sqrt{x^2 + a^2}}
\tag{2}
\]と修正する( soft-core potential と呼ぶそうです)

原点の特異性が取り除かれて、
x >> a の遠方では (1) のクーロンポテンシャルに漸近する。
こうすることによって、原点を通過する電子の動きが計算できるとのこと。
確かに、便利そうですね!

さっそく、シュレディンガー方程式を立てる。\[
-\frac{\hbar^2}{2m} u'' - \frac{e^2}{4\pi\varepsilon_0 \sqrt{x^2+a^2}} u = Eu
\tag{3}
\]原子単位系を用いて、長さの単位をボーア半径\[
x_0 = \frac{4\pi\varepsilon_0 \hbar^2}{me^2}
\tag{4}
\]エネルギーの単位をハートリーエネルギー\[
E_0 = \frac{e^2}{4\pi\varepsilon_0 x_0}
= \frac{me^4}{(4\pi\varepsilon_0)^2 \hbar^2}
\tag{5}
\]としてスケーリングすると、\[
-\frac{1}{2} u'' - \frac{1}{\sqrt{x^2+a^2}} u = Eu
\tag{6}
\]と簡単化できる。u'' について解いて、\[
u'' = -2\left( E + \frac{1}{\sqrt{x^2+a^2}} \right) u
\tag{7}
\]を4次ルンゲクッタで計算する。


参考文献
[1] J. H. Eberly et al. "Numerical Experiments in Strong and Super-Strong Fields"
in M. Gavrila Ed. "Atoms in Intense Laser Fields"

数値計算>量子力学 | コメント(0) | 2015/01/27 13:04

調和振動子の固有状態の数値計算 (7)

ちょっと本題からずれるのですが、
ついでに、もう少し高い励起状態までグラフにしてみました。
確率密度がどういう形になってるんだろう?と興味があったので。

数値計算ではなく、以下の解析解で計算しています。
規格化も入れています。\[
u_n(x) = (\pi^{1/2} 2^n n!)^{-1/2} H_n(x) e^{-x^2/2}
\tag{1}
\]

こちらが波動関数 $u(x)$ そのもの。
1DTISEHarmonicOscAnalyticEigenFunctions01.png

こちらが波動関数の2乗の確率密度 $|u(x)|^2$ をプロットしたもの。
1DTISEHarmonicOscAnalyticEigenFunctions02.png

なんとなく、イメージわきました(笑)

参考文献
[1] シッフ 「量子力学」(上)
数値計算>量子力学 | コメント(0) | 2015/01/27 12:38

調和振動子の固有状態の数値計算 (6)

固有値を求めるだけでなく、固有関数の方もちゃんと計算できてるか?
ということを確認。

初めの6個の固有状態を計算した結果。
青が偶関数で、赤が奇関数。規格化はしていません。
折り返した方 (x < 0) も書くと分かりやすいんでしょうけど、面倒なので・・・

1DTISEHarmonicOscEigenFunction01.png

固有エネルギーの分だけ下駄をはかせるお決まりの方法で表示してみました。

白い部分が古典的に許容される領域で、
影の部分は古典的には存在しえないが、
量子論的には波動関数が少し浸み出しているというのがこの表示だとよくわかります。

解析解\[
u(x) = a_n H_n(x) e^{-x^2/2}
\]の方も同様に書いてみて、比較してみました。
1DTISEHarmonicOscEigenFunction02.png

数値計算の方は規格化していないので、こちらの方の係数 $a_n$ を
数値計算の結果に合うように、適当に調整しています。

まあ、見た目ほぼ一致してるようですね。
どの程度の精度で一致しているかまでは評価してないけど、
まあこんなところでいいってことで・・・

ところで、実は調和振動子の計算をしたかったわけではありません。

解析解が分かるもので、正しく計算できることをテストしておいて、
次に、解析的には解けないようなポテンシャルで計算してみたいのです。

というわけで、次回からいよいよ解けないポテンシャルで計算してみようかと思います。
数値計算>量子力学 | コメント(0) | 2015/01/22 20:15

調和振動子の固有状態の数値計算 (5)

前に発散してしまうからという理由で、x = 2√E を無限遠に設定していたが、
前回の結果で、意外にも多少発散していても、零点は見つけられるということが分かったので、
無限遠を x = 30 に固定して、固有値の計算をやり直してみた。

これがその結果。

1.0000000001061
5.00000000687847
9.00000003899501
13.0000001166616
17.0000002600491
21.0000004893016
25.0000008245405
29.0000012858667
33.0000018933623
37.0000026670919
41.000003627103
45.0000047934266
49.0000061860783
53.0000078250579
57.0000097303504
61.0000119219255
65.0000144197385
69.0000172437302
73.0000204138268
77.0000239499403
81.0000278719684
85.0000321997946
89.0000369532885
93.0000421523053
97.0000478166863


固定したので、最初の数個もちゃんと精度が出ている。
ちなみに、x = 40 にしてみたら、さすがにやりすぎで、
低エネルギーの方では固有値が求められなかった。
(エネルギー掃引の刻みをもっと細かくすれば行けるかもしれないが・・・)

前回のように、固有値をグラフにしてみると、
1DTISEHarmonicOscEvaluesLinInterpInf30-01.png

相対誤差は、
1DTISEHarmonicOscEvaluesLinInterpInf30-02.png

今回は、すべての点で6桁以上の精度が得られている。

数値計算>量子力学 | コメント(0) | 2015/01/22 12:39

調和振動子の固有状態の数値計算 (4)

自動的に直線補間して、固有値を一気に求めるプログラムができた。
エネルギーを掃引しながら、無限遠での零点を見つけては、
直線補間して、零点の正確な値を出力していく。

零点の求め方は、まず、符号が反転したところを見つける。
順次エネルギーを掃引していき、$E_k$ と$E_{k+1}$ に対して、
$u_k u_{k+1} < 0$ となったら、符号が反転していると判定できる。

その時に、直線補間で零点の正確な値を求める。
その値を E とすると、\[
E = \frac{u_{k+1} E_k - u_k E_{k+1}}{u_{k+1} - u_k}
\]
1DTISEHarmonicOscEValuesLinInterpConcept01.png

以上のアルゴリズムで、E = 0.5~100 の範囲で求めた固有値一覧がこちら。

1.06019474008538
5.00001293981931
9.00000003919728
13.0000001157887
17.00000025828
21.0000004862375
25.0000008197467
29.0000012788781
33.0000017062101
37.0000026542532
41.0000036105981
45.0000047727627
49.0000054600211
53.0000077946515
57.0000085103253
61.000011879947
65.0000143713303
69.0000171885402
73.0000203514719
77.0000205266203
81.0000277945322
85.0000273957349
89.0000368597403
93.0000420505301
97.0000477065965


エネルギーの刻みは前と同様、0.01に設定し、
無限遠も前と同様、 x = 2√E に設定した。

僕の使っているコンパイラでは、double 型は仮数部が 52 bit らしいので、
$\log_{10} 2^{52} \simeq 15.7$ より
一応、15 桁ぐらいが有効だろうということで、
出力ストリームで precision(15) と指定して出力している。
(なぜか、15桁になっていないところがあるのは謎)

E > 81 で発散していると思っていたが、ちゃんと E = 100 まで計算できている。
前は、振動の振幅が大きすぎて、グラフソフトがうまく表示できていなかっただけらしい。

求めた固有値をグラフにすると、
1DTISEHarmonicOscEvaluesLinInterp01.png
表を見ても明らかであるが、E = 4n + 1 の理論直線上に並んでいるのが分かる。
(偶関数のみで計算しているので、4n + 3 の系列は現れない)

次に、理論値からの相対誤差 $|\Delta E|/E$ を調べてみた。
1DTISEHarmonicOscEvaluesLinInterp02.png
始めの2点の誤差が大きいのは、無限遠の値 2√E が小さすぎるためだろう。
それ以外は、6桁以上の精度が出ているので、とりあえず、十分。
数値計算>量子力学 | コメント(0) | 2015/01/21 19:20

調和振動子の固有状態の数値計算 (3)

調和振動子の固有状態の数値計算の続き。

予告した通り、無限遠を古典軌道の振幅($\sqrt{E}$)の2倍の点、
すなわち、 $x = 2\sqrt{E}$ のところに設定してみた。

x の刻みは h = 0.01。
偶関数の初期条件 $u(0) = 1$、$u'(0) = 0$ で、
各 E に対して、$x = 2 \sqrt{E}$ の時の u の値 $u_{\rm inf}$ を
4次ルンゲ・クッタで計算して、プロットしたものがこちら。

1DTISEHarmonicOscUinf01.png

E の刻みは、前回よりも細かくして、0.01 刻みにした。
E = 0.5 から 100 までおよそ 10000 個の値に対して計算しているが、
ノートパソコンで数十秒ほどで計算できてしまう(速いね!^^)

見てるウィンドウ |u|≦1 の中にデータ点がほとんど入ってないので、
ダメじゃんって思ったけど、もう少し細かく見てみることにする。
(ちなみに、E > 81 のところでは完全に発散してしまった)

まずは、E = 13 付近を拡大。

1DTISEHarmonicOscUinf02.png

計算したのは、0.01 刻みだが、グラフソフトが線形補間してくれる。
(E = 13.0000000 のところにデータ点がないのは、点を表示する設定を忘れてしまっただけ)

このグラフでゼロクロスしているところが固有値になる。
つまり、E = 13.00000012 が計算された固有値。

意外と、いい精度出てる!!!

真の理論値が E = 13 だから、
ΔE / E = 0.00000012 / 13 ≒ 10-8
8桁の精度が出てることになる。
こんな適当な計算で8桁も出るとはびっくり(笑)

次に、E = 61 の付近を見てみた。
1DTISEHarmonicOscUinf03.png

こちらは先ほどよりは精度が悪くなっているが、それでも7桁近く出ている。
これぐらいなら使えそう。

最後に、前回の無限遠を x = 5 に設定したケースで同じ部分を拡大して比較してみた。
1DTISEHarmonicOscUinf04.png
公平を期すために、こちらも E の刻みを 0.01 刻みで計算し直した。
2 √13 ≒ 7.2 なので、5 とそれほど変わらないけど、それでもだいぶ精度が違うようである。
(もっと E の大きいところで比較してみればよかったなあ・・・)

というわけで、線形補間すれば、結構な精度が出ることが分かった。

次は、線形補間を自動化して、固有値の近似値をすべて一気に計算できるように
プログラミングしてみようと思う。
数値計算>量子力学 | コメント(0) | 2015/01/20 22:35

調和振動子の固有状態の数値計算 (2)

前回、あらかじめ分かっている固有値の付近で E の値を動かして、
無限遠で発散しない解を見つけたが、
あらかじめ答が分かっていない場合を想定して、
E を広範囲で掃引して、解を見つけたい。

ここで悩ましいのが無限遠の設定

前回の計算では、とりあえず、x = 5 に設定していた。
x の値は、ポテンシャルが $V = x^2$ になるように無次元化されている。

1DTISEHarmonicOsc03.png

古典的に許容される領域は、$E \geq V = x^2$ だから、
古典軌道の端点は、$x = \pm \sqrt{E}$。

前回のように、 E = 1 や 3 であれば、x = 5 を無限遠としてもよさそう。
E を大きくしていって、E = 25 になると、x = 5 が古典軌道の端点になってしまい、
さらに E を大きくすると、古典的に許容される領域に入ってしまう。

無限遠を初めからもっと大きい値に設定すればよいのだが、
そうすると、途中の誤差が蓄積してしまうせいなのか、
E が固有値になっている場合でも、発散してしまう。

とりあえず、E = 0 から E = 50 まで E を 0.5 刻みで掃引して、
x = 5 における値を計算してみた。

偶関数の場合。
1DTISEHarmonicOsc04.png
奇関数の場合。
1DTISEHarmonicOsc05.png

x = 5 における値がゼロになる点(このグラフの零点)が固有値。

たとえば、偶関数の方で見ると、初めのうちは、E = 1, 5, 9 と 4 ずつ増えているが、
そのうち、少しずつずれてくる(赤丸の位置が零点とずれている)
これは、x = 5 が無限遠とみなせなくなるから。

たとえば、無限遠を古典領域の端点の2倍のところ、
すなわち、$ x = 2\sqrt{E}$ などと取ると、どうなるだろう?
やってみたいと思います。
数値計算>量子力学 | コメント(0) | 2014/12/29 11:28
 | HOME | Next »

上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。