粒子法による量子波束の数値解析
著者 廣野 史明, 岩沢 美佐子, 狩野 覚, 善甫 康成
出版者 法政大学情報メディア教育研究センター
雑誌名 法政大学情報メディア教育研究センター研究報告
巻 34
ページ 46‑51
発行年 2019‑07‑18
URL http://doi.org/10.15002/00022803
法政大学情報メディア教育研究センター研究報告 Vol.34 （2019）
粒子法による量子波束の数値解析
Quantum Double-slit Interference using SSPH
廣野 史明１） 岩沢 美佐子２） 狩野 覚３） 善甫 康成３）
Fumiaki Hirono, Misako Iwasawa, Satoru Kano, and Yasunari Zempo
１）法政大学情報科学研究科情報科学専攻
２）法政大学理工学部
３）法政大学情報科学部
Symmetric smoothed particle hydrodynamics (SSPH) method is applied to the electric structure calculation in real space, which is described by the time dependent Schrödinger equation. In the process of time evolution, we employed the Bohmian expression. Since the propagation is represented by trajectories, we can directly recognize the propagation as the particle movement. We have been studying this technique to be applied to the analysis of interference and diffraction of the double-slit. We compared the solution by SSPH with that of analytical one. Both results have a good agreement, and SSPH has enough accuracy and stability to use practical applications. We report our present status of the development SSPH for the electronic structure calculation.
Keywords : Electric structure, Wave packet, Double-slit, SSPH, Trajectory, Bohmian
1. はじめに
電子と原子核は物質を構成する基本的な粒子であ り、物質の構造や性質は電子の運動状態によって決 定される。この電子状態を知ることによって、物質 の力学的・熱的・電磁気的・光学的性質などを理解 することができる。現在この電子状態を求める計算 では密度汎関数法が最も広く使われている。
これらの電子状態計算において実空間法が広く使 われている。実空間法は大規模並列計算に適してお り、特に多数の原子を含む大規模な系において利用 される。実空間法の例として実空間格子と有限差分 法を用いた方法や有限要素法を用いた手法がある。
複雑な形状に対しても、この実空間のメッシュを構 築すれば計算を行うことは可能ではある。しかしな がらこの手法では格子間の接続関係を考慮して格子 やメッシュを構築する必要があるため格子の構築そ のものに困難さや時間的コストが存在する。
一方、粒子法は格子やメッシュを一切用いること なく有限個の粒子（計算点）のみを用いて偏微分方 程式を離散化する手法である。前述の格子の接続関
係が不要であるため粒子間の接続関係を与える必要 がない。そのため非一様な計算点の分布を用いて局 所的に分解能を向上させることを容易に実現できる と期待される。また解析対象の形状に合わせて計算 点を配置することも可能である。粒子法以外のメッ シュフリー法などを電子状態計算へ応用する研究は 行われているが、粒子法を電子状態計算へ応用する 試みはほとんど行われていない。代表的な粒子法の 一つにSmoothed particle hydrodynamics（SPH）があ る。SPHは宇宙物理学[1]や電磁流体学[2]など様々 な分野で流体の方程式を解くために利用されている 手法である。我々は粒子法の中でも安定に計算でき るSymmetric SPH（SSPH）[3, 4]を使用した。我々 は電子状態計算に粒子法を用いる研究を行ってきた [5-7]。
本論ではBohm形式で表現された量子的な波束に 粒子法を適用し二重スリット実験を題材として時間 発展を行った[8-10]。二重スリット実験では特に、
干渉によってこの計算点が密になる領域と疎になる 領域が存在し計算の安定性が損なわれてしまう。ま た波動関数 ψ が非常に小さくなる点（節など）で
法政大学情報メディア教育研究センター研究報告 Vol.34
量子ポテンシャルの計算が発散してしまうことも知 られている。我々は計算点の追加と削除を解析中に 動的に行うことにより均等に計算点を分布させ精度 を保ちつつ安定に計算する手法を考案した。
本論の構成は以下の通りである。まず1章では概 要を述べた。2章では計算手法について述べる。3 章では数値計算の方法と結果を示す。最後に4章で は考察を行う。
2. 手法
2.1 Bohm 形式
Bohm形式で表現された量子波束を用いて電子状 態計算を行う。波動関数を振幅 R と位相 S の積に 分けると時間依存のSchrödinger方程式は以下の通 りに表される。
式(2)の波動関数を(1)に代入して整理すると以 下のように古典力学におけるHamilton–Jacobiの運 動方程式と連続の式となり量子効果は Q として表 される。本稿では ℏ= e = m =1とする原始単位系 を用いている。
式(3)は古典の運動方程式に量子ポテンシャル Q を加えた形である。Q は R = ecの置き換えにより以 下のように表される。
なおBohm形式では量子ポテンシャルQ の計算 上ゼロ割の問題が生じる。上記の置き換え（R = ec） をするために形式上この問題を回避できるが数値計 算上は無視できない問題である。
我々は式(3)と式(5)の微分 を 求めるのに差分ではなく粒子法を用いる。粒子法を 用いると計算点の分布に自由度を持たせることがで きるという特徴がある。粒子法ではこれらの値は一 般に弱形式から求められる。我々が用いる粒子法は SSPHである。次の説で詳しく説明するが、これは数 値的に安定に計算できるという特徴を持つ。
上述の式(3)~(5)を用いて時間発展させ計算点の
粒子の軌跡を解析していく方法をトラジェクトリー 法という。
2.2 SSPH
SSPHは空間に分布する、ある点の微分係数を周 辺に分布する計算点から重み関数を使用し求める手 法である。まず、ある空間関数 ψ(x)を考える。こ のとき地点r'での空間量はデルタ関数を使用し以 下の恒等式(7)で表される。
この地点r'にて重み関数を使用し近似する。重 みWはδ関数的な性質があれば何でもよい。我々
はWendland関数を使用した[11]。これは関数が値
を表す領域が有限な範囲に限られており計算量が少 なくてすむという特徴があるからである。
SSPHではこの恒等式 (7) とそのモーメントの
Talyor展開から を導く。具体的には式
(7)にk次のモーメントをかけて展開する。
これらの項をモーメントの順に並べていくと式 (10)のような一次連立方程式になる。式(9)の左辺 をP、右辺をMDとすると、最終的に以下のMD = P を解きDを求めればよい。また実対称行列Mの逆 行列M -1は必ず存在することが知られている[3]。
図 1 Wendland関数 Table 1 Wendland function
𝑖𝑖𝑖𝜕𝜕
𝜕𝜕𝜕𝜕 � � �
�𝑖�
2𝑚𝑚 𝜵𝜵�� �� ��
𝜓𝜓�𝒙𝒙𝒙 𝒙𝒙� � ��𝒙𝒙𝒙 𝒙𝒙�𝑒𝑒���𝒙𝒙𝒙��ℏ
�𝜕𝜕𝜕𝜕�𝒙𝒙𝒙 𝒙𝒙�
𝜕𝜕𝒙𝒙 � 1
2𝑚𝑚�𝛁𝛁𝜕𝜕��� ��𝒙𝒙𝒙 𝒙𝒙� � ��𝜌𝜌𝜌 𝒙𝒙𝒙 𝒙𝒙�
𝜕𝜕𝜕𝜕�𝒙𝒙𝒙 𝒙𝒙�
𝜕𝜕𝒙𝒙 � � � �𝜕𝜕1
𝑚𝑚 ��� � �
� � � 1 2𝑚𝑚
𝜵𝜵�𝑅𝑅�𝒙𝒙𝒙 𝒙𝒙�
𝑅𝑅�𝒙𝒙𝒙 𝒙𝒙�
� � �1
2��𝛻𝛻𝛻𝛻��� 𝛻𝛻�𝛻𝛻�
𝜓𝜓�𝒓𝒓�� � � 𝜓𝜓�𝒓𝒓�𝛿𝛿�𝒓𝒓 � 𝒓𝒓��
� 𝑑𝑑𝒓𝒓
� � 𝜓𝜓�𝒓𝒓�𝑊𝑊�𝒓𝒓 � 𝒓𝒓�, ℎ�
� 𝑑𝑑𝒓𝒓
𝑊𝑊�𝑟𝑟𝑟 𝑟𝑟� �
⎩⎪
⎨
⎪⎧ 3 4𝑟 �1 �
1 2 𝑟𝑟 𝑟�
�
�2 �𝑟𝑟 𝑟�
��5𝑟𝑟
2𝑟 � 1� 𝑟 �0 � 𝑟𝑟 𝑟 � 2�
𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟 𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟0𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟 �𝑟𝑟
𝑟 � 2𝑟� 𝑟𝑟𝑟𝑟
� �𝒓𝒓 � 𝒓𝒓� ���𝜓𝜓�𝒓𝒓�𝑊𝑊�𝒓𝒓 � 𝒓𝒓�� ��𝑑𝑑𝒓𝒓
� �1
𝑙𝑙𝑙
∂� 𝜓𝜓
∂𝒓𝒓��
𝒓𝒓�
�
���
� �𝒓𝒓 � 𝒓𝒓�����𝑊𝑊
� 𝑑𝑑𝒓𝒓 �� � ��1� � � ��
(1) (2)
(3) (4) (5)
(6)
(7)
(8)
𝜵𝜵𝜵𝜵、𝜵𝜵𝟐𝟐𝜵𝜵、𝜵𝜵𝜵𝜵、𝜵𝜵𝟐𝟐𝜵𝜵 (9)
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕���, 𝜕𝜕�𝜕𝜕
𝜕𝜕𝜕𝜕��
��
3. 数値解析
本手法ではBohm形式を用いて時間発展させた。
一方、粒子法は時間発展を伴う諸問題の効率的な解 法として近年注目されている。また電子状態の時間 的な変化という意味ではBohm形式を用いることが 再度注目されている。そこで電子状態の時間発展を 表す粒子法をBohm形式と組み合わせることを考案 した。その中で本稿では二重スリット実験に注目し た。なぜなら二重スリット実験では計算点で表され た波束同士が衝突し密になる中心の地点と、拡散し 疎になる外側の地点が存在するからである。そのた め粒子法では扱いづらい問題であり検証に最適であ るため採用した。
二重スリット実験では同一の波源から放出された 光子や電子が2つのスリットを通過し波束となり自 身と干渉し観測面に縞模様をつくる。これは光学的 には光路差から求めることができる。本稿では電子 を量子的な波として表しているため、その位相差が 観測面に縞模様の分布を示す。
図2のような解析を行う。シミュレーションでは 同一の波源から電子が放出されたものとしつつもス リットを通過した時点から実際の計算を開始した。
図2では、波源から放出された電子はスリットを 通過し観測面に縞模様の像を作る。まず2次元の自 由粒子の波動関数を以下のように定義する。スリッ ト通過直後（t = 0）を考える。またスリット方向を x軸、進行方向をy軸とする。
これを解析の初期状態の波束とし、2つのスリッ トからそれぞれ入射させる。また、X、Yがそれぞ れ独立であることからSSPHでは干渉の起こる x 方 向のみ解析すればよい。一方 y 方向は一定速度で干 渉のない進行波として考えることができる。そのた め y 方向は時間軸でもある。
また、このX(x,t), Y(x,t)はそれぞれ解析解を持つ。
これらが一次元の自由粒子のSchrödinger方程式を 満たすとすると解は以下の通りである。
この解析解を先ほどのSSPHの x 方向の干渉の数 値計算と比較を行った。数値計算は式(3)~(5)につ いてスリットでの波束(11)~(13)を初期条件として SSPHを用いて時間発展させた。
その結果、数値計算の波束の密度を三次元でプ ロットしたのが図3である。図の奥側から手前側に 向かって波束は進行する。2つのスリットを通過し た波束は干渉しながら拡散していく様子がわかる。
またこのときの量子ポテンシャルをプロットしたも のが図4である。
この図3の横軸はx座標、奥方向は時刻、縦軸は 図 2 二重スリット実験の状況
Figure 2 Model of double-slit interference
𝜓𝜓�𝑥𝑥𝑥 𝑥𝑥𝑥 𝑥� � ��𝑥𝑥𝑥 𝑥�𝑌𝑌�𝑥𝑥𝑥 𝑥�
𝑋𝑋�𝑥𝑥𝑥 𝑥� � 1
�𝜋𝜋𝜋𝜋����� ��� ��𝑥𝑥� 2𝜋𝜋���
𝑌𝑌�𝑦𝑦𝑦 𝑦� � 1
��𝜎𝜎�����
��� ��𝑦𝑦�
2𝜎𝜎�� � ���𝑦𝑦�
𝑋𝑋�𝑥𝑥𝑥 𝑥𝑥� � 1
�𝜋𝜋𝜋𝜋�����
1
�1 � 𝑖𝑖𝑖𝑚𝑚𝜋𝜋��𝑥𝑥�
��
� ��� � �𝑥𝑥�
2𝜋𝜋���1 � 𝑖𝑖𝑖𝑚𝑚𝜋𝜋��𝑥𝑥��
𝑌𝑌�𝑦𝑦𝑦 𝑦𝑦� � 1
��𝑚𝑚�����
1
�1 � 𝑖𝑖𝑖𝑚𝑚𝑚𝑚��𝑦𝑦�
��
� ��� ��𝑦𝑦�� 𝑖𝑖𝑚𝑚���2𝑘𝑘�𝑦𝑦 � 𝑖𝑘𝑘𝑚𝑚 𝑦𝑦��� 2𝑚𝑚���1 � 𝑖𝑖𝑖𝑚𝑚𝑚𝑚��𝑦𝑦� � (10)
(11) (12)
(13)
(14)
(15)
法政大学情報メディア教育研究センター研究報告 Vol.34
波束の密度を表している。図の奥から入射した2つ の波束が干渉しながら拡散する。
図4ではスリットの中心から放射状に量子ポテン シャルが大きい点が広がっていることが分かる。大 きい値（黄）は計算点に対して斥力であり低い値（緑）
には引力である。波束が進むにつれ量子ポテンシャ ルQに縞状に強弱が表れていることがわかる。こ の図4の縦軸は x 座標、横軸は時刻（ステップ数）
を表している。図の左側の部分（0 Step）がスリッ トの位置を表す。
4. 考察
4.1 計算点としての粒子
本手法は量子波束を計算点として表し粒子法を用 いて時間発展を行うことが特徴である。粒子法では 粒子的な取り扱いをするが、本手法では電子状態を 表す計算点である。図5は波束を計算点として表し た際の模式図である。それぞれの計算点は物理量を もつ。今回は量子的な波束を表しているので、ある
図 3 波束の強度 Figure 3 Magnitude of wave packet
図 4 量子ポテンシャル Figure 4 Quantum Potential
図 5 計算粒子と計算点 Figure 5 Smoothed particles and
numerical calculation points
計算点iは波動関数の値 ψi 、座標xi 、位相Si、計算 点の移動速度 vi 、仮想質量miをそれぞれ持つ。
4.2 解析解との比較
粒子法での計算結果を詳細に検討するため、解析 解との比較を行いつつ各時刻での波束をプロットし た。図6～図8で上部の図が粒子法であり下部の図 が解析解の式(14)(15)をプロットしたものである。
各時刻で解析解と一致していることが確認できた。
この図6の波束の密度 ρ（青）と量子ポテンシャ ル Q（緑）のプロット。上側は解析結果で下側は同 時刻での解析解を表す。横軸は x 座標、縦軸は波束
図 6 初期状態（t = 0 a.u.）
Figure 6 Initial condition（t =0 a.u.）
密度（青）とエネルギー（緑）を表している。2つ のピークはそれぞれスリットを通過した直後の波束 である。ピーク内部では量子ポテンシャルが正であ ることから斥力が働き波束を拡散させていくことが 分かる。一方、原点付近では量子ポテンシャルが負 であるため引力が働き、波束同士が衝突する方向に 力が働いている。この原点付近では密度が0に近い ため、計算が破綻しやすくなってしまう。その対策 手法については後で述べる。
次の図7では の変化に従いQも変化している 様子が分かる。衝突した波束が干渉し始めているこ とがわかる。2つ存在する量子ポテンシャルの負方 向への強いピークでは計算点同士が強く引かれ合う こととなる。一方、波束の周辺部分ではゆるやかに 拡散していく。
図8はさらに時間が経過し拡散と干渉が起こり5 つの密度のピークが確認できた。このピークの数は 初期条件によって変化する。この密度や量子ポテン シャルの形状が解析解と一致していることを確認で きた。
4.3 計算点の追加と削除
粒子法では計算の精度はある注目している粒子
（計算点）とその周辺の粒子の分布にある。周辺の 粒子が多ければ精度は上がるが接近しすぎると計算
上の破綻が生じる。一方、周辺の粒子が少ないと十 分な精度が得られない。この対策として計算点を動 的に追加・削除を行うようにした。また節となり計 算が破綻してしまう点では計算点を取り除くことに した。また今回の二重スリットの実験では計算点が 密な地点と疎な地点が存在する。計算点が密な地点 では削除をし、疎な地点では追加する処理を行い計 算の安定性を確保した。
図9では350,000 Stepごろから計算点が動的に追
加されていることが分かる。縦軸は x 軸、横軸はス テップ数を表す。この波束は拡散しているため全体 的に計算点同士の間隔が広がってしまっている。こ 図 7 波束の時間発展（t = 2 a.u.）
Figure 7 Time evolution of wave packets（t = 2 a.u.） 図 8 波束の時間発展と干渉（t = 4 a.u.）
Figure 8 Time evolution and interference of wave packets
（t = 4 a.u.）
図 9 動的な計算点数の変更
Figure 9 Dynamical change of the number of calculation points
𝜵𝜵�𝑅𝑅
法政大学情報メディア教育研究センター研究報告 Vol.34
の間隔が一定の値を超えたところで、その計算点同 士の中心位置に計算点を自動で追加することとし た。
5. 結論
二重スリットによる電子の干渉現象を題材とし て、粒子法とBohm形式でのトラジェクトリー法の 組み合わせによる解析の有用性を確認した。この解 析では解析解とほぼ一致するという結果を得た。ま た計算点を自由に追加や削除をすることでR≅0の 数値的な不安定な点を除去する手法も実装した。
今後は解析例を追加し本手法の検証をさらに進め ていく予定である。
[1] P. Laguna,
参考文献
“Smoothed Particle Interpolation”, Astrophys. J., Vol.439, pp.814-821, 1995.[2] G. Ala, et al., “Smoothed particle electromagnetics:
A mesh-free solver for transients”, J. Comput. Appl.
Math., Vol.191, pp.194-205, 2006.
[3] G. M. Zhang and R. C. Batra, “Symmetric smoothed particle hydrodynamics (SSPH) method and its application to elastic problems”, Comput Mech, Vol.43, pp.321-340, 2009.
[4] R. C. Batra and G. M. Zhang, “SSPH basis functions for meshless methods, and comparison of solutions with strong and weak formulations”, Comput Mech, Vol.41, pp.527-545, 2008.
[5] S. Soichiro and Y. Zempo, “Smoothed particle method for real-space electronic structure calculations”, J.
Phys: Conf. Ser., Vol.510, 012037, 2014.
[6] S. Soichiro, M. Toogoshi and Y. Zempo, “Device Simulation using Symmetric Smoothed Particle Hydrodynamics”, J. Phys: Conf. Ser., Vol.905, 012011, 2017.
[7] 杉本宗一郎, “粒子法を用いた実空間における電 子状態計算”, 情報科学研究科 修士論文, 2014.
[8] P. R. Holland, “Quantum Theory of Motion”, Cambridge University Press, 1993.
[9] 廣野史明, 岩沢美佐子, 狩野覚, 善甫康成“粒子法
による波束の干渉の数値解析”, 日本コンピュータ 化学会2018年度秋年会, 2P18.
[10] 廣野史明, 岩沢美佐子, 狩野覚, 善甫康成, “粒子 法とBohm形式による波束の干渉の数値解析”, 日 本応用数理学会2018年度年会6.
[11] H. Wendland, “Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree”, Adv. Comput. Math., Vol.4, p.389, 1995.
