固体物理関係のグラフ（その1）

2020年度から固体物理の講義を持つようになったので、その授業のためにグラフをいくつか作りました。その備忘録も兼ねて、使ったスクリプトを掲載しておきます。

フェルミ分布関数

Fermi分布関数は \[ f(T,E) = \frac{1}{\exp[ (E - \mu) / k_{\mathrm{B}}T ] + 1} \] で与えられます。化学ポテンシャルを\(\mu/k_{\mathrm{B}}=10000~\mathrm{K}\)に固定し、様々な温度でのFermi分布関数をプロットするスクリプトです。

## HSV色空間からRGB色空間への変換。Wikipedia参照
#変換の補助関数群
HSV_Hi(h)=sgn(h)*floor(abs(h)/60)%6+(h>=0 ? 0 : 5)
HSV_f(h)=(sgn(h)*(abs(h)-(floor(abs(h))/360)*360.0)/60.0+(h>=0 ? 0 : 6))-HSV_Hi(h)
HSV2R(h,s,v)=( \
     HSV_Hi(h)==0 || HSV_Hi(h)==5 ? v                    : (\
     HSV_Hi(h)==1                 ? v*(1-HSV_f(h)*s)     : (\
     HSV_Hi(h)==4                 ? v*(1-(1-HSV_f(h))*s) : (\
     HSV_Hi(h)==2 || HSV_Hi(h)==3 ? v*(1-s)              : 0.0) )))
HSV2G(h,s,v)=( \
     HSV_Hi(h)==1 || HSV_Hi(h)==2 ? v                    : (\
     HSV_Hi(h)==3                 ? v*(1-HSV_f(h)*s)     : (\
     HSV_Hi(h)==0                 ? v*(1-(1-HSV_f(h))*s) : (\
     HSV_Hi(h)==4 || HSV_Hi(h)==5 ? v*(1-s)              : 0.0) )))
HSV2B(h,s,v)=( \
     HSV_Hi(h)==3 || HSV_Hi(h)==4 ? v                    : (\
     HSV_Hi(h)==5                 ? v*(1-HSV_f(h)*s)     : (\
     HSV_Hi(h)==2                 ? v*(1-(1-HSV_f(h))*s) : (\
     HSV_Hi(h)==0 || HSV_Hi(h)==1 ? v*(1-s)              : 0.0) )))
     
#HSV色空間のパラメータh,s,vをRGBの色を表す16進数文字列"#rrggbb"に変換する。
#3番目の引数の係数が200になっているのは、Greenを255にすると明るすぎて紙面上や画面上で見づらいため。
HSV2RGB(h,s,v)=sprintf("#%02x%02x%02x", 255*HSV2R(h,s,v), 200*HSV2G(h,s,v), 255*HSV2B(h,s,v))

#整数nをRGB文字列"#rrggbb"に変換する。
period=-6.0   #グラデーションの周期（360で色が一周する）
start=270.0     #グラデーションをスタートする色。0ならn=0は赤に対応する。
num2color(n)=HSV2RGB(360.0*n/period+start, 1.0, 1.0)

set grid
set term wxt font "Verdana,12"

set xlabel "{/:Italic E} / {/:Italic k}_B (K)" font ",15"
set ylabel "{/:Italic f}({/:Italic E})" font ",15"

# Chemical Potential (in Kelvin)
mu = 10000.0

# Temperatures to be plotted
array T[5] = [10, 60, 300, 2000, 10000]

# plot (1)
set xrange [0:20000]
set samples 1000
plot for [i=1:5] (1.0/(exp( (x-mu)/T[i] ) + 1)) w l lw 2 lc rgb num2color(i-1)  title sprintf("T = %.0f K", T[i])

pause -1

# plot (2)
set logscale y
set yrange [0.001:10]
set xrange [0:40000]
plot for [i=1:5] (1.0/(exp( (x-mu)/T[i] ) + 1)) w l lw 2 lc rgb num2color(i-1)  title sprintf("T = %.0f K", T[i]),\
     for [i=1:5] (exp(-(x-mu)/T[i])) w l lw 2 lc rgb num2color(i-1) dt "."  notitle

pause -1

# plot (3)
T=800.0
unset logscale xy
set xrange [0:20000]
set yrange [0:1]
set object 1 rect from first mu-2*T, graph 0 to first mu+2*T, graph 1 fc rgb "green" fillstyle solid 0.1 noborder
plot (1.0/(exp( (x-mu)/T ) + 1)) w l lw 2 lc rgb "green"  title sprintf("{/:Italic T} = %.0f K", T)

これの実行例は以下になります。温度上昇によって、Fermi分布関数がほぼ階段関数の形状から、なだらかな変化の関数へと変わっていくことがわかります。化学ポテンシャルと同程度の温度になると、もはや階段関数とは言えなくなります。（ただし、\(\mu\)の温度変化が考慮されていないことに注意; 化学ポテンシャルと同程度の温度では、化学ポテンシャル自体が元のFermiエネルギーから大きくずれるはずです。）

OKをクリックすると、以下のような縦軸をLogスケールにしたグラフができます。化学ポテンシャル\(\mu\)より十分高エネルギー側では、Fermi分布がボルツマン分布\(f = \exp[-(E-\mu)/k_{\mathrm{B}}T]\)に漸近していくことがわかります。

さらにOKをクリックすると、以下のようなグラフができます。Fermi分布関数は化学ポテンシャル\(\mu\)の周辺\(\mu\pm 2k_{\mathrm{B}}T\)の範囲で大きく変化しますが、その範囲を薄い背景色で表現しました。

フェルミ分布関数 その2

multiplotを用いて、いくつかの温度のFermi分布関数を縦に並べて表示します。

## HSV色空間からRGB色空間への変換。Wikipedia参照
#変換の補助関数群
HSV_Hi(h)=sgn(h)*floor(abs(h)/60)%6+(h>=0 ? 0 : 5)
HSV_f(h)=(sgn(h)*(abs(h)-(floor(abs(h))/360)*360.0)/60.0+(h>=0 ? 0 : 6))-HSV_Hi(h)
HSV2R(h,s,v)=( \
     HSV_Hi(h)==0 || HSV_Hi(h)==5 ? v                    : (\
     HSV_Hi(h)==1                 ? v*(1-HSV_f(h)*s)     : (\
     HSV_Hi(h)==4                 ? v*(1-(1-HSV_f(h))*s) : (\
     HSV_Hi(h)==2 || HSV_Hi(h)==3 ? v*(1-s)              : 0.0) )))
HSV2G(h,s,v)=( \
     HSV_Hi(h)==1 || HSV_Hi(h)==2 ? v                    : (\
     HSV_Hi(h)==3                 ? v*(1-HSV_f(h)*s)     : (\
     HSV_Hi(h)==0                 ? v*(1-(1-HSV_f(h))*s) : (\
     HSV_Hi(h)==4 || HSV_Hi(h)==5 ? v*(1-s)              : 0.0) )))
HSV2B(h,s,v)=( \
     HSV_Hi(h)==3 || HSV_Hi(h)==4 ? v                    : (\
     HSV_Hi(h)==5                 ? v*(1-HSV_f(h)*s)     : (\
     HSV_Hi(h)==2                 ? v*(1-(1-HSV_f(h))*s) : (\
     HSV_Hi(h)==0 || HSV_Hi(h)==1 ? v*(1-s)              : 0.0) )))
     
#HSV色空間のパラメータh,s,vをRGBの色を表す16進数文字列"#rrggbb"に変換する。
HSV2RGB(h,s,v)=sprintf("#%02x%02x%02x", 255*HSV2R(h,s,v), 255*HSV2G(h,s,v), 255*HSV2B(h,s,v))

#整数nをRGB文字列"#rrggbb"に変換する。
period=-4.3  #グラデーションの周期（360で色が一周する）
start=270.0     #グラデーションをスタートする色。0ならn=0は赤に対応する。
num2color(n)=HSV2RGB(360.0*n/period+start, 1.0, 1.0)
num2color_thin(n)=HSV2RGB(360.0*n/period+start, 0.2, 1)

set grid
set term wxt size 600,800 font "Verdana,12"

set xlabel "{/:Italic E} / {/:Italic k}_B (K)" font ",15"
set ylabel "{/:Italic f}"

# Chemical potential (in Kelvin)
mu = 10000.

array T[4] = [50, 200, 800, 3200]

set multiplot layout 4,1 downwards margins 0.15, 0.75, 0.15, 0.95 spacing 0.0, 0.05

set xrange [0:20000]
set yrange [0:1.]
set samples 1000
set arrow 10 from first mu, 0 to first mu, 1  nohead lw 1 dt "." front

set key outside samplen 0.5
do for [i=1:4] {
     unset xlabel
     set format x ""
     set object 1 rect from first mu-2*T[i], graph 0 to first mu+2*T[i], graph 1 fc rgb num2color_thin(i-1) fillstyle solid 1 noborder
     if (i==4){
	  set xlabel "{/:Italic E} / {/:Italic k}_B (K)" font ",15"
	  set format x
	  }
     plot (1.0/(exp( (x-mu)/T[i] ) + 1)) w l lw 2 lc rgb num2color(i-1)  title sprintf("{/:Italic T} = %.0f K", T[i]),\
	  (exp(-(x-mu)/T[i])) w l lw 2 lc rgb num2color(i-1) dt (2,2)  notitle
}

unset multiplot

実行例は以下の通り。一つ前の例と同じく、薄く背景色がついているのは化学ポテンシャル\(\mu\)の周辺\(\pm 2k_{\mathrm{B}}T\)の範囲です。この範囲内で分布関数が大きく変化することがわかります。また、点線はボルツマン分布\(f = \exp[-(E-\mu)/k_{\mathrm{B}}T]\)です。色のついた範囲よりも高エネルギー側ではボルツマン分布が良い近似になっていることがわかります。

状態密度とフェルミ分布関数

三次元電子系の状態密度は\(D(E) \propto \sqrt{E}\)で与えられます。それとFermi分布関数を掛け合わせると、有限温度の下での占有状態数の変化がわかります。それを模式的に表した図を作ります。（化学ポテンシャルの温度変化は考慮に入れていない）

set border 0

unset xtics
unset ytics
set format xy ""

set xlabel "Energy" font ",15"
set ylabel "Density of states" font ",15"

set term wxt size 800,400
set xrange [0:25]
set yrange [0:5]
set arrow 1 from graph 0, first 0 to graph 1, first 0 head filled lw 1 front
set arrow 2 from graph 0.,0 to graph 0.,1 head filled lw 1 front

x0 = 4.4
EF = 16
f(T, E, mu) = 1.0 / ( exp( ( E - mu ) / T ) + 1.0)

set arrow 10 from first EF, 0 to first EF, sqrt(EF) nohead lw 1 dt "." front
set samples 1000

plot  sqrt(x)*f(0.3, x, EF) w filledcurves y = 0 fs solid 0.5 fc "cyan" notitle,\
      sqrt(x)*f(0.3, x, EF) w l lc rgb "blue" lw 1 notitle,\
      sqrt(x) w l lw 2 lc 1 notitle ,\

pause -1

plot  sqrt(x)*f(0.1, x, EF) w l lc rgb "blue" lw 1 notitle,\
      sqrt(x)*f(0.2, x, EF) w l lc rgb "green" lw 1 notitle,\
      sqrt(x)*f(0.3, x, EF) w l lc rgb "orange" lw 1 notitle,\
      sqrt(x)*f(0.4, x, EF) w l lc rgb "red" lw 1 notitle,\
      sqrt(x) w l lw 2 lc 1 notitle

実行例は以下の通りです。有限温度では、化学ポテンシャル以上の状態が占有されたり、化学ポテンシャル以下の状態が非占有になったりすることがわかります。

OKをクリックすると、下のようにいくつかの温度での分布を重ねて示した図になります。

化学ポテンシャルの温度変化

化学ポテンシャルは、多くの場合わずかに温度変化します。3次元自由電子モデルの場合、その温度変化は、ゾンマーフェルト展開を用いて、 \[ \mu = E_{\mathrm{F}}^0 -\frac{\pi^2}{12 E_{\mathrm{F}}^0}k_{\mathrm{B}}^2T^2 \] と求められます。 この式を\(E_{\mathrm{F}}^0/k_{\mathrm{B}} = 10000~\mathrm{K}\)の場合にプロットすると以下のようになり、室温程度の温度範囲内では、\(\mu\)の温度変化はわずかであることがわかります。

この図は以下のスクリプトで作ったものです:

c_  = 2.99792458e+8  # [m/s]
me_ = 9.1093897e-31  # [kg]
e_  = 1.60217733e-19 # [C]
mu_  = pi*4.0e-7     # [N/A^2]
eps_ = 1.0/mu_/c_**2 # [N*m^2/C^2]
h_ = 6.6260755e-34   # [J*sec]
hbar_= h_/(2.0*pi)
Na_=6.0221367e+23    # [mol^(-1)]
kB_=1.380662e-23     # [J/K]
Phi0_=h_/(2*e_)

aB_ = (4*pi*hbar_**2*eps_ / (me_ * e_**2))

EF = 10000*kB_

mu(T) = EF/kB_ - pi**2 / (12*EF) * kB_ * T**2

set encoding utf8
set term wxt size 600,400 font "Verdana,12"

set xlabel "{/:Italic T} (K)" font ",15"
set ylabel "{/:Italic μ}/{/:Italic k}_B (K)" font ",15"
set grid
set xrange [0:500]
plot mu(x) w l lc rgb "dark-red" lw 2 not