光學三層膜反射率曲線推導－早安，苦命工程師的胡言亂語

前言

這個部落格當初的規劃，有希望放一些物理心得筆記的。但是因為忙，大概都沒有時間寫！現在比較閒一點，來個第一篇吧。來一篇一個光學上，三層薄膜的光學反射率曲線公式的推導吧。最早是因為有這方面的需求，所以翻了之前大學三年級所學的光學講義，然後再延伸所推導出來的。這個公式有實證過，所以應該沒有什麼問題。因為非常的有用，所以來經驗分享一下。假如你正在學習光學的話，這篇blog會是不錯的參考。後面有利用 R 語言，使用這個物理公式，來展示一下單層光學薄膜的反射率曲線，更進階的範例就放下一篇 blog 囉！

這個推導的過程，其實也是物理系學生對馬克士威爾方程式在光學上的應用的很標準推導方式啦！剛好可以看到這組包含一切大自然電磁現象的神的方程式，馬克士威爾方程式的神奇之處囉。

馬克士威爾方程式

這個鼎鼎大名的方程式真的是太重要囉，自然界所有的電磁現象都在這個方程式中。所以，它對人類的影響是巨大的。所有的光學現象，例如相機，照相...。所有的電磁波現象，例如你的 iPhone 可以跟基地台連線，AM/FM 收音機... 也都來自於這個方程式。電路學，輸電系統，電動馬達，發電機說穿了也是這個方程式來的喲。沒了這個方程式，恐怕大家現在還在蠻荒時代吧！ (愛因斯坦還利用它，最後思考出狹義相對論所規範出來的另一個我們自然界所存在的不可思議的現象，據此建構了全新的時間空間的架構。)

這四個非常優美如神般的方程式如下 (取自維基百科)

這四個方程式，要回到當時的歷史。對於當時電磁學的各種現象，庫倫已經發現的庫侖定律所描述的靜電現象。法拉第也發現了電磁感應的現象，磁移動生電，電生磁，還有就是重要的電場跟磁場的創見。再來不同的人眾多實驗也發現了穩定電流跟磁場的關係！

這麼多豐碩的成果，因為當時大家都是實驗物理學家，數學不太好啦！所以非常的紛亂，見霧是霧，看花是花，無法統一成一個漂亮美麗的理論。一直到馬克士威爾，他是個精通數學的數學家，最後把這些實驗出來的定律完全以法拉第的物理直覺，電場磁場的「場」的概念跟簡潔的向量分析數學給統合起來，成為四個漂亮的方程式（最後一個方程式，馬克士威爾在安培定律補上一項他的創見，位移電流)。

這裡有很棒的介紹，有興趣的別錯過了！

這樣的統合，成果也是豐碩的，馬克士威爾率先從他的方程式中預測出一種未知的電磁波動的存在，並率先計算出這種波動移動的速度就是光速。所以馬克士威爾認為「光」就是他方程式裡所預測出來的這種電磁波動，光就是這種未知的電磁波動的其中一種形式。

現在，我們就來操作這神一般的馬克士威爾方程式，來推導出我們要的，光學三層膜下的反射率曲線吧！

平面電磁波

來介紹一下平面電磁波的數學形式。

馬克士威爾的這四個方程式，經過化簡合併，最後都可以得到電場E 跟磁場H的二階偏微分波動方程式，物理學家或數學家對這個二階偏微分波動方程式是非常熟悉的，自然就預測出電場E 跟磁場H 會像水波或繩波一樣，會以波動的方式在空間中傳遞出去，形成所謂的電磁波。

弦波也是這種類型二階偏微分波動方程式的基本解的，所以以弦波為基礎的平面電磁波也會是這個二階偏微分波動方程式的解。

底下是平面電磁波的示意圖

一個平面法線向量 k 來指示這些平面移動的方向，跟這個向量k 垂直的是一層又一層的等位面(所以叫平面波)。等位面隨時間的改變沿著向量 k 所指示的方向移動著！

等位面的數學方程式為 f(x,y,z) = const, 讓 r = (x, y, z) 為一個隨意空間向量, 法向量 k, 所以等位面公式為

今天讓這個等位面沿著 k向量的方向隨時間移動, 所以隨時間移動的等位面公式為

所以隨時間移動的平面弦波可以寫成這樣

其中,

這個就是組成更複雜平面電磁波的基礎囉！

k 向量不只是指示方向而已，它也是弦波波動裡面波數的三維向量版本，omiga 是弦波波動裡的角速度參數，一個波動的速度u 由由角速度跟波數(或是頻率跟波長)來決定。所以這個弦波平面波的速度u 滿足下面這個式子

因為 cos(x) 不太容易運算，數學上我們喜歡用尤拉公式 exp(ix) = cos(x) + i sin(x) 將 cos(x) 換成有複數的 exp(ix)。因為 exp(ix) 微分後還是它自己乘上裡面提出來，這樣會比用一堆 sin, cos 來化簡簡單許多。最後只要捨棄複數部分，取實數部分就是我們要的！

所以最後電場 E跟磁場H 平面電磁波的形式如下

要注意的，這裡是最廣義的形式了， E0 是個複數向量，除了提供震盪的方向還提供相位差，所以可以組成各種的偏振波。(直線, 橢圓, 園... 偏振, 這裡懶得展示了，直接寫結果!)

同樣的，H0 也是個複數向量，除了提供震盪的方向還提供相位差，所以可以組成各種的偏振波。

滿足平面電磁波的馬克士威爾方程式

我們的平面電磁波是滿足馬克士威爾方程式，在這個條件下，馬克士威爾方程式可以改寫成很簡單的形式

我們發現，把平面電磁波對時間微分會變成這樣

把平面電磁波對空間微分會變成這樣

所以滿足平面電磁波的馬克士威爾方程式就可以很簡單的把旋度算符跟對時間的微分替換掉，形成下面四個很容易操作的形式：

(假定空間中無自由電荷及自由電流流動，所以這兩項為零喲!)

其中

B 跟 H間的係數是磁導率，假如電磁波不是在磁性物質的話，這個磁導率會跟在真空中一樣。(所以最後可以從方程式中消去)

D 跟 E間的係數是電介質的電容率，在不同的物質下，例如薄膜，這個值都會不太一樣。(習慣把它納入折射率中，方便計算)

滿足馬克士威爾方程式電磁波的特性

根據這四組平面電磁波的馬克士威爾方程式，可以告訴我們什麼呢?

這四個方程式，最重要應該是第二組法拉第定律跟第四組安培-馬克士威爾定律囉，因為這兩組是共軛的！

要同時滿足這兩組的外積向量的方向，k X E = H, k X H = -E，發現只有下面這種可能性！ k, E, H 兩兩互相垂直, 方向如下圖所示！(下面這張圖假定E跟H不會旋轉，是線性偏振的電磁波，事實上是也是可以隨時間改變轉動的，只要旋轉時依舊兩兩相互垂直，並不會違反馬克士威爾方程式！)

行進方向 k, 和電場波 E, 跟磁場波 H 兩倆相互垂直, 且 k 的方向等於 EXH (記住, 外積方向採用右手原則)

電磁波的能量等於 S= EXH, 也被稱為坡印廷向量. 電磁波能量流的流動由此定義.

因為這 k, E, H 兩兩相互垂直，所以它們的向量外積 sin(90) = 1, 是可以直接拆掉，變成純量方程式。

將第二式法拉第定律，跟第四式安培-馬克士威爾定律，這兩個物理來源截然不同的方程式彼此共軛。拆掉向量外積，相互帶入後可以得到

化簡後可以得到這電磁波動的速度u, 發現只跟磁性物質的磁導率和電介質的電容率有關。而真空下，電磁波的速度為光速c

上面的式子還有很重要的結論，電磁波(光)的速度只跟介質有關，只要是相同的介質下，速度就會一樣。另外磁性物質的磁導率和電介質的電容率在介質中會比在真空中還高，所以介質中電磁波(光)的速度會比光速還慢！介質的磁導率與電容率愈高，速度愈慢！大部分的物質都不是磁性物質，磁導率會跟在真空中一樣，所以介質中電磁波(光)的速度主要由電容率來決定。

接下來很重要的，我們知道電磁波動裡電場E, 會跟磁場H 相互垂直。那兩個之間的大小關係會是怎樣呢？

根據第四式，安培法拉第定律，去掉向量外積跟方向符號後可以得到如下的式子

再來假定我們的電磁波只會在電介質上運行，不會在磁導性物質上傳遞，這樣可以把折射率 n 給關連進來 (n=c/u)，同時定義一個常數阻抗 Z0

所以整個E 跟H 的大小關係式如下

很重要跟簡單的關係， H = n/Z0 E

也就是說，它們之間的大小比例，跟所在介質的折射率成正比。

基本關係式差不多介紹完畢了，可以開始我們的多層膜反射率曲線的推導囉！

單層膜反射率公式推導

要推導多層膜前，先來推導單層膜吧！因為多層膜公式是從單層膜來的！

示意圖如下

假設一開始光的介質，折射率為 n0, 薄膜折射率 n, 厚度 d, 最後面的物質折射率為 nt

入射光垂直薄膜，只考慮不旋轉直線偏振光，因為入射光垂直薄膜，所有反射光都是同一個入射光所分出來的，所以不用考慮什麼 TE波，TM波啦！

假設ㄧ開始入射光 E 平行紙面朝上，H 穿出紙面朝上，k垂直薄膜，所以所有的入射/反射/穿透光如下圖所示

電磁波動函數是需要有連續性的，不可以斷掉。所以在界面的地方，譬如說在 n0 的介質，瞬間進入薄膜 n 的介質的時候，在這個介面上，這兩個波動函數的數值要一樣。

所以根據介面邊界條件的連續性:

在第一個介面:

對E,

對H,

利用折射率關係式，將H 化成 E

在第二個介面:

對E,

對H,

利用折射率關係式，將H 化成 E

所以我們得到四組 E 的關係式，我們對薄膜中的電場強度 E1 不感興趣，消去它最後化簡得到底下兩組方程式

得到我們要的方程式了，這個方程式描述了反射率 r= E'0/E0 跟穿透率 t= ET/E0 的關係，把參數帶入直接就可以解出單層薄膜的反射率r 跟穿透率t。

多層膜的公式

但我們的興趣是在多層膜下的反射率。而剛剛那組方程式就是重要的起點。

有沒有注意到，因為對稱性，所以多層膜素乎只是單層膜的推廣。我們只要用遞歸的方式不斷的套用上面的方程式，就可以得到多層膜的方程式囉！

把上面的方程式改成矩陣的形式

因為對稱性，光每進入每一層薄膜，所遭受的的交互作用是一樣的，只差在不同的折射率和厚度。所以每一層可以看作是個線性算符。這個算符的矩陣跟關係如下

多層膜前後反射率關係式

多層膜關係式的最後主矩陣，假定矩陣係數如下

帶入剛剛的多層膜前後反射率關係式中，可以得到反射率 r 跟穿透率 t 的聯立方程組的兩個方程式。所以可以解出r 跟 t 如下

注意， r 跟 t 是複數。真實的反射率R 跟穿透率T 等於這個複數的絕對值大小。

三層薄膜反射率的推導

有了主矩陣跟所有對應的公式了，讓我們來推導三層薄膜的反射率曲線函數吧！

三層薄膜符號跟示意圖如下，光由折射率n0的介質入射三層薄膜, 每層薄膜依序為 (折射率 n1, 厚度d1), (折射率 n2, 厚度d2), (折射率 n3, 厚度d3), 最後抵達折射率 nt 的介質.

所以此三層薄膜的主矩陣如下,

呵呵！看起來很複雜對不對?

但對物理系的學生而言，這樣的符號運算真的是一塊小蛋糕！所以三個矩陣乘開後化簡，最後結果 A, B, C, D 如下 (有興趣的把它當作線上遊戲打副本，玩看看！)

有注意到 B, C 是純虛數，所以為了方便電腦運算，重新定義一下 a, b, c, d

所以反射率公式 r 可以化簡為

所以最後的反射率公式為

其中

這是完整三層薄膜的反射率曲線公式

有了上面這個公式，我們就可以開始計算各種薄膜組合的光學特性囉！

R 語言

來把上面的公式利用 R 語言給實作出來吧！上面這個公式還蠻複雜的，所以 Excel 非常的不適合啦. 用交談式的程式語言來實作會比較方便！

說來也有點巧遇，大概在 2006年的時候，筆者剛加入一個新公司。新公司畢露襤褸的，一開始連個像樣的統計軟體都沒有。可是報告還是得做呀，所以最後上網找了一下，發現了R 這個完全開源免費的統計語言。稍微地學了一下，就開始使用了。(蠻高階的，不難學！) 壓軸沒想到10年後，這個語言現在如此火紅，變成大數據處理專用的語言惹！

R 很高階，有很好的交談性，有很強的內建數學函數跟統計圖表的功能，所以非常適合用在這個三層薄膜反射率曲線的公式上。

程式很直接啦，就根據上面我們推導的公式，直接定義 aa(), bb(), cc(), dd() 四個函數來計算公式裡面 a, b, c, d 的四個值，最後利用ref 將 a, b, c, d 的四個值進一步統合成最後的反射率的數值。

要注意的細節，我們的 k 是波數，我們真實世界給的是光的波長，所以要先轉換一下波數 k = 2*pi / 波長 len

另外一個細節是當我們談到波長，是以我們在空氣中 n0=1 的觀察者來說的，可是介質裡面的波長會因為光移動速度比較慢而縮減，這個縮減率就是折射率。所以公式裡面的膜厚 d 需要因為這個波長的縮減修正，所以公式裡真實的膜厚d1 會是 d1 * n1/n0, d2 會是 d2 * n2/n0, d3 會是 d3 * n3/n0 。

還一個細節是粗糙度，公式是完美無粗糙的表面，可是真實的世界沒這麼完美。假如是非光滑平面的話，就有機會造成多重反射，光有機會再被反射一次。這裡是很約略的引進一個 rough 參數，讓反射率以 reflection^rough 次方的計算來估計這樣的效果。實際測試對粗糙表面的矽晶圓而言， rough 的數值大概在 1.1 ~ 1.5 之間。

我們先來單層膜的公式驗證吧！

這個公式相當地戎長，我們怎麼確認它是正確的呢？方法很簡單，我們可以利用一個單層膜來驗證。選定 n=1.5 / d= 85nm 的一個單層薄膜！

假如我們故意用

n1 = 1.5, d1 = 85 nm

n2 = nt, d2 = 0 nm

n3 = nt, d3 = 0 nm

假如結果跟用

n1 = n0, d1 = 0 nm

n2 = 1.5, d2 = 85 nm

n3 = nt, d3 = 0 nm

及跟用

n1 = n0, d1 = 0 nm

n2 = n0, d2 = 0 nm

n3 = 1.5, d3 = 85 nm

這三個結果都會一致的話，那我們的公式就是正確無誤的囉。因為 1.5 / 85nm 不管填入哪個位置，都得到正確的結果！

底下是程式跑出來的結果，結果正確無誤！以上面所設定的參數所計算出來的三條反射率曲線，完美的疊合在一起囉！

再稍微解說一下這個 n=1.5, d=85 nm 的結果

這個程式會自動標註波長 400nm - 700nm 這個可見光的區域，然後根據紅澄黃綠藍靛紫各單色光的波長，把每個顏色色光的位置標註上去。這樣就可以很容易判斷每一種色光反射的程度，進而判斷出肉眼所見可能的顏色呈現。

像上圖，反射最強的是紫色，其次是紅色。所以這個 1.5 / 85nm 的薄膜，肉眼會感覺是紫紅色的！

根據薄膜破壞性干涉的公式，當 n*d = 波長 / 4 時會有最強的破壞性干涉，這時候反射率會有最小值！所以這個 1.5 / 85nm 的薄膜，反射率最小值的位置其波長會是 4 n*d = 4*1.5*85 = 510nm, 這完全符合我們公式所畫出來的結果！所以所有的驗證，公式正確無誤！

剩下更多範例的解說，就留到下一篇 blog 囉！

(下一篇 blog, 光學三層膜反射率曲線公式應用篇，請點入閱讀)

2019.3.23 補充:

因為有網友留言，希望知道單層光學抗反射薄膜時，最佳薄膜折射率 n = sqrt(n0*nt) 這個公式怎麼來的？

我把推導過程放到這篇 BLOG 囉！(按我進入單層看反射膜反射率極值公式推導)

給有興趣的人參考！

R 語言原始程式列表：

# triple layers thin films

# 2010.1.28 Frank Lin

# reflection curve for 3 layers thin-film system

# 2018.12.13 - some small modifications for the use of my blog

# formula verification

aa=function(k,n1,n2,n3,d1,d2,d3)

{

aa=(cos(k*d1)*cos(k*d2)-(n2/n1)*sin(k*d1)*sin(k*d2))*cos(k*d3)-(cos(k*d1)*sin(k*d2)/n2+sin(k*d1)*cos(k*d2)/n1)*n3*sin(k*d3)

return(aa)

}

bb=function(k,n1,n2,n3,d1,d2,d3)

{

bb=-(cos(k*d1)*cos(k*d2)-(n2/n1)*sin(k*d1)*sin(k*d2))*sin(k*d3)/n3-(cos(k*d1)*sin(k*d2)/n2+sin(k*d1)*cos(k*d2)/n1)*cos(k*d3)

return(bb)

}

cc=function(k,n1,n2,n3,d1,d2,d3)

{

cc=-(n1*sin(k*d1)*cos(k*d2)+n2*cos(k*d1)*sin(k*d2))*cos(k*d3)-(-n1/n2*sin(k*d1)*sin(k*d2)+cos(k*d1)*cos(k*d2))*n3*sin(k*d3)

return(cc)

}

dd=function(k,n1,n2,n3,d1,d2,d3)

{

dd=-(n1*sin(k*d1)*cos(k*d2)+n2*cos(k*d1)*sin(k*d2))*sin(k*d3)/n3+(-n1/n2*sin(k*d1)*sin(k*d2)+cos(k*d1)*cos(k*d2))*cos(k*d3)

return(dd)

}

reflection=function(len,n1,d1,n2,d2,n3,d3,n0,nt,rough)

{

k=2*pi/len

d1=d1*n1/n0

d2=d2*n2/n0

d3=d3*n3/n0

a=aa(k,n1,n2,n3,d1,d2,d3)

b=bb(k,n1,n2,n3,d1,d2,d3)

c=cc(k,n1,n2,n3,d1,d2,d3)

d=dd(k,n1,n2,n3,d1,d2,d3)

ref=((a*n0-d*nt)^2+(b*nt*n0-c)^2)/((a*n0+d*nt)^2+(b*nt*n0+c)^2)

return(ref^rough)

}

NN = 500

xmin=300

xmax=1000

ymin=0

ymax=0.3

# rough is a factor to represent the roughness level of the surface

# rough=1 is perfect smooth surface. for Si wafers, it's experimentally about 1.1 ~ 1.5

rough=1

n0=1

nt=3.4

# curve 1

n1=1.5

d1=85

n2=nt

d2=0

n3=nt

d3=0

x=seq(xmin,xmax, by=(xmax-xmin)/NN)

y=reflection(x,n1,d1,n2,d2,n3,d3,n0,nt,rough)

par(mfrow=c(1,1))

plot(y ~ x,type="l", col="yellow", lwd=8, xlim=c(xmin,xmax), ylim=c(ymin,ymax), main="Reflection Curve for formula verification",xlab=c("Wave Length(nm)"), ylab=c(" Reflection"))

# curve 2

n1=n0

d1=0

n2=1.5

d2=85

n3=nt

d3=0

x=seq(xmin,xmax, by=(xmax-xmin)/NN)

y=reflection(x,n1,d1,n2,d2,n3,d3,n0,nt,rough)

i <- seq(NN-1)

segments(x[i],y[i], x[i+1], y[i+1], col="green", lwd=4)

# curve 3

n1=n0

d1=0

n2=n0

d2=0

n3=1.5

d3=85

x=seq(xmin,xmax, by=(xmax-xmin)/NN)

y=reflection(x,n1,d1,n2,d2,n3,d3,n0,nt,rough)

i <- seq(NN-1)

segments(x[i],y[i], x[i+1], y[i+1], col="red", lwd=2)

# add legends

legend(xmin,ymax, legend=c("1.5 / 85nm", "1.5 / 85nm", "1.5 / 85nm"), col=c("yellow", "green", "red"), lty=1:1:1, cex=0.8)

# draw visible light bondary

lab_Y = 0.05

lines(x=c(400,400),y=c(ymin,ymax*0.7), col="blue")

lines(x=c(400,700),y=c(lab_Y,lab_Y), col="blue")

lines(x=c(700,700),y=c(ymin,ymax*0.7), col="blue")

# draw color labels of visible light

# vilot 400 nm

arrow_len = (ymax-ymin)*0.1

xa=400