Part 1: Waveguide Basic Concept

ECE 5358/6358                                 han q le - copyrighted
Use solely for students registered for UH ECE 6358/5368 during courses - DO NOT distributed (copyrighted materials)

Load package

TE mode calc

IndexprofSSW[x_,nindex_List,a_]:=If[x≤ -a,nindex[[1]],If[x≤ a,nindex[[2]],nindex[[1]] ] ]  ;
ResolvTE0[nindex_List,a_,λ_,neff_]:= Module[{ϵ1,ϵ2,κ1, k2,k0},
		{ϵ1,ϵ2}=nindex^2 ; k0=2*Pi/λ;
    	κ1=k0*Sqrt[neff^2-ϵ1  ]; k2=k0*Sqrt[ϵ2-neff^2 ] ;
			{ (Cos[a k2] κ1-Sin[a k2] k2) 
,(Cos[a k2] k2+Sin[a k2] κ1) 
,(Cos[a k2] κ1-Sin[a k2] k2) (Cos[a k2] k2+Sin[a k2] κ1)}  ] ;
ResolvTE[nindex_List,a_,λ_,neff_]:= Module[{ϵ1,ϵ2,κ1, k2,k0},
		{ϵ1,ϵ2}=nindex^2 ; k0=2*Pi/λ;
    	κ1=k0*Sqrt[neff^2-ϵ1  ]; k2=k0*Sqrt[ϵ2-neff^2 ] ;
			(Cos[a k2] κ1-Sin[a k2] k2) (Cos[a k2] k2+Sin[a k2] κ1) ] ;

FieldTE[nindex_List,a_,λ_,neff_,x_]:= Module[{ϵ1,ϵ2,κ1, k2,k0,Norm},
		{ϵ1,ϵ2}=nindex^2 ; k0=2*Pi/λ;
    	κ1=k0*Sqrt[Abs[neff^2-ϵ1 ] ]; k2=k0*Sqrt[Abs[ϵ2-neff^2] ] ;
		Norm=(((-cos(4 a k2) κ1+κ1+2 a ((k2^2)+((κ1)^2))))/(2 (k2^2)))+
			(((k2-κ1) (k2+κ1) sin(4 a k2))/(4 (k2^3)))+(1/(κ1)) ;
(1/(Sqrt[Norm])) If[Abs[x]<a,{ cos((a+x)k2)+((sin((a+x) k2) κ1)/k2)          (* this is Ey  *)
,-I ( cos((a+x) k2) κ1-sin((a+x) k2) k2)/k0   (* this is Hz  *)
,- neff ( cos((a+x)k2)+((sin((a+x) k2) κ1)/k2))}   (* this is Hx  *)
	, If [x<0, {(E^(κ1 (x+a) ))                                      (* this is Ey  *) 
,-I κ1 (E^(κ1 (x+a) ))/k0                        (* this is Hz  *)
 ,- neff  (E^(κ1 (x+a) ))}                            (* this is Hx  *)
, ((sin(2 a k2)  ((k2^2)+((κ1)^2)))/(2 k2 κ1))  {(E^(-κ1 (x-a)     ))(* this is Ey  *) 
, I κ1 (E^(-κ1 (x-a) ))/k0  (* this is Hz  *)
,- neff  (E^(-κ1 (x-a) ))}  (* this is Hx  *) ] 
]
] ;

ParmTE[nindex_List,a_,λ_,neff_]:= Module[{ϵ1,ϵ2,κ1, k2,k0,Norm},
		{ϵ1,ϵ2}=nindex^2 ; k0=2*Pi/λ;
    	κ1=k0*Sqrt[Abs[neff^2-ϵ1 ] ]; k2=k0*Sqrt[Abs[ϵ2-neff^2] ] ;
		Norm=(((-cos(4 a k2) κ1+κ1+2 a (k0^2) (ϵ2-ϵ1) ))/(2 (k2^2)))
			+(((k2-κ1) (k2+κ1) sin(4 a k2))/(4 (k2^3)))+(1/(κ1)) ;
{a, neff, κ1, k2,k0,Norm}  ] ;
FieldTE2[parm_List,x_]:= Module[{a, neff, κ1, k2,k0,Norm} ,
{a, neff, κ1, k2,k0,Norm}= parm  ; 
(1/(Sqrt[Norm])) If[Abs[x]<a, 
		{ cos((a+x)k2)+((sin((a+x) k2) κ1)/k2)          (* this is Ey  *)
				,-I ( cos((a+x) k2) κ1-sin((a+x) k2) k2)/k0   (* this is Hz  *)
				,- neff ( cos((a+x)k2)+((sin((a+x) k2) κ1)/k2))}   (* this is Hx  *)
	, If [x<0, {(E^(κ1 (x+a) ))                                        (* this is Ey  *)
					          ,-I κ1 (E^(κ1 (x+a) ))/k0                   (* this is Hz  *)
					            ,- neff  (E^(κ1 (x+a) ))}                     (* this is Hx  *)
				 , ((sin(2 a k2)  ((k2^2)+((κ1)^2)))/(2 k2 κ1))  {(E^(-κ1 (x-a)     ))(* this is Ey  *)
					                                               , I κ1 (E^(-κ1 (x-a) ))/k0  (* this is Hz  *)
						                                ,- neff  (E^(-κ1 (x-a) ))}  (* this is Hx  *)
			 ] ]
					
		      ] ;

TM mode calc

ResolvTM[nindex_List,a_,λ_,neff_]:= Module[{ϵ1,ϵ2,κ1, k2,k0},
		{ϵ1,ϵ2}=nindex^2 ; k0=2*Pi/λ;
    	κ1=k0*Sqrt[Abs[neff^2-ϵ1 ] ];
		 k2=k0*Sqrt[Abs[ϵ2-neff^2] ] ;
			 (cos(a k2) ϵ2 κ1-sin(a k2) k2 ϵ1) (cos(a k2) k2 ϵ1+sin(a k2) ϵ2 κ1)  ] ;

FieldTM[nindex_List,a_,λ_,neff_,x_]:= Module[{ϵ1,ϵ2,κ1, k2,k0,Norm},
		{ϵ1,ϵ2}=nindex^2 ; k0=2*Pi/λ;
    	κ1=k0*Sqrt[Abs[neff^2-ϵ1 ] ]; k2=k0*Sqrt[Abs[ϵ2-neff^2] ] ;
Norm=((((neff/(ϵ2))))^2)((1/(2 (k2^2) ((ϵ1)^2))) ϵ1 ϵ2 κ1-ϵ1 ϵ2 Cos[4 a k2] κ1+2 a ((k2^2) ((ϵ1)^2)+((ϵ2)^2) ((κ1)^2))  
+(1/(4 (k2^3) ((ϵ1)^2)))(k2 ϵ1-ϵ2 κ1) (k2 ϵ1+ϵ2 κ1) Sin[4 a k2])
+((((neff/(ϵ1))))^2)(1/( κ1))  ;
If[Abs[x]≤ a,{ Cos[(a+x)k2]+((Sin[(a+x) k2]  ϵ2  κ1)/(k2 ϵ1)) 
,I  (((Cos[(a+x)k2] ϵ2  κ1)/(ϵ1))-Sin[(a+x) k2]k2)/(ϵ2 k0)
,neff (Cos[(a+x)k2]+((Sin[(a+x) k2] ϵ2  κ1)/(k2 ϵ1)))/ϵ2 }
, If [x<0, {(E^(κ1 (x+a) )),I κ1 (E^(κ1 (x+a) ))/(ϵ1*k0 ) , neff  (E^(κ1 (x+a) ))/ϵ1 }
,{(E^(-κ1 (x-a) )), -I κ1 (E^(-κ1 (x-a) ))/(ϵ1*k0 ) , neff  (E^(-κ1 (x-a) ))/ϵ1   }]]
]

Segment 1

0. Introduction

Waveguides are simply devices that we use to spatially confine light for transportation and processing, because spatially freely propagating and diffracting light is generally more difficult to control and process for applications.

1. Review: evanescent waves

1.1 Evanescent wave revisited

We have seen this before:

Animate[λ = 1; k0=((2 π)/λ); ω =2 π ;k =((2 π * 1.2)/λ);
 κ=(Sqrt[(k^2)-(k0^2) ]) ;
Plot3D[Re[ (E^(I (k x - ω t)-κ z))] ,{x,0,5},{z,0,2}
,BoxRatios->{5,2,2},PlotPoints->{50,10},ViewPoint->{-3,1,2}, Mesh->False, PlotRange->{-1,1}
	, ImageSize -> {300,200}] , {t,0,1}, AnimationRunning-> False]

What does this wave mean?

Light wave power does not flow in the evanescent direction but only parallel to it.

In real world situation, no medium is infinite and an evanescent wave may eventually be in contact with another medium that allows it to propagate with power (or radiate). But if the exponentially decaying wave amplitude is insignificantly small at the other side of the medium, it can be approximated as a completely evanescent wave in a semi-infinite medium for all practical purposes.

1.2 Evanescent wave in total internal reflection

Where can we find example of evanescent wave? Essentially, wherever we find a field wavefunction f[ξ] that is f[ξ]-> 0 as ξ -> ± ∞ (without single-direction propagation component (e^(i k ξ)) - Dual direction such as distributed Bragg reflector or photonic crystal is OK). Most often, this function f[ξ] is a purely exponential decaying term (e^(-κ ξ)). The evanescent wave can be in one or more spatial coordinate ξ. For example, (e^(-(Subscript[κ, 1]) x -(Subscript[κ, 2]) y)) would indicate evanescence in both x and y dimension.  It means that light does NOT penetrate into the medium in those dimension, or that it is “trapped” in the remaining freely propagating dimension. The evanescent wave can have other form, for example,  (e^(-a (ξ^2))) for parabolic graded index media, for example.

This evanescent effect can happen under many circumstances, from as familar as total internal reflection to periodic structures such as photonic bandgap crystals.

Let’s go back to the two-media reflection and refraction. Sometimes, we have cases like the below:

Let’s take a look at the problem again for TE mode. We already obtained the solution last LectSet4.
Summary of the solution: (see LectSet4):
                        (Subscript[E, 1])= (Subscript[E, inc]) +(Subscript[E, refl])                                     4(2.1.7)  (1.2.1a)
where:             (Subscript[E, inc])= (OverHat[y])    (e^(i  (Subscript[k, 1]) x )) (e^(i  (β z -ω t )))                       4(2.1.8)     (1.2.1b)
and:                 (Subscript[E, refl])= (OverHat[y]) r  (e^(-i  (Subscript[k, 1]) x ))(e^(i  (β z -ω t )))                    4(2.2.12)       (1.2.1c)
                        (Subscript[E, 2])=    (OverHat[y])  t  (e^(i  (Subscript[k, 2]) x ))(e^(i  (β z -ω t )))                     4(2.2.13)     (1.2.2)

where:             β=((2 π)/λ)(Sqrt[(Subscript[ϵ, 1]) (Subscript[μ, 1])]) Sin[θ] =(Subscript[k, 0]) (Subscript[n, 1]) Sin[θ]     4(1.1.1b)  (1.2.3)
                        (Subscript[k, 1])=(Sqrt[(Subscript[ϵ, 1]) (Subscript[μ, 1]) ((Subscript[k, 0])^2) -(β^2) ])=(Subscript[k, 0]) (Subscript[n, 1]) Cos[θ]               (1.2.4) 
Note that the difference now is that we assume medium 1 has index of refraction (Subscript[n, 1]) (and no gain or absorption) instead of 1 for air. But that is all. The rest is the same.
And:               (Subscript[k, 2])= (Sqrt[(Subscript[ϵ, 2]) (Subscript[μ, 2]) ((Subscript[k, 0])^2) -(β^2) ])= (Subscript[k, 0])(Sqrt[(Subscript[ϵ, 2]) (Subscript[μ, 2])-((((Subscript[n, 1]) Sin[θ]))^2) ])      (1.2.5)
The coefficients are:     r=(((Subscript[k, 1])-(Subscript[k, 2]))/((Subscript[k, 1])+(Subscript[k, 2])))   and   t=((2(Subscript[k, 1]))/((Subscript[k, 1])+(Subscript[k, 2])))                    (1.2.6)

We can plot the solutions:

Manipulate[   ω=1;  β=n1 Sin[θ]  ;
k1=n1 Cos[θ] ;k2=(Sqrt[(((n2+I κ2))^2) -(n1^2) ((Sin[θ])^2) ])  ;
w1=Min[8 Cos[θ],3] ; w2=(w1/(Cos[θ]))Re[k2]   ;
dGH= Re[I  Sin[θ]/(π k2)] ;
r=  ((k1-k2)/(k1+k2)) ; tr=   ((2 k1)/(k1+k2)) ;

U=DensityPlot[ If[y≥0, 
Re[ UnitBox[ (x  Cos[θ]+y Sin[θ])/w1] * Exp[I*2 π (x  β-y k1 -ω t)]
                   +If[Re[k2]>0,
UnitBox[(x  Cos[θ]-y Sin[θ])/w1] ,UnitBox[(x  Cos[θ]-y Sin[θ]-dGH)/w1]]
* r  Exp[I*2 π (x  β+y k1 -ω t ) ]]
,If[Re[k2]>0,
UnitBox[ (x  Re[k2]+y n1 Sin[θ])/w2 ] ,  UnitBox[ (x-0.5*dGH/Cos[θ]) Cos[θ]/w1 ]]*Re[tr* Exp[I*2 π (x  β-y k2 -ω t)] ]
] ,{x,-4,4},{y,-2.5,5.5}
,PlotPoints->{60,80},PlotRange->{{-4,4},{-2.5,5.5}, All},Mesh-> False, ColorFunction->"Rainbow"] ;
Show[U,Graphics[{White,Opacity[0.25],Line[{{-4,0},{4,0}}] }] ]
, {t,0,1}, {θ, 0,80*Pi/180} ,{{n1,2},1.33,2},{{n2,1},1,1.5},{κ2,0,0.1}]

Note what happens when θ is larger than some value: threre is no wave in medium 2. Let’s zoom in the region near the boundary.

Manipulate[  ω=1;  β=n1 Sin[θ]  ;
k1=n1 Cos[θ] ;k2=(Sqrt[(((n2+I κ2))^2) -(n1^2) ((Sin[θ])^2) ])  ;
w1=Min[8 Cos[θ],3] ; w2=(w1/(Cos[θ]))Re[k2]   ;
dGH= Re[I  Sin[θ]/(π k2)] ;
r=  ((k1-k2)/(k1+k2)) ; tr=   ((2 k1)/(k1+k2)) ;

U=DensityPlot[ If[y≥0, 
Re[ UnitBox[ (x  Cos[θ]+y Sin[θ])/w1] * Exp[I*2 π (x  β-y k1 -ω t)]
                   +If[Re[k2]>0,
UnitBox[(x  Cos[θ]-y Sin[θ])/w1] ,UnitBox[(x  Cos[θ]-y Sin[θ]-dGH)/w1]]
* r  Exp[I*2 π (x  β+y k1 -ω t ) ]]
,If[Re[k2]>0,
UnitBox[ (x  Re[k2]+y n1 Sin[θ])/w2 ] ,  UnitBox[ (x-0.5*dGH/Cos[θ]) Cos[θ]/w1 ]]*Re[tr* Exp[I*2 π (x  β-y k2 -ω t)] ]
] ,{x,-3,3},{y,-0.8,1.2}, AspectRatio->2/6.
,PlotPoints->{50,50},PlotRange-> {{-3,3},{-0.8,1.2},All},Mesh-> False, ColorFunction->"Rainbow"] ;
Show[U,Graphics[{White,Opacity[0.25],Line[{{-3,0},{3,0}}] }] ]
, {t,0,1}, {θ, 0,80*Pi/180} ,{{n1,2},1.33,2},{{n2,1},1,1.5},{κ2,0,0.1}]

It appears the wave penetrates just a little into 2nd medium and it travels along the surface instead of into the medium.
The reason is for that incident angle θ: 
                             (Subscript[ϵ, 2]) (Subscript[μ, 2])-((((Subscript[n, 1]) Sin[θ]))^2) <0                        (1.2.7)
Hence:    (Subscript[k, 2])=  (Subscript[k, 0])(Sqrt[(Subscript[ϵ, 2]) (Subscript[μ, 2])-((((Subscript[n, 1]) Sin[θ]))^2) ]) = i κ                 (1.2.8a)
where       κ=  (Subscript[k, 0])(Sqrt[((((Subscript[n, 1]) Sin[θ]))^2)-(Subscript[ϵ, 2]) (Subscript[μ, 2]) ])                          (1.2.8b)
Thus:      (Subscript[E, 2])=    (OverHat[y])  t  (e^(i  (Subscript[k, 2]) x ))(e^(i  (β z -ω t ))) =    (OverHat[y])  t  (e^(-κ x ))(e^(i  (β z -ω t )))   (1.2.2) (1.2.9)
The component (e^(-κ x )) is called evanescent, as we have studied above.
If that is the case, what can we expect of reflection? there is no power flowing into medium 2, then where does it go?
We know that        r=(((Subscript[k, 1])-(Subscript[k, 2]))/((Subscript[k, 1])+(Subscript[k, 2])))=(((Subscript[k, 1])-i κ)/((Subscript[k, 1])+i κ))    and              (1.2.10)
                             ((|r|)^2)=((|(((Subscript[k, 1])-i κ)/((Subscript[k, 1])+i κ))|)^2)=((((Subscript[k, 1])^2)+ (κ^2))/(((Subscript[k, 1])^2)+ (κ^2)))=1          (1.2.11)
which means that the reflected power is 100%: this is why we call total internal reflection.
The transmission coefficient 
                              t=((2(Subscript[k, 1]))/((Subscript[k, 1])+(Subscript[k, 2])))=((2(Subscript[k, 1]))/((Subscript[k, 1])+i κ))                    
 is NOT zero, but it doesn’t mean it carries power, because it is not the coefficient t that indicates power, which is only a ratio of the transmitted field amplitude at the boundary vs the incident field. The field decays to zero as x-> ∞, hence there is no power flowing into the medium.

Notice that  ((|r|)^2)=1 , but it is complex, hence it has a phase:

Manipulate[  β=n1 Sin[θ]  ;
k1=n1 Cos[θ] ;k2=(Sqrt[(((n2+I κ2))^2) -(n1^2) ((Sin[θ])^2) ])  ;
w1=Min[8 Cos[θ],3] ; w2=(w1/(Cos[θ]))Re[k2]   ;
r=  ((k1-k2)/(k1+k2)) ; tr=   ((2 k1)/(k1+k2)) ;
Graphics[{Yellow,Disk[{0,0},0.25,{0,Arg[r]}],Red, Arrowheads[0.05],Arrow[{{0,0},{Re[r],Im[r]}}],Blue,Dashed,Circle[{0,0},1]}, PlotRange-> {{-1,1},{-1,1}}, Frame-> True, GridLines-> Automatic]
, {t,0,1}, {θ, 0,90*Pi/180} ,{{n1,2},1.33,2},{{n2,1},1,1.5},{κ2,0,0.1}]

Note the negative phase shift of the reflection coefficient. This is related to Goos-Hanschen shift: the field does penetrate a little into the second medium before returning, hence, it experiences some delay, resulting in the negative phase shift.

1.3 Distributed reflected wave in Bragg reflector (1-D photonic crystal)

A very similar wave is the distributed reflected wave that we learn in previous Lect Set 4, Section 7. Although the form of this wave is  (e^(-γ x)) ( a (e^(i q x))+b  (e^(-i q x))), which is not exactly like TIR evanescent wave (e^(-κ x)), it has the same essential evanescent property in the envelop term (e^(-γ x)): the net wave power does not flow into the evanescent direction, nor is it lossy along the evanescent direction, but it is simply confined and reflected back as illustrated below.

1.3.1 Modules for calculation

Here, we have a slight modification from previous chapter.

Sp[n_List,L_List,λ_,θ_]:=Module[{ dH, dL,kH, kL , k0, nH, nL, d1,d2},
{nH, nL}=n  ; { d1,d2}=L  ;
k0=2 Pi/λ  ; kH=k0 (Sqrt[(nH^2)-(((Sin[θ]))^2)])  ; kL=k0 (Sqrt[(nL^2)-(((Sin[θ]))^2)])  ;
dH=kH*d1  ; dL=kL*d2  ;
(({{cos(dL), ((sin(dL))/kL)}, {-sin(dL) kL, cos(dL)}})) . (({{cos(dH), ((sin(dH))/kH)}, {-sin(dH) kH, cos(dH)}})) ];
Reflect[n_List,L_List,ns_,λ_,θ_,M_]:=Module[{SM,ki,ks} ,
ki=(2 Pi/λ) Cos[θ]  ; ks=  (2 Pi/λ) (Sqrt[(ns^2)-(((Sin[θ]))^2)])   ;
       SM=MatrixPower[Sp[n,L,λ,θ],M] ;
1-((2 (SM[[1,1]] ks+I SM[[2,1]]))/(ks (SM[[1,1]]-I ki SM[[1,2]])+ki SM[[2,2]]+I SM[[2,1]]))  ] ;

Exinc[ λ_,k1_,x_,t_]:=  If[x≤ 0, (E^(I  ( k1 x-2 π  t)))  ,]  ;
Exref[BraggParm_,λ_,k1_, x_,t_]:=Module[ {n,L,ns,M,θ,R } ,
{n,L,ns,M}=BraggParm  ; θ=ArcCos[k1 λ/(2 π)] ;
R= Reflect[n,L,ns,λ, θ,M]  ;
         If[x≤ 0,  R  (E^(-I  ( k1 x + 2 π  t)))   ,   ]
]  ;
ExBragg[BraggParm_,λ_,k1_, x_,t_]:=Module[ {n,L,ns,M, θ,R, Λ
                       ,nH, nL, d1,d2, k0,kH, kL,S, u,m, coef, F} ,
{n,L,ns,M}=BraggParm  ;
{nH, nL}=n  ; { d1,d2}=L  ;Λ=d1+d2   ;
k0=2 Pi/λ  ;    θ=ArcCos[k1 /k0] ;
 kH=k0 (Sqrt[(nH^2)-((Sin[θ])^2)])  ; kL=k0 (Sqrt[(nL^2)-((Sin[θ])^2)])  ;
R= Reflect[n,L,ns,λ,θ,M]   ;
S=(({{cos(kL d2 ), ((sin(kL d2))/kL)}, {-sin(kL d2) kL, cos(kL d2)}})). (({{cos(kH d1), ((sin(kH d1))/kH)}, {-sin(kH d1) kH, cos(kH d1)}}))   ;
If[x> 0 && x< M Λ, 
u=Mod[x,Λ]  ;
m=Floor[x/Λ]  ;
If[u<d1 ,
coef=(({{(1/2), -(I/(2 kH))}, {(1/2), (I/(2 kH))}})) . MatrixPower[S,m] . (({{R+1}, {I k1 (1-R)}})) ;
F={coef[[1,1]] (E^(I  (kH u -2 π t))) , coef[[2,1]]  (E^(-I  (kH u + 2 π t)))  }
,
coef=(({{(1/2), -(I/(2 kL))}, {(1/2), (I/(2 kL))}})) . (({{cos(kH d1), ((sin(kH d1))/kH)}, {-sin(kH d1) kH, cos(kH d1)}})) . MatrixPower[S,m] . (({{R+1}, {I  k1 (1-R)}})) ;
F={coef[[1,1]] (E^(I (kL (u-d1) -2 π t)))  , coef[[2,1]]  (E^(-I  (kL (u-d1)+2 π t)))}
],  ] ;
F
];

1.3.2 Illustration

nb={1.6,1.4}; ns=1.5; λ=1 ;
kc=2 π(Sqrt[(nb^2)-(((Sin[π/2-0.1]))^2)])/λ;
L=π/(2.  kc)   ;
BraggParm={nb, L, ns,20} ;
Manipulate[  θ=π/2-φ  ;k1=2 π* Cos[θ] /λ;
Animate[Plot[{Re[Exinc[λ, k1,x,t]],Re[Exref[BraggParm,λ,k1, x,t]],
Re[Release[ExBragg[BraggParm,λ,k1, x,t]] ]}, {x,-8,10}
,PlotStyle->{Blue,Red,Hue[0.35,1,0.5,1],Purple}
,Filling-> {{1-> {Axis,Hue[0.5,1,1,0.5]} },{2-> {Axis,Hue[0.,1,1,0.5]}}
,{3-> {Axis,Hue[0.35,1,1,0.5]}} (* ,{4-> {Axis,Hue[0.7,1,1,0.5]}} *)}
, PlotRange-> {-1.2,1.2}]
, {t,0,1}, AnimationRunning-> False]
,{{λ,1},0.8,1.2},{{φ,0.325,"axial angle"},0.05,0.95*Pi/2}]

Note how the Bragg reflector relects the wave at the designed axial angle with its distributed reflected wave. (The reflection does not come from a particular interface like a metal mirror with nm-scaled skin depth, but is distributed throughout the Bragg stack (~10 μm in the illustration above), such that the incoming wave exponentially decays as it gives power to the exponentially growing reflecting wave).

nb={1.6,1.4}; ns=1.5; λ=1 ;
kc=2 π(Sqrt[(nb^2)-(((Sin[π/2-0.1]))^2)])/λ;
L=π/(2.  kc)   ;
BraggParm={nb, L, ns,50} ;
Manipulate[  θ=π/2-φ  ;k1=2 π* Cos[θ] /λ;
β=2 π* Sin[θ] /λ ;
Animate[Uplot=DensityPlot[Re[(E^(I β z)) (Exinc[λ, k1,-x,t]+ Exref[BraggParm,λ,k1,-x,t])],{z,-4,4}, {x,0,8}
,PlotPoints->{30,30}
, PlotRange-> {-2,2}, ColorFunction->"Rainbow"] ;
Lplot=DensityPlot[
Re[(E^(I β z)) (ExBragg[BraggParm,λ,k1,-x,t] . {1,1} )],{z,-4,4}, {x,-8,0}
,PlotPoints->{30,30}
, PlotRange-> {-2,2}, ColorFunction->"Rainbow" ] ;
Show[ {Uplot,Lplot},PlotRange-> {{-4,4},{-8,8}}, AspectRatio-> 2 ]
, {t, 0 ,1 }, AnimationRunning-> False]
,{{λ,1},0.8,1.2},{{φ,0.1,"axial angle"},0.05,0.95*Pi/2}]

The left 2 figures shows that when the beam axial angle (angle parallel to the horizontal axis) is small (a glancing incident) such that the vertical wave component is within the Bragg stack reflection band (transmission band gap), the wave has an evanescent envelop (the exponentially decaying portion in the lower half of the figure) and the light is 100% reflected. This is seen from the 100% contrast interference between incident and reflection light.
The right figure shows that at larger angle, the wave is in the transmission band and can propagate into the Bragg stack (the light escapes and is not trapped).
In the left figure situation, if we have BR on all sides, the wave would be trapped and can only travel along the z (horizontal) axis.

2. Waveguiding from evanescent wave confinement

2.1 Concept illustration and discussion

Over all concept:
1- Waveguides are structures that can “conduct” light, which means light can travel in certain reduced dimension(s), while being confined in the other dimension(s) with evanescent field. 
2- The main region where the light travels is called the “core” of the waveguide, it can be a dielectric, or even vacuum or air. The region where the light is evanescent, and effects the confinement of the light is called “cladding,” which can also be a dielectric, or structure that does not allow the light to propagate through and escape.
3- A waveguide allows only certain optical solutions (called optical modes) that can propagate, resulting in certain unique relationship between the frequency of light and it’s propagation and field profile. Optical modes are the key property of a waveguide and determines its application.

2.2 Ray optics concept

In ray optics, the behavior light is modeled with straight lines of rays, which can be used to explain waveguide to a certain extent. However, only wave optics treatment is trule rigorous. Ray optics does offer a the intuitively useful concepts of acceptance angle and numerical aperture.

As illustrated in the left figure, a ray enters into the perpendicular facet of a planar waveguide or fiber can be trapped (or guided) or escape into the cladding depending on the incident angle. If the incident angle is larger than certain value (Subscript[θ, a]), the ray will escape into the core because the condition for total internal reflection is not met. This angle is called the acceptance angle. 
This angle can be calculated as follow. Let (Subscript[θ, c]) be the waveguide internal total reflection angle:
                                                  (Subscript[θ, c])=ArcSin[((Subscript[n, 2])/(Subscript[n, 1]))] 
 Then, the refracted angle at the facet is:
                                                 (Subscript[θ, r])=(π/2)-(Subscript[θ, c])=(π/2)-ArcSin[((Subscript[n, 2])/(Subscript[n, 1]))] 
 The acceptance angle is give by Snell's law:
                                                (Subscript[n, 0])Sin[(Subscript[θ, a])]=(Subscript[n, 1]) Sin[(Subscript[θ, r])] 
                                                (Subscript[n, 0])Sin[(Subscript[θ, a])]=(Subscript[n, 1]) Cos[(Subscript[θ, c])]  
 We can approximate:            (Subscript[n, 0])Sin[(Subscript[θ, a])]=(Subscript[n, 1]) (Sqrt[1-((((Subscript[n, 2])^2) )/((Subscript[n, 1])^2))])=(Sqrt[((Subscript[n, 1])^2)-((Subscript[n, 2])^2) ])≃(Sqrt[(Subscript[ϵ, 1])-(Subscript[ϵ, 2])])
 Since (Subscript[n, 0]) of air is ~ 1, we can also write:   Sin[(Subscript[θ, a])]=≃(Sqrt[(Subscript[ϵ, 1])-(Subscript[ϵ, 2])])

However, ray optics theory of waveguide is only an approximate and is rarely used where accurate and critical engineering design of waveguide and applications are involved.

2.3 Wave optics framework and an example: symmetric slab waveguide

We approach the waveguide problem from the general Helmholtz wave quation 
         (∇^2)E + ((ω^2)/(c^2))ϵ(x,y,z) E = 0           (2.3.1a)                  
         (∇^2)H + ((ω^2)/(c^2))ϵ(x,y,z) H = 0           (2.3.1b)
where ϵ(x,y,z) is a function of space.

In a waveguide, the dielectric should be constant along one dimension, say z. Then, the solution for the electric field, called mode, or transverse optical mode, can be written:
                      E(x,y)  (e^(i(β z - ω t))).    (2.3.2)
The wave equation becomes:
                          (((∂^2)E)/(∂(x^2)))+(((∂^2)E)/(∂(y^2))) -((β^2)-((ω^2)/(c^2)) ϵ(x,y)) E=0   (2.3.3)

The essence is to solve (2.3.3) which usually have infinite number of solutions. But we are only interested in those solutions that satisfy:
                     E(x,y)->0    as  (x^2)+(y^2)-> ∞       (2.3.4)
Of course, we assume that E(x,y)->0  is not due to lossy dielectric  ϵ(x,y), which is assumed to be real. In practice, it may have a small complex component, but this lossy component is not the reason for bounded mode condition in (2.3.4).

In practice,  ϵ(x,y) is often not an analytic or smoothly continuous function of x and y over the entire space, but can be piecewise with different profiles in different domains. How do we approach such a case? The key steps are:
- Solve for solutions at the furthest out domain as (x^2)+(y^2)-> ∞ . This domain is usually called the cladding.
- Select only those solutions that also satisfy (2.3.4):  E(x,y)->0    as  (x^2)+(y^2)-> ∞  
- Join solutions from the outer most domain (cladding) to the inner domain, usually called “core”.

It is most helpful to gain insight by trying to solve a simple problem, and the simplest would be a symmetric dielectric slab waveguide.

3. Symmetric slab waveguide

3.1 Slab waveguide

It is illustrative to try to solve Eq. (2.3.3) for a specific example: symmetric slab wave guide. As the name implies, it is just a slab of dielectric with 2 parallel surfaces. But it doesn’t have to be in air. It can be cladded on both sides with another dielectric. Hence, the profile is shown below.

In fact, it doesn’t look any difference from what we studied from the previous chapter:

Except that here, we are not interested in the solution involving traveling wave incidenting on and exiting the slab in its orthogonal direction. These waves have real (Subscript[k, 1]): (e^(±i (Subscript[k, 1])x )).
We are interested in those solutions that are bounded to the planar layers of the slab, in other words, evanescent both top and bottom (e^(± κ x)) (Re κ)
But mathematically, the problem setup is exactly the same:

But with some difference: the basis solutions we are look for in the cladding layers are (e^(+± κ x)) with κ being real, given by:
                i κ=(Sqrt[(Subscript[ϵ, 1]) ((Subscript[k, 0])^2)-(β^2) ])  or:   κ=(Sqrt[(β^2) -(Subscript[ϵ, 1]) ((Subscript[k, 0])^2)])     (3.1.1)
 Thus:                          (({{(Subscript[F, 1])[x]}, {(Subscript[F, 1])'[x]}}))=(({{(Subscript[ψ, 1])[x], (Subscript[ψ, 2])[x]}, {(Subscript[ψ, 1])'[x], (Subscript[ψ, 2])'[x]}})) . (({{(Subscript[A, left])}, {(Subscript[B, left])}}))        4(3.2.14)  (3.1.2)
                                   (({{(Subscript[F, 1])[x]}, {(Subscript[F, 1])'[x]}}))=(({{(e^(+κ (x+a))), (e^(+- κ (x+a)))}, {((κe)^(+κ (x+a))), -((κe)^(+-κ (x+a)))}})) . (({{(Subscript[A, left])}, {(Subscript[B, left])}}))          (3.1.3)
 Also, we note that we can’t have (Subscript[B, left])≠0 because the solution would go to ∞ as x-> -∞, hence:
                                   (({{(Subscript[F, 1])[x]}, {(Subscript[F, 1])'[x]}}))=(({{(e^(+κ (x+a))), (e^(+- κ (x+a)))}, {((κe)^(+κ (x+a))), -((κe)^(+-κ (x+a)))}})) . (({{(Subscript[A, left])}, {0}}))          (3.1.4) 
  We choose basis function (e^(+±κ (x+a))) instead of  (e^(+± κ x)) is just for convenience when we set x=-a:
                           (({{(Subscript[F, 1])[-a]}, {(Subscript[F, 1])'[-a]}}))=(({{1, 1}, {κ, -κ}})) . (({{(Subscript[A, left])}, {0}}))            (3.1.5)

For the slab, we define:
                                          (Subscript[k, 2])=(Sqrt[(Subscript[ϵ, 2]) ((Subscript[k, 0])^2) -(β^2) ])                     (3.1.6)
One question is, how do we know (Subscript[k, 2]) is real? We don’t care it is real or complex. We simply define it to be that way and will find out later whether it is real or complex after obtaining the solution.
The transfer matrix in the slab is:
                S=(({{(Subscript[ψ, 1])[a], (Subscript[ψ, 2])[a]}, {(Subscript[ψ, 1])'[a], (Subscript[ψ, 2])'[a]}})) (((({{(Subscript[ψ, 1])[-a], (Subscript[ψ, 2])[-a]}, {(Subscript[ψ, 1])'[-a], (Subscript[ψ, 2])'[-a]}})))^(-1))             4(3.2.19) (3.1.7)
                S=(({{(e^(i (Subscript[k, 2])a )), (e^(-i (Subscript[k, 2])a ))}, {i (Subscript[k, 2])(e^(i (Subscript[k, 2])a )), -i (Subscript[k, 2])(e^(-i (Subscript[k, 2])a ))}})) (((({{(e^(-i (Subscript[k, 2])a )), (e^(i (Subscript[k, 2])a ))}, {i (Subscript[k, 2])(e^(-i (Subscript[k, 2])a )), -i (Subscript[k, 2])(e^(i (Subscript[k, 2])a ))}})))^(-1))    (3.1.8)

FullSimplify[(({{(E^(I (Subscript[k, 2])a )), (E^(-I (Subscript[k, 2])a ))}, {I (Subscript[k, 2])(E^(I (Subscript[k, 2])a )), -I (Subscript[k, 2])(E^(-I (Subscript[k, 2])a ))}})) . Inverse[(({{(E^(-I (Subscript[k, 2])a )), (E^(I (Subscript[k, 2])a ))}, {I (Subscript[k, 2])(E^(-I (Subscript[k, 2])a )), -I (Subscript[k, 2])(E^(I (Subscript[k, 2])a ))}}))]]

We have seen this before:
                      S= (({{cos(2 a (Subscript[k, 2])), ((sin(2 a (Subscript[k, 2])))/(Subscript[k, 2]))}, {-sin(2 a (Subscript[k, 2])) (Subscript[k, 2]), cos(2 a (Subscript[k, 2]))}}))            (3.1.9)
 which is the same as:
        S=(({{(Subscript[ψ, 1])[l], (Subscript[ψ, 2])[l]}, {(Subscript[ψ, 1])'[l], (Subscript[ψ, 2])'[l]}})) (((({{(Subscript[ψ, 1])[0], (Subscript[ψ, 2])[0]}, {(Subscript[ψ, 1])'[0], (Subscript[ψ, 2])'[0]}})))^(-1))=(({{cos(l (Subscript[k, 2])), ((sin(l (Subscript[k, 2])))/(Subscript[k, 2]))}, {-sin(l (Subscript[k, 2])) (Subscript[k, 2]), cos(l (Subscript[k, 2]))}}))   4(4.1.3) 
except that here:             l =2 a

Finally:       
                 (({{(Subscript[F, 3])[x]}, {(Subscript[F, 3])'[x]}}))=(({{(e^(+κ (x-a))), (e^(+- κ (x-a)))}, {((κe)^(+κ (x-a))), -((κe)^(+-κ (x-a)))}})) . (({{(Subscript[A, right])}, {(Subscript[B, right])}}))          (3.1.10a)
Notice that we have to choose (Subscript[A, right])=0, otherwise we have an infinite solution at x->∞.
                 (({{(Subscript[F, 3])[x]}, {(Subscript[F, 3])'[x]}}))=(({{(e^(+κ (x-a))), (e^(+- κ (x-a)))}, {((κe)^(+κ (x-a))), -((κe)^(+-κ (x-a)))}})) . (({{0}, {(Subscript[B, right])}}))          (3.1.10b)       
and:          (({{(Subscript[F, 3])[a]}, {(Subscript[F, 3])'[a]}}))=(({{1, 1}, {κ, -κ}})) . (({{0}, {(Subscript[B, right])}}))      (3.1.11)

We can put all together:
    (({{(Subscript[F, 3])[a]}, {(Subscript[F, 3])'[a]}}))=(({{1, 1}, {κ, -κ}})) . (({{0}, {(Subscript[B, right])}}))=(({{cos(2 a (Subscript[k, 2])), ((sin(2 a (Subscript[k, 2])))/(Subscript[k, 2]))}, {-sin(2 a (Subscript[k, 2])) (Subscript[k, 2]), cos(2 a (Subscript[k, 2]))}})) . (({{(Subscript[F, 1])[-a]}, {(Subscript[F, 1])'[-a]}}))    (3.1.12a)     
        (({{1, 1}, {κ, -κ}})) . (({{0}, {(Subscript[B, right])}}))=(({{cos(2 a (Subscript[k, 2])), ((sin(2 a (Subscript[k, 2])))/(Subscript[k, 2]))}, {-sin(2 a (Subscript[k, 2])) (Subscript[k, 2]), cos(2 a (Subscript[k, 2]))}})) .(({{1, 1}, {κ, -κ}})) . (({{(Subscript[A, left])}, {0}}))       (3.1.12b)
Or:
     (({{0}, {(Subscript[B, right])}}))=(((({{1, 1}, {κ, -κ}})))^(-1)) . (({{cos(2 a (Subscript[k, 2])), ((sin(2 a (Subscript[k, 2])))/(Subscript[k, 2]))}, {-sin(2 a (Subscript[k, 2])) (Subscript[k, 2]), cos(2 a (Subscript[k, 2]))}})) .(({{1, 1}, {κ, -κ}})) . (({{(Subscript[A, left])}, {0}}))     (3.1.12c)

(({{0}, {(Subscript[B, right])}}))==FullSimplify[Inverse[ (({{1, 1}, {κ, -κ}}))].(({{cos(2 a (Subscript[k, 2])), ((sin(2 a (Subscript[k, 2])))/(Subscript[k, 2]))}, {-sin(2 a (Subscript[k, 2])) (Subscript[k, 2]), cos(2 a (Subscript[k, 2]))}})). (({{1, 1}, {κ, -κ}}))]. (({{(Subscript[A, left])}, {0}}))

This equation is the key result we are looking for:
                   (({{0}, {(Subscript[B, right])}}))==(({{(Subscript[A, left]) (cos(2 a (Subscript[k, 2]))+((sin(2 a (Subscript[k, 2])) ((κ^2)-(Subscript[k, 2]^2)))/(2 κ (Subscript[k, 2]))))}, {((sin(2 a (Subscript[k, 2])) (Subscript[A, left]) ((κ^2)+(Subscript[k, 2]^2)))/(2 κ (Subscript[k, 2])))}}))   (3.1.13)

3.2 Characteristic equation

Expand Eq. (3.1.13) above into 2 equations:
   0=(Subscript[A, left]) (cos(2 a (Subscript[k, 2]))+((sin(2 a (Subscript[k, 2])) ((κ^2)-(Subscript[k, 2]^2)))/(2 κ (Subscript[k, 2]))))             (3.2.1a)
   (Subscript[B, right]) =((sin(2 a (Subscript[k, 2])) ((κ^2)+(Subscript[k, 2]^2)))/(2 κ (Subscript[k, 2])))(Subscript[A, left])                                (3.2.1b)
Eq. (3.2.1b) is merely a proportional relationship between 2 constant (Subscript[A, left]) , (Subscript[B, right]). But the only way (3.2.1a) has a solution that is non-trivial is:

cos(2 a (Subscript[k, 2]))+((sin(2 a (Subscript[k, 2])) ((κ^2)-(Subscript[k, 2]^2)))/(2 κ (Subscript[k, 2])))==0            (3.2.2a) 
Or:                         2 (Subscript[k, 2]) κ cos(2 a (Subscript[k, 2]))+((κ^2)-(Subscript[k, 2]^2)) sin(2 a (Subscript[k, 2]))== 0     (3.2.2b) 
This is called the characteristic equation. What quantity is undetermined in this equation?
- First, we note that κ and (Subscript[k, 2]) are defined:           κ=(Sqrt[(β^2) -(Subscript[ϵ, 1]) ((Subscript[k, 0])^2)]) ;  (Subscript[k, 2])=(Sqrt[(Subscript[ϵ, 2]) ((Subscript[k, 0])^2) -(β^2) ])
- Next, we note that parameters (Subscript[ϵ, 1]) , (Subscript[ϵ, 2]), and a are given. ω is of course also given.

The only thing we can’t choose is β. This result says that we can’t just choose any arbitrary value of β, but must find one that satisfies the characteristic equation (3.2.2b) above.

In fact, this is the most important concept in waveguide. All waveguides have some characteristic equation. Even if some cannot be expressed in analytic form like here, it is an equation (numerical) that dictates what the allowed propagation constant should be, and determines the optical modes that can be transmitted in the waveguide.

There are several way to express the same characteristic equation above that can give us some insight of the problem:

TrigExpand[cos(2 a (Subscript[k, 2]))+((sin(2 a (Subscript[k, 2])) ((κ^2)-(Subscript[k, 2]^2)))/(2 (Subscript[k, 2]) κ))]

Factor[%]

Hence, the equation can also be written:
                   (cos(a (Subscript[k, 2])) κ-sin(a (Subscript[k, 2])) (Subscript[k, 2])) (cos(a (Subscript[k, 2])) (Subscript[k, 2])+sin(a (Subscript[k, 2])) κ)==0    (3.2.3a) 
which are really 2 equations:
       (even modes)           cos(a (Subscript[k, 2])) κ-sin(a (Subscript[k, 2])) (Subscript[k, 2]) ==0   (3.2.3b)  
       (odd modes)            cos(a (Subscript[k, 2])) (Subscript[k, 2])+sin(a (Subscript[k, 2])) κ==0    (3.2.3c)

What is the meaning of this? It means that there are two classes of modes and solutions, one for each equation. But why? Because of symmetry. Like in any quantum problem, if we have symmetry, we have a subset of eigensolutions that belong to an eigenspace of that geometry. Similarly here, we have a symmetric waveguide, which has the reflection symmetry that should have even solution: F[x]=F[-x] or odd solution F[x]=-F[-x]: We can guess that each "sub-characteristic" equation is for each class of modes and solutions.

If fact, (3.2.3b) gives us the even solutions of the form Cos[(Subscript[k, 2]) x]  in the core.  Eq. (3.2.3c) gives us the odd solutions with the form  Sin[(Subscript[k, 2]) x].
For these equations, if (Subscript[k, 2]) is imaginary, it has no solution. Hence, we now know that it has to be real, which makes sense for a traveling wave.

3.3 Symmetric slab waveguide solution

Once we obtain β from (2.3.18), the modes are given:
          Left (or bottom):    (e^(κ (x+a)))                                         (3.3.1a)
          Core:                   c (e^(i (Subscript[k, 2])x ))+d  (e^(-i (Subscript[k, 2])x ))                        (3.3.1b)
         Right (or top):       (e^(-κ (x-a)))                                         (3.3.1c)
Coefficients c and d are given from:
      (({{(e^(-i (Subscript[k, 2])a )), (e^(i (Subscript[k, 2])a ))}, {i (Subscript[k, 2])(e^(-i (Subscript[k, 2])a )), -i (Subscript[k, 2])(e^(i (Subscript[k, 2])a ))}})) . (({{c}, {d}}))=(({{1, 1}, {κ, -κ}})) . (({{1}, {0}}))         (3.3.1d)

Let’s do some illustration.

3.4 Solve the characteristic equation for β

Execute the code package at the top first for computation in this section.

This is the index profile. If core index is < cladding index:

indexprof={3.5,3.4}; a=2.5 ; λ=1.5 ;
Plot[IndexprofSSW[x,indexprof,a] ,{x,-5,5}, PlotRange-> All, Filling-> {1-> {Axis,Hue[0.3,1,1,0.5]}}, PlotStyle-> {Thick,Red}] 
Plot[Release[ResolvTE[indexprof,a,λ,x]],{x,3.5,3.7}, PlotStyle-> {Red,Blue,{Dashed,Black}}]

We see that the characteristic function has no solution. The left hand side is complex.

This is the index profile. Now, let core index > cladding index:

indexprof={3.4,3.5}; a=2.5 ; λ=1.5 ;
idpplot=Plot[IndexprofSSW[x,indexprof,a] ,{x,-5,5}, PlotRange-> All, Filling-> {1-> {Axis,Hue[0.3,1,1,0.5]}}, PlotStyle-> {Thick,Red}] ;
resfplot=Plot[Release[ResolvTE0[indexprof,a,λ,x]],{x,3.4,3.499}, PlotStyle-> {Red,Blue,{Dashed,Black}}] ;
Manipulate[
{Show[idpplot, Graphics[{Thick,Black,Line[{{-5,neff},{5,neff}}]}] ]
,Show[resfplot, Graphics[{Thick,Hue[0.3,1,0.5],Line[{{neff,-10},{neff,10}}]}] ]}
,{neff,indexprof[[1]], indexprof[[2]]} ]

indexprof={3.4,3.5}; a=2.5 ; λ=1.5 ;
idpplot=Plot[IndexprofSSW[x,indexprof,a] ,{x,-5,5}, PlotRange-> All, Filling-> {1-> {Axis,Hue[0.3,1,1,0.5]}}, PlotStyle-> {Thick,Red}, ImageSize->{150,150}] ;
resfplotV=ParametricPlot[{ResolvTE0[indexprof,a,λ,x][[3]],x} ,{x,3.4,3.499}, PlotStyle-> {Thick,Purple}
,AspectRatio-> 1, ImageSize->{150,150}] ;
Manipulate[
{Show[idpplot, Graphics[{Thick,Black,Line[{{-5,neff},{5,neff}}]}] , ImageSize->{150,150}]
,Show[resfplotV, Graphics[{Thick,Hue[0.3,1,0.5],Line[{{-10,neff},{10,neff}}]}] ]}
,{neff,indexprof[[1]], indexprof[[2]]} ]

We define              (Subscript[n, eff])=(β/(Subscript[k, 0])) ;                       (3.4.1)
which is also called effective index (of refraction) of a mode

Let’s collect all the mode effective index:

indexprof={3.4,3.5}; a=2.5 ; λ=1.5 ;
segment={3.499,3.495,3.48,3.47,3.45,3.42,3.405}; neff={} ;
For[ i=1, i≤ 6, i++,
	AppendTo[neff, 
	 x /. FindRoot[ResolvTE[indexprof,a,λ,x]==0,{x,segment[[i+1]],segment[[i]]}]
                            ] ];
neff

The effective indices of the 6 modes are:

neff

resfplotVsm=ParametricPlot[{0.25*ResolvTE0[indexprof,a,λ,x][[3]],x} ,{x,3.4,3.499}, PlotStyle-> {Thick,Purple}
] ;
Show[idpplot, resfplotVsm,Graphics[
Table[{Thick,Hue[(i-1)/(0.1+Length[neff]),1,1],Line[{{-5,neff[[i]]},{5,neff[[i]] }}]}
,{i,1,Length[neff]}] ] ]

3.5 Plotting E field (TE polarization)

3.5.1 (Subscript[E, y]) profile

Let’s calculate and plot (Subscript[E, y]) spatial profile, which is F[x] function in our discussion above:

indexprof={3.4,3.5}; a=2.5 ; λ=1.5 ;
allp={};
For[i=1,i≤ Length[neff],i++,
	         AppendTo[allp,  Plot[FieldTE[indexprof,a,λ,neff[[i]],x][[1]],{x,-7,7}
	, PlotRange->{-1.,1}, PlotStyle->{Hue[(i-1)/(0.1+Length[neff]),1,0.4]}
,Filling-> {1-> {Axis, Hue[(i-1)/(0.1+Length[neff]),1,1,0.35]}} ]
]
             ] ;
allp

{,,,,,}

3.5.2 Traveling electric field

Now, we plot the whole (Subscript[E, y]) as a function of time. Note that it is essentially a combination of a standing wave in the x direction (because of waveguide) and a traveling wave in the z direction:

Manipulate[
k0=2 π/λ ;
Animate[Plot3D[FieldTE[indexprof,a,λ,neff[[i]],x][[1]]Cos[(z neff[[i]]*k0-2 π t)]
,{x,-5,5},{z,0,2}, Mesh-> False, PlotPoints-> {25,40}, PlotRange-> {-0.8,0.8}, BoxRatios->{1,1,0.5}],{t,0,1,0.025}, DisplayAllSteps->True, AnimationRunning->False]
,{i,1,6,1}]

k0=2 π/λ ;allp0={}; i=4 ;
For[j=0,j≤ 40,j++,t=0.025 j ;
AppendTo[allp0,Plot3D[FieldTE[indexprof,a,λ,neff[[i]],x][[1]]Cos[(z neff[[i]]*k0-2 π t)]
,{x,-5,5},{z,0,1.2}, Mesh-> False, PlotPoints-> {25,40}, PlotRange-> {-0.8,0.8}, BoxRatios->{1,1,0.5}]
]
];

Animate[Show[allp0[[i]] ],{i,1,Length[allp0],1}, DisplayAllSteps->True, AnimationRunning->False]

Manipulate[
k0=2 π/λ ;
Animate[DensityPlot[FieldTE[indexprof,a,λ,neff[[i]],x][[1]]Cos[(z neff[[i]]*k0-2 π t)]
,{x,-5,5},{z,0,1.5}, Mesh-> False, PlotPoints-> {25,30}, PlotRange-> {-0.8,0.8}, ColorFunction->"Rainbow"],{t,0,1,0.02}, DisplayAllSteps->True, AnimationRunning->False]
,{i,1,6,1}]

k0=2 π/λ ;allp1={}; i=5 ;
For[j=0,j≤ 40,j++,t=0.025 j ;
AppendTo[allp1,DensityPlot[FieldTE[indexprof,a,λ,neff[[i]],x][[1]]Cos[(z neff[[i]]*k0-2 π t)]
,{x,-4,4},{z,0,1.5}, Mesh-> False, PlotPoints-> {25,30}, PlotRange-> {-0.8,0.8}, ColorFunction->"Rainbow"]
]
];

Animate[Show[allp1[[i]],AspectRatio->1.5/8 ],{i,1,Length[allp1],1}, DisplayAllSteps->True, AnimationRunning->False]

3.6 Vector graphics

Remember that EM field are vector field. Not scalar field. If it has only one component like Ey above, then we can look at the amplitude of that single component just as a scalar plot. But if it has more than one component, to have a sense of what really happens, we must look at vector field. This is the case for the magnetic induction field.

This is a plot of the amplitude of Ey, along with Hz and Hx for mode 3. You can change the mode parameter to plot other modes.

mode=3 ; nep=neff[[mode]] ;
Plot[ {FieldTE[indexprof,a,λ, nep,x][[1]]
		       ,Im[FieldTE[indexprof,a,λ, nep,x] [[2]] ]
				,-FieldTE[indexprof,a,λ, nep,x][[3]]},
		        {x,-5,5}
	, PlotRange->All
	, PlotStyle->{{Hue[0,1,0.8]}
			                      , {Hue[0.28,1,0.5]}
			                       ,{Hue[0.6,1,0.8]}  }]

3.6.1 Ey in the x-z plane

This is an illustration of E field (not to scale x, z) in the x-z plane (y=0)

Below is the code:

indexprof={3.4,3.5}; a=2.5 ; λ=1.5 ;
segment={3.499,3.495,3.48,3.47,3.45,3.42,3.405}; neff={} ;
For[ i=1, i≤ 6, i++,
	AppendTo[neff, 
	 x /. FindRoot[ResolvTE[indexprof,a,λ,x]==0,{x,segment[[i+1]],segment[[i]]}]
                            ] ]

indexprof={3.4,3.5}; a=2.5 ; λ=1.5 ;k0=2 π/λ ;
 mode=1 ; 
pz=40; px=10 ; Δz=1./pz  ;Δx=8/(2*px ) ;
eyvec={} ;
For[j=0,j< 40,j++,t=j/40. ;
AppendTo[ eyvec,
Graphics3D[Flatten[ 
Table[{Hue[1.3*FieldTE[indexprof,a,λ,neff[[mode]],Δx m][[1]],0.7,1,0.8]
,Arrowheads[0.02],Thickness[0.0025],Arrow[{{ Δz l,0,Δx m},{Δz l
,10 FieldTE[indexprof,a,λ,neff[[mode]],Δx m][[1]]Cos[(Δz l neff[[mode]]*k0-2 π t)]
,Δx m}} ]  },{m,-px,px},{l,0,pz}]
,1]
, BoxRatios-> {1.5,0.7,1} (* , PlotRange-> {{0,2},{-1,1},{-1,1}}  *)]
]
]

Manipulate[eyvec[[j]],{j,1,Length[eyvec],1}]

3.6.2 Magnetic field

Remember that EM field are vector field. Not scalar field. If it has only one component like Ey above, then we can look at the amplitude of that single component just as a scalar plot. But if it has more than one component, to have a sense of what really happens, we must look at vector field. This is the case for the magnetic induction field.

Linked to Fig. 3.1.3.1: Magnetic field lines of mode 0

The vector field varies smoothly at the center and "turns" at the boundary

Linked to Fig. 3.1.3.2

Linked to Fig. 3.1.3.3

Where we have interference with higher kx, we see "twirl" or vortex, or  whirpool: this is in higher order mode

Linked to Fig. 3.1.3.4

Linked to Fig. 3.1.3.5

3.6.3 Field lines

Field line is a concept that "imagines" line that a vector field tries to follow. If we have a line {x[s],y[s],z[s]}. Then its tangential vector is:
{(∂_s) x[s],(∂_s) y[s],(∂_s) z[s]}=α {(Subscript[F, x]),(Subscript[F, y]),(Subscript[F, z])}   where α is just a proportional constant, which can be 1. Thus, the system of D. E. is:
(∂_s) (({{x[s]}, {y[s]}, {z[s]}}))=(({{(Subscript[F, x])[x[s],y[s],z[s]]}, {(Subscript[F, y])[x[s],y[s],z[s]]}, {(Subscript[F, z])[x[s],y[s],z[s]]}})) . It is quite simple to integrate this eq. numerically. In simple case, we can even do it analytically. For example in this case: 
(∂_s) (({{x[s]}, {y[s]}, {z[s]}}))=Re[(({{-i β Cos[k x]}, {0}, {-k Sin[k x]}}))(e^(i β z)) ] =(({{β Cos[k x]Sin[β z]}, {0}, {-k Sin[k x]Cos[β z]}}))

TrigReduce[(({{β Cos[k x[s]]Sin[β z[s]]}, {-k Sin[k x[s]]Cos[β z[s]]}}))]

(({{k x(s)+β z(s)}, {k x(s)-β z(s)}}))==(({{k, β}, {k, -β}})).(({{x}, {z}}))

FullSimplify[(({{k, β}, {k, -β}})).(({{(1/2) (β Sin[u[s]]-β Sin[v[s]] )}, {(1/2) (-k Sin[v[s]]-k Sin[u[s]] )}}))]

(∂_s) (({{u[s]}, {v[s]}}))==(({{- sin(v(s))}, { sin(u(s))}}))

u' Sin[u]+v' Sin[v]==0
(∂_s) (Cos[u]+Cos[v])==0
Cos[u]+Cos[v]= Constant

k=0.25; β=1.;fl={};
For[i=1,i<21,i++,AppendTo[fl,
ContourPlot[Cos[k x+β y]+Cos[k x-β y]==(i-10)*0.2,{y,-10,10},
	                        {x,-10,10}, PlotPoints->50,ContourStyle->{Hue[i*0.045,1,0.9], Thickness[0.005]}, DisplayFunction->Identity]]] ;
Show[fl, DisplayFunction->$DisplayFunction]

3.7 Modal dispersion

How does the β depend on ω?

nindex={1.4,1.5}; a=1.25 ;
c=2.99792458 ; 
plot1=Plot3D[ResolvTE[nindex,a,c/f,x],
		          {x,nindex[[2]]-0.005,nindex[[1]]+0.005},
                 {f,2,4.}, PlotPoints->{70,30}, DisplayFunction->Identity] ;
plot2=Plot3D[0,
		           {x,nindex[[2]]-0.005,nindex[[1]]+0.005},
                 {f,2,4.}, PlotPoints->2,DisplayFunction->Identity] ;
Show[plot1,plot2,DisplayFunction->$DisplayFunction,
ImageSize->{500,500}, BoxRatios-> {1,1,1},ViewPoint->{2,-5,3}]

Linked to Fig. 4.4.1

Plots of β[ω] or (Subscript[n, eff])[ω] or β[f] or (Subscript[n, eff])[f]  or  β[λ] or (Subscript[n, eff])[λ]  as a function of frequency are called dispersion curves.

nindex={1.4,1.5}; a=1.25 ;
c=2.99792458 ; 
plot1=Plot3D[ResolvTE[nindex,a,λ,x],
		          {x,nindex[[2]]-0.005,nindex[[1]]+0.005},
                 {λ,0.75,1.5}, PlotPoints->{70,30}, DisplayFunction->Identity] ;
plot2=Plot3D[0,
		           {x,nindex[[2]],nindex[[1]]},
                {λ,0.75,1.5}, PlotPoints->2,DisplayFunction->Identity] ;
Show[plot1,plot2,DisplayFunction->$DisplayFunction,
ImageSize->{500,500}, BoxRatios-> {1,1,1},ViewPoint->{2,-5,3}]

We would expect that β[ω]=((2 π (Subscript[n, eff])[ω])/λ) =((ω (Subscript[n, eff])[ω])/c)≈((ω (Sqrt[⟨ϵ⟩]))/c) where ⟨ϵ⟩ is some sort of averaged dielectric constant of the waveguide, because (Subscript[n, eff])[ω] is approximately the modal average of the entire dielectric constant profile:
                (Subscript[n, eff;m])[ω] ≈ (Sqrt[(((∫_(-∞)^∞)ϵ[x] (((|(Subscript[E, m])[x]|))^2)dx)/((∫_(-∞)^∞) (((|(Subscript[E, m])[x]|))^2)dx))])               (3.7.1)
However, as the frequency change, the mode profile changes, and (Subscript[n, eff])[ω] is not constant for each mode.
We define dispersion terms:

InverseGVDterm=(∂_ω) ((ω (Subscript[n, eff])[ω])/c) ;
			
dispersion=(∂_(ω,ω))((ω (Subscript[n, eff])[ω])/c)

InverseGVDterm

We see that the FIRST order derivative of (Subscript[n, eff])[ω] already contributes to the dispersion term. This term ((2 (Subscript[n, eff]^′)(ω))/c)+((ω (Subscript[n, eff]^″)(ω))/c) is called modal dispersion, if we exclude the contribution of material dispersion from (Subscript[n, eff]^′)(ω) and (Subscript[n, eff]^″)(ω). Thus, by virtue of the change of the mode profile vs. ω, we have a change of (Subscript[n, eff])[ω] that gives modal dispersion.