2次元六方晶系モデル結晶における非線形波の伝搬に関する計算機実験

伊多波 正徳, 湊 淳, 小澤 哲, 比企 能夫


Return

1 はじめに

ソリトンは、速度を変えずに波形を安定に保ちながら伝播する、線形の波に無い特異な性質を持つ。フォノンのような線形の波とソリトンのような非線形の波は、結晶中の種々の欠陥に衝突した際に、異なった反応が予想される。それぞれの波の性質を利用すること、つまり、物理的な応答の違いを積極的に利用することが考えられる。具体的には、フォノンおよびソリトンを新しいプローブとして応用することである。ゆえに、ソリトン問題についての研究に強い関心が持たれた。著者のグループは、ソリトンを応用した新しいプローブやデバイス等の開発を念頭に、それを実現するために必要な基礎データの収集に関する研究を行ってきた。具体的には、分子動力学(MD)法を使って、非線形質量−バネモデル結晶中の原子の励起機構についての研究を広範囲にわたって行った。
初めに、1次元格子において、入力変位の大きさに応じて、線形波およびソリトンの生成に成功し、また、格子欠陥との衝突実験を行った[1 - 3]。観測された線形波は、波数ベクトルを持ち、フォノンの性質に近い。次に、2次元正方格子について、多重ソリトンの生成および格子欠陥との衝突実験を行った[4, 5]。更に、3次元立方格子について、線形波およびソリトンの生成に成功した。前述の計算機実験では、相互作用は最近接原子についてのみ考慮された。しかしながら、2次元正方格子と3次元立方格子の実験は本質的には1次元的とみなすことができるため、第2近接原子による作用を考慮した実験を行った[6 - 9]。これら一連の研究から得られた知見のうちの一つは次のとおりである。第2近接原子の作用が大きい場合では、伝搬したソリトンのエネルギーの減少が最も少ないのは2次元正方格子であった。このことから、実際にソリトンを生成することに適した物質は面間相互作用が少ない準2次元的な結晶であることが予測された。そのような条件に近い物質であるグラファイト単結晶はソリトンを生成するための実験室実験で使用される試料の有力な候補であると考えている。本研究では、実験室実験を実現することに向けた基礎データの収集を行う目的で、グラファイト単結晶を想定した2次元六方格子モデル結晶について、分子動力学計算機実験を行った。
結晶の2つの方向について、第2近接原子の作用を変化させた場合におけるソリトンの速度とエネルギーについて調査した。そして、1次元格子、2次元正方格子、3次元立方格子および2次元六方格子について、ソリトンの速度とエネルギーを比較した。

2 モデルと計算

グラファイト単結晶は層状構造をもつ物質である。面間の原子間相互作用は、面内のそれに比べて十分小さいので、近似的にこれを無視して考えられるであろう。Figure 1に示すように2次元六方格子上に原子を配置し、それらが互いにバネでつながっている質量−バネモデルが採用される。個々の原子の運動は次の運動方程式で記述できる。

ここでmは質量、Ri, j は位置ベクトル、ijは原子を識別するための添え字、Tは時間、fi, j は位置エネルギーである。次に、原子間作用は中心力を考えており、原子の位置エネルギーはべき級数展開の4次の項まで考慮して次式で表現する。

ここでC(j)は原子間力の強さを決めるj次のパラメータ定数、Rp,qは位置ベクトル、pqは原子を識別するための添え字、(e)は基底状態(平衡状態)であること、またDi, jDp, qは原子の平衡位置からの変位を意味している。式(2)の表現は現実の物質のポテンシャル係数を代入できるという点で優れている。前述の運動方程式について、4次のルンゲ・クッタ法による数値積分を行い、原子の運動を時間の関数として追跡する。MD単位[10]が採用され、原子質量m = 1、格子パラメータL = 1000、時間T = 1である。また、原子間力定数についてはC(1)=1、C(2)= -0.1、C(3)= (C(2))2の関係を仮定する。理由は次のとおりである。3つの定数の比はC(1) : C(2) : C(3)= 1 : −10/L : 100/L2と書くことができる。原子ごとの位置エネルギーは、4次まで考慮して、次のように与えられる。

ここで、Δxは原子の変位である。3つの定数の比を考慮すると、次式となる。

ここで、結晶の高次の弾性定数について考えると、n次の弾性定数は次のように定義される。

2次、3次、そして、4次の弾性定数は歪が零の場合、
CIJ = C(1)L2, CIJK = -C(1) (10/L)L3
CIJKL = C(1) (100/L2) L4
である。ゆえに、CIJ : CIJK : CIJKL = 1 : -10 : 100である。これらの比は、高次の弾性理論の結論、また、実験的に評価された種々の結晶の高次の弾性定数について、大きさはオーダーで一致する[11]。
次に、式(2)のpqに関する和は、最近接原子対(nn : nearest neighbor)と第2近接原子対(nnn : next nearest neighbor)について考慮される。すなわち、

である。ここで、それらの相互作用の比k

と定義する。
DPは2つの方向から原子面に与えられる。
(A) x方向入力の場合、原子の数は86×32である。 DPを与える原子の座標は、Xが42と43、Yが0から31については[A]モード入力、Xが0、Yが0から31については[B]モード入力と呼ぶ。原子面の全ての原子にx方向のDPが与えられる。
(B) y方向入力の場合、原子の数は31×148である。 DPを与える原子の座標は、Xが0から31、Yが74と75については[A]モード入力、Xが0から31、Yが0については[B]モード入力と呼ぶ。原子面の全ての原子にy方向のDPが与えられる。
原子の位置は、格子パラメータLによって正規化されたXY座標値で表される。


Figure 1. Two - dimensional hexagonal model crystal used in the computer simulation.

3 結果と考察

3. 1 ソリトンの伝搬の様相

実験では、モデル結晶に平面状に入力変位DPを与え、縦波の伝搬を観測した。DPの方向に対し垂直な方向に周期境界条件が課された。平面状にDPが与えられるために、原子の変位の対称性が保たれる。そのため、DPの方向に垂直な方向の系のサイズは、DPの方向のそれに比べて小さく扱える。モデル結晶の端面は自由端である。
初めに、実験の様子を客観的に明らかにするために、波の伝搬の様子をFigure 2に示す。これはx方向についての[A]モードの場合である。横軸は原子位置Xである。縦軸は、原子のx方向の変位Dが横波的に任意単位で、また、全(運動+位置)エネルギーEが任意単位で示されている。時間経過はmds(MD steps)である。モデル結晶の中央にDPを与えた後、時間経過を経るに従い、正および負の方向に励起された波が伝搬する様子が示されている。この図は線形波とソリトンが混在している状態である。


Figure 2. Displacement D and energy E(arbitrary units) of atoms vs position X in long time for the x direction case.

非線形効果が現れないほど小さなDPを与えると、線形波のみが生成される。一方、それより大きなDPを与えると、線形波とソリトンが生成される。DPを種々変えた実験結果から、線形波の速度を知ることができる。つまり、速度を知ることによって、複数の励起の中から線形波の位置を特定することができる。また、ソリトンは、Dの波形がシャープになること、そして、エネルギーが空間的に集中することから判断される。結果の詳細は正方向に伝搬する波について述べる。
Figure 3は、結晶の中心にDPを与えた([A]モード)後、右へ進む波の先頭が結晶の長さの90%の位置まで伝搬した時の原子位置X(x方向の場合)またはY(y方向の場合)に対する原子のxまたはy方向の変位DEである。位置が130の原子についてのDおよびEは零であり、DおよびEのそれぞれの原点となる。これらのFigureは、振幅の最大値が収まるように描かれている。x方向の場合にYは9行、y方向の場合にXは21列で観測を行った。実験は、fi, j(nn)を一定とし、k = 0,0.1,...,1.0の範囲で行った。Figureはその一部を示している。k = 0の場合、DP = 10(MD units)である。DP(k≠0)は全エネルギーEtotal(k≠0)とEtotal(k = 0)が等しくなるように調整された。つまり、同じ入力エネルギーの条件のもとで、実験は行われた。


Figure 3. Displacement D and energy E(arbitrary units) of atoms vs atomic position x or y for an input pulse DP(k=0)=10(MD units) applied to two central atomic planes([A]mode input) in a Hexagonal crystal. The results for cases of different values of k are given.

x方向の場合。k = 0ではXが110付近に緩やかに変化するDの波形が現れた。また、その付近では、ほぼ平らなバックグランドにエネルギーのパルスが一つ存在しており、これはソリトンの特徴である。また、Xが80付近では線形波が現れた。次に、kが大きくなると、先頭を走るソリトンが2つに分裂して、エネルギーのピークの減少が認められた。
y方向の場合。先頭を走るソリトンのEが小さく表示されているが、これはエネルギーのピーク値で全体をスケーリングしているためである。全てのkの場合について、Yが65から80の範囲について、エネルギーの大きな局在ソリトンあるいは線形波が現れた。一方、この様な局在した波が現れた現象は、x方向については、k = 0.1とk = 0.2の場合において、現われた。


Figure 4. Displacement D and energy E(arbitrary units) of atoms vs atomic position x or y for an input pulse DP(k=0)=10(MD units) applied to one left atomic plane([B]mode input) in a Hexagonal crystal. The results for cases of different values of k are given.

Figure 4は、結晶の左端にDPを与えた場合であり、[B]モードについて示している。結果の傾向はFigure 3とほぼ同様である。先頭を走るソリトンは、[A]モードより長い距離を伝搬しているため、その特徴がより鮮明に現れた。ソリトン励起の特有な性質[1]であるDのステップ状の波形は、xとy方向ともに、k = 0で認められる。y方向の場合。[A]モードで見られた入力位置付近の大きな局在ソリトンや線形波はk = 0以外では現れなかった。

3. 2 1D、2Dおよび3Dとの比較

前述の[A]モードと同種の計算機実験は1次元格子(1D)、2次元正方格子(2D)および3次元立方格子(3D)について行われた[9]。1Dの場合。128個の原子からなる格子(Xは1から128)について、x方向の変位DP(k = 0) = 10(MD units)が中央(Xは63と64)の原子に与えられた。2Dの場合。128×128個の原子からなる正方格子(Xは1から128、Yは1から128)について、x方向のDP(k = 0) = 10(MD units)が中央(Xは63と64、Yは1から128)の原子列に与えられた。3Dの場合。128×16×16個の原子からなる立方格子(Xは1から128、Yは1から16,Zは1から16)について、x方向のDP(k = 0) = 10(MD units)が中央(Xは64と65、Yは1から16、Zは1から16)の原子列に与えられた。
Figure 5は、1D、2D、3D、2次元六方格子のx方向および2次元六方格子のy方向について、パラメータkに対する、先頭を走るソリトンの伝搬速度VnとそのエネルギーEnを正規化して示している。添え字nは正規化を意味する。次の現象が認められた。


Figure 5. Normalized propagation velocity Vn and normalized strength En of soliton vs parameter k for the case of [A] mode input for 1D, 2D, 3D, and two-dimensional hexagonal(x and y direction) crystals.

  1. kが増加するにつれ、Vnは早くなり、Enは減少した
  2. kが増加するにつれ、Enの減少が最も少ないのは2次元六方格子のx方向の場合であった。
kが増加することは物質が硬くなることを意味する。硬くなるほど音速が早くなることは妥当である。一方、kの増加によるEnの減少については次のような理由が考えられる。
  1. 先頭を走る波は、kが増加すると、より強く散乱される。そのために、エネルギーの分布が広がり、後追いソリトンが生成される。先頭を走るソリトンのエネルギーの一部は後追いソリトンの生成に費やされる。
  2. 局在ソリトンが大きなエネルギーを抱えた場合、先頭を走るソリトンのエネルギーが少なくなる。
  3. 線形波が大きなエネルギーを抱えた場合、先頭を走るソリトンのエネルギーが少なくなる。
x方向では、殆どのkで、(2)と(3)の要因による大きな影響が無かった。先頭を走るソリトンは主に(1)の要因による影響を受けるが、ソリトンを維持するために必要なエネルギーを保持していたために、Enの減少が少なかったと考えられる。
y方向では、ほぼ全てのkで、Yの中央付近に、エネルギーの大きな局在ソリトンや線形波が生成された。(2)と(3)の要因による影響を大きく受けたため、先頭を走るソリトンのエネルギーが弱くなった。
そのために、ソリトンを維持するために必要なエネルギーを保持できず、さらに、(1)の要因が影響してEnが大きく減少したと考えられる。Figure 5に示されている実験対象となった系は(1)から(3)の要因による影響を組み合わせてもっている。本研究で調査した範囲において、最もEnの減少が少なかったのは2次元六方格子のx方向であった。

4 まとめ

2次元六方格子質量−バネモデル結晶について、分子動力学計算機実験を行った。結晶の2つの方向について、ソリトンは原子の励起として観測された。ソリトンのエネルギーは、場合によっては、平らなバックグランド上に単一のパルスとして生成された。この様なパルスはプローブとして利用することにとって都合が良い。次に、先頭を走るソリトンの伝搬方向の違いによる際立つ差異は、[A]モードについて、ソリトンのエネルギーの減少の差として認められた。ソリトンのエネルギーの減少はx方向の方がy方向より少なかった。これは入力変位を与えた位置付近に生じたエネルギーの大きな局在ソリトンや線形波が影響したことが考えられる。次に、最隣接原子間ポテンシャルと第2近接原子間ポテンシャルの比に対する生成されたソリトンの伝搬速度とエネルギーを、1D、2D、3Dおよび2次元六方格子のそれぞれのモデル結晶について、比較した。共通して現れた現象は、第2近接原子間力が増加するにつれ、ソリトンの伝搬速度は早くなり、ソリトンのエネルギーは減少したことである。また、第2近接原子間力が増加するにつれ、ソリトンのエネルギーの減少が最も少ないのは2次元六方格子のx方向の場合であった。
本論文では、2次元六方格子について、ソリトンが生成されることが示された。現実の物質でソリトンを生成する場合、準2次元的な物質であるグラファイトは試料の有力な候補であると考えている。

参考文献

[ 1] S. Ozawa and Y. Hiki, Jpn. J. Appl. Phys., 34, 2590 (1995).
[ 2] S. Ozawa and Y. Hiki, Physica B Condensed Matter, 219&220, 473 (1996).
[ 3] S. Ozawa and Y. Hiki, Jpn. J. Appl. Phys., 35, 3205 (1996).
[ 4] H. Tanaka, S. Ozawa and Y. Hiki, Jpn. J. Appl. Phys., 36, 2964 (1997).
[ 5] H. Tanaka, S. Ozawa and Y. Hiki, Jpn. J. Appl. Phys., 37, 2794 (1998).
[ 6] S. Ozawa, R. Komuro and Y. Hiki, Physica B Condensed Matter, 263-264, 90 (1999).
[ 7] Y. Hiki, R. Komuro, A. Minato, S. Ozawa, Jpn. J. Appl. Phys., 38, 3020 (1999).
[ 8] Md. M. A. Joarder, A. Minato, S. Ozawa, and Y. Hiki, Jpn. J. Appl. Phys., 39, 2941 (2000).
[ 9] Md. M. A. Joarder, A. Minato, S. Ozawa, and Y. Hiki, Jpn. J. Appl. Phys., 40, 3501 (2001).
[10] S. Ozawa, Md. M. A. Joarder, A. Minato,Y. Hiki, Journal of Alloys and Compounds, 10, 97 (2000).
[11] Y. Hiki, Ann. Rev. Mater. Sci., 11, 51 (1981).


Return