有限元矩陣的激光諧振腔模式分析

陳 凱，朱東旭，張平才

(中國衛(wèi)星海上測控部，江陰214431)

式中，λ為激光波長;k為波數，表示在2π長度上出現的全波數目;U1(r1，φ1)為腔鏡1上任意一點P的場分布;ρ為兩腔鏡間任意兩點之間的距離;r1∈(0，a1)，φ1∈(0，2π)，r2∈(0，a2)，φ2∈(0，2π);S1為腔鏡1的積分面元。根據積分需要，作如下參量

引 言

激光光束的質量決定了激光的聚焦特性和傳輸特性，它與激光諧振腔內的光場模式密切相關。經典諧振腔理論已經給出了部分典型諧振腔的模式分析方法和解析解。但是，在激光器的不斷發(fā)展過程中，出現了各種新型的腔型和光學元件，往往用解析方法難以計算，必須轉而采用數值的方法求解，如Fox-Li法［1］、Prony 法［2］、快速傅里葉變換法(fast Fourier transform，FFT)［3］、有限差分法(finite differential method，FDM)［4］和有限元矩陣法(finite element matrix，FEM)［5］。由于 FEM 在設置單元的大小和形狀時很方便，因而在處理復雜幾何曲面時更加有效，具有較強的適應性、較高的精度等特點。FEM被提出以來，一直受到國內外學者的廣泛關注和研究［6-8］，并在波導光學領域中得到了非常廣泛的應用。作者的研究目的就是在FEM的基礎上，做進一步的研究和改進，使激光諧振腔模式的分析更加精確和方便。

1 理論分析

FEM的基本思想是將一個連續(xù)的求解區(qū)域離散化為有限個單元，用求解域上待求的未知場函數在各個單元結點的數值及其差值函數來分片地表示每個單元的場函數［9］。隨著單元數目的增加或者單元自由度的增加及插值函數精度的提高，解的近似程度將不斷改進。如果單元是滿足收斂要求的，近似解最后將收斂于精確解。

Fig.1 General spherical cavity

式中，λ為激光波長;k為波數，表示在2π長度上出現的全波數目;U1(r1，φ1)為腔鏡1上任意一點P的場分布;ρ為兩腔鏡間任意兩點之間的距離;r1∈(0，a1)，φ1∈(0，2π)，r2∈(0，a2)，φ2∈(0，2π);S1為腔鏡1的積分面元。根據積分需要，作如下參量

由于反射鏡繞腔軸旋轉對稱，故腔內的光場分布也應具有繞腔軸的旋轉對稱性［10］。不妨設:

式中，E2為關于復振幅的函數，l為角向橫模序數。

將(2)式～(3)式代入(1)式，可得:

式中，E1也為關于復振幅的函數。對于計算過程中的積分項，式中作出如下簡化:

現考慮腔鏡1上的場分布對腔鏡2上的場分布的影響。將腔鏡1按照一定的順序劃分為單元1～M。當劃分數M足夠大時，可認為腔鏡1上每個單元的復振幅起伏不大，即復振幅均勻分布。在這種情況下，可用復振幅U1［m］來代替腔鏡1上第m個單元的復振幅分布。同理，將腔鏡2劃分為單元1～N。對于腔鏡 2上的第n個單元，其復振幅U2［n］可認為是腔鏡1上單元1～M作用并疊加的結果:

由于諧振腔的腔長L通常遠大于腔鏡的孔徑a，并且貝塞爾函數Jl具有強收斂性［12］。故在a1m～a1m+Δa1m的區(qū)間里進行積分時，可將被積函數代入具體數值a1m后看作常數，提到積分號外面。為了方便計算，對圓形腔鏡只進行徑向劃分，不考慮角向劃分，可得:

式中，Δa1m=a1/M，則腔鏡1第m個單元的半徑為a1m=ma1/M;Δa2n=a2/N，則腔鏡2第n個單元的半徑為a2n=na2/N。簡化(7)式，可得:

式中，A12［n，m］為腔鏡1到腔鏡2的渡越公式，與兩腔鏡上的離散劃分位置有關:

于是，腔鏡2上的場分布可寫成矩陣形式:

(10)式表征了光場由腔鏡1到腔鏡2的一次渡越，簡記為U2=A12×U1。

同理，腔鏡2上的場分布作用于腔鏡1上可有:

不難看出，A12與A21互為轉置矩陣。對于腔鏡1上的場分布，簡記為U1=A21×U2。

腔鏡上的光場在腔內一次往返后，其復振幅分布Ui'可表示為:

根據自再現原理［13-15］，經過足夠多次渡越以后，在腔內往返一次后能夠“再現”出發(fā)時的場分布。唯一可能的變化是，鏡面上各點的場振幅按同樣的比例衰減，各點的相位發(fā)生同樣大小的滯后。于是，有如下關系:

式中，γ表征諧振腔各階模式的本征值;本征向量Ui表示不同γ對應的諧振腔各階模式分布;矩陣A為傳輸矩陣，描述光場復振幅在腔內往返一次的變化，即有A=Aji×Aij。當諧振腔為對稱諧振腔時，渡越矩陣Aij=Aji=A，單程渡越就能實現模的自再現。

2 實例分析及討論

下面以圓形鏡諧振腔為研究對象，利用FEM方法對諧振腔模式進行分析。設定激光器波長為632.8nm，兩腔鏡的曲率半徑分別為1m和1000m(可近似看作平面鏡)，孔徑分別為0.5mm和1mm，劃分格數均為50，腔長為0.2m。

根據(13)式，可得到每個自再現模在整個鏡面上的場分布，包括振幅分布和相位分布。根據計算的結果，圖2表征了TEM00模到TEM03模的本征值γ、單程損耗)和振幅分布。

由圖2可以看出，TEM00模在腔內的模式競爭中具有明顯的優(yōu)勢。因此，FEM方法的計算結果可以很好地用來分析諧振腔的模式鑒別能力，該計算結果與用 Fox-Li法求得的解析解完全相符［16-20］。事實上，FEM方法可用于任何結構諧振腔的模式計算或光束傳輸模擬。

Fig.2 Amplitude distribution of different modes in round mirror cavity

3 單元格數的合理劃分

菲涅耳數N'是諧振腔的一個十分重要的參量，它能夠表征諧振腔的衍射損耗，其表達式為:

菲涅耳數N'對模式增益的影響是物理意義上的，N'越大，衍射損耗越小，故而模式增益越大;而單元格數G對模式增益的影響只是數值意義上的，G越大，劃分越精細，模式增益越大也越趨于穩(wěn)定。相比之下，單元格數G對模式增益的影響較小，是在衍射損耗的基礎上使結果進一步精確化，更逼近真實值。因此，在菲涅耳數N'相同的條件下確定合理的單元格數Gb更為恰當。

3.1 單元格數的合理劃分

為了方便分析單元格數對模式計算的影響，這里以對稱圓形鏡諧振腔作為研究對象，即兩腔鏡的曲率半徑R1=R2=R，兩腔鏡的孔徑a1=a2=a。兩腔鏡面上選擇相同的單元劃分方案，單元格數均設為G。為了能將菲涅耳數N'考慮進來，不妨將劃分格數G除以菲涅耳數N'，即G/N'，作為單元格精細程度。以此作為自變量，并不影響單元格數G變化而菲涅耳數N'不變的前提。

運用多重循環(huán)的方法對波長λ、孔徑a、腔長L進行取值，擬合得到不同單元格精細程度G/N'與相對誤差d的關系式(如圖3所示):

從而不難得出合理單元格數Gb與相對誤差d的關系:

式中，d為分析結果隨單元格數G增加的相對誤差值。

該方法的優(yōu)勢就在于可根據不同的相對誤差d要求確定合理的單元格數Gb。

Fig.3 Curve fitting between the cell sophistication G/N'and the relative error of mode gain d

3.2 實例計算與分析

對于一個腔長為0.6m、兩腔鏡的孔徑均為1mm、激光波長為632.8nm的諧振腔，要求計算的相對誤差d=0.1%。根據(16)式計算得合理單元格數Gb應為26.91。分別將單元劃分格數G設置為10，20，25，30，35和40，計算相應本征值絕對值的最大值(見表1)以及TEM01模的光場分布(見圖4)。

Table 1 The maximum absolute value of the eigenvalue

Fig.4 The optical field distribution under different number of elements G(TEM01mode)

由表1可見，當單元格數G達到合理單元格數Gb后再增加單元格數計算的相對誤差分別為0.08%，0.04%和0.03%，均小于所要求的計算相對誤差0.1%。同樣不難從圖4可以發(fā)現，單元格數G達到合理劃分格數Gb后，TEM01模的光場分布變化也較小，基本趨于穩(wěn)定。

4 結論

從基爾霍夫-菲涅耳衍射積分方程出發(fā)，采用FEM劃分的思想對腔鏡進行徑向單元劃分，并將單元節(jié)點處的光場復振幅變量分離出來，可一次求解多個模式，與初值的選擇無關，適用于復雜腔型的計算。對腔鏡單元的合理劃分進行了有益的探索，結果符合較好。

［1］ BOYD G D，KOGELNIK H.Generalized confocal resonator theory［J］.Bell System Technical Journal，1962，41(4):1347-1369.

［2］ GORDON J P，KOGELNIK H. Equivalencerelationsamong spherical mirror optical resonators［J］.Bell System Technical Journal，1964，43(6):2873-2886.

［3］ FLAMMER C.Spheroidal wave functions［M］.Stanford，USA:Stanford University Press，1957:251-273.

［4］ HEURTLEY J C，STREIFER W.Optical resonator modes-circular reflectors of spherical curvature［J］.Journal of the Optical Society of America，1965，55(11):1472-1479.

［5］ FOX A G，LI T.Resonant modes in a maser interferometer［J］.Bell System Technical Journal，1961，40(2):453-488.

［6］ SIEGMAN A E，MILLER H Y.Unstable optical resonator loss calculations using the prony method［J］.Applied Optics，1970，9(12):2729-2736.

［7］ QIN Y X，TANG X H，ZHONG R T，et al.Study on mode calculation of transmission matrix for laser resonator［J］.Chinese Journal of Lasers，2008，35(10):1463-1468(in Chinese).

［8］ YONG F Q，WANG Y Q.Numerical matrix method for calculation of resonator mode［J］.Optics ＆ Optoelectronic Technology，2005，3(5):40-42.

［9］ SIEGMAN A E，SZIKLAS E A.Mode calculations in unstable resonators with flowing saturable gain.Ⅰ.Hermite-Gaussian expansion［J］.Applied Optics，1974，13(12):2775-2791.

［10］ ZHUANG G L，DING T Q.A model for numerical analysis of cavity-modes and far-fields of circular-hole and ring-coupling laser resonators［J］.Laser Journal，1987，8(3):188-192(in Chinese).

［11］ SZIKLAS E A，SIEGMAN A E.Mode calculations in unstable resonators with flowing saturable gain.Ⅱ.Fast Fourier transform method［J］.Applied Optics，1975，14(8):1874-1889.

［12］ LIU Z Y，WU X Y，FAN D Y.Collins formula in frequency-domain and fractional Fourier transform［J］.Optics Communications，1998，155(1/3):7-11.

［13］ CHEN J B.Laser principle and application［M］.Beijing:Electronic Industry Press，2004:261-267(in Chinese).

［14］ CHENG Y Y，WANG Y Q，HU J.A novel eigenvector method for calculation of optical resonator modes and beam propagation［J］.Acta Physica Sinica，2004，53(8):2576-2580(in Chinese).

［15］ GAO Zh H.Wave pattern property of self-reproductive mode in laser resonator［J］.Acta Photonica Sinica，2000，29(8):726-729(in Chinese).

［16］ HELFERT S F，PREGLA R.A finite difference beam propagation algorithm based on generalized transmission line equations［J］.Optical and Quantum Electronics，2000，32(6/8):681-690.

［17］ STOFFER R，HOEKSTRA H J W M.Efficient interface conditions based on a 5-point finite difference operator［J］.Optical and Quantum Electronics，1998，30(5/6):375-383.

［18］ VASSALLO C.1993-1995 optical mode solvers［J］.Optical and Quantum Electronics，1997，29(2):95-114.

［19］ POLYCHRONOPOULOS S J，ATHANASOULIAS G B，UZUNOGLU N K.Advanced mode solver using an integral equation technique and entire domain plane wave basis functions［J］.Optical and Quantum Electronics，1997，29(2):127-137.

［20］ SELLERI S，VINCETTI L，CUCINOTTA A，et al.Complex FEM modal solver of optical waveguides with PML boundary conditions［J］.Optical and Quantum Electronics，2001，33(4):359-371.