地域密着型エリア広告配信リクルートの地域密着型広告ドコイク?アドネットワーク Ads by ドコイク?[無料でホームページを作成] [通報・削除依頼]
[無料でホームページを作成] [通報・削除依頼]

フォノン

  ここではフォノンの計算に関する情報について纏める。
-----------------------------------------------------------------------
■ 定積比熱と熱伝導率
 温度Tで熱平衡にあるフォノン系の内部エネルギーは、調和近似で、
 U = ΣqΣj(h/2π)*ωj(q)*[nj(q)+1/2]= ΣqΣj(h/2π)*ωj(q)*[fBE +1/2]
で与えられる。ここでnj(q)は(q,j)で指定されるフォノンの平均の個数であり、ボース分布(ボース・アインシュタイン分布)関数
 fBE = nj(q) = 1/(e^(h/(2π*ωj(q)*kB*T))-1)
である。定積比熱を求めるためには、体積一定のもとで内部エネルギーUを温度で微分すればよい。
 温度が変化しても体積が一定であれば、ωj(q)に変化はなく一定に保たれるから、定積比熱Cvはnj(q)の温度変化によって決定される。すなわち、上記の2つの式から定積比熱Cvは、
 Cv = (dU/dT)v = kBqΣj((h/2π)*ωj(q)/(kB*T))^2*nj(q)*[nj(q)+1]
= kBqΣj((h/2π)*ωj(q)/(kB*T))^2*fBE*[fBE+1]
で与えられる。この式を用いて具体的にCvを求めようとする場合、Debye近似が用いられる。
 この近似は、絶縁体の比熱の温度依存性を説明するのに十分なものである。
 Debye近似ではCvに寄与するフォノンは、音響的分枝のフォノンが主なものであると仮定し、また音響的分枝の分散関係に対しては次式を仮定する。
 ωj(q) = cj*q
 さらにΣは第一Brillouin帯内で(q,j)のとりうる3rN個(Nは基本単位胞の数、rは基本単位胞内の原子の数)の値について和をとることを表すが、この和を振動数に関する次のような積分で置き換えて近似する。
 すなわちqについての和を積分に置き換え、これを[ωj(q) = cj*q]の関係を用いてωについての積分で置き換える。
 qとjに関する和を積分に書き換えると下記のようになる。
 ΣqΣj → ∫0ωDg(ω)dω
ただし、g(ω)=9*r*N*(1/ωD^3)*ω^2
 g(ω)dωはωからω+dωの中に振動数を持つような固有振動の状態の数を与えるから、g(ω)を状態密度と呼ぶ(0からωDまで積分すると全状態数の3rNになる)。またωDはこの近似でとりうる最大の振動数を与え、Debye振動数と呼ばれる。
 以上のDebye近似を用いると、定積比熱Cvは、次のように書き直すことができる。
 Cv = 9*r*N*kB*(T/ΘD)^3*∫0ΘD/T[(x^4*exp(x))/(exp(x)-1)^2]dx
 ただしΘDは、
 ΘD = ((h/2π)*ωD)/(kB*T)
で定義され、これをDebye温度と呼ぶ。十分低温では上記の積分の上限は∞と近似してよいから、
 Cv = 9*r*N*kB*(T/ΘD)^3, T =< ΘD/25
となり、比熱はT^3に比例することがわかる。また十分高温では、ΘD/Tのべきに展開して、
 Cv=3*r*N*kB[1-(1/20)*(ΘD/T)^2+...], T > ΘD
が得られるので、高温での比熱はDulong-Petitの値3*r*N*kBに漸近していくことがわかる。
 熱伝導度Kはフォノンの群速度vと緩和時間τを用いれば次のように表される(phon3pyのoutput fileに書かれている式とほぼ同じ)。
 K = (1/3)*Cv*v^2*τ
群速度vはフォノンの分散関係の傾きから、緩和時間は非調和振動項(例えば3次)を加味して得られる。
日本金属学会編「金属物性基礎講座」丸善

◇ 私の勝手な見解
 フォノンの分散関係を見ると、多くの場合、光学フォノンは音響フォノンよりも傾きが緩やかであることから群速度は遅い(また傾きが緩やかであることから光学フォノンのDOSは大きい)から、緩和時間τが同じであれば、熱伝導は小さくなることが予測される。
 十数%までの重元素置換では、多くの光学フォノンが音響フォノンに重なるようにはならないが、音響および光学フォノンの群速度は低下し、低いエネルギー側へとシフトした形になる。そのため、ΘDより十分に高温側では、定積比熱Cvの温度による増加量は小さくなり((ΘD/T)で展開すると Cv=3*r*N*kB[1-(1/20)*(ΘD/T)^2+...])、低くなった群速度と相まって熱伝導率が低下する。
 定積比熱Cvや熱伝導率Kの温度依存性はボース分布関数fBEによって影響を受ける。フォノンの分散またはDOSにfBE*[fBE+1]、つまり(fBE2+fBE)を掛けて積分すればCvが、さらに群速度の2乗と緩和時間も掛けて積分すれば熱伝導度が得られるからである。
 低いエネルギーは主に音響フォノンが主である。低いエネルギー側に光学フォノンが入るような方法を探すか、重元素によって群速度が遅くなるようにすることが必要になる。
 平均自由行程 l = v*τとなるので、もし重元素置換で、緩和時間がそれほど変わらず、群速度vが0.9倍になったら、粒界散乱の効果を得るためには単純計算で置換前の0.9倍だけ距離を縮めなければならない。
※ 焼結で粒界拡散の距離を縮める場合、得られた相対密度から電気抵抗率を算出するための体積の補正が必要なのではないかと思っています。気孔がある場合にはその点を補正してみたらと思っています(多くの先生方や研究者の報告はもっとよい値が出ているのではないでしょうか?)。

◇ 状態密度分布の合計(フォノンのDOSの合計)
 各原子からx,y,z方向にバネが出ていると簡単に考える。そうすると、原子数Nにx,y,zの合計3のバネがあるから、全状態数は3Nとなる。実際、フォノンDOSの面積を合計すると3Nとなる。バンド分散の数(固有振動のモードの数=分枝)も3Nとなる。結晶中の全イオンの運動の自由度も3Nである。
※ Debye近似では、Debye近似での状態密度と実際のphononの状態密度は形状が異なるが、総面積が3Nであることは変わらないこと、そこそこ定積比熱の振る舞いを説明できる点などがある。

◇ MEMO
・relaxation time approximation
   -∂n/∂t = Tr( (n-equilibrium(n))/t )
 n = phonon occupation number]
 equilibrium(n) =1/(eℏw/kBT-1)
・Single mode relaxation time approximation (SMRT)
 Single modeでは対角(Tr)を取る。よって下記となる。
 -∂n/∂t = Tr( (n-equilibrium(n))/t )
・Fermi’s golden rule > Im (self energy) = Gamma
 Phonon lifetime t = 1/2Gamma
・ボース分布関数は50 Kで0から30 meV(=7.25 THz)程度まで広がりを持つ。
 1 THz = 4.13567 meV
-----------------------------------------------------------------------------
◆ レクチャー
Tutorial 3a: Materials Simulation by First-Principles Density Functional Theory I : https://www.youtube.com/watch?v=rx-pL3fkS5o
-----------------------------------------------------------------------------
◆ 各計算コードの特徴
□ Phonopy
  フォノン分散やフォノンDOSを計算可能。V_Simで原子の動きが見える(ANIMEで行う。これは未確認)。計算しやすい。
□ PWscf
  計算が速い。フォノン分散やフォノンDOSを計算可能。計算しやすい。
□ VASP
 計算が速い。フォノン分散やフォノンDOSを計算可能。かなり計算しやすい。
□ PHON
  Γ点におけるフォースと原子の動きを知ることができる。フォノン分散やフォノンDOSを計算可能。私は計算に成功していない。
□ PHASE
  Γ点のみ計算可能。IR及びRamanのどちらに活性があるかを議論できる。私は計算に成功していない。
-----------------------------------------------------------------------------
◆ 格子振動を計算することで以下の情報を得ることができる。
・比熱
・熱膨張
・熱伝導
・強誘電体
・誘電体の光物性
・相転移
-----------------------------------------------------------------------------
・フォノン分散を描き、周波数が負になっているところが、振動数が虚数であることを意味する。
・虚数の振動数が計算結果で現れた場合には構造の変化が予測される。
・虚数の振動数が現れるのはエネルギー的に不安定な構造になっているため。
※ ポテンシャルVの二回微分であることを思い浮かべると良い。下に凸の場合(エネルギー的に安定であることが分かる)は、二回微分すると正の値になる。一方、上に凸の場合(は他の位置の方がエネルギーが低くなるポテンシャルになるので、エネルギー的に不安定であることが分かる)は、二回微分すると負の値になる。
-----------------------------------------------------------------------------
・振動数が得られることにより「振動のエントロピー」や「ヘルツホルムの振動のエネルギー」、「比熱」を計算することが可能となる。
-----------------------------------------------------------------------------
◆ フォノンの計算法
・ The small displacement method 
  supercell で計算することが必要になる。しかし、MgOなどではsupercellを大きくしても実験値と合わない[PH1]。
・ Linear response [LR1, LR2]
  primitive cell ( or conventional cell ) のみの計算で比較的精度の高い計算が出来る(私の誤りで無ければだが・・・・・・)。PWscfはこの計算が可能。VASPも ver.5.1以降で計算可能になった。IBIRON=8 では対称性を考慮して計算できるために対称性がある系では計算速度が速い。
[PH1] http://www.homepages.ucl.ac.uk/~ucfbdxa/pubblicazioni/COMPHY3792.pdf
[PH2] http://icms3.weebly.com/uploads/3/5/9/0/3590130/version1.pdf
[LR1] http://www.mcc.uiuc.edu/summerschool/2007/qmc/tutorials/Tutorial_phonons_alfe.pdf
[LR2] http://cms.mpi.univie.ac.at/vasp/vasp/IBRION_7_IBRION_8.html
[LR3] http://www.equilibriumtrix.net/2011rioworkshop/Palumbo_lecture_harmQHA.pdf
-----------------------------------------------------------------------------
◆ Frozen Phonon 法
1. ポテンシャルを計算し、力の定数を計算する。
2. 調和振動近似を導入(調和振動のもとで運動方程式を解く)して、基準モードと固有振動数を計算する。
* 最近は、quasi-harmonic approximation (the prefix “quasi” is there to indicate that the force constant matrix may depend on volume) [PH1-PH3]というのもある。恐らく、準調和振動近似とでも呼ぶのであろうか。フォース(力)の定数(より正確には行列にしたもの)が体積に依存するとして調和振動のもとで運動方程式を解く。
[PH3] http://phonopy.sourceforge.net/qha.html
-----------------------------------------------------------------------------
◆ 温度を考慮した全エネルギーの探索
・温度を考慮した全エネルギー = 通常の第一原理計算による全エネルギー + ヘルツホルムの振動の自由エネルギー
・上記の式で、ある格子定数での通常の第一原理計算による全エネルギーとヘルムホルツの振動の自由エネルギーで絶対温度を変化させたときの全エネルギーの和をプロット(横軸を格子定数、縦軸をエネルギー)する。この計算を格子定数を変えて行う。そして、一番エネルギーの低い状態(各温度と格子定数の組)を探す。
※ ヘルツホルムの振動の自由エネルギー = -kBT ln Z
  ここでkBはボルツマン定数、T は絶対温度、Zは分配関数。[1]
・フォノンのDOSから求める方法
  温度を考慮した全エネルギーE(T) = 通常の第一原理計算による全エネルギー Eel(E) + Eph(T)
  Eph(T)=(1/2)Σq (h/2π)Wq + kBq ln(nq+1)
ここでWqは振動数、nqはWqに対する状態数。実際には縦軸にフォノンのDOS(nq)、横軸にWqの図を描いて、上記の式に当て嵌めて振動数Wqを0から無限までにして計算すればよい。[2]
[1] http://www.google.co.jp/url?q=http://repository.kulib.kyoto-u.ac.jp/dspace/bitstream/2433/159093/1/Keisan%2520Kogaku_14_2_995.pdf&sa=U&ei=qwQ4UbvUBYPdkAWb9oCIBw&ved=0CCIQFjAC&usg=AFQjCNEHS8BeNf3eKRIV8zyD9NNPY7ioEw 
[2] 「計算機マテリアルデザイン入門」、大阪大学出版会
-----------------------------------------------------------------------------
「1 eVの平均運動エネルギーをもつ気体の温度は11 604 K」であることを理解しておくと良い。常温だと約0.03 eVに相当する。そのため、大抵のフォノンに関する計算では、Frozen Phonon法(フローズンフォノン法)で十分であろう。
[1] http://ja.wikipedia.org/wiki/%E9%9B%BB%E5%AD%90%E3%83%9C%E3%83%AB%E3%83%88 
-----------------------------------------------------------------------------
◆ 力の定数の計算
Fi = ∂ Etot / ∂ R =  − ∂ <φ|H|φ> / ∂ R = − <φ|∂H|φ> / ∂ R
ここで、波動関数に平面波を用いてる場合には、平面波では波動関数φが原子位置に依存しないため、∂ φ / ∂ R = 0 となる。平面波で無い場合はpulay 補正[1]が必要となる。
[1] http://www.lehigh.edu/~kk04/Theory_of_Periodic_Systems-DFT-FP-LAPW.pdf
-----------------------------------------------------------------------
■ 調和近似
 ポテンシャルの展開式において、変位uについて2次の項のみにとどめる近似を調和近似と呼ぶ。
 この近似では、イオンは平衡点のまわりで調和振動を行うからである。調和近似の範囲内では、結晶の膨張係数は0であり、定圧比熱と定積比熱は等しく、弾性定数は温度依存性を示さないなど、実際の結晶の性質とは必ずしも合致しない結果が導かれる。それにもかかわらず、調和近似は格子振動を理解するために最も基本的で重要な概念を与えるものである。
 さらに付け加えるべきことは、調和近似の形式はそのままにして概念を拡張することにより、この近似の適用範囲を広げることができることである。たとえば、調和項における微分を、結晶の温度における格子点で微分した値をとることにすると、力の定数は温度に依存して値が変化すると考え直すことができるから、調和近似の枠内でも物理的に有用な考察が可能となる。
 このような近似を本来の近似と区別するために、特に準調和近似と呼ぶこともある。
日本金属学会編「格子振動・非線形光学」丸善
-----------------------------------------------------------------------
■ 非調和項
 非調和項が、調和近似で求められた固有振動モード間の非調和相互作用エネルギーを表す項であるとみなしてもよい。 
 フォノンの描像では、非調和項はフォノン同士の相互作用を表すものである。
 3次の非調和項は3個のフォノンの衝突相互作用エネルギーを表し、運動量保存則が成り立つ。
  (h/2π)*(±q1 ±q2 ±q3) = (h/2π)*K
ここで、(h/2π)*q1,(h/2π)*q2,(h/2π)*q3は3個のフォノンの運動量をそれぞれ表し、Kは0または任意の逆格子ベクトルである。±の符号はそれぞれフォノンの消滅および生成に対応する。
 K=0の場合はふつうの運動量保存則で、この場合の衝突過程をN過程(Normal process)と呼ぶ。
 K≠0の場合をU過程(Umklapp process)と呼ぶ。(h/2π)*Kを格子系が担うべき運動量(擬運動量)と見なせば、上式を擬運動量保存則と呼んでも良い。U過程は結晶の併進対称性の結果可能となった過程である。
 4次以上の非調和項も同様に4個以上のフォノンの衝突過程に対応し、関与するフォノンの運動量の間には、擬運動量保存則が存在する。
 フォノンの描像は元来調和近似に基づいてつくられたものであり、フォノンが互いに独立な粒子として存在するかぎり一定のエネルギーを持っている。
 しかしながら、非調和項を通してフォノン同士の衝突が可能となると、フォノンはもはやもとの形のままで存在し続けるとは限らず、ある時は消滅して2個のフォノンに姿を変え、その後また2個のフォノンの姿を消して再びもとのフォノンに立ち返るというな種々の変貌を繰り返すことが量子力学的に可能である。この場合もとのフォノンの姿に立ち返っている時間が十分長ければ、一定の振動数を持った粒子として観測することが可能であり、したがってやはり量子数で指定されるフォノンという”粒子”像をもち続けることができる。
 このときは結晶の全格子振動エネルギーは、このような新しいフォノンの個数にそのエネルギーを掛けて、すべてのとりうる値について総和したもので与えられると考えてよい。新しいフォノンは非調和項の影響を受けたものであるから、一般にはそのエネルギーは調和近似でのフォノンのエネルギーとは異なった値になっていることが期待される。
 一方、非調和項の影響が、もとのフォノンの姿を変貌させる効果を持つということは、新しいフォノンが有限の寿命を持つ効果として現れる。
 この新しい”準粒子”はその寿命がある程度長いかぎり”粒子”としての意味を失わない。
 一般に格子振動を表す結晶のハミルトニアンは非調和項を含んでいるから、実際に観測されるフォノンは調和近似でのフォノンではなく、むしろ非調和の影響のため有限の寿命をもった準粒子としてのフォノンであると考えるのが正しい。
日本金属学会編「格子振動・非線形光学」丸善
-----------------------------------------------------------------------
■ 調和近似における各イオンの変位に対する運動方程式
フックの法則から、力Fと変位xを用いて下記のように表される。
 F = -k*x
ここでkはばね定数(または力の定数)と呼ばれる。ばねに蓄えられたエネルギーU(=ポテンシャルエネルギー)は
 U = (1/2)*k*x^2
と表される。
 各イオンの変位に対する運動方程式は、変位を単純なxから一般化してuとすると下記のようになる。
 Mi*∂2uα(i)/∂t2 = -ΣjΣβ Φα,β(i, j)*uβ(j)
 Φα,β(i, j) ≡ ∂2U/(∂uα(i)∂uβ(j))
ここでMiはi番目の原子の質量、Φα,β(i, j)は上式で定義される力の定数である。
日本金属学会編「格子振動・非線形光学」丸善
-----------------------------------------------------------------------
■ 分散関係
 波数ベクトルqと振動数との関係ωj(q)を固有振動モードの分散関係と呼ぶ。
 この分散関係を理論的に求めるためには、ダイナミカル・マトリックス Dが与えられなければならない。
日本金属学会編「格子振動・非線形光学」丸善
-----------------------------------------------------------------------
■ ダイナミカル・マトリックス
 調和近似において、各イオンの変位に対する運動方程式は、一般に平面波の合成で表すことができる。
 Mi*∂2uα(i)/∂t2 = -ΣjΣβ Φα,β(i, j)*uβ(j)
 そのような一般解を求めるために、まず特解を求めてみる。波数ベクトルq,振動数ωを持った波が解の一つであると仮定すると、下記のような永年方程式を得ることができる。
 (ω^2*1-D)U=0
ここで1は単位行列、Uは1列3r行の行列であり、ダイナミカル・マトリックスDは3r x 3rの行列である。
 この式が意味のある解を持つための必要条件は、永年方程式、
 |ω^2*1-D|=0
が満たされることである。左辺は3r行3r列の行列式であるから、ω^2について解くと、与えられたqに対して一般に3r個の異なる解が得られる。この3r個の振動数wj(q)を固有振動数と呼ぶ。
 分散関係を理論的に求めるためには、ダイナミカル・マトリックスが与えられなければならない。
日本金属学会編「格子振動・非線形光学」丸善
-----------------------------------------------------------------------
■ ハミルトニアン
 全系のハミルトニアンは、フォノンの粒子像で描き表されることが分かった。
 この量子化の手続きは、電磁場を量子化して格子像を得たのとまったく並行的であることを注意しておく。
 したがって光子がBose統計に従う粒子であるのと同じく、フォノンもBose統計に従い、個数演算子の固有値として許される値は0,1,2,3...などの正の整数である。
 このことは物理的に、同じ量子数(q,j)で指定される一つの状態を占めるフォノンの個数は、整数ならば何個でも許されることを意味する。それゆえ有限の温度Tで熱平衡状態にあるフォノン系では、(q,j)状態を占めるフォノンの平均個数は、
 n(q,j) = 1/[exp((h/2π)*ωj(q)/kbT)-1]
で与えられる。このような分布をBose分布と呼ぶ。このようにフォノンと格子は、統計的な性質では類似しているが、次の点で明らかに区別される。
 すなわち光子の場合、エネルギーは(h/2π)*ω(q)=(h/2π)*c*q (cは光の速さ)であるが、フォノンのエネルギー(h/2π)*ω(q)は、結晶特有の分散関係で与えられるq依存性を持つ。
 また当然のことながらフォノンは真空中では存在しえない。さらにフォノンには横波のほかに縦波フォノンも存在することを忘れてはならない。
日本金属学会編「格子振動・非線形光学」丸善
-----------------------------------------------------------------------
■ フォノンの運動量
 ここでフォノンの持つ運動量について次のことを注意しておく。結晶の並進対称性の結果、任意の逆格子ベクトルをKとして、
 ωj(q)=ωj(q+K)
が成り立つから、(q+K,j)で指定される固有振動モードと(q,j)で指定されるそれとは互いに独立なものではない。
 これに対応して(q+K,j)フォノンと(q,j)フォノンは同じエネルギーを持つ、本来別種類のフォノンであるとは言えない。言い換えれば、フォノンを指定する波数ベクトルqは、本来任意の逆格子ベクトルKだけの任意性を持っている。
日本金属学会編「格子振動・非線形光学」丸善
-----------------------------------------------------------------------
■ 縦波と横波分枝
 3個のフォノンの関与するN過程を考えると、運動量保存則
 q1 + q2 = q3
およびエネルギー保存則
 ω1 + ω2 = ω3
を同時に満足しなければならないという条件から、3個のフォノンはどの分枝に属するものでもよいというわけにはいかない。
 すなわち、上式を満足するには三つのベクトルが三角形を作ることを表すから
 q1 + q2 >= q3
が要求される。一方、縦波分枝の方が横波分枝よりも同じqに対してω/qの値が大きいことから、運動量およびエネルギー保存則と矛盾しないためには、横波分枝をT, 縦波分枝をLと書いて、
 T + T ←→ L
 T + L ←→ L
の過程しか許されないことがわかる。したがって、横波と縦波のフォノンの緩和時間は必然的に異なったωと温度依存性を持つことが期待される。
日本金属学会編「格子振動・非線形光学」丸善
-----------------------------------------------------------------------
■ イオン結晶の分散関係
 アルカリ・ハライドのようなイオン結晶の分散関係について、一般的に次のことがいえる。
 基本単位胞内に二つのイオンを含むから、固有振動モードは6個存在し、そのうち3個はq→0, ω→0となる。
 この3個のモードは小さいqの値に対して、振動数ωがqに比例するから音響的分枝と呼ばれる。
 またこの分枝では正負のイオンが同位相であり、電気的中性の連続媒質内の振動波と同じようにふるまうから、正負イオンの重心の運動が振動数を決めていると考えてよい。一方、qの大きい値に対する振動数は、重いイオンの質量の平方根に逆比例することが予測される。
 音響的分枝は一つの縦波と二つの横波の分枝とから成るが、逆格子空間で回転対称性を持つqの方向では、この二つの横波分枝は縮退している。音響的分枝以外の三つの分子は、正負イオンがπだけずれた位相で相対的に運動するモードである。
 振動数はq→0で有限の値にとどまり、かつ正負イオンの還元質量の平方根に逆比例する。
 このことはイオンのふるまいが、正負イオンからなるプラズマの振動と同じであることから理解される。
 このモードはq〜0で正負イオンが電気的双極子能率をつくるような振動であり、光学的に観測されるから光学的分枝と呼ばれる。qの大きな値に対しては、軽い方のイオンが振動数の値を決定すると考えてよい。
 なお光学的縦波分枝と光学的横波分枝の振動数は、q→0の極限で一致せず、Layddane-Sachs-Tellerの関係式
 ωLO/ωTO = (ε0/ε∞)^(1/2)
が成り立つ。ここでε0は静電的誘電率、ε∞は光学誘電率である。縦波に対しては、媒質の長波長の分極が付加的な力として作用するから、その分だけ力の定数が大きくなり、縦波振動数ωLOが横波振動数ωTOより高い振動数となるわけである。
 なおこれらの一般的な性質は、一次元格子のモデルを用いて容易に理解することができる。
日本金属学会編「格子振動・非線形光学」丸善
-----------------------------------------------------------------------
■ デバイ公式
 デバイ公式は定量的には、全温度領域にわたって実測値と必ずしも一致するというわけではない。実際にどれだけの差異を生じるかを見るために普通採用される方法は、Debye公式を内挿式のように考え、ΘDを経験的なパラメーターと見直して、これを実験的に求めるやり方である。
 すなわちΘDのかわりにパラメーターΘを導入し、実測の比熱が、
 Cvexp(T) = Cvdebye(Θ/T) ≡ 9*r*N*kB*(T/Θ)^3*∫0Θ/T[(x^4*exp(x))/(exp(x)-1)^2]dx
で表されるとして、逆にΘを測定温度Tの関数として求める。もしDebye近似が正しければ、ΘはTによらず一定の値ΘDに等しく求まるはずであるが、実際にはΘ(T)がTに依存して変化する結果が得られる。
 一般的にいえることは、このようにして実験的に求めたΘ(T)の値は、低温の極限においては、音速の測定値を用いて求めたDebye温度ΘDと良い一致を示している。しかしながら、ほぼデバイ温度/50くらいの温度以上で実測値と異なってくる。これは実際のフォノン分散関係が(状態密度も含めて)デバイ近似のものから外れてくるためであると考えられる。
 実測値と一致しなくなる原因としては、以上のほかに非調和項の影響が考えられる。
 デバイ公式は調和近似で求めたものであるから、非調和項の影響が無視できなくなると思われる高温の領域に対して、
デバイ公式を理論的に補正しなければならない。
日本金属学会編「格子振動・非線形光学」丸善
-----------------------------------------------------------------------
■ 熱伝導度
 熱伝導現象では結晶全体が熱平衡にあるわけではなく、温度勾配が存在するから、一定の量子数(q,j)で指定されるフォノン(対応する固有振動は平面波で表される)のかわりに、ある場素xのまわりに小さい広がり|凅|をもって局在したフォノンの塊り(波束)を考えなければならない。
 波束のもつ波数ベクトルをqとする。これはqに分布の極大を持ち、そのまわりに|冫|の広がりを持つような重みをつけて、固有振動を重ね合わせたものでつくられていると考えてよい。|冫|は|冫|・|凅|〜1を満たす程度の広がりである。
 波束に対応するフォノンの分布関数Nj*q,x)を導入すると、熱流Qは
 Q = ΣqΣj Nj(q,x)*(h/2π)*ωj*∂ωj(q)/∂q
で与えられる。ここで∂ωj(q)/∂qは波束の群速度である。これはqを-qで置き換えると符号を変えるから、分布関数Nがq→-qの置き換えで不変であるような形であれば上式のqについての和はゼロ、すなわちQ=0となることがわかる。特にNが熱平衡分布
 Nj0(q) = 1/[exp((h/2π)*ωj(q)/kB*T)-1]
のときは明らかにQ=0である。このように分布関数Nが熱平衡分布N0から非中心対称的にずれた分だけ熱流に寄与するわけであるが、このずれは一つには結晶内の温度勾配によって生じる。一方、有限な定常熱流を得るためには、分布関数Nを熱平衡分布に戻そうとする緩和機構が存在しなければならない。
 完全結晶の場合のフォノン同士の衝突過程には、N過程とU過程が存在する。もしN過程だけが許されるとすると、これによってはフォノン系の全運動量は保存されるから、温度勾配によってフォノン系が一度運動量を得れば、いつまでも保存され、熱伝導度は発散することになる。実際に、任意の速度λで流れるフォノンの定常分布
 Njλ(q)= 1/[exp((h/2π)*(ωj(q)+λ*q)/(kB*T))-1]
は運動量を保存する分布であり、フォノンがこの分布関数をとれば、N過程によって熱平衡分布に緩和させることはできないので、この状態のままでは熱伝導度は無限大となる。
 しかしながら、だからといってN過程は熱伝導度の計算で捨て去ることはできない。なぜならN過程は、種々の状態を占めるフォノンの分布関数を変える働きをするから、たとえそれ自身は熱抵抗を与えなくても、熱抵抗を与えるような他の散乱過程と組み合わさって緩和に影響するからである。熱抵抗を与える散乱過程としては、U過程のほかに試料表面によるフォノンの散乱、不純物、格子欠陥などの格子の不完全性による散乱過程が考えられる。
 これらが熱伝導にどのように寄与するかについて以下で述べることにする。
 一般に定常状態でのフォノン分布関数を決めるBoltzmann方程式は、
 vq,j*∇Nj(q,x)=[∂Nj(q,x)/∂t]
で与えられる。vq,jは群速度である。左辺は温度勾配があるためのフォノンの輸送による分布関数の変化を表し、温度勾配について線形な近似では∇Nj(q,x)は
 ∇Nj(q,x)/ = ∇T*[∂Nj0(q,x)//∂T]
で置き換えてよい。上式の右辺はフォノンの散乱過程による分布関数の変化の割合を表し、フォノン衝突項とよばれる。この項は、フォノンどうしの散乱過程を表す非線形積分項を含み複雑な形をしているので、ふつうこれを線形化する近似が用いられる。
 分布関数が熱平衡分布に戻る緩和は、この衝突項によって決まるから、緩和現象を特徴づける緩和時間パラメーターτjeff(q)を導入して、形式的に
 [∂Nj(q,x)/∂t] = -(Nj(q,x)-Nj0(q))/τjeff(q)
と書けるものと仮定する。このときBoltzmann方程式の解Nj(q,x)は求まるから、これを熱流の式(上記のQを求める式)に代入し、熱伝導度テンソル
 Kα,β = ΣqΣ (vq,j)α* (vq,j)βjeff(q)*Cj(q)
を得る(phon3pyのoutput fileに書かれている式と同じ)。ここで、
 Cj(q) ≡ ∂[(h/2π)*ωj(q)*∂Nj0(q)]/∂T
とおいてあるが、これは(q,j)フォノンの比熱への寄与を表す。
日本金属学会編「格子振動・非線形光学」丸善
-----------------------------------------------------------------------
■ 熱伝導度の測定結果を解析するためによく使われる近似
 実際に熱伝導度の測定結果を解析するためには、次のような近似が用いられる。すなわち、群速度v(q,j)を平均の速度vで置き換え、かつフォノン・スペクトルはすべての分枝に共通にω(q,j)=v*qで表されるとする。
 ただしqのとりうる範囲はちょうど運動の自由度で許される0からqD(デバイ波数)までの範囲とする。また結晶は等方的であるとして、vα vβの空間平均をとりδα,β*v^2/3で置き換える。qについての和は積分で置き換え、jについての和は単に分枝の数3を乗ずることで置き換えることにすると、熱伝導度は、
 K = kb/(2π^2)*(2π*kb*T/h)^3*(1/v)*∫τ(x,T)*[(x^4*exp(x))/(exp(x)-1)^2]dx
 x = ω(q)/(2π*kb*T)
ここで積分は0からデバイ温度/Tまで取る。
と変形される。この式から分かるように、重要なkとは熱伝導度はもっぱら緩和時間パラメーターτによって決まることである。
 したがって緩和時間パラメーターがフォノンの散乱過程とどのような関係にあるかを知らなければならない。
※ フォノンの群速度(dω/dq)は電子の時と同様に、分散関係の傾きから求めることができる。
日本金属学会編「格子振動・非線形光学」丸善
-----------------------------------------------------------------------
■ N過程と散乱に直接寄与する過程
 フォノンの散乱過程の中には、U過程、境界、不純物、格子欠陥などのよる散乱のように熱抵抗を与える散乱過程と、それ自身では熱抵抗を与えないN過程がある。したがってこれらの役割を分類することができれば、フォノンの熱伝導についてさらに深く知ることができる。Callawayに従って、フォノン衝突項を二つの部分に分け
 (∂N/∂t)=-(N-N0)/τR - (N-Nλ)/τN
で表せると仮定する。右辺第一項は熱抵抗を与えるような散乱過程を通して、分布関数が緩和時間τRで平衡に近づくことを表し、第二項はN過程にによって下記の式
 Njλ(q)= 1/[exp((h/2π)*(ωj(q)+λ*q)/(kB*T))-1]
で与えられる分布に近づく割合を表す。熱抵抗を与える散乱過程がいくつかある場合は、その過程による緩和時間をτiとして、τR
 (1/τR) = Σi (1/τi)
で与えらえる。これらの式とDebye近似を使って熱伝導度の式を表すと、
 K = K1 + K2
となる(K1とK2の詳細は参考文献を参照してください)。
1)熱抵抗を与える散乱過程(R過程と呼ぶ)が支配的な場合
 この場合はR過程による緩和時間がN過程のそれより速いから、τN >> τRであり、K2を無視ることができる。
 結局これはτをτRで置き換えたものになる。
2)N過程が支配的で、τN << τRの場合
 K1を無視することができる。
 この場合、τでは理解が難しいため、熱抵抗 W ≡ 1/Kで考える方が理解しやすい。
 すなわち、(1/τR) = Σi (1/τi)に注意して、
 W = Σi Wi
と表すことができる。結果的に式の中に、τNは入らないが、N過程が十分速いために、熱抵抗は各R過程による抵抗の和で与えられる。
 このように熱伝導度がR過程でどのように決まるかはN過程の緩和時間τNの速さに依存し、N過程が十分速く、フォノン系の中でのエネルギーの再配分が速やかに行われるときは、熱伝導度のかわりに熱抵抗に着目し、各R過程がどのような抵抗の寄与をするかを考えればよい。
 逆の極限でN過程が遅いときは、各R過程がどのように緩和時間τRに寄与するかによって熱伝導度が決まる。
 τNは非平衡フォノン分布が、任意の速度λをもった定常分布Nλに緩和する時間の目安を与える。
日本金属学会編「格子振動・非線形光学」丸善
-----------------------------------------------------------------------
■ 熱伝導度の温度依存性
 N過程は種々の状態を占めるフォノンの分布状態を変える働きをするから、たとえそれ自身は熱抵抗を与えなくても、熱抵抗を与えるような他の散乱過程と組み合わさって緩和に影響する。
 熱抵抗を与える散乱過程としては、U過程のほかに試料表面によるフォノンの散乱、不純物、格子欠陥などの格子の不完全性による散乱過程が考えられる。
(省略)
 重要なことは熱伝導度はもっぱら緩和時間によって決まることである。したがって、緩和時間がフォノンの散乱過程とどのような関係にあるかを知らなければならない。

1)十分低温でフォノンの平均自由行程が長くなると、試料表面での散乱が熱伝導度を決めるようになる
 緩和時間^-1 = フォノンの群速度の平均 / Lc
 ここで、LcはCasimirの長さと呼ばれ、半径Rの円筒試料ではLc=2R, また一辺Dの角形の棒状試料ではLc=1.12Dととられる。
 この過程での緩和時間は角振動するによらないから、十分低温であるT<<デバイ温度では、熱伝導度はT^3に比例する。

2)同位元素による散乱
 同位元素は、力の定数が等しく質量のみが異なる不純物と考えてよいから、フォノンはこれによってRayleigh散乱を受ける。

3)U過程による散乱
 簡単のため3個のフォノンの関与する散乱過程を考える。いまq'とq"の2個のフォノンが消滅して1個のqフォノンが生まれる過程を例にとると、保存則は、
 q-q'-q"=K≠0
で与えられる。定温領域においては、フォノンは低振動数領域に多く分布し、したがって主要なフォノンの|q|は十分小さい。
 それゆえ上式を満たすようなU過程が可能であるためには、q'+q"が第1 Brillouin帯の外に出るようなものでなければならないから、|q'|,|q"|ともに|q|程度の小さい値をとる場合は除外される。このことは3個のフォノンのうちいずれかは、ある程度大きなqあるいはせいぜいBrillouin帯の領域の1/4くらいのqはもっていなければならないことを意味する。
 したがって、デバイ波数をqDとすると、このようなフォノンのqの大きさはq〜qD/2程度と考えてよい。
 T << デバイ温度の場合
  散乱確率はあるいは緩和時間の逆数の主な温度依存性は、
  緩和時間^-1 = A*exp(-デバイ温度/(a*T))
  ここでA〜2である。
 高温領域 T>デバイ温度の場合
  分布関数はTに比例するから、緩和時間^-1 ∝ Tとなり、熱伝導度はT^-1に比例することが示される。

4)N過程による散乱の寄与
 N過程だけでは熱抵抗を生じないが、異なる振動モード間にエネルギーを再分配させる役割を果たすから、他の散乱過程による緩和に影響を与える。

5)点欠陥による散乱
 空格子点や格子間不純物のような点欠陥は、母体結晶との質量の違いあるいは原子間力の違いによる格子の乱れの他に、点欠陥のまわりに弾性的なひずみを生じ、これがフォノン散乱の主因になる。ひずみの場を最も簡単に
 U(r) = A*(r/r)*(1/r^2) (r>r0)
 U(r) = 0 (r<r0)
で表す。ここでr0は原子間距離程度の大きさに選ぶ。このような点対称の場による散乱過程に対しては、
 τp^-1 〜2*10^2*Np*(γ*A)^2*ω^4/vav^3
となることが示される。ここでNpは点欠陥の濃度、γはGruneisen定数である。
 同位元素による散乱の場合と同じく、τpはω^-4に比例し、ω→0で強い発散を示すが、N過程を同時に考えるとこの発散は弱められ、結果的には、点欠陥による散乱はTに比例する項を熱伝導度に寄与することが分かる。

6)転位による散乱
 Gruneisen定数が、刃状転位による散乱での緩和時間に影響する。
 積層欠陥については、角周波数の2乗に反比例する。
 小傾角境界は結晶粒の境界に垂直なBurgersベクトルをもつ刃状転位が、等間隔dで並び、境界の両側の結晶の傾き角Θ〜b/dが小さいものをいう。dがフォノン波長に比べて小さい場合には、熱伝導度にはそれほど重要な寄与を与えない。
 
 以上、主な緩和過程について述べたが、高純度で完全に近い絶縁体結晶の場合には、熱伝導度の温度依存性について一般的に次のことがいえる。
 すなわち低温側では、試料表面の境界散乱によって決まるT^3依存性を持つち、中間温度領域で極大を示す。さらに高温側では、U過程が重要となり、温度依存性は(T^n)*exp(デバイ温度/(a*T))で表される。 また十分高温領域では1/Tに比例するようになる。
 熱伝導度が極大を示す温度領域では、同位元素、不純物散乱などによる緩和機構に敏感である。
 熱伝導度の測定においては、高温側から試料を通過してくるフォノンの流れは単色ではなく、幅の広いスペクトルを持っていること、第二にこのスペクトル分布の極大を示す振動数が温度に依存するから、試料の温度を変えることはスペクトル分布の極大を示す振動数を変えることに相当していることである。
 仮にフォノン分布を熱平衡分布で近似すると、振動数がω〜ω+dω内にあるフォノンのもつ平均のエネルギーは
 (h/2π)*ω^3*[exp((h/2π)*ω)/kbT-1]^-1*dω
で与えられ、(h/2π)*ω〜2.8kbTに極大を持っていることが分かる。
日本金属学会編「格子振動・非線形光学」丸善
-----------------------------------------------------------------------
■ 金属の熱伝導度の温度依存性
 フォノン―フォノン散乱のN過程による緩和が十分速やかな場合には、フォノンの熱伝導度Kよりもむしろ熱抵抗W=1/Kで考えた方がよい。
 Wは種々の散乱過程からの熱抵抗の寄与の和で与えられる。すなわち、T < ΘDに対して、
 (1/K) = Wge + WB + Wp + Wu + W
と表される。ここで、
 電子による散乱のためのフォノンの熱抵抗: Wge=E*T^-2
 境界散乱による熱抵抗: WB=B*T^-3
 点欠陥あるいは同位元素によるもの: Wp=P*T
 U過程: Wu=U*T^n*exp(-ΘD/(α*T))
 転位によるもの: WD=D*T^-2
のような温度依存性を持つ。
日本金属学会編「格子振動・非線形光学」丸善
-----------------------------------------------------------------------
■ 非調和結晶とフォノン
 これまでに述べてきたフォノンの描像は、原子変位の大きさが格子定数に比して十分小さいような格子振動に関するものであった。一方、非調和項の影響が決して小さい摂動と考えられないような場合、すなわち非調和性がむしろ本質的であるような現象については、格子振動に対してどのようなモデルを考えればよいか、またフォノンの概念を拡張してこれに適用しようとすると、どのような新しい描像が可能であるのかなどの興味ある問題が残されている。
 一般に結晶の相転移においては、温度(または圧力など)を変化させたとき、転移点を境にして原子配置に変化が生じ、結晶の対称性が変化する。このとき、ともの原子配置に固有なフォノンは、いつまでも安定なモードのままであるとはかぎらず、何らかの変化が生じることが期待される。
 特に興味があるのは、結晶の相転移に際して特定のフォノン・モードが不安定になる現象であり、これは実際以下に述べるように、ある種の結晶相転移において観測されている。
 ここではフォノン・モードの不安定化に起因する結晶の相転移を特に構造相転移とよぶことにし、構造相転移とフォノンとの関係について詳述する。
 元来、結晶がある圧力、温度の下で定まった構造をもつとき、イオンはその結晶構造で定まる格子点で最小のポテンシャル・エネルギーをもつ。それゆえイオンがその格子点からわずかに変位するとき、イオンをもとの格子点に引き戻そうとする復元力が常に働いており、結果として結晶固有の格子振動モードが存在することになる。
 重要なことは、この場合、格子の固有振動モードは結晶の空間対称性を反映していなければならず、群論の言葉でいえば、結晶の属する空間群またはその部分群の規約表現の基底ベクトルになっていなければならない。
 固有振動モード(あるいは同じことであるがフォノン)は波数ベクトルqと、偏りおよび分枝を区別する量子数jとで指定されるから、(q,j)の組は結晶の対称性で許されるものに制限されていることがわかる。
 いま結晶が別の構造に変化する場合を考えると、一般にはもとの格子点はイオンにとって安定なポテンシャル最小の点ではなくなり、これに伴ってもとの格子点から見た復元力は小さくなり、極端な場合は符号さえ変えることがありうる。そのような場合には、もとの結晶構造に固有な格子振動モードのうち、少なくとも一つは振動数がゼロに近づき、ついには虚数にさえなってしまう(力の定数が負の値になれば振動数は虚数となる)。
 実際、SrTiO3の高温相は立方晶であるが、温度を下げて行くと約105 Kで正方晶に相変化する。
 このとき波数ベクトル(q=1/2,1/2,1/2)を持ち、かつ規約表現Γ25の対称性を持つモードが不安定化し、その振動数は相転移点でゼロとなるような温度変化を示す。
 このように構造相転移の転移点に近づくにしたがって、振動数がゼロに近づいて行くフォノン・モードをソフト・モードと呼んでいる。
 ある特定の(q,j)で指定されるフォノン・モードがソフト化し、振動数ゼロの極限で結晶構造の相転移が起こることを概念的に理解するために、一次元格子の振動を考えてみる。
 格子定数dの一次元格子において、隣り合うイオンが互いに逆位相で振動する波数q=2π/2dの縦波フォノン・モードを表す。このような固有振動モードがソフト化すると、振動数が小さくなるにつれてその振幅[基準座標Q(ω(q))で表される]は増大し、同時にQの微分値は小さくなって、振動数ゼロの極限で凍結してしまう。
 結局、一次元格子は格子定数が2dとなり、単位胞内に2個のイオンを含む構造に変化したことになる。
 イオンの安定位置の変位をもたらすような結晶変態を変位型相転移とよぶこともある。
 このときは振動数ゼロの極限でこの変位を表す基準座標(厳密には基準座標の統計力学的平均値)をQ0とすると、転位点を境にして一方の温度領域でQ0=0,他の領域でQ0≠0であるから、Q0を秩序パラメタ―に選んで相転移を現象論的に扱うことができる。
 一般に結晶のHelmholtzの自由エネルギーをFとすると、これはイオンの変位ベクトルについて2次以上の関数になっている。このことを考慮すると、FをあらわにQ0の関数として表したとき、転移点近傍では次のようにQ0のべきに展開することができる。この表式が、
 F = F0 + (1/2)*A*Q0^2 + (1/4)*C*Q0^4 + ・・・ (C > 0)
のように与えられる場合をまず考える。これはSrTiO3のように反転対称を持つ結晶構造の場合に対応し、対称性からFはQ0の奇数べきの項を含まない。転移点よりも高温の相でQ0=0の状態が安定相であるためには、熱力学的安定条件
(∂F/∂Q0)= 0 (Q0=0) (一回微分)
(∂2F/∂Q02)> 0 (Q0=0) (2回微分)
から、この相でA>0であることが要求される。もしA<0となればQ0=0の状態はもはや安定相ではなくなり、安定な相ではQ0≠0となることが期待される。したがって相転移が起こるのは、Aが温度の関数として変化してちょうどA=0となるときであり、その温度が転移点である。
 転移点Tc以下の温度になると、Q0≠0であるような安定相が連続的に現れてくることがわかる。
 一方、FがQ0の奇数べきの項を含みうる場合には、Fの展開式に(1/3)*B*Q0^3の項を付け加えればよい。このときもしB<0でかつC>0であれば、自由エネルギーは温度の低下とともに下に凸の曲線を描くようになる。
 転移点Tcにおいては、有限のQ0をもつ状態が高温相(Q0=0の状態)と同じ自由エネルギーの値を持つことになり、したがって、そこで相転移が不連続に起こる。もちろん、T<Tcの温度ではFの極小値に対応するQ0≠0の状態が安定相である。そのような振る舞いをもつ場合は一次相転移であり、そうでない場合には二次相転移である。熱力学的にはGibbs自由エネルギーG(T,P)の一次微分が転移点で不連続な変化を示すもの、すなわち潜熱の伴うものを一次相転移と定義し、
二次微分が不連続な変化を示すもの、たとえば比熱の不連続変化が観測されるものを二次相転移と定義する。上に述べたSrTiO3の105Kでの構造相転移は、変位型相転移であり、かつ二次相転移である。
 場合によっては、結晶の対称性からFの展開式がQ0の奇数べきの項を含まないにも関わらず、一次相転移を示すような物質がある。この場合には、さらに6次の項である(1/6)*D*Q0^6を付け加え、C<0, D>0と仮定することによって、一次相転移を現象論的に説明することができる。
 またいくつかの強誘電体の構造相転移に見られるように、ソフト・モードに対応するQ0のほかに自発分極Pその他の秩序パラメーターを導入する方が記述に便利な場合がある。
 係数Aが転移点Tcで符号を変えるような温度依存性を持つというモデルは、ソフト・モードの振動数がTcでゼロとなることと密接に関係している。したがって構造相転移の構造を解明するためには、まず格子の不安定化の表れとして特定のフォノン・モードがソフト化する原因を探らなければならない。
 特定のフォノン・モードがソフト化する原因は結晶によって異なり、種々のモデルが考えられているが、大きく分けて電子状態の不安定性が引き金になるものと、格子振動の非調和項が主因と考えられるものがある。
 金属においては、電子の不安定化に伴って電子の運動がイオンの運動を完全に遮蔽してしまうためであると考えることができる。遷移金属カルコゲナイドのような層状物質で見られるように、波数空間でのFermi面の構造がかなりの面積にわたって互いに平行な面をもつような場合には、電子系は低温で新しい電子状態が安定になることが論じられる。
 絶縁体においては、伝導電子が存在しないから、フォノンのソフト化はフォノン同士の非調和相互作用によってのみ生じうる。
日本金属学会編「格子振動・非線形光学」丸善

 結晶ポテンシャルは、絶縁体の場合はイオン間の相互作用によるポテンシャルエネルギーであり、金属の場合には伝導電子によるポテンシャルがこれに加えられる。

 調和近似の範囲を超えて、イオンの動的ふるまいを正しく求めるためには、結晶ポテンシャルを構成するイオン自身が振動していることを反映して、これによる効果を結晶ポテンシャルの中に自己無撞着的に取り入れておかなければならない。
 この効果としては、一般に次の三つが考えられる。すなわち
1)伝導電子とフォノンの相互作用による媒質の動的偏極効果
2)着目するフォノンの振動数に依存した伝導電子の遮蔽効果
3)イオン同士の非調和相互作用による媒質の動的偏極効果
-----------------------------------------------------------------------
■ 分散関係の実験的解析
 フォノンの分散関係を知るために最も有効に利用される現象は、X線散漫散乱、熱中性子散乱、赤外吸収およびRaman散乱(Brillouin散乱を含む)などである。
 これらの現象を利用して得られる分散関係の情報は、フォノンの波数ベクトルの全領域にわたるわけではなく、またエネルギー領域も限られている。
 熱中性子散乱の実験で得られる情報は、フォノンのエネルギーおよび波数ベクトルの大きさについてそれぞれ99%の領域を網羅している。しかし熱中性子散乱では、エネルギーおよび波数ベクトルの測定精度が約1%であり、かつより小さい波数ベクトルのフォノンについては高い測定精度は望めない。一方、赤外吸収およびRaman散乱の実験は、小さい波数ベクトルのフォノンに関する情報を得るための有力な手段となっている。
 このように熱中性子散乱およびX線散漫散乱の実験と光学的な実験は、フォノンの分散関係に関する情報を得るために、互いに相補的な関係にあることがわかる。

◇ X線散漫散乱
 格子振動の振動数は、X線の振動数に比して無視できるほど小さいので、X線の結晶による散乱は、結晶内の電子電荷の瞬間的な空間分布で決まると考えて差し支えない。
 格子振動に伴ってゆらぐ結晶内電子の電荷分布を、空間的な場所の関数としてen(r)と書くことにする。入射X線は単色であるとし、その波数ベクトルをk0と書く。また空間電荷分布en(r)によってある方向に散乱されたX線に着目し、その波数ベクトルをkとする。上に述べたように空間電荷分布の時間依存性は無視できるので、この散乱は弾性散乱と見なしてよく、波数ベクトルの大きさはほとんど変わらない。すなわち
 |k0| = |k| = 2π/λ
である。λはX線の波長である。
 Q ≡ k - k0
を散乱ベクトルと呼ぶ。空間電荷密度en(r)を平面波に分解して考えたとき、散乱ベクトルQの波とちょうど同じ位相の平面波成分だけが散乱に寄与するから、k方向への散乱波の振幅は、数学的な言葉でいえばen(r)のフーリエ変換
 ∫[en(r)*exp(-iQ・r)]d^3r
に比例する。したがって、k方向への散乱強度は、
 I = I0 * C^2 *|<∫[n(r)*exp(-iQ・r)]d^3r>|^2
で与えられる。ここでI0は入射X線強度であり、C^2は、
C^2 = ( (1+cosΘ*cosΘ) / 2) * ( (e/c)^2 / (m*R) )^2
で与えられる。Rは結晶から観測点までの距離、Θはkとk0のなす角である。<>は、結晶の統計学的状態についての平均を表す。以上のようにX線の散乱強度の測定は、電荷密度の空間的なゆらぎのスペクトルを観測することに相当する。

 時間的な原子位置の変位を反映した強度をIとすると、変位を調和振動モードの1次結合と表せば、
 I = I0 + I1 + I2 + ...
と幾つかの項に分けた表現をすることができる。この第1項はQ=Kに相当する方向の回折像強度を与えるが、I1以下の項はこの回折像の広がりを与えるものであり、散漫散乱と呼ばれている。
 I1の項では、
Q ≡ k - k0 = K ± q
を満たす散乱方向に有限の散乱強度が見出せることを示している。±qはフォノンの波数ベクトルであり、この散乱に1個のフォノンの吸収または放出が関与していることを物語っている。
 同様にI2の項は、
 Q = K ± q1 ± q2
で与えられる散乱ベクトルの方向への散乱を表す。これは、q1、q2の二つのフォノンの関与する散乱過程によるものである。
 このように格子振動の影響は、Debye-Waller因子だけ散乱強度を弱めるとともに、回折像にはLaue条件Q = Kで与えられる散乱方向のまわりの広がり(ぼやけ)として現れる。これが散漫散乱(diffuse scattering)とよばれる理由であるが、特にI1を1次の散漫散乱、I2を二次の散漫散乱などと呼んで区別することもある。
 フォノンの分散関係を適当に仮定すれば下記の2つの式を使って、
  <|Q|^2> =( (h/2π)/ωj(q) )*(n(q,j) + 1/2)
  n(q,j) = 1/[exp( (h/2π)*ωj(q)/(kb*T) )-1]
<|Q|^2>を計算で求めることができるので、I1またはI2までの項だけで散漫散乱の強度が与えられると仮定すれば、実測される強度分布と理論から予測される分布を一致させることによって、フォノン分散関係あるいはフォノンの振動数分布を求めることができる。
 なお実験的には、測定された強度から、Compton散乱その他の寄与を注意深く補正しておかなければならないことはいうまでもない。なお実験の詳しい解析については、一例として文献[C. B. Walker: Phys. Rev. 103 (1956) 547.]を参照されたい。
 熱中性子の可干渉散乱断面積が小さいような物質では、X線散漫散乱の実験は、フォノンの分散関係を得るための有効な実験手段となっている。たとえば、Cdは熱中性子に対して強い吸収を示すので、熱中性子散乱ではフォノンの分散関係を得ることは困難であるが、X線散漫散乱の実験からフォノン分散関係が得られている。

◇ 熱中性子散乱
 エネルギーEをもつ中性子のde Brogile波長は、
 λ = h / √(2*mN*E)
で与えられる。ここでhはPlank定数、mNは中性子の質量である。試みに10meVのエネルギーをもった中性子を考えると、このエネルギーをE=kb*Tで温度にすると、T=120Kに相当するので熱中性子と呼ばれる。これは1.38km/sの速度をもち、de Brogile波長は2.86Åである。したがって、このような熱中性子はX線の場合と同様に、結晶解析のみならず、格子振動に関する情報を得るために用いられる。しかしながら、結晶による中性子散乱において特徴的なことは、熱中性子の振動数がフォノンのそれと同程度であり、そのために熱中性子の結晶による散乱においては、X線の場合とは異なって、弾性散乱のみならず非弾性散乱効果も顕著に現れることである。
 実はこの事実が、フォノンに関するより詳細な情報を直接得るための有力な手段を提供する。中性子は格子点近傍にある原子核とのに直接相互作用をして散乱される。位置r(lk)にある原子核と熱中性子との相互作用をふつう、ポテンシャル
 U(r) = (h/(2π*mN))*b(lk)*δ(r-r(lk))
で表す。これをFermiの擬ポテンシャルと呼んでる。ここでb(lk)は,lk番目イオンの原子核の散乱長(scattering length)と呼ばれ、X線散乱の場合の原子散乱因子f0に相当するものである。
 散乱の瞬間における原子核の位置r(lk)は、格子点R0(lk)からu(lk)だけ変位したところにあるとする。
 変位u(lk)の時間的な変化の速さは、上述したように熱中性子の振動数とほぼ同程度の大きさであるから、X線の場合とはことなってU(r,t)の時間依存性を無視することはできない。言い換えれば、X線散漫散乱においては電子電荷密度の空間的ゆらぎを観測するのに対して、熱中性子散乱においては、原子核の空間分布の時間的ならびに空間的ゆらぎを観測することができる。
 したがって熱中性子の入射波の波数ベクトルをk0としたとき、波数ベクトルkをもった散乱波の振幅は、X線の場合の拡張として、原子核の空間密度分布n(r,t)の時間、空間に関するフーリエ変換
 n(Q,ω)≡∫(∫[n(r,t)*exp(-iQ・r+iωt)]d^3r)dt
に比例することが期待される。ここでn(r,t)は原子核の散乱長b(lk)を含めて、
 n(r,t) = Σ[b(lk)*δ(r-R0(lk)-u(lk))]
とおかれる。また散乱ベクトルQはX線の場合と同じく
 Q ≡ k - k0
で与えられる。ωは散乱中性子と入射中性子のエネルギー差に相当する振動数で
 (h/2π)*ω ≡ (h/(2π*mN))*(k^2 - k0^2)
で与えられる。上記で与えられるn(Q,ω)、すなわち、原子核の空間的密度分布n(r,t)の時間的ならびに空間的ゆらぎのスペクトルを実験的に観測するためには、散乱中性子の散乱角と同時にエネルギーを測定しなければならない。このような測定で得られる量は、散乱波方向の単位立体角、単位エネルギーあたりの散乱微分断面積であり、この量は<|n(Q,ω)|^2>に比例する。<>の記号は結晶の状態についての統計学的平均を表す。
 中性子はスピン1/2をもつから、散乱体である原子核の核スピンSが0でないときは、両方のスピンが平行か反平行かで散乱長の値は異なる。これをそれぞれb+、b-と書くことにする。
 結晶中の原子核は、中性子との衝突の瞬間において、あるものは衝突する中性子と平行の核スピンをもち、あるものは反平行の核スピンをもっているから、これはちょうどb+、b-の2種類の散乱長を付与された原子核が、格子点に無秩序に配置されているとした場合と等価である。そのうえ一般に同位元素が結晶内に存在することをも考慮すると結晶内のすべての原子核が同じ散乱長bを付与されているとせず、異なる値の散乱長が結晶内の各原子核に存在することをも考慮すると、結晶内の全ての原子核が同じ散乱長bを付与されているとせず、異なる値の散乱長が結晶内の各原子核に無秩序に付与されていると考えなければならない。したがって、このような結晶から散乱される中性子線は、一部は可干渉的であり、残りは非干渉的になっている。前者を可干渉散乱(coherent scattering)、後者を非干渉散乱(incoherent scattering)と呼んで区別する。これに対応して、散乱微分断面積は、可干渉および非干渉散乱微分断面積の二つの項に分けることができる。
 結晶の状態について統計学的平均をとるときに、散乱長の平均は原子核のスピン状態について平均すれば十分である。可干渉散乱では、原子核のスピン状態についての平均をとるとき<b>av≠0であることから、可干渉散乱断面積は<b>av^2に比例する。
  変位に対して、調和振動近似を用いると、可干渉ならびに非干渉散乱微分断面積を見やすい形に書き換えることができる。そこで得られる可干渉散乱微分断面積は、中性子のエネルギーが変化しない可干渉弾性散乱とフォノンの関与する可干渉的な非弾性散乱で表現される。
 フォノンの分散関係を得るために重要なのは、可干渉的な非弾性散乱の部分となる。この中には、1個のフォノンの放出または吸収を伴った散乱過程に対応する部分が存在する。そこでの因子凵iQ±q)は
 Q ≡ k - k0 = ∓q + K (Kは逆格子ベクトル)
が成り立つ方向に可干渉散乱が生じることを物語っている。波数ベクトルに(h/2π)を掛けたものが運動量であるから、上記の1音子散乱過程での運動量保存則であると見なすことができる。また同項の因子
 δ(ω±ω(j,q))は、
エネルギー保存則、
 (h/2π)*ω ≡ (h*k/2π)^2/2mN - (h*k0/2π)^2/2mN = ∓(h/2π)*ω(j,q)
がこのときに同時に成り立っていなければならないことを示している。このように1個のフォノンの関与する可干渉非弾性散乱においては、上記の二つの条件が同時に成り立っていなければならないことから、散乱に関与するフォノンの分散関係を求めることが原理的に可能である。しかも同じ項の係数はX線散漫散乱の場合と同様、|(Q・e(j,q))|^2 なる因子を含んでいるので、散乱ベクトルQの方向と、フォノンの偏りベクトルe(j,q)が平行である場合にその係数の大きさは最大になる。このことから逆に、着目するフォノンの偏りを決定することができる。
 非干渉散乱微分断面積を与える式では、弾性散乱と非弾性散乱に項に分けられる。さらに1個のフォノンの放出または吸収を伴う項が存在する。この因子δ(ω±ω(j,q))は、エネルギー保存則が成り立つことを示してるが、運動量に対する保存則は存在しない。同項において、(q,j)についての和を積分に置き換えると、1音子の関与する非干渉散乱微分断面積は[1+振動数ωのフォノンの個数]*フォノンの振動数分布(調和振動モードの数の分布)の形が含まれている。非干渉非弾性散乱の測定から、フォノンの振動数分布を直接求めることができる。
 以上のように熱中性子の散乱は、フォノンの分散関係あるいは振動数分布を得るための直接的でかつ有力な手段を提供するが、物質によっては、非干渉散乱断面積が可干渉散乱断面積に比べて圧倒的に大きいものがあるので、いつもフォノン分散関係の測定が可能であるというわけではない。また中性子は原子核によって散乱されるばかりでなく、吸収されることもある。一般に原子核による中性子の全散乱断面積は、1〜100meVのエネルギーをもつ熱中性子に対しては波長によらないが、吸収断面積は波長に比例して増大する。吸収断面積の値は元素によってことなるが、吸収断面積が特に大きい値をもつものは、少数の元素の同位元素に限られているので、中性子散乱の実験は、大部分の物質についてフォノンに関する直接的な情報を得るための有力な手段であることには変わりはない。断面積の単位としてはふつうbarn(=10^-24cm2)が用いられる。
 熱中性子のエネルギー(または波長)を決定するためには、結晶によるBragg散乱を利用した分光法か、またはパルス状の中性子線が既知の距離を飛翔するのに要する時間を測定するいわゆるTOF(time of flight)の方法が用いられる。フォノン分散関係を決定するために、1音子の関与する熱中性子の可干渉散乱に着目する。この散乱過程では、波数ベクトルに対する保存則およびエネルギー保存則が同時に成り立たなければならない。
 入射波数ベクトルk0の方向に対して、試料の選ばれた結晶軸方向がなす角ψ(時計回りのとき正とする)、散乱波数ベクトルkがなす角をΦ(時計回りのとき正とする)ととる。運動量保存則を満たすため
 -k*sin(Φ-ψ)-k0*sin(ψ) = Qx = qx + Kx
 -k*cos(Φ-ψ)-k0*cos(ψ) = Qy = qy + Ky
またエネルギー保存則は、
 ( (h/2π)^2 * (k^2 - k0^2) )/2mN= ∓(h/2π)*ω(j,q)
で与えられる。一方、実験的に動かしうるパラメーターはk0, k, Φ, ψの四つである。したがって、そのうちの一つ例えばkの値を固定しておいて、上の3式を同時に満たすような他の三つのパラメーターの組を見出すことは可能である。この原理に基づいて工夫されたのが、Brockhouseらの考案による三軸結晶分光器(triple-axis crystal spectrometer)である。
単結晶MによるBragg反射によって単色化された熱中性子は、単結晶試料(ふつう数cm^3)Sによって
可干渉非弾性散乱を受けたのち、単結晶AによるBragg反射によって再び分光され、BF3または3Heを使った検出器Dで検出される。このとき独立に動かしうるパラメーターは、k0、あるいは同じことであるが入射エネルギーE0=((h*k0/2π)^2)/2mNを決めるための角度ΘM、散乱中性子エネルギーE'((h*k/2π)^2)/2mNを決めるための角度ΘA、および前に述べたΦとψを同時に動かす。
 Qを一定に選ぶことはフォノンのqを指定することであり、ΘAを変えたときの中性子のカウント数最大のところを見いだせば、E'したがってω(j,q)が決まる。保存則はE0とE'について対称であるから、E0を固定するかわりにE'を固定し、ΘMを変える方法も用いられる。このようにQ=k-k0すなわちqを一定に保ち、E'-E0を変化させて保存則を満たすようなE'-E0の値を見出す方法をconstant Q methodと呼んでる。
 フォノン分散関係でω(j,q)がqに対して鋭く立ち上がるような領域に対しては、Qを一定にするかわりに( (h/2π)^2 * (k^2 - k0^2) )/2mN を一定に保ち、保存則を満たすようなQを見出す方法を用いることもある。これをconstant ω methodと呼んでいる。
日本金属学会編「格子振動・非線形光学」丸善
-----------------------------------------------------------------------------
■ 計算条件
・PAW法 + LDA で 分極計算にベリー位相法、フォンノン分散の計算にフローズンフォノン法
■ 実験条件
・ ラマン分光法によるラマン散乱の観測(典型的なソフトモードが観測されている例がある[2])。転移温度前後でソフトモードの振動数が変化している。
■ 重要な解説
・通常の実数の振動数を有する格子振動の場合は原子変位に伴って復元力が働き、元の原子位置に戻す力が働くことにより、ある振動数で振動することになるが、虚数の振動数を有する振動の場合には原子変位に伴う復元力が働かず、原子が元の位置には戻らずに新しい位置に変位してしまうことを意味しており、ソフトモードフォノンと呼ばれる。
・結晶格子体積を数%(2%)変化させて、ソフトモードフォノンがフォノン分散曲線に現れるかを調べている。ソフトモードフォノンが存在するときには相転移が予想される。
・ソフトモードを有する物質にソフトモード有しない物質を置換固溶すれば、置換量に従って、ソフトモードが抑制され、強誘電性相転移が抑制される可能性がある。
[1] 森分 博紀ら、「第一原理計算と高精度実験の連携による強誘電体材料研究」、応用物理、81 (2012) 760.
[2] H. Taniguchi et al., Phys. Rev. B 76 (2007) 212103.
-----------------------------------------------------------------------------
◆ フォノン計算用のスクリプト
-------------------------------------------------------------------
□ phonopy-1.6.4 (VASP用スクリプト)
コメント:
  下記では ENCUTを 25 Ry (= 350 eV)程度にしてあります。INCARはINCAR_backupとして保存して、新しいINCARを作ります。IALGOはVASP.4.5, 4.6, 5.2 ではディフォルトで設定されますが念のために記しておきます。vasp.5.2はホームフォルダーの下のVASPという名称のフォルダで解凍しています。下記のスクリプトは単位胞の構造をPOSCARに入力しておくと、2x2x2のスーパーセルとなるようにしています。簡単なスクリプトですので、色々と作り変えてみてください。

#!/bin/bash

cp POSCAR POSCAR-unitcell
$HOME/phonopy-1.6.4/bin/phonopy -d --dim="2 2 2" -c POSCAR-unitcell
mv SPOSCAR POSCAR

mv INCAR INCAR_backup
echo PREC = Accurate >> INCAR
echo ENCUT = 350 >> INCAR
echo IBRION = 8 >> INCAR
echo IALGO = 38 >> INCAR
echo EDIFF = 1.0e-08 >> INCAR
echo ISMEAR = 0 >> INCAR
echo SIGMA = 0.1 >> INCAR
echo LREAL = .FALSE. >> INCAR
echo ADDGRID = .TRUE. >> INCAR
echo LWAVE = .FALSE. >> INCAR
echo LCHARG = .FALSE. >> INCAR

mv KPOINTS KPOINTS_backup
echo phonon calculation supercell 2x2x2 >> KPOINTS
echo  0  >> KPOINTS
echo Monkhorst Pack >> KPOINTS
echo 4 4 4 >> KPOINTS
echo 0 0 0 >> KPOINTS

mpirun -np 4 $HOME/VASP/vasp.5.2/vasp
$HOME/phonopy-1.6.4/bin/phonopy --fc vasprun.xml

$HOME/phonopy-1.6.4/bin/phonopy --dim="2 2 2" -c POSCAR-unitcell band.conf
$HOME/phonopy-1.6.4/bin/bandplot band.yaml
-------------------------------------------------------------------
□ band.conf
コメント:band.conf などでは、FORCE_CONSTANTS = READ を書き込みます。http://www.cryst.ehu.es/ でのKVECから各空間群に対するk点の名称と位置が掛かれます。Fm-3mですと下記はGamma-X-W-K-Gamma-Lになります。
ATOM_NAME = XX XX XX
DIM =  2 2 2
BAND = 0.0 0.0 0.0  0.5 0.5 0.0  0.5 0.75 0.25  0.375 0.75 0.375  0.0 0.0 0.0  0.5 0.5 0.5
FORCE_CONSTANTS = READ
-------------------------------------------------------------------
□ PWscf
DYNで発散したり、no convergence となる場合は、alpha_mix(1) = 0.7 を 0.1 にする。
nmix_ph や tr2_ph を変えたり、occupation = "fixed" にしたりすることも検討する。
nmix_ph で100を超えた値としたい場合はコードを書き換えてコンパイルし直します。
-----------------------------------------------------------------------------
アクセス数
ページビュー数