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

陳 凱，朱東旭，張平才

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

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

引 言

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

1 理論分析

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

Fig.1 General spherical cavity

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

由于反射鏡繞腔軸旋轉(zhuǎn)對稱，故腔內(nèi)的光場分布也應(yīng)具有繞腔軸的旋轉(zhuǎn)對稱性［10］。不妨設(shè):

式中，E2為關(guān)于復(fù)振幅的函數(shù)，l為角向橫模序數(shù)。

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

式中，E1也為關(guān)于復(fù)振幅的函數(shù)。對于計(jì)算過程中的積分項(xiàng)，式中作出如下簡化:

現(xiàn)考慮腔鏡1上的場分布對腔鏡2上的場分布的影響。將腔鏡1按照一定的順序劃分為單元1～M。當(dāng)劃分?jǐn)?shù)M足夠大時(shí)，可認(rèn)為腔鏡1上每個(gè)單元的復(fù)振幅起伏不大，即復(fù)振幅均勻分布。在這種情況下，可用復(fù)振幅U1［m］來代替腔鏡1上第m個(gè)單元的復(fù)振幅分布。同理，將腔鏡2劃分為單元1～N。對于腔鏡 2上的第n個(gè)單元，其復(fù)振幅U2［n］可認(rèn)為是腔鏡1上單元1～M作用并疊加的結(jié)果:

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

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

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

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

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

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

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

腔鏡上的光場在腔內(nèi)一次往返后，其復(fù)振幅分布Ui'可表示為:

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

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

2 實(shí)例分析及討論

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

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

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

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

3 單元格數(shù)的合理劃分

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

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

3.1 單元格數(shù)的合理劃分

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

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

從而不難得出合理單元格數(shù)Gb與相對誤差d的關(guān)系:

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

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

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

3.2 實(shí)例計(jì)算與分析

對于一個(gè)腔長為0.6m、兩腔鏡的孔徑均為1mm、激光波長為632.8nm的諧振腔，要求計(jì)算的相對誤差d=0.1%。根據(jù)(16)式計(jì)算得合理單元格數(shù)Gb應(yīng)為26.91。分別將單元?jiǎng)澐指駭?shù)G設(shè)置為10，20，25，30，35和40，計(jì)算相應(yīng)本征值絕對值的最大值(見表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可見，當(dāng)單元格數(shù)G達(dá)到合理單元格數(shù)Gb后再增加單元格數(shù)計(jì)算的相對誤差分別為0.08%，0.04%和0.03%，均小于所要求的計(jì)算相對誤差0.1%。同樣不難從圖4可以發(fā)現(xiàn)，單元格數(shù)G達(dá)到合理劃分格數(shù)Gb后，TEM01模的光場分布變化也較小，基本趨于穩(wěn)定。

4 結(jié)論

從基爾霍夫-菲涅耳衍射積分方程出發(fā)，采用FEM劃分的思想對腔鏡進(jìn)行徑向單元?jiǎng)澐?，并將單元?jié)點(diǎn)處的光場復(fù)振幅變量分離出來，可一次求解多個(gè)模式，與初值的選擇無關(guān)，適用于復(fù)雜腔型的計(jì)算。對腔鏡單元的合理劃分進(jìn)行了有益的探索，結(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，F(xiàn)AN 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.