ガラスの非平衡ダイナミクス
山 本 量 一
h 京都大学大学院工学研究科化学工学専攻 615-8510 京都市西京区京都大学桂 [email protected]小 貫 明
h 京都大学大学院理学研究科物理学・宇宙物理学専攻 606-8502 京都市左京区北白川追分町 [email protected] 液体を融点以下に急冷するとアモルファス状態のまま固化してガラスになる．ガラス転移と呼ばれるこの現象は古く から知られているが，転移が起こる本質的な機構については未だに解明されていない．これまで平衡に近い状態での線 形応答に関心が集中したが，最近になってエイジング過程や流動下など平衡から大きく離れた状態でのガラスの非線形 応答に関心が集まり出した．ガラスのように平衡状態さえ満足に分かっていない系の非平衡・非線形の応答など手も足 も出ないと思いがちであるが，最近の計算機シミュレーションは意外なほど単純な性質の存在を明らかにした． 1. はじめに 液体を徐冷すると融点付近で固化して結晶になるが，急 冷すると融点以下でも過冷却の液体相が準安定に存在する． さらに温度を下げると液体的なアモルファス構造が凍結さ れた固体，つまりガラスと呼ばれる状態に到達する．液体 を高圧・高密度状態にしても同様の現象が起こる．複雑な 形状の分子系・サイズの異なる分子の混合系などでガラス 化が起こり易い. ガラス転移と呼ばれる液体・固体変化は， 気液転移や結晶化などの相転移現象と同様に普遍的な現象 であり，様々な立場から研究がなされてきた.1–11) しかし その本質的な機構については明確には分かっておらず物性 物理学の主要な未解決問題の１つである． 図 1 では9)数種の物質の粘性率 η をガラス転移温度 Tgで 規格化された温度の逆数 Tg/T の関数として示す*1．過冷 却状態では原子スケールの粒子の配置換えが熱揺らぎでは 起こり難い. したがって活性化過程を熱揺らぎで乗り越え るための構造緩和時間 τ が温度低下とともに急激に長くな る. τ はストレスの緩和時間でもあり線形粘性率は η = Gτ の関係で表される. ここで G は τ より早い微小変形に対 するずり弾性率で温度に対し緩やかにしか変化しないので， η は τ に比例すると考えて良い. 図では低温になるに従い 粘性率が飛躍的に増大している．ガラス状物質は粘性率の 温度依存性によって”strong” と”fragile” という２つの特徴 的なカテゴリーに分類できる．ほぼ直線で近似できるもの は strong ガラスと呼ばれ，η はアレニウス型の温度依存性 η ∝ exp(E/KBT ) を示す．そこでは液体状態であろうと ガラス転移温度近傍であろうと活性化エネルギー E の値に 大きな変化はない．SiO2(シリカ) や GeO2のように原子の ネットワーク構造がしっかりと保持されている物質の多く がこのカテゴリーに属する． 指向性の強い共有結合を持たず，分子間相互作用が等方 的な低分子からなる物質の多くは fragile ガラスとなる．こ の場合，左端の液体状態からグラフの中ば付近までは通常 のアレニウス則に従い指数関数的な増加を示す．しかしさ らに低温になるとグラフの傾きが急激に大きくなり粘性率 *1Tg は明確な相転移温度を意味しない. 粘性率が 1013ポアズに達し た温度が Tgとして定義されている．また図１は線形，つまり周波数零の 極限でのずり粘性率を図示している． は Vogel-Fulcher 式 η ∝ exp(E/KB(T − T0)) で近似できる ようになる．つまり Tgより低いある温度 T0に向かって粘 性率が強く発散するかのように振る舞う．しかし T0に深い 意味があるのかどうかも含めて，fragile ガラスの粘性率に 急激な増加をもたらす機構はわかっていない*2．このよう な劇的な温度変化は粘性率のみならず種々の動的性質（自 己拡散係数・分子の回転緩和時間など）に見られるが，静 的な性質（圧力・体積・液体構造など）にはごく僅かな変 化しか現れない．静的構造を見る限りマクロからミクロな スケールまで構造に目立った変化が見られず，ガラス状態 は液体状態と大差がないままである．通常の相転移現象と は大きく異なる点である． 僅かな静的構造変化と劇的緩慢化を定量的に結びつける ことができればガラス転移の解明へ向けて大きく前進する． アモルファス構造では長距離秩序が存在しないため，分子 間距離のミクロな液体構造の変化を取り上げる必要がある． この考え方に沿って 1980 年代の前半，臨界現象のダイナ ミクスで成功を収めたモード結合理論がガラス転移に応用 され，ミクロスケールの密度揺らぎの減衰が温度の低下と ともに緩慢化し，ある温度 Tcで凍結するという結果が導 かれた12)．このガラスのモード結合理論では，静的な構造 関数の極大のわずかの増大が動的緩慢化を引き起こしてい る. この理論はその後のガラス転移の研究に大きな影響を 与えたが，理論の予測する Tcは実際のガラス転移温度 Tg よりずっと高い温度であることが後に明らかになった．図 1 で言えば左端から 2/3 ぐらいの範囲はよいが，それより 右側の挙動を説明できない．また数学的近似の物理描像が はっきりしないことや，長波長揺らぎが支配的な臨界現象 に対して有効な近似をガラス転移のようなミクロな現象に 適用することに対し批判がある．このようなガラスのモー ド結合理論については川崎による解説3)がある．液体論に 立脚した取り扱いとは異なる立場の理論的アプローチとし て，トラッピング拡散モデルなどガラス転移を確率過程と して定式化する議論2, 4)，現象論的に複数のオーダーパラ メータを用いる議論13)，エネルギーランドスケープに基づ *2他にもボゾンピークと呼ばれる低エネルギー励起の起源や Kauzmann のパラドックスと呼ばれるアモルファス構造の残余エントロピーに関する 問題など，ガラス転移に付随した未解決問題は多い．く議論4, 9, 11)，スピングラスのレプリカ対称性の破れとの 類似性を論じる議論14)などが提案されている．また液体か ら迫るのでなく逆に格子欠陥が大量に存在する固体に空間 時間に依存した「自由体積」（隙間）を導入しアモルファス 構造の動的性質を理解しようとする試みもなされている.15) 2. ガラス化に伴う動的不均一性と非平衡ダイナミクス 外力や流動がなく平衡に近い場合，粒子は熱揺らぎによっ て運動し拡散する．液体では個々の粒子の拡散は無相関に 近いが，ガラス状態に近づくにつれ高くなるエネルギー障 壁のため個別運動は凍結される．ガラス転移の近傍では強 い相関をもった稀に起こる共同運動だけが生き残り，系の動 的性質を決定する．満員電車内で１人だけが移動しようと すると無理があるが，周りの人々と協調すると動き易いこ とを連想するとよい．ガラス転移に近づくに従って顕在化 する個別→共同という運動性の変化を解析することで，ガ ラス転移の本質に近づけるかもしれない．最近分子動力学 (MD) シミュレーションにより液体からガラス状態に近づ くにつれ粒子の共同運動に起因する動的不均一性が顕著と なり，空間的に拡大することが確認されてきた7, 16–19)．本 稿では我々が行った「動的不均一性」の可視化と定量化の 試みを紹介する．スピングラス20, 21) でも転移点近傍で似 たような長距離相関が出現するという報告がある．前出の ガラスのモード結合理論は均一平均場理論であり長距離の 不均一性は顕わに考慮されていない． ガラス転移に関する従来の物理サイドの主関心は平衡に 近い状態における線形応答にあった．外力が作用して非平 衡状態にあるガラスの応答は工学的な問題*3と直結してい るにもかかわらず，理解はあまり進んでいない. ガラス転 移点近傍において，例えば剪断（シア）流を与えることで 非平衡状態を実現できる.22, 23) 構造緩和時間の逆数 1/τ よ りもシア率 ˙γ を大きくすると， ˙γ の増加とともに粘性率が 顕著に小さくなることが知られている (非ニュートン性）.24) 非平衡状態にあるガラスの研究は，ガラス転移を広い視点 で理解することであり，粉体や泡・液滴などでできた複雑な 散逸系のレオロジーの理解にも寄与する. 図 2 に示される ように10)，温度や密度といった熱力学状態変数を変えると， 赤い平面内での系の液体的 ↔ 固体的変化が問題になる．一 方，Load(剪断流をも含ませた意味の荷重）を与えると，荷 重が小さい時には弾性的（elastic) にあるいは固体的に振る 舞うが，荷重を強くすると塑性的 (plastic) に，ついには液 体的に振る舞う．粉体などでは熱揺らぎは無視できるので， 青い面内での固体的 ↔ 液体的変化（jamming 転移）が問 題になる．ガラスの非平衡ダイナミクスの研究は，これま で別々に研究されてきた赤い面で起こる熱的ガラス転移と 青い面で起こる jamming を統合的に捉えることになる． 本稿では議論しないが，ガラスにおける特徴的な非平衡 現象としてエイジング (aging) 効果が昔から知られている． 高温で平衡状態にある系をガラス転移温度以下に急冷した *3例えば無結晶合金のアモルファス金属は破壊され難く強度が高い. 際に出現する非平衡状態では，時々刻々と遅くなる特殊な 緩和現象が観測される． 3. シミュレーションによる解析 MD シミュレーションとは，分子や原子など物質を構成 する粒子間の相互作用を近似的なポテンシャル関数で置き 換え，粒子の運動方程式をコンピューターで数値的に解く ことによって注目する系の静的・動的な性質を調べる手法 である25)．MD シミュレーションでは物質のミクロな時系 列情報が得られるのですべての物理量が計算可能であると いう利点がある．反面，近似的なポテンシャルを用いる以 上選んだモデル系が現実系をどの程度正確に再現している のかという問題が常にある. しかしガラス転移などのよう に物質の細かい個性には依存しないがそれでいて直感を超 えた普遍的現象を研究するには非常に有効である．今回は モデルとして 2 次元 (2D) 及び 3 次元 (3D) でのソフトコ ア粒子の 2 成分混合物を用いた*4．粒子間ポテンシャルは α, β(∈ 1, 2) 成分の粒子の大きさ σα, σβを用いて vαβ(r) = ² µ σαβ r ¶12 , σαβ= 1 2(σα+ σβ), (1) と粒子間距離 r のべき乗で表される. 粒子数 N1 = N2 = 5000，サイズ比 σ2/σ1 = 1.4 (2D) 又は 1.2 (3D)，質量比 m2/m1= 2，全数密度 n1+ n2= ρ = 0.8/σd1 とし（d は次 元），通常液体からガラスに近い過冷却状態まで (0.337 ≤ T ≤ 2.45 (2D)， 0.234 ≤ T ≤ 0.772 (3D)) 温度を変えて シミュレーションを行った．典型的な粒子配置を示した図 3 から推察されるようにサイズ比は高密度状態での粒子配 置の乱れ（frustration）の程度を決める. 本稿では長さと 時間の単位として σ1と τ0= (m1σ21/²)1/2を用いる．温度 の単位は ²/kBである. 静止状態でのシミュレーションに加 え，平均流がシア流 ˙γyex(ex:x 方向の単位ベクトル) であ る流動状態でもシミュレーションを行った． 3.1 動的不均一性の可視化と定量化：ボンド切断 局所的な粒子の再配置と拡散を定量化するために，初期 時刻 t0で隣接していた粒子対が時間の経過とともに離れて いくイベントを追跡する．時刻 t0で粒子対 j, k の位置が次 式を満たせば隣接していると判断し，粒子対に仮想的なボ ンドを与える． rjk(t0) = |rj(t0) − rk(t0)| ≤ A1σαβ, (2) ここで rj(t) は粒子 j の時刻 t における位置を表す. 2 次元 での典型的な粒子配置と t0におけるボンドの様子を図 3 に 示す．この粒子対 j, k が時間 ∆t の後に次の式 (3) の条件を 満たせば，両者は ∆t の間に離れた，つまりボンドは消滅 したと判断する． rjk(t0+ ∆t) > A2σαβ (3) 今回は隣接・解離の判定に A1= 1.1，A2= 1.6 (2D)，A1= A2 = 1.5 (3D) を用いたが，これらの数値を多少変えても *4内部構造のない単純な１成分系では容易に結晶化が起こってしまう． ２成分系でサイズ比を大きくし過ぎると相分離する.
結果は大きく変わらない．時間間隔 ∆t が大きくなるとと もに生き残るボンドの数はほぼ指数関数的に減少し，初期 の値の e−1になる ∆t をボンドの平均寿命 τb= τb(T, ˙γ)（T と ˙γ の関数）と定義する. 前出の熱的緩和時間 τ は τb(T, 0) と同一視してよい. それではガラスの動的不均一性の可視化に移ろう．図 4 は時間間隔 [t0, t0+ 0.05τb] の間に消滅したボンドの空間分 布を (a) 通常液体と (b) ガラスに近い過冷却液体について 示したものである．どちらもシアはかかっていない．各ドッ トは消滅したボンドの中点 Rjk= (rj+ rk)/2 が時刻 t0に 存在した場所であり，∆t = 0.05τbの間にその近傍で粒子 配置の変化（構造緩和）が起こっている．通常液体 (a) で は配置変化はほぼ均一であるが，過冷却液体 (b) ではドッ トの分布に強い相関が生じ，特徴的な長さ ξ を持った構造 緩和の起こり易いクラスター状の領域が不均一に分布して いる．これがガラス化に伴い顕在化する動的不均一性の実 態である． 図示した動的不均一性を定量的に解析するために切れた ボンドの中間位置 Rjk に対する構造因子を次式で求める． Sb(q) = ¿¯¯ ¯ ¯ X broken bonds exp(iq · Rjk) ¯ ¯ ¯ ¯ 2À (4) 図５に示すように計算された Sb(q) は磁性体の臨界点近傍 での秩序変数の空間相関を表す Ornstein-Zernike 型23)に Sb(q) = Sb(0)/(1 + ξ2q2). (5) 非常によく合致する. このため不均一性の相関長 ξ = ξ(T, ˙γ) が求められる. さらなる臨界現象との類似点として，波数 q が大きくなると Sb(q) の値が温度に余り依存しなくなるこ と，Sb(0) ∝ ξ2であることも挙げられる. 動的不均一性は，実験的にも共焦点顕微鏡を用いたコロ イドガラスの直接観察によって確認されている26, 27)．また， 種々の実験結果を説明する際にしばしば使われたアイディ アに共同運動領域 (CRR) というものがある5, 6)．実験技術 の進歩によってその実態が明らかになりつつあり28)，本稿 の動的不均一性と共通点が多い． 3.2 動的不均一性と粒子の自己拡散 動的不均一性が顕著に影響を与える例として，自己拡散 がある．通常液体状態では粒子（直径 σ）の自己拡散係数 D と粘性率 η との間に Stokes-Einstein 則 D = kBT /2πση (6) がよい精度で成立している．しかし多くのガラス形成物質 についてガラス状態に近づくにつれ D は上記則よりも何桁 も大きいことが観測されている．即ち D は 1/η ∝ 1/ταほ どには減少しない. 我々は 3 次元の過冷却液体中での 1 粒 子運動を調べ，そのメカニズムについて調べた29)． 小さい粒子（直径 σ1の粒子 1) の散乱関数を導入する. Fs(q, t) = N1−1hX j eiq·∆rj(t)i (7) 粒子数 N1で割ることで Fs(q, 0) = 1 と規格化する. 変位ベ クトル ∆rj(t) はシア流下で次のように定義する. ∆rj(t) = rj(t) − rj(0) − ˙γ Z t 0 dt0yj(t0)e x (8) ここでは平均流による移送効果が第三項により除去されて いる. 構造緩和時間 ταは Fs(q, t = τα) = e−1 (9) で決める（ただし q = 2π/σ1）. するとシア下の非ニュー トン領域までをも含めた広い T と ˙γ の範囲で ταと τbは τα∼= 0.1τbで，η と ταは η = (2πkBT /σ13)ταでほぼ結びつ くことを確認した．後者の関係は実験的にも確かめられて いる30)．これを用いると Stokes-Einstein 則 (6) は Dτα= σ2 1/(2π)2 となり Dταが T と ˙γ によらず一定値をとること になる．しかし計算結果では図 6 に示されるように，ガラ ス化に近づくに従い Dταの値は増大する．この原因を考 える．時間領域 t ≥ 0.1ταにおいては Einstein の式 6Dt = h(∆r(t))2i が成立しているため29)，van Hove 自己相関関 数 Gs(r, t) ≡ h PN1 j=1δ(∆rj(t) − r)i/N1を用いて Dτα=1 6h(∆r(τα)) 2i = 1 6 Z ∞ 0 dr 4πr4Gs(r, τα) (10) を得る．図 7 では過冷却状態で 4πr4Gs(r, τα) を図示して いる．時間 t = ταにおいて高温 (T = 0.473) ではガウス分 布になっているのに対し，低温 (T = 0.267) では r > 1 の 部分がガウス分布より著しく増加しており Dταの増大をも たらしている．過冷却状態では動き易い粒子と動き難い粒 子の差がガウス分布の予測よりずっと大きい. ガラス転移 近傍における自己拡散運動の非ガウス的挙動は実験的にも 確認された31)． 図 8 は (a) 通常液体と (b) ガラスに近い過冷却液体につ いて，時間間隔 ∆t = 0.1τα に個々の粒子が実際に移動し た変位を三角錐で表示したものである．通常液体では粒子 の変位は分散が小さく一様であるが，過冷却状態では分散 が大きく著しく不均一である．複数の三角錐が連結してス トリングを形成し，さらにそれが集まり凝集している様に 見える. この変位の不均一空間分布は前節のボンド切断の 動的不均一性の空間分布とぴったり一致することが確認さ れており，D の Stokes-Einstein 則からのずれは動的不均一 性に由来すると結論できる．なお不均一構造は時間空間で 揺らいでおり，時間間隔が動的不均一性の寿命 (∼ 3τα) よ り長くなると徐々に消滅する． 1 粒子運動についてさらに考察する．液体状態での粒子 運動は単一の時定数を持つブラウン運動でよく近似できる． 過冷却状態では時定数に空間分布が出現すると考えると， Gs(r, t) は局所的な拡散係数 D(x, t) を持つガウス分布 Gs(x, r, t) = [4πD(x, t)t]−3/2exp[−r2/4D(x, t)t] (11) の重ね合わせで表現できるだろう．図 9 にスケールされた 分布関数√6Dt4πr2Gs(r, t) を r∗ = r/√6Dt に対してプ
ロットした．カーブの下の面積は１に規格化されている. ス ケールされた分布関数は，動的不均一性の寿命より短い時 間 t . 3ταでは r & 1 (r∗& (6Dt)−1/2) の領域でマスター カーブに漸近している．これは予測 (11) と整合している. しかし r < 1 の領域では曲線は互いに大きくずれる．つま り粒子間距離より大きなスケールでは局所拡散係数の描像 はよいが，それより小さなスケールでは適用できない．動 的不均一性の寿命より長い時間 t & 10τα（図 9 挿入図）で は不均一性がならされ，分布関数はガウス分布に漸近する． 3.3 ガラス化近傍の非平衡ダイナミクス シア下の非平衡状態を考える. 我々のシミュレーションで はボンドの減衰率（ボンドの平均寿命 τb∼= 10ταの逆数）は 1/τb(T, ˙γ) = 1/τb(T, 0) + Ab˙γ (12) と，単純に熱的寄与とシアによる寄与の和で書けた. Abは 定数である (今回用いたモデルでは 0.57 (2D)，0.80 (3D)). (9) 式以下に述べたように粘性率は η ∝ ταであるから η(T, ˙γ) = η(T, 0)/(1 + τη˙γ) (13) ここで τη= Abτb(T, 0) ∝ τα(T, 0) である. ガラス化近傍で は τb(T, 0) は極めて長く非ニュートン領域 1/τb(T, 0) < ˙γ ¿ 1/τ0は非常に広い (τ0はミクロな時間). そこでは τα∼ 1/ ˙γ であり平均ストレス η ˙γ は一定値を取る. これらの関係は 実験でも確認されている24, 32)．この振る舞いは次のように も説明できる. 粒子配置変化が時定数 1/ ˙γ で起こり ² 程 度のエネルギーを放出するなら，単位体積当たり熱発生率 は ˙Q ∼ n² ˙γ となる（n は平均密度）．一方流体力学的関係 ˙ Q = η ˙γ2から η ∼ n²/ ˙γ が得られる. シア下で動的不均一性はどうなるか？図 4(c) は (b) と同 じ温度でシアが与えられているが，クラスターのサイズが 小さくなり分布が均一化していることがわかる．3 次元で も同様の結果が得られている．図 10 に示すが広い T と ˙γ の範囲で不均一性の相関長 ξ とボンドの寿命 τb の間に簡単 な動的スケーリング則が近似的に見い出された． τb∝ ξz (14) ここで 2 次元で z = 4，3 次元で z = 2 である．非ニュート ン領域では ξ ∝ ˙γ−1/zである. 低次元の方が準安定配置の 解消に時間が費されるであろうことから z はより大きいは ずであるが，z の具体値の導出はまだない. 3.4 シア下での等方性と有効温度 1 つの粒子の周りの局所構造を考える．熱揺らぎなどで 粒子の配置換えが発生する時間間隔は ταであるが，配置換 え自体はミクロなスケールの出来事であり完了するまでに ミクロな時間 τ0しかかからない. そのため配置換えがシア で引き起こされるにしても，シア率 ˙γ が 1/τ0よりはるかに 小さい限り，局所構造の異方性は小さい. 実際に非線形領 域でも図３のような粒子配置や図４のような切断ボンドの スナップショットを見るだけではシアの存在は判然としな い. 例えば同時刻対相関関数 gαβ(r−r0) = hnα(r)nβ(r0)i の ピークの高さの異方性は高々5% 程度しかない*5. 式 (7) の １粒子運動に対応する散乱関数 Fs(q, t) においては式 (8) に より対流効果が排除されており，波数べクトル q の方向依 存性はほとんどない. 臨界溶液などの遅く緩和する系に流動 を与えると散乱強度は容易に非等方になるのに対し22, 23)， 低分子ガラスの等方性はユニークである． モード結合理論にシアを取り入れる試みがなされ，シアに よって緩和が促進されるという我々のシミュレーション結果 と一致する結論が得られている.33, 34) 成分 1 の密度時間相 関関数 hn1(r, t)n1(r0, t0)i を考えよう. 定常状態でこの関数 は G(R, t−t0) とおける. ただし R = r−r0− ˙γ(t−t0)ye xで exは流動方向の単位ベクトル).23) 動的構造関数 (G(R, t) のフーリエ変換） F (q, t) = hX jk
eiq·(rk(0)−rj(t))+iqx˙γtyj(t)i (15) の数値計算がなされたが，図 10 に示すようにこの量は波数 べクトル q の方向にほとんどよらない35)．即ちシアによる 構造緩和の促進がほとんど等方的に起こる． これらのの驚くほど簡明なシアの効果はシア下での有効 温度 Teff の存在を示唆する. 即ち τb(Teff, 0) = τb(T, ˙γ) と 置くことによりシア下の定常状態をより高い温度 Teffの静 止状態にマップできると考える. 簡単のためアレニウス型 の温度依存性 τb(T, 0) = C exp(E/kBT ) を仮定すると， E/Teff = E/T − kBln · 1 + AbC ˙γ exp(E/kBT ) ¸ (16) となり，シアは臨界現象における磁場のような役割を果た す.18) 最近注目されているトピックの一つに非平衡状態におけ る揺動散逸定理 (FDT) の破れの問題がある．エイジングの 過程にあるスピングラスについて議論が始まり38, 39)，ガラ スの分野でも考察された．ガラスではエイジング実験の待 ち時間 twの代わりにシアによっても非平衡度を制御でき る．非平衡では通常の FDT に出てくる温度 T を有効温度 Teffに置き換えて R(t) = − 1 kBTeff dC(t) dt (17) とする．ここで R(t) は任意の物理量の応答関数，C(t) は 対応する時間相関関数である．粒子密度の揺らぎに関する 応答関数と相関関数の関係をミュレーションにより調べ実 際に Teffを決めた例がある36)．シア下の系にトレーサー粒 子を入れ，その質量を変えて粒子の運動エネルギーを測定 したものもある36)．トレーサー質量が軽い場合は早い運動， 重い場合は遅い運動をとらえることになるが，驚くべきこ とに質量を大きくしていくと運動エネルギー，つまりその トレーサーの感じる系の温度が上昇する．しかもその上昇 した温度が先の FDT で出てきた Teff に一致するというの *5平均ストレス η ˙γ は g αβ(r) を使った空間積分で書けるので23)，その 異方性は重要な意味がある.
である．その後他の物理量についても同様の手法で有効温 度が決定されたが，それらの間の一致は良いようである37)． 5. まとめと将来の問題 MD シミュレーションによってガラス転移近傍における ボンド切断事象から動的不均一性を可視化した．長い緩和 時間のスケールで観察すると，低温になるに従い粒子の配 置換えが活発な領域と不活発な領域とのコントラストが際 立ってくる．それは臨界揺らぎと酷似している．この動的 不均一性は自己拡散運動に重大な影響を与えている．しか しガラスのモード結合理論には長距離相関が現れないこと に注意しよう．ガラス転移はミクロな機構で起こるといえ る．ガラスの動的不均一性は磁性体における臨界揺らぎと のような第一義的意義を持つわけではないが，それだから こそ起源・役割は依然良く分からない． シア流を与えることでガラス研究の内実が飛躍的に深ま り広がることも示した．アモルファス合金の変形機構や粉 体などのレオロジーと共通の問題に直面する．シア下でも 動的不均一性の相関長 ξ と寿命 τbの間には臨界現象に類 似した簡単な動的スケーリング則が成立していることがわ かった．また流動状態にあるガラスのダイナミクスは非常 に等方的であり，より高い有効温度 Teff で静止状態にある 系でのダイナミクスとほとんど同じように見える． あるガラス形成物質においては実に 20 − 200nm にも及 ぶ密度の不均一性が観測されている.40) 類似の不均一性が ゲルにおいて観測されている.23) ゲル生成時の網目構造 の不均一性が弾性相互作用により，温度変化の後に密度変 動を引き起こす. ガラスでも似たようなシナリオが想定で きる. また不均一のゲルに一軸変形を与えると abnormal butterfly 呼ばれる異方的散乱が観測される. ガラスでのこ のような実験は興味深い． 高分子ではガラス化に伴い粒子間構造緩和時間 ταと鎖の 緩和時間 τR = N2ταがともに長くなる.41) N は鎖上の ビーズの数であり絡み合いはないとした. 散乱強度は非常 に小さいシア ˙γ < 1/τRのもとで異方的になる．シア下で 散乱と誘電測定を組み合わせた実験が興味深い． 本研究の一部は文部科学省科学研究費の補助の下に行わ れた．シミュレーションは京都大学化学研究所スーパーコ ンピュータラボラトリ，及び東京大学医科学研究所ヒトゲ ノム解析センターの並列計算機を用いて行なわれた． 参考文献 1) 江上毅: 日本物理学会誌 39 (1984) 505. 2) 樋渡保秋, 宮川博夫, 小田垣孝: 日本物理学会誌 46 (1991) 90. 3) 川崎恭治: 日本物理学会誌 48 (1993) 869; 「非平衡と相転移」(朝 倉書店, 2000). 4) 小田垣孝: 固体物理 33 (1998) 489. 5) 金谷利治: 高分子 52 (2003) 773. 6) 山室修: 高圧力の化学と技術 9 (1999) 117. 7) 村中正, 樋渡保秋: 熱物性 13 (1999) 92. 8) 山本量一, 小貫明: 高圧力の化学と技術 9 (1999) 134. 9) P.G. Debenedetti, F.H. Stillinger: Nature 410 (2001) 259. 10) A.J. Liu, S.R. Nagel: Nature 396 (1998) 21.
11) S. Sastry: Nature 409 (2001) 164.
12) W. G¨otze: in Liquids, Freezing, and the Glass Transition, ed. J.P. Hansen, D. Levesque, and J. Zinn-Justin, Elsevier, Ams-terdam, 1991, pp.289-503.
13) H. Tanaka: J. Chem. Phys. 111 (1999) 3163.
14) B. Coluzzi, G. Parisi, P. Verrocchio: Phys. Rev. Lett. 84 (2000) 306; R. Di Leonardo, L. Angelani, G. Parisi, G. Ruocco: Phys. Rev. Lett. 84 (2000) 6054.
15) A. Onuki, A. Furukawa, A. Minami: in the proceedings of Stat-phys 22, Bangalore, 4-9 July, 2004.
16) M.M. Hurley, P. Harrowell: Phys. Rev. E 52, (1995) 1694: D.N. Perera: J. Phys.: Condens. Matter 10, (1998) 10115.
17) R. Yamamoto, A. Onuki: J. Phys. Soc. Jpn. 66, (1997) 2545. 18) R. Yamamoto, A. Onuki: Phys. Rev. E 58, (1998) 3515. 19) W. Kob et al.: Phys. Rev. Lett. 79, (1997) 2827; C. Donati et
al.: Phys. Rev. Lett. 80, (1998) 2338.
20) 高山一: 「スピングラス」(丸善, 1991).
21) 西森秀稔: 「スピングラス理論と情報統計力学」(岩波書店, 1999). 22) 小貫明: 日本物理学会誌 50 (1995) 20.
23) A. Onuki: “Phase Transition Dynamics” (Cambridge University Press, 2002).
24) J.H. Simmons, R.K. Mohr, C.J. Montrose: J. Appl. Phys. 53, (1982) 4075. J.H. Simmons, R. Ochoa, K.D. Simmons, J.J. Mills: J. Non-Cryst. Solids 105, (1988) 313.
25) 上田顕: 「分子シミュレーション」 (裳華房, 2003) 26) W.K. Kegel, A. van Blaaderen: Science 287 (2000) 290. 27) E.R. Weeks, J.C. Crocker, A.C. Levitt, A. Schofield, D.A.
Weitz: Science 287 (2000) 627.
28) M.D. Ediger: Annu. Rev. Phys. Chem. 51 (2000) 99. 29) R. Yamamoto, A. Onuki: Phys. Rev. Lett. 81, (1998) 4915. 30) F. Mezei, W. Knaak, B. Farago: Phys. Rev. Lett. 58, (1987)
571; D. Richter, R. Frick, B. Farago: Phys. Rev. Lett. 61, (1988) 2465.
31) R. Zorn: Phys. Rev. B 55, (1997) 6249; T. Kanaya, I. Tsukushi, K. Kaji: Supplement to Prog. Theor. Phys., 126, (1997) 133. 32) R. Di Leonardo, F. Ianni, G. Ruocco: Phys. Rev. E 71 (2005)
011505.
33) M. Fuchs, M. E. Cates: Phys. Rev. Lett. 89, 248304 (2002). 34) K. Miyazaki, D.R. Reichman: Phys. Rev. E 66 (2002) 050501. 35) K. Miyazaki, D.R. Reichman, R. Yamamoto: Phys. Rev. E 70
(2004) 011501.
36) L. Berthier, J.-L. Barrat: J. Chem. Phys. 116 (2002) 6228. 37) I.K. Ono, C.S. O’Hern, D.J. Durian, S.A. Langer, A. Liu, S.R.
Nagel: Phys. Rev. Lett. 89 (2002) 095703; C.S. O’Hern, A. Liu, S.R. Nagel: Phys. Rev. Lett. 93 (2004) 165702.
38) 川村光: 日本物理学会誌, 59 (2004) 9.
39) 佐々木志剛，根本幸児: 日本物理学会誌 58 (2003) 580. 40) E.W. Fischer: Physica A 201, 183 (1993).
図 1: 種々の物質の粘性率 η をガラス転移温度 Tgで規格化 した温度の逆数 Tg/T の関数としてプロットしたもの．温 度依存性により”strong” と”fragile” という２つの特徴的な カテゴリーに分類されている．文献 9 より． 図 2: 一般化されたガラス/jamming 転移の概念的な相図． 面で囲まれた内側の原点近傍の領域がガラス/jamming 状 態を表す．x，y，z 軸はそれぞれ荷重，密度の逆数，温度 を表し，各々原点から離れるに従ってガラス/jamming 状 態から遠ざかる．文献 10 より． 図 3: 2 次元系における典型的な粒子配置．濃い灰色の円は 小さい粒子，薄い灰色の円は大きい粒子，粒子の中心を結 んだ直線は隣接ボンドを表す．
(a)
(b)
(c)
図 4: (2), (3) 式で定義される ∆t = 0.05τbの間に消滅した 隣接ボンドの空間分布．(a) は T = 2.54（液体状態），(b) は T = 0.337（過冷却状態）でそれぞれシアがない場合．低 温になりガラス化するに従って，構造緩和が不均一になる ことがわかる．(c) は T = 0.337 で ˙γ = 0.25 × 10−2のシア がかかっている場合．シアの効果で均一性が回復されるこ とがわかる．矢印は (5) 式にフィットして得られた相関長 ξ の長さを表す．Tξ 6E T 6 E 7 図 5: 各温度・各シア率の下で得られた消滅したボンドの 構造因子（(4) 式で定義）のスケーリングプロット（3 次元 系）．実線は Ornstein-Zernike 型の関数 1/(1 + (qξ)2) を表 すが，シミュレーションデータとの一致は大変よい． 10–1 100 101 102 103 104 105 10–2 10–1 100 τα D τα T = 0.234 0.267 0.306 0.352 0.473 0.772 図 6: 各温度で得られた Dταを ταの関数としてプロット． 直線は Einstein-Stokes 則を表すが，ταが大きいつまり低 温側ではシミュレーションデータと合わなくなる． 0 1 2 3 4 0.1 0.2 0.3 0.4 0.5 4 π r 4 G s (r, τα ) r T = 0.234 T = 0.473 Gaussian 図 7: 4πr4Gs(r, τα) を r の関数としてプロットしたもの． 波線は T = 0.473（液体状態），実線は T = 0.267（過冷却 状態）のシミュレーション結果，点線は単純なブラウン運 動の結果である．それぞれ曲線と x 軸で囲まれた部分の面 積が 6Dταを与える． 図 8: 時間間隔 ∆t = 0.1ταに個々の粒子が実際に移動した 変位を三角錐で表示．(a) は T = 0.473（液体状態），(b) は T = 0.267（過冷却状態）である．過冷却状態では三角 錐の長さの分散が大きくなり，空間的にも不均一になる． 1 2 3 4 0 0.5 1 1.5 1 2 0 0.5 1 t=10τα t=100τα (6Dt) 1/2 4 π r 2 G s (r,t) r(6Dt)–1/2 t=0.5τα t=3τα Gaussian 図 9: 動的不均一性の寿命よりも短い時間領域 0.5τα< t < 3ταでの van Hove 自己相関関数 Gs(r, t) のスケーリングプ ロット．テールの部分に Gauss 型とは異なる漸近形の存在 が確認できる．挿入図は動的不均一性の寿命より長い時間 領域でのプロット．長時間の後に Gauss 型に漸近していく ことが分かる．
100 101 102 100 101 102 103 104 105 L τb ( T, γ ) ξ( T, γ ) . . T = 0.337 1.43 0.85 0.526 2.54 (a) 100 101 102 100 101 102 103 104 105 L τb ( T, γ ) ξ( T, γ ) . . T = 0.234 0.352 0.306 0.267 0.473 (b) 0.772 図 10: 各温度・各シア率の下で得られた τbと ξ に見いだ されるユニバーサルな関係．(a) は 2 次元，(b) は 3 次元の 結果である．実線は (14) 式で表させる動的スケーリング則 であるが，2 次元では z = 4，3 次元では z = 2 とするとシ ミュレーションデータとよい一致が見られる． 10−2 10−1 100 101 102 103 104 105 0 0.5 1 Fρρ (q ,t )/S ρρ (q ) t 0 γ = 10−1 10−4 (a) (b) (c) (d) 10−3 10−2 . qm=5.8 図 11: (15) 式で定義された密度相関関数 F (q, t) のプロッ ト．2 次元系 T = 0.526 の過冷却状態の結果である．|q| = qm= 5.8．(a)，(b)，(c)，(d) は異なる方向を表しており， それぞれ q={qx, qy} = qm{1, 0}， √ 2 2 qm{1, 1}，qm{0, 1}， √ 2 2 qm{−1, 1} である．
