《1 引言》

1 引言

光通信的迅速发展促使波导光电子器件的需求不断激增,掩埋矩形光波导及脊形光波导是构成波导光电子器件的基本结构,其理论分析有助于优化器件结构。当光波导芯层与包层折射率反差较大时,标量本征模不能反映模场本质,而必须求解矢量本征模;由于这两类光波导所承载的矢量本征模没有解析解,近似方法如Mareatilli法^[1]及有效折射法(EIM)^[2]只适合于求解标量本征模,因此数值方法是求解光波导矢量本征模不可缺少的工具,常用的数值方法主要有有限差分法(FDM)^[3]、有限元法(FEM)^[4]等。但这两类方法导致计算大矩阵,对计算机内存要求高,求解时间长。级数展开法(SEM)将待求电场履成级数形式,如正弦级数^[5]、厄米-高斯级数^[6]等,最终归结为求解标准本征值方程,待求矩阵小,计算效率高,广泛用来求解光波导的标量及矢量本征模问题^[5~7]。

由于光波导是开放式结构,其承载的本征模场应消逝在无穷远处,而实际计算时只能取有限计算窗口,由此产生边界截断,引起非物理反射,影响计算精度。变量变换将无限z,y平面映射成单位平面,使单位平面边界上的电磁场自然为零,等效于自然边界条件,因而避免了边界截断问题,非物理反射自动消除,计算精度因而得到提高,目前仅应用于分析光波导的标量本征模^[8]。笔者运用变量变换级数展开法(VTSEM)求解掩埋矩形光波导及脊形光波导所承载的矢量本征模,获得矢量本征模场分布及其传播常数,供优化器件结构参考。

《2 方法描述》

2 方法描述

单色光电场矢量波方程可简化为偏微分方程:

\(\frac{\partial^{2} u_{x}}{\partial x^{2}}+\frac{\partial^{2} u_{x}}{\partial y^{2}}+k_{0}^{2}\left(n^{2}-\bar{n}^{2}\right) u_{x}+ \\ 2 \frac{\partial}{\partial x}\left[u_{x} \frac{\partial \ln (n)}{\partial x}+u_{y} \frac{\partial \ln (n)}{\partial y}\right]=0\) , (1a)

\(\frac{\partial^{2} u_{y}}{\partial x^{2}}+\frac{\partial^{2} u_{y}}{\partial y^{2}}+k_{0}^{2}\left(n^{2}-\bar{n}^{2}\right) u_{y}+ \\ 2 \frac{\partial}{\partial y}\left[u_{x} \frac{\partial \ln (n)}{\partial x}+u_{y} \frac{\partial \ln (n)}{\partial y}\right]=0\) ,   (1b)

式中k₀为自由空间波数,\(\bar{n}\)为有效折射率,β=k_0\(\bar{n}\)即为传播常数,n(x,y)描述了光波导横截面的折射率分布，u_x , u_y为x, y方向的电场分量。引入如下变量变换：

\(x=\sigma_{x} \tan [\pi(\xi-0.5)]\),   (2a)

\(y=\sigma_{y} \tan [\pi(\eta-0.5)]\),    (2b)

其中σ_x,σ_y为转换因子，无限平面x∈[-∞, + ∞],   y∈ [-∞, ＋∞]则映射成单位平面ξ∈[0, 1], η ∈[0, 1],方程组(la)、(1b)变为

\(\left(\frac{\mathrm{d} \xi}{\mathrm{d} x}\right)^{2} \frac{\partial^{2} u_{x}}{\partial \xi^{2}}+\frac{\mathrm{d}^{2} \xi}{\mathrm{d} x^{2}} \frac{\partial u_{x}}{\partial \xi}+\left(\frac{\mathrm{d} \eta}{\mathrm{d} y}\right)^{2} \frac{\partial^{2} u_{x}}{\partial \eta^{2}}+ \\ \frac{\mathrm{d}^{2} \eta}{\mathrm{d} y^{2}} \frac{\partial u_{x}}{\partial y}+k_{0}^{2}\left(n^{2}-\bar{n}^{2}\right) u_{x}+ \\ 2\left[\left(\frac{\mathrm{d} \xi}{\mathrm{d} x}\right)^{2} \frac{\partial}{\partial \xi}\left(u_{x} \frac{\partial \ln (n)}{\partial \xi}\right)+\frac{\mathrm{d}^{2} \xi}{\mathrm{d} x^{2}} u_{x} \frac{\partial \ln (n)}{\partial \xi}+\right. \\ \left.\frac{\mathrm{d} \xi}{\mathrm{d} x} \frac{\mathrm{d} \eta}{\mathrm{d} y} \frac{\partial}{\partial \xi}\left(u_{y} \frac{\partial \ln (n)}{\partial \eta}\right)\right]=0\) ,   (3a)

\(\left(\frac{\mathrm{d} \xi}{\mathrm{d} x}\right)^{2} \frac{\partial^{2} u_{y}}{\partial \xi^{2}}+\frac{\mathrm{d}^{2} \xi}{\mathrm{d} x^{2}} \frac{\partial u_{y}}{\partial \xi}+\left(\frac{\mathrm{d} \eta}{\mathrm{d} y}\right)^{2} \frac{\partial^{2} u_{y}}{\partial \eta^{2}}+ \\ \frac{\mathrm{d}^{2} \eta}{\mathrm{d} y^{2}} \frac{\partial u_{y}}{\partial y}+k_{0}^{2}\left(n^{2}-\bar{n}^{2}\right) u_{y}+ \\ 2\left[\left(\frac{\mathrm{d} \eta}{\mathrm{d} y}\right)^{2} \frac{\partial}{\partial \eta}\left(u_{y} \frac{\partial \ln (n)}{\partial \eta}\right)+\frac{\mathrm{d}^{2} \eta}{\mathrm{d} y^{2}} u_{y} \frac{\partial \ln (n)}{\partial \eta}+\right. \\ \left.\frac{\mathrm{d} \xi}{\mathrm{d} x} \frac{\mathrm{d} \eta}{\mathrm{d} y} \frac{\partial}{\partial \eta}\left(u_{x} \frac{\partial \ln (n)}{\partial \xi}\right)\right]=0\),    (3b)

将电场u_x(ξ, η), u_y(ξ, η)在区间ξ∈[0, 1], η ∈[0, 1]展成正弦级数，为正交完全集：

\(u_{x}(\xi, \eta)=\sum_{k=1}^{N_{x} \times N_{y}} c_{k}^{x} \varphi_{k}(\xi, \eta)= \\ \sum_{p=1}^{N_{x}} \sum_{q=1}^{N_{y}} c_{p q}^{x} 2 \sin (p \pi \xi) \sin (q \pi \eta)\),     (4a)

\(u_{y}(\xi, \eta)=\sum_{k=1}^{N_{x} \times N_{y}} c_{k}^{y} \varphi_{k}(\xi, \eta)= \\ \sum_{p=1}^{N_{x}} \sum_{q=1}^{N_{y}} c_{p q}^{y} 2 \sin (p \pi \xi) \sin (q \pi \eta)\),     (4b)

其中正整数N_x,N_y分别为ξ及η方向的级数项,c_k^x ,c_k^y分别为u_x,u_y的展开系数,且有如下关系

\(\varphi_{k}=2 \sin (p \pi \xi) \sin (q \pi \eta)\),    (5)

\(p=\operatorname{ent}\left(\frac{k-1}{N_{y}}\right)+1\),               (6a)

\(q=\bmod \left(k-1, N_{y}\right)+1\),    (6b)

式中ent表示取整,mod表示求余。将方程式(4a)、式(4b)分别代人方程式(3a)、式(3b),并应用伽辽金法,得到如下矩阵方程组

\(\begin{aligned} \boldsymbol{A}^{x} \boldsymbol{C}^{x}+\boldsymbol{B}^{y} \boldsymbol{C}^{y} &=\beta^{2} \boldsymbol{C}^{x} \end{aligned}\),     (7a)

\(\begin{aligned} \boldsymbol{A}^{y} \boldsymbol{C}^{y}+\boldsymbol{B}^{x} \boldsymbol{C}^{x} &=\beta^{2} \boldsymbol{C}^{y} \end{aligned}\)。   (c7b)

式中各矩阵元为

\(A_{l, k}^{x}=\int_{\xi=0}^{1} \int_{\eta=0}^{1}\left[\left(\frac{\mathrm{d} \xi}{\mathrm{d} x}\right)^{2} \frac{\partial^{2} \varphi_{k}}{\partial \xi^{2}}+\frac{\mathrm{d}^{2} \xi}{\mathrm{d} x^{2}} \frac{\partial \varphi_{k}}{\partial \xi}+\right. \\ \left.\left(\frac{\mathrm{d} \eta}{\mathrm{d} y}\right)^{2} \frac{\partial^{2} \varphi_{k}}{\partial \eta^{2}}+\frac{\mathrm{d}^{2} \eta}{\mathrm{d} y^{2}} \frac{\partial \varphi_{k}}{\partial y}+k_{0}^{2} n^{2} \varphi_{k}\right] \varphi_{l} \mathrm{~d} \xi \mathrm{d} \eta+ \\ \int_{\xi=0}^{1} \int_{\eta=0}^{1} 2\left[\left(\frac{\mathrm{d} \xi}{\mathrm{d} x}\right)^{2} \frac{\partial}{\partial \xi}\left(\varphi_{k} \frac{\partial \ln (n)}{\partial \xi}\right)+\right. \\ \left.\frac{\mathrm{d}^{2} \xi}{\mathrm{d} x^{2}} \varphi_{k} \frac{\partial \ln (n)}{\partial \xi}\right] \varphi_{l} \mathrm{~d} \xi \mathrm{d} \eta\)，     （8a）

\(B_{l, k}^{y}=\int_{\xi=0}^{1} \int_{\eta=0}^{1} 2\left[\frac{\mathrm{d} \xi}{\mathrm{d} x} \frac{\mathrm{d} \eta}{\mathrm{d} y} .\right. \\ \left.\frac{\partial}{\partial \xi}\left(\varphi_{k} \frac{\partial \ln (n)}{\partial \eta}\right)\right] \varphi_{l} \mathrm{~d} \xi \mathrm{d} \eta\),                                         (8b)

\(\begin{array}{c} A_{l^{\prime}, k^{\prime}}^{y}=\int_{\xi=0}^{1} \int_{\eta=0}^{1}\left[\left(\frac{\mathrm{d} \xi}{\mathrm{d} x}\right)^{2} \frac{\partial^{2} \varphi_{k^{\prime}}}{\partial \xi^{2}}+\frac{\mathrm{d}^{2} \xi}{\mathrm{d} x^{2}} \frac{\partial \varphi_{k^{\prime}}}{\partial \xi}+\right. \\ \left.\left(\frac{\mathrm{d} \eta}{\mathrm{d} y}\right)^{2} \frac{\partial^{2} \varphi_{k^{\prime}}}{\partial \eta^{2}}+\frac{\mathrm{d}^{2} \eta}{\mathrm{d} y^{2}} \frac{\partial \varphi_{k^{\prime}}}{\partial y}+k_{0}^{2} n^{2} \varphi_{k^{\prime}}\right] \varphi_{l^{\prime}} \mathrm{d} \xi \mathrm{d} \eta+ \\ \int_{\xi=0}^{1} \int_{\eta=0}^{1} 2\left[\left(\frac{\mathrm{d} \eta}{\mathrm{d} y}\right)^{2} \frac{\partial}{\partial \eta}\left(\varphi_{k^{\prime}} \frac{\partial \ln (n)}{\partial \eta}\right)+\right. \\ \left.\frac{\mathrm{d}^{2} \eta}{\mathrm{d} y^{2}} \varphi_{k^{\prime}} \frac{\partial \ln (n)}{\partial \eta}\right] \varphi_{l^{\prime}} \mathrm{d} \xi \mathrm{d} \eta \end{array}\),         (8c)

\(\begin{array}{c} B_{l^{\prime}, k^{\prime}}^{x}=\int_{\xi=0}^{1} \int_{\eta=0}^{1} 2\left[\frac{\mathrm{d} \xi}{\mathrm{d} x} \frac{\mathrm{d} \eta}{\mathrm{d} y} \cdot\right. \\ \left.\frac{\partial}{\partial \eta}\left(\varphi_{k^{\prime}} \frac{\partial \ln (n)}{\partial \xi}\right)\right] \varphi_{l^{\prime}} \mathrm{d} \xi \mathrm{d} \eta \\ \end{array}\) ,                                        (8d)

\(\beta^{2}=k_{0}^{2} n^{2} \\ \),                                                                       (8e)

\(\boldsymbol{C}^{x}=\left[c_{1}^{x}, c_{2}^{x}, c_{3}^{x}, \cdots, c_{N_{x} \times N_{y}}^{x}\right]^{\mathrm{T}} \\ \),                                   (8f)

\(\boldsymbol{C}^{y}=\left[c_{1}^{y}, c_{2}^{y}, c_{3}^{y}, \cdots, c_{N_{x} \times N_{y}}^{y}\right]^{\mathrm{T}}, \\ k, l, k^{\prime}, l^{\prime}=1,2,3, \cdots, N_{x} \times N_{y}。\)                                 (8g)

方程式(7a)、式(7b)有相同的本征值,因此可合并成一个矩阵本征值方程

\(\left[\begin{array}{ll} \boldsymbol{A}^{x} & \boldsymbol{B}^{y} \\ \boldsymbol{B}^{x} & \boldsymbol{A}^{y} \end{array}\right]\left[\begin{array}{l} \boldsymbol{C}^{x} \\ \boldsymbol{C}^{y} \end{array}\right]=\beta^{2}\left[\begin{array}{l} \boldsymbol{C}^{x} \\ \boldsymbol{C}^{y} \end{array}\right]\)。                              （9）

因此,获得矩阵A^x、B^y、A^y、B^x的矩阵元后,通过求解方程式(9)的本征值,可得矢量本征模的有效折射率,相应的本征矢即为方程式(4)的展开系数,矢量本征模自然获得。

方程式(8a)及式(8c)中前4个二重积分利用正弦函数的正交性可得其解析表达式^[8];方程式(8a)、式(8b)、式(8c)及式(8d)中含有折射率自然对数的导数项可利用分步积分法消除^[7]。求得x、y方向的电场分量后,由电磁场的各分量关系式可求得电磁场的剩余分量,此时求得的模场为准TE模;方程式(1a)、式(lb)换成磁场方程,可求得准TM模,方法类似。

《3 数值结果》

3 数值结果

方程式(2a)、式(2b)提供的是非线性变换,矢量本征模场将随σ_x,σ_y变化而变化,为此必须优化选择σ_x,σ_y,以使方程式(9)快速收敛并保证矢量本征模场不变形。根据文献[8],如图1及图2所示的掩埋矩形光波导及脊形光波导分别选取σ_x=ω/2, σ_y=h/2即可满足要求。

《图1》

图1 掩埋短形光波导的横向尺寸

Fig.1 Cross section of buried rectangular optical waveguide

《图2》

图2 脊形光波导的横向尺寸

Fig.2 Cross section of rib optical waveguide

按文献[7]的数值计算实例和模型,给出EH₁₁模有效折射率的精确值为1.44795;用伽认金法算出的EH₁₁模有效折射率为1.44862,与精确值相差达0.046%;本文方法为1.448157,与精确值相差仅0.015%,是一种高精度的数值方法。

表1给出了图1所示的掩埋矩形光波导归一化传播常数B随归一化频率,的变化关系,其中

-\(B=\frac{\bar{n}^{2}-n_{\mathrm{s}}^{2}}{n_{\mathrm{c}}^{2}-n_{\mathrm{s}}^{2}}\)     (10)

\(\nu=\frac{h}{\lambda} \sqrt{n_{\mathrm{c}}^{2}-n_{\mathrm{s}}^{2}}\)      (11)

ω/h=2, n_c=1.15pm,级数项N:=Ny=12。本文结果与采用其它方法如非线性追代有限差分法(NLIFDM)^[9]、半矢量有限差分法(SVFDM)^[10]及傅立叶算子变换法(FOTM)^[11]计算的结果吻合较好。但上述方法的计算窗口随归一化频率ν增大,导出的矩阵随之增大,为此必须反复选取计算窗口,非常繁琐。本文的计算窗口固定为单位平面,待求矩阵的大小不变,较为方便。另外,当级数项N_x=N_y=12,方程式(9)导出的矩阵阶数为288x288;采用以上方法,计算窗口必须为芯层区的3倍以上,剖分成144x144点阵(最后也须求解288x288矩阵)所得的结果将非常粗糙,因此必须进一步细分计算窗口,导出的矩阵远大于本文方法的矩阵。

《表1》

表1 掩埋矩形光波导的归一化传播常数B随归一化频率ν的变化关系
Table 1 Variation of the normalized propagation constants B for the buried rectangular optical waveguide with the normalized frequency ν

表2中ν=1.0,其他参数与表1相同。由表可见标量分析得出的归一化传播常数B始终大于矢量分析得出的结果。只有当芯层折射率n_c与包层折射率n_s相近时两者才近似相等。

表2 掩埋矩形光波导的归一化传播常数B随折射率n_s的变化关系
Table 2 Variation of the normalized propagation constants B for the buried rectangular optical wave guide with the refractive index n_s

图3给出了掩埋矩形光波导所承载的矢量基模电场分布,其中ω=4μm,h=2μm,n_c=1.5,n_s=1.45,λ=1.15μm。可见u_x与u_y的模场形状差异很大,u_y在角上畸变。

《图3》

图3 掩埋矩形光波导承载的矢量基模电场分布
Fig.3 Field intensity distributions of vectorial fundamental electric mode supported by a buried rectangular optical waveguide

表3为图2所示的脊形光波导归一化传播常数卫随!的变化关系,其中ω=1.5μm,h=1μm,n₀=1.0,n_c=3.44,n_s=3.40,λ=1.15μm, N_x=N_y=15。由于脊形光波导结构较矩形光波团导复杂,因此级数项较大时才能达到收敛解。由表可见,归一化传播常数B随t增大而增大。本文结果接近有限元法(FEM)获得的结果^[12]。

《表3》

表3 脊形光波的归一化传播常数B随t的变化关系

Table 3 variation of the normalized propagation constants B for the rib optical wave guide with t

图4分别给出了t=0.2μm,t=1.0μm时脊

《图4》

图4 脊形光波导承载的矢量基模电场分布

Fig. 4 Field intensity distributions of vectorial fundamental electric mode supported by a rib optical waveguide 

形光波导所承载的矢量基模电场分布,其中ω=2μm,h=1.3μm,n₀=1.0,n_c=3.44,n_s=3.34,λ=1.55μm。 由图4可见t较大时电场向光波导两侧渗透明显,而t较小时电场向两侧渗透小,较集中在芯区。

《4 结论》

4 结论

运用变量变换级数展开法,获得了掩埋短形光波导及脊形光波导的矢量本征模场分布及其传播常数,并观察到电场在角上的畸变。变量变换使无限平面映射成单位平面,从而使单位平面边界上的电磁场自然为零,等效于自然边界条件,因而避免了边界截断,非物理反射自动消除,提高了计算精度。另外,由于这种方法所导致的矩阵阶数小,有较高的计算效率。本文分析的结果与其它方法的结果吻合较好,可以用来优化波导光电子器件的结构。