Optical Fiber Intro - Part 1
Waveguide concepts

Overview and study guide

Utility package

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 κ1-ϵ1 ϵ2 Cos[4 a k2] κ1+2 a ((k2^2) ((ϵ1)^2)+((ϵ2)^2) ((κ1)^2)))/(2 (k2^2) ((ϵ1)^2)))
                     +          (((k2 ϵ1-ϵ2 κ1) (k2 ϵ1+ϵ2 κ1) Sin[4 a k2])/(4 (k2^3) ((ϵ1)^2))))
					+((((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 }
   , ((Sin[2 a k2]  ((k2^2)((ϵ1)^2)+((κ1)^2)((ϵ2)^2)))/(2 k2 κ1 ϵ1 ϵ2))  *   {(E^(-κ1 (x-a) ))
					                    , -I κ1 (E^(-κ1 (x-a) ))/(ϵ1*k0 ) , neff  (E^(-κ1 (x-a) ))/ϵ1   }] ]
					
]

ParmTM[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=((((neff/(ϵ2))))^2)((( ϵ1 ϵ2 κ1-ϵ1 ϵ2 cos(4 a k2) κ1+2 a ((k2^2) ((ϵ1)^2)+((ϵ2)^2) ((κ1)^2)))/(2 (k2^2) ((ϵ1)^2)))
+    (((k2 ϵ1-ϵ2 κ1) (k2 ϵ1+ϵ2 κ1) sin(4 a k2))/(4 (k2^3) ((ϵ1)^2))))
		+((((neff/(ϵ1))))^2)(1/( κ1)) ; 
{a,ϵ1,ϵ2, neff, κ1, k2,k0,Norm}  ] ;
FieldTM2[parm_List,x_]:= Module[{a,ϵ1,ϵ2, neff, κ1, k2,k0,Norm} ,
{a,ϵ1,ϵ2, neff, κ1, k2,k0,Norm}= parm  ; 
(1/(Sqrt[Norm])) 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 }
			, ((sin(2 a k2)  ((k2^2)((ϵ1)^2)+((κ1)^2)((ϵ2)^2)))/(2 k2 κ1 ϵ1 ϵ2))  
	{(E^(-κ1 (x-a) )) , -I κ1 (E^(-κ1 (x-a) ))/(ϵ1*k0 ) , neff  (E^(-κ1 (x-a) ))/ϵ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  *)
			 ] ]
					
		      ] ;

ResolvTE[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) κ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  *)
			 ] ]
					
		      ] ;

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 κ1-ϵ1 ϵ2 Cos[4 a k2] κ1+2 a ((k2^2) ((ϵ1)^2)+((ϵ2)^2) ((κ1)^2)))/(2 (k2^2) ((ϵ1)^2)))
                     +          (((k2 ϵ1-ϵ2 κ1) (k2 ϵ1+ϵ2 κ1) Sin[4 a k2])/(4 (k2^3) ((ϵ1)^2))))
					+((((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 }
   , ((Sin[2 a k2]  ((k2^2)((ϵ1)^2)+((κ1)^2)((ϵ2)^2)))/(2 k2 κ1 ϵ1 ϵ2))  *   {(E^(-κ1 (x-a) ))
					                    , -I κ1 (E^(-κ1 (x-a) ))/(ϵ1*k0 ) , neff  (E^(-κ1 (x-a) ))/ϵ1   }] ]
					
]

ParmTM[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=((((neff/(ϵ2))))^2)((( ϵ1 ϵ2 κ1-ϵ1 ϵ2 cos(4 a k2) κ1+2 a ((k2^2) ((ϵ1)^2)+((ϵ2)^2) ((κ1)^2)))/(2 (k2^2) ((ϵ1)^2)))
+    (((k2 ϵ1-ϵ2 κ1) (k2 ϵ1+ϵ2 κ1) sin(4 a k2))/(4 (k2^3) ((ϵ1)^2))))
		+((((neff/(ϵ1))))^2)(1/( κ1)) ; 
{a,ϵ1,ϵ2, neff, κ1, k2,k0,Norm}  ] ;
FieldTM2[parm_List,x_]:= Module[{a,ϵ1,ϵ2, neff, κ1, k2,k0,Norm} ,
{a,ϵ1,ϵ2, neff, κ1, k2,k0,Norm}= parm  ; 
(1/(Sqrt[Norm])) 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 }
			, ((sin(2 a k2)  ((k2^2)((ϵ1)^2)+((κ1)^2)((ϵ2)^2)))/(2 k2 κ1 ϵ1 ϵ2))  
	{(E^(-κ1 (x-a) )) , -I κ1 (E^(-κ1 (x-a) ))/(ϵ1*k0 ) , neff  (E^(-κ1 (x-a) ))/ϵ1
					                }
]]
  ] ;

1. Introduction waveguide concept

1.1 Discussion

If we let a light beam freely propagates in free space, what happens to the beam?

Source code: Demo Free space diffraction

Manipulate[
scalelist={"μm","10μm","100μm","mm","cm","m","10m","100m","km"};
wmin =0.5 ; λmax=3. ;
If[RUN==1, 
Manipulate[
pos=Position[scalelist,zR ] ; 
ipos=pos[[1,1]] ;zmax=(10^ipos) ;
z0=π (w0^2)/λ ; z0lim=π (wmin^2)/λmax ;
bwmax=0.75*(10^(-1))wmin (Sqrt[1+((10^(2(ipos-1)+2))/(z0lim^2))]) ; bwexp=Floor[Log[10,bwmax] ] ;
bwR= bwmax*(10^(-bwexp));
 Plot[{-w0 (Sqrt[1+(((z^2)*(10^(2(ipos-1))))/(z0^2))])*(10^(-bwexp)),w0 (Sqrt[1+(((z^2)(10^(2(ipos-1))))/(z0^2))])*(10^(-bwexp))},{z,0,10},AspectRatio->0.25,PlotRange->{-bwR,bwR},Filling->{1->{{2}, Hue[0.42,1,1,0.7]}}
,Frame-> True,FrameLabel->{"1 unit= "<>zR,"Beam size\n 1 unit="<>scalelist[[bwexp+1]]}
,LabelStyle->{18,FontFamily-> "Calibri"}, ImageSize->{600,200}]
,{{ w0,10,"Beamwaist (μm)"},0.5,20},{{ λ,1.5,"wavelength (μm)"},1,3}
,{{zR,"10μm","Range scale"},{"μm","10μm","100μm","mm","cm","m","10m","100m","km"}, ControlType->SetterBar}]
,  ]
,{RUN,{0,1}}]

Demo Free space diffraction

The beam grows big. It can grow small like when we focus, but becomes big again after the minimum spot.
This is known as diffraction.

ParametricPlot3D[{(Sqrt[1+0.1 (z^2)])Cos[ϕ],(Sqrt[1+0.1(z^2)]) Sin[ϕ], z},{z,-3,8},{ϕ,-Pi,Pi}]

If we like to propagate a beam while keeping it the same shape (cross section) everywhere, how do we do it? Example:

A pair of parallel planar mirrors!
But we don't have to use mirrors, there are other ways, the main concept is that we can engineer a medium in such a way that we can use to "guide" the beam along the desired direction of propagation, this is the concept of waveguide. Optical waveguide is simply a waveguide for light. This is to distinguish it from waveguides for other EM waves such as microwave, RF...

1.2 Intro to optical waveguides

There are several ways to achieve optical waveguide. It is relevant, from the practical point of view, to be aware of technologies that are well established in applications vs. those that are theoretically feasible but might not be practical. Mirror or reflecting-wall optical waveguide, for example can be achieved, but are not as practical as dielectric waveguides. Holey fibers or photonic bandgap (PBG) waveguides are likewise promising for special cases but have not been as widespread as dielectric waveguides.

Here we'll concentrate on dielectric waveguides because this is still the principal technology in fiber optic communications.

1.2.1 Concept of total internal reflection

Snell's law (Subscript[n, 1]) Sin[(Subscript[θ, 1])]=(Subscript[n, 2]) Sin[(Subscript[θ, 2])]
or:                     Sin[(Subscript[θ, 2])]=((Subscript[n, 1])/(Subscript[n, 2]))Sin[(Subscript[θ, 1])] 
If  (Subscript[n, 1]) >(Subscript[n, 2]) , there will be values of (Subscript[θ, 1]) such that Sin[(Subscript[θ, 2])]>1. Whent this happens, total internal reflection occurs, and the critical angle is (Subscript[θ, c])=ArcSin[((Subscript[n, 2])/(Subscript[n, 1]))] .

1.2.2 Waveguiding effect from total internal reflection

Summary:
1- Waveguides are structures that can "conduct" light, having certain spatial profile of the optical dielectric constant
2- Because the optical dielectric constant is not spatially uniform, waveguides allow only certain optical solutions (called optical modes) that: 
	- can be finite in space: confined or bounded modes
	- can extend to infinity, but only with certain relationship of the fields over the inhomogenous dielectric region
3- At the boundary between two regions of different ϵ, the boundary condition of the EM field critically determines the optical modes.

1.2.3 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])])

What if (Sqrt[(Subscript[ϵ, 1])-(Subscript[ϵ, 2])])>1?

Manipulate[
If[RUN==1, 
ϵ2=1.3  ;
Manipulate[Plot[{0.5-((Sqrt[ϵ1-ϵ2])/(Sqrt[1+ϵ2-ϵ1])) x,-0.5+((Sqrt[ϵ1-ϵ2])/(Sqrt[1+ϵ2-ϵ1])) x}, {x,-3,0}, Filling-> {1->{{2},{Hue[0.6,0.4,1]}}}
,PlotRange-> {{-3,0},{-3,3}}]
,{ϵ1 , 1.4,2.5}, AnimationRunning-> False]
,  ]
,{RUN,{0,1}}]

The quantity (Subscript[n, 0])Sin[(Subscript[θ, a])] is a measure of low large the the angular opening of the waveguide to accept (or to radiate) the beam. Thus, a concept is used to describe the quantity: numerical aperture or NA.

The light 
Skew rays propagate without passing through the center axis of the fiber.

The acceptance angle for skew rays is larger than the acceptance angle of meridional rays. This condition explains why skew rays outnumber meridional rays. Skew rays are often used in the calculation of light acceptance in an optical fiber. The addition of skew rays increases the amount of light capacity of a fiber. In large NA fibers, the increase may be significant. 
The addition of skew rays also increases the amount of loss in a fiber. Skew rays tend to propagate near the edge of the fiber core. A large portion of the number of skew rays that are trapped in the fiber core are considered to be leaky rays. Leaky rays are predicted to be totally reflected at the core-cladding boundary. However, these rays are partially refracted because of the curved nature of the fiber boundary. Mode theory is also used to describe this type of leaky ray loss.

The main result of skew ray is that because it can bounce around with a larger angle (glancing blow), with the skew angle γ, it has a larger acceptance angle:
                                           (Subscript[n, 0])Sin[(Subscript[θ, as])]=((Subscript[n, 1])/(Cos[γ])) (Sqrt[1-((((Subscript[n, 2])^2) )/((Subscript[n, 1])^2))])=(1/(Cos[γ]))(Sqrt[((Subscript[n, 1])^2)-((Subscript[n, 2])^2) ])
In other words, it increases by the term (1/(Cos[γ])).

2. Planar dielectric waveguide or slab waveguide

2.1 General concept of planar dielectric waveguide

A stack of layers of dielectric materials (arbitrary number of layers) can form a waveguide as long as it has EM solutions that propagate in the plane of the layers.

How do we know it functions as a waveguide? we must solve the wave equation (from the Maxwell's equations) and see if a propagating mode exists.

The general wave equation is:
(∇^2)E + ((ω^2)/(c^2))ϵ(x,y,z) E = 0              (2.1a)
(∇^2)H + ((ω^2)/(c^2))ϵ(x,y,z) H = 0            (2.1b)
In a waveguide, the dielectric is a constant along one dimension, say z. If the waveguide solution exists, then the electric field can be written:
E(x,y)  (e^(i (β z - ω t))), where β, called the propagation constant is unknown. ω is the frequency of light, which is given. The wave equation becomes:
(((∂^2)E)/(∂(x^2)))+(((∂^2)E)/(∂(y^2))) -((β^2)-((ω^2)/(c^2)) ϵ(x,y)) E=0             (2.2)
1- If a solution of Eq. (2.2) exists, it is called a mode of the structure 
2- Eq. (2.2) may have many solutions or none
3- A solution may have a unique value of β, or a continuous range of value of β
3- The mode is called a propagating wave if Re[β]≠0

Ex. 2.1 What does the wave (e^(i(β z - ω t))) look like for:
1. β=1
2. β=1+ i 0.05
2. β= i 0.05

Source code: Demo propagation constant

Manipulate[
If[RUN== 1,
ω =2 π  ;
Manipulate[Plot[ Re[(E^(I ((βR+I βI) z - ω t))) ],{z,0,50}, ImageSize-> {600,100}, PlotRange-> {{0,50},{-1,1}}, Frame-> True, AspectRatio-> 0.1, Filling-> Axis] 
,{{βR,1},{-1,0,1}},{βI,0,0.05}, {{t,0,Button["time",BaseStyle->{16, Bold,FontFamily-> "Calibri"}
,Background->Yellow]},0,1}  ]
,TextCell["Not running"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]]
,Row[{Button["Toggle box to start, and turn off when not running",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Black],Control[{{RUN,0,"->"},{0,1}}]
}]
]

Demo propagation constant

Ex. 2.2  Plot 3D the wave (e^(-(x^2)/(w^2)))  (e^(i(β z - ω t))) for:
β=1 and w=2.

Source code: Demo Gaussian confined propagation

Manipulate[
If[RUN== 1,ω =2 π  ;
Manipulate[
 
Animate[Plot3D[ (E^(-(x^2)/(w^2))) Re[(E^(I ((βR+I Sign[βR+0.00000001]βI) (z-(1-Sign[βR])25) - ω t))) ],{x,-10,10},{z,0,50}, (*  ImageSize-> {400,400},  *)Mesh-> False, PlotRange-> {{-10,10},{0,50},{-1.2,1.2}} 
,PlotPoints-> {21,150}, BoxRatios-> {2,5,1}, ViewPoint-> {5,8,4}]
, {{t,0,Button["time",BaseStyle->{16, Bold,FontFamily-> "Calibri"}
,Background->Yellow]},0,1}  ,AnimationRunning->False]
,{{βR,1},{-1,0,1}},{βI,0,0.05},{w,2,4}  ]
,TextCell["Not running"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]]
,Row[{Button["Toggle box to start, and turn off when not running",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Black],Control[{{RUN,0,"->"},{0,1}}]
}]
]

Demo Gaussian confined propagation

Manipulate[
DensityPlot[(E^(-(x^2)/(w^2))),{x,-5,5},{y,-10,10}, ColorFunction->GrayLevel, PlotRange->All], {w,0.5,3}]

Manipulate[
DensityPlot[(E^(-(x^2)/(wx^2))) (E^(-(y^2)/(wy^2))),{x,-5,5},{y,-5,5}, ColorFunction->GrayLevel,PlotPoints-> {30,30}, PlotRange->All], {wx,0.5,3}, {wy,0.5,3}]

2.2 An example: 3-layer symmetric slab waveguide

In practice, for complex design, Eq. (2.2) above can be solved by numerical method. Analytical approach is useful to illustrate the qualitative aspect of the dielectric waveguide behavior. Here, we will look at the simplest 3-layer symmetric slab waveguide

Consider the trivial case of symmetric 3 layers with indices (Subscript[ϵ, 1]) (cladding) and (Subscript[ϵ, 2]) (core).

2.2.1 Solutions in each layer

Step 1: We can find the general solutions in each layer

Equations (2.2) becomes
(((∂^2)E)/(∂(x^2))) -((β^2)-((ω^2)/(c^2)) ϵ(x)) E=0             (2.3)
which has two cases:
1- If x is in the core:
 (((∂^2)E)/(∂(x^2))) -((β^2)-((ω^2)/(c^2)) (Subscript[ϵ, 2])) E=0               (2.4.a)
2- If x is in the cladding:
 (((∂^2)E)/(∂(x^2))) -((β^2)-((ω^2)/(c^2)) (Subscript[ϵ, 1])) E=0               (2.4.b)
 Define: (Subscript[((k^2)), 1 or 2])=((ω^2)/(c^2)) (Subscript[ϵ, 1 or 2]) -(β^2)

DSolve[(∂_(x,x))Ef[x]+((Subscript[k, 2])^2) Ef[x]==0, Ef[x],x ]

TrigToExp[{{Ef(x)->(Subscript[c, 1]) sin((Subscript[k, 2]) x)+(Subscript[c, 2]) cos((Subscript[k, 2]) x)}}]

DSolve[(∂_(x,x))Ef[x]+((Subscript[k, 1])^2) Ef[x]==0, Ef[x],x ]

DSolve[(∂_(x,x))Ef[x]-((Subscript[k, 1])^2) Ef[x]==0, Ef[x],x ]

The solutions are:
                                (Subscript[ψ, core;A])[x]=(e^(i (Subscript[k, 2]) x));  (Subscript[ψ, core;B])[x]=(e^(-i (Subscript[k, 2]) x))      (2.5) 
in the core, where                   (Subscript[k, 2]) = (Sqrt[(Subscript[ϵ, 2]) ((Subscript[k, 0])^2)- (β^2)])                  (2.5a)
where                                                 ((Subscript[k, 0])^2)≡((ω^2)/(c^2))  
and
                                 (Subscript[ψ, clad;A])[x]=(e^(i (Subscript[k, 1]) x));  (Subscript[ψ, clad;B])[x]=(e^(-i (Subscript[k, 1]) x))       (2.6)
for radiating wave in the cladding;  (Subscript[k, 1]) = (Sqrt[(Subscript[ϵ, 1]) ((Subscript[k, 0])^2)- (β^2)])          (2.6a)
OR
                                 (Subscript[ψ, clad;A])[x]=(e^( (Subscript[κ, 1]) x));  (Subscript[ψ, clad;B])[x]=(e^(- (Subscript[κ, 1]) x))       (2.7)
for evanescent wave (bounded) in the cladding 
                                        (Subscript[κ, 1]) = (Sqrt[(β^2)-(Subscript[ϵ, 1]) ((Subscript[k, 0])^2)]) .                        (2.7a)
For waveguide mode (bounded), we expect only solutions Eq. (2.7)      
The E or H field (TE or TM):
                                          (E or H)=F[x](e^(i(β z-ω t))) (OverHat[y]).          (2.8a)
where:
                                      F[x]=(Subscript[c, A]) (Subscript[ψ, A])[x]+(Subscript[c, B]) (Subscript[ψ, B])[x]               (2.8b)
or:                                 F[x] = (({{(Subscript[ψ, A])[x], (Subscript[ψ, B])}})[x]) . (({{(Subscript[c, A])}, {(Subscript[c, B])}}))         (2.8c)
Remember that we don't have to keep the term (e^(i(β z-ω t) ))because it is the same everywhere.

2.2.2 Boundary condition

Step 2: A wave solution must satisfy boundary condition

We have the solutions (e^(±i (Subscript[k, 2]) x)) and and (e^(±i (Subscript[k, 1]) x)) or (e^(± (Subscript[κ, 1]) x)) , does this mean that the problem is solved? Obviously not, because we still don't have specific solutions. These are general solutions in EACH region. They must still match each other at the boundary.

For dielectric media, surface charge and surface current are zero. The boundary conditions are:
                                                           ∇. (OverVector[D]) = 0 ; ∇. (OverVector[B]) = 0 ;                    (2.9)
For normal component of E field: (Subscript[ε, 1]) (Subscript[E, 1]) = (Subscript[ε, 2]) (Subscript[E, 2]);  (D=ε E)              (2.10)
and H field:                                      (Subscript[μ, 1])(Subscript[H, 1]) = (Subscript[μ, 2]) (Subscript[H, 2])  (B=μ H)           (2.11)
For most nonmagnetic media, μ=1, so rarely do we have discontinuity of H across an interface.
For tangential E field, the component is continuous:
                                                            (Subscript[E, 1]) = (Subscript[E, 2])                           (2.12)
and for zero surface current, so is the H field:
                                                             (Subscript[H, 1]) = (Subscript[H, 2])                         (2.13)

Any polarization can be split into 2 components:  TE: transverse electric or TM: transverse magnetic

2.3 Boundary conditions for different polarizations

Each polarization, TE or TM has unique boundary condition

2.3.1 TE mode

Condition for boundary condition for TE polarization

Electric field is parallel to the interface: (Overscript[E, ->]) = (e^(i (β z-ω t))) F[x] (OverHat[y]):  there is only tangential component: that's why we call it transverse electric (TE) mode. What about the magnetic field? 
-(1/c)μ((∂ (Overscript[H, ->]))/(∂ t))= ∇×(Overscript[E, ->])  =  (({{(OverHat[x]), (OverHat[y]), (OverHat[z])}, {(∂_x), (∂_y), (∂_z)}, {(Subscript[E, x]), (Subscript[E, y]), (Subscript[E, z])}})) = (({{(OverHat[x]), (OverHat[y]), (OverHat[z])}, {(∂_x), (∂_y), (∂_z)}, {0, (Subscript[E, y]), 0}})) = (({{(OverHat[x]), (OverHat[y]), (OverHat[z])}}))(({{-i β F[x]}, {0}, {(∂_x)F[x]}}))(e^(i (β z-ω t)))  (2.14)
The normal boundary condition of H (x direction) yields the same as that of E tangential boundary condition: no surprise here. The tangential boundary condition of H (z direction) yields (∂_x)F[x] continuity. So if we have a general solution for F[x]: 
                                        F = (({{(Subscript[ψ, A]), (Subscript[ψ, B])}})) . (({{(Subscript[c, A])}, {(Subscript[c, B])}})) ;             (2.8c) above
the boundary conditions, which can be written: 
                                    (({{F}, {F'}}))= (({{(Subscript[ψ, A]), (Subscript[ψ, B])}, {(Subscript[ψ, A])', (Subscript[ψ, B])'}})) (({{(Subscript[c, A])}, {(Subscript[c, B])}}))              (2.15)
have to be continuous, or: 
                 (({{(Subscript[ψ, A]^((1))), (Subscript[ψ, B]^((1)))}, { (Subscript[ψ, A]^((1)))', (Subscript[ψ, B]^((1)))'}})) (({{(Subscript[c, 1A])}, {(Subscript[c, 1B])}})) = (({{(Subscript[ψ, A]^((2))), (Subscript[ψ, B]^((2)))}, { (Subscript[ψ, A]^((2)))', (Subscript[ψ, B]^((2)))'}})) (({{(Subscript[c, 2A])}, {(Subscript[c, 2B])}})).           (2.16)
From this, we can think of S-matrix approach again: 
                  (((({{(Subscript[ψ, A]^((2))), (Subscript[ψ, B]^((2)))}, { (Subscript[ψ, A]^((2)))', (Subscript[ψ, B]^((2)))'}})))^(-1))(({{(Subscript[ψ, A]^((1))), (Subscript[ψ, B]^((1)))}, { (Subscript[ψ, A]^((1)))', (Subscript[ψ, B]^((1)))'}})) (({{(Subscript[c, 1A])}, {(Subscript[c, 1B])}})) = (({{(Subscript[c, 2A])}, {(Subscript[c, 2B])}}))        (2.17)
or:                                  (({{(Subscript[S, 11]), (Subscript[S, 12])}, {(Subscript[S, 21]), (Subscript[S, 22])}})) (({{(Subscript[c, 1A])}, {(Subscript[c, 1B])}})) = (({{(Subscript[c, 2A])}, {(Subscript[c, 2B])}})) ;                 (2.18)
where: 
                               (({{(Subscript[S, 11]), (Subscript[S, 12])}, {(Subscript[S, 21]), (Subscript[S, 22])}})) =  (((({{(Subscript[ψ, A]^((2))), (Subscript[ψ, B]^((2)))}, { (Subscript[ψ, A]^((2)))', (Subscript[ψ, B]^((2)))'}})))^(-1))(({{(Subscript[ψ, A]^((1))), (Subscript[ψ, B]^((1)))}, { (Subscript[ψ, A]^((1)))', (Subscript[ψ, B]^((1)))'}}))    (2.19)  
This is called interface S-matrix.
(There are several approaches to defining S matrix: do not get confused. S matrix is a generic concept, not a specific definition)

2.3.2 TM mode

Condition for boundary condition for TE polarization

In this case, the magnetic field is parallel to the interface: (Overscript[H, ->]) = (e^(i (β z-ω t))) F[x] (OverHat[y]):  we call it transverse magnetic (TM) mode. What about the electric field? 
              (1/c)ϵ ((∂ (Overscript[E, ->]))/(∂ t))  = ∇×(Overscript[H, ->])  
               =  (({{(OverHat[x]), (OverHat[y]), (OverHat[z])}, {(∂_x) , (∂_y), (∂_z)}, {0, F[x](e^(i (β z-ω t))), 0}})) = (({{(OverHat[x]), (OverHat[y]), (OverHat[z])}}))(({{-i β F[x]}, {0}, {(∂_x)F[x]}}))(e^(i (β z-ω t)))  (2.20)
                                       ((-i ω)/c)ϵ (Overscript[E, ->])  =  (({{(OverHat[x]), (OverHat[y]), (OverHat[z])}}))(({{-i β F[x]}, {0}, {(∂_x)F[x]}}))(e^(i (β z-ω t)))                 (2.21)
The continuity of D normal gives F[x] continuity (but NOT electric field); the continuity of E tangential gives:
                               (1/(Subscript[ϵ, 1]))(∂_x)F[x](Subscript[|, 1]) = (1/(Subscript[ϵ, 2]))(∂_x)F[x](Subscript[|, 2]) 
                               F (x) derivative is discontinuous                                               (2.22)
Thus, given
                                     (({{F}, {F'}}))= (({{(Subscript[ψ, A]), (Subscript[ψ, B])}, {(Subscript[ψ, A])', (Subscript[ψ, B])'}})) (({{(Subscript[c, A])}, {(Subscript[c, B])}}))  ;                                        (2.15)  
The boundary value that we want to match is:
                                   (({{F}, {(1/ϵ)F'}})) = (({{1, 0}, {0, (1/ϵ)}})) (({{(Subscript[ψ, A]), (Subscript[ψ, B])}, {(Subscript[ψ, A])', (Subscript[ψ, B])'}})) (({{(Subscript[c, A])}, {(Subscript[c, B])}}))                            (2.23)
To write interface S-matrix:
                     (({{(Subscript[S, 11]), (Subscript[S, 12])}, {(Subscript[S, 21]), (Subscript[S, 22])}})) =  (((({{(Subscript[ψ, A]^((2))), (Subscript[ψ, B]^((2)))}, { (Subscript[ψ, A]^((2)))', (Subscript[ψ, B]^((2)))'}})))^(-1))(({{1, 0}, {0, (Subscript[ϵ, 2])}}))(({{1, 0}, {0, (1/(Subscript[ϵ, 1]))}}))(({{(Subscript[ψ, A]^((1))), (Subscript[ψ, B]^((1)))}, { (Subscript[ψ, A]^((1)))', (Subscript[ψ, B]^((1)))'}}))
                                        = (((({{(Subscript[ψ, A]^((2))), (Subscript[ψ, B]^((2)))}, { (Subscript[ψ, A]^((2)))', (Subscript[ψ, B]^((2)))'}})))^(-1))(({{1, 0}, {0, ((Subscript[ϵ, 2])/(Subscript[ϵ, 1]))}}))(({{(Subscript[ψ, A]^((1))), (Subscript[ψ, B]^((1)))}, { (Subscript[ψ, A]^((1)))', (Subscript[ψ, B]^((1)))'}}))               (2.24)

2.3.3 Equations for solutions

The solution is obtained by satisfying all boundary conditions. Each boundary condition gives a set of equation. Various coefficients will be determined by the Eqs.

2.4 Bounded mode TE

Characteristic equation for TE polarization: the equation that determines allowed (valid) propagation constants β for each polarization

2.4.1 Characteristic equation

The S matrix for the core:

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

Notice that the S matrix has translational symmetry.

This connect the left and the right of the core region:  
(Subscript[S, 2]) = (({{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]))}}))
For TE mode, the interface M matrix is unity. The eigen equation is:   
  (({{((Subscript[ψ, A])^(right most))[(Subscript[x, 1])], ((Subscript[ψ, B])^(right most))[(Subscript[x, 1])]}, {(∂_x) ((Subscript[ψ, A])^(right most))[(Subscript[x, 1])], (∂_x) ((Subscript[ψ, B])^(right most))[(Subscript[x, 1])]}})) . (({{(Subscript[c, An])}, {(Subscript[c, Bn])}})) =
	                                                               (Subscript[S, Total]) (({{((Subscript[ψ, A])^(left most))[(Subscript[x, 1])], ((Subscript[ψ, B])^(left most))[(Subscript[x, 1])]}, {(∂_x) ((Subscript[ψ, A])^(left most))[(Subscript[x, 1])], (∂_x) ((Subscript[ψ, B])^(left most))[(Subscript[x, 1])]}})) . (({{(Subscript[c, A1])}, {(Subscript[c, B1])}})) 
Here is a way to choose neat solution:
 ((Subscript[ψ, A])^(right most))[x] =(e^( (Subscript[κ, 1])(x-L))) ; ((Subscript[ψ, B])^(right most))[x] =(e^(- (Subscript[κ, 1])(x-L)  ))
 ((Subscript[ψ, A])^(left most))[x] =(e^( (Subscript[κ, 1])x)) ; ((Subscript[ψ, B])^(left most))[x] =(e^(- (Subscript[κ, 1])x ))
 (({{1, 1}, {(Subscript[κ, 1]), -(Subscript[κ, 1])}})).(({{0}, {(Subscript[c, right])}}))=(({{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]))}})). (({{1, 1}, {(Subscript[κ, 1]), -(Subscript[κ, 1])}})). (({{(Subscript[c, left])}, {0}}))

(({{0}, {(Subscript[c, right])}}))==FullSimplify[Inverse[ (({{1, 1}, {(Subscript[κ, 1]), -(Subscript[κ, 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]))}})). (({{1, 1}, {(Subscript[κ, 1]), -(Subscript[κ, 1])}}))]. (({{(Subscript[c, left])}, {0}}))

There are 2 equations. For the first one, the only way we have a non-trivial solution is:
                                                                 cos(L (Subscript[k, 2]))+((sin(L (Subscript[k, 2])) ((Subscript[κ, 1]^2)-(Subscript[k, 2]^2)))/(2 (Subscript[k, 2]) (Subscript[κ, 1])))=0.
This is the characteristic equation: it says that we can't just choose arbitrary value of β to satify the wave equation. There is a unique propagation constant for each mode.
 What does the second equation give us?
                                                    (Subscript[c, right])=(Subscript[c, left]) ((sin(L (Subscript[k, 2]))  ((Subscript[k, 2]^2)+(Subscript[κ, 1]^2)))/(2 (Subscript[k, 2]) (Subscript[κ, 1])))
a proportional relation between (Subscript[c, right]) and (Subscript[c, left]): we can choose an arbitrary (Subscript[c, left]), then (Subscript[c, right]) is fixed or vice versa. The wave amplitude is fixed relatively between sections, there is ONLY ONE arbitrary amplitude for the whole wave, as it should be.
There are several way to express the same characteristic equation:

((Cos[L (Subscript[k, 2])])/(Sin[L (Subscript[k, 2])]))==((((Subscript[k, 2]^2)-(Subscript[κ, 1]^2)))/(2 (Subscript[k, 2]) (Subscript[κ, 1])))

cos(L (Subscript[k, 2]))+((sin(L (Subscript[k, 2])) ((Subscript[κ, 1]^2)-(Subscript[k, 2]^2)))/(2 (Subscript[k, 2]) (Subscript[κ, 1])))==0 /. L->2 a

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

Factor[%]

(cos(a (Subscript[k, 2])) (Subscript[κ, 1])-sin(a (Subscript[k, 2])) (Subscript[k, 2])) (cos(a (Subscript[k, 2])) (Subscript[k, 2])+sin(a (Subscript[k, 2])) (Subscript[κ, 1]))==0

We notice that this means this single characteristic equation gives us 2 separate characteristic sub-equations: 
                                    cos(a (Subscript[k, 2])) (Subscript[κ, 1])-sin(a (Subscript[k, 2])) (Subscript[k, 2]) = 0 
AND                           cos(a (Subscript[k, 2])) (Subscript[k, 2])+sin(a (Subscript[k, 2])) (Subscript[κ, 1])=0
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: 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 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, the even solutions are: Cos[(Subscript[k, 2]) x]  in the core. At the core boundary:
(({{Cos[(Subscript[k, 2]) a]}, {-(Subscript[k, 2]) Sin[(Subscript[k, 2]) a]}})). In the cladding, the solution is C (e^(-(Subscript[κ, 1])(x-a) )): At the cladding boundary: (({{C}, {-C (Subscript[κ, 1])}})). Choose C=Cos[(Subscript[k, 2]) a], the equation is: (({{Cos[(Subscript[k, 2]) a]}, {-(Subscript[k, 2]) Sin[(Subscript[k, 2]) a]}}))= (({{Cos[(Subscript[k, 2]) a]}, {-Cos[(Subscript[k, 2]) a] (Subscript[κ, 1])}}))
Obviously, the charac. eq. is: (Subscript[k, 2]) Sin[(Subscript[k, 2]) a]=Cos[(Subscript[k, 2]) a] (Subscript[κ, 1]) or Cos[(Subscript[k, 2]) a] (Subscript[κ, 1])- (Subscript[k, 2]) Sin[(Subscript[k, 2]) a]=0 . We can similarly verify the same for odd solutions.

Numerical illustration of TE propagation constants

ResolvTE[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) κ1-sin(a k2) k2) (cos(a k2) k2+sin(a k2) κ1)  ] ;

indexprof={3.4,3.5 }  ; a=2.5 ; λ=1.5 ;
Plot[ResolvTE[indexprof,a,λ,x],{x,3.4,3.5}]

The roots of the equation are where the function crosses the horizontal axis, which is ==0 . These are discrete modes of the waveguide.

FindRoot[ResolvTE[{3.4,3.5 },2.5,1.5,n]==0,{n,3.42,3.44}]

Summary of key concepts:

1- This value n (example is 3.43683 above) is called the MODE EFFECTIVE INDEX.
                                             (Subscript[n, eff])=(β/(Subscript[k, 0]))
In other words, the mode acts as if it sees this index value in the waveguide.
2- Each mode has its own index, different from each other
3- The speed of propagation is:  v=(c/(Subscript[n, eff])) , hence different modes move at different speeds. This causes INTERMODAL DISPERSION. If a light pulse is a linear combination of many modes, each modal component will move at its own speed, and after awhile, be out of step with others.

2.4.2 EM field

From the propagation constant, the E field profile of the mode is determined.

So we find the mode. How do we find the E field?
Left region:  (({{F[x]}, {(∂_x) F[x]}}))=(({{(e^((Subscript[κ, 1]) (x+a))), (e^(-(Subscript[κ, 1]) (x+a)))}, {(Subscript[κ, 1])(e^((Subscript[κ, 1]) (x+a))), -(Subscript[κ, 1])(e^(-(Subscript[κ, 1]) (x+a)))}})). (({{(Subscript[c, left])}, {0}})) 
Middle or core region: (({{F[x]}, {(∂_x) F[x]}}))=(({{(e^(i (Subscript[k, 2]) x)), (e^(-i (Subscript[k, 2]) x))}, {i (Subscript[k, 2])(e^(i (Subscript[k, 2]) x)), -i (Subscript[k, 2])(e^(-i (Subscript[k, 2]) x))}})). (({{(Subscript[c, m1])}, {(Subscript[c, m2])}}))
(({{(e^(i (Subscript[k, 2]) x)), (e^(-i (Subscript[k, 2]) x))}, {i (Subscript[k, 2])(e^(i (Subscript[k, 2]) x)), -i (Subscript[k, 2])(e^(-i (Subscript[k, 2]) x))}})). (({{(Subscript[c, m1])}, {(Subscript[c, m2])}}))(Subscript[|, x→-(SuperPlus[a])]) =(({{(e^((Subscript[κ, 1]) (x+a))), (e^(-(Subscript[κ, 1]) (x+a)))}, {(Subscript[κ, 1])(e^((Subscript[κ, 1]) (x+a))), -(Subscript[κ, 1])(e^(-(Subscript[κ, 1]) (x+a)))}})). (({{(Subscript[c, left])}, {0}}))(Subscript[|, x→-(SuperMinus[a])]) 
(({{(Subscript[c, m1])}, {(Subscript[c, m2])}}))=Inverse[ (({{(e^(i (Subscript[k, 2]) x)), (e^(-i (Subscript[k, 2]) x))}, {i (Subscript[k, 2])(e^(i (Subscript[k, 2]) x)), -i (Subscript[k, 2])(e^(-i (Subscript[k, 2]) x))}}))(Subscript[|, x→-(SuperPlus[a])])] .(({{(e^((Subscript[κ, 1]) (x+a))), (e^(-(Subscript[κ, 1]) (x+a)))}, {(Subscript[κ, 1])(e^((Subscript[κ, 1]) (x+a))), -(Subscript[κ, 1])(e^(-(Subscript[κ, 1]) (x+a)))}}))(Subscript[|, x→-(SuperMinus[a])]) . (({{(Subscript[c, left])}, {0}}))
Remember interface S matrix?
S=  (((({{(e^(i (Subscript[k, 2]) (Subscript[x, n]))), (e^(-i (Subscript[k, 2]) (Subscript[x, n])))}, {i (Subscript[k, 2])(e^(i (Subscript[k, 2]) (Subscript[x, n]))), -i (Subscript[k, 2])(e^(-i (Subscript[k, 2]) (Subscript[x, n])))}})))^(-1)) .(({{(e^((Subscript[κ, 1]) (Subscript[x, n]))), (e^(-(Subscript[κ, 1]) (Subscript[x, n])))}, {(Subscript[κ, 1])(e^((Subscript[κ, 1]) (Subscript[x, n]))), -(Subscript[κ, 1])(e^(-(Subscript[κ, 1]) (Subscript[x, n])))}}))
So here is how to get the full solution:
Chose (Subscript[c, left]) = 1 (or anything we like, except for 0)
Use S interface matrix to get to the next region, and the next, and so on... But: use the "second" equation to get the last region.
Finally: E=F[x](e^(i(β z-ω t)))
((i ω)/c)μ(Overscript[H, ->])=  (({{(OverHat[x]), (OverHat[y]), (OverHat[z])}}))(({{-i β F[x]}, {0}, {(∂_x)F[x]}}))(e^(i (β z-ω t)))
(({{(Subscript[H, x])}, {(Subscript[H, z])}}))=(1/(Subscript[k, 0]))(({{-β F[x]}, {-i (∂_x)F[x]}})) (e^(i (β z-ω t)))

FullSimplify[Inverse[ (({{(E^(I (Subscript[k, 2]) x)), (E^(-I (Subscript[k, 2]) x))}, {I (Subscript[k, 2])(E^(I (Subscript[k, 2]) x)), -I (Subscript[k, 2])(E^(-I (Subscript[k, 2]) x))}}))] .(({{(E^((Subscript[κ, 1]) (x+a))), (E^(-(Subscript[κ, 1]) (x+a)))}, {(Subscript[κ, 1])(E^((Subscript[κ, 1]) (x+a))), -(Subscript[κ, 1])(E^(-(Subscript[κ, 1]) (x+a)))}}))/. x->-a]

(({{(((E^(I a (Subscript[k, 2]))) ((Subscript[k, 2])-I (Subscript[κ, 1])))/(2 (Subscript[k, 2]))), (((E^(I a (Subscript[k, 2]))) ((Subscript[k, 2])+I (Subscript[κ, 1])))/(2 (Subscript[k, 2])))}, {(((E^(-I a (Subscript[k, 2]))) ((Subscript[k, 2])+I (Subscript[κ, 1])))/(2 (Subscript[k, 2]))), (((E^(-I a (Subscript[k, 2]))) ((Subscript[k, 2])-I (Subscript[κ, 1])))/(2 (Subscript[k, 2])))}})).(({{1}, {0}}))

FullSimplify[ (({{(E^(I (Subscript[k, 2]) x)), (E^(-I (Subscript[k, 2]) x))}, {I (Subscript[k, 2])(E^(I (Subscript[k, 2]) x)), -I (Subscript[k, 2])(E^(-I (Subscript[k, 2]) x))}})).(({{(((E^(I a (Subscript[k, 2]))) ((Subscript[k, 2])-I (Subscript[κ, 1])))/(2 (Subscript[k, 2])))}, {(((E^(-I a (Subscript[k, 2]))) ((Subscript[k, 2])+I (Subscript[κ, 1])))/(2 (Subscript[k, 2])))}}))]

So the solutions are:  
 If[Abs[x]<a, (* in the core region *) 
		{ 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  *)
Left side: (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  *) 
Right side: x>0
Sign[ sin(2 a k2)   ] {(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  *)
For practical use, we still need to normalize:

FullSimplify[Coreintegration=(∫_(-a)^a)(((cos((a+x)k2)+((sin((a+x) k2) κ1)/k2)))^2)\[DifferentialD]x]

FullSimplify[Cladintegration=(∫_a^∞)((((E^(-κ1 (x-a) ))))^2)\[DifferentialD]x]

N factor: (((-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))

2.4.3 Illustration: what does the "profile" looks like

We can plot F[x] which is directly proportional to Ey profile for various mode

nindex={3.4,3.5}; a=2.5 ; λ=1.5 ;
Plot[ResolvTE[nindex,a,λ,x],{x,nindex[[1]],nindex[[2]]}]

nindex={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<7, i++,
	AppendTo[neff, 
	 x /. FindRoot[ResolvTE[nindex,a,λ,x]==0,{x,segment[[i+1]],segment[[i]]}]
                            ] ]

neff

allp={};
nindex={3.4,3.5}; a=2.5 ; λ=1.5 ;
For[i=1,i<Length[neff]+1,i++,
parm=ParmTE[nindex,a,λ,neff[[i]] ]  ;
AppendTo[allp,  Plot[FieldTE2[parm,x][[1]],{x,-7,7}
	, PlotRange->{-1.,1}
, PlotStyle->{Thickness[0.02],Hue[i*0.1,1,0.8]}
,Filling-> {1->{ Axis,{Hue[i*0.1,1,0.8,0.5]}}}
, Ticks-> None, AspectRatio-> 1.] 
]
             ] ;
Show[GraphicsRow[allp] ]

2.4.4 Traveling electric field

The field is essentially a standing wave in the x direction (because of waveguide) and a traveling wave in the z direction
Link to 3.1.2.1

imode=1 ;  β=2 π *neff[[imode]]/λ  ; ω= 2 π ; 
parm=ParmTE[nindex,a,λ,neff[[imode]] ]  ;
allp={} ;
For[i=1,i≤ 40, i++, t=0.025*i ;
PrintTemporary["i=   ",i];
AppendTo[allp,Plot3D[FieldTE2[parm,x][[1]]  Cos[β z - ω t] ,{x,-4,4},{z,0,2},   ImageSize-> {400,200}
,PlotPoints->{15,40},  Mesh-> False, PlotRange-> {{-4,4},{0,2},{-1,1}} 
, BoxRatios-> {4,2,1}, ViewPoint-> {5,8,6}] ]
];
Animate[allp[[i]]  ,{i,1, Length[allp],1}, AnimationRunning->False]

Link to Fig .1.2.1

2.4.5 Vector graphics

Remember that EM field are vector field. Not scalar field. We can look at the amplitude of a single component, but to have a sense of what really happens, we must look at vector field.

Linked to Fig. 3.1.3.1 Magnetic field lines of mode 0

What does the Hz component look like (compared with Hy?)

Plot[ {FieldTE[nindex,a,λ,neff[[1]],x][[1]]
		       ,Im[FieldTE[nindex,a,λ,neff[[1]],x]] [[2]]
				,-FieldTE[nindex,a,λ,neff[[1]],x][[3]]},
		        {x,-7,7}
	, PlotRange->All
	, PlotStyle->{{Hue[0,1,0.8]}
			                      , {Hue[0.4,1,0.8]}
			                       ,{Hue[0.8,1,0.8]}  }]

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

Link to Fig. 3.1.3.2

Link to 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

2.4.6 Field lines

Field line is a concept that "draws" the 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->100,ContourStyle-> {Hue[i*0.045,1,0.9], Thickness[0.005]}, DisplayFunction->Identity]]]
Show[fl, DisplayFunction->$DisplayFunction]

Link to section 3.1.1

Profile of all components: Ey, Hx, and Hz
Link to section 3.1.1 fig 3.1.3.a

2.5 Bounded mode TM

Characteristic equation for TM polarization: the equation that determines allowed (valid) propagation constants β for each polarization

2.3.1 Characteristic equation

The S matrix for the core: (Subscript[S, 2]) = (({{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]))}}))

For TE mode, the interface M matrix is unity. The eigen equation is:  
  M.(({{((Subscript[ψ, A])^(right most))[(Subscript[x, 1])], ((Subscript[ψ, B])^(right most))[(Subscript[x, 1])]}, {(∂_x) ((Subscript[ψ, A])^(right most))[(Subscript[x, 1])], (∂_x) ((Subscript[ψ, B])^(right most))[(Subscript[x, 1])]}})) . (({{(Subscript[c, An])}, {(Subscript[c, Bn])}})) =
	                                                                (Subscript[S, Total]) .M.(({{((Subscript[ψ, A])^(left most))[(Subscript[x, 1])], ((Subscript[ψ, B])^(left most))[(Subscript[x, 1])]}, {(∂_x) ((Subscript[ψ, A])^(left most))[(Subscript[x, 1])], (∂_x) ((Subscript[ψ, B])^(left most))[(Subscript[x, 1])]}})) . (({{(Subscript[c, A1])}, {(Subscript[c, B1])}})) 
Here is a way to choose neat solution:
 ((Subscript[ψ, A])^(right most))[x] =(e^( (Subscript[κ, 1])(x-L))) ; ((Subscript[ψ, B])^(right most))[x] =(e^(- (Subscript[κ, 1])(x-L)  ))
 ((Subscript[ψ, A])^(left most))[x] =(e^( (Subscript[κ, 1])x)) ; ((Subscript[ψ, B])^(left most))[x] =(e^(- (Subscript[κ, 1])x ))
 (({{1, 1}, {(Subscript[κ, 1]), -(Subscript[κ, 1])}})).(({{0}, {(Subscript[c, right])}}))=(({{1, 0}, {0, ((Subscript[ϵ, 1])/(Subscript[ϵ, 2]))}}))(({{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]))}})).(({{1, 0}, {0, ((Subscript[ϵ, 2])/(Subscript[ϵ, 1]))}})). (({{1, 1}, {(Subscript[κ, 1]), -(Subscript[κ, 1])}})). (({{(Subscript[c, left])}, {0}}))

(({{0}, {(Subscript[c, right])}}))==FullSimplify[Inverse[ (({{1, 1}, {(Subscript[κ, 1]), -(Subscript[κ, 1])}}))].(({{1, 0}, {0, ((Subscript[ϵ, 1])/(Subscript[ϵ, 2]))}})).(({{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]))}})).(({{1, 0}, {0, ((Subscript[ϵ, 2])/(Subscript[ϵ, 1]))}})). (({{1, 1}, {(Subscript[κ, 1]), -(Subscript[κ, 1])}}))]. (({{(Subscript[c, left])}, {0}}))

There are 2 equations. For the first one, the only way we have a non-trivial solution is:
                                              cos(L (Subscript[k, 2]))+sin(L (Subscript[k, 2]))((( (Subscript[ϵ, 2]) (Subscript[κ, 1]))/(2 (Subscript[k, 2]) (Subscript[ϵ, 1])))-(((Subscript[k, 2]) (Subscript[ϵ, 1]))/(2 (Subscript[ϵ, 2]) (Subscript[κ, 1]))))=0.
This is the characteristic equation: it says that we can't just choose arbitrary value of β. There is a unique propagation constant for each mode.
 What does the second equation give us? 
                                                (Subscript[c, right])=(Subscript[c, left]) sin(L (Subscript[k, 2])) ((((Subscript[k, 2]^2) (Subscript[ϵ, 1]^2)+(Subscript[ϵ, 2]^2) (Subscript[κ, 1]^2)))/(2 (Subscript[k, 2]) (Subscript[ϵ, 1]) (Subscript[ϵ, 2]) (Subscript[κ, 1])))

Similar as above, we can decompose the characteristic equation

TrigExpand[cos(L (Subscript[k, 2]))+sin(L (Subscript[k, 2]))((( (Subscript[ϵ, 2]) (Subscript[κ, 1]))/(2 (Subscript[k, 2]) (Subscript[ϵ, 1])))-(((Subscript[k, 2]) (Subscript[ϵ, 1]))/(2 (Subscript[ϵ, 2]) (Subscript[κ, 1])))) /. L->2 a]

Factor[%]

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)  ] ;

nindex={3.4,3.5 }; a=2.5 ; λ=1.5 ;
Plot[ResolvTM[nindex,a,λ,x],{x,3.4,3.5}]

FindRoot[ResolvTM[{3.4,3.5 },2.5,1.5,x]==0,{x,3.405,3.42}]

2.3.2 EM field

From the propagation constant, the E field profile of the mode is determined.

How do we find the EM field?
(Overscript[H, ->]) = (e^(i (β z-ω t))) F[x] (OverHat[y])
((-i ω)/c)ϵ (Overscript[E, ->])  =  (({{(OverHat[x]), (OverHat[y]), (OverHat[z])}}))(({{-i β F[x]}, {0}, {(∂_x)F[x]}}))(e^(i (β z-ω t))) 

Left region:  (({{F[x]}, {(∂_x) F[x]}}))=(({{(e^((Subscript[κ, 1]) (x+a))), (e^(-(Subscript[κ, 1]) (x+a)))}, {(Subscript[κ, 1])(e^((Subscript[κ, 1]) (x+a))), -(Subscript[κ, 1])(e^(-(Subscript[κ, 1]) (x+a)))}})). (({{(Subscript[c, left])}, {0}})) 
Middle or core region: (({{F[x]}, {(∂_x) F[x]}}))=(({{(e^(i (Subscript[k, 2]) x)), (e^(-i (Subscript[k, 2]) x))}, {i (Subscript[k, 2])(e^(i (Subscript[k, 2]) x)), -i (Subscript[k, 2])(e^(-i (Subscript[k, 2]) x))}})). (({{(Subscript[c, m1])}, {(Subscript[c, m2])}}))
(Subscript[M, 2]).(({{(e^(i (Subscript[k, 2]) x)), (e^(-i (Subscript[k, 2]) x))}, {i (Subscript[k, 2])(e^(i (Subscript[k, 2]) x)), -i (Subscript[k, 2])(e^(-i (Subscript[k, 2]) x))}})). (({{(Subscript[c, m1])}, {(Subscript[c, m2])}}))(Subscript[|, x→-(SuperPlus[a])]) =(Subscript[M, 1]). (({{(e^((Subscript[κ, 1]) (x+a))), (e^(-(Subscript[κ, 1]) (x+a)))}, {(Subscript[κ, 1])(e^((Subscript[κ, 1]) (x+a))), -(Subscript[κ, 1])(e^(-(Subscript[κ, 1]) (x+a)))}})). (({{(Subscript[c, left])}, {0}}))(Subscript[|, x→-(SuperMinus[a])]) 
(({{(Subscript[c, m1])}, {(Subscript[c, m2])}}))=Inverse[ (({{(e^(i (Subscript[k, 2]) x)), (e^(-i (Subscript[k, 2]) x))}, {i (Subscript[k, 2])(e^(i (Subscript[k, 2]) x)), -i (Subscript[k, 2])(e^(-i (Subscript[k, 2]) x))}}))(Subscript[|, x→-(SuperPlus[a])])] .(Subscript[M, T]).(({{(e^((Subscript[κ, 1]) (x+a))), (e^(-(Subscript[κ, 1]) (x+a)))}, {(Subscript[κ, 1])(e^((Subscript[κ, 1]) (x+a))), -(Subscript[κ, 1])(e^(-(Subscript[κ, 1]) (x+a)))}}))(Subscript[|, x→-(SuperMinus[a])]) . (({{(Subscript[c, left])}, {0}}))
Remember interface S matrix?
S=  (((({{(e^(i (Subscript[k, 2]) (Subscript[x, n]))), (e^(-i (Subscript[k, 2]) (Subscript[x, n])))}, {i (Subscript[k, 2])(e^(i (Subscript[k, 2]) (Subscript[x, n]))), -i (Subscript[k, 2])(e^(-i (Subscript[k, 2]) (Subscript[x, n])))}})))^(-1)).(Subscript[M, T]) .(({{(e^((Subscript[κ, 1]) (Subscript[x, n]))), (e^(-(Subscript[κ, 1]) (Subscript[x, n])))}, {(Subscript[κ, 1])(e^((Subscript[κ, 1]) (Subscript[x, n]))), -(Subscript[κ, 1])(e^(-(Subscript[κ, 1]) (Subscript[x, n])))}}))
So here is how to get the full solution:
Chose (Subscript[c, left]) = 1 (or anything we like, except for 0)
Use S interface matrix to get to the next region, and the next, and so on... But: use the "second" equation to get the last region.  
Again:  
(Overscript[H, ->]) = (e^(i (β z-ω t))) F[x] (OverHat[y])
((-i ω)/c)ϵ (Overscript[E, ->])  =  (({{(OverHat[x]), (OverHat[y]), (OverHat[z])}}))(({{-i β F[x]}, {0}, {(∂_x)F[x]}}))(e^(i (β z-ω t))) 
(({{(Subscript[E, x])}, {(Subscript[E, z])}}))=(1/(ϵ (Subscript[k, 0])))(({{ β F[x]}, {i (∂_x)F[x]}})) (e^(i (β z-ω t)))

FullSimplify[Inverse[ (({{(E^(I (Subscript[k, 2]) x)), (E^(-I (Subscript[k, 2]) x))}, {I (Subscript[k, 2])(E^(I (Subscript[k, 2]) x)), -I (Subscript[k, 2])(E^(-I (Subscript[k, 2]) x))}}))] .(({{1, 0}, {0, ((Subscript[ϵ, 2])/(Subscript[ϵ, 1]))}})).(({{(E^((Subscript[κ, 1]) (x+a))), (E^(-(Subscript[κ, 1]) (x+a)))}, {(Subscript[κ, 1])(E^((Subscript[κ, 1]) (x+a))), -(Subscript[κ, 1])(E^(-(Subscript[κ, 1]) (x+a)))}}))/. x->-a]

(({{(((E^(I a (Subscript[k, 2]))) ((Subscript[k, 2]) (Subscript[ϵ, 1])-I (Subscript[ϵ, 2]) (Subscript[κ, 1])))/(2 (Subscript[k, 2]) (Subscript[ϵ, 1]))), (((E^(I a (Subscript[k, 2]))) ((Subscript[k, 2]) (Subscript[ϵ, 1])+I (Subscript[ϵ, 2]) (Subscript[κ, 1])))/(2 (Subscript[k, 2]) (Subscript[ϵ, 1])))}, {(((E^(-I a (Subscript[k, 2]))) ((Subscript[k, 2]) (Subscript[ϵ, 1])+I (Subscript[ϵ, 2]) (Subscript[κ, 1])))/(2 (Subscript[k, 2]) (Subscript[ϵ, 1]))), (((E^(-I a (Subscript[k, 2]))) ((Subscript[k, 2]) (Subscript[ϵ, 1])-I (Subscript[ϵ, 2]) (Subscript[κ, 1])))/(2 (Subscript[k, 2]) (Subscript[ϵ, 1])))}})).(({{1}, {0}}))

FullSimplify[ (({{(E^(I (Subscript[k, 2]) x)), (E^(-I (Subscript[k, 2]) x))}, {I (Subscript[k, 2])(E^(I (Subscript[k, 2]) x)), -I (Subscript[k, 2])(E^(-I (Subscript[k, 2]) x))}})).(({{(((E^(I a (Subscript[k, 2]))) ((Subscript[k, 2]) (Subscript[ϵ, 1])-I (Subscript[ϵ, 2]) (Subscript[κ, 1])))/(2 (Subscript[k, 2]) (Subscript[ϵ, 1])))}, {(((E^(-I a (Subscript[k, 2]))) ((Subscript[k, 2]) (Subscript[ϵ, 1])+I (Subscript[ϵ, 2]) (Subscript[κ, 1])))/(2 (Subscript[k, 2]) (Subscript[ϵ, 1])))}}))]

Normalization of Ex

FullSimplify[Coreintegration=(∫_(-a)^a)(((cos((a+x)k2)+((sin((a+x) k2) ϵ2  κ1)/(k2 ϵ1))))^2)\[DifferentialD]x]

((((neff/(ϵ2))))^2)((( (ϵ1 ϵ2 κ1-ϵ1 ϵ2 cos(4 a k2) κ1+2 a ((k2^2) ((ϵ1)^2)+((ϵ2)^2) ((κ1)^2)))/(2 (k2^2) ((ϵ1)^2)))+(((k2 ϵ1-ϵ2 κ1) (k2 ϵ1+ϵ2 κ1) sin(4 a k2))/(4 (k2^3) ((ϵ1)^2))))

FullSimplify[Cladintegration=(∫_a^∞)((((E^(-κ1 (x-a) ))))^2)\[DifferentialD]x]

nindex={1,3.5}; a=0.5 ; λ=1.5 ;
Plot[ResolvTM[nindex,a,λ,x],{x,nindex[[1]],nindex[[2]]}]

2.3.3 The most important feature of TM mode: the E field discontinuity

What does the trasnverse magnetic field Hy look like?
Linked to Fig. 3.2.1.1

But remember the normal electric field Ex boundary condition:
(Subscript[ϵ, 1]) (Subscript[E, x;1])=(Subscript[ϵ, 2]) (Subscript[E, x;2 ])
Linked to Fig. 3.2.1.2

An important concept (for all waveguide, not just slab): confinement factor: How much of the electric field is really in the core region? We'll see about this in later section.

F[x] is directly proportional to Hy profile for various mode

nindex={1,3.5}; a=0.5 ; λ=1.5 ;
Plot[ResolvTM[nindex,a,λ,x],{x,nindex[[1]],nindex[[2]]}]

segment2={3.49,3.4,3,2.5,1.5,1.01}; neff={} ;
For[ i=1, i<Length[segment2], i++,
	AppendTo[neff, 
	 x /. FindRoot[ResolvTM[nindex,a,λ,x]==0,{x,segment2[[i+1]],segment2[[i]]}]
                            ] ]

neff

Let's look at the transverse magnetic Hy

allp={};
For[i=1,i≤ Length[neff],i++,
parm=ParmTM[nindex,a,λ,neff[[i]] ]  ;
AppendTo[allp,  Plot[FieldTM2[parm,x][[1]],{x,-2,2}
	, PlotRange->{-10.,10}
, PlotStyle->{Thickness[0.02],Hue[i*0.1,1,0.8]}
,Filling-> {1->{ Axis,{Hue[i*0.1,1,0.8,0.5]}}}
, Ticks-> None, AspectRatio-> 0.75] 
]
             ] ;
Show[GraphicsRow[allp] ]

Let's look at the Electric field Ex

allp={};
For[i=1,i≤ Length[neff],i++,
parm=ParmTM[nindex,a,λ,neff[[i]] ]  ;
AppendTo[allp,  Plot[FieldTM2[parm,x][[3]],{x,-2,2}
	, PlotRange->{-2.,2} 
, PlotStyle->{Thickness[0.012],Hue[i*0.1,1,0.8]}
,Filling-> {1->{ Axis,{Hue[i*0.1,1,0.8,0.5]}}}
, Ticks-> None, AspectRatio-> 0.75] 
]
             ] ;
Show[GraphicsRow[allp] ]

Fig. 3.2.1.2

Compare Hy and Ex of the lowest order mode:

Plot[{Re[FieldTM[nindex,a,λ,neff[[1]],x][[1]] ],
	2*Re[FieldTM[nindex,a,λ,neff[[1]],x][[3]] ]},{x,-1,1}
,PlotStyle->{{RGBColor[1.,0.,0]}, {RGBColor[0.,0.,1]}}]

Compare Hy and Ex of the highest order mode:

Plot[{Re[FieldTM[nindex,a,λ,neff[[5]],x][[1]] ],
	2*Re[FieldTM[nindex,a,λ,neff[[5]],x][[3]] ]},{x,-4,4}
	, PlotRange->{{-4,4},All}
,PlotStyle->{{RGBColor[1.,0.,0]}, {RGBColor[0.,0.,1]}}
,AspectRatio-> 0.4]

Linked to reviewed materials

3. Key concepts of waveguide

Although we study a particular waveguide geometry above, the slab waveguide, several important concepts are applicable to any waveguide, and can be illustrated with the slab waveguide.

3.1 Intensity and energy flow

The use of waveguide is to conduct lightwave: energy flow

3.1.1 Energy flow: Poynting vector - TE mode

Energy flow for TE mode

Intensity inside a waveguide is obtained by evaluating the Poynting vector. The time-averaged Poynting vector is:  
                                            ⟨(Overscript[S, →])⟩=(c/(8 π))Re[ (Overscript[E, ->]) ×(Overscript[(SuperStar[H]), ->])]                         (3.1)

Consider the example of slab waveguide.
For TE mode:
                                  (Overscript[E, ->]) =  E[x,z] (OverHat[y ])(e^(-i ω t ))                                                  (3.2)
                                (Overscript[H, ->])=-(i/(μ (Subscript[k, 0]))) ∇×(Overscript[E, ->])  = (({{(OverHat[x]), (OverHat[y]), (OverHat[z])}, {(∂_x), (∂_y), (∂_z)}, {0, (Subscript[E, y]), 0}})) 
=-(i/(μ (Subscript[k, 0])))(({{(OverHat[x]), (OverHat[y]), (OverHat[z])}}))(({{-(∂_z)(Subscript[E, y])[x,z]}, {0}, {(∂_x)(Subscript[E, y])[x,z]}}))(e^(-i ω t))                    (3.3)
                   ⟨(Overscript[S, →])⟩=(c/(8 π))Re[ (Overscript[E, ->]) ×(Overscript[(SuperStar[H]), ->])] 
           =(c/(8 π)) (1/(μ (Subscript[k, 0])))Re|({{(OverHat[x]), (OverHat[y]), (OverHat[z])}, {0, (Subscript[E, y]), 0}, {-i(∂_z)(SuperStar[(Subscript[E, y])])[x,z], 0, i (∂_x)(SuperStar[(Subscript[E, y])])[x,z]}})|       (3.4)

Det[(({{(OverHat[x]), (OverHat[y]), (OverHat[z])}, {0, (Subscript[Ef, y])[x,z], 0}, {-I(∂_z)(SuperStar[(Subscript[Ef, y])])[x,z], 0, I (∂_x)(SuperStar[(Subscript[Ef, y])])[x,z]}}))]

⟨(Overscript[S, →])⟩=c/(8 π)Re[ (Overscript[E, ->]) ×(Overscript[(SuperStar[H]), ->])]
        =(c/(8 π μ (Subscript[k, 0])))Re[ (OverHat[x]) i  (Subscript[E, y])(x,z)  ((∂ (SuperStar[(Subscript[E, y])]))/(∂ x)) +(OverHat[z]) i (Subscript[E, y])(x,z) ((∂ (SuperStar[(Subscript[E, y])]))/(∂ z))]         (3.5)
Is it possible for a net energy to flow along x direction? Intuitively, no. The energy can bounce back and forth from within the waveguide, but the net roundtrip (or the total space average) must be zero. 
If (Subscript[E, y])(x,z) is real, then so is ((∂ (SuperStar[(Subscript[E, y])]))/(∂ x))  and Re[ (OverHat[x]) i  (Subscript[E, y])(x,z) ((∂ (SuperStar[(Subscript[E, y])]))/(∂ x))]= 0: Obviously the intensity along x direction = 0. But the principle must be true much more generally, not for just this case of TE mode and real (Subscript[E, y])(x,z).
The intensity of the z-traveling component is given by:
                       |⟨(Overscript[S, →])⟩.(OverHat[z])| =(c/(8 π μ (Subscript[k, 0]))) |Re[ i (Subscript[E, y])(x,z) ((∂ (SuperStar[(Subscript[E, y])]))/(∂ z))]|                            (3.6)
and the power is given by:  
            P=(∫_(-∞)^∞)|⟨(Overscript[S, →])⟩.(OverHat[z])|dx =(c/(8 π μ (Subscript[k, 0]))) (∫_(-∞)^∞)|Re[ i (Subscript[E, y])(x,z) ((∂ (SuperStar[(Subscript[E, y])]))/(∂ z))]|dx        (3.7)
For mode m:   
                                   ((∂ (SuperStar[(Subscript[E, y])]))/(∂ z))= -i (Subscript[β, m]) (e^(-i (Subscript[β, m]) z)) (SuperStar[(Subscript[E, m])])[x]                                      (3.8)
The intensity is:
                            |⟨(Overscript[S, →])⟩.(OverHat[z])| =(c/(8 π μ (Subscript[k, 0])))(Subscript[β, m]) (|(Subscript[E, m])[x]|)^2 
                  =((c (Subscript[n, eff;m]))/(8 π μ))  (|(Subscript[E, m])[x]|)^2                                 (3.9)
                                       P= ((c (Subscript[n, eff;m]))/(8 π μ))(∫_(-∞)^∞)(((|(Subscript[E, m])[x]|))^2)dx                               (3.10)
These expressions are very similar to planewave, with the difference being the effective index, which is mode-dependent.
This is the reason why normalization with respect to power is slightly different from normalization with respect to wavefunction:
For just wavefunction:             (∫_(-∞)^∞)(|(Subscript[E, m])[x]|)^2 = 1.                            (3.11a)
For power:                         (Subscript[n, eff;m])(∫_(-∞)^∞)(|(Subscript[E, m])[x]|)^2dx = 1                    (3.11b)
(all other constants are dropped for convenient).
Either way is fine as long as one remembers what to use in each calculation.

3.1.2 For TM mode

Energy flow for TM mode

⟨(Overscript[S, →])⟩=(c/(8 π)) (1/(ϵ  (Subscript[k, 0]))) Re[-i (OverHat[x]) (SuperStar[(Subscript[H, y])]) (Subscript[H, y]^((1,0)))(x,z) -i (OverHat[z]) (SuperStar[(Subscript[H, y])]) (Subscript[H, y]^((0,1)))(x,z)]       (3.12)
The intensity is given by:  
                                        |⟨(Overscript[S, →])⟩.(OverHat[z])| =(c/(8 π ϵ  (Subscript[k, 0]))) |Re[ -i (SuperStar[(Subscript[H, y])])(x,z) ((∂ (Subscript[H, y]))/(∂ z))]|            (3.13)
Now we need to convert to Ex:
                                  (Subscript[E, x]) =- (i/(ϵ  (Subscript[k, 0])))(∂_z)(Subscript[H, y])[x,z] .                           (3.14)
 For a given mode m: 
                                     (Subscript[E, m;x]) =- (i/(ϵ  (Subscript[k, 0])))i (Subscript[β, m]) (Subscript[H, y])[x,z] 
            =((Subscript[β, m])/(ϵ  (Subscript[k, 0]))) (Subscript[H, y])[x,z]  = ((Subscript[n, eff;m])/ϵ) (Subscript[H, y])[x,z]          (3.15)
Thus:                                    (Subscript[H, y])[x,z]= (ϵ/(Subscript[n, eff;m]))(Subscript[E, m;x ]).                               (3.16)
Then: 
                                    |⟨(Overscript[S, →])⟩.(OverHat[z])| =(c/(8 π )) (ϵ/(Subscript[n, eff;m]))|Re[ (SuperStar[(Subscript[E, m;x])]) (Subscript[E, m;x])]|
            = (c/(8 π )) (ϵ/(Subscript[n, eff;m]))(((|(Subscript[E, m;x])|))^2)               (3.17)
Indeed that the intensity for TM mode is proportional to the transverse component of the electric field, as expected. 
The power of the mth mode is given by:  
                       (Subscript[P, m])=(∫_(-∞)^∞)|⟨(Overscript[(Subscript[S, m]), →])⟩.(OverHat[z])|dx = (c/(8 π (Subscript[n, eff;m ])))(∫_(-∞)^∞)ϵ[x](((|(Subscript[E, m;x])|))^2)dx       (3.18)
Notice that (∫_(-∞)^∞)ϵ[x](((|(Subscript[E, m;x])|))^2)dx can be written as:  
                         (∫_(-∞)^∞)ϵ[x](((|(Subscript[E, m;x])|))^2)dx =(((∫_(-∞)^∞)ϵ[x](((|(Subscript[E, m;x])|))^2)dx)/((∫_(-∞)^∞)(((|(Subscript[E, m;x])|))^2)dx))(∫_(-∞)^∞)(((|(Subscript[E, m;x])|))^2)dx 
                                          = ⟨ϵ⟩(∫_(-∞)^∞)(((|(Subscript[E, m;x])|))^2)dx     (3.19)
where  ⟨ϵ⟩ is the spatial average of ϵ. Thus:
                                                    (Subscript[P, m])=(c/(8 π  ))((⟨ϵ⟩)/(Subscript[n, eff;m]))(∫_(-∞)^∞)(((|(Subscript[E, m;x])|))^2)dx .               (3.20)
Notice the suttle difference in the power expression of TE and TM mode. The ratio ((⟨ϵ⟩)/(Subscript[n, eff;m]))can be considered as a modal average index.

3.1.3: In the lab: intensity plot

Intensity profile of guided modes

In the lab, if we look at a waveguide facet, what to we see? We see the intensity.

nindex={3.4,3.5}; a=2.5 ; λ=1.5 ;
segment={3.499,3.495,3.48,3.47,3.45,3.42,3.405}; neffTE={} ;
For[ i=1, i<Length[segment], i++,
	AppendTo[neffTE, 
	 x /. FindRoot[ResolvTE[nindex,a,λ,x]==0,{x,segment[[i+1]],segment[[i]]}]
                            ] ];
alldpTE={};
For[i=1,i<Length[neffTE]+1,i++,
	AppendTo[alldpTE,DensityPlot[neffTE[[i]]*Abs[FieldTE[nindex,a,λ,neffTE[[i]],x][[1]]]^2,{y,0,20},{x,-4,4}
	 ,PlotPoints->{2,202}, Mesh->False, ColorFunction->GrayLevel
	,ImageSize->{700,200}, AspectRatio-> 2/7.] ]
             ]

alldpTE

Intensity plot

We recall earlier that for some TM mode, if the index difference is large, the E field leakage is also quite large. We can see how it is in this example.
Note: this waveguide is NOT the same as what we show for TE above.

nindex={1,3.5}; a=0.5 ; λ=1.5 ;
Plot[ResolvTM[nindex,a,λ,x],{x,nindex[[1]],nindex[[2]]}]

segment2={3.49,3.4,3,2.5,1.5,1.01}; neffTM={} ;
For[ i=1, i<Length[segment2], i++,
	AppendTo[neffTM, 
	 x /. FindRoot[ResolvTM[nindex,a,λ,x]==0,{x,segment2[[i+1]],segment2[[i]]}]
                            ] ];
alldpTM={};
Eps[x_]:= If[Abs[x]<=a, nindex[[2]]^2,nindex[[1]]^2]; 

For[i=1,i≤ Length[neffTM],i++,
PrintTemporary[i];
AppendTo[alldpTM,DensityPlot[(Eps[x]/neffTM[[i]] )*Abs[FieldTM[nindex,a,λ,neffTM[[i]],x][[3]]]^2,{y,0,20},{x,-1.5,1.5} ,PlotPoints->{2,200}
, ColorFunction->GrayLevel
,ImageSize->{700,200}, AspectRatio-> 2/7.] ]
             ]

alldpTM

TM mode intensity: note the discontinuity

3.2 Confinement factor

An important concept (for all waveguide, not just slab): confinement factor: How much of the electric field is really in the core region? The fraction of the field in the core is called the confinement factor. (in complex structure, we need to define what "core" region is the region of our interest) 
                                (((∫_core^□)(((|E[x,y]|))^2)dxdy  )/((∫_all^□)(((|E[x,y]|))^2)dxdy))            (3.21)

3.3 Multi-mode propagation

3.3.1 Concept

If we launch an optical beam into a multimode waveguide, what do we get out at the other end? Again we apply the principle of linear superposition:
                                                     Φ[x]=(∑_(m=1)^N)(Subscript[a, m])(Subscript[F, m])[x] .                         (3.22)
We can represent the input signal as a finite-rank vector: 
Φ=(({{(Subscript[a, 1])}, {(Subscript[a, 2])}, {.}, {(Subscript[a, N])}})). The full expression of the electric field is: 
                    (Overscript[E, →])=(∑_(m=1)^N)(Subscript[a, m])(Subscript[(Overscript[E, →]), m])[x](e^(i ((Subscript[β, m])z-ω t)))  =(e^(-i ω t))(∑_(m=1)^N)(Subscript[a, m])(Subscript[(Overscript[E, →]), m])[x](e^(i (Subscript[β, m])z))    (3.23)
We notice one very important feature: we cannot separate the x and z dependent component of the net field (Overscript[E, →]) anymore. Thus: while the z-dependent and the x-dependent component can be separated for every eigenmode of the waveguide, they can not be so for any arbitrary wave.
The reason is modal dispersion as discussed above: each mode has a different propagation constant (Subscript[β, m]).

What does this entail? Remember how when you add waves together, not only amplitude but also phase is very important? Let's consider the input signal at z=0: 
                                        (Overscript[E, →])[x,z=0]=(e^(-i ω t))(∑_(m=1)^N)(Subscript[a, m])(Subscript[(Overscript[E, →]), m])[x] .              (3.24)
At z=L, we have:
                                       (Overscript[E, →])[x,z=L]=(e^(-i ω t))(∑_(m=1)^N)(Subscript[a, m])(e^(i (Subscript[β, m])L))(Subscript[(Overscript[E, →]), m])[x] 
                       =(e^(-i ω t))(∑_(m=1)^N)(Subscript[a, m])'(Subscript[(Overscript[E, →]), m])[x]        (3.25)
where                                                (Subscript[a, m])'=(Subscript[a, m])(e^(i (Subscript[β, m])L)).
Thus, the profile (Overscript[E, →])[x,z=L] can be very different from that at z=0.
This very important aspect is the foundation for waveguide circuit: modulate and control of light path through waveguide.

3.3.2 Example

Consider a slab waveguide with 2 modes:

nindex={1.4,1.5}; a=1.25 ; λ=1.5 ;
Plot[ResolvTE[nindex,a,λ,x],{x,nindex[[1]],nindex[[2]]}]

segment={1.49,1.46,1.42}; neff={} ;
For[ i=1, i<Length[segment], i++,
	AppendTo[neff, 
	 x /. FindRoot[ResolvTE[nindex,a,λ,x]==0,{x,segment[[i+1]],segment[[i]]}]
                            ] ]

neff

Let's look at the waveguide mode

Plot[{FieldTE[nindex,a,λ,neff[[1]],x][[1]]
		   ,FieldTE[nindex,a,λ,neff[[2]],x][[1]]},{x,-6,6}
	, PlotRange->{-1.,1}, PlotStyle->{{Thickness[0.01],Hue[0.35,1,0.5]},{Thickness[0.01],Hue[0.8,1,0.8]}}
]

At z=0, let's have a wave that consists of half and half of the low mode and the high mode:

w=0.5;
Plot[{w*FieldTE[nindex,a,λ,neff[[1]],x][[1]]
		   +(1-w)*FieldTE[nindex,a,λ,neff[[2]],x][[1]]},{x,-4,4}
	, PlotRange->{-1.,1}, PlotStyle->{{Thickness[0.01],Hue[0.35,1,0.5]},{Thickness[0.01],Hue[0.8,1,0.8]}}]

What does this wave look like as it propagates along the waveguide?
At z= 0:
(Overscript[E, →])[x,z=0]=(1/2)(e^(-i ω t))((Subscript[(Overscript[E, →]), 1])[x]+(Subscript[(Overscript[E, →]), 2])[x] )
(Overscript[E, →])[x,z=L]=(1/2)(e^(-i ω t))((e^(i (Subscript[β, 1])L)) (Subscript[(Overscript[E, →]), 1])[x]+(e^(i (Subscript[β, 2])L)) (Subscript[(Overscript[E, →]), 2])[x] )
                  =(1/2)(e^(-i ω t))(e^(i (Subscript[β, 1])L))( (Subscript[(Overscript[E, →]), 1])[x]+(e^(i ((Subscript[β, 2])-(Subscript[β, 1]))L)) (Subscript[(Overscript[E, →]), 2])[x] )
What happens at z=L when (e^(i ((Subscript[β, 2])-(Subscript[β, 1]))L)) = (e^(i π)) ? 
(Overscript[E, →])[x,z=L]=(1/2)(e^(-i ω t))(e^(i (Subscript[β, 1])L))( (Subscript[(Overscript[E, →]), 1])[x]- (Subscript[(Overscript[E, →]), 2])[x] )
Let's compare the field profile at z=0 and z=L, dropping the oscilling term (e^(-i ω t))(e^(i (Subscript[β, 1])L)) :

w=0.5;
Plot[{w*FieldTE[nindex,a,λ,neff[[1]],x][[1]]
		   +(1-w)*FieldTE[nindex,a,λ,neff[[2]],x][[1]]
	     , w*FieldTE[nindex,a,λ,neff[[1]],x][[1]]
		   -(1-w)*FieldTE[nindex,a,λ,neff[[2]],x][[1]]},{x,-4,4}
	, PlotRange->{-1.,1}, PlotStyle->{{Thickness[0.01],Hue[0.35,1,0.5]},{Thickness[0.01],Hue[0.8,1,0.8]}} ]

Remarkable: this tells us that the profile changes quite a bit as the wave propagates: it is not always a constant profile wave in a waveguide! What if it keeps on going?
(Overscript[E, →])[x,z=2L]=(1/2)(e^(-i ω t))(e^(i 2 (Subscript[β, 1])L))( (Subscript[(Overscript[E, →]), 1])[x]+(e^(i 2 ((Subscript[β, 2])-(Subscript[β, 1]))L)) (Subscript[(Overscript[E, →]), 2])[x] )
and (e^(i ((Subscript[β, 2])-(Subscript[β, 1]))L)) = (e^(i π)) : (Overscript[E, →])[x,z=2L]=(1/2)(e^(-i ω t))(e^(i 2 (Subscript[β, 1])L))( (Subscript[(Overscript[E, →]), 1])[x]+(Subscript[(Overscript[E, →]), 2])[x] ): we are back to square one again.
We can see all this by looking at its intensity. From the above expression of the Poynting vector: 
|⟨(Overscript[S, →])⟩.(OverHat[z])| = (c/(8 π μ (Subscript[k, 0])))|Re[ i (Subscript[E, y])(x,z) ((∂ (SuperStar[(Subscript[E, y])]))/(∂ z))]| 
(Subscript[E, y])(x,z)= w (Subscript[E, y;1])(x)(e^(i (Subscript[β, 1])z)) +(1-w) (Subscript[E, y;2])(x)(e^(i (Subscript[β, 2])z ))
((∂ (SuperStar[(Subscript[E, y])]))/(∂ z))=-i (Subscript[β, 1]) w (Subscript[E, y;1])(x)(e^(-i (Subscript[β, 1])z)) -i (Subscript[β, 2])(1-w) (Subscript[E, y;2])(x)(e^(-i (Subscript[β, 2])z ))
i (Subscript[E, y])(x,z) ((∂ (SuperStar[(Subscript[E, y])]))/(∂ z))=(w (Subscript[E, y;1])(x)(e^(i (Subscript[β, 1])z)) +(1-w) (Subscript[E, y;2])(x)(e^(i (Subscript[β, 2])z)) )
		                             ((Subscript[β, 1]) w (Subscript[E, y;1])(x)(e^(-i (Subscript[β, 1])z)) + (Subscript[β, 2])(1-w) (Subscript[E, y;2])(x)(e^(-i (Subscript[β, 2])z)) )

w=.;

Expand[(w (Subscript[Ef, y;1])(x)(E^(I (Subscript[β, 1])z)) +u (Subscript[Ef, y;2])(x)(E^(I (Subscript[β, 2])z)) )
		                             ((Subscript[β, 1]) w (Subscript[Ef, y;1])(x)(E^(-I (Subscript[β, 1])z)) + (Subscript[β, 2])u (Subscript[Ef, y;2])(x)(E^(-I (Subscript[β, 2])z)) )]

i (Subscript[E, y])(x,z) ((∂ (SuperStar[(Subscript[E, y])]))/(∂ z))=(w^2) (Subscript[β, 1]) (((Subscript[Ef, 1])(x))^2)+(u^2) (Subscript[β, 2]) (((Subscript[Ef, 2])(x))^2)
			       +((e^(i z( (Subscript[β, 2])- (Subscript[β, 1])))) (Subscript[β, 1]) +(e^(-i z( (Subscript[β, 2])- (Subscript[β, 1]))))  (Subscript[β, 2])) u w (Subscript[Ef, 2])(x) (Subscript[Ef, 1])(x)
There is a real term: (w^2) (Subscript[β, 1]) (((Subscript[Ef, 1])(x))^2)+(u^2) (Subscript[β, 2]) (((Subscript[Ef, 2])(x))^2) which is the individual sum of the intensity of each component. But more importantly, there is the interference term:
((e^(i z( (Subscript[β, 2])- (Subscript[β, 1])))) (Subscript[β, 1]) +(e^(i z( (Subscript[β, 2])- (Subscript[β, 1]))))  (Subscript[β, 2])) u w (Subscript[Ef, 2])(x) (Subscript[Ef, 1])(x). 
What is the effect of this interference term? Let's plot it out:

w=0.5; β=((2π)/λ)neff;
Plot3D[Re[(w*FieldTE[nindex,a,λ,neff[[1]],x][[1]](E^(I β[[1]]z))
		   +(1-w)*FieldTE[nindex,a,λ,neff[[2]],x][[1]](E^(I β[[2]]z)))*
		(w*β[[1]]*FieldTE[nindex,a,λ,neff[[1]],x][[1]](E^(-I β[[1]]z))
		   +(1-w)*β[[2]]*FieldTE[nindex,a,λ,neff[[2]],x][[1]](E^(-I β[[2]]z)))
				 ]
		,{z,0,30},{x,-2.5,2.5}, PlotPoints->{120,40}
	, PlotRange->All, Mesh->False, BoxRatios->{30,5,5},
ViewPoint->{0,0,5}]

Notice what happens: the wave seems to move from right side of the waveguide to left and back to right again and so on... : this is the behavior of a ray that bounces back and forth between the two boundary surface. Where is the shift maximum?

L=(π/(Abs[(β[[2]]-β[[1]])]))

What is a Δβ coupler? why is it called Δβ?

We have a waveguide with this property: cladding index= 1.4 ; thickness 2.5 um, the core index can be tuned with an external electric field from 1.43 to 1.5. A beam with 1.5 um wavelength is coupled to one side of the waveguide (mode1 + mode2), how does the external electric field can change its propagation property?

nindex={1.4,ncore}; a=1.25 ; λ=1.5 ;allp={};
For[i=1,i<12, i++, {ncore=1.45+(i-1)*(1.55-1.45)/10;
		AppendTo[allp,
Plot[ResolvTE[{1.4,ncore},a,λ,x],{x,1.4,ncore}, DisplayFunction->Identity]]
	              }]

Show[allp, DisplayFunction->$DisplayFunction]

nindex={1.4,1.45}; a=1.25 ; λ=1.5 ;
Plot[ResolvTE[nindex,a,λ,x],{x,nindex[[1]],nindex[[2]]}]

segment={1.44,1.43,1.401}; neff={} ;
For[ i=1, i<Length[segment], i++,
	AppendTo[neff, 
	 x /. FindRoot[ResolvTE[nindex,a,λ,x]==0,{x,segment[[i+1]],segment[[i]]}]
                            ] ]

neff

w=0.5; β=((2π)/λ)neff;
plot1=Plot3D[Re[(w*FieldTE[nindex,a,λ,neff[[1]],x][[1]](E^(I β[[1]]z))
		   +(1-w)*FieldTE[nindex,a,λ,neff[[2]],x][[1]](E^(I β[[2]]z)))*
		(w*β[[1]]*FieldTE[nindex,a,λ,neff[[1]],x][[1]](E^(-I β[[1]]z))
		   +(1-w)*β[[2]]*FieldTE[nindex,a,λ,neff[[2]],x][[1]](E^(-I β[[2]]z)))
				 ]
		,{z,0,30},{x,-2.5,2.5}, PlotPoints->{120,40}
	, PlotRange->All, Mesh->False, BoxRatios->{30,5,5},
ViewPoint->{0,0,5}]

nindex={1.4,1.55}; a=1.25 ; λ=1.5 ;
Plot[ResolvTE[nindex,a,λ,x],{x,nindex[[1]],nindex[[2]]}]

segment={1.54,1.5,1.46}; neff={} ;
For[ i=1, i<Length[segment], i++,
	AppendTo[neff, 
	 x /. FindRoot[ResolvTE[nindex,a,λ,x]==0,{x,segment[[i+1]],segment[[i]]}]
                            ] ]

neff

w=0.5; β=((2π)/λ)neff;
plot2=Plot3D[Re[(w*FieldTE[nindex,a,λ,neff[[1]],x][[1]](E^(I β[[1]]z))
		   +(1-w)*FieldTE[nindex,a,λ,neff[[2]],x][[1]](E^(I β[[2]]z)))*
		(w*β[[1]]*FieldTE[nindex,a,λ,neff[[1]],x][[1]](E^(-I β[[1]]z))
		   +(1-w)*β[[2]]*FieldTE[nindex,a,λ,neff[[2]],x][[1]](E^(-I β[[2]]z)))
				 ]
		,{z,0,30},{x,-2.5,2.5}, PlotPoints->{120,40}
	, PlotRange->All, Mesh->False, BoxRatios->{30,5,5},
ViewPoint->{0,0,5}]

3.4 Modal dispersion

The propagation constant β of each mode is a function of frequency ω: this is known as modal dispersion relation.

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))])
However, as the frequency change, the mode profile changes, and (Subscript[n, eff])[ω] is not constant for each mode. Recalling what we define as dispersion:

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.

Summary basic concepts

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.