修士学位論文
2成分レナード・ジョンズ
コロイド分散系の
ガラス転移近傍における
ブラウン動力学シミュレーション
平成 19 年度
(平成 20 年 1 月 31 日提出)
東北大学大学院工学研究科
ナノメカニクス専攻
木村祐人
Brownian-dynamics Simulation of Binary Colloidal Suspensions
with Lennard-Jones Potentials near the Glass Transition
Yuto Kimura
Abstract
The glass is one of the familiar materials for us. It is known as a solid state with disordered
structure. The process that glass forming melt changes into a glassy state is called glass
transition and observed in various systems, for example, soda-lime glass, polymer blend,
colloidal suspension and metallic glass. Although glass is used in the widely-rangig field
of technologies, the mechanism of glass transition is not fully understood. As the glass
transition point approaches, the viscosity of the system increases rapidly, and finally
the system loses the fluidity. Near the glass transition, several time scales appear and
the relaxation process becomes complex. β- and α-relaxation is typical example. We can
observe the change of relaxation process in multi-step decay of self-intermediate scattering
function. Recently, the glass transition has attracted a remarkable attention in the field
of statistical physics and material sciences. It is obvious that the study of the glass
transition has significant meanings scientifically and technologically.
In the present thesis, we study both one component system and binary mixture of colloidal suspensions with Lennard-Jones potential by employing an extended scale Browniandynamics simulation. The control parameter of our system is termperature T and simulations were executed at various termperature ranged from 5.00 to 0.44 in a dimensionles
value. In order to research the static and dynamical behavior of the system, some physical quantities such as effective pressure, radial distribution function, static structure
factor, coarse grained density, mean square displacement, mean fourth displacement,
non-Gaussian parameter and self-intermediate scattering function were investigated into
details. From those results, we observed liquid-crystal transition in the one component
system. On the other hand, the binary mixture showed characteristic behaviors of supercooled liquid. Large number of particles suppressed the finite size effect and resulted
the logarithmic growth in the mean square displacement, and logarithmic decay in the
self-intermediate scattering function. We analyzed the mean square displacement, a typical physical quantity which show the dynamical behavior of the system, using the mean
field theory (MFT) proposed by Tokuyama. The simulation results showed the very well
match with MFT.
From the universal point of view of MFT, we compared our results with the data obtained from molecular dynamics simulations of binary glass forming mixture with same
interaction. The static properties, such as pressure and radial distribution function of the
two system were completely identical. Additionally, taking into account the short-time
self-diffusion coefficient DSS whose effect for the dynamics was neglected in the present
simulation, we obtained the same behavior of not only the long-time diffusion coefficient but also the time development of mean square displacement after early β-process.
These facts indicate the importance of the hydro dynamic interaction for the dynamics of
colloidal suspensions near the glass transition and validity of the the standpoint of MFT.
目次
第1章
1.1
1.2
1.3
1.4
緒言
研究背景 . . . .
これまでの研究
本研究の目的 .
論文の構成 . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
第 2 章 モデル
2.1 ブラウン運動とランジュバン方程式 .
2.2 相互作用 . . . . . . . . . . . . . . . .
2.2.1 流体力学的相互作用 . . . . .
2.2.2 レナード・ジョンズ相互作用
2.3 システムの導入 . . . . . . . . . . . .
2.4 時空間スケール . . . . . . . . . . . .
第3章
3.1
3.2
3.3
3.4
3.5
3.6
計算手法
計算式 . . . . . . . . . .
無次元化 . . . . . . . . .
差分式の導出 . . . . . .
周期境界条件 . . . . . .
初期配置 . . . . . . . . .
計算プログラムの高速化
3.6.1 粒子登録法 . . .
3.6.2 サブセル法 . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
第 4 章 解析方法
4.1 準備 . . . . . . . . . . . . . . .
4.2 密度の導入 . . . . . . . . . . .
4.2.1 密度と密度揺らぎ . . . .
4.2.2 粗視化された密度揺らぎ
4.3 動径分布関数と静的構造因子 .
4.4 有効圧力 . . . . . . . . . . . . .
4.5 平均二乗変位 . . . . . . . . . .
4.6 ノンガウシアンパラメタ . . . .
4.7 散乱関数 . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
i
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1
1
3
5
5
.
.
.
.
.
.
7
7
8
8
8
10
11
.
.
.
.
.
.
.
.
13
13
13
15
17
17
19
19
21
.
.
.
.
.
.
.
.
.
23
23
24
24
25
27
29
30
31
33
4.7.1 van Hove 関数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.7.2 自己中間散乱関数 . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.8 平均場理論 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
第 5 章 計算結果と考察
5.1 単成分系 . . . . . . . . . . . . . . . . . . .
5.1.1 計算条件 . . . . . . . . . . . . . . .
5.1.2 有効圧力 . . . . . . . . . . . . . . .
5.1.3 動径分布関数 . . . . . . . . . . . .
5.1.4 平均二乗変位 . . . . . . . . . . . .
5.1.5 ノンガウシアンパラメータ . . . . .
5.1.6 自己中間散乱関数 . . . . . . . . . .
5.2 2成分系 . . . . . . . . . . . . . . . . . . .
5.2.1 計算条件 . . . . . . . . . . . . . . .
5.2.2 有効圧力 . . . . . . . . . . . . . . .
5.2.3 動径分布関数と静的構造因子 . . .
5.2.4 粗視化された密度揺らぎ . . . . . .
5.2.5 平均二乗変位 . . . . . . . . . . . .
5.2.6 ノンガウシアンパラメタ . . . . . .
5.2.7 自己中間散乱関数 . . . . . . . . . .
5.2.8 平均場理論による解析 . . . . . . .
5.2.9 時間スケール . . . . . . . . . . . .
5.3 分子動力学法による計算結果との比較 . .
5.3.1 圧力と動径分布関数 . . . . . . . .
5.3.2 自己長時間拡散係数と平均二乗変位
第6章
6.1
6.2
6.3
6.4
6.5
6.6
6.7
付 録A
A.1
A.2
A.3
A.4
A.5
結言
単成分系と2成分系の比較 . . . . . . . .
緩和過程の複雑化 . . . . . . . . . . . . .
A 粒子と B 粒子の比較 . . . . . . . . . .
平均場理論による解析 . . . . . . . . . .
特性時間の変化 . . . . . . . . . . . . . .
分子動力学法による計算機実験との比較
最後に . . . . . . . . . . . . . . . . . . .
ランジュバン方程式の性質
ガウス積分 . . . . . . . .
一様乱数と正規乱数 . . .
初期配置の作製法 . . . . .
平均場理論について . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
ii
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
33
33
36
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
39
39
39
39
40
41
42
43
44
44
44
44
48
49
51
55
58
60
64
64
65
.
.
.
.
.
.
.
71
71
72
72
72
73
73
74
.
.
.
.
.
75
75
77
78
79
81
謝辞
87
iii
iv
第 1 章 緒言
1.1
研究背景
我々の住む世界には様々な物質が存在している. それらは周囲の環境の温度によって
大きくその性質を変える. 低温から温度を上昇させると固体, 液体, 気体, さらにはプラ
ズマ状態が得られる. 有史以来, 我々はこれらの物質状態を様々な形で生活に利用してき
た. その中で固体に注目すると, 身近なものでは鉄やアルミニウム, 銅などの金属, シリコ
ンなどの半導体, ポリエチレンやポリプロピレンなどのプラスチック, さらには木材やセ
ラミックなどが挙げられる. もう一つ, 身近な固体の例としてガラスを挙げることが出来
る. ガラスは窓ガラス等の建築材, コップ等の生活用品, 瓶等の保存容器, カメラのレンズ
等の光学部品, 白熱電球や蛍光灯などの照明装置の材料として利用されて来た [1]. また,
近年では液晶ディスプレイに代表されるディスプレイの基板, 光ファイバー等の情報伝達
装置, 光学デバイスの材料として, さらに利用の場を広げつつある.
前置き無しにガラスという言葉を用いたが, 我々がガラスと聞いて思い浮かべるのは窓
ガラスやガラス瓶であろう. これらの材料は珪酸 (Na2 O-CaO-SiO2 ) を主原料としたソー
ダ石灰ガラスである [1]. 伝統的なガラスの製法は高温で原料を融解し, それを冷却する
というものである. 冷却が進むにつれてガラス形成物質の液体はある温度の付近で粘度を
急激に増し, 流動性を失い, やがて機械的に固体と見なせる状態となる. これがガラスで
あり, 物質がガラスに変化する現象はガラス転移と呼ばれている. ソーダ石灰ガラスの他
にも一般的に用いられているガラスには様々な組成を持ったものが存在するが, ほとんど
は酸化物を混ぜ合わせた組成を持っており, 酸化物ガラスと呼ばれる. 特に, SiO2 , B2 O3 ,
P2 O5 , GeO2 等は単体でガラス化することが知られており, ガラス形成酸化物と呼ばれる.
ここで SiO2 のガラスである石英ガラスは, SiO2 の結晶である水晶と同じ組成を持ってい
る. この例からも分かるように, ガラスは組成によらない物質の状態の一つであると言う
ことが出来る. 水晶と石英ガラスを比較すると, 水晶は原子配列が規則的であるのに対し,
石英ガラスは原子が不規則に並んでいるという違いがある. この不規則な構造というの
1
第 1 章. 緒言
1.1 研究背景
は, ガラスの一つの特徴である.
ガラス溶液は固化する付近の温度で急激に粘性が増加するが, 同じような現象はガラス
以外の物質でも見られる [2]. ガラス転移点近傍の粘性の挙動に注目すると, 1013 ポアズ
の温度をガラス転移温度と便宜的に定義することが出来る. この方法で Angell は様々な
物質のガラス転移温度付近での粘性の挙動を整理した [3]. また, ガラス転移近傍ではい
くつかの特性時間が表れ, それが急激に変化して緩和過程が複雑になることが知られてい
る. α 緩和と β 緩和はその顕著な例である. 緩和過程の複雑化は自己中間散乱関数の多段
階緩和等に見ることが出来る. ガラス転移現象は酸化物ガラスのみならず, 近年注目を集
めている金属ガラス, 高分子系, 糖等の有機分子, さらにはコロイド分散系等, 様々な物質,
系で観測されることが知られている. その中でもコロイド分散系は原子・分子系と比較し
て時空間スケールが非常に大きく, ガラス転移現象が系や物質によらない普遍的な現象で
あることを伺わせる.
コロイド分散系とは分散媒となる物質中に, 10nm から 1µm 程度の大きさを持った異な
る物質が分散している物質状態である [4, 5]. これは我々の身近に多く存在している. 特
に分散媒が液相, 分散質が固相の場合をサスペンションと呼ぶ. メタクリル酸メチル樹脂
(poly-(methyl methacrylate), PMMA) の表面に poly-(12-hydroxy stearic acid) の薄い層
を持つ球形コロイド粒子を用いて中性コロイド分散系を実現することが出来る. 中性コロ
イド分散系において液体状態と結晶状態が存在し, さらにはガラス状態も観察されること
が報告されている [6]. このコロイド粒子のガラス状態はコロイドガラスと呼ばれる. コ
ロイドガラスの構成粒子は数 µm 程度であり, ガラスの構成要素である原子の大きさと比
較すると 1000 倍程度大きい. それに伴い, 系の時間変化, すなわちダイナミクスも原子・
分子系に比べて遅くなる. これらの特徴から, コロイド分散系ではガラス転移点付近で系
のダイナミクスの直接的な観察が比較的容易に行われる. 実際, 各種散乱実験や中性コロ
イド分散系を共焦点顕微鏡等を用いて観察する実験が行われており, 多くの有益な情報が
提供されている [7, 8, 9, 10, 11].
近年, ガラス転移現象が統計物理学の分野で大きな注目を集めている [12, 13, 14, 15].
その契機となったのが先に述べた中性コロイド分散系におけるガラス状態の発見 [6], そ
してモード結合理論のガラス転移への応用である [16, 17]. コロイドガラスの発見はガラ
ス転移現象の持つ普遍性を示し, 物理学の研究対象としての価値を高めたと言えるだろ
う. また, モード結合理論の登場はそれまで部分的な時間スケールでしか扱うことの出来
2
第 1 章. 緒言
1.2 これまでの研究
なかったガラス転移点近傍の物質のダイナミクスを統一的に扱うことが出来る可能性を
示した [12].
ガラス転移現象の研究には計算機実験も多く用いられる. バルク金属ガラスの応用研究
が進められているということもあり, 金属ガラスを模擬した計算機実験も積極的に行なわ
れている [18]. 理論的研究や実験と比較すると, 以下のような計算機実験の利点を挙げる
ことが出来る. 一つ目は系のモデル化を行ない, 現象の本質を抜き出すことが出来るとい
う点である. ガラス転移現象には原子・分子, あるいはコロイド粒子等, 相の構成粒子の
相互作用の多体相関をどのように取り扱うかという問題が横たわっている. しかし, 実際
の系では相互作用自体が複雑であり, 全てを厳密に取り扱うことは困難である. 二つ目が,
モデル化した系について第一原理から計算することにより, 理論の検証を行なうことが出
来るという点である. ここで言う第一原理とは量子力学のことではなく, 統計的手法を用
いて理論を構築する際の出発式という意味で, コロイド分散系ではランジュバン方程式が
それに当たる. 例えば, モード結合理論に関しては近年詳細な数値計算が可能となり, そ
れを計算機実験の結果と比較することで問題点が明らかになりつつある [19, 20, 21]. 三
つ目が, 計算過程において系に関する情報を全て取り出すことが出来るという点である.
これにより, 実験に比較して詳細な解析を行なうことが出来る. さらには, 粒子ごとの大
きさ, 密度, 系の温度, 初期状態等を意図した通りに設定出来るのも計算機実験の特長で
ある. また, コロイド分散系に関して言えば, 実験で懸念されるコロイド粒子の沈殿や凝
集が起こらないという点も利点と言える.
以上のように, ガラス転移現象は様々な物質, 分野に見られる普遍的な現象であり, そ
の研究は工学的, 物理的に意義のあるものである. 本質的なガラス転移現象の解明には画
期的な理論的アプローチが必要であるが, その開発には計算機実験の提供する情報が実験
と同等に必要不可欠なものとなっている.
1.2
これまでの研究
ガラスは五千年以上昔から人類に利用されて来たことが知られている [1]. 長い歴史の
中で洗練されて来た成分, 製法等は現在のものにも通じるものがある. しかしガラス転移
現象の理論的な取り扱いに関しては, ここ数十年で発展したと言っても過言ではない.
ガラス転移点近傍の動的特性の研究は, モード結合理論がガラス転移現象に応用されて
以来発展を遂げた [16, 17]. 初期のモード結合理論はガラス転移点で長時間自己拡散係数
3
第 1 章. 緒言
1.2 これまでの研究
が 0 になるという特異性を含んでいたが, 改良が施され現在に至っている [22, 23, 24]. コ
ロイド分散系におけるモード結合理論は Szamel らによって構築された [25]. Fuchs らに
開発されたアルゴリズム [26] に基づいて近年 Voigtmann ら, および Flenner らによりモー
ド結合理論の詳細な数値解が計算され, 対応する系における計算機実験と比較が行われた
が [19, 20, 21], 近年徳山によってモード結合理論は β 緩和過程を正しく記述出来ていない
ことが示唆されている [27].
ガラス転移点近傍での動的特性をよく説明する理論として, 徳山による平均場理論が知
られている. 1994 年に徳山らはコロイド分散系について詳細な理論計算を行い, 短時間自
己拡散係数, および長時間自己拡散係数の振る舞いを予測した [28]. 後に平均二乗変位の
時間発展方程式を平均場理論として求め, コロイド分散系のみならず原子分子系にも拡張
が行われた [29]. 平均場理論はコロイド分散系, 原子分子系と異なるその動的性質を説明
出来ることが報告されている [30].
Pusey らは半径に 5%の分散を持つ中性コロイド分散系を用いて液体や結晶等の相変化,
さらにはガラス状態が実現されることを発見した [6]. 以来中性コロイド分散系は原子分
子系の模擬系として注目を浴び, 自己中間散乱関数等の様々な物理量が測定され, ガラス
転移の研究に必要な情報が提供されてきた [7]. コロイド分散系を研究する大きな利点は,
原子分子系に比較して時空間スケールが大きいため実験において粒子の運動を直接追跡
出来るという点である. この利点を生かして, 過冷却液体の特徴である不均一な動的特性
が盛んに研究されている. Weeks らは共焦点顕微鏡を用いてコロイド粒子の三次元位置
情報の計測を実現し, 過冷却状態において粒子が協調運動を示すことを報告した [8]. 彼ら
は後に過冷却状態において移動度の大きいコロイド粒子が移動度の小さい粒子と異なり,
線状のクラスタを形成していること等を実験で観測している [9]. 中性コロイド分散系以
外にもモデル系が考案され, 実験が行われている. König らは磁性コロイドを重力により
サンプルの気液界面に安定させた擬二次元磁性コロイド分散系を用いて流体力学的相互
作用の少ないモデル系を構築し, 過冷却状態における協調運動等を観察した [10, 11].
希ガス分子の相互作用のモデルとして知られているレナード・ジョンズ相互作用は, ガ
ラス転移現象の研究においても積極的に用いられている. Kob と Andersen は2成分金
属ガラスである Ni80 P20 の相互作用をレナード・ジョンズ相互作用を用いてモデル化し,
1000 粒子の計算機実験で自己中間散乱関数, α 緩和時間, 長時間自己拡散係数等を測定し
た [31, 32, 33]. 以降, 彼らが用いたパラメータでその後多くの計算機実験が行われている
4
第 1 章. 緒言
1.3 本研究の目的
[20, 21, 34, 35, 36].
1.3
本研究の目的
本研究ではレナード・ジョンズ相互作用を持った2成分コロイド分散系について, 10976
粒子の大規模な計算機実験をブラウン動力学法を用いて行う [20, 21]. 参考とするため, 一
成分のみが存在する単成分系についても計算を行い, 比較を行う. この系を選んだ理由は
以下の通りである.
• 相互作用が単純であり, 安定した過冷却状態を示すことが知られている.
• 2成分レナード・ジョンズコロイド分散系において粒子数が 10000 程度の大規模な
計算は未だ行われていない.
• 同一の相互作用を用いた分子動力学法による計算が行われており, 計算結果を比較
することによりガラス転移現象の普遍性について考察出来る.
計算を通じて有効圧力, 動径分布関数等の静的な物理量, 平均二乗変位やノンガウシア
ンパラメータ, 自己中間散乱関数等の動的な物理量を計測することにより, 系の静的, 動的
な特性の温度依存性を詳細に調査する. また, 平均二乗変位に着目し, 徳山により提案さ
れた平均場理論を用いて, 長時間自己拡散係数, 系の特性時間等について解析を行う. さ
らにガラスは系によらない普遍的な現象であるという視点に基づき, 得られた結果を同一
の相互作用を持った分子動力学法による計算 [37] と比較し, 考察を行う.
1.4
論文の構成
本論文は 6 章から構成される. 次章では本研究における最も基礎的な方程式であるラ
ンジュバン方程式について簡単に触れた後, 計算機実験に用いた具体的なシステムを導入
する. 粒子間相互作用について, 流体力学的相互作用, レナード・ジョンズ型相互作用の
検討を行い, その後注目する時空間スケールを考慮して実際に計算に用いる方程式を導出
する.
その次の 3 章では主に計算機実験を行う上での手法を述べる. まず, 2 章で求めた方程
式について無次元化, 差分化を行う. 次いで, 境界条件, 初期条件 (初期配置), 計算プログ
ラムの高速化の手法について, 実装方法を念頭に置いて述べる.
5
第 1 章. 緒言
1.4 論文の構成
4 章では解析手法を導入する. 本研究では様々な物理量を計算機実験により計測するが,
その計算方法がここに記してある. 最初に粒子密度を導入し, それを用いて他の物理量を
導く, という手順で記述した.
5 章では計算結果と考察を述べる. 各種物理量について単成分系, 2成分系それぞれの
結果を示し, 比較と考察を行っている. 前半では静的な物理量を議論し, 後半ではダイナ
ミクスに関連した動的な物理量の結果を示す.
6 章では本論文の総括を行う.
6
第 2 章 モデル
2.1
ブラウン運動とランジュバン方程式
溶液中に分散したコロイド粒子はランダムな運動, ブラウン運動をすることが知られて
いる [38]. Langevin は 1902 年にブラウン運動のモデルとしてランジュバン方程式を提案
した. コロイド粒子の質量を m, 時刻 t における位置ベクトルを X(t), 速度を V (t) とす
ると, 粒子の運動は以下の運動方程式で記述される.
d
X(t) = V (t),
dt
d
m V (t) = −γV (t) + R(t).
dt
(2.1)
(2.2)
ここで (2.2) 式の右辺はどちらも周囲の溶媒から受ける力である. 右辺第一項は速度に比
例した粘性抗力であり, 係数 γ はストークスの法則より γ = 6πηa となる. ここに η は溶媒
の粘性率, a は粒子の半径である. 第二項は溶媒の分子の衝突に起因した揺動力と呼ばれ
るランダム力であり, 溶媒分子の運動の特性時間 tm 程度の時間スケールで変動する. R(t)
の時間変化はコロイド粒子の運動に比べて非常に速いため, ブラウン運動を考える際には
同時刻のみ相関を持つと見なせる. これはデルタ関数を用いて表現出来る. 実際, R(t) は
次の関係を満たす.
hR(t)R(t0 )i = 2γkB T δ(t − t0 )1.
(2.3)
ここで kB はボルツマン定数, T は溶媒の温度, 1 は単位テンソルである. (2.3) 式は溶媒の
持つ熱エネルギーとブラウン運動の駆動力である揺動力の間の関係を示しており, 揺動散
逸関係と呼ばれている.
ランジュバン方程式の特徴的な点は, 式 (2.2) の右辺に溶媒を流体と見なした第一項と
溶媒を粒子と見なした第二項を同時に含んでいる点である. 粒子が大きすぎては第一項
の効果が大きいためランダムな運動は起こらず, 逆に粒子が小さすぎると第二項の効果が
大きくなり運動が適度に緩和しなくなってしまう. ブラウン運動は二つの項が同程度の
時に観察される運動であると言える.
7
第 2 章. モデル
2.2
2.2.1
2.2 相互作用
相互作用
流体力学的相互作用
1994 年に徳山と Oppenheim は流体力学的相互作用を考慮したランジュバン方程式を提
案した [28, 39].
d
Xi (t) = Vi (t),
dt
N
N
X
X
d
m Vi (t) = −γ
Gij Vi (t) +
Fij (t) + Ri (t)
dt
j
j6=i
(2.4)
(2.5)
揺動散逸関係は,
hRi (t)Ri (t0 )i = 2γkB T δ(t − t0 )δij 1
(2.6)
となる. ここで Fij は粒子 i が粒子 j から受ける力, δij はクロネッカーのデルタである.
Gij は繰り込まれた摩擦係数テンソルであり,
Gij = [(1 + g)−1 ]ij
(2.7)
と書ける. また, gij は流体力学的相互作用を表す Oseen テンソルであり,

 0
µ
¶ for i = j,
gij =
3 ai
Xij Xij

1+
for i 6= j
4 Xij
Xij Xij
(2.8)
と書ける. 流体力学的相互作用とは, コロイド粒子の運動により周囲の流体が受けた力が
時間とともに伝搬し, 他の粒子に力を与えるというものである. 今日では, 流体力学的相
互作用はコロイド分散系におけるガラス転移点を決定する要因の一つとして非常に重要
であることが知られている.
流体力学的相互作用がコロイド粒子のダイナミクスに与える影響は大きく分けて二つ
ある. 一つ目が, 短時間自己拡散係数 DSS の生成, 二つ目が長時間自己拡散係数 DSL の低下
である. これら二つの影響により, 実際のコロイド分散系では流体力学的相互作用を考慮
しない計算機実験に比較してガラス転移体積分率が低くなることが知られている.
2.2.2
レナード・ジョンズ相互作用
希ガスの分子間力のモデルとして, 下記のレナード・ジョンズ相互作用がよく知られて
いる.
Φ(r) = 4²
·³ ´
σ 12
r
8
−
³ σ ´6 ¸
r
.
(2.9)
第 2 章. モデル
2.2 相互作用
図 2.1: Kob-Andersen 型の相互作用レナード・ジョンズ相互作用.
σ, ² はそれぞれ長さ, エネルギーの次元を持ったパラメータであり, アルゴンの場合では
²/kB = 120K, σ = 3.4 × 10−10 m と見積もられている. 第一項は粒子同士の斥力を表して
いる. べき指数の 12 という値には意味は無いが, 第二項の倍数であるため計算機内での扱
いは容易になっている. 第二項はファンデルワールス引力と呼ばれ, 斥力よりも長い範囲
で作用する. Kob と Andersen は2成分金属合金の一つである Ni80 P20 をレナード・ジョ
ンズ相互作用によりモデル化し, 計算機実験を遂行した [31]. 距離 r だけ離れた位置にあ
る粒子 α と粒子 β(α, β ∈ (A, B)) の間の相互作用は,
·³
¸
σαβ ´12 ³ σαβ ´6
αβ
Φ (r) = 4²αβ
−
.
r
r
(2.10)
で記述される. A 粒子が Ni, B 粒子が P に対応しており, 粒子数の比は NA : NB = 4 : 1,
数密度は約 1.2 である.
彼らが選んだ計算条件は安定した過冷却状態を再現することが知られており, 本研究で
もこの相互作用, 粒子数比, 数密度を用いることにする. Kob-Andersen モデルのパラメー
タを以下に示す.
σAA = 1.0, ²AA = 1.0,
σAB = 0.8, ²AB = 1.5,
σBB = 0.88, ²BB = 0.5.
(2.11)
このパラメメータ値でポテンシャル Φαβ (r) の概形を描いたのが図 2.1 である. 図 2.1 より,
9
第 2 章. モデル
2.3 システムの導入
異粒子間の相互作用が最も強く, 安定距離も短いことが分かる. これにより系は低温にお
いても結晶化や相分離を起こすこと無く, 過冷却状態へと移行する.
2.3
システムの導入
N 個の球形コロイド粒子が, 溶媒中に分散している系を考える. 粒子は質量 mA の A 粒
子と質量 mB の B 粒子の二種類がそれぞれ NA , NB 個存在しているとする. 溶媒の温度
は T , 粘性率は η で一定であるとする. i(= 1 . . . N ) 番目の粒子の運動は, 次のランジュバ
ン方程式で記述出来る. いま tαB = mα /γ α 程度の時間スケールに着目すると,
d α
X (t) = Viα (t),
dt i
(2.12)
X X αβ
d α
m
Vi (t) = −γ α Viα (t) +
Fij (t) + Riα (t).
dt
j
β
Nβ A,B
α
(2.13)
ここで Xiα (t) は時刻 t の i 番目の粒子 α の位置ベクトル, Viα (t) は時刻 t の i 番目の粒子 α
の速度ベクトル, γ α は粒子 α の速度に比例した抗力係数, Fijαβ (t) は時刻 t に i 番目の粒子
が j 番目の粒子 β から受ける力, Riα (t) は時刻 t に i 番目の粒子が周囲の溶媒から受ける
力であり, 以下の揺動散逸関係を満たす.
D
E
β 0
α
Ri (t)Rj (t ) = 2γ α kB T δ(t − t0 )δij δαβ .
(2.14)
また, 粒子間の相互作用は Kob-Andersen 型のレナード・ジョンズ相互作用を用いる.
Fijαβ (t) = −∇αi Φαβ (|Xiα (t) − Xjβ (t)|).
(2.15)
また, 相互作用に合わせて A 粒子と B 粒子の比, および数密度も Kob らの系と同じ条件
を設定した.
本研究で扱う系では粒子の体積分率は 0.59 と大きく, 実際のコロイド分散系では流体力
学的相互作用が系のダイナミクスに大きく影響を与える. しかし, 複雑な長距離力である
流体力学的相互作用を効率的に計算する方法は未だ開発されておらず, 我々が興味を持っ
ている系の構造が変化するような時間領域まで計算機実験を行うことは非常に困難であ
る. 従って, 流体力学的相互作用については計算を通じて gij = 0 とし, 働かないものとし
た. しかし, コロイド過冷却液体状態の本質的な振る舞いを調べるには十分な系であると
考えられる.
図 2.2 に今回計算を行なう2成分系のスナップショットを示す. 色の薄い粒子が A 粒子,
濃い粒子が B 粒子に対応している.
10
第 2 章. モデル
2.4 時空間スケール
図 2.2: 計算モデルのスナップショット. T = 0.50 において t = 10000 まで計算した粒子
配置を用いて作成した.
2.4
時空間スケール
ランジュバン方程式で記述されるブラウン運動の特性時間 tB の間に粒子が移動する距
離は lB という長さで特徴づけられる. lB はブラウン粒子の半径の千分の一程度の長さで
ある.本研究ではブラウン粒子同士の相互作用による運動, それによる系の構造が変化と
いったブラウン粒子を単位とした現象を取り扱う.このような現象の時間スケール, 空間ス
ケールはtB や lB よりも大きいものである. 空間スケールとしてはブラウン粒子の大きさ
を基準とするのが妥当であると考えられる. そこで, ブラウン粒子がその半径 a 程度移動
する時間,
tD =
a2
aγ
=
D0
kB T
(2.16)
を考える.tD は拡散時間と呼ばれる. tD は tB と比較すると 106 程度大きい値になる. tD
程度の時間スケールでランジュバン方程式 (2.2) 式の各項のオーダを考えてみると, 左辺
の慣性項は右辺の二項に比べて小さくなる. このとき, 2.13 式から慣性項を消去すると,
d α
1
Xi = α Fiα (t) + fiα (t).
dt
γ
11
(2.17)
第 2 章. モデル
2.4 時空間スケール
を得る. ここで新たな揺動力の揺動散逸関係は,
D
fiα (t)fjβ (t0 )
E
kB T
δ(t − t0 )δij δαβ
γα
= 2D0α δ(t − t0 )δij δαβ
= 2
となる. D0α = kB T /γ は粒子 α の初期拡散係数である.
本研究における計算機実験では (3.7) を用い, 系の時間発展を計算する.
12
(2.18)
第 3 章 計算手法
3.1
計算式
2.4 節で述べた通り, 本研究の計算機実験では以下の方程式を計算に用いる.
d α
1
Xi = α Fiα (t) + fiα (t).
dt
γ
(3.1)
ここで fiα (t) は, 以下の揺動散逸関係を満たす.
< fiα (t)fjβ (t0 ) >= 2D0α δ(t − t0 )δij δαβ 1.
(3.2)
D0α = kB T /γ α は粒子 α の初期拡散係数である. γ α は一般的には粒子の半径に依存し, 粒
子種によって異なるのだが, 本研究では計算時に γ A = γ B = γ = 1 とおく. これにより,
A 粒子と B 粒子の初期拡散係数は同じ値になる.
さて, 粒子間の相互作用 Fijαβ について, 式 (2.10) を用いて実際に計算してみると,
Fijαβ = −∇αi Φ(Xijαβ )
Ã
! αβ
6
6
24²αβ σαβ
σαβ
Xij
=
−1
αβ 7
αβ 6
(Xij )
(Xij )
Xijαβ
(3.3)
となる. 距離 r にある二粒子 α と β の間に働く力の大きさ F αβ (r) を図 3.1 に示す. 図 3.1
より r ≥ 2.5 に存在する粒子同士にはほとんど力が働かないことが分かる. よって計算に
際しては力のカットオフ距離 Rf を設定し, Rf の距離内に存在する粒子同士の相互作用の
みを考えることにする.
3.2
無次元化
ここでは計算中の単位量について考える.
まず長さであるが, これは σAA を単位量とする. これは主成分である A 粒子の直径に対
応した量であり, 粒子同士の相互作用による運動や, 系の構造の変化といった粒子を単位
とした現象を記述するには適当な量であると思われる.
13
第 3 章. 計算手法
3.2 無次元化
図 3.1: Kob-Andersen 型の相互作用レナード・ジョンズ相互作用で記述される粒子間に
作用する力. 値が正の距離では斥力, 負の距離では引力が作用する.
2
A
A
次に時間であるが, これは tA
D = σAA /D0 = γ σAA /kB T を用いる. tD は粒子 A が自由
拡散によって σAA 程度移動する量である. ここで, tA
D は, 半径を用いて定義した通常の拡
散時間に比べ 4 倍大きくなる点に注意が必要である. 特に他の計算とデータを比較する様
な場合には気をつけねばならない.
エネルギーは kB T , 力は kB T /σAA を用いて無次元化する.
次に, 3.1 節で確認した式について無次元化した変数を導入する. 無次元量をハット (. ˆ. .)
を用いて表すものとする. まず 3.7 式は,
dX̂iα
d(Xiα /σAA )
=
d(t/tA
dt̂
D)
α
tA
D dXi
=
σAA dt
½
¾
γ A σAA 1 α
α
F + fi
=
kB T
γA i
Fiα
fiα
1
+
= α A
γ /γ kB T /σAA kB T /γ A σAA
1
= α F̂iα + fˆiα ,
γ̂
14
(3.4)
第 3 章. 計算手法
3.3 差分式の導出
となる. 揺動散逸関係 3.8 式は,
µ
<
fˆiα (t)fˆjβ (t0 )
> =
µ
γ A σAA
kB T
¶2
< fiα (t)fjβ (t0 ) >
¶2
γ A σAA
2kB T
=
·
δ(t − t0 )δij δαβ 1
kB T
γα
γ A δ(t − t0 )
δij δαβ 1
= 2 α
γ
tA
D
1
= 2 α δ(t̂ − tˆ0 )δij δαβ 1,
γ̂
(3.5)
となる. この無次元化により, 揺動力から温度の情報が消去されている点に注意しよう.
同様に粒子間相互作用 Fijαβ についても計算すると,
F̂ijαβ
Fijαβ
=
kB T /σAA
"
#
6
6
2σαβ
σAA 24²αβ σαβ
=
− 1 (riα − rjβ )
αβ
αβ
6
kB T
rij
(rij )
"
#
12
24²̂αβ σ̂αβ
2
1
1
(r̂iα − r̂jβ ),
−
=
αβ
αβ 6
6
(σ̂αβ )
r̂ij (r̂ij )
T̂
(3.6)
となる. 揺動力とは対照的に, 相互作用には無次元化によって温度の情報が現れている.
無次元化の前は温度を変化させても相互作用の大きさは変わらず, 揺動力, すなわち粒子
の動き易さが変化する事により系の状態が変化するという描像であった. これに対し, 無
次元化後は揺動力は常に一定で, 相互作用の強さが変化する事により系の状態が変化する
という描像になっている. 適当に単位量を選べば前者の描像のまま無次元化を行うこと
も可能である. 今回はデータ整理の都合上後者を採用したが, 両者は全く等価であり, 無
次元化の如何で計算結果は変化しない.
3.3
差分式の導出
ここでは計算機で方程式を数値的に解くために必要となる差分式を導出する. 無次元
化の結果と, γ α = 1 という条件から計算式は,
dX̂iα (t)
= F̂iα (t) + fˆiα (t),
dt̂ E
D
β
α
fˆi (t)fˆj (t0 ) = 2δ̂(t − t0 )δij δαβ 1
15
(3.7)
(3.8)
第 3 章. 計算手法
3.3 差分式の導出
となる. (3.7) の両辺を t ∼ t + ∆t で積分する. ここで δt は計算機実験を行なう際の時間
刻みに対応しており, 時間の単位量 tA
D に比べて十分小さい値である. 左辺と右辺第一項
については,
Z
t+∆t
ds
t
Z
dX̂iα (s)
= X̂iα (t + ∆t) − X̂iα (t),
ds
(3.9)
t+∆t
dsF̂iα (s) = F̂iα (t)∆t + O((∆t)2 )
(3.10)
t
となる. 右辺第二項について, 新たな確率変数 Uiα (t, t + ∆t) を導入して,
Z
Uiα (t, t
t+∆t
dsfˆiα (s)
+ ∆t) =
(3.11)
t
とおく. Uiα (t, t + ∆t) について, 揺動散逸関係 3.8 式から,
*½Z
¾ (Z
D
E
t+∆t
β
dsfˆiα (s)
Uiα (t, t + ∆t)Uj (t0 , t0 + ∆t) =
Z
t
Z
t+∆t
=
t0 +∆t
ds
Z
t
t+∆t
=
Z
t0
Z
t0
t0 +∆t
ds
Z
t
t+∆t
=
t0 +∆t
ds
t0
0
t
t0 +∆t
t0
)+
ds0 fˆjβ (s0 )
D
E
ds0 fˆiα (s)fˆjβ (s0 )
ds0 2δ̂(s − s0 )δij δαβ 1
ds0 2δ̂(s − s0 )δ(t − t0 )δij δαβ 1
= 2∆tδ(t − t )δij δαβ 1
(3.12)
と計算できる. 従って分散 1 の Gauss 乱数 u を用いて,
Uiα (t, t + ∆t) =
√
2∆t(uex + uey + uez )
(3.13)
と書ける. 改めて u = uex + uey + uez とおけば求める差分式は,
X̂iα (t + ∆t) = X̂iα (0) + F̂iα (t)∆t +
√
2∆tu + O((∆t)2 )
(3.14)
と求められる.
実際に計算する際には Gauss 乱数を発生させるのは計算負荷が大きいため, これを一様
√ √
乱数で置き換える. その際に乱数の分散が一致するように, 平均 0, 大きさ (− 3, 3) の
一様乱数を用いる (付録 A.3 参照).
16
第 3 章. 計算手法
3.4
3.4 周期境界条件
周期境界条件
我々が興味を持つ系は液体, 過冷却液体状態, およびガラス状態ののバルクにおける性
質であるが, これを計算機実験で調査する際にはいくつかの困難が存在する. 例として,
一辺 1cm の立方体中に半径 1µm の球形コロイド粒子が体積分率 0.59 で分散している場
合を考えると, 粒子数 N は
N=
0.59 × (10−2 )3
= 1.4 × 1011 .
4
−6 )3
π
×
(10
3
(3.15)
となり, 膨大な数の粒子が含まれていることが分かる. 現実的に計算機実験が出来るのは
せいぜい数千から数万程度の粒子数までであり, この系の全ての粒子の運動を計算機実験
で追跡することは不可能である. また粒子数が数万の計算機実験において単位セル内に
粒子を配置すると, x, y, z それぞれの方向には高々数十の粒子しか存在しないことにな
る. これでは粒子の運動は単位セルの壁面の運動を大きく受けることになり, バルクの性
質を再現出来ないということが懸念される.
これらの問題への対策として, 周期境界条件を用いるという方法がある. 周期境界条件
の概念図を図 3.2 に示す. 実際の計算は三次元空間内で行なわれるが, 図では分かり易い
ように二次元で示した. 図中には九つの正方形が描いてあるが, このうち中心の正方形が
基本セルと呼ばれるものであり, 実際に粒子が存在する空間である. 色のついた球が粒子
に対応する. 周期境界条件では境界上に壁面は存在しないが, 代わりに向かい合う境界が
同一視される. 例えば, 基本セル内のある粒子が左側の境界からセルの外に飛び出すと,
右側の境界の同じ位置から同じ速度を持った同種の粒子が入ってくるという寸法である.
これにより壁で囲まれていないにも関わらず, 系の粒子数, 運動量は保存される. 粒子同
士の相互作用を考える際にも同様であり, 粒子は基本セル内の粒子だけでなく, 周囲の仮
想セル (必要であればさらに外側にある仮想セルを含む) から力を受ける. これにより計
算する粒子の少なさによる壁面からの影響を取り除き, よりバルクに近い状態で計算機実
験が可能となる.
3.5
初期配置
ブラウン動力学法における計算機実験では初期条件としてそれぞれの粒子の位置を与
える必要がある. 今回は粒子の初期配置として, 主にランダム配置を用いた. この他に結
17
第 3 章. 計算手法
3.5 初期配置
図 3.2: 周期的境界条件の概念図.
晶格子上に粒子を配置する方法もある. 一般的には系が平衡状態にあれば最終状態は初
期状態に依存しないため, どちらを用いても同じ結果が得られるはずである. しかし, 相
転移点付近では系の緩和時間が非常に長くなり, 実行可能な計算時間内では系の状態に初
期配置依存性が現れることがある. 固液転移においては双方の初期配置とも結晶化する
最高の温度を凝固点, 融解する最低の温度を融点とし, この初期配置依存性が現れる領域
を共存領域と呼ぶ. 近年, この共存領域はシミュレーションで考える系の大きさに依存し,
剛体球系では粒子数を多くするほど共存領域が縮小し, 融点と凝固点が接近することが報
告されている [40].
粒子数 10976, 体積分率 0.59 の単成分系において作成した面心立方格子型, ランダム配
置による初期配置のスナップショットを図 3.3, 図 3.4 に示す. ランダム配置の作成方法に
は Jodrey らによって開発された方法を用いた [41]. 面心立方格子型配置では格子点上に
粒子を配置したため, 周期的に同じ配置が繰り返されていることが分かる. これに対しラ
ンダム配置では粒子の配置に秩序が見られない. 次に, 図 3.3 および図 3.4 と同様の初期
配置について動径分布関数 g(r) を計算した結果を図 3.5 と図 3.6 に示す. g(r) はある粒子
を中心とした半径 r, 厚さ dr の球殻における粒子の数密度を, 系全体の数密度で規格化し
18
第 3 章. 計算手法
3.6 計算プログラムの高速化
図 3.3: 単成分系における面心立方格子
型配置のスナップショット.
図 3.4: 単成分系におけるランダム配置
のスナップショット.
図 3.5: 面心立方格子型配置における動
径分布関数.
図 3.6: ランダム配置における動径分布
関数.
たものである. 面心立方格子型配置では粒子が相対的に存在する距離が限定されており,
動径分布関数のピーク値が大きいものとなっている. また, 直径の 10 倍程度の非常に長
い距離を経ても相関が存在することが分かる. これに対し, ランダム配置では 7σAA 程度
で g(r) の値はほぼ 1 に収束しており, 長距離相関は存在しない. なお, 付録 A.4 で, 初期
配置の作成方法を簡単に述べた.
3.6
3.6.1
計算プログラムの高速化
粒子登録法
ブラウン動力学法の計算の主な流れは,
19
第 3 章. 計算手法
3.6 計算プログラムの高速化
1. 粒子位置から各粒子ごとの相互作用を計算
2. 揺動力の計算
3. 粒子位置の差分を計算
4. 粒子位置の更新と周期境界条件の適用
5. 計算結果の出力
6. 繰り返し
というものである. 計算のオーダーを考えると, 2 から 5 の計算は O(N ) であるが, 1 は粒
子の組を考えて和を取る必要があるため, 計算負荷は O(N 2 ) になる. これは単純に考え
ると粒子数を 10 倍にした際計算に要する時間が 100 倍になるいうことであり, 本研究の
ように緩和時間が長い系での計算機実験では大きな問題になる.
本研究で用いるレナード・ジョンズ相互作用は短距離力であり, ある粒子に力を与える
のはカットオフ距離内にある粒子のみである. 力のカットオフ距離 Rf を 2.5 として体積
分率 0.59 の系においてこの距離内にある粒子の数 Ncutoff を見積もってみると,
Ncutoff
0.59 × 34 π × 2.53
=
= 73.75.
4
π × 0.53
3
(3.16)
となりせいぜい 80 個程度であることが分かる. このため, 力の計算の際に全ての粒子の
組について計算を行なうのは大きな無駄である. 粒子登録法は粒子の位置が時空間的に
ある程度局在化していることを利用して, この不要な計算を少なくする手法である. その
概念図を図 3.7 に示す. 粒子登録法では計算開始時に全ての粒子の組について Rf よりも
長い距離 Rc 内にあるものを抽出し, リストに登録する. リストは十分な大きさを持った
二次元配列等を用いればよい. そして登録後 tbkm の間はリストにある粒子の組について
のみ Rf 内にあるかどうかの判定を行なう. ある時間内に粒子の移動する距離は有限であ
るから, Rc および tbkm の値を適当に選べば計算結果を粒子登録法を用いない場合と完全
に一致させることが出来る. Rc の値を小さくするとリストに登録される粒子の組は少な
くなるが, tbkm を小さくする必要がある. 計算する系や条件によって最適な Rc と tbkm の
組み合わせは異なるので, ある程度経験的に決定する必要がある.
20
第 3 章. 計算手法
3.6 計算プログラムの高速化
図 3.7: 粒子登録法の概念図.
図 3.8: サブセル法の概念図.
3.6.2
サブセル法
基本セルをさらに小さなセル, サブセルに分割することによって計算を高速化すること
が出来る [42]. この概念図を図 3.8 に示す. 計算機内では, 粒子位置をサブセルの大きさで
割り整数型の変数に格納することで容易に値を離散化できる. 離散化された粒子位置を
二次元配列に格納すれば, 任意の距離内にある粒子の組を取り出す際の計算手続きが容易
になる. 例えば, 力のカットオフ距離がサブセルの大きさ三個分に等しければ, 力の計算
は注目する粒子が入っているサブセルの周囲 26 個のセル内にある粒子を考えれば十分で
21
第 3 章. 計算手法
3.6 計算プログラムの高速化
ある. 周囲のサブセル内に存在する粒子は配列の番号を適当に操作すれば得られるから,
粒子間距離を計算する必要が無くなる. これにより, カットオフ距離内に存在する粒子を
抽出する手続きの計算手間を格段に減らすことが出来る.
本研究では基本的に粒子登録法を用い, tbkm ごとの粒子登録の際にサブセル法を利用し
た. これにより粒子登録法, サブセル法を単体で用いた場合よりも計算時間を短縮するこ
とに成功した.
22
第 4 章 解析方法
4.1
準備
2 章で導入された系について, コロイド粒子の持つ全エネルギー, すなわちハミルトニ
アンがどうなるか考えよう. 全エネルギーは運動エネルギーとポテンシャルエネルギー
の和として与えられるから, ハミルトニアン H(r N , pN ) は次のように書ける.
A,B Nα
A,B Nα
X pα 2
X
X
X
i
H(r , p ) =
+
Φαβ (Xijαβ )
α
2m
α
α
i
i
N
N
(4.1)
©
ª
©
ª α
A
A
B
B
ここで r N = X1A , . . . , XN
, X1B , . . . , XNBB , pN = pA
1 , . . . , pNα , p1 , . . . , pNB , pi は i
α
番目の粒子 α の運動量, Xijαβ = |Xiα − Xjβ | である. ハミルトニアンの変数として r N , pN
の計 6N 個の変数を用いたが, これはこれらの変数の値の組が系の 1 つの状態に対応して
いると考えることが出来るからである. r N と pN の作る 6N 次元の空間は位相空間 (phase
space) と呼ばれる. 位相空間全体の積分は以下のように定義される.
Z
Z
Z
N
dΓ =
dr
dpN
=
A,B Nα Z
YY
Zα
=
Z
ここで積分
Z
dXiα
dpαi
Z
Z
i
Z
. . . dXNBB
...
Z
Z
Z
Z
A
B
A
dp1 . . . dpNA dp1 . . . dpB
NB .
dX1A
dXNAA
dX1B
(4.2)
Z
dXiα
および
dXiα は座標空間および運動量空間全体の積分を表す. 例え
ば三次元直交座標系 (Xiα = (Xiα x , Xiα y , Xiα z ), pαi = (pαi x , pαi y , pαi z )) では,
Z
Z ∞
Z ∞
Z ∞
α
α
α
dXi =
dXi x
dXi y
dXiα z ,
−∞
−∞
−∞
Z ∞
Z ∞
Z ∞
Z
α
α
α
dpαi z ,
dpi y
dpi x
dpi =
−∞
−∞
である.
23
−∞
(4.3)
第 4 章. 解析方法
4.2 密度の導入
考えている系は体積一定で熱浴に接していることから, 力学量 A(r N , pN ) の平衡平均
hAi をカノニカル分布を用いて以下のように定義する.
Z
Z
­
®
N
N
N
A(r , p ) = dr
dpN feq (r N , pN )A(r N , pN ).
ここで feq (r N , pN ) はカノニカル分布であり, 以下のように与えられる.
Z
Z
1
N
N
N
feq (r , p ) = 3N
dr
dpN exp[−βH(r N , pN )].
h N !ZN
ここで h は Planck 定数, β は逆温度 (β = 1/kB T ) である. ZN は規格化条件
Z
Z
N
dr
dpN feq (r N , pN ) = 1
を満足するための関数であり,
1
ZN = 3N
h N!
Z
Z
dr
£
¤
dpN exp −βH(r N , pN )
N
(4.4)
(4.5)
(4.6)
(4.7)
で与えられる.
4.2
4.2.1
密度の導入
密度と密度揺らぎ
粒子 α(∈ (A, B)) の時刻 t, 位置 r における密度を次のように導入する.
α
n (r, t) =
Nα
X
δ(r − Xiα (t)).
(4.8)
Nα
=: nαeq
V
(4.9)
i=1
nα (r, t) の平衡平均は,
hnα (r, t)i =
である.
これらを用いると粒子 α の密度揺らぎ δnα (r, t) は次のようにかける.
δnα (r, t) = nα (r, t) − nαeq .
(4.10)
全く同様に, 全粒子の密度, 平均密度, 密度揺らぎを導入しておく.
n(r, t) =
A,B
X
nα (r, t),
(4.11)
α
N
,
V
δn(r, t) = n(r, t) − neq .
neq =
24
(4.12)
(4.13)
第 4 章. 解析方法
4.2 密度の導入
本研究のシステムでは運動量が緩和した後の時間スケールを議論するので, この密度, お
よび揺らぎはシステムを記述するすべての情報を持っていると考えられる. 実際, 他の物
理量はこれらから導出することが出来る.
4.2.2
粗視化された密度揺らぎ
ここでは粗視化された密度揺らぎを導入する. これによりある特定の長さスケール以
上の揺らぎを取り出して議論する事が可能となる.
まず, 実空間で定義された物理量 A(r) のフーリエ変換 F [A(r)] を以下のように定義
する.
Z
F [A(r)] :=
dre−ik·r A(r) = Ak.
(4.14)
ここで i は虚数単位, k は波数ベクトルである. 逆変換は以下のように定義される,
F
−1
1
[Ak] :=
(2π)d
Z
dkeik·r Ak.
(4.15)
ここで d は空間次元であるり, 今回は d = 3 である.
密度のフーリエ変換を計算すると,
Z
dre−ik·r n(r, t)
nk :=
=
=
A,B N Z
X
X
α i=1
A,B N
XX
α
dre−ik·r δ(r − Xiα (t))
e−ik·Xi (t) .
α
(4.16)
i=1
ここで波数カットオフ kc を導入して, 波数空間における密度 nk を二つの部分に分ける.
すなわち,
nk = nk≤kc (k, t) + nk>kc (k, t).
(4.17)
着目する空間パターンが rc = 2π/kc 程度以上の長さスケールで記述できる場合, 第二項
は不要な情報である. 第一項は系の情報を粗視化した密度揺らぎと考えることが出来る.
ここで粒子の大きさを考慮した, 粗視化された密度揺らぎを以下のように定義する.
nc (r, t) :=
A,B Nα
X
X
α
i
Z
1
dΩαi
Z
dΩαi
25
1
(2π)3
Z
α (t)
dkeik·(r−Yi
k≤kc
.
(4.18)
第 4 章. 解析方法
4.2 密度の導入
ここで, Yiα (t) は原点から粒子 i の任意の表面上の点へ延ばしたベクトルである. 粒子 i の
中心から任意の表面上の点へ延ばしたベクトルを aαi = Yiα − Xiα , 粒子 i の中心周りに張
る立体角を Ωαi とおくと, 粒子の大きさ程度以下の長さスケールについて, 角度依存性を
以下のように消去して計算出来る.
nc (r, t) =
A,B Nα
X
X
α
Z
dΩαi
i
Z
Z
1
dΩαi
k≤kc
dk ik·(r−Yiα (t)
e
(2π)3
Z
A,B Nα
X 1 Z
X
dk ik·(r−Xiα (t))
α
e
=
dΩαi e−ik·ai (Ωi )
3
4π k≤kc (2π)
α i=1
A,B Nα Z
sin(kai )
1 XX
α
dkeik·(r−Xi (t)) 4π
=
4
32π α i=1 k≤kc
kai
Z 2π Z π
A,B Nα Z kc
1 XX
sin(kai )
α
dk
dφ
dθeik|r−Xi (t)| cos θ k 2 sin θ
=
3
8π α i=1 0
kai
0
0
Z
A,B Nα Z kc
sin(kai ) 2 1
1 XX
α
dk
=
k
d(cos θ)eik|r−Xi (t)| cos θ
2
2π α i=1 0
kai
−1
A,B Nα Z kc
sin(k|r − Xiα (t)|) sin(kai )
1 XX
dk
=
2π 2 α i=1 0
ai |r − Xiα (t)|
A,B Nα
1 XX
=
f (kc |r − Xiα (t)|).
2π 2 α i=1
(4.19)
ここで,
ci cos(ci ) sin x − x sin(ci ) cos x
ci x(x2 c2i )
ci = kc ai .
f (x) =
である. 計算過程において,
Z π Z
dθ
0
2π
dφ sin θe
−ikr cos θ
Z
= −2π
0
·
= 2π
1
dk sin(kc a) sin(kc b) =
0
a2
(4.21)
dpe−ikrp
−1
−ikrp
e
ikr
sin(kr)
= 4π
kr
および
Z kc
(4.20)
¸1
−1
1
{−a cos(kc a) sin(kc b) + b sin(kc a) cos(kc b)} ,
− b2
を用いた.
26
(4.22)
(4.23)
第 4 章. 解析方法
4.3
4.3 動径分布関数と静的構造因子
動径分布関数と静的構造因子
ここでは動径分布関数と静的構造因子を導入する. これらは共に 2 粒子の位置の同時刻
相関に対応し, 実空間, および波数空間における系の構造を記述することが出来る.
まず, 粒子 α に関して, 密度 - 密度の自己相関関数を次式で定義する.
Z
1
αα
dr 0 hnα (r 0 + r)nα (r 0)i .
G (r) :=
N
(4.24)
密度の定義を代入して計算すると,
*( N
)( N
)+
Z
α
α
X
X
1
αα
0
0
α
0
α
G (r) =
dr
δ(r + r − Xi )
δ(r − Xj )
Nα
i=1
j=1
+
*N N
Z
α
α X
X
1
δ(r 0 + r − Xiα )δ(r 0 − Xjα )
=
dr 0
Nα
i=1 j6=i
+
*N
Z
α
X
1
+
δ(r 0 + r − Xiα )δ(r 0 − Xiα )
dr 0
Nα
i=1
Nα
Nα X
­
®
1 X
δ(r − [Xiα − Xjα ]) + δ(r).
=
Nα i=1 j6=i
ここで,
nαeq g αα (r)
Nα
Nα X
­
®
1 X
δ(r − [Xiα − Xjα ])
:=
Nα i=1 j6=i
(4.25)
(4.26)
とおくと次式を得る.
Gαα (r) = nαeq g αα (r) + δ(r).
(4.27)
gαα (r) は無次元の関数であり, 相対位置がr の位置に二つの粒子 α が存在する確率を表す.
系が等方的であれば g αα (r) は r = |r| のみに依存すると考えられる. r のみに依存した
関数 g αα (r) は角度に関する情報を消去して得られる. すなわち,
Z
1
αα
dΩi g αα (r(r, Ωi )).
g (r) = Z
dΩi
(4.28)
この g αα (r) が粒子 α の動径分布関数 (radial distribution function) である. g αα (r)dr はあ
る粒子 α の中心から見たときに, 他の粒子 α の中心が半径 r ∼ r + dr の球殻の数密度を
系全体の数密度で規格化したものである. その意味から, 式 4.26 より次のように表現する
ことも可能である.
®
­
α
α
|)
−
X
δ(r
−
|X
1
j
i
g αα (r) =
4πr2
nαeq
27
(4.29)
第 4 章. 解析方法
4.3 動径分布関数と静的構造因子
同様に, 波数空間における密度 nαk の自己相関 S(k) を考える.
S αα (k) :=
=
=
=
=
=
=
=
1 ­ α α ®
n n
Nα k −k
)+
*( N
)( N
α
α
X
X
1
α
α
e−ik·Xi
eik·Xj
Nα
i=1
j=1
Nα D
Nα X
E
α
α
1 X
e−ik·[Xi −Xj ]
Nα i=1 j=1
Nα D
Nα D
Nα X
E
E 1 X
1 X
−ik·[Xiα −Xjα ]
−k·[Xiα −Xiα ]
e
e
Nα i=1
Nα i=1 j6=i
*N N
+
Z
Z
α
α X
X
1
α
α
1+
dr dr 0
eik·[Xi −Xj ] δ(r − Xiα )δ(r 0 − Xjα )
Nα
i=1 j6=i
*N N
+
Z
Z
α
α X
X
1
0
dr dr 0eik·[r−r ]
δ(r − Xiα )δ(r 0 − Xjα )
1+
Nα
i=1 j6=i
Z
Z
1
α
α
dr dr 0eik·[Xi −Xj ] g αα(2) (r, r 0)
1+
Nα
Z
α
(4.30)
1 + neq dre−ik·r (g αα (r) − 1),
なる関係があることが分かる. 逆変換を考えれば次の関係が得られる.
Z
1
α
αα
neq (g (r) − 1) =
dkek·r [S αα (k) − 1] .
(2π)3
(4.31)
さて, 系の状態が等方的であれば, g(r) と同様に S(k) は k = |k| のみの関数となる.
g(r)dr = g(r)r2 sin θdrdθdφ であるから,
Z ∞ Z π Z 2π
eq
dr
dθ
dφr2 e−ikr cos θ (g(r) − 1)
S(k) = 1 + nα
0
0
Z ∞ 0
sin(kr)
= 1 + 4πneq
drr2 (g(r) − 1)
,
α
kr
0
(4.32)
なる関係が得られる. この S(k) は静的構造因子 (static structure factor) と呼ばれ, 散乱
実験から測定可能である. 同様に S αα (k)dk = S αα (k)k 2 sin θdkdθdφ であるから,
Z
sin(kr)
1
αα
α
.
neq (g (r) − 1) = 2 dkk 2 (S αα (k) − 1)
2π
kr
(4.33)
が成り立つ.
(4.32) 式および (4.33) 式は S αα (k) と g αα (r) の関係を与える. 互いに変換が可能であり,
両者が等価な情報を持っていることが分かる.
28
第 4 章. 解析方法
4.4
4.4 有効圧力
有効圧力
原子・分子系における計算機実験ではビリアル定理を用いて圧力を計算し, 系の相を判
定する指標とすることが出来る. ここではコロイド分散系でも同様の物理量を考え, 圧力
と同様に計算を行なう.
まず, コロイド分散系についてビリアル定理を導出する. 導出過程は文献 [44] を参考に
した. i 番目の粒子 α ∈ (A, B) の運動方程式は
mα
d α
Vi = −γ α Vi + Fiα + Wiα + Riα .
dt
(4.34)
ここで Fiα は他の粒子との相互作用, Wiα は壁から受ける力. Riα は周囲の溶媒から受け
る力であり, 以下の揺動散逸関係を満足する.
Riα (t)Rjβ (t0 ) = 2γ α kB T δ(t − t0 )δij δαβ 1.
(4.35)
粒子 i の位置ベクトル Xiα を内積して時間平均を取ると,
m
­
®
d
hViα · Xiα it − mα Viα 2 t = −γ α hViα · Xiα it + h(Fiα + Wiα ) · Xiα it + hRiα · Xiα it .
dt
(4.36)
となる. Vi · Xi は有限の値であるから左辺第一項は十分長い時間平均を取れば,
¿
À
iT
d
1h
Vi · Xi = 0.
(Vi · Xi ) = lim
T →∞ T
dt
0
t
(4.37)
となる. 他の項, hViα · Xiα it , hRiα · Xiα it も今回は tA
D 程度の時間に注目するので無視する.
すべての粒子について和をとると, 左辺は系の平均運動エネルギーになるので,
+
* A,B N
α
XX
1 α2
V
hKit =
2 i
α i=1
+
* A,B N
α
1 XX
(Fiα + Wiα ) · Xiα .
= −
2
α i=1
(4.38)
t
この, 系の運動エネルギーの長時間平均が粒子の座標と力の平均に等しいという関係をビ
リアル定理と呼ぶ. ここでビリアルとは式 (4.38) の右辺の量である. 全粒子の運動エネル
ギーの平均は,
3
hKit = N kB T,
2
29
(4.39)
第 4 章. 解析方法
4.5 平均二乗変位
となるので,
* A,B N
α
XX
α
+
= −3N kB T −
Xiα · Wiα
* A,B N
α
XX
α
i
+
Xiα · Fiα
(4.40)
i
となる. 左辺は全粒子が容器の壁から受ける力のビリアルの和になっている. ここで, 系
の圧力は容器の壁が粒子から単位面積辺りに受ける力の平均として与えられる. よって,
コロイド粒子が作る圧力を P とするとガウスの積分定理を用いて左辺を次のように計算
できる.
* A,B N
α
XX
α
+
Xiα
·
Z
=−
Wiα
P n · rdS = −3P V.
(4.41)
S
i
以上により,
PV
1
=1+
N kB T
3N kB T
* A,B N
α
XX
α
+
Xiα · Fiα
.
(4.42)
i
ここで, 本研究では粒子間相互作用はレナード・ジョンズ型相互作用を用いているので,
Fi =
A,B Nα
X
X
α
= −
j6=i
A,B Nα
X
X
α
=
Fijαβ
j6=i
A,B Nα
X
X
α
∇αi Φ(Xijαβ )
j6=i
6
24²αβ σαβ
"
8
(Xijαβ )
#
6
2σαβ
(Xijαβ )6
− 1 Xijαβ .
(4.43)
より, 圧力は次のように計算出来る.
A,B Nα Nβ
6
X 24²αβ σαβ
N kB T
1 XX
P =
−
V
6V α,β i j6=i (X αβ )6
ij
4.5
"
6
2σαβ
(Xijαβ )6
#
−1
Xijαβ · Xijαβ
Xijαβ Xijαβ
.
(4.44)
平均二乗変位
系のダイナミクスを観測する際, 最も直感的な方法はそれぞれの粒子の運動がどうなっ
ているかを観察するというものである. 粒子の運動を追跡する方法としては観測開始時
から各時刻において粒子がどれだけ動いたか, すなわち初期位置からの変位 ∆Xiα (t) を観
測するというものが考えられる.
∆Xiα (t) := Xiα (t) − Xiα (0).
30
(4.45)
第 4 章. 解析方法
4.6 ノンガウシアンパラメタ
粒子の変位の時間発展, すなわち粒子の軌跡を見れば直感的に粒子がどのような運動をし
ているか知ることが出来る. しかし粒子ごとに運動の様子は異なるため, 一つの粒子の軌
跡を追いかけても普遍的な情報を効率的に取り出すことが出来るとは言いがたい. 粒子
の平均的な振る舞いを知るためには平均を取ればよいが, 単純に初期位置からの変位の平
均を取ると 0 になってしまう. そこで, 次の平均二乗変位 (mean square displacement) が
用いられる.
®
® ­
­
M2α (t) := |Xiα (t) − Xiα (0)|2 = ∆Xiα (t)2
(4.46)
式 (4.46) は粒子 α の平均二乗変位の定義である. 系が十分一様であれば, 以下の式でも計
算出来ると考えられる.
M2α (t)
Nα
1 X
'
∆Xiα (t)2 .
Nα i
(4.47)
本研究では A 粒子の平均二乗変位 M2A (t) について主に解析を行う. 計算機実験から平
均二乗変位を計測する際には式 (4.47) を用いる.
4.6
ノンガウシアンパラメタ
ここでは 4.5 で導入した平均二乗変位を一般化しよう. 実際, 変位の n 次モーメントを
考えれば平均 n 乗変位が定義出来る.
Mnα (t) := h|Xiα (t) − Xiα (0)|n i = h∆Xiα (t)n i
(n = 1, 2, . . . )
(4.48)
系が等方的であれば, ∆Xiα (t) の分布も等方的になる. この時, 明らかに
α
M2m−1
(t) = 0 (m = 1, 2, . . . )
(4.49)
である.
∆Xiα (t) の分布関数 P (∆Xiα ) の形が分かれば, Mnα は直ちに計算出来る. ここでは Mnα
のおおまかな性質を調べるために, 例として ∆Xiα (t) が等方的で分布関数がガウス分布,
すなわち
PGauss (∆Xiα , t)
¸
·
∆Xiα 2
4π
exp −
=
(4πD0 t)d/2
4D0 t
31
(4.50)
第 4 章. 解析方法
4.6 ノンガウシアンパラメタ
に従う場合を考えよう. 実際, 自由ブラウン運動においては変位の分布は式 (4.50) に従う.
α
M2n
(t) を計算すると,
α
[M2n
(t)]Gauss =
=
=
=
=
®
∆Xiα (t)2n Gauss
Z
d(∆Xiα )PGauss (∆Xiα , t)∆Xiα (t)2n
Z ∞
4π
α
α
2(n+1) −∆∆Xiα (t)2 /4D0 t
d(∆X
)∆X
(t)
e
i
i
(4πD0 t)d/2 0
√
4π
1
(2n + 1)!!(2D0 t)n+1 40 t (付録 A.2 参照.)
d/2
(4πD0 t) 2
(2n + 1)!!
(M2α (t))n .
(4.51)
dn
­
となる.
ここで, M2 と M2n を用いて定義された以下の量を考える [45].
αn (t) =
α
(2n + 1)!! M2n
(t)
− 1.
α
n
3
(M2 (t))n
(4.52)
(4.51) 式に注意すると, 式 (4.52) の右辺第一項は, ∆Xiα の分布がガウス分布の時は 1 であ
り, このとき αn (t) = 0 となることが分かる. 逆に粒子の変位がガウス分布からずれる程,
すなわち拡散運動から外れる程, αn (t) の値は大きくなると考えられる.
この αn (t) をノンガウシアンパラメタ (non-Gaussian parameter) と呼ぶ. 本研究では以
下の α2A (t) を主に計測する.
α2A (t) =
3 M4A (t)
− 1.
5 (M2A (t))2
(4.53)
液体状態では α2A はピークを持つ形の関数となる. 短時間領域では粒子の運動は短時間自
己拡散運動であり, 変位はガウス分布に従うから, α2A (t) の値は 0 である. また, 長時間領
域では粒子同士の相互作用の結果, 粒子の運動過程は長時間自己拡散に移行する. この場
合も変位はガウス分布に従うため, 長時間領域でも α2A (t) の値は 0 となる. α2A (t) が値を
持つ領域は, 粒子の初期拡散運動が周囲の粒子のよって阻害され, 単純な拡散運動では記
述出来なくなっている領域である. 粒子が受ける力は他の粒子の配置に大きく依存する
から, α2A (t) は系の構造の情報を反映すると考えられる.
32
第 4 章. 解析方法
4.7
4.7.1
4.7 散乱関数
散乱関数
van Hove 関数
ここでは, 再び粒子密度の相関を考える. 4.3 節では同時刻の相関を議論したが, ここで
考えるのは時間が異なる密度の相関である. 以下の関数を考えよう.
Gαα (r, r 0; t) =
1
hnα (r 0 + r, t)nα (r 0 , 0)i .
Nα
(4.54)
これは, 初期時刻における r の密度と, 時刻 t において r 0 だけ離れた点における密度の相
関である. この相関から座標系の原点の取り方の影響を排除するために r 0 に関して空間
積分を行うと, van Hove により導入された以下の関数が得られる [46].
+
*N N Z
α
α X
X
1
dr 0 δ(r 0 + r − Xiα (t))δ(r 0 − Xjα (0)) ,
Gαα (r, t) =
Nα i=1 j=1
*N N
+
α
α X
1 X
=
δ(r − [Xiα (t) − Xjα (0)]) .
Nα i=1 j=1
(4.55)
この関数は van Hove 関数 (van Hove function) と呼ばれる. van Hove 関数ににおいて
t = 0 とおくと (4.27) 式と一致するのは明らかである. 同時刻相関の場合と同様に自己相
関と多相関の部分に分けることが出来る.
αα
Gαα (r, t) = Gαα
s (r, t) + Gd (r, t).
ここで,
Gαα
s
Gαα
d
*N
+
α
X
1
=
δ(r − [Xiα (t) − Xiα (0)]) ,
Nα i=1
*N N
+
α
α X
1 X
α
α
δ(r − [Xi (t) − Xj (0)])
=
Nα i=1 j6=i
(4.56)
(4.57)
である.
4.7.2
自己中間散乱関数
密度-密度相関を議論する際, 実空間で議論するよりも波数空間で議論した方が有用で
ある場合が多い. それらは散乱実験の結果と対応するからである. ここでは, 中間散乱関
数 (intermediate scattering function) と呼ばれる以下の関数を導入する.
F α (k, t) :=
®
1 ­ α
nk (t)nα−k(0) .
Nα
33
(4.58)
第 4 章. 解析方法
4.7 散乱関数
密度の定義を代入して計算すると,
)( N
)+
*( N
α
α
X
X
1
α
α
e−ik·Xj (0)
e−ik·Xi (t)
F α (k, t) =
Nα
j=1
i=1
*N N
+
α
α X
¡
£ α
¤¢
1 X
exp −ik · Xi (t) − Xjα (t)
=
Nα i=1 j=1
+
*N N Z
α
α X
X
1
dre−ik·r δ(r − [Xiα (t) − Xjα (0)])
=
Nα i=1 j=1
Z
=
dre−ik·r Gαα (r, t)
= F [Gαα (r, t)] .
(4.59)
となり, van Hove 関数のフーリエ変換であることが分かる. 中間散乱関数は波数空間にお
ける密度の初期値との相関であり, 時間発展を調べることにより初期値からの緩和過程を
見ることが出来る.
さて, 中間散乱関数も自己相関と他相関の部分に分けることが出来る.
F α (k, t) = Fsα (k, t) + Fdα (k, t).
(4.60)
ここで,
1
Fsα (k, t) :=
Nα
1
Fdα (k, t) :=
Nα
*N
α
X
+
−ik·[Xiα (t)−Xiα (0)]
e
i=1
*N N
α
α X
X
.
+
−ik·[Xiα (t)−Xjα (0)]
e
.
(4.61)
i=1 j6=i
である.
Fsα (k, t) を特に自己中間散乱関数 (self-intermediate scattering function) と呼ぶ. 系が
等方的であれば Fsα は |k| = k のみの関数になり,
*N Z
+
Z π
α
2π
X
1
α
−ik∆Xiα (t) cos θ
dφ
Fs (k, t) =
dθ sin θe
4πNα i=1 0
0
+
*N
α
sin(k∆Xiα (t))
1 X
, (4.62)
=
Nα i=1 k∆Xiα (t)
と書ける. ここで k が小さいと仮定して sin x の Taylor 展開,
sin x = x −
1 3 1 5
x + x − ...
3!
5!
34
(4.63)
第 4 章. 解析方法
4.7 散乱関数
より k 4 の項までを計算すると,
Fsα (k, t)
¸
Nα ·
1 α
1 X
1
2
α
4
6
1 − M2 (t)k +
=
M (t)k + O(k ) .
Nα i=1
6
120 4
(4.64)
となることが分かる.
さらに, ln(1 + x) の Taylor 展開,
1
1
1
ln(1 + x) = x − x2 + x3 − x4 + . . .
2
3
4
を利用して (4.64) 式を展開し直すと,
µ
¶
1 α
1
α
2
α
4
6
ln(Fs (k, t)) = ln 1 − M2 (t)k +
M (t)k + O(k )
6
120 4
¶
µ
1
1 α
α
2
4
4
M (t)k + O(k )
=
− M2 (t)k +
6
120 4
µ
¶2
1
1 α
1
α
2
4
6
−
− M2 (t)k +
M (t)k + O(k ) + . . .
2
6
120 4
1
1
1
= − M2α (t)k 2 +
M4α (t) −
(M2α (t))2 + O(k 6 )
6
120
72
1 α
1
2 4
2
α
= − M2 (t)k +
(M2 (t)) k · α2α (t) + O(k 6 ).
6
72
より,
Fsα (k, t)
¸
·
1
1 α
2 4
2
α
α
4
(M2 (t)) k · α2 (t) + O(k ).
= exp − M2 (t)k +
6
72
(4.65)
(4.66)
(4.67)
となり, Fsα (k, t) とノンガウシアンパラメタ α2α (t) が関係づけられる. O(k 6 ) を無視して
考えると, α2α (t) が十分小さいとき, 例えば粒子が拡散運動に近い運動をしている時には
Fsα (k, t) は M2α (t) の増加とともに指数関数的に減少することが分かる. 実際, 温度が高く
(あるいは体積分率が低く)α2α (t) の値が大きく増加しない様な場合には Fs (k, t) は指数関
数的に減少する. また温度が低い領域でも, 運動が短時間自己拡散係数 DSS に支配される
短時間領域, 長時間自己拡散係数 DSL に支配される長時間領域では α2α (t) の値は 0となり
,Fs (k, t) も部分的に指数関数に近い形で緩和する. これに対し, 粒子の運動が周囲の粒子
によって制限される中間時間領域では α2α (t) の値が 1 のオーダまで増大する. これにより,
(4.67) 式の指数関数の中の第二項が第一項と拮抗し, Fsα (k, t) の値はこの時間領域ではあ
まり減少しなくなる.
35
第 4 章. 解析方法
4.8
4.8 平均場理論
平均場理論
1994 年に徳山と Oppenheim は自身は流体力学的相互作用を考慮したランジュバン方
程式 ((2.5) 式) を出発式として中性コロイド分散系における非線形確率拡散方程式を導出
し, 中性コロイド分散系の長時間自己拡散係数と体積分率 φ の関係を導いた [28, 39].
DS (φ) = DSS (φ)
DS (φ)
1+ S
D0
µ
1−
φ
φTc O
9
φ
32
¶µ
¶−2 .
φ
1 − TO
φc
(4.68)
ここで φTc O は拡散係数の値が仮想的に 0 になる特異点である. 流体力学的相互作用にっ
て生成される短時間自己拡散係数 DSS も体積分率の関数として以下のように与えられる
[28].
D0
,
1 + H(φ)
2b2
c
bc(2 + c)
H(φ) =
−
−
.
1 − b 1 + 2c (1 + c)(1 − b + c)
DSS (φ) =
(4.69)
(4.70)
ここで b = (9φ/8)1/2 , c = 11φ/16 である. 式 (4.70) に今回計算するモデルの体積分率を
代入すると, DSS ' 0.27 となる.
また, 式 (4.68) を一般化した次の式は様々な系において長時間拡散係数の振る舞いを記
述出来ることが知られている [30].
DS (φ) =
1+
DSS (p)
µ
D0
D0
¶µ
¶−2 .
p
p
1−
pc
pc
(4.71)
ここで p はコントロールパラメータであり, 中性コロイド分散系では体積分率, 原子・分
子系や本研究で扱う系では逆温度 1/kB T に対応する.
平均場理論はコロイド分散系のダイナミクスを記述する理論として, 以下の平均二乗変
位の時間発展方程式として与えられる [47].
¤
£
d
2
M2 (t) = 2dDSL + 2d DSS − DSL e−M2 (t)/lf .
dt
(4.72)
ここで lf は自由長 (free length) と呼ばれるフィッティングパラメータであり, ケージの大
きさに対応する. 式 (4.72) は, コロイド分散系において自己中間散乱関数が式 (4.67) の k 2
までの式, すなわち,
¸
· 2
k M2 (t)
Fs (k, t) = exp −
6
36
(4.73)
第 4 章. 解析方法
4.8 平均場理論
図 4.1: lf = 0.1, DSS = 1.0, DSL = 1.0 × 10−3 , a = 0.5 における平均場理論による平均二
乗変位とその特性時間の比較. 図中の点が左から順に tf , tγ , tβ , tL に対応している. 破線
は式 (4.76) の曲線.
に従うと仮定し, これと非線形確率拡散方程式を用いて導出された [28].
(4.72) 式は容易に解くことが出来て (付録 A.5 参照),
·
´¸
DSS ³ 2dDSL t/lf2
s
−1 ,
M2 (t) = lf ln 1 + L e
DS
(4.74)
となり, 平均二乗変位を解析的に表すことが出来る.
式 (4.74) の振る舞いについて考えてみる. まず, t が lf2 /DSL に比べて十分小さいとして,
L
2
e2dDS t/lf ' 1 +
2dDSL t
lf2
(4.75)
を式 (4.74) へ代入すると次式を得る.
#
S
2dD
t
S
M2 (t) ' lf2 ln 1 +
.
lf2
"
(4.76)
tβ := lf2 /DSL と定義すれば, t ¿ tβ で M2 (t) は対数関数的に増加することが分かる. tβ の
37
第 4 章. 解析方法
4.8 平均場理論
他にも, 平均場理論の平均二乗変位について以下の様な特性時間が定義出来る [29].
tf :=
lf2
DSS
tγ := (tf tβ )1/2
lf2
tβ :=
DSL
a2
tL :=
.
DSL
(4.77)
(4.78)
(4.79)
(4.80)
これらの特性時間と式 (4.74) 式の比較を図 4.1 に示す.
式 (4.74) は DSS により支配される短時間自己拡散運動から DSL に支配される長時間自己
拡散運動へのクロスオーバーを表現していることが分かる. 図 4.1 からそれぞれの特性時
間について, tf は粒子の運動が制限され始める時間, tγ は平均二乗変位の振る舞いが対数
的な増加から離れ始める時間, tβ は粒子の運動が長時間自己拡散へと移行する過程を特徴
づける時間, tL は粒子の運動が長時間自己拡散運動に移行する時間, とその意味をとらえ
ることが出来る.
38
第 5 章 計算結果と考察
5.1
5.1.1
単成分系
計算条件
単成分系における計算条件を以下に示す.
粒子数
差分法
時間刻み
システムサイズ
初期配置
境界条件
力のカットオフ距離
粒子登録半径
粒子登録ステップ数
サブセル分割数
:
:
:
:
:
:
:
:
:
:
10979
Euler 法
10−4
21.35
ランダム
周期境界条件
2.5
8.0
50
16
なお, システムサイズに関してはレナード・ジョンズ相互作用のパラメータである σαβ
(α, β ∈ (A, B)) を粒子の直径とみて, Kob らの計算機実験 [31] と体積分率が等しくなるよ
うに設定した.
5.1.2
有効圧力
単成分系における有効圧力の計算結果を図 5.1 に示す. 横軸は温度 T 縦軸が有効圧力で
ある. 温度が低くなるほど有効圧力の値が減少していることが分かる.
T = 5.00 から T = 1.50 までは有効圧力の値が連続的に変化しており, 一本の曲線上に
乗っている. これに対し, T = 1.50 と T = 1.20 の間の変化はそれより高い温度領域のも
のと比較して不連続であり, より低い温度である T = 0.50 でもその傾向は変わらない. こ
の原因として, T = 1.50 と T = 1.20 の間で何らかの相転移が起こっていることが推測さ
39
第 5 章. 計算結果と考察
5.1 単成分系
図 5.1: 有効圧力 P の温度依存性. 温度は T = 5.00, 4.50, 4.00, 3.50, 3.00, 2.50, 2.00, 1.50,
1.20, 1.00, 0.50 で計算した.
れる. 次に具体的に系の構造がどのように変化したのか調べるため, 次節で動径分布関数
を測定する.
5.1.3
動径分布関数
単成分系における動径分布関数 g AA (r) の計算結果を図 5.2 および図 5.3 に示す. 有効圧
力の結果において異なった振る舞いが見られた高温領域と低温領域がそれぞれ図 5.2 と図
5.3 に対応している.
二つの図の動径分布関数関数は明らかに異なった傾向を示している. まず, 高温領域
(T ≥ 1.50) では r ∼ 5.0 付近までの短距離相関のみが存在し, それより大きい距離では
g AA (r) の値はほぼ 1 である. 液体の特徴的な動径分布関数である. これに対し低温領域
(T ≤ 1.20) では r ≥ 10 の領域でも相関が存在している. 図中の点線は面心立方格子型配
置の動径分布関数である. さらに T = 1.20, 1.00, 0.50 についてピークの位置が面心立方
格子型配置とよく対応していることが見て取れる. このことから, T = 1.20, 1.00, 0.50 の
三つの温度では系が結晶化を起こしているということが分かる.
40
第 5 章. 計算結果と考察
図 5.2: 単成分系における動径分布関数
g AA (r). 下から順に T = 5.00, 4.50, 4.00,
3.50, 3.00, 2.50, 2.00, 1.50.
5.1 単成分系
図 5.3: 単成分系における動径分布関数
g AA (r). 下から順に T = 1.20, 1.00 0.50.
破線は面心立方格子の動径分布関数.
図 5.4: 平均二乗変位 M2A (t) の時間発展. 左から T = 5.00, 4.00, 3.00, 2.00, 1.00.
5.1.4
平均二乗変位
単成分系における平均二乗変位 M2A (t) の時間発展を図 5.4 に示す. いずれの温度でも時
間が小さい領域では M2A (t) の値が傾きが 1 に従って増加していることが見て取れる. こ
41
第 5 章. 計算結果と考察
5.1 単成分系
図 5.5: ノンガウシアンパラメータ α2A (t) の時間発展. 左から T = 5.00, 4.00, 3.00, 2.00,
1.00. 温度が低いほど細い実線で表示されている.
れは短時間自己拡散係数運動であり, 流体力学的相互作用を無視した本研究では,
M2A (t) ' 6D0 t.
(5.1)
である.
短時間領域の振る舞いは一致するが, それ以降の時間領域では T ≥ 2.00 と T = 1.00 の
間で M2A (t) の振る舞いが大きく異なる. T ≥ 2.00 では一時的に M2A (t) の傾きが減少する
ものの値は増加続け, 十分時間が経過すると再び傾き 1 の拡散運動, 長時間自己拡散を示
している. 対照的に, T = 1.00 では M2A (t) の値がある程度以上は増加せず, 粒子の運動が
制限されていることが分かる. これは結晶格子が安定で, 粒子を捕捉しているためと考え
られる.
5.1.5
ノンガウシアンパラメータ
単成分系におけるノンガウシアンパラメータ α2A (t) の時間発展を図 5.5 に示す. 高温領
域 (T ≥ 2.00) では α2A (t) は t=0.1 ∼ 1.0 程度の時間領域で 0.1 ∼ 0.2 程度のピークを示し
ている. ピークの大きさは温度が低くなるほど大きくなっているのが見て取れる. これに
対し, T = 1.00 では値が大きく, 安定しないものとなっている. 高温領域では十分時間が
経つと α2A (t) の値は 0 になっているが, T = 1.00 では少なくとも今回計算を行なった時間
42
第 5 章. 計算結果と考察
5.1 単成分系
図 5.6: 自己中間散乱関数 M2A (t) の時間発展. 左から T = 5.00, 4.00, 3.00, 2.00, 1.00.
範囲においては α2A (t) の値は数十程度のままであった. この原因としては, 平均四乗変位
の振る舞いがまだ落ち着いていないということが考えられる.
5.1.6
自己中間散乱関数
単成分系における自己中間散乱関数 FsA (k, t) の時間発展を図 5.6 に示す. 波数について
は, 単成分系における静的構造因子 S(k)AA を計算していないため, 便宜的に2成分系に
おける S AA (k) が最大値を取る波数, k = 7.25 で計算を行なった. 高温領域では t = 1.0 程
度までに FsA (k, t) はほぼ 0 に緩和しているが, T = 1.00 では t = 1000 においても FsA (k, t)
の値は 0.8 程度となっている. T = 1.00 では平均二乗変位と同様に粒子が結晶格子に束縛
されているため, 粒子は格子内の運動しか出来ず, 構造緩和が制限されているものと考え
られる.
43
第 5 章. 計算結果と考察
5.2
5.2.1
5.2 2成分系
2成分系
計算条件
2成分系における計算条件を以下に示す.
粒子数
差分法
時間刻み
システムサイズ
初期配置
境界条件
力のカットオフ距離
粒子登録半径
粒子登録ステップ数
サブセル分割数
5.2.2
:
:
:
:
:
:
:
:
:
:
10976(A 粒子 8781, B 粒子 2195)
Euler 法
10−4
20.89
ランダム
周期境界条件
2.5
8.0
50
16
有効圧力
図 5.7 に様々な温度における有効圧力の計算結果を示す. 今回は t = 10000 まで計算し
た粒子配置を用いた. 温度が低下するにつれ, 有効圧力の値が減少するのが見て取れる.
この点においては単成分系の計算結果と同じである. しかし, 今回計算した温度領域にお
いて単成分系で観察されたような圧力の飛びは観察されなかった. このことから2成分系
ではこの温度領域で結晶化や相分離等の一次相転移は起こっていないと考えられる. 圧
力は系の構造と密接に関係していることから, 温度が低くなるにつれて系の構造も連続的
に変化していると予測される. 次節では, 実際に系の構造を調べるため動径分布関数につ
いて議論する.
5.2.3
動径分布関数と静的構造因子
図 5.8 と図 5.9 はそれぞれ T = 5.00, 0.50 における動径分布関数 g(r), g AA (r), g AB (r),
g BB (r) を比較したものである. 尚, 計算には t = 10000 まで計算した位置データを用い,
長さの刻み ∆r = 0.02 として計算した. まず, T = 5.00 の動径分布関数について考察を行
う. いずれの動径分布関数も r ≥ 5.0 ではほぼ値が 1 となっており, 長距離相関が存在し
44
第 5 章. 計算結果と考察
5.2 2成分系
図 5.7: 2成分系における有効圧力 P の温度依存性. 温度は T = 5.00, 3.00, 2.00, 1.50,
1.00, 0.90, 0.80, 0.70, 0.60, 0.55, 0.50, 0.47, 0.45, 0.44.
図 5.8: T = 5.00 における動径分布関数, g(r), g AA (r), g AB (r), g BB (r) の比較.
ないことが見て取れる. A 粒子の動径分布関数 g AA (r) に注目すると, 単分散における高温
領域の結果と似た端的な液体の形になっていることが分かる. A 粒子と B 粒子の相関で
ある g AB (r) も g AA (r) と同じような形が見られるが, 第一ピークはより小さい r の点に出
45
第 5 章. 計算結果と考察
5.2 2成分系
図 5.9: T = 0.50 における動径分布関数, g(r), g AA (r), g AB (r), g BB (r) の比較.
現している. 全粒子の動径分布関数は第一ピークの内側に g AB (r) の第一ピークの影響が
表れている. B 粒子の動径分布関数 g BB も短距離相関を持ち長距離相関は存在しないと
いう点では同じであるが, 値が他の動径分布関数に比較して 1 に近くなっている. このこ
とから B 粒子はセル中にほぼ一様に分散していると考えられる.
次に T = 0.50 の動径分布関数について T = 5.00 と比較しながら考察を行う. 今回計算
を行った温度の中では低い温度であるが, 単成分系で見られたような長距離相関はいずれ
の動径分布関数にも見られない. これは有効圧力の結果から得られた予想と対応してい
る. しかし高い温度と同じという訳ではなく, 同じ2成分系の T = 5.00 における結果 (図
5.8) と比較すると, 違った特徴が見られる. まず, 全体的にピークの値が大きくなってい
る. これは粒子の運動が相互作用により束縛されるようになっていることの現れであり,
単成分系でも観察されている. 特に異粒子同士の相関 g AB (r) の第一ピークが大きくなっ
ており, それに伴って全粒子の動径分布関数 g(r) の第一ピークに割れが見られる. 相互
作用の構造に与える影響が大きくなったため, もっとも安定な相互作用を持つ A-B 粒子
が隣り合う確率が大きくなっているためと考えられる. B 粒子の動径分布関数では逆に
T = 5.00 で見られた第一ピークが消失し, 代わりに r = 1.38 付近にピークが出現してい
る. g AB (r) で見られたように, B 粒子は A 粒子と隣合う確率が高くなっているため B 粒
子同士が隣合う確率が小さくなっていることを表している. r = 1.38 に現れたピークは A
46
第 5 章. 計算結果と考察
5.2 2成分系
図 5.10: 様々な温度における動径分布関数 g(r), g AA (r), g AB (r), g BB (r) の計算結果. 下か
ら順に, T = 5.00, 3.00, 2.00, 1.50, 1.00, 0.90, 0.80, 0.70, 0.60, 0.50, 0.55, 0.50, 0.47, 0.45,
0.44 に対応している.
粒子を一つ挟んだ位置に存在する B 粒子との相関であると考えられる. また, g AA (r) は
T = 5.00 と同じような形をしているが, 第二ピークに注目すると分裂していることが分
かる. 同様の分裂が g AB (r) にも見られる. この第二ピークの分裂は過冷却液体の動径分
47
第 5 章. 計算結果と考察
5.2 2成分系
布関数に現れる特徴の一つとして知られているものである.
図 5.10 に様々な温度における動径分布関数 g(r) の計算結果を示す. 温度は T = 5.00,
3.00, 2.00, 1.50, 1.00, 0.90, 0.80, 0.70, 0.60, 0.55, 0.50, 0.47, 0.45, 0.44 である. 曲線が重
ならないように 1 ずつずらして表示してある. 温度が低くなるにつれて, 図 5.8 と図 5.9 の
間に見られるような違いが徐々に現れていることが分かる. 高い温度と低い温度で動径
分布関数の特徴は異なるが, 単成分系のように温度変化に対する急激な構造の変化は観察
されない. これは有効圧力に飛びが見られないという点と対応している. 次に, ピークの
分裂する温度について考察する. g(r) の第一ピークについては T = 3.00 で既に割れが見
られる. g AA (r), g AB (r) の第二ピークの分裂は T = 0.70 ∼ 0.60 付近から起こっているこ
とが分かる. この温度付近で g BB (r) の第一ピーク (高温における第二ピーク) も分裂して
いることが分かる.
図 5.11 に様々な温度における静的構造因子 S(k) の計算結果を示す. 温度は T = 5.00,
3.00, 2.00, 1.50, 1.00, 0.90, 0.80, 0.70, 0.60, 0.55, 0.50, 0.47, 0.45, 0.44 である. 曲線が重
ならないように 1 ずつずらして表示してある. 計算には式 (4.32) を用い, 図 5.10 に示した
各温度における動径分布関数をフーリエ変換することにより求めた. A 粒子同士の静的
構造因子 S AA (k) に注目すると, 最大値をとる波数は kmax = 7.25 であることが分かる. あ
とで計算する A 粒子の自己中間散乱関数 FsA (k, t) はこの波数を用いて計算する.
5.2.4
粗視化された密度揺らぎ
ここでは粗視化された密度揺らぎの計算結果を示す. カットオフ長は rc = 2.5 とし,
t = 5000 における粒子配置を用いて計算した. 図 5.12 に T = 1.00, 0.50 における nc = 0.2
の等値面を可視化したもの示す. 等値面の作成にはマーチン・テトラヘドラ法を用い, フ
リーの三次元レイトレーシングライブラリである Pov-Ray を用いて可視化を行った. 二
つの図に大きな違いは見られず, ある時刻のみを見ても温度による違いをとらえられない
ことが分かる.
今回, 密度揺らぎのスナップショットのみからは, 温度による違いをとらえることが出
来なかった. しかし, 空間を粗視化してとらえるという意味でこの手法は強力なものであ
り, 動的な解析を組み合わせれば過冷却やガラス状態の特徴を捉えるのに役立つと考えら
れる.
48
第 5 章. 計算結果と考察
5.2 2成分系
図 5.11: 様々な温度における静的構造因子 S(r), S AA (r), S AB (r), S BB (r) の計算結果. 下
から順に, T = 5.00, 3.00, 2.00, 1.50, 1.00, 0.90, 0.80, 0.70, 0.60, 0.50, 0.55, 0.50, 0.47,
0.45, 0.44 に対応している.
5.2.5
平均二乗変位
A 粒子および B 粒子の平均二乗変位 M2A (t) および M2B (t) の時間変化を図 5.13 および
図 5.14 にそれぞれ示す. t = 5000 から計測したものである.
49
第 5 章. 計算結果と考察
5.2 2成分系
T = 1.00
T = 0.50
図 5.12: T = 1.00, 0.50 における密度揺らぎの計算結果.
図 5.13: 様々な温度における平均二乗変
位 M2A (t) の時間発展. 左から温度 T =
5.00, 3.00, 2.00, 1.50, 1.00, 0.90, 0.80,
0.70, 0.60, 0.55, 0.50, 0.47, 0.45, 0.44.
図 5.14: 様々な温度における平均二乗変
位 M2B (t) の時間発展. 詳細は 5.13 に同
じ.
t ≤ 10−3 程度の短時間領域ではすべて温度で M2A (t), M2B (t), M2 (t) の時間変化が傾き
1 の直線にほぼ一致している. この直線は, M2A (t) = 6D0 t である. 本研究の計算では流体
力学的相互作用を考慮していないため, 短時間自己拡散係数 DSS は存在せず, 初期運動は
D0 = kB T /γ のみに支配され, 初期の運動の様子はすべての温度で同じである. t ≥ 10−3
程度から, M2A (t) の傾きが 1 よりも段々と小さくなってくる. この領域は, 粒子が周りの
50
第 5 章. 計算結果と考察
5.2 2成分系
粒子から力を受けることにより運動を阻害されている時間領域である. この時間領域では
丁度周囲の粒子が粒子をかごのようにとらえているようであり, ケージ効果と呼ばれてい
る. 温度が低くなるほどケージ効果が始まるのが早く, 長い時間領域にわたって粒子の運
動が制限されていることが分かる. このケージ効果の強さが長時間自己拡散係数の低下の
度合いと対応している. 単成分系の計算でも結晶状態でも粒子の運動が制限されるのが
観察されたが, 結晶状態では平均二乗変位の値が中間時間領域以降変化しないのに対し,
図 5.13 の M2A (t) は少しずつ値が増加している. これは平均場理論で予測される対数関数
に従ったものになっている [48]. 長時間領域では M2A (t) の変化が再び傾き 1 の直線になっ
ている. これは長時間自己拡散運動であり, 長時間自己拡散係数 DSL で特徴づけられる.
M2A (t) ∼ 6DSL t.
(5.2)
単分散系では液体状態と結晶状態で平均二乗変位が劇的に変化したが, 2成分系では温
度の変化に対する M2A (t) の変化は連続的である. しかし, 高温と低温の M2A (t) を比較す
るとその形が大きく異なっており, 系のダイナミクスが大きく変化していることがわかる.
次に図 5.14 について図 5.13 と比較しながら考察する. M2B (t) の基本的な振る舞いは
M2A (t) と同じである. 初期には D0 で拡散し, 次いで β 緩和が観察され, その後長時間自
己拡散へと移行する. より詳しく両者を比較するため, 図 5.15 に T = 1.00 および 0.50 に
おける M2A (t) と M2B (t) の比較を示す. 同じ温度で両者を比較すると B 粒子の方が A 粒
L(A)
子よりも大きく移動していることが分かる. A 粒子の長時間係数 DS
L(B)
間係数 DS
と B 粒子の長時
の比較を図 5.16 に示す. 温度が低いほど長時間自己拡散係数の差が大きく
なっていることが分かる.
5.2.6
ノンガウシアンパラメタ
まず, 図 5.17 および 5.18 に様々な温度における平均 4 乗変位 M4A (t) および M4B (t) の時
間発展を示す.
平均 4 乗変位は短時間領域, 長時間領域でそれぞれ,
(
120(D0 t)2 (短時間領域)
M4 (t) ∼
120(DSL t)2 (長時間領域)
(5.3)
に従う. A 粒子よりも B 粒子の長時間拡散係数が大きくなることを平均二乗変位の結果
から確認したが, M4A (t) と M4B (t) を比較するとその差がさらに顕著に現れていることが
見て取れる. 特に中間時間領域での振る舞いが大きく異なる.
51
第 5 章. 計算結果と考察
5.2 2成分系
L(A)
とB
図 5.15: T = 1.00, 0.50 における M2A (t)
と M2B (t) の比較.
図 5.16: A 粒子の長時間係数 DS
L(B)
粒子の長時間係数 DS の比較.
図 5.17: 様々な温度における平均 4 乗
変位 M2A (t) の時間発展. 詳細は 5.13 に
同じ.
図 5.18: 様々な温度における平均 4 乗
変位 M2B (t) の時間発展. 詳細は 5.13 に
同じ.
式 (4.52) から分かるように, 5(M2 (t))2 /3 と比較すると初期と後期で振る舞いが一致
する. この二つの比較を示したのが図 5.19 である. T = 1.00, 0.50 において M4A (t) と
5(M2A (t))2 /3 の時間発展を示した. 左の二つの曲線が T = 1.00, 残りが T = 0.50 であり,
5(M2A (t))2 /3 は細い曲線で示してある. 初期と後期で振る舞いが一致しているのに対し,
中間時間領域ではそれぞれの温度で M4A (t) の方が値が上回っている. そして, その差は粒
52
第 5 章. 計算結果と考察
5.2 2成分系
図 5.19: M4A (t) と M2A (t) の比較.
子の運動が後期拡散過程に移行するクロスオーバ点付近で最大になっており, T = 0.50 で
より顕著に値の差が現れている. M2A (t) のケージ効果が顕著に現れている時間領域から
DSL に支配される後期拡散過程に移行する付近で差が大きくなっていることが分かる. こ
の図における M4 (t) と 5(M2 (t))2 /3 の差 (両対数目盛りであるから, 正確には比) がノンガ
ウシアンパラメタに対応している. 従って, α2A (t) は運動過程が後期拡散過程に移行する
時間を特徴づける量であるととらえることが出来る.
図 5.20, 図 5.21 に様々な温度における A 粒子のノンガウシアンパラメタ α2A (t) および
B 粒子のノンガウシアンパラメタ α2B (t) の時間発展を示す. 図 5.19 で見た通り, 短時間領
域と長時間領域では値が 0 であり, ある時刻において最大値をとるという振る舞いが見ら
れる. 温度が低くなると最大値が大きくなることが分かる. これは図 5.19 から得られた知
見と対応する. また, α2A (t), α2B (t) の最大値が現れる時刻が温度が低くなるほど遅くなっ
ている. これは平均二乗変位において, 運動が長時間拡散に移行するのが温度が低くなる
ほど遅くなるという点と対応していると考えられる.
次に, 図 5.20 と図 5.21 を比較すると, 同じ温度で見た場合 A 粒子よりも B 粒子の α2 (t)
の方が最大値の値が大きく, 最大値が表れる時間が早くなっている. これらを具体的に比
較するために, α2A (t) と α2B (t) の最大値と最大値が現れる時刻の温度依存性をプロットし
たものを図 5.22, 図 5.23 に示す. 尚, 図の横軸は逆温度である. 図 5.22 より, α2 (t) の最大
53
第 5 章. 計算結果と考察
5.2 2成分系
図 5.20: 様々な温度におけるノンガウシ
アンパラメタ α2A (t) の時間発展. 詳細は
5.13 に同じ.
図 5.21: 様々な温度におけるノンガウシ
アンパラメタ α2B (t) の時間発展. 詳細は
5.13 に同じ.
図 5.22: ノンガウシアンパラメタ α2A (t),
α2B (t) の最大値の温度依存性. 実線が粒
子 A, 破線が粒子 B.
図 5.23: ノンガウシアンパラメタ α2A (t),
α2B (t) が最大値を取る時刻の温度依存性.
実線が粒子 A, 破線が粒子 B.
値は T ≥ 1.00 ではごく小さい値であるが温度が低くなるにつれ増大し, T = 0.60 付近か
ら急激に大きくなっていることが分かる. α2A (t) の最大値と α2B (t) の最大値を比較すると,
T ≥ 1.00 の高い温度ではほぼ同じ値であるが, 温度が低くなるにつれ α2B (t) の方が大きく
なっているのが分かる. T = 0.44 では (α2B (t))max /(α2B (t))max = 2.28 程度である.
54
第 5 章. 計算結果と考察
5.2 2成分系
図 5.24: 対数正規分布によるフィッティングで得られた tNGP と元の α2A (t) が最大値を取
る時刻の比較.
同様に図 5.23 より, α2A (t) と α2B (t) が最大値をとる時刻も温度によりそのオーダーが大
きく変化していることが確認出来る. 最大値を取る時刻は高い温度では同程度の値であ
るが, 温度が低くなると α2B の方が早く最大値を取るようになる. T = 0.44 では α2A (t) は
t = 203.636, α2B (t) は t = 27.5587 で最大値を取っており, 7.39 倍程度の差がある.
図 5.23 から α2A (t) が最大値をとる時刻について大まかな傾向はつかめたが, α2A (t) は大
きな揺らぎを含むためそれぞれの温度での平均的な振る舞いをとらえていない可能性が
ある. これを解決する方法としてはサンプル平均を取るのが一般的であるが, 本研究のよ
うな大規模計算では多数のアンサンブルを計算するのに非常に長い時間が必要となる. そ
こで今回は α2A (t) を便宜的に対数正規分布で近似し, 近似曲線の最大値を α2A (t) の極大点
tNGP とすることにする. 対数正規分布のフィッティングにより得られた tNGP と 5.23 に示
した α2A (t) の極大点の比較を図 5.24 に示す. 対数正規分布を用いた近似により得られた
結果, 破線が元の極大点である. フィッティングを行なったことによりばらつきが抑えら
れていることが分かる. 対数正規分布を用いる理由には物理的根拠は無いが, 以降は tNGP
を使用する.
5.2.7
自己中間散乱関数
図 5.25 に様々な温度における自己中間散乱関数の時間発展を示す. 波数に関しては
S AA (k) の第一ピークが現れる波数, k = 7.25 を用い, t = 5000 まで予備計算を行った後に
計測した. 尚, 計算に際しては (4.62) 式を用いた.
55
第 5 章. 計算結果と考察
5.2 2成分系
図 5.25: 様々な温度における自己中間散乱関数 FSA (k, t) の時間発展. 実線が計算機実験の
結果であり, 左から順に温度 T = 5.00, 3.00, 2.00, 1.50, 1.00, 0.90, 0.80, 0.70, 0.60, 0.55,
0.50, 0.47, 0.45, 0.44. 破線は自由拡散運動における解析解 FSA (k, t) = exp(−k 2 D0 t)
ごく初期の時間領域 (t ≤ 10−3 程度) では, FSA (k, t) の振る舞いは全ての温度においてほ
とんど等しいものとなっている. いま (4.66) 式において k 4 以上の項がごく小さいと仮定
すると,
¶
M2A (t) 2
k .
' exp −
6
µ
FsA (k, t)
(5.4)
と書ける. t ≤ 10−3 においては M2A (t) ' 6D0 t であったから,
¢
¡
FsA (k, t) ' exp −D0 tk 2 .
(5.5)
となる. これは一粒子の自由ブラウン運動における FSA (k, t) に対応している. (5.5) 式に
基づく曲線を図 5.25 中の破線で示した. M2A (t) が 6D0 t と一致する時間領域程度では, よ
く振る舞いが一致していることが分かる.
10−3 ≤ t ≤ 10−1 では, FsA (k, t) の現象が図 5.25 の上で直線的になっていることが分か
る. これは平均場理論により予測されている対数関数的な減衰である [48].
(
次に, FSA (k, t) に対して (5.4) 式がどの程度適用出来るか検討する. FS k, t) の時間発展
と今回計算した M2A (t) を (5.4) 式に代入して描いた曲線を比較したのが図 5.26 である. 図
中の実線が各温度において得られた M2A (t) のデータを用いて描いた近似曲線であり, 実際
56
第 5 章. 計算結果と考察
5.2 2成分系
図 5.26: 自己中間散乱関数の計算結果と M2A (t) を用いた曲線の比較.
の FSA (k, t) の値を中抜きの点で示してある. (5.5) 式に比べ, より長い時間領域で FSA (k, t)
と振る舞いが一致していることが分かる. 振る舞いが一致している時間領域は温度によっ
て異なるが, FSA (k, t) = 0.8 となる付近までは (5.4) 式でよく一致している. FSA (k, t) の値
が 0.8 程度よりも小さくなると, 曲線の値が FSA (k, t) よりも小さくなるのが観察される.
この傾向は低温になるほど顕著である. この近似が成り立つ程度の時間領域においては
k 4 以上の項の値は小さく, 以降の時間領域では温度が低いほどその影響が大きいことが
分かる. これら α2A (t) の最大値が低温になるほど大きく, 最大値が表れる時間が遅くなる
という結果と対応する.
次に, 次数を上げた k 4 までの式,
"
M A (t) 2 1
k +
FSA (k, t) ' exp − 2
6
2
µ
M2A (t)
6
¶2
#
α2A (t)k 2 .
(5.6)
がどの程度まで FSA (k, t) と一致するか検証する. 図 5.27 は M2A (t) と α2A (t) の計算結果を
(5.6) 式に代入して描いた曲線と FSA (k, t) を比較したものである. (5.6) 式による曲線を
実線, FSA (k, t) を中抜きの点で示した. 図 5.26 と比べると, 曲線が一致する時間領域がよ
り長くなっていることが見て取れる. 特に T ≥ 1.0 については FSA (k, t) が 0 になるまで
両者はよく一致しており, より k の高次の項からの緩和過程への寄与が小さいことが分
かる. FSA (k, t) と近似曲線が一致しなくなる点については温度によって大きくことなる.
57
第 5 章. 計算結果と考察
5.2 2成分系
図 5.27: 自己中間散乱関数の計算結果と M2A (t), α2A (t) を用いた曲線の比較.
FSA (k, t) の値でみると FSA (k, t) = 0.5 ∼ 0.8 程度であり, 温度が低くなるほど値が高い所
から一致しなくなっている. 時間領域は温度が低いほど長い領域で一致する. これらの傾
向は図 5.26 と同じであり, 近似が成り立つ領域は FSA (k, t) の何らかの特性時間と関係し
ていると予想される. 図 5.26 と異なるのは, 時間が経過すると曲線の値が FSA (k, t) よりも
大きくなる点である. この傾向は温度が低くなるほど顕著になり, T ≤ 0.55 では曲線が 0
に緩和しなくなっている. α2A (t) は最大値を取った後 0 に収束するが, 低温では最大値を
とる時間および 0 に収束するまでの時間が遅くなるため, T ≤ 0.55 では (M2A (t))2 の値が
これを打ち消すほど大きくなっており, 曲線のこの様な振る舞いを引き起こしていると考
えられる. k 6 以降の項はこれを打ち消すような形で FSA (t) に入っていると予測される.
5.2.8
平均場理論による解析
M2A (t) について, 平均場理論を用いた解析を行う. 図 5.28 に, 様々な温度の M2A (t) と,
それぞれの温度のデータに対して (4.74) 式によるフィッティングを行って得られた平均場
理論の曲線を示す. 広い温度領域, 時間領域で平均場理論の式と計算機実験の結果が一致
していることが分かる. しかしながら, 温度が低くなると初期の拡散運動からケージ効果
により M2A (t) の傾きが減少する付近で平均場理論の値よりも計算機実験の値が大きくな
58
第 5 章. 計算結果と考察
5.2 2成分系
図 5.28: 平均場理論による平均二乗変位 M2A (t) のフィッティング.
るのが観察される. 特に, T = 0.45, 0.44 の低い温度のデータは t ≤ 1 程度の領域で明らか
なずれが見られる. これらの温度では自己中間散乱関数の結果からも分かるように緩和
時間が非常に長くなっているため, 系がまだ非平衡状態である可能性が高い. これが短時
間領域でのずれを生み出しているのではないかと考えられる.
L(A)
次に, 平均場理論のフッティングパラメータである長時間自己拡散係数 DS
と自由長
(A)
lf の温度依存性を図 5.29, 図 5.30 にそれぞれ示す. どちらも横軸は温度の逆数 1/T であり,
L(A)
DS
については値の変化が非常に大きいため常用対数を取ったものを示した. 図 5.29 中の
実線は平均場理論による拡散係数の理論式 (特異関数型) である. 後に詳しく述べるが, この
理論式のパラーメータは同じ相互作用を持った分子動力学法による原子系のシミュレーショ
ンにおける DSL のフィッティングパラメータを適当に変換して得られたものである. 図 5.29
では, 温度が低くなるにつれ長時間拡散係数の値が連続的, かつ急激に減少している様子が
見てとれる. 例えば T = 1.00 と T = 0.50 の値を比較すると (4.930×10−2 )/(7.291×10−4 ) =
59
第 5 章. 計算結果と考察
5.2 2成分系
図 5.29: A 粒子の長時間自己拡散係数
L(A)
DS の温度依存性.
図 5.30: 自由長 lf の温度依存性.
67.61, T = 0.90 と T = 0.45 の値では (3.405 × 10−2 )/(7.852 × 10−5 ) = 433.7 と, 温度が半
分になる間に実に数十分の一から数百分の一のオーダーにまで減少している. T > 1.00
の温度が高い領域では DSL が理論式を下回っており, ずれが生じている. (4.71) 式はガラ
ス転移点近傍に特異点を仮定し, その周りで展開して得られたということを考えると, 展
開した点から離れているために近似が成り立たなくなっているのではないかと推測され
L(A)
る. また, T ≤ 0.50 の低い温度では理論線よりも DS
の値が大きくなっている. 図 5.28
からも分かる通り, これらの温度はまだ非平衡の可能性があり, さらに計算を継続すると
拡散係数が減少すると考えられる. T = 0.45, 0.44 の二つの温度に関しては特異点よりも
L(A)
低い温度であるから理論線と対応するには DS
が 0 になる必要がある. しかし, 今回の
計算で得られた M2A (t) の結果から T = 0.45, 0.44 においても長時間拡散運動は確認され
ており, これが実現されるとは考えにくい.
(A)
図 5.30 でも, 温度が低くなるにつれ lf
L(A)
かる. 値が変化する範囲は DS
5.2.9
の値が連続的かつ急激に変化していることが分
ほど大きくは無く, 全温度を通じて 0.1 程度の値である.
時間スケール
図 5.31 に平均場理論の平均二乗変位の特性時間である tf , tβ , tL , および A 粒子のノン
ガウシアンパラメタ α2A (t) が最大値を取る時間 tNGP の温度依存性を示す. それぞれの特
性時間の計算には 5.29 と図 5.30 の値を用いた. tNGP については, 図 5.23 に示したもので
60
第 5 章. 計算結果と考察
5.2 2成分系
図 5.31: 平均場理論の特性時間とノンガウシアンパラメタの最大値をとる時刻の比較. 上
から順に tL , tNGP , tγ , tβ , tf である.
ある. まず M2A (t) の特性時間についてであるが, tf は温度が低くなると減少するのに対
し, tγ , tβ , tL は増加している. 値の変化の幅は tβ が最も大きく, その次が tL であり, tf は
最も変化の幅が小さい. このため, 温度が低くなると tL と tβ の比が小さくなっている. こ
れは温度が低いほど β 緩和領域から長時間拡散領域への移行が短い時間で行われている
ということを意味する. T = 5.0 では tf と tL の比は 102 程度であるが, = 0.44 では 106 程
度まで広がっている. tf は ケージ効果が見られ始める時間であり,tβ は運動が長時間拡散
に移行し始める時間であるから, 温度の低下によりケージ効果が見られる領域, すなわち
β 緩和領域が温度の低下とともに非常に長くなっていることが確認出来る.
tβ と tNGP は粒子の運動が長時間拡散過程に移行する時間スケールを特徴づけるもので
ある. 図 5.31 において, 両者は別々の物理力から測定したものであるにも関わらず, その
値および振る舞いが近いことが見て取れる. これを直接比較したものが図 5.32 である. 図
中の破線は tβ = tNGP である. 高い温度では両者はほぼ等しく, 温度が低くなると tNGP の
61
第 5 章. 計算結果と考察
5.2 2成分系
図 5.32: tβ と tNGP の比較. 縦軸が tβ , 横
軸が tNGP である. 破線は tβ = tNGP を
示す.
図 5.33: A 粒子の長時間自己拡散係数
L(A)
DS と tNGP の比較. 破線はべき乗則に
よる近似曲線
方が tβ よりも大きくなるのが見て取れる.
L(A)
tNGP について A 粒子の長時間自己拡散係数 DS
との関係も調べた. これを図 5.33 に
示す. 図中の破線はべき乗則,
tNGP = AxB
(5.7)
による近似曲線であり, 最小二乗法による解析より A = 0.019523, B = −0.906933 と求め
られた. 広い温度でべき乗則が成り立っていることが分かる.
さらに, 各特性時間の自由長 lf に対する依存性を調べた. これを図 5.34 に示す. 横軸は
1/lf , 縦軸は特性時間の対数である. lf は温度が低くなるほど小さくなるから, 図の右側
の点ほど温度が低い状態に対応している. 温度が低くなるにつれ tf は減少し, tf 以外の時
間は大きくなっているのが見て取れる.
tf の lf に対する依存性はその定義から,
µ
¶
lf2
log10 (tf ) = log10
D0
= −2 log10 (lf−1 ) − log10 (D0 )
(5.8)
である.
他の特性時間については T = 0.70 付近からはそれより高い温度に比較して tf 以外の時
間の変化が大きくなっていることが分かる. 特に, T = 0.70, 0.60, 0.55, 0.50 における特
62
第 5 章. 計算結果と考察
5.2 2成分系
図 5.34: 平均場理論の特性時間とノンガウシアンパラメタの最大値をとる時刻の自由長
lf 依存性. 上から順に tL , tNGP , tγ , tβ , tf である. 温度は左から順に T = 5.00, 3.00, 2.00,
1.50, 1.00, 0.90, 0.80, 0.70, 0.60, 0.55, 0.50, 0.47, 0.45, 0.44 に対応している.
性時間の増加が顕著である. T = 0.70 は温度は A 粒子の動径分布関数 g AA (r) の第二ピー
クに割れが出現する温度であり, 過冷却液体は通常の液体に比べて緩和過程が複雑である
ということを考えると, この特性時間の lf 依存性の変化も過冷却を特徴づけるものであ
ると推測することが出来る.
63
第 5 章. 計算結果と考察
5.3 分子動力学法による計算結果との比較
図 5.35: ブラウン動力学法と分子動力学法における圧力, 有効圧力 P の比較. 黒丸がブラ
ウン動力学法における有効圧力, 白丸が分子動力学法における圧力である.
5.3
5.3.1
分子動力学法による計算結果との比較
圧力と動径分布関数
ここでは同一の相互作用を持った分子動力学法での計算機実験との比較を行なう. 図
5.35 は本研究で得られた有効圧力を鳴海らによる Kob-Andersen 型レナード・ジョンズ系
における計算機実験から得られた圧力と比較したものである [37]. 横軸が逆温度 1/T , 縦
軸が圧力および有効圧力である. コロイド分散系と原子・分子系というように系の時空間
スケールや初期運動が全く異なる系における計算であるにも関わらず, 両者の温度依存性
は非常に良く一致している. ブラウン動力学における粒子の初期運動はブラウン運動, 分
子動力学法では力学運動であるが, それらの違いは圧力には現れていないことが分かる.
単分散系における考察で見たように圧力は系の構造と対応していることを考えると, 二つ
の系の静的な構造は温度 T が等しければ一致しているのではないかと予測出来る. これ
64
第 5 章. 計算結果と考察
5.3 分子動力学法による計算結果との比較
T = 1.00
T = 0.50
図 5.36: T = 1.00, 0.50 におけるブラウン動力学法, 分子動力学法の動径分布関数の比較.
実線がブラウン動力学法の計算結果, 白丸が分子動力学法の計算結果.
を検証するために, 二つの系について同じ温度における動径分布関数を比較したのが図
5.36 である. 先の予測の通り, T = 1.00, 0.50 のどちらもブラウン動力学法と分子動力学
法の計算結果が非常に良く一致した. わずかなずれに関しても誤差の範囲内であると考
えられる.
有効圧力と圧力, および動径分布関数の計算結果がほぼ完全に一致したことから, 同じ
相互作用を持ったブラウン動力学法と分子動力学法において静的な構造は温度のみによっ
て決定され, 初期運動の違いは影響を与えないということが分かった. 続いて次節では,
系の動的特性を特徴づける自己長時間拡散係数と平均二乗変位を比較し, 二つの系の違い
についてさらに考察を行う.
5.3.2
自己長時間拡散係数と平均二乗変位
二つの計算は異なる系について行われており, 適当に単位量をスケールし直す必要があ
る. 長さの単位量 d0 はどちらも共通で,
d0 =
σAA
[BD, MD],
2
65
(5.9)
第 5 章. 計算結果と考察
5.3 分子動力学法による計算結果との比較
L(A)
図 5.37: 分子動力学法とブラウン動力学法における粒子 A の長時間自己拡散係数 DS
の計算結果の比較. 黒丸がブラウン動力学法の結果, 白丸が分子動力学法の結果である.
時間の単位量 t0 はそれぞれ,
 2
σAA


[BD]

 4D0
t0 =
r
2


mσAA


[MD],
12²AA
(5.10)
を用いて無次元化した.
図 5.37 にブラウン動力学法と分子動力学法の A 粒子の長時間自己拡散係数の比較を示
す. 二組のデータが示されているが, 黒丸がブラウン動力学法, 白丸が分子動力学法の計
算結果である. 図中の曲線は平均場理論によるフィッテングである. まず, 分子動力学法
の結果については, パラメータは κ = 46.77, 1/Tc = 2.1285 となっており, 以下の理論線
でフィッティングを行った.
DSL MD (T ) =
µ
1+κ
66
D0
¶µ
¶−2 .
Tc
Tc
1−
T
T
(5.11)
第 5 章. 計算結果と考察
5.3 分子動力学法による計算結果との比較
L(A)
図 5.38: 分子動力学法における粒子 A の長時間自己拡散係数 DS と, ブラウン動力学
L(A)
法から得られた DS を短時間自己拡散係数 DSS で割ったものの比較. 詳細は図 5.37 に
同じ.
1/T ' 0.5 付近から特異点近傍まで理論線と計算結果がよく一致していることが分かる.
次に, 本研究による計算結果に対する平均場理論の式は,
DSL BD (T ) =
1+κ
DSS
µ
D0
D0
¶µ
¶−2 .
Tc
Tc
1−
T
T
(5.12)
で与えられる. ここで DSS は仮想的に導入された短時間自己拡散係数で, 本研究のモデル
の体積分率が 0.59 であることから DSS = 0.27 としている. 残り二つのパラメータについ
ては, 分子動力学法の場合と同じ値, κ = 46.77, 1/Tc = 2.1285 を用いている. 1/T ' 1.0
付近から特異点近傍まで理論線と計算結果がよく一致していることが見て取れる.
驚くべきことに二つの計算機実験は全く違うスケールの現象を取り扱っているにも関
わらず, κ と 1/Tc の値が同じ曲線でフィッティングすることが出来ることが分かる. しか
も違いは DSS という初期運動に関係した因子にのみ依存している. DSS による違いを取り
67
第 5 章. 計算結果と考察
5.3 分子動力学法による計算結果との比較
図 5.39: ブラウン動力学法と分子動力学法における A 粒子の平均二乗変位の比較.
除くために, ブラウン動力学法の結果と理論線をを DSS = 0.27 でスケールして比較したも
のが図 5.38 である. 温度が高い領域での振る舞いは異なるが, T ≥ 1.0 では両者がほぼ一
致していることが分かる. これは, 温度が高い領域では粒子同士の相互作用に加え初期運
動が長時間拡散運動に影響しており, 逆に温度が低い領域では粒子同士の相互作用が長時
間拡散を支配しているためと考えられる.
次に, 系の動的特性は長時間自己拡散係数が支配しているという平均場理論の視点に
基づき, 同じ拡散係数を持った分子動力学法の平均二乗変位と, ブラウン動力学法の平
均二乗変位の比較を行う. 具体的には T = 0.70[BD] と T = 2.00[MD], T = 0.50[BD] と
T = 0.714[MD] の二組の温度について比較を行う. ここで [BD] はブラウン動力学法, [MD]
は分子動力学法における値を示す. 比較した結果を図 5.39 に示す. 黒丸がブラウン動力
学法における平均二乗変位, 白抜きの点が分子動力学法における平均二乗変位である. ブ
√
ラウン動力学法は t0 , 分子動力学法は t0 T でそれぞれ時間スケールを無次元化してある.
68
第 5 章. 計算結果と考察
5.3 分子動力学法による計算結果との比較
初期運動はブラウン動力学法はブラウン運動, 分子動力学法は力学運動であり,
(
t [BD]
M2A (t) ∝
(5.13)
t2 [MD],
にそれぞれ従っており, 全く異なる運動である. にも関わらず, 二組の温度において, β 緩
和過程以降の両者の振る舞いは非常に良く一致していることが分かる. これは初期運動
を除いた系の動的特性は長時間自己拡散係数の値により支配されるという平均場理論の
主張と一致するものである.
圧力と動径分布関数は温度が同じであれば二つの系でほぼ完全に一致したが, 動的特性
を支配する長時間自己拡散係数は DSS の分だけブラウン動力学法の方が大きかった. これ
らは徳山らにより報告された剛体球流体による計算機実験と対応するものである [49, 50].
ガラス転移現象は系の動的特性と密接に関わった現象であるから, 二つの計算を比較する
と本研究の系の方がガラス転移温度が低くなるのではないかと考えられる.
69
第 5 章. 計算結果と考察
5.3 分子動力学法による計算結果との比較
70
第 6 章 結言
本研究ではガラス転移現象の解明に資するべく, レナード・ジョンズコロイド分散系に
おいて大規模なブラウン動力学法による計算機実験を遂行した. 粒子数が 10000 程度か
つ長時間の計算機実験は行われている数も少なく, 有限サイズ効果も小さいので非常に有
用なデータが得られたと考えている. 計算は単成分系と2成分系の二つの系で行われた.
計算結果としては大きく分けて有効圧力や動径分布関数等の静的物理量, 平均二乗変位や
自己中間散乱関数等の動的な物理量を測定し, 考察を行った. また, 計算結果を同一の相
互作用を持った分子動力学法の計算機実験と平均場理論の立場から比較することにより,
様々な知見が得られた. 以下にそれらをまとめる.
6.1
単成分系と2成分系の比較
単成分系
単成分系では温度が低くなると有効圧力に急激な変化が観察された. 動径分布関数か
ら, 系の構造は有効圧力に急激な変化が見られた領域よりも高温では液体, 低温では結晶
になっていることが確認された. すなわち, 単成分系では液体-結晶転移が観測された.
2成分系
2成分系では有効圧力, 動径分布関数ともに温度に対する変化は連続的であり, 単成分
系に見られたような急激な変化は観測されなかった. 動径分布関数も温度変化に対し連
続的に変化してた. このことから2成分系では成分を増やすことにより結晶化を避ける
ことが出来たと考えられる.
2成分系では液体-結晶転移を避けたことにより, 全ての温度で平均二乗変位から長時
間自己拡散が観察された. 自己中間散乱関数も単成分系のように平坦になることはなく,
温度が低くなるにつれ緩和に要する時間は増大するものの, 1 から 0 への緩和が見られた.
71
第 6 章. 結言
6.2
6.2 緩和過程の複雑化
緩和過程の複雑化
2成分系では温度が低くなるにつれ動的特性が変化し, 緩和過程が複雑になる様子が観
察された. 温度が低くなるとノンガウシアンパラメタの値も 1∼5 程度と非常にその値が
大きくなり, 粒子の変位のガウス分布からのずれが大きくなることが見られた. 自己中間
散乱関数から観察される系の緩和過程も温度が低くなるにつれ二段階緩和が見られるよ
うになった. β 緩和過程においては, 平均場理論により予測されていた対数関数的な減衰
が観察された. これは大規模な系で計算を行ったことにより有限サイズ効果が比較的少
なかったためと推測される.
6.3
A 粒子と B 粒子の比較
2成分系において, A 粒子と B 粒子で静的特性, 動的特性に違いが見られた.
主成分である A 粒子の動径分布関数は全ての温度で液体の様な特徴が見られたが, T =
0.70 ∼ 0.60 付近から第二ピークに割れが出現した. 対して, 高い温度では B 粒子はほぼ
一様に分布していたのに対し, 低温では B 粒子同士が隣合わない傾向が強くなるという
変化が見られた. また, 低い温度では A 粒子と B 粒子の短距離相関が強く現れた.
粒子の運動に関しては, B 粒子の方が A 粒子よりも平均二乗変位および拡散係数が大
きい値をとり, その差は温度が低くなるほど大きくなることが分かった. ノンガウシアン
パラメータに関しては B 粒子の方が最大値が大きく, 最大値をとる時刻も早いことが分
かった.
6.4
平均場理論による解析
2成分系において, 平均場理論による解析を行った.
平均二乗変位について平均場理論でフィッティングを行ったところ, 両者はよく一致し
た. T = 0.45, 0.44 に関しては β 緩和過程の初期において平均場理論よりも計算結果の方
が値が大きくなったが, これは系が非平衡であるためと考えられる.
長時間自己拡散係数は温度が低くなるにつれ急激に減少した. その温度特性は平均場
理論の理論式と T = 0.50 付近までよく一致することが確かめられた.
72
第 6 章. 結言
6.5
6.5 特性時間の変化
特性時間の変化
平均場理論の特性時間とノンガウシアンパラメータが最大値をとる時間 tNGP について
比較を行い, 考察を行った.
平均場理論の特性時間の温度依存性から, 温度が低くなると系の β 緩和領域が劇的に増
大することが確認された. 特に β 緩和過程から長時間拡散過程への移行の指標である tβ
の増大が顕著であった.
tβ と tNGP は高い温度ではほぼ一致し, 温度が低くなると tNGP の方が大きくなる傾向
L(A)
が見られた. tNGP と A 粒子の長時間ん自己拡散係数 DS
の関係をべき乗則で調べたと
ころ,
L(A) −0.90
tNGP ∼ (DS
)
,
(6.1)
でよく近似出来るという結果が得られた.
また, 各種特性時間の自由長 lf に対する依存性を調べたところ, T = 0.70 ∼ 0.50 の温
度に対応する領域でで tγ , tβ , tNGP , tL の値が急激に上昇することが分かった. T = 0.70
は g AA の第二ピークに割れが出現する温度に対応しており, この特性時間の lf 依存性の
変化も過冷却液体を特徴づけるものである可能性がある.
6.6
分子動力学法による計算機実験との比較
同じ相互作用, 数密度を持つ分子動力学法の計算機実験と結果を比較した.
静的な物理量の温度依存性について有効圧力と動径分布関数を用いて比較を行ったと
ころ, 両者ほぼ完全に一致した. このことから系の静的な性質に初期運動の違いは表れな
いということが予測される.
平均場理論の視点に基づいて A 粒子の自己長時間拡散係数を比較した所, 特異点が一
致し, さらにもう一つのパラメータの差は短時間自己拡散係数 DSS の値として表れた. ま
た, 長時間自己拡散係数が同じ値を持つ温度で平均二乗変位を比較したところ, β 緩和過
程後期以降の振る舞いが非常に良く一致した. これは系の動的な性質は長時間自己拡散
係数が支配しているという平均場理論の主張と対応している. これはもし初期拡散 D0 を
DSS で置き換えて計算を置き換えれば, 二つの系の結果は全く一致するということを意味
しており, 流体力学的相互作用のコロイド分散系の動的特性ににおける重要性, およびガ
ラス転移が普遍的な現象であるという視点の正当性を示すものである.
73
第 6 章. 結言
6.7
6.7 最後に
最後に
大規模なブラウン動力学法による計算機実験によって, 低温液体状態における緩和過
程の複雑化, 温度低下に伴う特性時間の急激な変化等, 過冷却液体の特徴的な振る舞いと
されているものを計算機上で再現することが出来た. また, 同一の相互作用を持った分子
動力学法による計算機実験との比較から, ガラス転移現象, 特に β 緩和過程後期以降のの
普遍性, および流体力学的相互作用の重要性を確認することが出来た. 本研究を通じて今
後のガラス転移現象の解明のために非常に有用なデータを蓄積することが出来たと言え
よう.
74
付 録A
A.1
ランジュバン方程式の性質
ランジュバン方程式,
d
X(t) = V (t)
dt
d
m V (t) = −γ α V (t) + R(t).
dt
hR(t)R(t0 )i = 2γkB T δ(t − t0 )1.
(A.1)
(A.2)
(A.3)
を解き, ブラウン運動の基本的な性質を確認する [51]. (A.2) 式は一階の線形微分方程式
と見ることが出来るから, 0 から t まで積分して,
−t/tB
V (t) = e
1
V (0) +
m
Z
t
dse−(t−s)/tB R(s).
(A.4)
0
を得る. ここで,
tB =
m
,
γ
(A.5)
とおいた.(A.2) 式の形から明らかであるように tB はブラウン運動の特性時間の一つであ
り, ブラウン緩和時間と呼ばれる. 平均速度 hV (t)i は,
1
hV (t)i = e
hV (0)i +
m
−t/tB
= e
hV (0)i .
−t/tB
Z
t
dse−(t−s)/tB hR(s)i
0
(A.6)
となる. ある時刻における速度の情報は tB 程度時間が経過すると失われてしまうことが
分かる.
次にブラウン粒子の運動エネルギー,
®
1 ­
K(t) = m V (t)2 ,
2
75
(A.7)
付録 A.
A.1 ランジュバン方程式の性質
の時間特性を考える.K(t) を時間で微分すると,
À
¿
d
dV (t)
K(t) = m V (t) ·
dt
dt
­
®
2
= −γ V (t) + hV (t) · R(t)i
2
= − K(t) + hV (t) · R(t)i .
tB
(A.8)
となる. ここで右辺第二項は,
Z
1 t
< V (s)R(s) =
dse−(t−s)/tB hR(s) · R(t)i
m 0
Z
kB T t
dse−2(t−s)/tB δ(s − t)
=
tB 0
kB T
= d
.
tB
(A.9)
と計算出来るので結局,
d
kB T
2t
.
K(t) = − K(0) + d
dt
tB
tB
(A.10)
となる. これを 0 から t まで積分すれば,
−2t/tB
K(t) = e
Z
t
K(0) +
dse−2(t−s)/tB d
0
kB T
tB
d
= e−2t/tB K(0) + (1 − e−2t/tB ) kB T.
2
(A.11)
右辺第一項は tB /2 程度で緩和するため, tB より十分大きい時間では K(t) = d kB2T となる
ことが分かる.
最後にブラウン粒子の平均二乗変位,
®
­
M2 (t) = [X(t) − X(0)]2 .
(A.12)
の振る舞いを考える. まず, 時刻 t におけるブラウン粒子の位置 X(t) は V (t) を積分すれ
ば得られる.
Z
t
X(t) = X(0) +
"
−s0 /tB
ds0 e
0
1
+
m
Z
s0
#
−(s0 −s)/tB
dse
R(s)
0
= X(0) + (1 − e−t/tB )V (0)tB
Z s1
Z
1 t
+
ds2 e−(s1 −s2 )/tB R(s2 ).
ds1
m 0
0
(A.13)
故に,
hX(t)i = hX(0)i + (1 − e−t/tB )tB hV (0)i .
76
(A.14)
付録 A.
A.2 ガウス積分
である. また M2 (t) は,
­
®
M2 (t) = (1 − e−t/tB )2 t2B V (0)2
Z s1
Z t
2
−t/tB
+ (1 − e
ds2 e−(s1 −s2 )/tB R(s2 ) · V (0)
ds1
)tB
m
Z t
Z s1
Z0 t
Z0 u1
1
ds1
ds2
du1
du2 e1(s1 +u1 −s2 −u2 )/tB R(s2 ) · R(u2 ).(A.15)
2
m 0
0
0
0
となる. (A.15) 式の右辺第一項は hV (0)i = dkB T /m を用いると,
®
­
dkB T
(1 − e−t/tB )2 t2B V (0)2 = (1 − 2e−t/tB + e−2t/tB )t2B
m
¡
¢
−t/tB
−2t/tB
+e
= dD0 tB 1 − 2e
.
(A.16)
となる. ここで D0 := kB T /γ とした. (A.15) 式の右辺第二項は hR(s) · V (0)i = 0 より 0
となる. 第三項は,
=
=
=
=
Z t
Z s1
Z t
Z u1
1
ds1
ds2
du1
du2 e−(s1 +u1 −s2 −u2 )/tB R(s2 )R(u2 )
m2 0
Z t
Z0 t
Z 0t
Z 0t
1
ds1
ds2
du1
du2 e−(s1 +u1 −s2 −u2 )/tB Θ(s1 − s2 )Θ(u1 − u2 )2dB T δ(s2 − u2 )
2
m 0
Z t 0 Z t 0 Z t 0
2dkB T
du1 e−(s1 +u1 −2s2 )/tB Θ(s1 − s2 )Θ(u1 − s2 )
ds2
ds1
mtB 0
Z t
Z0 t
Z0 t
2dkB T
ds2
ds1
du1 e2s2 /tB e−(s1 +u1 )/tB
mtB 0
s2
s2
¡
¢
−t/tB
− tB e−2t/tB .
(A.17)
dD0 2t − 3tB + 4tB e
となる. ここで Θ(t) はステップ関数である. 以上により平均二乗変位 M2 (t) は,
£
¤
M2 (t) = 2dD0 (t − tB ) + tB e−t/tB .
(A.18)
となる.
A.2
ガウス積分
ここではガウス積分,
Z
I2n =
∞
dxx2n e−ax
2
(n = 1, 2, . . . )
−∞
の一般的な形を求める. ここで a は正の定数である.
77
(A.19)
付録 A.
A.3 一様乱数と正規乱数
まず, I0 については,
½Z
I02
=
∞
−ax2
¾ ½Z
¾
∞
dxe
dye
−∞
Z
Z ∞
2
2
dx
dye−a(x +y )
−∞
−∞
Z 2π Z ∞
2
dθ
dr re−ar
0
0
#∞
"
−ar2
e
2π −
2a
0
π
a
−ay 2
−∞
∞
=
=
=
=
より,
r
I0 =
π
a
(A.20)
(A.21)
と計算出来る.
次に,
∂I0
=
∂a
Z
∞
∂ −ax2
e
∂a
−∞
Z ∞
∂
2
= −
dx x2 e−ax
∂a
−∞
= −I2
dx
である. 一方,
∂I0
∂
=
∂a
∂a
r
π
1
=−
a
2
であるから,
1
I2 =
2
となる. 同様に計算すれば,
I2n
r
r
π
a3
π
.
a3
(2n − 1)!!
=
2n
r
(A.22)
(A.23)
(A.24)
π
a2n+1
(A.25)
と求まる.
A.3
一様乱数と正規乱数
ランジュバン方程式の差分式に現れる確率変数 u(t) は平均 0 分散 1 の正規乱数である.
実際に計算を行う際には毎ステップ正規乱数を発生させるのは計算負荷が大きいので一
78
付録 A.
A.4 初期配置の作製法
様乱数 ξ で代用する. ここで ξ の確率分布 pξ (ξ) はパラメータ hξi, wξ を用いて,
(
1
|ξ − hξi | ≤ wξ ,
2wξ
ξ=
0
|ξ − hξi | > wξ .
(A.26)
と書ける. hξi は ξ の平均値, wξ は分布の広がりを表すパラメータである. その際に ξ の
分散が 1 になるようにパラメータを決定せねばならない.
まず, 明らかに hξi = 0 である. ξ の分散 σξ2 を計算してみると,
Z ∞
2
σξ =
dξξ 2 p(ξ)
Z−∞
wξ
1
dξξ 2
=
2wξ
−wξ
· 3 ¸wξ
1 ξ
=
2wξ 3 −wξ
wξ2
=
.
3
(A.27)
となる. σξ2 = 1 となれば良いので結局,
wσ =
√
3.
を得る. ランジュバンを数値的に解く際は正規乱数 u を幅
(A.28)
√
3, 平均 0 の一様乱数で置き
換えて計算出来ることが分かる.
A.4
初期配置の作製法
希薄な場合 (一様乱数を用いた方法)
ランダム配置の作成方法として最も単純なのが一様乱数を用いた方法である. 直感的
にも分かり易く, 比較的希薄な初期配置 (φ ∼ 0.30 程度以下) を作る場合はこの方法で十
分である. 以下に大まかな手順を紹介する.
例として, 立方体の単位セル内に半径 ai (i = 1,) の粒子 N 個, 体積分率 φ の初期配置を
作成する場合を考えよう. 各粒子は球状で, セルには周期境界条件が適用されるものと
する.
1. システムサイズ L の決定
体積分率の定義式,
PN
φ=
79
i
4
πa3i
3
L3
(A.29)
付録 A.
A.4 初期配置の作製法
より, L を決定する. すなわち,
Ã
L=
4π X 3
a
3φ i i
N
!1/3
.
(A.30)
2. i 番目の粒子位置の仮決定.
i 番目の粒子の仮の位置 Xtmpi を乱数を用いて与える.ξ を (0, 1) の一様乱数とすると,
Xtmpi = (L − 2a)ξex + (L − 2a)ξey + (L − 2a)ξez
(A.31)
とする. ここで ex , ey , ez はそれぞれ x, y, z 方向の単位ベクトルである. これによ
り x, y, z それぞれの方向において一様な確率分布で粒子位置を決定する事が出来
る. 直径分だけ減じているのは, 周期境界条件を考慮した場合, セルの境界付近で密
度が高くなってしまうのを防ぐためである.
3. j(< i) 番目の粒子との重なり判定粒子を一様に分布させる事が出来ても, 単に 2 の
手順を繰り返すだけでは粒子の重なりが生じる. これでは現実的な配置とは言いが
たい. そこで既に配置が決定している j(< i) 番目の粒子との重なり判定を行う. 具
体的には, 粒子間距離が ai + aj よりも大きいか判定する. すなわち判定式は,
|Xtmpi − Xj | ≥ (ai + aj ).
(A.32)
判定式が,
• すべての j(< i) について真.→4 へ移動.
• いずれかの j(< i) について偽.→2 へ移動.
4. i を i + 1 に置き換えて, 2 ∼ 3 を i = N となるまで繰り返し
体積分率が高い場合
過冷却状態, ガラス状態をの計算機実験を行う場合, 体積分率の高い (φ ≥ 0.50) 初期配
置を作る必要がある.A.4 の方法では,φ が 0.30 程度より大きくなると粒子の重なりが頻繁
に発生し, 初期配置を作るのが困難になる. 高い体積分率における初期配置の作成には,
A.4 の様な単純な方法ではなく, 工夫したアルゴリズムが必要である.
本研究では体積分率が高い場合の初期配置の作製法として, Jodrey と Tory によって開
発された方法 [41] を用いて, 高い体積分率におけるランダムな粒子配置を作成した. この
方法は徐々に半径を大きくしながら, 粒子の重なりを除去して行くのが特徴である.
80
付録 A.
A.5
A.5 平均場理論について
平均場理論について
(4.72) 式を M2 (t) について解いてみよう.
2
f (t) = eM2 (t)/lf
(A.33)
とおくと,
d
1 d
f (t) = f (t) · 2 M2 (t)
dt
lf dt
´
³
f (t)
L t/l2
L
S
L −2dDS
f
=
2dD
+
2d[D
−
D
]e
S
S
S
lf2
=
2dDSL
2d(DSS − DSL )
f
(t)
+
.
lf2
lf2
これは 1 階の線形微分方程式であるから,
(
Z
L t/lf 2
2dDS
f (t) = e
t
f (0) +
L t/l2
−2dDS
f
dse
0
2d(DSS − DSL )
lf2
(A.34)
)
·
¸t )
2
l
−
L
2
L
2
f
−
= e2dDS t/lf 1 +
e−2dDS s/lf
2
L
lf
2dDS
0
³
´
S
L
D −D
L
2
L
2
= e2dDS t/lf + S L S e2dDS t/lf − 1
DS
´
S ³
D
L
2
= 1 + SL e2dDS t/lf − 1 .
DS
(
2d(DSS
DSL )
(A.35)
故に,
M2 (t) = lf2 ln(f (t))
µ
½
¶¾
L
2
DSS
e2dDS t/lf
2
= lf ln 1 + L e
−1
,
DS
となり, (4.74) 式が得られる.
81
(A.36)
付録 A.
A.5 平均場理論について
82
参考文献
[1] 咲花済夫, トコトンやさしいガラスの本, 日刊工業新聞社 (2004).
[2] C.A. Angell, Nature 267, 1924 (1995).
[3] C. A. Angell, J. Non-Cryst. Soliids 102, 205 (1988).
[4] 日本化学学会, コロイド科学 I 基礎および分散・吸着, 東京化学同人 (1995).
[5] 日本化学学会, コロイド科学 III 生体コロイドおよびコロイドの応用, 東京化学同人
(1995).
[6] P. N. Pusey and W. van Megen, Nature 320, 340 (1986).
[7] W. van Megen, T. C. Mortensen and S. R. Williams, Phys. Rev. E 58, 6073 (1998).
[8] E. R. Weeks, J. C. Crocker, A. Schofield and D. A. Weitz, Science 287, 627 (2000).
[9] E. R. Weeks, D. A. Weitz, Phys. Rev. Lett 89, 095704-1 (2002).
[10] M. Kollmann, R. Hund, B. Rinn, G. Nägele, K. Zahn, H. König, G. Maret, R. Klein
and J. K. G. Dhont, Eulophys. Lett. 58, 919 (2002).
[11] H. König, R. Hund, K. Zahn, and G. Maret, Eur. Phys. J. E 18, 287 (2005).
[12] W. Götze, J. Phys. Condens. Matter 11, A1 (1999).
[13] M. D. Ediger, Annu. Rev. Chem 51, 99 (2000).
[14] C. A. Angell, K. L. Ngai, G. B. McKenna, P. F. McMillan, and S. W. Martin, J.
Appl. Phys. 88, 3113 (2000).
[15] P. G. Debenedetti and F. H. Stillinger, Nature 410, 259 (2000).
83
参考文献
[16] U. Bentzelius, W. Götze and A. Sjölander, J. Phys. C 17, 5915 (1984).
[17] E. Leutheusser, Phys. Rev. A 29, 2765 (1984).
[18] X. J. Han and H. Teichler, Phys. Rev. E 75, 061501 (2007).
[19] T. Voigtmann, A. M. Puertas, and M. Fuchs, Phys. Rev E 70, 061506 (2004).
[20] E. Flenner and G. Szamel, Phys. Rev. E 72, 031508 (2005).
[21] E. Flenner and G. Szamel, Phys. Rev. E 72, 011205 (2005).
[22] G. Szamel, Phys. Rev. Lett 90, 228301 (2003).
[23] K. Miyazaki, B. Bagchi, and A. Yethirai, J. Chem. Phys. 121, 8120 (2004).
[24] G. Biroli, J. Bouchaud, K. Miyazaki, and D. R. Reichman, Phys. Rev. Lett 97,
195701 (2006).
[25] G. Szamel and H. Löwen, Phys. Rev. A 44, 8215 (1991).
[26] M. Fuchs, W. Götze, I. Hofacker, and A. Latz, J. Phys.: Condens. Matter 3, 5047
(1991).
[27] M. Tokuyama, Physica A 387, 1926 (2008).
[28] M. Tokuyama and I. Oppenheim, Phys. Rev. E 50, R16 (1994).
[29] M. Tokuyama, Physica A 364, 23 (2006).
[30] M. Tokuyama, T. Narumi and E. Kohira, Physica A 385, 439 (2007).
[31] W. Kob and H. C. Andersen, Phys. Rev. Lett. 73, 1376 (1994).
[32] W. Kob and H. C. Andersen, Phys. Rev. E 51, 4626 (1995).
[33] W. Kob and H. C. Andersen, Phys. Rev. E 52, 4134 (1995).
[34] P. Gallo, R. Pellarin, and M. Rovere, Phys. Rev. E 67, 041202-1 (2003).
[35] G. Szamel, and E. Flenner, Phys. Rev. E 74, 021507 (2006).
84
参考文献
[36] E. Flenner, and G. Szamel, J. Phys. Condens. Matter 19, 205125 (2007).
[37] T. Narumi and M. Tokuyama, Rep. Inst. Fluid Sci. 19, 73 (2007).
[38] 米沢富美子, ブラウン運動, 共立出版 (1986).
[39] M. Tokuyama and I. Oppenheim, Physica A 216, 85 (1995).
[40] E. Kohira and M. Tokuyama, Rep. Inst. Fluid Sci. 19, 91 (2007).
[41] W. S. Jodrey and E. M. Tory, Phys. Rev. A 32, 2347 (1985).
[42] 青山幸男, チューニング技法入門, (2006).
[43] J. P. Hansen and I. R. McDonald, Theory of Simple Liquids (Academic, London,
1986).
[44] 上田顯, 分子シミュレーション, 裳華房 (2003).
[45] A. Rahman, Rhys. Rev. 136, A405 (1964).
[46] L. van Hove, Rhys. Rev. 95, 249 (1954).
[47] M. Tokuyama, Phys. Rev. E 62, R5915 (2000).
[48] M. Tokuyama, Mean-Field Theory of Glass Transitions (lecture note), (2007).
[49] M. Tokuyama, H. Yamazaki, and Y. Terada, Phys. Rev. E 67, 062403-1 (2003).
[50] M. Tokuyama, H. Yamazaki, and Y. Terada, Physica A 328, 367 (2003).
[51] 徳山道夫, ナノ物性物理学講義ノート, (2006).
85
参考文献
86
謝辞
論文の終わりに臨み, 日頃より懇切丁寧な御指導と御鞭撻を賜りました指導教官である
東北大学原子分子材料科学高等研究機構 東北大学流体科学研究所徳山道夫教授に謹んで
感謝の意を表します. 研究の方向性を導いて下さいましたことに関しても深く感謝の意
を表します.
修士論文中間審査会においてコロイド分散系の計算機実験に関して大変貴重なご意見
を賜りました東北大学流体科学研究所 小原拓教授に深く感謝の意を表します.
修士論文中間審査会において今後の研究の方向性に関して大変貴重なご意見を賜りま
した東北大学流体科学研究所 西山秀哉教授に深く感謝の意を表します.
計算機実験中の温度範囲の設定に関して貴重なご意見を賜りました東北大学工学研究
科 佐多教子准教授に深く感謝の意を表します.
計算機実験中の物理量に関しまして貴重なご意見を賜りました東北大学流体科学研究
所 米村茂准教授に深く感謝の意を表します.
計算機実験の手法および物理量に関して貴重なご意見を賜りました東北大学流体科学
研究所 徳増崇准教授に深く感謝の意を表します.
計算プログラム, 計算機実験の実施方法等につきまして多くのご助言を賜りました東北
大学流体科学研究所 寺田弥生助教に深く感謝の意を表します. 研究室に関わる事項全般
につきましてお世話になりましたことに大しても合わせて深く感謝の意を表します.
研究室での生活を支えて下さった, 東北大学流体科学研究所 高橋正嘉技官に深く感謝
の意を表します.
研究室での生活を支えて下さった, 東北流体科学研究所 工藤琢技官に深く感謝の意を
表します.
研究室での生活を支えて下さった, 徳山・寺田研究室秘書 湯村彩子さんに深く感謝の
意を表します.
統計物理学の理論や解析手法に関する疑問にいつも丁寧に答えて下さった徳山・寺田
研究室の日高邦昌先輩に深く感謝の意を表します.
87
謝辞
計算機実験の手法やコンピュータ全般の質問に対して親切にに教えて頂いた徳山・寺
田研究室の古林敬顕先輩に深く感謝の意を表します.
理論や解析手法, ガラス転移全般についての質問にいつも快く応じて下さった徳山・寺
田研究室の鳴海孝之先輩に深く感謝の意を表します. 本論文の執筆に際し, 貴重なデータ
を提供して下さったことに対して, 重ねて深く感謝の意を表します.
良い研究環境を作ってくれた徳山・寺田研究室の佐々木一誠君, 藤井宏之君, 鈴木洋平
君, 中西慎也君に深く感謝します. また, 本年度で徳山・寺田研究室を卒業する鈴木洋平
君, 中西慎也君の今後の活躍を期待します.
徳山・寺田研究室の皆様には公私にわたりお世話になり, 本論文の執筆に当たっても多
大な激励を頂きました. 重ねて感謝の意を申し上げます.
これまで不肖な私を励まし, 支えてくれた家族, 友人達に心より感謝の意を表します.
最後になりましたが, 本論文の数値計算は東北大学流体科学研究所未来流体情報創造
センターのスカラー並列計算システム SGI Altix3700Bx2 を用いて行なわれました. スー
パーコンピューター無しでは本研究の大規模かつ長時間にわたる計算は実現出来ません
でした. この場を借りまして, 未来流体情報創造センターの方々に感謝の意を表します.
88
ダウンロード

2成分レナード・ジョンズ コロイド分散系の ガラス転移近傍における