粒子法による量子波束の数値解析

Abstract

The particle method has no restriction on the particle arrangement, where the calculation is performed. The purpose is, by using this feature, to develop the method to analyze the dynamics of the electronic state. In order to describe the time evolution of the electronic state, we have developed a new method to solve the time-dependent wave equation, using the Bohmian that is very compatible with the particle method. In this form, there is also numerical instability in the region, where the value of the wave function is very small. We have applied our technique to the wave packet dynamics on harmonic potential and the analysis of interference by double slits as a simple system with analytical solution and easy to compare accuracy, and realized that both results are in good agreement.

1 はじめに

粒子法はメッシュを用いない計算手法であるため計算点の配置の自由度が高い．粒子法を電子状態に適用した場合，高精度な計算が必要な領域へ集中的に計算点を配置することで効率的な計算を行うことが期待できる．粒子法の中でも比較的計算精度が高いSymmetric Smoothed Particle Hydrodynamics (SSPH) が開発されている [1, 2]．これを用いて固有状態を算出することを行ってきた [3, 4]．これらの場合，粒子の位置は任意である．一方で動的な時間依存の電子状態を解析するためには計算点である粒子を電子状態の変化とともに動かす必要がある．そこで我々は電子状態を表す波動関数をBohm形式で記述することにより [5]，この課題を解決することができた．解析的な解がある単純な系として，調和ポテンシャル上での波束の運動の解析および二重スリットによる干渉の解析結果について報告する．

2 計算手法

電子状態計算にSPHを用いて解析するとき，波動関数は積分形式で表現される．有限数のサンプリング点を用いて表現するため粒子法と呼ばれる．波動関数は次のような恒等式で表す．   

$\psi \left(r^{\prime }\right)=\underset{\text{Ω}}{\overset{}{\int }}\psi \left(r\right)\delta \left(r-r^{\prime }\right)dr \cong \underset{\text{Ω}}{\overset{}{\int }}\psi \left(r\right)W\left(r-r^{\prime }, h\right)dr$(1)

このδ(r)はKernel関数W(r,h) = W(|r|,h)として近似する．ここで h はsmoothing lengthでありWの広がりを表す指標である．Kernel関数としては，δ 関数的な性質があれば何でもよい．我々はWendland関数を採用している [6]．局所的に値を持つ関数であり数値計算上の負荷が少ない．

我々が用いているSSPHでは，式(1)とその k-次のモーメント式(2)についてTaylor展開を行うとMD=P の形で表される1次連立方程式となる．最終的に電子状態計算で必要になる ψ,∇ψ,∇²ψ などの項はMD = Pを解きDを求めればよい．   

$\begin{array}{c}\underset{\text{Ω}}{\overset{}{\int }}{\left(r-r_i\right)}^k\psi \left(r\right)W\left(r-r_i,h\right)dr\\ \cong \sum _{l=0}^m\left(\underset{\Omega }{\overset{}{\int }}{\left(r-r_i\right)}^{k+l}Wdr\right) \left(\frac{1}{l!}{\frac{{\partial }^l\psi }{\partial r^l}|}_{r_i}\right) \\ \left(k=0,1,\dots ,m\right)\end{array}$(2)

波動関数に合わせ粒子を動的に配置する場合，時間依存のSchrödinger方程式を$\psi \left(r,t\right)=R\left(r,t\right)\exp \left[iS\left(r,t\right)\right]$のBohm形式で表すと次のようなHamilton–Jacobiの運動方程式と連続の式となる．連続の式により粒子についてLagrange描像での解析が可能になる．   

$-\frac{\partial S\left(r,t\right)}{\partial t}=\frac{1}{2m}{\left(\nabla S\right)}^2+V\left(r,t\right)+Q\left(\rho ;r,t\right)$(3)

  

$\frac{\partial \rho \left(r,t\right)}{\partial t}+\nabla ·\left(\rho \frac{1}{m}\nabla S\right)=0$(4)

ここで Q は量子ポテンシャルであり次式で表され   

$Q=-\frac{1}{2m}\frac{{\nabla }^{2}R\left(r,t\right)}{R\left(r,t\right)}$(5)

量子効果を与えるものである．なお後述するが，数値計算上はゼロ割に注意する必要がある．

3 結果と考察

SPHの効果がよくわかる2次元調和ポテンシャル上での波束の運動と，波束の広がりと干渉が端的にわかる単純な系である二重スリットについて波束の動的な解析を行った．これらは何れも解析的な解がよく知られており [7]，精度の確認が容易である．

Figure 1は2次元調和ポテンシャルV(x²+y²)/2に基底状態のGaussian波束$\psi =\exp \left[-\left(x^2+y^2\right)/2\right]$を置き，その動的な解析を行ったものである．計算領域に粒子を配置し，そこでの波束の様子を動的に観測できる．それ以外に電子状態を計算する粒子(計算点)はない．調和ポテンシャルの中心にGaussian波束を置いた場合，時間発展をさせてもGaussian波束とそれを表す粒子の移動は生じず，安定した計算であることがわかる．

Figure 1.

 The particle arrangement for the Gaussian wave packet (dot), and the density is expressed in a color contour.

Figure 2は初期の波束に初速度を与え時間発展をさせたものである．Gaussianが変形することなく十分な時間にわたり時間発展を行うことが可能となった．粒子を等間隔な直交メッシュ位置に配置しても，またランダムに配置しても時間発展を行うことが可能であった．波束の移動に従って粒子が移動していくので，計算上の精度にも問題が生じない．これは粒子配置に制約がないという粒子法の特徴によるものである．これらの結果は解析解と非常によく一致している．

Figure 2.

 The oscillation of the Gaussian wave packet.

Figure 3は二重スリットによる波束の干渉を粒子法で解析したものである．この系では二重スリットの遠方上(この図では紙面の裏側)に1つの波源があると仮定し，スリットに到達したときは，ほぼ平面波とみなせるとする．その結果スリットからは波源を同一とするGaussianがそれぞれ放出されることになる．そこから更にGaussian波束が放出されるとして得られた電子密度を表示したものである．Figure 4は二重スリットから出てくる波束が時間とともに干渉していく様子がわかる．その干渉している波束の様子(左)を電子密度と量子ポテンシャルの時間変化を解析解(右)と比較したものである．非常によく一致していることがわかる．

Figure 3.

 The charge density profile for two Gaussian slits separated by 8.0 a.u. and the slit width is 4.0 a.u..

Figure4.

 The density profiles and quantum potentials at the time 0, 2.4, and 3.6 a.u. obtained by our particle method (left) and analytical solutions (right).

この系は粒子法で表現する場合に非常に難しい例である．式(3)のポテンシャルVがなく，波束干渉により生じる量子ポテンシャルQのみが存在する．二重スリットの外側で波束は拡散し，中央付近では波束が重なりQによる量子的な効果により干渉が生じていることがわかる．このような波束の干渉や波束がR(r,t) ∼ 0となる．このような場合Qの計算が数値的に不安定となる．これを解決するために次の2つの手法を導入している．

非常に大きな1つのGaussian波束とスリット放出された干渉する波束と線形結合させた波束について時間発展を行い，別途計算しておいた1つのGaussian波束の結果を差引くことで，干渉を安定に表現することができた．R ∼ 0近傍での数値的な不安定を取り除くことが可能となった．もちろん別途，大きな1つのGaussian波束の時間発展の計算を行っておく必要がある．

干渉が発生する場所付近では粒子が非常に疎あるいは密に偏ることが起こる．粒子法では粒子が接近し過ぎる場合や粒子間の距離が極端に広がる場合，微分の精度を電子状態計算レベルに保つことは難しい．このような状況を避け精度を上げるため，粒子の分布が疎な場合は粒子の追加を行い，粒子の分布が集中している場合には粒子の消去を動的に行う手法を導入した．これは粒子があくまで計算点であるという特徴を用いたものである．

4 まとめ

実空間における電子状態計算の空間離散化手法として粒子法の1つであるSSPHを適用する方法について示した．また動的な時間依存の電子状態を解析するためには計算点である粒子を電子状態の変化とともに動かす必要がある．そこで電子状態を表す波動関数をBohm形式で記述する手法を開発した．解析的な解がある単純な系として試験的に，調和ポテンシャル上での波束の運動の解析および二重スリットによる干渉の解析を行った．我々の手法と解析解は非常によく一致することが分かった．

謝辞

なお，この研究はJSPS科研費16K05047の助成支援を受けたものである．

参考文献

• [1]   R. C. Batra, G. M. Zhang, Comput. Mech., 41, 527 (2007). doi:10.1007/s00466-007-0209-3
• [2]   S. Sugimoto, Y. Zempo, J. Phys. Conf. Ser., 510, 012037 (2014). doi:10.1088/1742-6596/510/1/012037
• [3]   S. Sugimoto, Y. Zempo, J. Comput. Chem. Jpn., 13, 187 (2014). doi:10.2477/jccj.2014-0027
• [4]   K. Kitayama, M. Toogoshi, Y. Zempo, J. Phys. Conf. Ser., 905, 012011 (2017). doi:10.1088/1742-6596/905/1/012011
• [5]   C. L. Lopreore, R. E. Wyatt, Phys. Rev. Lett., 82, 5190 (1999). doi:10.1103/PhysRevLett.82.5190
• [6]   H. Wendland, Adv. Comput. Math., 4, 389 (1995). doi:10.1007/BF02123482
• [7]   e.g., L. I. Schiff, 3rd ed., "Quantum Mechanics", McGraw Hill (1968).