SITE ZHANG,1,2,*CHRISTIAN HELLMANN,2 AND FRANK WYROWSKI1
1Applied Computational Optics Group, Institute of Applied Physics, Friedrich Schiller University Jena, Max-Wien-Platz 1, 07743 Jena, Germany
2Wyrowski Photonics UG, Kahlaische Straße 4, 07745 Jena, Germany
*Corresponding author: site.zhang@uni‑jena.de
Received 8 March 2017; revised 29 April 2017; accepted 3 May 2017; posted 3 May 2017 (Doc. ID 290298); published 19 May 2017
Received 8 March 2017; revised 29 April 2017; accepted 3 May 2017; posted 3 May 2017 (Doc. ID 290298); published 19 May 2017
通過使用平面波譜分析以及S矩陣方法，我們研究了一般電磁場經過由各向同性介質或者各向異性介質構成的光學層狀結構的傳播。我們還開發了一種基于快速傅里葉變換技術的算法，具有數值高效的采樣規則。通過將此算法與其它系統建模技術相結合，我們展示了一些仿真案例，如經過一個各向同性法珀標準具的光場傳輸以及具有任意方位和光軸方向的單軸晶體平板的光場傳輸。
OCIS codes: (260.0260) Physical optics; (260.2110) Electromagnetic optics; (260.1440) Birefringence; (230.4170) Multilayers.
https://doi.org/10.1364/AO.56.004566
1.引言
由平行平面構成的光學層在光學中廣泛應用。層狀結構可以用作許多不同情況的模型，像平板和標準具�；谶@個事實，光與層狀結構相互作用的主題一直引起大家的注意并且對此已經進行了大量的研究。
在這類研究中，大多數觀點都側重于平面波，然而僅僅少數的研究使用了平面波譜方法（SPW）來考慮一般的電磁場。例如，參考文獻[1-6]中研究了各向同性-各向同性的界面上，高斯光束的反射率和透射率；在參考文獻[7-11]中研究了各向同性層或者平板的情況；參考文獻[12-22]討論了各向同性-各向異性界面的情況，在參考文獻[23-26]中則討論了各向異性層或者平板的情況。
上面所提到的許多研究都用于特定的研究主題，像[1,3,5]中研究了高斯光束全內反射的橫向偏移，并且他們常常關注于具體的配置。因此，將這些方法推廣到更一般的情況的可能性受到了限制。
在這篇文章中，我們從一個更一般的觀點來考慮此問題。光學層幾乎不會單獨使用；相反，他們常常是一個光學系統的一部分并且和其他的元件一起使用，如圖1中所示�；诖耸聦�，我們遵循場追跡的概念[27]，并使用不同的場追跡算子組合[28-32]，如圖1中所示，以對一個包含了層介質元件的系統進行物理光學模擬�？紤]到模擬是對整個系統而不是單個元件，仿真層結構必須與系統的前后部分相連接。這要求我們傳播步驟（圖1中的P）進行適當的考慮，將前一個元件的輸出連接到當前元件的輸入，并將當前元件的輸出傳遞到下一個元件。一般情況下，這樣的傳輸步驟會出現在平行或者非平行平面之間。在參考文獻[28,29]中已經提到了平行平面間幾種有效的傳輸方法，在參考文獻[33]中則可以找到對非平行平面間傳輸的一個詳細的討論。在這篇文章中，我們不會研究傳輸步驟，但會關注層狀結構的元件算子C。
此外，從數值計算的觀點出發，為了執行一個連續且有效的系統模擬，要求元件算子C
正確地處理采樣場數據并和其他的算子以一種統一的格式傳遞場數據；
優化數值計算的效率。
考慮到上述兩個標準，我們開發了一種具有自動數值采樣規則的SPW方法。與之前一些利用積分方法對空間和角譜相關的傅里葉變換進評估的研究相比（如參考文獻[23]中的二維中點規則和參考文獻[12-14，20，25]中的Stamnes–Spjelkavik–Pedersen方法[34]），我們使用了快速傅里葉變換（FFT）技術，此技術在大部分數值軟件包中容易訪問并且效率高。再加上在角譜域中經過深入考慮的數值采樣規則，我們的方法具有一般適用性，對層元件和入射場沒有任何限制。因此，此算法可以直接包含在一個物理光學系統模擬之中。
圖1.結合使用不同的場追蹤算子來模擬光學系統： C是元件算子，P是相鄰元件之間的傳輸算子。
2.理論
如圖2所示，層狀結構分別由兩個位于  和  的平行平面構成。  和  的區域充滿了復折射率為  和  的均勻各向同性介質。參考文獻[27]中表明使用橫向分量Ex和Ey已足夠表征均勻各向同性介質中電磁場了。因此，我們可以使用以下表達式來描述此問題：
其中，分別在平面 和 處定義輸入和輸出橫向電場矢量，（兩者位于界面的數學位置，但總是認為在均勻介質的一側），由下式給出
其中  。方程（1）中的元件算子是一個2x2的矩陣形式，
圖2.層狀結構分別由兩個位于 和 的平行平面構成。 和 的區域由均勻各向同性介質填充，其折射率分別是 和 。輸出場和輸出場在層表面進行定義，但總是在相應的各向同性介質的一側。
在這章節，我們的目標是找到C的精確的形式，以連接層介質元件的輸入和輸出場。為了研究與層結構的相互作用，我們對輸入橫向場分量進行了一個傅里葉變換，并獲得了
其中  , F表示二維傅里葉變換，
 。逆傅里葉變換定義如下
方程（6）中的積分可以解釋為將  分解為具有不同橫向波矢分量κ的平面波。因此，在我們的情況下，每個輸入平面波都可以單獨處理——我們首先計算每個輸入平面波的輸出，然后進行求和從而獲得輸出場。
此外，根據邊界條件對電磁場施加的連續性要求，可以顯示出一個給定的輸入平面波在與層結構相互作用的過程中其橫向波矢分量κ必定保持不變。同樣可以顯示出，通過疊加原理的有效性，不同的κ之間沒有耦合。因此，對于輸出角譜，我們可以寫下
其中
公式（8）中分別是透射和反射系數矩陣。為了計算T(κ) 或者R(κ)，我們選擇使用數值穩定S矩陣方法。為了計算S矩陣，首先必須確定每個各向異性層的平面波�；谖墨I[35]中Berreman的4x4矩陣公式，Landry和Maldonado開發并展示了一種數值友好的形式，見參考文獻[23]。我們采用了他們的方法，對于每一層，求解了參考文獻[23]中由方程（28）所描述的特征系統的特征值和特征向量。
不同于[23，25]中直接使用本征解來構建一個轉換矩陣，另外，我們還需要根據他們的傳輸方向整理出平面波，這是為計算S矩陣所做的一個必需的準備。為此，我們遵循[36]中4.3部分由Li所提出的標準。
然后可以應用遞歸S矩陣公式。我們在這篇文中不再重復給出已發展成熟的S矩陣方法，讀者可以參考文獻[37]中的方程（5）-（8）以獲得更多的信息。在我們的情況中，由于沒有反向傳輸輸入場，我們僅對正向透射或者反向反射感興趣，因此這篇文章中的矩陣系數T(κ) 和R(κ)對應于[37]中方程（5c）或者（5d）中的子矩陣T_uu 或者R_du。
一旦獲得了矩陣系數，通過方程(7)即可獲得輸出角譜。對輸出角譜進行一個逆傅里葉變換，我們獲得了輸出橫向場矢量
通過聯合方程（4），（7），（9），我們可以寫出從輸入場到輸出場整個計算流程，如果
我們以透射的情況作為例子，則
因此，方程（1）中元件算子C的精確形式如下
通過使用系數矩陣R(κ)代替 T(κ)，可以獲得反射情況的表達式。
3.算法
按照方程（10）的順序，我們可以應用一種數值算法以計算場經過分層介質元件的傳輸。讓我們從橫向輸入場矢量  和  ，以圖3（a）中的均勻網格進行采樣。這種網格定義為  ,其中  和  作為索引數，δx和δy是x方向和y方向的采樣距離。初始采樣參數應該受到合適的控制以使他們符合先前算子的奈奎斯特-香農采樣定理。方程（10）中的F為聯系空間域和頻率域的算子，可以使用不同的數值方法實現，像廣泛使用的快速傅里葉變換（FFT）技術，以及包含了能夠進一步提高數值效率的半解析傅里葉變換[38]和啁啾z變換[39-41]的更高級的方法。在此篇文章中，我們使用了FFT技術，并以此獲得輸入角譜。但我們的算法不受限于該技術。
圖3.在角譜域定義均勻采樣網格。（a），初始網格定義，5x5采樣點，  和  為采樣距離；（b），使用5x5采樣點，沿垂直方向生成一個測試網格， 和  =0.5 為采樣距離；（c）中使用9x5個采樣點定義了沿水平方向的一個測試網格，  =0.5 和 為采樣距離。（（b）和（c）中的實心點是在（a）的初始網格中出現的共同的采樣點位置，而空心圓環與初始網格中不一樣的采樣點。
其中角譜域中的均勻網格定義為  ，其中 和 為Kx和Ky的采樣距離。利用FFT計算的結果，在兩個域中的采樣點是一樣的，因此我們有  。
接下來，將輸入角譜乘以系數矩陣，然后我們可以獲得輸出角譜場，例如，對于透射
角譜  的采樣通過傅里葉變換關系自動確保。然而，  的采樣沒有必要進行確保，因為他們包含了系數和輸入角譜之間的點乘。這個問題將在第四部分的案例中進行清楚地說明。
一般來說，采樣距離 ， 和帶寬  必須進行適當調整以確保正確的采樣。在我們的情況中，方程（13）中的運算沒有改變譜范圍。因此，我們僅僅只需要找到合適的采樣距離即可。此問題并沒有解析解，因為系數矩陣T或者R是通過遞歸S矩陣方法數值上獲得的。確保 采樣的唯一方法是實行數值測試。
讓我們在角譜域定義一個測試采樣網格，如如3（b）或者3（c）所示。測試網格定義為  ,其中  為指數；為了不改變角譜范圍，要求  。與初始網格相比，它也需要進行細化，這意味著  。
一方面，我們以一種嚴謹的方式來計算測試網格上的輸出橫向角譜分量。由于輸入橫向角譜分量的合適采樣，根據需要可以進行插值，例如，在測試網格上，我們獲得
其中“INTRPL”代表的是數值插值運算。為了在一個更精細的網格上嚴格地獲得輸入場，如方程（14）中，我們總是采用基于FFT的Sinc插值方法，把它們代入到方程（13）中，可以嚴格地獲得輸出橫向角譜量
另一方面，通過對  的插值，我們也可以獲得測試網格上的輸出橫向角譜量，并將插值結果表示為  。我們將此過程描述為
值得強調的是，在方程（16）中無需再使用嚴格的Sinc插值。因為這些結果僅僅用于評估一下情況。在這篇文章中，我們使用三次插值。上面的插值沒必要給出正確的輸出值，因此，在方程（16）的左邊，我們對那些量使用上標“intrpl”。下一步，我們定義
作為插值結果和嚴格結果之間的相對偏差。只要  在初始網格上的采樣滿足奈奎斯特-香農采樣定理，基于他們的插值結果就不會對嚴格的插值方法表現出一個大的偏差。在這種情況下，根據所選的插值方法，兩種結果的偏差應該在數值誤差級次，由σ_0表示。使用上述標準，通過一步步的減小采樣間距來測試場采樣，直至σ<σ_0。
以上的測試過程對應著一個循環，在每個測試循環中都會執行方程（14）到（17）的計算。為了在每個測試循環中充分地使用測試數據，如方程（15）中的結果，我們總是使用圖3（b）和3（c）中的測試網格。這樣規則的測試網格定義為
對每個循環內，需要對先前的采樣距離減半。使用這個測試網格具有如下好處：當計算  時，僅僅需要計算那些在空心圓圈位置處的值，而實心點位置處的值在之前的測試循環中已經計算過了，并且讓這些值簡單地進行傳輸即可。在這種方式下，這些基于測試目的所進行嚴格計算的值不會被丟掉，而會作為下一次循環中的起始點以重新使用。因此，在一個完整的計算中所執行的是嚴格的計算，從而驗證合適的采樣，并且這些計算的值能夠有效地用于構建最終輸出場。我們將上面的方法總結在算法1中。
算法1：通過合適的采樣控制來進行輸出角譜計算的數值流程
1)在初始網格嚴格地計算
2)初始化相對誤差值σ=+∞
3)判斷σ>σ_0， ⊳如果為真，則采樣不合適
4)沿Kx（或者Ky）二等分采樣距離
5)應用新的采樣距離定義測試網格 
6)在測試網格中對輸入角譜插值，根據方程（14）獲得 
7)根據方程（15），在測試網格上嚴格地計算輸出角譜 
8)根據方程（16），通過插值獲得 
9)根據方程（17），計算誤差σ
10)如果σ>σ_(0 )  ，則⊳將目前的輸出場設置為下一次循環的起始點
11)返回
為了有效地處理非對稱情況，例如，光束在x方向和y方向有不同的發散角，算法1中的測試需要按順序沿兩個方向進行。開始方向的選擇是任意的，在我們的情況中，我們是沿y方向開始測試的。
4.示例
在VirtualLab Fusion[42]軟件中，我們將第三節中提出的算法實現在“可編程元件”的編程界面中。這個元件可以與VirtualLab Fusion中其它的物理光學仿真技術進行聯合仿真。接下來，我們展示了四個案例：前兩個主要關注元件本身并以一種嚴格的數值方式檢查此算法；后兩個案例中，元件將用于光學系統中，例如，此算法與其它仿真技術一起使用。
在進入實際的案例之前，設置方程（17）中的迭代終止標準σ_0很重要。對于在此文中所使用的三次插值，我們預先檢查了它在一般情況下的表現并在我們的數值環境中找到了一個0.01的基準值。
A.各向同性標準具
第一個案例模擬了一個線性偏振高斯光束經過一個標準具的傳輸。通過這個案例，我們將清楚的指出第三節中所說的采樣問題并描述了算法1的工作原理。標準具由熔融石英制成，兩側有多層薄膜，如圖4所示。關于標準具的光學參數和結構的更多信息，請見表1。輸入場為波長633nm，x方向線偏光的高斯光束。其在元件的輸入平面定義為E_⊥^in,且束腰半徑為（2um，2um）。在經過傅里葉變換后，我們獲得其角譜，同樣具有高斯輪廓，如圖5所示。按照方程（10）中的操作算子序列，輸入角譜將乘以透射或者反射率系數。我們仍以透射作為例子，并且對于線性偏振輸入場，我們使用t_xx和t_yx乘以E ̃_x^in，以獲得輸出角譜分量。
圖4.由熔融石英制成，兩側有多層膜的標準具。其結構和光學參數如表1中所示。
表1 標準具的結構和光學參數
圖5 .（a）輸入高斯場分量的振幅；（b）對應的角譜分量。由于輸入光場為沿x方向的線性偏振光，因此僅顯示Ex分量。
如第三節中所指出，乘積  的采樣不能自動得到保證，此案例中將顯示該現象。標準具由于其頻率選擇功能（頻譜或角頻率）而得到廣泛的使用。在我們的案例中，角頻率選擇可以解釋為系數 txx和 tyx以一種方式調制輸入角譜，以使特定的角頻率加強而其它的減弱。這種調制可以出現在一種非常精細的頻率水平上。因此，需要使用更精細的采樣以在輸出角譜中解析這樣一個精細的調制。為了獲得需要的采樣間距，我們遵循算法1，圖6中顯示了部分結果。
通過圖6和表2，我們根據算法1中的步驟描述了工作流程，如下：
第一步：從圖6（a）中所示的輸入角譜開始，計算各個系數并乘上元件矩陣以生成初始化的輸出角譜；
第二步：初始化相對偏差σ=+∞；
第三步：開始測試循環；
第四&五步：將采樣距離沿κ_x或者κ_y方向減半，以定義測試網格 ，對應的采樣點如表2中所示；
第六步：對輸入角譜在測試網格上插值；
第七步：在測試網格上嚴格的計算輸出角譜 ，在此例中，對應圖6（e）；
第八步：執行插值以獲得 ，在此例中對應圖6（b）-6（d）；
第九步：比較嚴格仿真和插值結果，并計算相對誤差；
第十步：對于較差的插值結果，如6（b）和6（c），表2中的0-3行，其結果是σ>σ_0，嚴格的結果 將會傳遞給下一個循環并用于輸入；否則，返回當前的結果。
圖6.算法1中不同步驟時沿κ_x方向的一維提取結果：（a）輸入角譜振幅，（b）-（d），在測試循環中的插值角譜振幅以及（e）在最后的循環中嚴格地計算輸出角譜。所有的子圖中的值都縮放到相同的范圍內。
從表2中我們也可以看到測試首先是沿y方向，之后再沿x方向執行的，如第三節的最后所提到的。從第0輪到第4輪測試，采樣距離δκx并沒有改變，因此采樣點Mx^，保持不變；在第四輪測試后σ<σ0，沿y方向的測試終止，意味著場數據已經可以從先前一輪的結果中恢復。因此，沿y方向的采樣點數是705，在第四輪中額外的704個數據僅僅是用于測試目的，對最終的輸出并沒有貢獻。然后以一種類似的方式沿x方向開始測試并在第11輪終止，同樣，場數據可以從先前一輪的結果中恢復。因此，最終輸出的采樣點數固定在2817x705。在表2的測試輪中，包含的總的采樣點數是5633x705+45x704，數據45x704來自于沿y方向的最后測試輪。再次強調一下，為了在角譜域進行合適的采樣，有必要采用如此大的采樣點數。并且，除了沿x方向或者y方向最后一輪測試，用于測試目的而所有嚴格計算地值都用于構建最終輸出場。
表2 Etalon模擬中每個測試輪次中的采樣參數和誤差
文章中的模擬，使用的是一臺Intel Core i7-4910MQ處理器,2.9GHz，32Gb物理內存的電腦。表2中所顯示的示例所需的總的計算時間是118s。注意計算時間最多的部分是花在S矩陣計算上，而最后步驟的傅里葉逆變換僅花了約0.5s。
從表2中也可以看出沿x方向和沿y方向的測試輪次數也不一樣，即沿兩個方向所需的采樣間距不同。因此，我們在算法中更傾向沿每個方向分別執行采樣測試。
通過這個算法，我們獲得了標準具的透過場，如圖7所示。除了圖7中透射場的尺寸遠大于圖5中的輸入場，非零輸出場分量  也值得注意，盡管其振幅要遠小于  。Ex和Ey之間出現的偏振串擾是由于系數矩陣T或者R的非對角形式造成的。從 的振幅分布中，還可以看出精細的同心圓，這是由于標準具多層結構中發生的多重反射造成的。
圖7.以任意單位給出的標準具輸出場分量的振幅，（a）：|E_x^out |，（b）：|E_y^out |
B.單軸晶體板
我們的方法對各向同性介質層和各向異性介質層都同樣適用。在這個案例中，我們演示了聚光干涉儀的原理，其可以用于精確的測量單軸晶體的光軸傾斜角度[43]。當使用會聚單色光照射正交偏振器間的晶體時，即可產生聚光干涉。圖8中展示了聚光干涉儀的簡要原理。在我們的模擬中，光源是一個Ex線性偏振，波長為633nm，NA=0.25的球面波，在晶體板前20mm。光源將輸入場  傳遞到晶體板的前表面，后表面的輸出場將被分析。在這個子節中，我們沒有討論晶體板元件外部區域的傳播步驟。
圖8.使用聚光干涉儀測量光軸的角度。偏振片P1沿x方向，偏振片P2沿y方向，光軸表示為“o.a”。
我們使用聚光干涉儀，測試了一個6mm厚的向列液晶板，其光軸與z軸成θ角。其中n_e=1.7，n_0=1.5[43]，設θ=0°，2°以及5°，我們觀察了晶體板后y方向的輸出場，如圖9所示。當θ=0°，晶體的光軸垂直于晶體板表面，同心環的原點位于中心，如圖9（a）所示；當θ≠0°時，如此例中的圖9（b）和
圖9（c），同心環的原點相應的發生了橫向偏移。我們在模擬中獲得的值與參考文獻[43]中分析的結果一致。
圖9.針對不同的θ值，在晶體板后表面的輸出場分布（振幅平方）。使用沿y方向的偏振片P2分析輸出場，因此僅顯示了E_y分量。紅色的十字標出了同心圓環的原點，為（0，0），（210um，0）以及（525um，0），分別對應子圖（a），（b）和（c）。
為了獲得圖9中的透射場，算法1再次被用于確定采樣距離。為了完整性，表3中給出了θ=5°時此例中采樣距離的數值測試的細節。
在表3所顯示的測試輪次中，包含了3969x1985+63x1984個采樣點，對于每個采樣點，計算了各向異性S矩陣，耗時416s。
表3.光軸為5°的單軸晶體板在每個測試輪次中的采樣距離和偏差
C.傾斜的各向異性晶體板系統
如第一節所述，我們開發的算法是系統建模的一部分。在這個例子中，我們模擬了一個包含方解石晶體板的光學系統以研究偏振轉換以及渦旋的生成[44]，如圖10所示。因為是針對整個系統進行建模，因此在處理過程使用了幾種不同的場追跡算子[27]：偏振片使用瓊斯矩陣算子，透鏡被當作是理想透鏡，到傾斜方解石板表面的光場傳輸使用參考文獻[33]中的技術，方解石板中的光場傳輸則使用了此文中的方法。
根據參考文獻[44]，給出了方解石晶體折射率，其中
其中λ單位為微米。
圖10.包含傾斜晶體板的系統。晶體板材料為方解石，厚度為6mm，光軸在圖中以“o.a.”表示，垂直于晶體板表面。偏振片P1用于生成沿x方向的偏振光，而偏振片P2用于分析。透鏡L1和L2焦距都為30mm。
此例中我們從α=0°開始，在系統中幾個位置的場分布如圖11所示。在偏振片P1前的輸入場為波長633nm，x方向和y方向束腰半徑1.5mm的高斯光束。經過偏振片P1之后，獲得了一個沿x方向的線性偏振場。透鏡L1后的聚焦場是方解石晶體板的輸入場。根據算法1，我們獲得了晶體板后表面的輸出場，如圖11（a）和11（b）。最后，在偏振片P2之后的平面上，我們獲得了系統的輸出場，如圖11（c）和11（d）所示。
圖11.圖10中方解石晶體板及P2偏振片后的場分布（振幅平方），此例中α=0°。
讓我們來仔細地觀察一下圖11中所展現的結果。數值測試細節在表4中給出。與4A部分的情況相似，由于T或者R的非對角線形式，沿x方向的線性偏振輸入場為輸出場增加了一個非零y偏振方向。但此例中單軸晶體板的偏振串擾相對更強，E(x )≈E(y ),而在標準具的例子中，E(x )和E(y )的比例大概是1000：1，如圖7（a）和7（b）。這是由于方解石晶體的雙折射性質，導致E(x )和E(y )分量之間具有一個強的耦合。
表4 1.5mm高斯光束入射到單軸晶體板，在每個測試輪次中的采樣距離和偏差
P2后的場，如圖11（e）和11（f）所示，與參考文獻[44]中的測量結果吻合的很好。除了由雙折射引起的四極結構，我們同樣也看到了疊加在其上的同心圓環結構。與圖11的分析類似，這些精細結構是由于板內的多重反射造成的�？梢栽趨⒖嘉墨I[44]中的圖4中看到同樣的影響，盡管實驗測量的對比度較低。
下一步，如參考文獻[44]中所展示，我們將晶體傾斜一定的角度α并將輸入高斯束腰設置為0.5mm，以生成單極光渦旋。模擬中，我們使用[33]中的方法處理非平行板之間的傳輸。通過這種方式，可以獲得傾斜晶體板表面的輸入場，且傾斜板表面的輸出場可以繼續傳輸至下一個元件。模擬結果如圖12和13所示。當α=1.2°時，可以在表5中獲得采樣距離的數值測試細節�？梢钥吹�，沿x方向和y方向的偏差σ分別在輪次2和5中滿足停止判據。我們沿兩個方向繼續往前迭代一次（標記有*）以獲得一個過采樣，以為了更好的與之后的非平行平面間的傳播相結合[33]，這之間的傳播通常需要過采樣因子為2.
圖12.當α=0.8°時，方解石晶體板和偏振片P2之后的輸出場（振幅平方）
圖12和13再次與參考文獻[44]中的結果一致。與之前分析α=0°時的情況類似，光渦旋結構是由方解石晶體的雙折射產生的。如圖13（c）和13（d）所示，同樣也可以看到相位錯位引起的光渦旋。除了渦旋，圖13（b）中還可以看到環形的精細結構，這是由于晶體板內的多重反射造成的。
將表5中的數據測試與之前的進行比較，我們可以看到一個更快的收斂，特別是沿x方向的測試。包括表5中的所有測試循環，經過傾斜的方解石晶體板的傳輸僅使用了大約4s的計算時間就完成了，而整個系統的模擬，包含透鏡，大約是30s的時間。由于方解石板繞y軸傾斜1.2°，若將晶體板表面作為參考，則對應著一個輕微的傾斜入射場。根據傅里葉變換偏移理論，沿y方向的傾斜將導致角譜沿此方向有一個偏移。因此，透射場的角譜同樣存在沿x方向和y方向的不同調制。
圖13.當α=1.2°時，方解石晶體板和偏振片P2后的輸出場，其中（a）和（b）為振幅平方的分布，（c）和（d）為相位分布
表5 0.5mm高斯光束入射到傾斜晶體板時，每輪測試中的采樣參數和偏差
D.在空氣-單軸表面的反射
對于反射情況，同樣能以類似于算法中的透射情況的方法進行處理，我們將在下文展示。我們選用文獻[45]的工作，即空氣-單軸晶體表面發生的自旋霍爾效應，來作為一個例子。系統原理圖如圖14所示，且自旋相關納米量級偏移可以通過檢查反射場的中心來進行測量。
圖14.在空氣-單軸晶體表面的反射。晶體元件使用  坐標系統來表示，其光軸（o.a.）沿  方向，輸入和反射場分在  中給出
在我們的模擬中，使用了波長為632.8nm的線性偏振高斯場作為輸入場，在此波長下，對于o光和e光，LiNbO3晶體樣品的折射率分別為n0=2.232和ne=2.156。在文獻[45]中并沒有給出輸出光場的光束尺寸，我們將其設置為20um。由于反射場的橫向偏移是由空間頻率域中的局部線性相位決定的，因此不改變光束尺寸的大小。
為了研究此效應，反射矢量場投影在圓偏單位矢量坐標系，而后我們測量了右旋圓偏振分量的中心。此外，對Ex和Ey兩個線性偏振入射場，我們旋轉晶體元件并監控位移隨著入射角度的變化。模擬結果如圖15中所示。
我們使用參數掃描獲得圖15中的結果。對于每個入射角度，我們使用我們的算法計算了反射場，且在一個大的角度范圍內顯示出了較好的適應性。模擬結果與[45]中的結果匹配的很好。
圖15.入射角度與右圓偏振量反射場的橫向偏移之間的關系，（a）Ex和（b）Ey分別為線性偏振輸入。
總結
我們從系統仿真的角度研究了電磁場經過各向同性或各向異性層介質元件的傳輸。利用SPW分析并將以方程（10）的形式概括的結果，以為后續的數值實現做準備。我們清楚地討論了在很多實際問題中會遇到的采樣問題，此外，我們提出了一種數值測試算法，以更有效地在角譜域確定采樣參數。通過標準具的案例，我們詳細地描述了我們的算法的工作流程，顯示了我們方法的一般有效性。此方法已經成功地應用到了對激光晶體封裝技術中的應力誘導雙折射的分析[46]。
致謝
我們感謝Olga Baladron-Zorita女士對此文章的校正以及其日常的幫助。
參考文獻
B. R. Horowitz and T. Tamir, “Lateral displacement of a light beam at a dielectric interface*,” J. Opt. Soc. Am. 61, 586–594(1971).
M. McGuirk and C. K. Carniglia, “An angular spectrum representation approach to the Goos–Hänchen shift,” J. Opt. Soc. Am. 67, 103–107(1977).
S. Kozaki and H. Sakurai, “Characteristics of a Gaussian beam at a dielectric interface,” J. Opt. Soc. Am. 68, 508–514 (1978).
C. C. Chan and T. Tamir, “Angular shift of a Gaussian beam reflected near the Brewster angle,” Opt. Lett. 10, 378–380 (1985).
F. I. Baida, D. V. Labeke, and J.-M. Vigoureux, “Numerical study of the displacement of a three-dimensional Gaussian beam transmitted at total internal reflection near-field applications,” J. Opt. Soc. Am. A17, 858–866 (2000).
S. Zhang, D. Asoubar, F. Wyrowski, and M. Kuhn, “Efficient and rigorous propagation of harmonic fields through plane interfaces,”Proc. SPIE 8429, 84290J (2012).
M. Tanaka, K. Tanaka, and O. Fukumitsu, “Transmission and reflection of a Gaussian beam at oblique incidence on a dielectric slab,” J. Opt. Soc. Am. 67, 819–825 (1977).
S. Kozaki and H. Harada, “Beam displacement of a reflected beam at an interface between an inhomogeneous medium and free space,”J. Opt. Soc. Am. 68, 1592–1596 (1978).
C. W. Hsue and T. Tamir, “Lateral displacement and distortion of beams incident upon a transmitting-layer configuration,” J. Opt.Soc. Am. A 2, 978–988 (1985).
R. P. Riesz and R. Simon, “Reflection of a Gaussian beam from a dielectric slab,” J. Opt. Soc. Am. A 2, 1809–1817 (1985).
T. Tamir, “Nonspecular phenomena in beam fields reflected by multilayered media,” J. Opt. Soc. Am. A 3, 558–565 (1986).
J. J. Stamnes and D. Jiang, “Focusing of electromagnetic waves into a uniaxial crystal,” Opt. Commun. 150, 251–262 (1998).
J. J. Stamnes and V. Dhayalan, “Transmission of a two-dimensional Gaussian beam into a uniaxial crystal,” J. Opt. Soc. Am. A 18,1662–1669 (2001).
J. J. Stamnes and V. Dhayalan, “Double refraction of a Gaussian beam into a uniaxial crystal,” J. Opt. Soc. Am. A 29, 486–497(2012).
L. I. Perez, “Reflection and non-specular effects of 2d Gaussian beams in interfaces between isotropic and uniaxial anisotropic media,” J. Mod. Opt. 47, 1645–1658 (2000).
L. I. Perez, “Nonspecular transverse effects of polarized and unpolarized symmetric beams in isotropic-uniaxial interfaces,” J. Opt. Soc.Am. A 20, 741–752 (2003).
S. Stallinga, “Axial birefringence in high-numerical-aperture optical systems and the light distribution close to focus,” J. Opt. Soc. Am.A 18, 2846–2859 (2001).
S. Stallinga, “Light distribution close to focus in biaxially birefringent media,” J. Opt. Soc. Am. A 21, 1785–1798 (2004).
M. Jain, J. Lotsberg, J. Stamnes, and Ø. Frette, “Effects of aperture size on focusing of electromagnetic waves into a biaxial crystal,” Opt.Commun. 266, 438–447 (2006).
M. Jain, J. K. Lotsberg, J. J. Stamnes, Ø. Frette, D. Velauthapillai, D.Jiang, and X. Zhao, “Numerical and experimental results for focusing of three-dimensional electromagnetic waves into uniaxial crystals,”J. Opt. Soc. Am. A 26, 691–698 (2009).
S. Zhang and F. Wyrowski, “Simulations of general electromagnetic fields propagation through optically anisotropic media,” Proc. SPIE 9630, 96300A (2015).
D. Asoubar, S. Zhang, and F. Wyrowski, “Simulation of birefringenceeffects on the dominant transversal laser resonator mode caused by anisotropic crystals,” Opt. Express 23, 13848–13865(2015).
G. D. Landry and T. A. Maldonado, “Gaussian beam transmission and reflection from a general anisotropic multilayer structure,” Appl. Opt.35, 5870–5879 (1996).
V. Dhayalan and J. J. Stamnes, “Focusing of electromagnetic waves into a dielectric slab, I: exact and asymptotic results,” Pure Appl. Opt.7, 33–52 (1998).
V. Dhayalan and J. Stamnes, “Focusing of electromagnetic waves into a dielectric slab, II: numerical results,” J. Eur. Opt. Soc. 6,11036 (2011).
A. Turpin, Y. V. Loiko, T. K. Kalkandjiev, and J. Mompart, “Light propagation in biaxial crystals,” J. Opt. 17, 065603 (2015).
F. Wyrowski and M. Kuhn, “Introduction to field tracing,” J. Mod. Opt.58, 449–466 (2011).
D. Asoubar, S. Zhang, F. Wyrowski, and M. Kuhn, “Parabasal field decomposition and its application to non-paraxial propagation,”Opt. Express 20, 23502–23517 (2012).
D. Asoubar, S. Zhang, F. Wyrowski, and M. Kuhn, “Efficient semianalytical propagation techniques for electromagnetic fields,”J. Opt. Soc. Am. A 31, 591–602 (2014).
H. Zhong, S. Zhang, and F. Wyrowski, “Parabasal thin-element approximation approach for the analysis of microstructured interfaces and freeform surfaces,” J. Opt. Soc. Am. A 32, 124–129 (2015).
F. Wyrowski, H. Zhong, S. Zhang, and C. Hellmann, “Approximate solution of Maxwell’s equations by geometrical optics,” Proc. SPIE9630, 963009 (2015).
F. Wyrowski and C. Hellmann, “Geometrical optics reloaded,”Optik Photonik 10, 43–47 (2015).
S. Zhang, D. Asoubar, C. Hellmann, and F. Wyrowski, “Propagation of electromagnetic fields between non-parallel planes: a fully vectorial formulation and an efficient implementation,” Appl. Opt. 55, 529–538 (2016).
34. J. J. Stamnes, B. Spjelkavik, and H. M. Pedersen, “Evaluation of diffraction integrals using local phase and amplitude approximations,”Opt. Acta 30, 207–222 (1983).
D. W. Berreman, “Optics in stratified and anisotropic media:4 × 4-matrix formulation,” J. Opt. Soc. Am. 62, 502–510 (1972).
L. Li, “Reformulation of the Fourier modal method for surface-relief gratings made with anisotropic materials,” J. Mod. Opt. 45,1313–1334 (1998).
L. Li, “Note on the s-matrix propagation algorithm,” J. Opt. Soc. Am. A20, 655–660 (2003).
Z. Wang, “Analytical handling of optical wavefront,” Master’s thesis(Friedrich Schiller University, 2015).
L. Rabiner, R. Schafer, and C. Rader, “The chirp z-transform algorithm,” IEEE Trans. Audio Electroacoust. 17, 86–92 (1969).
J. L. Bakx, “Efficient computation of optical disk readout by use of the chirp z transform,” Appl. Opt. 41, 4897–4903 (2002).
M. Leutenegger, R. Rao, R. A. Leitgeb, and T. Lasser, “Fast focus field calculations,” Opt. Express 14, 11277–11291 (2006).
Physical optics design software, “Wyrowski virtuallab fusion,” developed by Wyrowski Photonics UG, distributed by LightTrans GmbH. http://www.lighttrans.com.
B. L. V. Horn and H. H. Winter, “Analysis of the conoscopic measurement for uniaxial liquid-crystal tilt angles,” Appl. Opt. 40, 2089–2094(2001).
Y. Izdebskaya, E. Brasselet, V. Shvedov, A. Desyatnikov, W.Krolikowski, and Y. Kivshar, “Dynamics of linear polarizationconversion in uniaxial crystals,” Opt. Express 17, 18196–18208(2009).
Y. Qin, Y. Li, X. Feng, Z. Liu, H. He, Y. Xiao, and Q. Gong, “Spin Hall effect of reflected light at the air-uniaxial crystal interface,” Opt.Express 18, 16832–16839 (2010).
P. Ribes-Pleguezuelo, S. Zhang, E. Beckert, R. Eberhardt, F.Wyrowski, and A. Tünnermann, “Method to simulate and analyse induced stresses for laser crystal packaging technologies,” Opt.Express 25, 5927–5940 (2017).
|
