7
非平衡系での構造形成
7.1
非保存 vs. 保存
前章で導出した TDGL 方程式では，各相の物質量（粒子数）が保存されない．そのため， 相転移温度以下では，多くの場合，ほぼ完全にどちらかの相だけになってしまう．これは，第 4 章で紹介した Metropolis Monte Carlo 法における，原子種の転換 に 対応する． 磁性転移（スピンの向きの揃い方が変わる）のように 非保存の物理量の相転移を調べるに はよいモデル であるが，合金系の相分離のような場合には少し不自然に見える． 保存系を調べたい場合には，連続の式equation of continuity を出発点にする必要がある： ∂s(⃗r, t) ∂t + ⃗∇ ⃗J (⃗r, t) = 0 (79) ここで， ⃗J は 秩序変数 s の流束flux である．前節の式 (67) に代わって，J が外力に線形⃗ に依存することを仮定する： ⃗ J (⃗r, t) =−k⃗∇δF δs (80) 式 (79) と組み合わせると ∂s(⃗r, t) ∂t = k∆ δF δs (81) 自由エネルギーモデルは前節と同じもの，式 (59)，を仮定すると，前章と同様の無次元化 定係数 k は無次元化により消去して しまった．あるいは１としたと考えて もよい． により次の保存系 TDGL 方程式が得られる： ∂s ∂t =−∆ [ −T∗s− s3+ ξ∆s] (82) 図 7–14 は，リスト 7.1 のプログラムで求めた１次元系の計算結果である．非保存系 (図 6–12) と比較してみよう．時間発展とともにやがてどちらかの相に統一されてしまう非保 存系の場合とは異なり，特徴的な空間周期構造spatial periodic structure が発達するのがわ かる．また，緩和のしかたは非保存の場合と比べてゆっくりとしている．従って，空間構 造が発達する様子を調べるためには計算ステップ数を増やす必要がある． なお，保存系の TDGL 方程式では，空間の４階微分演算子 ∆2≡ ∂4 ∂x4+ ∂4 ∂y4+· · · が現れ る．このとき，陽解法では r≡ ∆t (∆x)4 について，比 r :−4r : 1 + 6r : −4r : r の空間平均を とることになり，一般的には２階微分の場合に比べて小さな ∆t で刻む必要があるだろう． リスト 7.2 と図 7–15 は２次元系のプログラムと計算例である．#include <stdio.h> #include <stdlib.h> #define NUM_CELL 100 #define NUM_STEP 80000 #define NUM_SAVE 200 double state[NUM_CELL]; void initial( ) { int icell; for (icell=0;icell<NUM_CELL;icell++) { state[icell]=0.1*(rand( )/(double)RAND_MAX-0.5); } }
void data_out(int iout) { int icell; char filen[100]; FILE *filep; sprintf(filen,"T-1.0/%5.5d.dat",iout); printf("%s\n",filen); filep=fopen(filen,"w"); for (icell=0;icell<NUM_CELL;icell++) { fprintf(filep,"%7.4f\n",state[icell]); } fclose(filep); } int main() { int istep,icell,icell0,icell1,iout; double temp,st,st0,st1; double dstate[NUM_CELL]; double xi=1.0; double dtime=0.05; scanf("%lf",&temp); initial(); istep=0; data_out(istep); for (istep=1;istep<=NUM_STEP;istep++) { for (icell=0;icell<NUM_CELL;icell++) {
icell0=icell-1; if (icell0==-1) icell0=NUM_CELL-1; icell1=icell+1; if (icell1==NUM_CELL) icell1=0; st=state[icell];
dstate[icell] = -temp*st - st*st*st + xi*(state[icell0]+state[icell1]-2*st); }
for (icell=0;icell<NUM_CELL;icell++) {
icell0=icell-1; if (icell0==-1) icell0=NUM_CELL-1; icell1=icell+1; if (icell1==NUM_CELL) icell1=0;
state[icell] -= dtime*(dstate[icell0]+dstate[icell1]-2*dstate[icell]); } if (istep%NUM_SAVE==0) data_out(istep); } return 0; } List 7.1: １次元保存系TDGLモデルの数値計算プログラム例 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0 20 40 60 80 100 O rd e r Pa ra me te r s Spatial Coordinate T= -0.1 0 200 2000 10000 20000 30000 50000 80000 sqrt(-T) -1 -0.5 0 0.5 1 0 20 40 60 80 100 O rd e r Pa ra me te r s Spatial Coordinate T= -1.0 0 200 2000 10000 20000 21000 25000 30000 sqrt(-T) 図 7–14: １次元保存系 TDGL 数値計算結果の例
#include <stdio.h> #include <stdlib.h> #include <math.h> #define NUM_CELL 100 #define NUM_STEP 200000 #define NUM_SAVE 1000 #define dtime 0.02 // 非保存系よりも小さくとる必要がある double state[NUM_CELL][NUM_CELL]; FILE *fplot; void initial() //---{ int ic,jc; for (ic=0;ic<NUM_CELL;ic++) for (jc=0;jc<NUM_CELL;jc++) state[ic][jc]=1.0*(rand()/(double)RAND_MAX-0.5); }
void data_out(int istep) //---{
int ic,jc; char cfile[100]; static int flag=0; FILE *filep; sprintf(cfile,"%6.6d.dat",istep); filep=fopen(cfile,"w"); for (ic=0;ic<NUM_CELL;ic++) { for (jc=0;jc<NUM_CELL;jc++) { fprintf(filep," %7.4f",state[ic][jc]); } fprintf(filep,"\n"); } fclose(filep); // gnuplot コマンドファイル出力 fprintf(fplot,"splot \’%s\’ matrix\n",cfile); if (flag==0) { fprintf(fplot,"pause -1\n"); flag=1; } else { fprintf(fplot,"pause 0.1\n"); } } int main() //---{ int istep,ic,jc,ic0,ic1,jc0,jc1; double temp,st,grad;
static double dstate[NUM_CELL][NUM_CELL]; static double ddstate[NUM_CELL][NUM_CELL];
double xi=1.0; // 表面張力係数 ξ は 1.0 に固定してある
printf("Input T: "); scanf("%lf",&temp);
fplot=fopen("tdgl.plt","w"); // gnuplot (Ver 4.6) コマンドファイル作成 fprintf(fplot,"set term windows font \"Arial,20\"\n");
fprintf(fplot,"set size square\n"); fprintf(fplot,"unset surface\n"); fprintf(fplot,"set style data pm3d\n"); fprintf(fplot,"set view map\n");
fprintf(fplot,"set pm3d implicit at b\n");
fprintf(fplot,"set palette rgbformulae 33, 13, 10\n"); fprintf(fplot,"set key at %3d,%3d\n",NUM_CELL,NUM_CELL+10);
fprintf(fplot,"set cbrange [%6.1f:%6.1f] noreverse \n",-1.2*sqrt(fabs(temp)),1.2*sqrt(fabs(temp))); fprintf(fplot,"set size square\n");
initial(); istep=0; data_out(istep); for (istep=1;istep<=NUM_STEP;istep++) { for (ic=0;ic<NUM_CELL;ic++) { for (jc=0;jc<NUM_CELL;jc++) { st=state[ic][jc];
ic0=ic-1; if (ic0==-1) ic0=NUM_CELL-1; ic1=ic+1; if (ic1==NUM_CELL) ic1=0; grad =state[ic0][jc]+state[ic1][jc]-2*st;
jc0=jc-1; if (jc0==-1) jc0=NUM_CELL-1; jc1=jc+1; if (jc1==NUM_CELL) jc1=0; grad+=state[ic][jc0]+state[ic][jc1]-2*st; dstate[ic][jc]= -temp*st - st*st*st + xi*grad; }}
for (ic=0;ic<NUM_CELL;ic++) { for (jc=0;jc<NUM_CELL;jc++) {
st=dstate[ic][jc];
ic0=ic-1; if (ic0==-1) ic0=NUM_CELL-1; ic1=ic+1; if (ic1==NUM_CELL) ic1=0; grad =dstate[ic0][jc]+dstate[ic1][jc]-2*st; jc0=jc-1; if (jc0==-1) jc0=NUM_CELL-1; jc1=jc+1; if (jc1==NUM_CELL) jc1=0; grad+=dstate[ic][jc0]+dstate[ic][jc1]-2*st; ddstate[ic][jc]=grad; }} for (ic=0;ic<NUM_CELL;ic++) { for (jc=0;jc<NUM_CELL;jc++) { state[ic][jc] -= dtime*ddstate[ic][jc]; }} if (istep%NUM_SAVE==0) { printf("%7d step\n",istep); data_out(istep); } } fprintf(fplot,"set out\n"); fclose(fplot); return 0; }
List 7.2: ２次元保存系 TDGL モデルの数値計算プログラム例．やはり gnuplot 4.6 (or later) 用のコマンドファイルを生成する． ’000000.dat’ matrix 0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100 -1.5 -1 -0.5 0 0.5 1 1.5 ’005000.dat’ matrix 0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100 -1.5 -1 -0.5 0 0.5 1 1.5 ’010000.dat’ matrix 0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100 -1.5 -1 -0.5 0 0.5 1 1.5 ’050000.dat’ matrix 0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100 -1.5 -1 -0.5 0 0.5 1 1.5 ’100000.dat’ matrix 0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100 -1.5 -1 -0.5 0 0.5 1 1.5 ’200000.dat’ matrix 0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100 -1.5 -1 -0.5 0 0.5 1 1.5 図 7–15: ２次元保存系 TDGL 方程式の数値計算（List 7.2）の可視化例：T∗=−1.0 の場合．図 6–13 と比べてみよう．
7.2
空間構造の解析
非平衡状態からの緩和が生み出す空間構造やその時間変化を定量的に調べる標準的な方 法の１つが，空間自己相関関数 spatial autocorrelation function G(⃗r, t) を求めることであ る．この関数は，各時刻における，秩序変数 s(⃗r, t) の空間的相関を調べるものであり，次 式で定義される： G(⃗r, t)≡ ⟨s(⃗r0+ ⃗r, t)· s(⃗r0, t)⟩0 ⟨s(⃗r0, t)· s(⃗r0, t)⟩0 (83) ここで，⟨. . .⟩0は，いろいろな場所 ⃗r0についての統計平均をとることを表す．２次元以上 で，等方的な系であれば，当然あらゆる方向についても平均をとることができるので，G は距離のみの関数 G(r, t) となる．なお分母は単なる規格化因子 [r = 0 で G(0, t) = 1 とな るように] である． プログラム例（別途求めた TDGL 方程式の数値解を解析するツール）を List 7.3 に，ま た結果の例を図 7–16 に示した． (補足) 自己相関関数 autocorrelation function 前期の「熱物理工学」で次の形の自己相関関数を紹介したが，受講者の方は憶えている だろうか？ V AF (t) =⟨⃗v(t0)· ⃗v(t0+ t)⟩t 0 (84)
これは，速度自己相関関数 velocity autocorrelation function とよばれ，粒子がいつまで自 分の速度を憶えているか，を大まかに示すものであった．例えば，Brown 運動をする粒子 の場合は， V AF (t)∝ exp [ −m γt ] (85) のように減衰することが示された（m は粒子質量，γ は抵抗係数）．これから，γ m程度の 時間がたつと粒子は自分が持っていた速度を忘れることがわかった．ここで取りあげる空 間自己相関関数 は，これと原理的には同じことを，空間について行うものであり，ある場 所での s (⃗r) の影響がどの程度遠くまで残っているかを示すものと考えればよい． また，同じく前期の「原子系の動力学セミナー」で紹介した 動径分布関数 radial distri-bution function も類似の空間相関関数の一種と考えてよい．（この場合，厳密には，「自己」 とは言えないが…）
#include <stdio.h> #include <math.h> #define NUM_CELL 200 // システムサイズ #define NUM_COR 100 // 相関関数を求める範囲 #define DEL_COR 1.0 // 相関関数の dr int main( ) { int i,j,i0,j0,is,id,jd,dist; int nstep; int num[NUM_COR];
double sum[NUM_COR], sum0=0.0; double state[NUM_CELL][NUM_CELL]; char cdum[100]; FILE *fdata,*fout; scanf("%d",&nstep); sprintf(cdum,"%6.6d.dat",nstep); fdata=fopen(cdum,"r");
if (fdata==NULL) {printf("File NOT Exist\n"); return -9;} sprintf(cdum,"%6.6d.cor",nstep); fout=fopen(cdum,"w"); for (i=0;i<NUM_CELL;i++) { for (j=0;j<NUM_CELL;j++) { fscanf(fdata,"%lf",&state[i][j]); } } fclose(fdata); for (is=0;is<NUM_COR;is++) { num[is]=0; sum[is]=0.0; } for (i0=0;i0<NUM_CELL;i0++) { for (j0=0;j0<NUM_CELL;j0++) { sum0+=state[i0][j0]*state[i0][j0]; }} sum0/=(NUM_CELL*NUM_CELL); for (i0=0;i0<NUM_CELL;i0++) { printf("%3d\n",i0);
for (j0=0;j0<NUM_CELL;j0++) { // まず自分の位置(i0, j0)を決めて for (i=0;i<NUM_CELL;i++) {
for (j=0;j<NUM_CELL;j++) { // 次に相手の位置(i, j)を決める id=i-i0;
if (id <-NUM_CELL/2) id+=NUM_CELL; // 周期境界条件 else if (id > NUM_CELL/2) id-=NUM_CELL;
jd=j-j0;
if (jd <-NUM_CELL/2) jd+=NUM_CELL; else if (jd > NUM_CELL/2) jd-=NUM_CELL;
dist=sqrt((double)(id*id+jd*jd))/DEL_COR+0.5; // 距離を求めて if (dist<NUM_COR) { num[dist]++; // ヒストグラム積算 sum[dist] +=state[i0][j0]*state[i][j]; } }} }} for (is=0;is<NUM_COR;is++) { fprintf(fout,"%7.3f %10.7f %10d\n",is*DEL_COR,sum[is]/sum0/num[is],num[is]); } return 0; } List 7.3: 空間自己相関関数を求めるプログラムの例． -0.50 -0.20 0 0.25 0.50 0.75 1 0 20 40 60 Sp a ti a l C o rre la ti o n Distance 2D TDGL (Conserved System) 000000 000100 000800 002700 006400 012500 021600 034300 051200 072900 100000 133100 172800 219700 図 7–16: 保存系 TDGL モデルの空間自己相関関数の例．
7.3
空間構造の時間変化
この節の詳細は，川崎 (参考文献 5) や 藤阪 (同 6) の教科書などを参照され たい． TDGL モデルは，T∗ < 0 の場合に，相分離が進行し，構造形成が行われる．この節で は，その時間変化の特徴を，簡単な解析（線形解析 や 次元解析）によって調べる． 7.3.1 前期過程：線形解析 秩序変数 s = 0 (すなわち高温で出現することの多い 無秩序状態disordered state) とい う初期状態から，温度を下げることによって出現する秩序状態ordered state のパターン の時間変化を考えよう．初めのうちは，s の空間的な揺らぎの成長が，パターンの変化となって現れる．線形モード解析linear mode analysis はこうした変化を調べるために利用で
きるツール tool である． s = 0 のまわりの揺らぎ δs(⃗r) を波数空間で考える： 波数空間で考える（＝ Fourier 変換す る）のは，ラプラシアン ∆ の扱いが 容易になるからである． δs(⃗r) =∑ k δ˜skexp[i⃗k· ⃗r] (86) TDGL 方程式 (69) または (82) に代入して，δ˜skの１次までの項を考える (これが linear
analysis である) と，それぞれの k について ∆ exp[i⃗k· ⃗r] = −k2exp[i⃗k· ⃗r] とな
ることに注意． ∂δ˜sk ∂t ≃ ( −T∗− ξk2)δ˜s k 非保存系 k2(−T∗− ξk2)δ˜s k 保存系 (87) となる．右辺 δ˜skの係数が揺らぎの成長速度を支配することになる．相変化を起こすのは これはもちろん，dxdt = ax の解が x(t)∝ exp[at] だからである． T∗< 0 であることに注意すると，その成長速度は線形解析の範囲内では図 7–17 のように なると予想される． (1) 非保存系では，k = 0 のモードの揺らぎ が最も速く成長する．つまり，長波長の揺 らぎが起こりやすく，その成長速度は， δ˜sk=0(t)∝ e−T ∗t (88) である．逆に，短波長 k > √ −T∗ ξ の揺らぎは急速に減衰することが予想される． -2 -1 0 1 2 0 0.5 1 1.5 2 G ro w th R a te Wave Number k Non-conserved Conserved 図 7–17: 波数 k のモードの線形近似での成長速度：T∗=−1.0，ξ = 1.0 の例．
(2) 保存系 では，k2についての２次式の形から，波数 k = √ −T∗ 2ξ の揺らぎが最も速く x(a− bx) は x = a 2b において最大値 a2 4b をとるからである． 成長し，その速度は δ˜sk(t)∝ e (T ∗ )2 4ξ t (89) である．これが，保存系において，周期的な空間構造（縞模様）が発達する理由で ある． 7.3.2 後期過程：次元解析 十分に時間が経過すると，パターンの変化は揺らぎの成長ではなく，界面の移動が主な 原因となる．非線形性（TDGL モデルなら，s3の項の影響）が強くあらわれるので，線形 解析は使えず，一般には解析は面倒になる．ここでは，簡単な次元解析dimension analysis による定性的な考察のみを紹介する． (1) 非保存系： 図 7–18 (左) のように，界面張力のために凸の部分が後退していく．そ の速度は大まかにいって曲率に比例（曲率半径に反比例）すると予想されるので，系 「曲率に比例」は大きな根拠があるわ けではない．曲率がゼロ（＝直線，平 面）なら界面は動かないはずので，曲 率に比例するのではないか，という単 純な予想である． の特徴的な長さcharacteristic length L に対して，L t ∝ 1 L となる．すなわち， L∝ t1/2 (90) となることが予想される．L としては，例えば空間相関長を考えればよい． (2) 保存系： 図 7–18 (右) のように，界面が変形するためには 凸部から凹部への物質移 動が必要である．物質量が保存されるために，界面の移動速度は曲率そのものでは なく，曲率の場所的な違い に支配されると期待できる．その結果，L t ∝ 1 L2，すな 次元解析から，曲率∝ L−1であり， その空間微分により，曲率の空間変化 率∝ L−2というわけである． わち L∝ t1/3 (91) となることが予想される．こうした議論により，保存系の緩和が非保存系よりも遅い ことが説明できる． 非保存系 保存系 図 7–18: TDGL モデルによる界面の運動（後期過程）の模式図．
