光ファイバーとベッセル関数

初めにコアとクラッドからなるステップ型の光ファイバーは比較的に解析的に計算することができます。

Maple を使った光ファイバーの解析例を紹介します。

[Maple Plot] [Maple Plot]

最初に一般的な準備を行っておきます。

>    restart;
with(plots):
with(plottools):
DO:=Pi/180:

Warning, the name changecoords has been redefined

Warning, the assigned name arrow now has a global binding

ステップ型光ファイバーの構造

ステップ型光ファイバーの構造をシリンダー・オブジェクトを使って表示してみましょう。

>    c1:=rotate(cylinder([0,0,1],3,10,color=blue),0,90*DO,0):

>    c2:=rotate(cylinder([0,0,0],1,10,color=red),0,90*DO,0):

>    display({c1,c2},scaling=constrained);

[Maple Plot]

ファイバーの断面構造

通常シングルモードの光ファイバーはコアの屈折率とクラッドの屈折率は近いものが使われます。

コアの屈折率が少し高くなっています。n0はクラッド、n1はコアの屈折率です。aはコアの半径で5ミクロンとします。

>    n0:=1.4443;n1:=1.4487;a:=5;

n0 := 1.4443

n1 := 1.4487

a := 5

ステップ状の屈折率はpiecewiseを使って定義することができます。

>    n:=(r,theta)->piecewise(r<a,n1,n0);

n := proc (r, theta) options operator, arrow; piecewise(r < a,n1,n0) end proc

>    evalf(n(6,2));

1.4443

断面の屈折率分布をプロットしてみましょう。

>    plot3d(n(sqrt(x^2+y^2),arctan(y,x)),x=-10..10,y=-10..10,grid=[60,60]);

[Maple Plot]

ベッセル関数

このような光ファイバーを光が伝搬するとき、特定のモードと呼ばれる光フィールドすなわち電界と

磁界のベクトル場として伝わります。もっとも次数の少ないモードを基本モードと呼びます。

ここで扱う単純な構造の光ファイバーで、コアとクラッドの屈折率が近接している場合シングルモード

と呼ばれこの基本モードのみ伝搬します。この基本モードの電界Exに対する方程式は次のようになります。

>    1/r*Diff(r*diff(E[x](r),r),r)+(k^2*n^2-beta^2)*E[x](r)=0;

1/r*Diff(r*diff(E[x](r),r),r)+(k^2*n^2-beta^2)*E[x](r) = 0

ここで k は波数、n(r) は屈折率分布、 beta  は伝搬定数と呼ばれます。

コアとクラッドの屈折率が近接していることを仮定したモードの近似を LP モードと呼びます。

この方程式の解はベッセル関数で与えられます。Maple にはベッセル関数が用意されていますのでまず J0, J1, J2 のベッセル関数をプロットしてみましょう。

>    J0:=x->BesselJ(0,x);J1:=x->BesselJ(1,x);J2:=x->BesselJ(2,x);

J0 := proc (x) options operator, arrow; BesselJ(0,x) end proc

J1 := proc (x) options operator, arrow; BesselJ(1,x) end proc

J2 := proc (x) options operator, arrow; BesselJ(2,x) end proc

>    plot({J0(x),J1(x),J2(x)},x=0..10);

[Maple Plot]

同様に K0, K1, K2 もプロットしてみましょう。

>    K0:=x->BesselK(0,x);K1:=x->BesselK(1,x);K2:=x->BesselK(2,x);

K0 := proc (x) options operator, arrow; BesselK(0,x) end proc

K1 := proc (x) options operator, arrow; BesselK(1,x) end proc

K2 := proc (x) options operator, arrow; BesselK(2,x) end proc

>    plot({K0(x),K1(x),K2(x)},x=0..5,0..5);

[Maple Plot]

分散方程式と伝搬定数

波長を lambda  とすると波数は次のようになります。実際的に波長を1.55ミクロンとしましょう。

>    k:=2*Pi/lambda;

k := 2*Pi/lambda

>    lambda:=1.55;

lambda := 1.55

Delta  をつぎのように定義しておきます。

>    Delta:=(n1^2-n0^2)/(2*n1^2);

Delta := .3032593461e-2

kappa = sqrt(k^2*n1^2-beta^2) , sigma = sqrt(beta-k^2*n0^2)  と置くと Ex に関する波動方程式は次のようになります。

>    eq_1:=1/r*diff(r*diff(g(r),r),r)+kappa*g(r)=0;

eq_1 := 1/r*(diff(g(r),r)+r*diff(g(r),`$`(r,2)))+kappa*g(r) = 0

>    dsolve(eq_1,g(r));

g(r) = _C1*BesselJ(0,kappa^(1/2)*r)+_C2*BesselY(0,kappa^(1/2)*r)

このように解は 0 次ベッセル関数 J0 で表されます。どうようにEyに関しても解けます。

>    kappa:=sqrt(k^2*n1^2-beta^2);

kappa := (3.494246312*Pi^2-beta^2)^(1/2)

>    sigma:=sqrt(beta-k^2*n0^2);

sigma := (beta-3.473053055*Pi^2)^(1/2)

u = kappa*a , w = sigma*a  と置くと、伝搬定数 beta  は次のようになります。u と w の間には

>    sqrt(u^2+w^2)=kappa*a*sqrt(n_1^2-n_0^2);

(u^2+w^2)^(1/2) = 5*(3.494246312*Pi^2-beta^2)^(1/2)*(n_1^2-n_0^2)^(1/2)

の関係がなりたち、これを v とおきます。f を光の周波数とすると

>    v:=2*Pi*f*a/c*sqrt(n_1^2-n_0^2);

v := 10*Pi*f/c*(n_1^2-n_0^2)^(1/2)

の量で規格化周波数と呼びます。

>    v:='v';

v := 'v'

>    beta:=k*n1*sqrt(1-2*Delta*(u/v)^2);

beta := 1.869290323*Pi*(1-.6065186922e-2*u^2/v^2)^(1/2)

b を ((beta/k)^2-n0^2)/(n1^2-n0^2)  とすると分散方程式は次のようになります。

>    b:='b';

b := 'b'

>    eq:=J0(v*sqrt(1-b))/(sqrt(1-b)*J1(v*sqrt(1-b)))=K0(v*sqrt(b))/(sqrt(b)*K1(v*sqrt(b)));

eq := BesselJ(0,v*(1-b)^(1/2))/(1-b)^(1/2)/BesselJ(1,v*(1-b)^(1/2)) = BesselK(0,v*b^(1/2))/b^(1/2)/BesselK(1,v*b^(1/2))

v を 2 としたとき、この式を b に関して fsolve で数値解を解くと b が求められます。

>    v:=2;

v := 2

>    u:=v*sqrt(1-b);

u := 2*(1-b)^(1/2)

>    b:=fsolve(eq,b=0.5);

b := .4161633927

>    u:=v*sqrt(1-b);

u := 1.528184030

>    w:=v*sqrt(b);

w := 1.290214544

b が決まり、u, w が定まり伝搬定数 beta  も決まります。

>    beta;

1.865977736*Pi

分散曲線

v を変化させるとbが変化します。これを分散曲線と呼びます。このグラフをプロットしてみましょう。

>    m:=25:pp:=NULL:dv:=(5-0.4)/m:v:=0.4:b:='b':

>    for k to m do
v:=evalf(v+dv);eq:=J0(v*sqrt(1-b))/(sqrt(1-b)*J1(v*sqrt(1-b)))=K0(v*sqrt(b))/(sqrt(b)*K1(v*sqrt(b))):
bb:=fsolve(eq,b=0.5);
pp:=pp,[v,bb];
od:

>    listplot([pp]);

[Maple Plot]

電界の表示

>    A:=1.0;

A := 1.0

>    q:=u^2*w^2/(2*v^2)*J0(u)/(u*J1(u));

q := .4492989200e-1

A は光の大きさで振幅に相当する定数です。q を上のように定義すると電界の Ex, Ey はそれぞれ次のようになります。

>    Ex:=(r,theta)->A*(J0(u/a*r)-q*Delta*J2(u/a*r)*cos(2*theta));

Ex := proc (r, theta) options operator, arrow; A*(J0(u/a*r)-q*Delta*J2(u/a*r)*cos(2*theta)) end proc

>    Ey:=(r,theta)->-A*q*Delta*J2(u/a*r)*sin(2*theta);

Ey := proc (r, theta) options operator, arrow; -A*q*Delta*J2(u/a*r)*sin(2*theta) end proc

これをベクトル場表示してみましょう。コアの径がわかるように円を同時に表示します。

この電界はコア内部の電界を表しているのでコアの外側は正しいベクトル場ではありません。

>    p1:=fieldplot([Ex(sqrt(x^2+y^2),arctan(y,x)),Ey(sqrt(x^2+y^2),arctan(y,x))],x=-10..10,y=-10..10):

>    p2:=implicitplot(x^2+y^2=a^2,x=-10..10,y=-10..10):

>    display({p1,p2},scaling=constrained);

[Maple Plot]

piecewise 関数を使ってコア内部とクラッドの電界を関数で定義することができます。

>    Ex:=(r,theta)->piecewise(r<a,A*J0(u/a*r),A*J0(u)/K0(w)*K0(w/a*r));

Ex := proc (r, theta) options operator, arrow; piecewise(r < a,A*J0(u/a*r),A*J0(u)/K0(w)*K0(w/a*r)) end proc

これをプロットしてみましょう。コアの径が分かるようにコア径を tubeplot を使ってリング状に表示します。

>    p1:=plot3d(Ex(sqrt(x^2+y^2),arctan(y,x)),x=-10..10,y=-10..10,shading=zhue):

>    p2:=tubeplot([0,0,t],t=0.5..0.53,radius=a,tubepoints=20):

>    display({p1,p2});

[Maple Plot]

>   

参考文献:

光導波路の基礎   日本電信電話(株)岡本勝就  コロナ社(1992)フォトニクスシリーズ13