Mathematicaで人参Beta-caroteneの量子力学的考察事始めVol.Ⅱ
Mathematicaで人参Beta-caroteneの量子力学的考察事始めVol.Ⅰからの続き1.2.5.一次元散乱 数値計算法は、束縛状態の場合と同様に、一次元散乱状態に対しても有効です。ここでは、定常状態に基づくアプローチと波動パケットに基づくアプローチという、全く異なる2つの散乱解析手法について説明します。1回照射法によるエネルギー固有状態1次元空間に局在化したポテンシャル障壁またはポテンシャル井戸があると仮定し、明確な運動量κを持つ右向きの粒子を入射させる。ポテンシャルは、確率Rで運動量 −κの粒子を反射するか、確率T = 1 − Rで運動量kの粒子を透過すると予想される。（簡略化の為、ポテンシャルには「段差」が無いと仮定する。従って、透過運動量は入射運動量と同じです。但し、この仮定を緩和するのは容易です。）要約すると、ポテンシャルの非自明な部分から離れた領域における波動関数は次のようになる。ψ(x) = (Aeiκx+ Be−iκ xas x →-∞ Ceiκxas x →+∞ (9)反射確率はR = |B|²/|A|²となります。 この波動関数はあらゆる場所でTISE（time-independent Schrödinger equation）時間依存しないシュレーディンガー方程式を満たす必要がある為、未知の振幅は、任意のκに対してTISEを数値積分する事ができる。コツは、右側から始める事です。右側では は単純な既知の形式を持ち、 は一度だけ"shoot"を実行し、ポテンシャルを左方向に積分して漸近領域に到達すると、係数AとBを実質的に読み取る事ができる。右側の適切な境界条件は、ψ= 1およびψ′ = κで、どちらも同じ （十分に大きな）xで、Cの大きさは任意に1とする。この方法を説明する為に、単純な長方形の障壁を使用する。高さは、任意のエネルギー単位で50です。In[]:v[x_] := If[x > 0 && x < 1, 50, 0];xMax = 3;Plot[v[x], {x, -xMax, xMax}, Exclusions -> None,AxesOrigin -> {-xMax, 0}]Out[]: TISEを解くには、適切な境界条件を設定したNDSolveを使用するだけで済む。複素波動関数はNDSolveにとって全く問題にならないが、それをプロットするには少し手間が掛かる。In[]:energy = 40;κ = Sqrt[2 energy];solution =NDSolve[{psi’’[x] == -2 (energy - v[x]) psi[x], psi[xMax] == 1,psi’[xMax] == I*κ}, psi, {x, -xMax, xMax}];Plot[Abs[psi[x]]^2 /. solution, {x, -xMax, xMax},PlotRange -> {0, All}, Filling -> Axis,ColorFunction -> Function[x, Hue[Arg[psi[x]/.solution]/(2Pi)]],ColorFunctionScaling -> False, PlotPoints -> 500]Out[]: 示された数値では、障壁の左側に大きな振幅の干渉パターンが得られ、古典的に禁制な障壁では指数関数的に減衰し、右側では（比較すると）振幅が無視できる程小さくなる。次に、障壁の幅を1ではなく0.1にしてもう一度試して見て下さい。すると、トンネル効果の確率が顕著になる事が分かる。左側の平均波振幅と比較して干渉縞が大きい程、反射確率は高くなる。この考え方を厳密にして式を導出するのは比較的簡単な演習です。R = |B|²/|A|²=((|ψ|max-/|ψ|min)/(|ψ|max+/|ψ|min))²ここで、「max」と「min」は、左側の領域における|ψ|のピークと谷を指す。Mathematicaスクリプトでは、これは次のようになる。In[]:max = FindMaximum[Abs[psi[x]] /. solution,{x, -.8 xMax, -.2 xMax}][[1]];min = FindMinimum[Abs[psi[x]] /. solution,{x, -.8 xMax, -.2 xMax}][[1]];reflection = ((max - min)/(max + min))^2Out[]:0.9996659893453111（Mathematicaの数値最適化関数に、やや恣意的な探索区間を指定した。これらの関数はそれぞれリストを返し、その最初の要素は極値の位置です。）1.2.6.Mathematicaに実装された関数を使用してTISE（time-independent Schrödinger equation）の解を可視化する。Ⅰ.シュレディンガー(Schrödinger)演算子を使用する。 シュレディンガー(Schrödinger)演算子をパラメータh，ポテンシャルVで指定する。In[]:h=1/2;V[x_]:=x2ζ=-h2*u''[x]+V[x]*u[x];細分化されたメッシュ上で小さい方から10個の固有値と固有関数を求める{vals,funs}=NDEigensystem[ζ,u[x],{x,-3,3},10,Method->{"PDEDiscretization"->{"FiniteElement",{"MeshOptions"->{MaxCellMeasure->0.05}}}}];valsOut[]:{0.10000016225988244`,0.30000113385045746`,0.5000040417329519`,0.7000101643270652`,0.9000207683332286`,1.1000371088172065`,1.3000604292934639`,1.5000919618072917`,1.7001329270156669`,1.9001845342669028`}hでスケールされ，それぞれの固有値でオフセットされた固有関数10個を可視化する。In[]:Show[Plot[Evaluate[h*funs+vals],{x,-3,3}],Plot[V[x],{x,-3,3}],PlotRange->{{-3,3},{0,2}},AxesOrigin->{-3,0},ImageSize->Large]Out[]:Ⅱ.SchrodingerPDEComponentを使用する。 SchrodingerPDEComponentは，ニュートンの第2法則に類似した非相対論的量子力学であり，量子力学における波動関数の時間発展を記述する。a.換算プランク定数でスケールされそれぞれの固有値でオフセットされた固有関数を可視化する。In[]:Show[Plot[Evaluate[1/10*funs+vals],{x,-3,3}],Plot[x^2,{x,-3,3}],Sequence[PlotRange->{{-3,3},{0,2}},AxesOrigin->{-3,0}],ImageSize->Large]Out[]:b.固有関数を可視化して固有値にラベルを付ける。In[]:Plot[funs,{x,-3,3},Sequence[PlotRange->All,AxesLabel->{x,Ψ},PlotLegends->en],ImageSize->Large]Out[]:Ⅲ.NDEigensystemを使用する。 NDEigensystemは，固有モードソルバとしても知られるものであるが，領域上で微分方程式の固有値と固有ベクトルを求める数値固有ソルバである．a.シュレディンガー(Schrödinger)演算子をパラメータ，ポテンシャル で指定する。In[]:h=1/2;V[x_]:=x2\[ScriptCapitalL]=-h2*u''[x]+V[x]*u[x];細分化されたメッシュ上で小さい方から10個の固有値と固有関数を求めるIN[]:{vals,funs}=NDEigensystem[\[ScriptCapitalL],u[x],{x,-3,3},10,Method->{"PDEDiscretization"->{"FiniteElement",{"MeshOptions"->{MaxCellMeasure->0.05}}}}];vals{vals,funs}=NDEigensystem[\[ScriptCapitalL],u[x],{x,-3,3},5];valsOut[]:{0.1,0.300001,0.500004,0.70001,0.900021,1.10004,1.30006,1.50009,1.70013,1.90018}b.hによってスケールされ，それぞれの固有値によってオフセットされた固有関数を可視化する。In[]:Show[Plot[Evaluate[h*funs+vals],{x,-3,3}],Plot[V[x],{x,-3,3}],PlotRange->{{-3,3},{0,5}},AxesOrigin->{-3,0},ImageSize->Large]Out[]:c.シフトを定義して固有値問題を解く。In[]:{en,funs}=NDEigensystem[{H[0],DirichletCondition[ψ[r]==0,True]},ψ[r],{r,0,100},4,Method->{"Eigensystem"->{"Arnoldi","Shift"->λ0},"SpatialDiscretization"->{"FiniteElement",{"MeshOptions"->{"MaxCellMeasure"->0.001}}}}];電子ボルトで入手された固有値を見る。In[]:en UnitConvert[Quantity[1,"RydbergEnergy"],"Electronvolts"]Out[]:{-0.850356eV,-1.51174eV,-3.40142eV,-13.6057eV}これらは最も負である固有値に対応する．しかし，順序は逆順であり，固有状態をソートし直す必要がある．SortByを使って固有状態をソートする。In[]:{en,funs}=SortBy[{en,funs}\[Transpose],Re]\[Transpose];確率密度をプロットする。In[]:Plot[Evaluate[funs^2],{r,0,30},PlotLegends->en UnitConvert[Quantity[1,"RydbergEnergy"],"Electronvolts"],PlotRange->All]Out[]:d.NDEigensystemを使って1次元時間独立シュレディンガー方程式を数値的に解く方法In[]:V[r_]=-(1/r)-(1.0415223038416566` E^(-0.9990999998788636` r))/rWith[{shift=10,d=30,n=5},{ev,ef}=NDEigensystem[{shift f[r]+V[r] f[r]-1/2 f''[r],DirichletCondition[f[r]==0,True]},f[r],{r,0,d},n,Method->{"SpatialDiscretization"->{"FiniteElement",{"MeshOptions"->{"MaxCellMeasure"->0.001}}},"Eigensystem"->{"Arnoldi"}}];evShifted=ev-shift]With[{n=4,d=30,amplitudes={-1,5,20,30}},Plot[Evaluate[Table[evShifted[[i]]+amplitudes[[i]] 1/r ef[[i]],{i,n}]],{r,0,d},PlotRange→{-4,4},Epilog→{Gray,Dashed,Table[Line[{{0,evShifted[[i]]},{d,evShifted[[i]]}}],{i,n}]}]]Out[]:{-1.33732035927342`,-0.18643394348523046`,-0.07104738501501373`,-0.03483728099528349`,-0.0011951714032178984`}e.NDEigensystemを使って一次元時間独立シュレディンガー方程式を解く。 この可視化の目的は，バリアの高さと井戸の長さが固有状態に与える影響を確認する事である。井戸の長さを伸ばすと各波動関数のエネルギーがどのように減少するかが分かる．一方で，バリアの高さを増すとエネルギーの固有値も大きくなる。ポテンシャルバリア，井戸の幅，あるいは固有状態の数の任意の値について固有状態を計算する関数を定義する。In[]:solver[V0_,d_,n_]:=NDEigensystem[{With[{V00=V0,d0=d},\[ScriptCapitalH]//Evaluate],DirichletCondition[Ψ[z]==0,True]},Ψ[z],{z,-4 d,4 d},n,Method->{"PDEDiscretization"->{"FiniteElement",{"MeshOptions"->{"MaxCellMeasure"->0.01 d}}}}];Manipulateを設定する。Manipulate[Block[{en,funs},{en,funs}=solver[V0,d,n];Row[{ListPlot[en->en,Sequence[Filling->Bottom,PlotLabel->"Energy eigenstates",ImageSize->Small,PlotRange->{0,1.}]],Show[Plot[With[{V00=V0,d0=d},V[z]//Evaluate],{z,-2 d,2 d},Sequence[PlotRange->All,AxesLabel->{"z","V0,Ψ"},Exclusions->None]],Plot[Evaluate[0.25 funs+en],{z,-2 d,2 d},PlotRange->All],ImageSize->Large]}]],{{V0,0.75},0.5,1},{{d,10},5,25},{{n,4},{2,3,4,5}},SaveDefinitions->True]Out[]:Ⅳ.Plotを使用する。ポテンシャル井戸内の固有関数IN[]:f[n_,x_]:=Abs[((1/Pi)^(1/4) HermiteH[n,x])/(E^(x^2/2) Sqrt[2^n n!])]^2Plot[Evaluate@Append[Table[f[n,x]+n+1/2,{n,0,7}],x^2/2],{x,-4,4},Filling->Table[n->n-1/2,{n,1,8}]]Out[]:1.3.TDSE（time-dependent Schrödinger equation） 時間依存シュレーディンガー方程式の数値解を出す。1.3.1.波束の時間発展散乱計算で定常状態を用いると、実際に何かが出入りしているようには見えない為、学習者にとっては混乱を招く可能性がある。より直感的な代替手段として、time-dependent Schr¨odinger equation時間依存シュレーディンガー方程式（TDSE）を、局所的な右向きの波束を初期状態として解く方法がある。それでは、その方法を説明する。自然単位系では、TDSEは次のようになる。i(∂ψ/∂t)=(-1/2)(∂2ψ/∂t2)+V(x)ψ (11)空間二階微分に対して「二階中心差分近似」を用いて、xを等間隔の値のリストに離散化する。i(∂ψ/∂t)=(-1/2)(ψ(x+dx)+ψ(x-dx)-2ψ(dx))/(dx)2+V(x)ψ(x) (12)ここで、dxは隣接するx値間の（非微小な）間隔です。この式を簡略化する為に、dx=1となる単位系を用いる。少し整理すると、TDSEは次のようになる。i(∂ψ/∂t)=-(ψ(x+1)+ψ(x-1))/2+(1+V(x))ψ(x) (13)時間を離散化する必要がある。つまり、∂ψ/∂tを離散近似で置き換える。最も簡単な方法は、次のように記述する事です。∂ψ/∂t≈ψ(x, t+dt) −ψ(x, t-dt)/dt (14)式13の右辺の各ψは、t（t+dtや他の時刻ではなく）で評価されるという前提で計算されます。しかし、この近似は未来時刻t+dtは考慮するものの、過去の時刻t-dtを考慮しない為、前方バイアスがかかっており、不正確です。中心差分近似を用いる方が遥かに正確です。∂ψ/∂t≈(ψ(x, t+dt) −ψ(x, t-dt))/(2dt) (15)これは近未来と直近の過去を対称的に扱うものです。この近似式を式13に代入し、ψ (x,t+dt)について解くと、完全に離散化されたTDSEのバージョンが得られる。ψ(x, t+dt) =ψ(x, t−dt) + i dt (ψ(x+1, t) +ψ(x−1,t) − 2(1+V(x))ψ(x, t) (16)式16は、我々が中心差分法と呼ぶTDSE（時間依存シュレディンガー方程式）の解法の核となる部分です。手順としては、離散化された格子上の全てのx点についてループ処理を行い、この式を用いてψ(x,t+dt)を計算し、時間を進めて同じプロセスを繰り返し、次の時間増分、更に次の時間増分、といった具合にψを計算していくす。あと少しだけ詳細を補足する必要がある。dtの値を大きくすると計算は速くなるが、dtを大きくし過ぎると、処理が不安定になる。小さな誤差が時間と共に指数関数的に増大し、最終的にψが巨大で意味不明な値になる。許容できる最大のdtの値はV(x)に依存する。我々の単位では、この最大値は常に0 .5未満です。そして実際には、dt=0.45が殆どの場合うまく機能する。式16は、計算対象の空間区間の両端では適用できない。何故なら、その場合、式16は両端からdxだけ離れたx値を参照する事になるからです。従って、通常のアプローチは、両端で常にψ=0と設定する事です。これにより、システムは実質的に無限に広がる正方形の井戸の中に埋め込まれる。そして、式16は内部の点でのみ使用する。 始める。勿論、初期波動関数ψ (x,0)を与える必要がある。しかし、式16を初めて使用する為には、時刻ゼロの1ステップ前のψ (x,-dt)も知る必要がある。簡単な解決策は、式14を用いて、dtを負の値として計算する事です。この近似はあまり正確ではないが、一度しか使用しない為、その影響は最小限です。それでは、中心差分法をMathematicaを用いて実装した例を示す。まず、空間格子のサイズを定義し、ポテンシャルエネルギー値の表を作成し、ポテンシャル関数をプロットする。In[]:xMax = 500;v = Table[If[x > 250 && x < 255, 0.05, 0], {x, 1, xMax}];ListPlot[v, Joined -> True]Out[]: Mathematicaのリストインデックスは0ではなく1から始まるため、xの範囲を1からxMax迄とするのが最も簡単です。我々の任意のポテンシャルエネルギー関数は、グリッド間隔が1距離単位の場合にうまく機能する幅と高さを持つ長方形の障壁です。次に、時間変数、時間ステップ、及び波動関数自体を初期化する（出力は無）。In[]:t = 0;dt = 0.45;x0 = 100;p0 = 0.25;a = 30;psi = Table[Exp[-(x-x0)^2/a^2] Exp[I p0 x], {x, 1, xMax}];psiLast = psiNext = Table[0, {x, 1, xMax}];For[x = 2, x < xMax, x++,psiLast[[x]] = psi[[x]]- (I*dt/2)*(psi[[x+1]] + psi[[x-1]] – 2(1+v[[x]])psi[[x]])];続編：Mathematicaで人参Beta-caroteneの量子力学的考察事始めVol.Ⅲ
