Optical Fiber Intro - Part 2
Optical fiber

Loading package

CharG[nindex_List,λ_,ρ_,n_,l_]:=Module[{ϵ1,ϵ2,δ1, δ2,ξ1, ξ2,k0, J,K,Jp,Kp},
		{ϵ1,ϵ2}=nindex^2 ; k0=2*Pi/λ;
    	δ1=Sqrt[Abs[ϵ1 -n^2]]; δ2=Sqrt[Abs[ϵ2-n^2] ];
		ξ1=δ1*ρ *k0 ; ξ2=δ2*ρ* k0 ;
		J=BesselJ[l,ξ1] ; K=BesselK[l, ξ2];
		Jp=0.5*(BesselJ[l-1,ξ1] -BesselJ[l+1,ξ1]);
		Kp=-0.5*(BesselK[l-1,ξ2] + BesselK[l+1,ξ2]) ;
     (δ1 J Kp+ δ2 K Jp) (δ1  ϵ2 J Kp+δ2 ϵ1 K Jp)
	-(( (n^2))/( (k0^2) (ρ^2) ((δ1)^2) ((δ2)^2)))(l^2) (((ϵ1-ϵ2))^2) K^2 J^2   ]    ;

Coeff[nindex_List,λ_,ρ_,n_,l_]:=Module[{ϵ1,ϵ2,δ1, δ2,ξ1, ξ2,k0, J,K,Jp,Kp},
		{ϵ1,ϵ2}=nindex^2 ; k0=2*Pi/λ;
    	δ1=Sqrt[Abs[ϵ1 -n^2]]; δ2=Sqrt[Abs[ϵ2-n^2] ];
		ξ1=δ1*ρ *k0 ; ξ2=δ2*ρ* k0 ;
		J=BesselJ[l,ξ1] ; K=BesselK[l, ξ2];
		Jp=0.5*(BesselJ[l-1,ξ1] -BesselJ[l+1,ξ1]);
		Kp=-0.5*(BesselK[l-1,ξ2] + BesselK[l+1,ξ2]) ;
		{ 1, ((I l J K n (ϵ1-ϵ2))/(ρ k0 δ1 δ2 (δ2  K   Jp+δ1 J Kp))),
			-((J δ1)/(K δ2)) , 
                             -((I l (J^2)  n (ϵ1-ϵ2))/(ρ k0 ((δ2)^2) (δ2 K Jp+δ1 J Kp)))} 
 ]

FieldFiber[coor_List, nindex_List,λ_,rad_,n_,l_,coef_List]:=Module[{ρ, ϕ, ϵ1,ϵ2,q1, q2,ξ1, ξ2,k0, J,K,Jp,Kp, field},
		{ρ, ϕ}=coor ;
			{ϵ1,ϵ2}=nindex^2 ; k0=2*Pi/λ;
    	q1=Sqrt[Abs[ϵ1 -n^2]]; q2=Sqrt[Abs[ϵ2-n^2] ];
		ξ1=q1*ρ *k0 ; ξ2=q2*ρ* k0 ;
		J=BesselJ[l,ξ1] ; K=BesselK[l, ξ2];
		Jp=0.5*(BesselJ[l-1,ξ1] -BesselJ[l+1,ξ1]);
		Kp=-0.5*(BesselK[l-1,ξ2] + BesselK[l+1,ξ2]) ;
		field=If[ Abs[ρ]<=rad,
			 { J  q1
				, coef[[2]]*  J  q1
				, I n Jp -coef[[2]]*l *J/ξ1
				, -l *n J/ξ1 -coef[[2]] * I *Jp
				, l *ϵ1  J/ξ1+coef[[2]]*I n Jp
				, I ϵ1 Jp-coef[[2]]*l*J*n/ξ1}
			,{-coef[[3]]*K q2
				,-coef[[4]]*K q2
				,coef[[3]]*I n Kp-coef[[4]]*l K/ξ2
				,-coef[[3]]*l n K/ξ2-coef[[4]]*I *Kp
                  ,coef[[3]]*l K ϵ2/ξ2+coef[[4]] I n Kp
	,coef[[3]] I ϵ2 Kp-coef[[4]]l K n/ξ2}]  
;
		{I field * (E^(I ϕ l)),I field *(E^(-I ϕ l))}
		          ]  ;
(*  Note: the field output six components are: (({{(Subscript[E, z]) }, {(Subscript[H, z])}, {(Subscript[E, ρ])}, {(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))    *)

Overview and study guide

4. The fiber cylindrical geometry

4.1 Introduction

Step 1: Use cylindrical coordinate

For fiber, we use the cylindrical coordinate. Let z is the direction of propagation:
Let                                           E=(e^(i(β z - ω t)))F[ρ,ϕ]                             (4.1)
The wave equation becomes:
                  (1/ρ)((∂ F)/(∂ ρ)) + (((∂^2) F)/(∂ (ρ^2))) + (1/(ρ^2))(((∂^2) F)/(∂ (ϕ^2)))  +(((ω^2)/(c^2))ϵ(ρ,ϕ)  -  (β^2) )F= 0          (4.2)
Let's assume further angular geometry: 
                            E=(e^(i(β z - ω t)))F[ρ,ϕ]   ->   (e^(i(β z - ω t))) (e^(i l ϕ)) F[ρ]           (4.3)
Substitute in Eq. (4.2)
                             (1/ρ)((∂ F)/(∂ ρ)) + (((∂^2) F)/(∂ (ρ^2)))  +(((Subscript[k, 0])^2) ϵ(ρ)  -((l^2)/(ρ^2)) - (β^2) )F= 0            (4.4)
What does the vector field look like? E, H  = (({{(Subscript[(E or H), ρ])}, {(Subscript[(E or H), ϕ])}, {(Subscript[(E or H), z])}}));

Similar to the planar waveguide, we solve the wave Eq. in each medium and make sure to match the boundary conditions.
To solve Eq. (4.4) in each medium, we need to know  ϵ(ρ). The profile of ϵ(ρ) can be non-trivial. If each medium is just a uniform dielectric, then ϵ(ρ)  is just a constant. In some case, people make graded-index fiber in which ϵ(ρ) varies as a function of ρ, such as linear and quadratic.
Regardless, the key is still the boundary condition as we can see.

4.2 Boundary condition for fiber geometry

Step 2: Set up boundary condition

Just like the planar waveguide, we need to connect D and B (E and H) across an interface. First, we need to find the relation between these fields from the Maxwell's Eqs. in the system of coordinate of choice.
The first is:
                                                -(1/c)μ((∂ (Overscript[H, ->]))/(∂ t))= ∇×(Overscript[E, ->])                                     (4.5)

Curl[{Eρ[ρ,ϕ,z],Eϕ[ρ,ϕ,z],Ez[ρ,ϕ,z]}, Cylindrical[ρ,ϕ,z] ]

Eq. (4.5) , -(1/c)μ((∂ (Overscript[H, ->]))/(∂ t))=((i ω)/c)μ (Overscript[H, ->]) = ∇×(Overscript[E, ->])      gives:

                ((i ω)/c)μ(({{(Subscript[H, ρ])}, {(Subscript[H, ϕ])}, {(Subscript[H, z])}})) = (({{(1/ρ)(∂_ϕ) (Subscript[E, z])[ρ,ϕ,z] - (∂_z) (Subscript[E, ϕ])[ρ,ϕ,z]}, {(∂_z) (Subscript[E, ρ])[ρ,ϕ,z]-(∂_ρ) (Subscript[E, z])[ρ,ϕ,z]}, {(∂_ρ) (Subscript[E, ϕ])[ρ,ϕ,z]+(1/ρ) (Subscript[E, ϕ])[ρ,ϕ,z]-(1/ρ)(∂_ϕ) (Subscript[E, ρ])[ρ,ϕ,z]}}))         
			
                           =  (({{((i l)/ρ)(Subscript[E, z])[ρ] - i β (Subscript[E, ϕ])[ρ]}, {i β (Subscript[E, ρ])[ρ]-(∂_ρ) (Subscript[E, z])[ρ]}, {(∂_ρ) (Subscript[E, ϕ])[ρ]+(1/ρ) (Subscript[E, ϕ])[ρ]-((i l)/ρ) (Subscript[E, ρ])[ρ]}}))                                        (4.6)
and the second Maxwell's Eq.
                                              ((-i ω)/c)ϵ  (Overscript[E, ->])  = ∇×(Overscript[H, ->])                                             (4.7)
gives:
				
              ((-i ω)/c)ϵ(({{(Subscript[E, ρ])}, {(Subscript[E, ϕ])}, {(Subscript[E, z])}})) =  (({{((i l)/ρ)(Subscript[H, z])[ρ] - i β (Subscript[H, ϕ])[ρ]}, {i β (Subscript[H, ρ])[ρ]-(∂_ρ) (Subscript[H, z])[ρ]}, {(∂_ρ) (Subscript[H, ϕ])[ρ]+(1/ρ) (Subscript[H, ϕ])[ρ]-((i l)/ρ)(Subscript[H, ρ])[ρ]}}))                           (4.8)

4.3 Reduction to simpler forms

We have six components. Not all boundary conditions for them are independent from each other, because of their relations as described by Maxwell's eqs. To reduce them, we first start out with tangentials: (Subscript[E, z]) and (Subscript[H, z]): they themselves are continuous across a boundary. We will find how other are related to them. We have a solution with one unknown for each: (Subscript[E, z])[ρ], (Subscript[H, z])[ρ]. 
In the follow, the blue terms are those that we consider as "known" with which we try to express the others.
        ((-i ω)/c)ϵ(({{(Subscript[E, ρ])}, {(Subscript[E, ϕ])}})) =  (({{((i l)/ρ)(Subscript[H, z])[ρ] - i β (Subscript[H, ϕ])[ρ]}, {i β (Subscript[H, ρ])[ρ]-(∂_ρ) (Subscript[H, z])[ρ]}}))= (({{- i β (Subscript[H, ϕ])[ρ]}, {i β (Subscript[H, ρ])[ρ]}}))+ (({{((i l)/ρ)(Subscript[H, z])[ρ] }, {-(∂_ρ) (Subscript[H, z])[ρ]}}))     (4.9)
which gives: 
                (({{(Subscript[E, ρ])}, {(Subscript[E, ϕ])}})) = (( i β)/(-i ϵ  (Subscript[k, 0]))) (({{- (Subscript[H, ϕ])[ρ]}, {(Subscript[H, ρ])[ρ]}}))+ ((i  )/((Subscript[k, 0]) ϵ))(({{((i l)/ρ)(Subscript[H, z])[ρ] }, {-(∂_ρ) (Subscript[H, z])[ρ]}}))
              = -(β/( ϵ  (Subscript[k, 0]))) (({{- (Subscript[H, ϕ])[ρ]}, {(Subscript[H, ρ])[ρ]}}))+((i  )/((Subscript[k, 0]) ϵ))(({{((i l)/ρ)(Subscript[H, z])[ρ] }, {-(∂_ρ) (Subscript[H, z])[ρ]}}))                   (4.10)
and from Eq. (4.6):
                         ((i ω)/c)μ(({{(Subscript[H, ρ])}, {(Subscript[H, ϕ])}, {(Subscript[H, z])}})) =  (({{((i l)/ρ)(Subscript[E, z])[ρ] - i β (Subscript[E, ϕ])[ρ]}, {i β (Subscript[E, ρ])[ρ]-(∂_ρ) (Subscript[E, z])[ρ]}, {(∂_ρ) (Subscript[E, ϕ])[ρ]+(1/ρ) (Subscript[E, ϕ])[ρ]-((i l)/ρ) (Subscript[E, ρ])[ρ]}})),
we obtain:
                         ((i ω)/c)μ(({{(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))   = (({{- i β (Subscript[E, ϕ])[ρ]}, {i β (Subscript[E, ρ])[ρ]}}))+ (({{((i l)/ρ)(Subscript[E, z])[ρ] }, {-(∂_ρ) (Subscript[E, z])[ρ]}}))             (4.11)
Which gives: 
                    (({{(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))   = ((i β )/( i  (Subscript[k, 0]) μ))(({{-  (Subscript[E, ϕ])[ρ]}, {(Subscript[E, ρ])[ρ]}}))+ (( 1)/( i (Subscript[k, 0])  μ))(({{((i l)/ρ)(Subscript[E, z])[ρ] }, {-(∂_ρ) (Subscript[E, z])[ρ]}}))             (4.12)
Thus, starting from the first equation:  
       (({{(Subscript[E, ρ])}, {(Subscript[E, ϕ])}})) =  -(( β)/(ϵ  (Subscript[k, 0]))) (({{0, -1}, {1, 0}}))(({{(Subscript[H, ρ])}, {(Subscript[H, ϕ])}})) +((i  )/((Subscript[k, 0])  ϵ)) (({{((i l)/ρ)(Subscript[H, z])[ρ] }, {-(∂_ρ) (Subscript[H, z])[ρ]}})) 
          =   -(( β)/(ϵ  (Subscript[k, 0])))(({{0, -1}, {1, 0}})) {((β )/(  (Subscript[k, 0]) μ))(({{0, -1}, {1, 0}}))(({{(Subscript[E, ρ])}, {(Subscript[E, ϕ])}})) + (( 1)/( i (Subscript[k, 0])  μ))(({{((i l)/ρ)(Subscript[E, z])[ρ] }, {-(∂_ρ) (Subscript[E, z])[ρ]}}))} 
		                                              +((i  )/((Subscript[k, 0])  ϵ))(({{((i l)/ρ)(Subscript[H, z])[ρ] }, {-(∂_ρ) (Subscript[H, z])[ρ]}}))        
            =  ((β^2)/(ϵ μ  (Subscript[k, 0]^2)))(({{(Subscript[E, ρ])}, {(Subscript[E, ϕ])}})) + ((i β )/( (Subscript[k, 0]^2)  ϵ μ))(({{(∂_ρ) (Subscript[E, z])[ρ] }, {((i l)/ρ)(Subscript[E, z])[ρ]}})) +((i  )/((Subscript[k, 0])  ϵ))(({{((i l)/ρ)(Subscript[H, z])[ρ] }, {-(∂_ρ) (Subscript[H, z])[ρ]}}))           (4.13)
Or:
           (1-((β^2)/(ϵ μ  (Subscript[k, 0]^2)))) (({{(Subscript[E, ρ])}, {(Subscript[E, ϕ])}})) =  ((i β )/(  ((Subscript[k, 0])^2) ϵ μ))(({{(∂_ρ) (Subscript[E, z])[ρ] }, {((i l)/ρ)(Subscript[E, z])[ρ]}})) +((i )/((Subscript[k, 0]) ϵ))(({{((i l)/ρ)(Subscript[H, z])[ρ] }, {-(∂_ρ) (Subscript[H, z])[ρ]}}))  ;    (4.14)
 and:
         (({{(Subscript[E, ρ])}, {(Subscript[E, ϕ])}})) =((ϵ μ ((Subscript[k, 0])^2))/((ϵ μ ((Subscript[k, 0])^2)- (β^2) )))((i β )/(  ((Subscript[k, 0])^2) ϵ μ))(({{(∂_ρ) (Subscript[E, z])[ρ] }, {((i l)/ρ)(Subscript[E, z])[ρ]}})) +((ϵ μ ((Subscript[k, 0])^2))/((ϵ μ ((Subscript[k, 0])^2)- (β^2) )))((i )/((Subscript[k, 0]) ϵ))(({{((i l)/ρ)(Subscript[H, z])[ρ] }, {-(∂_ρ) (Subscript[H, z])[ρ]}}))      (4.15)
         (({{(Subscript[E, ρ])}, {(Subscript[E, ϕ])}})) =((i β)/((ϵ μ ((Subscript[k, 0])^2)- (β^2) )))(({{(∂_ρ) (Subscript[E, z])[ρ] }, {((i l)/ρ)(Subscript[E, z])[ρ]}})) +i(((Subscript[k, 0]) μ)/((ϵ μ ((Subscript[k, 0])^2)- (β^2) ))) (({{((i l)/ρ)(Subscript[H, z])[ρ] }, {-(∂_ρ) (Subscript[H, z])[ρ]}}))    (4.16)

Another way to write it: 
        (({{(Subscript[E, ρ])}, {(Subscript[E, ϕ])}})) =(1/((ϵ μ ((Subscript[k, 0])^2)- (β^2) )))(({{0, i β, -(Subscript[k, 0]) μ(l/ρ), 0}, {-((β l)/ρ), 0, 0, -i (Subscript[k, 0]) μ}})). (({{(Subscript[E, z])[ρ]}, {(∂_ρ) (Subscript[E, z])[ρ]}, {(Subscript[H, z])[ρ]}, {(∂_ρ) (Subscript[H, z])[ρ]}}))             (4.17)
Do the same thing for H: From above:
               (({{(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))   = ((β )/(  (Subscript[k, 0]) μ))(({{-  (Subscript[E, ϕ])[ρ]}, {(Subscript[E, ρ])[ρ]}}))+ (( 1)/( i (Subscript[k, 0])  μ))(({{((i l)/ρ)(Subscript[E, z])[ρ] }, {-(∂_ρ) (Subscript[E, z])[ρ]}})) 
                           =  ((β )/(  (Subscript[k, 0]) μ))(({{0, -1}, {1, 0}}))(({{(Subscript[E, ρ])[ρ]  }, {(Subscript[E, ϕ])[ρ]}}))+ (( 1)/( i (Subscript[k, 0])  μ))(({{((i l)/ρ)(Subscript[E, z])[ρ] }, {-(∂_ρ) (Subscript[E, z])[ρ]}}))               (4.18)

         (({{(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))   = ((β )/((Subscript[k, 0]) μ))(({{0, -1}, {1, 0}})). (((i β)/((ϵ μ ((Subscript[k, 0])^2)- (β^2) )))(({{(∂_ρ) (Subscript[E, z])[ρ] }, {((i l)/ρ)(Subscript[E, z])[ρ]}})) +i(((Subscript[k, 0]) μ)/((ϵ μ ((Subscript[k, 0])^2)- (β^2) ))) (({{((i l)/ρ)(Subscript[H, z])[ρ] }, {-(∂_ρ) (Subscript[H, z])[ρ]}})) ) 
                                                      + (( 1)/( i (Subscript[k, 0])  μ))(({{((i l)/ρ)(Subscript[E, z])[ρ] }, {-(∂_ρ) (Subscript[E, z])[ρ]}}))                               (4.19)
 After simplification:

Simplify[((β )/((Subscript[k, 0]) μ))((I β)/((ϵ μ ((Subscript[k, 0])^2)- (β^2) )))(({{0, -1}, {1, 0}})). (({{dEz }, {((I l)/ρ)Ez}}))  + (( 1)/( I (Subscript[k, 0])  μ))(({{((I l)/ρ)Ez }, {-dEz }})) ]
((β )/((Subscript[k, 0]) μ))I(((Subscript[k, 0]) μ)/((ϵ μ ((Subscript[k, 0])^2)- (β^2) )))(({{0, -1}, {1, 0}})).(({{((I l)/ρ)Hz }, {-dHz}}))

Simplify[(({{((Ez l ϵ (Subscript[k, 0]))/(ϵ μ ρ (Subscript[k, 0]^2)-(β^2) ρ))}, {((I dEz ϵ (Subscript[k, 0]))/(ϵ μ (Subscript[k, 0]^2)-(β^2)))}}))*((ϵ μ (Subscript[k, 0]^2)-(β^2))/(-I ϵ (Subscript[k, 0])))]
Simplify[(({{((I dHz β)/(ϵ μ (Subscript[k, 0]^2)-(β^2)))}, {-((Hz l β)/(ρ (ϵ μ (Subscript[k, 0]^2)-(β^2))))}}))*((ϵ μ (Subscript[k, 0]^2)-(β^2))/(I  β))]

(({{(Subscript[H, ρ])}, {(Subscript[H, ϕ])}})) =((-i ϵ (Subscript[k, 0]))/(ϵ μ (Subscript[k, 0]^2)-(β^2)))(({{((i l)/ρ)(Subscript[E, z])[ρ] }, {-(∂_ρ) (Subscript[E, z])[ρ]}}))     + ((i  β)/(ϵ μ (Subscript[k, 0]^2)-(β^2)))(({{(∂_ρ) (Subscript[H, z])[ρ]}, {((i l)/ρ)(Subscript[H, z])[ρ]}}))        (4.20)
or:
              (({{(Subscript[H, ρ])}, {(Subscript[H, ϕ])}})) =(1/(ϵ μ (Subscript[k, 0]^2)-(β^2)))(({{ϵ (Subscript[k, 0])(( l)/ρ), 0, 0, i  β}, {0, i ϵ (Subscript[k, 0]), -((β l )/ρ), 0}})). (({{(Subscript[E, z])[ρ]}, {(∂_ρ) (Subscript[E, z])[ρ]}, {(Subscript[H, z])[ρ]}, {(∂_ρ) (Subscript[H, z])[ρ]}}))         (4.21)
So, we have solved for all components of E and H in terms of (Subscript[E, z]) and (Subscript[H, z]).

4.4 The use of S matrix

Step 3: Method to solve boundary condition problem and obtain characteristic equation

Can S-matrix formalism be used in this case for boundary condition? Yes.
Let:                                     (({{(Subscript[E, z ])}, {(∂_ρ) (Subscript[E, z])}}))= (({{(Subscript[ψ, A]), (Subscript[ψ, B])}, {(Subscript[ψ, A])', (Subscript[ψ, B])'}})).(({{(Subscript[c, A])}, {(Subscript[c, B])}})) ;                     (4.22a)
                                           (({{(Subscript[H, z ])}, {(∂_ρ) (Subscript[H, z])}}))= (({{(Subscript[ψ, A]), (Subscript[ψ, B])}, {(Subscript[ψ, A])', (Subscript[ψ, B])'}})). (({{(Subscript[d, A])}, {(Subscript[d, B])}})) .                     (4.22b)
Since the boundary condition involves mixes of E and H, we can combine them into a 4x4 matrix: 
                         (({{(Subscript[E, z ])}, {(∂_ρ) (Subscript[E, z])}, {(Subscript[H, z ])}, {(∂_ρ) (Subscript[H, z])}}))=(({{(Subscript[ψ, A]), (Subscript[ψ, B]), 0, 0}, {(Subscript[ψ, A])', (Subscript[ψ, B])', 0, 0}, {0, 0, (Subscript[ψ, A]), (Subscript[ψ, B])}, {0, 0, (Subscript[ψ, A])', (Subscript[ψ, B])'}})) (({{(Subscript[c, A])}, {(Subscript[c, B])}, {(Subscript[d, A])}, {(Subscript[d, B])}}))                     (4.23)
We write other boundary condition: From above: 
         (({{(Subscript[E, ρ])}, {(Subscript[E, ϕ])}})) =(1/((ϵ μ ((Subscript[k, 0])^2)- (β^2) )))(({{0, i β, -(Subscript[k, 0]) μ(l/ρ), 0}, {-((β l)/ρ), 0, 0, -i (Subscript[k, 0]) μ}})). (({{(Subscript[E, z])[ρ]}, {(∂_ρ) (Subscript[E, z])[ρ]}, {(Subscript[H, z])[ρ]}, {(∂_ρ) (Subscript[H, z])[ρ]}}))            (4.24)
         (({{ϵ (Subscript[E, ρ])}, {(Subscript[E, ϕ])}})) =(1/((ϵ μ ((Subscript[k, 0])^2)- (β^2) )))( ({{0, i β ϵ, -(Subscript[k, 0]) μ ϵ(l/ρ) , 0}, {-((β l)/ρ), 0, 0, -i (Subscript[k, 0]) μ }}))(({{(Subscript[E, z]) }, {(∂_ρ)(Subscript[E, z])}, {(Subscript[H, z]) }, {(∂_ρ)(Subscript[H, z])}}))       (4.25)
Also from above: 
             (({{(Subscript[H, ρ])}, {(Subscript[H, ϕ])}})) =(1/((ϵ μ ((Subscript[k, 0])^2)- (β^2) )))(({{ϵ  (Subscript[k, 0])(( l)/ρ), 0, 0 , i β}, {0, i ϵ  (Subscript[k, 0]), -(( l)/ρ) β, 0 }}))(({{(Subscript[E, z]) }, {(∂_ρ)(Subscript[E, z])}, {(Subscript[H, z]) }, {(∂_ρ)(Subscript[H, z])}}))             (4.26)
Combine them:
         (({{ϵ (Subscript[E, ρ])}, {(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))=(1/((ϵ μ ((Subscript[k, 0])^2)- (β^2) )))(({{0, i β ϵ, -(Subscript[k, 0]) μ ϵ(l/ρ) , 0}, {-(( β l)/ρ), 0, 0, -i (Subscript[k, 0]) μ }, {ϵ  (Subscript[k, 0])(( l)/ρ), 0, 0 , i β}, {0, i ϵ  (Subscript[k, 0]), -(( l)/ρ) β, 0 }})) (({{(Subscript[E, z]) }, {(∂_ρ)(Subscript[E, z])}, {(Subscript[H, z]) }, {(∂_ρ)(Subscript[H, z])}}))             (4.27)
With 2 more conditions for Ez and Hz: 
           (({{1, 0, 0, 0}, {0, 0, 1, 0}, {0, ((i β ϵ)/Q), ((-(Subscript[k, 0]) μ ϵ)/Q)(l/ρ) , 0}, {-(( β l)/(Q ρ)), 0, 0, -((i (Subscript[k, 0]) μ )/Q)}, {ϵ  (Subscript[k, 0])(( l)/(Q ρ)), 0, 0 , ((i β)/Q)}, {0, ((i ϵ  (Subscript[k, 0]))/Q), -(( l)/(ρ Q)) β, 0 }})) . (({{(Subscript[E, z]) }, {(∂_ρ)(Subscript[E, z])}, {(Subscript[H, z]) }, {(∂_ρ)(Subscript[H, z])}}))=(({{(Subscript[E, z]) }, {(Subscript[H, z])}, {ϵ (Subscript[E, ρ])}, {(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))  ;           (4.28)
where                                                          (1/((ϵ μ ((Subscript[k, 0])^2)- (β^2) )))  = (1/Q)                      (4.29)
       (({{(Subscript[E, z]) }, {(Subscript[H, z])}, {ϵ (Subscript[E, ρ])}, {-(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))=(({{1, 0, 0, 0}, {0, 0, 1, 0}, {0, ((i β ϵ)/Q), ((-(Subscript[k, 0]) μ ϵ)/Q)(l/ρ) , 0}, {(( β l)/(Q ρ)), 0, 0, ((i (Subscript[k, 0]) μ )/Q)}, {ϵ  (Subscript[k, 0])(( l)/(Q ρ)), 0, 0 , ((i β)/Q)}, {0, ((i ϵ  (Subscript[k, 0]))/Q), -(( l)/(ρ Q)) β, 0 }})) (({{(Subscript[ψ, A]), (Subscript[ψ, B]), 0, 0}, {(Subscript[ψ, A])', (Subscript[ψ, B])', 0, 0}, {0, 0, (Subscript[ψ, A]), (Subscript[ψ, B])}, {0, 0, (Subscript[ψ, A])', (Subscript[ψ, B])'}})) (({{(Subscript[c, A])}, {(Subscript[c, B])}, {(Subscript[d, A])}, {(Subscript[d, B])}}))         (4.30)

We have only 4 unknowns and 6 equations! Surely the eqs are not all independent of each other. We can check out:

Simplify[ Det[(({{0, ((I β ϵ)/Q), ((-(Subscript[k, 0]) μ ϵ)/Q)(l/ρ) , 0}, {(( β l)/(Q ρ)), 0, 0, ((I (Subscript[k, 0]) μ )/Q)}, {ϵ  (Subscript[k, 0])(( l)/(Q ρ)), 0, 0 , ((I β)/Q)}, {0, ((I ϵ  (Subscript[k, 0]))/Q), -(( l)/(ρ Q)) β, 0 }}))]] /. ϵ μ (Subscript[k, 0]^2)- (β^2)  ->Q  /.  (β^2)-ϵ μ (Subscript[k, 0]^2)->-Q

This set of 4 is linearly independent for l not =0

Simplify[ Det[(({{1, 0, 0 , 0}, {0, 0, 1, 0}, {ϵ  (Subscript[k, 0])(( l)/(Q ρ)), 0, 0 , ((I β)/Q)}, {0, ((I ϵ  (Subscript[k, 0]))/Q), -(( l)/(Q ρ)) β, 0 }}))]]

This set of 4 is also linearly independent. The sets below are not:

Det[(({{1, 0, 0 , 0}, {0, 0, 1, 0}, {0, ((I β ϵ)/Q), -(Subscript[k, 0]) μ ϵ(l/(Q ρ)) , 0}, {0, ((I ϵ  (Subscript[k, 0]))/Q), -(( l)/(Q ρ)) β, 0 }}))]
 Det[(({{1, 0, 0 , 0}, {0, 0, 1, 0}, {(( β l)/(Q ρ)), 0, 0, ((I (Subscript[k, 0]) μ)/Q) }, {ϵ  (Subscript[k, 0])(( l)/(Q ρ)), 0, 0, ((I β)/Q)}}))/Q]

These sets are NOT lin. indep. We can choose a set like this 
(({{(Subscript[E, z]) }, {(Subscript[H, z])}, {ϵ (Subscript[E, ρ])}, {-(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))=(({{1, 0, 0, 0}, {0, 0, 1, 0}, {0, ((i β ϵ)/Q), ((-(Subscript[k, 0]) μ ϵ)/Q)(l/ρ) , 0}, {(( β l)/(Q ρ)), 0, 0, ((i (Subscript[k, 0]) μ )/Q)}, {ϵ  (Subscript[k, 0])(( l)/(Q ρ)), 0, 0 , ((i β)/Q)}, {0, ((i ϵ  (Subscript[k, 0]))/Q), -(( l)/(ρ Q)) β, 0 }})) (({{(Subscript[ψ, A]), (Subscript[ψ, B]), 0, 0}, {(Subscript[ψ, A])', (Subscript[ψ, B])', 0, 0}, {0, 0, (Subscript[ψ, A]), (Subscript[ψ, B])}, {0, 0, (Subscript[ψ, A])', (Subscript[ψ, B])'}})) (({{(Subscript[c, A])}, {(Subscript[c, B])}, {(Subscript[d, A])}, {(Subscript[d, B])}}))          (4.31)
(({{(Subscript[E, z]) }, {(Subscript[H, ϕ])}, {(Subscript[H, z])}, {-(Subscript[E, ϕ])}}))=(({{1, 0, 0 , 0}, {0, ((i ϵ  (Subscript[k, 0]))/Q), -(( l)/(Q ρ)) β, 0 }, {0, 0, 1, 0 }, {(( β l)/(Q ρ)), 0, 0, ((i (Subscript[k, 0]) μ)/Q)}}))(({{(Subscript[ψ, A]), (Subscript[ψ, B]), 0, 0}, {(Subscript[ψ, A])', (Subscript[ψ, B])', 0, 0}, {0, 0, (Subscript[ψ, A]), (Subscript[ψ, B])}, {0, 0, (Subscript[ψ, A])', (Subscript[ψ, B])'}})) (({{(Subscript[c, A])}, {(Subscript[c, B])}, {(Subscript[d, A])}, {(Subscript[d, B])}}))           (4.32)
Since all these components are tangential, they are continuous across an interface.

Det[(({{1, 0, 0 , 0}, {0, ((I ϵ  (Subscript[k, 0]))/Q), -(( l)/(Q ρ)) β, 0 }, {0, 0, 1, 0 }, {(( β l)/(Q ρ)), 0, 0, ((I (Subscript[k, 0]) μ)/Q)}}))]

Inverse[  (({{1, 0, 0 , 0}, {0, ((I ϵ  (Subscript[k, 0]))/Q), -(( l)/(Q ρ)) β, 0 }, {0, 0, 1, 0 }, {(( β l)/(Q ρ)), 0, 0, ((I (Subscript[k, 0]) μ)/Q)}})) ]

The interface S matrix here is:
S= (((({{(Subscript[ψ, A1]), (Subscript[ψ, B1]), 0, 0}, {(Subscript[ψ, A1])', (Subscript[ψ, B1])', 0, 0}, {0, 0, (Subscript[ψ, A1]), (Subscript[ψ, B1])}, {0, 0, (Subscript[ψ, A1])', (Subscript[ψ, B1])'}})))^(-1))  (({{1, 0, 0, 0}, {0, -((i (Subscript[Q, 1]))/((Subscript[ϵ, 1]) (Subscript[k, 0]))), -((i l β)/((Subscript[ϵ, 1]) ρ (Subscript[k, 0]))), 0}, {0, 0, 1, 0}, {((i l β)/(μ ρ (Subscript[k, 0]))), 0, 0, -((i (Subscript[Q, 1]))/(μ (Subscript[k, 0])))}}))          (({{1, 0, 0 , 0}, {0, ((i (Subscript[ϵ, 2])  (Subscript[k, 0]))/(Subscript[Q, 2])), -(( l)/((Subscript[Q, 2]) ρ)) β, 0 }, {0, 0, 1, 0 }, {(( β l)/((Subscript[Q, 2]) ρ)), 0, 0, ((i (Subscript[k, 0]) μ)/(Subscript[Q, 2]))}}))(({{(Subscript[ψ, A2]), (Subscript[ψ, B2]), 0, 0}, {(Subscript[ψ, A2])', (Subscript[ψ, B2])', 0, 0}, {0, 0, (Subscript[ψ, A2]), (Subscript[ψ, B2])}, {0, 0, (Subscript[ψ, A2])', (Subscript[ψ, B2])'}}))                (4.33)
which connect the four fields (Subscript[E, z]),(Subscript[H, ϕ]),(Subscript[H, z]),(Subscript[E, ϕ]) from one side to the other.

Simplify[ (({{1, 0, 0, 0}, {0, -((I (Subscript[Q, 1]))/((Subscript[ϵ, 1]) (Subscript[k, 0]))), -((I l β)/((Subscript[ϵ, 1]) ρ (Subscript[k, 0]))), 0}, {0, 0, 1, 0}, {((I l β)/(μ ρ (Subscript[k, 0]))), 0, 0, -((I (Subscript[Q, 1]))/(μ (Subscript[k, 0])))}})).(({{1, 0, 0 , 0}, {0, ((I (Subscript[ϵ, 2])  (Subscript[k, 0]))/(Subscript[Q, 2])), -(( l)/((Subscript[Q, 2]) ρ)) β, 0 }, {0, 0, 1, 0 }, {(( β l)/((Subscript[Q, 2]) ρ)), 0, 0, ((I (Subscript[k, 0]) μ)/(Subscript[Q, 2]))}}))]

This is the transfer matrix we can use.

5. Step-index fiber

5.1 Introduction

There are several designs for fiber. Functionally, there are only 2 kinds: single mode and multi-mode. But there are many ways to design fibers for either function. Different designs are made to attain certain desirable properties for each type for various applications.

Here we study the step-index fiber, which is the most common case. This type of fiber can have many layers, but each layer is uniform and the index of refraction is constant.
Below we will study the simplest of all, a structure that consists of only a core and a cladding that is thick enough to be considered as " infinite". Although some real fibers have slightly more elaborated index-of-refraction profile, the case of this simple step-index fiber is very useful to illustrate the more general features common to all fibers.

5.2 Discussion

To solve for the optical modes, we need to do two things:
         1- Solve the wave equation (Helmholtz) for each region
         2- Apply the boundary conditions at various interfaces
To solve for the waveguide mode, we need solutions of the eq.: 
             (1/ρ)((∂ F)/(∂ ρ)) + (((∂^2) F)/(∂ (ρ^2)))  +(((Subscript[k, 0])^2) ϵ(ρ)  -((l^2)/(ρ^2)) - (β^2) )F= 0
          ρ((∂ F)/(∂ ρ)) +(ρ^2) (((∂^2) F)/(∂ (ρ^2)))  +((((Subscript[k, 0])^2) ϵ(ρ) - (β^2))(ρ^2) -(l^2)  )F= 0
For the step-index fiber, ϵ(ρ) is constant of ρ for each layer of the fiber. 
Let  ϵ(ρ) be constant, and define (q^2) = ((Subscript[k, 0])^2) ϵ -  (β^(2 ))
                                  ρ((∂ F)/(∂ ρ)) +(ρ^2) (((∂^2) F)/(∂ (ρ^2)))  +((q^2)(ρ^2) -(l^2)  )F= 0          (5.1)

DSolve[ ρ F'[ρ] + (ρ^2) F''[ρ] +(((q )^2)(ρ^2)  -(l^2)  )F[ρ ]== 0 ,
	                             F[ρ ], ρ ]

5.2.1 Bessel function solution

Let q=(Sqrt[ϵ ((Subscript[k, 0])^2) -(β^(2 ))])for core region (high index), ϵ ((Subscript[k, 0])^2) -(β^2)> 0 and q is real > 0
Define  ξ = ρ q .  Then ((∂ )/(∂ ρ)) ->q ((∂ )/(∂ ξ  )); ((∂^2)/(∂ (ρ^2))) ->(q^2)((∂^2)/(∂ (ξ^2)))
 ρ F'[ρ] + (ρ^2) F''[ρ] +(((q )^2)(ρ^2)  -(l^2)  )F[ρ ]== 0  becomes:

ρ F'[ρ] + (ρ^2) F''[ρ] +(((q )^2)(ρ^2)  -(l^2)  )F[ρ ] /. ρ-> ξ /. q->1

DSolve[Fa''(ξ) (ξ^2)+Fa'(ξ) ξ+((ξ^2)-(l^2)) Fa(ξ)==0,Fa[ξ],ξ]

The linearly indepent solutions are Bessel functions J and Y

In the cladding region where ϵ ((Subscript[k, 0])^2) -(β^2)< 0; define  q=(Sqrt[(β^2)- ϵ ((Subscript[k, 0])^2)]). 
                   ρ F'[ρ] + (ρ^2) F''[ρ] -( ((q )^2)(ρ^2)  +(l^2)  )F[ρ ]== 0                  
 or                      ξ F'[ξ] + (ξ^2) F''[ξ] -( (ξ^2)  +(l^2)  )F[ξ ]== 0
The solutions are Bessel I and K

DSolve[ξ F'[ξ] + (ξ^2) F''[ξ] -( (ξ^2)+(l^2)  )F[ξ ]== 0,
	     F[ξ],ξ]

Plot[ Release[Table[BesselJ[n, x] ,{n,0,5}]],{x,0,10}
, PlotStyle->{{Hue[0.2]},{Hue[0.3]},{Hue[0.4]},{Hue[0.5]},{Hue[0.6]},{Hue[0.7]},{Hue[0.8]},{Hue[0.9]}} 
,Frame-> True, GridLines->Automatic
,FrameLabel->{}
,LabelStyle->{14, FontFamily-> "Arial"}]

Bessel J is regular almost everywhere, so it a physically valid solution.

BesselY[2,0]

Because Bessel Y is singular at the origin, it cannot be used in the core region. We can similarly consider Bessel I and K:

Plot[ Release[Table[BesselI[n, x] ,{n,0,5}]],{x,0.,10}
, PlotStyle->{{Hue[0.2]},{Hue[0.3]},{Hue[0.4]},{Hue[0.5]},{Hue[0.6]},{Hue[0.7]},{Hue[0.8]},{Hue[0.9]}} ]

Plot[ Release[Table[BesselK[n, x] ,{n,0,5}]],{x,0.1,2}
, PlotStyle->{{Hue[0.2]},{Hue[0.3]},{Hue[0.4]},{Hue[0.5]},{Hue[0.6]},{Hue[0.7]},{Hue[0.8]},{Hue[0.9]}} ]

Again, we see that Bessel I is infinite at infinity (like positive exponential), it cannot be used for the cladding region.

5.2.2 S matrix

Now we can solve for step index fiber:
In the core, ρ < a: Solution is (Subscript[J, l])[ξ], because (Subscript[Y, l])[ξ] goes to infinity at 0. In the cladding, solution is K because I goes to infinity at inifinity. Our interface S matrix is 
S= (((({{(Subscript[K, l]), (Subscript[I, l]), 0, 0}, {(Subscript[q, 2]) (Subscript[K, l])', (Subscript[q, 2]) (Subscript[I, l])', 0, 0}, {0, 0, (Subscript[K, l]), (Subscript[I, l])}, {0, 0, (Subscript[q, 2]) (Subscript[K, l])', (Subscript[q, 2]) (Subscript[I, l])'}})))^(-1)).(({{1, 0, 0, 0}, {0, (((Subscript[Q, 1]) (Subscript[ϵ, 2]))/((Subscript[Q, 2]) (Subscript[ϵ, 1]))), ((i l β ((Subscript[Q, 1])-(Subscript[Q, 2])))/(ρ (Subscript[k, 0]) (Subscript[Q, 2]) (Subscript[ϵ, 1]))), 0}, {0, 0, 1, 0}, {-((i l β ((Subscript[Q, 1])-(Subscript[Q, 2])))/(μ ρ (Subscript[k, 0]) (Subscript[Q, 2]))), 0, 0, ((Subscript[Q, 1])/(Subscript[Q, 2]))}})).(({{(Subscript[J, l]), (Subscript[Y, l]), 0, 0}, {(Subscript[q, 1]) (Subscript[J, l])', (Subscript[q, 1]) (Subscript[Y, l])', 0, 0}, {0, 0, (Subscript[J, l]), (Subscript[Y, l])}, {0, 0, (Subscript[q, 1]) (Subscript[J, l])', (Subscript[q, 1]) (Subscript[Y, l])'}}))         (5.2)
But here is a case of overwork: the solution has no Y and I because of other boundary conditions. We can still use the full S matrix. But why do the full algebra while simpler algebra exists? We'll show that S-matrix in this case can be reduced significantly.
For future use, remember that: the solution are  (Subscript[J, l])[(Sqrt[(Subscript[ϵ , 1])((Subscript[k, 0])^2) -(β^(2 ))]) ρ] in the core region. We further define 
(Subscript[q, 1])=(Sqrt[(Subscript[ϵ, 1])((Subscript[k, 0])^2) -(β^(2 ))]) and (Subscript[q, 2])=(Sqrt[(β^(2 ))-(Subscript[ϵ, 2])((Subscript[k, 0])^2)]), denoting (Subscript[ϵ, 1]) for the core and (Subscript[ϵ, 2]) for cladding. Their relation with (Subscript[Q, 1]) and (Subscript[Q, 2]) are:
(Subscript[Q, 1])=((Subscript[q, 1])^(2 )); (Subscript[Q, 2])=-((Subscript[q, 2])^2);
Note:           (Subscript[Q, 1])-(Subscript[Q, 2])=((Subscript[q, 1])^(2 ))+((Subscript[q, 2])^2)=((Subscript[k, 0])^2) ((Subscript[ϵ, 1])-(Subscript[ϵ, 2])) 
Here are the equations:

(({{1, 0, 0 , 0}, {0, ((I (Subscript[ϵ, 1])  (Subscript[k, 0]))/(Subscript[Q, 1])), -(( l)/((Subscript[Q, 1]) ρ)) β, 0 }, {0, 0, 1, 0 }, {(( β l)/((Subscript[Q, 1]) ρ)), 0, 0, ((I (Subscript[k, 0]) μ)/(Subscript[Q, 1]))}})).(({{(Subscript[J, l]), (Subscript[Y, l]), 0, 0}, {q1 (Subscript[J, l])', q1 (Subscript[Y, l])', 0, 0}, {0, 0, (Subscript[J, l]), (Subscript[Y, l])}, {0, 0, q1 (Subscript[J, l])', q1 (Subscript[Y, l])'}})).(({{c1}, {0}, {c2}, {0}}))==
	(({{1, 0, 0 , 0}, {0, ((I (Subscript[ϵ, 2])  (Subscript[k, 0]))/(Subscript[Q, 2])), -(( l)/((Subscript[Q, 2]) ρ)) β, 0 }, {0, 0, 1, 0 }, {(( β l)/((Subscript[Q, 2]) ρ)), 0, 0, ((I (Subscript[k, 0]) μ)/(Subscript[Q, 2]))}})).(({{(Subscript[K, l]), (Subscript[BI, l]), 0, 0}, {q2 (Subscript[K, l])', q2 (Subscript[BI, l])', 0, 0}, {0, 0, (Subscript[K, l]), (Subscript[BI, l])}, {0, 0, q2 (Subscript[K, l])', q2 (Subscript[BI, l])'}})).(({{c3}, {0}, {c4}, {0}}))

We can drop any Y and I column:

(({{1, 0, 0 , 0}, {0, ((I (Subscript[ϵ, 1])  (Subscript[k, 0]))/(Subscript[Q, 1])), -(( l)/((Subscript[Q, 1]) ρ)) β, 0 }, {0, 0, 1, 0 }, {(( β l)/((Subscript[Q, 1]) ρ)), 0, 0, ((I (Subscript[k, 0]) μ)/(Subscript[Q, 1]))}})).(({{(Subscript[J, l]), 0}, {q1 (Subscript[J, l])', 0}, {0, (Subscript[J, l])}, {0, q1 (Subscript[J, l])'}})).(({{c1}, {c2}}))==
	(({{1, 0, 0 , 0}, {0, ((I (Subscript[ϵ, 2])  (Subscript[k, 0]))/(Subscript[Q, 2])), -(( l)/((Subscript[Q, 2]) ρ)) β, 0 }, {0, 0, 1, 0 }, {(( β l)/((Subscript[Q, 2]) ρ)), 0, 0, ((I (Subscript[k, 0]) μ)/(Subscript[Q, 2]))}})).(({{(Subscript[K, l]), 0}, {q2 (Subscript[K, l])', 0}, {0, (Subscript[K, l])}, {0, q2 (Subscript[K, l])'}})).(({{c3}, {c4}}))

We can write:
                  (({{(Subscript[J, l]), 0, -(Subscript[K, l]), 0}, {((i  q1 (Subscript[k, 0]) (Subscript[ϵ, 1]) (Subscript[J, l]^′))/(Subscript[Q, 1])), -(( l β (Subscript[J, l]))/(ρ (Subscript[Q, 1]))), -((i  q2 (Subscript[k, 0]) (Subscript[ϵ, 2]) (Subscript[K, l]^′))/(Subscript[Q, 2])), (( l β (Subscript[K, l]))/(ρ (Subscript[Q, 2])))}, {0, (Subscript[J, l]), 0, -(Subscript[K, l])}, {(( l β (Subscript[J, l]))/(ρ (Subscript[Q, 1]))), ((i  q1 μ (Subscript[k, 0]) (Subscript[J, l]^′))/(Subscript[Q, 1])), -(( l β (Subscript[K, l]))/(ρ (Subscript[Q, 2]))), -((i  q2 μ (Subscript[k, 0]) (Subscript[K, l]^′))/(Subscript[Q, 2]))}})).(({{c1}, {c2}, {c3}, {c4}}))==(({{0}, {0}, {0}, {0}}))   (5.3)

The determinant has to vanish for nontrivial solution to exist.

This is the ultimate result from boundary condition that gives us the characteristic equation

5.3 Characteristic equation (algebra)

5.3.1 Equation set up

Simplify[Det[ (({{(Subscript[J, l]), 0, -(Subscript[K, l]), 0}, {((I  q1 (Subscript[k, 0]) (Subscript[ϵ, 1]) (Subscript[J, l]^′))/(Subscript[Q, 1])), -(( l β (Subscript[J, l]))/(ρ (Subscript[Q, 1]))), -((I  q2 (Subscript[k, 0]) (Subscript[ϵ, 2]) (Subscript[K, l]^′))/(Subscript[Q, 2])), (( l β (Subscript[K, l]))/(ρ (Subscript[Q, 2])))}, {0, (Subscript[J, l]), 0, -(Subscript[K, l])}, {(( l β (Subscript[J, l]))/(ρ (Subscript[Q, 1]))), ((I  q1 μ (Subscript[k, 0]) (Subscript[J, l]^′))/(Subscript[Q, 1])), -(( l β (Subscript[K, l]))/(ρ (Subscript[Q, 2]))), -((I  q2 μ (Subscript[k, 0]) (Subscript[K, l]^′))/(Subscript[Q, 2]))}}))  ]]

5.3.2 Algebraic reduction

((Subscript[J, l]^2) ((Subscript[k, 0]^2) (Subscript[Q, 1]^2) (q2^2) (Subscript[ϵ, 2]) μ (ρ^2) ((((((Subscript[K, l])))'))^2)-(l^2) ((((Subscript[Q, 1])-(Subscript[Q, 2])))^2) (β^2) (Subscript[K, l]^2))+(Subscript[k, 0]^2) (Subscript[Q, 2]^2) (q1^2) (Subscript[ϵ, 1]) μ (ρ^2) (Subscript[K, l]^2) ((((((Subscript[J, l])))'))^2)-(Subscript[k, 0]^2) (Subscript[Q, 1]) (Subscript[Q, 2]) q1 q2 ((Subscript[ϵ, 1])+(Subscript[ϵ, 2])) μ (ρ^2) (Subscript[J, l]) (Subscript[K, l]) (((Subscript[J, l])))' (((Subscript[K, l])))') /. μ-> 1 /. ((Subscript[Q, 1])-(Subscript[Q, 2]))-> ((Subscript[ϵ, 1])-(Subscript[ϵ, 2]))((Subscript[k, 0])^2)

(Subscript[J, l]^2) (Subscript[k, 0]^2) (Subscript[Q, 1]^2) (q2^2) (Subscript[ϵ, 2]) (ρ^2) ((((((Subscript[K, l])))'))^2)-(Subscript[k, 0]^2) (Subscript[Q, 1]) (Subscript[Q, 2]) q1 q2 ((Subscript[ϵ, 1])+(Subscript[ϵ, 2])) (ρ^2) (Subscript[J, l]) (Subscript[K, l]) (((Subscript[J, l])))' (((Subscript[K, l])))'+(Subscript[k, 0]^2) (Subscript[Q, 2]^2) (q1^2) (Subscript[ϵ, 1]) (ρ^2) (Subscript[K, l]^2) ((((((Subscript[J, l])))'))^2)==(l^2) (β^2) (Subscript[k, 0]^4)  ((((Subscript[ϵ, 1])-(Subscript[ϵ, 2])))^2) (Subscript[K, l]^2) (Subscript[J, l]^2)

FullSimplify[(Subscript[k, 0]^2) (Subscript[Q, 2]^2) (q1^2) (Subscript[ϵ, 1]) (ρ^2) (Subscript[K, l]^2) ((((((Subscript[J, l])))'))^2)-(Subscript[k, 0]^2) (Subscript[Q, 1]) (Subscript[Q, 2]) q1 q2 ((Subscript[ϵ, 1])+(Subscript[ϵ, 2])) (ρ^2) (Subscript[J, l]) (Subscript[K, l]) (((Subscript[J, l])))' (((Subscript[K, l])))'+(Subscript[k, 0]^2) (Subscript[Q, 1]^2) (q2^2) (Subscript[ϵ, 2]) (ρ^2) (Subscript[J, l]^2) ((((((Subscript[K, l])))'))^2)]==(Subscript[k, 0]^4) (l^2) (β^2) ((((Subscript[ϵ, 1])-(Subscript[ϵ, 2])))^2) (Subscript[J, l]^2) (Subscript[K, l]^2)

FullSimplify[(Subscript[k, 0]^2) (ρ^2) ((Subscript[Q, 2]) q1 (Subscript[K, l]) (((Subscript[J, l])))'-(Subscript[Q, 1]) q2 (Subscript[J, l]) (((Subscript[K, l])))') ((Subscript[Q, 2]) q1 (Subscript[ϵ, 1]) (Subscript[K, l]) (((Subscript[J, l])))'-(Subscript[Q, 1]) q2 (Subscript[ϵ, 2]) (Subscript[J, l]) (((Subscript[K, l])))')
	/. (Subscript[Q, 1])-> (q1^2) /. (Subscript[Q, 2])-> -(q2^2) ]

(q1 (Subscript[J, l]) (((Subscript[K, l])))'+q2 (Subscript[K, l]) (((Subscript[J, l])))') (q1 (Subscript[ϵ, 2]) (Subscript[J, l]) (((Subscript[K, l])))'+q2 (Subscript[ϵ, 1]) (Subscript[K, l]) (((Subscript[J, l])))') ==(1/((Subscript[k, 0]^2) (q1^2) (q2^2) (ρ^2) ))(Subscript[k, 0]^4) (l^2) (β^2) ((((Subscript[ϵ, 1])-(Subscript[ϵ, 2])))^2) (Subscript[J, l]^2) (Subscript[K, l]^2)

((q2/(Subscript[k, 0])) (Subscript[K, l]) (Subscript[J, l]^′)+(q1/(Subscript[k, 0])) (Subscript[J, l]) (Subscript[K, l]^′)) ((q2/(Subscript[k, 0])) (Subscript[K, l]) (Subscript[ϵ, 1]) (Subscript[J, l]^′)+(q1/(Subscript[k, 0])) (Subscript[J, l]) (Subscript[ϵ, 2]) (Subscript[K, l]^′))==((Subscript[k, 0]^4)/((q1^2) (q2^2)))((β^2)/(Subscript[k, 0]^2))((((l^2)((Subscript[ϵ, 1])-(Subscript[ϵ, 2])))^2)/((Subscript[k, 0]^2) (ρ^2)))(Subscript[J, l]^2) (Subscript[K, l]^2)

Define: (Subscript[δ, 1])= (q1/(Subscript[k, 0]))=(1/(Subscript[k, 0]))(Sqrt[(Subscript[ϵ, 1])((Subscript[k, 0])^2) -(β^(2 ))])=(Sqrt[(Subscript[ϵ, 1])-((β^(2 ))/((Subscript[k, 0])^2))])=(Sqrt[(Subscript[ϵ, 1])-((Subscript[n, eff])^2)]) 
            (Subscript[δ, 2])= (Sqrt[((Subscript[n, eff])^2)-(Subscript[ϵ, 2])])

((Subscript[δ, 2]) (Subscript[K, l]) (Subscript[J, l]^′)+(Subscript[δ, 1]) (Subscript[J, l]) (Subscript[K, l]^′)) ((Subscript[δ, 2])(Subscript[ϵ, 1]) (Subscript[K, l])  (Subscript[J, l]^′)+(Subscript[δ, 1]) (Subscript[ϵ, 2])(Subscript[J, l])  (Subscript[K, l]^′))==(((Subscript[n, eff])^2)/(((Subscript[δ, 1])^2) ((Subscript[δ, 2])^2)))((((l^2)((Subscript[ϵ, 1])-(Subscript[ϵ, 2])))^2)/((Subscript[k, 0]^2) (ρ^2)))(Subscript[J, l]^2) (Subscript[K, l]^2)

5.3.3 Characteristic equation

The characteristic equation here is thus:
                ((Subscript[δ, 2]) (Subscript[K, l]) (Subscript[J, l]^′)+(Subscript[δ, 1]) (Subscript[J, l]) (Subscript[K, l]^′)) ((Subscript[δ, 2])(Subscript[ϵ, 1]) (Subscript[K, l])  (Subscript[J, l]^′)+(Subscript[δ, 1]) (Subscript[ϵ, 2])(Subscript[J, l])  (Subscript[K, l]^′))==(((Subscript[n, eff])^2)/(((Subscript[δ, 1])^2) ((Subscript[δ, 2])^2)))((((l^2)((Subscript[ϵ, 1])-(Subscript[ϵ, 2])))^2)/((Subscript[k, 0]^2) (ρ^2)))(Subscript[J, l]^2) (Subscript[K, l]^2)

5.3.4 Numerical characteristic Eq. calculation

Source code: mode effective index demo

Manipulate[
If[INIT== 1,
CharG[nindex_List,λ_,ρ_,n_,l_]:=Module[{ϵ1,ϵ2,δ1, δ2,ξ1, ξ2,k0, J,K,Jp,Kp},
		{ϵ1,ϵ2}=nindex^2 ; k0=2*Pi/λ;
    	δ1=Sqrt[Abs[ϵ1 -n^2]]; δ2=Sqrt[Abs[ϵ2-n^2] ];
		ξ1=δ1*ρ *k0 ; ξ2=δ2*ρ* k0 ;
		J=BesselJ[l,ξ1] ; K=BesselK[l, ξ2];
		Jp=0.5*(BesselJ[l-1,ξ1] -BesselJ[l+1,ξ1]);
		Kp=-0.5*(BesselK[l-1,ξ2] + BesselK[l+1,ξ2]) ;
     (δ1 J Kp+ δ2 K Jp) (δ1  ϵ2 J Kp+δ2 ϵ1 K Jp)
	-(( (n^2))/( (k0^2) (ρ^2) ((δ1)^2) ((δ2)^2)))(l^2) (((ϵ1-ϵ2))^2) K^2 J^2   ] ; 
INIT=0 
,TextCell["Not ready"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]] 
;


Manipulate[
Column[{Button["Mode Search for Step Index Fiber"
,BaseStyle->{20, Bold,Black,FontFamily-> "Arial"}
,Background->Hue[0.5,0.7,1] ]

,If[RUN== 1,nclad=1.45;
Manipulate[
Indexprof={ncore,1.45};
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
Plot[Log[Abs[CharG[Indexprof,λ, rad,x,l]]] , {x,Indexprof[[2]]+dn,Indexprof[[1]]-dn}
		,PlotStyle->{{RGBColor[1,0.,0.]},{RGBColor[1,0,0]},{RGBColor[0,0,1]}} 
, PlotRange->{{Min[{nplow,nphi}],Max[{nplow,nphi}]+0.001}
, All}
,Frame-> True, GridLines->Automatic
,FrameLabel->{"Effective mode index","Log CharFunction"}
,LabelStyle->{14, FontFamily-> "Arial"} ]
,Grid[{{Control[{{l,0,Button["angular mode",BaseStyle->{16, Bold,Black,FontFamily-> "Calibri"},Background->Hue[0.15,1,1,1]]},0,5,1,ControlType-> SetterBar}]
,
Control[{{λ,1.55,Button["wavelength λ (μm)",BaseStyle->{16, Bold,Black,FontFamily-> "Calibri"},Background->Hue[0.036,0.5,1,1]]},0.5,3}]}
,{Control[{{ncore,1.48,Button["core index",BaseStyle->{16, Bold,Black,FontFamily-> "Calibri"},Background->Hue[0.3,1,1,1]]},1.46,1.52}]
,Control[{{rad,2.75,Button["core radius (μm)",BaseStyle->{16, Bold,Black,FontFamily-> "Calibri"},Background->Hue[0.45,1,1,1]]},2.5,4}]}}]
,Control[{{nplow,nclad+0.0005,"lower limit"},nclad+0.0005,ncore-0.001,ControlPlacement-> Bottom}]
,Control[{{nphi,ncore-0.0005,"upper limit"},nclad+0.001,ncore-0.0005,ControlPlacement-> Bottom}]
]
,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}}]
}]
]
,Row[{Button["Click once to initialize",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Red],Control[{{INIT,0,"->"},{0,1}}]
}]
]

Demo mode effective index

5.3.5  Concept of degeneracy

How many mode do we see in this curve below?

One. Obviously. Are we sure? We know that light has polarization. Supposed we have a rectangular waveguide like this (left), we know we have a TE and a TM mode. But if the rectangular WG becomes a perfect square, what happens to the modes? Can we distinguish TE and TM?

We can't. Because of 90-degree rotation symmetry, the two modes are indistinguishable by a rotation. They must therefore have the same propagation constant. This is the case of degeneracy. Similar to quantum mechanics, symmetry leads to degeneracy: many quantum states occupy the same energy, although each has a different quantum number. Now we have a perfect rotational symmetry in cylindrical waveguide. Would it be imaginable that each mode must have generacy also? Let's look at the characteristic equation:
((q2/(Subscript[k, 0])) (Subscript[K, l]) (Subscript[J, l]^′)+(q1/(Subscript[k, 0])) (Subscript[J, l]) (Subscript[K, l]^′)) ((q2/(Subscript[k, 0])) (Subscript[K, l]) (Subscript[ϵ, 1]) (Subscript[J, l]^′)+(q1/(Subscript[k, 0])) (Subscript[J, l]) (Subscript[ϵ, 2]) (Subscript[K, l]^′))==((Subscript[k, 0]^4)/((q1^2) (q2^2)))((β^2)/(Subscript[k, 0]^2))((((l^2)((Subscript[ϵ, 1])-(Subscript[ϵ, 2])))^2)/((Subscript[k, 0]^2) (ρ^2)))(Subscript[J, l]^2) (Subscript[K, l]^2)
Does it show any sign of degeneracy?
Recall the equation:
                     ρ F'[ρ] + (ρ^2) F''[ρ] +(((q )^2)(ρ^2)-(l^2)  )F[ρ ]== 0 
What if we change l into -l ? The equation remains the same. It's similar to changing β into -β, the equation is the same. The solution (e^(i β z)) describes a forward moving wave and (e^(-i β z)), a backward moving wave. They should have the same propagation constant. So (e^(i l ϕ)) describes one type of wave (we'll see later), then (e^(-i l ϕ)) must describes a symmetric wave with the same wave number. Indeed, we have the following mathematical relation:
 (Subscript[J, l])[x] = (((-1))^l)(Subscript[J, -l])[x]; (Subscript[J, l])'[x] = (((-1))^l)(Subscript[J, -l])'[x];(Subscript[K, l])[x]=(Subscript[K, -l])[x]; and (Subscript[K, l])'[x]=(Subscript[K, -l])'[x]
 Thus: with reverse sign: 
 ((q2/(Subscript[k, 0])) (Subscript[K, -l]) (Subscript[J, -l]^′)+(q1/(Subscript[k, 0])) (Subscript[J, -l]) (Subscript[K, -l]^′)) ((q2/(Subscript[k, 0])) (Subscript[K, -l]) (Subscript[ϵ, 1]) (Subscript[J, -l]^′)+(q1/(Subscript[k, 0])) (Subscript[J, -l]) (Subscript[ϵ, 2]) (Subscript[K, -l]^′))==((Subscript[k, 0]^4)/((q1^2) (q2^2)))((β^2)/(Subscript[k, 0]^2))((((l^2)((Subscript[ϵ, 1])-(Subscript[ϵ, 2])))^2)/((Subscript[k, 0]^2) (ρ^2)))(Subscript[J, -l]^2) (Subscript[K, -l]^2) 
  (((-1))^(2l))((q2/(Subscript[k, 0])) (Subscript[K, l]) (Subscript[J, l]^′)+(q1/(Subscript[k, 0])) (Subscript[J, l]) (Subscript[K, l]^′)) ((q2/(Subscript[k, 0])) (Subscript[K, l]) (Subscript[ϵ, 1]) (Subscript[J, l]^′)+(q1/(Subscript[k, 0])) (Subscript[J, l]) (Subscript[ϵ, 2]) (Subscript[K, -l]^′))==((Subscript[k, 0]^4)/((q1^2) (q2^2)))((β^2)/(Subscript[k, 0]^2))((((l^2)((Subscript[ϵ, 1])-(Subscript[ϵ, 2])))^2)/((Subscript[k, 0]^2) (ρ^2)))(Subscript[J, l]^2) (Subscript[K, l]^2) 
 we get back the same Eq. Thus, we have degeneracy: each propagation constant has at least 2 modes, corresponding to (e^(i l ϕ)) and (e^(-i l ϕ))

5.4 Determine the mode profile

Calculation of E and H field inside the fiber

5.4.1 Field expressions

Once we have the mode propagation constant, we can obtain the mode profile. Remember the eq.: 
(({{(Subscript[E, z]) }, {(Subscript[H, z])}, {ϵ (Subscript[E, ρ])}, {-(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))=(({{1, 0, 0, 0}, {0, 0, 1, 0}, {0, ((i β ϵ)/Q), ((-(Subscript[k, 0]) μ ϵ)/Q)(l/ρ) , 0}, {(( β l)/(Q ρ)), 0, 0, ((i (Subscript[k, 0]) μ )/Q)}, {ϵ  (Subscript[k, 0])(( l)/(Q ρ)), 0, 0 , ((i β)/Q)}, {0, ((i ϵ  (Subscript[k, 0]))/Q), -(( l)/(ρ Q)) β, 0 }})) (({{(Subscript[ψ, A]), (Subscript[ψ, B]), 0, 0}, {(Subscript[ψ, A])', (Subscript[ψ, B])', 0, 0}, {0, 0, (Subscript[ψ, A]), (Subscript[ψ, B])}, {0, 0, (Subscript[ψ, A])', (Subscript[ψ, B])'}})) (({{(Subscript[c, A])}, {(Subscript[c, B])}, {(Subscript[d, A])}, {(Subscript[d, B])}})) 
What this means is that if we determine the coefficients, we have all 6 field components. So we can write for inner core: 
(({{(Subscript[E, z]) }, {(Subscript[H, z])}, {(Subscript[E, ρ])}, {-(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))=(({{1, 0, 0, 0}, {0, 0, 1, 0}, {0, ((i β)/(Subscript[Q, 1])), ((-(Subscript[k, 0]) μ )/(Subscript[Q, 1]))(l/ρ) , 0}, {(( β l)/((Subscript[Q, 1]) ρ)), 0, 0, ((i (Subscript[k, 0]) μ )/(Subscript[Q, 1]))}, {(Subscript[ϵ, 1])  (Subscript[k, 0])(( l)/((Subscript[Q, 1]) ρ)), 0, 0 , ((i β)/(Subscript[Q, 1]))}, {0, ((i (Subscript[ϵ, 1])  (Subscript[k, 0]))/(Subscript[Q, 1])), -(( l)/(ρ (Subscript[Q, 1]))) β, 0 }}))(({{(Subscript[J, l]), (Subscript[Y, l]), 0, 0}, {q1 (Subscript[J, l])', q1 (Subscript[Y, l])', 0, 0}, {0, 0, (Subscript[J, l]), (Subscript[Y, l])}, {0, 0, q1 (Subscript[J, l])', q1 (Subscript[Y, l])'}})) (({{(Subscript[c, A])}, {(Subscript[c, B])}, {(Subscript[d, A])}, {(Subscript[d, B])}})) 
or: 
(({{(Subscript[E, z]) }, {(Subscript[H, z])}, {(Subscript[E, ρ])}, {-(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))=(({{(Subscript[J, l]), (Subscript[Y, l]), 0, 0}, {0, 0, (Subscript[J, l]), (Subscript[Y, l])}, {((i q1 β  (Subscript[J, l]^′))/(Subscript[Q, 1])), ((i q1 β  (Subscript[Y, l]^′))/(Subscript[Q, 1])), -((l (Subscript[J, l]) (Subscript[k, 0 ]))/(ρ (Subscript[Q, 1]))), -((l (Subscript[k, 0]) (Subscript[Y, l ]))/(ρ (Subscript[Q, 1])))}, {((l β (Subscript[J, l]))/(ρ (Subscript[Q, 1]))), ((l β (Subscript[Y, l]))/(ρ (Subscript[Q, 1]))), ((i q1 (Subscript[k, 0]) (Subscript[J, l]^′))/(Subscript[Q, 1])), ((i q1 (Subscript[k, 0]) (Subscript[Y, l]^′))/(Subscript[Q, 1]))}, {((l (Subscript[J, l]) (Subscript[k, 0]) (Subscript[ϵ, 1]))/(ρ (Subscript[Q, 1]))), ((l (Subscript[k, 0]) (Subscript[Y, l]) (Subscript[ϵ, 1]))/(ρ (Subscript[Q, 1]))), ((i q1 β (Subscript[J, l]^′))/(Subscript[Q, 1])), ((i q1 β (Subscript[Y, l]^′))/(Subscript[Q, 1]))}, {((i q1 (Subscript[k, 0]) (Subscript[ϵ, 1]) (Subscript[J, l]^′))/(Subscript[Q, 1])), ((i q1 (Subscript[k, 0]) (Subscript[ϵ, 1]) (Subscript[Y, l]^′))/(Subscript[Q, 1])), -((l β (Subscript[J, l]))/(ρ (Subscript[Q, 1]))), -((l β (Subscript[Y, l]))/(ρ (Subscript[Q, 1])))}})).(({{(Subscript[c, 1])}, {0}, {(Subscript[c, 2])}, {0}})) 
or: 
(({{(Subscript[E, z]) }, {(Subscript[H, z])}, {(Subscript[E, ρ])}, {-(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))=(({{(Subscript[J, l]), 0}, {0, (Subscript[J, l])}, {((i q1 β  (Subscript[J, l]^′))/(Subscript[Q, 1])), -((l (Subscript[J, l]) (Subscript[k, 0 ]))/(ρ (Subscript[Q, 1])))}, {((l β (Subscript[J, l]))/(ρ (Subscript[Q, 1]))), ((i q1 (Subscript[k, 0]) (Subscript[J, l]^′))/(Subscript[Q, 1]))}, {((l (Subscript[J, l]) (Subscript[k, 0]) (Subscript[ϵ, 1]))/(ρ (Subscript[Q, 1]))), ((i q1 β (Subscript[J, l]^′))/(Subscript[Q, 1]))}, {((i q1 (Subscript[k, 0]) (Subscript[ϵ, 1]) (Subscript[J, l]^′))/(Subscript[Q, 1])), -((l β (Subscript[J, l]))/(ρ (Subscript[Q, 1])))}})).(({{(Subscript[c, 1])}, {(Subscript[c, 2])}})) . 
For convenience and numerical calculation purpose, we will redefine the coefs:  
(({{(Subscript[D, 1])}, {(Subscript[D, 2])}}))=(1/(Subscript[δ, 1]))(({{(Subscript[c, 1])}, {(Subscript[c, 2])}})); Notation: (Subscript[δ, 1])=(Sqrt[(Subscript[ϵ, 1])-((Subscript[n, eff])^2)]); (Subscript[δ, 2])=(Sqrt[((Subscript[n, eff])^2)-(Subscript[ϵ, 2])]) ;
 (({{(Subscript[E, z]) }, {(Subscript[H, z])}, {(Subscript[E, ρ])}, {-(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))=(({{(Subscript[J, l]) (Subscript[δ, 1]), 0}, {0, (Subscript[J, l]) (Subscript[δ, 1])}, {i (Subscript[n, eff]) (Subscript[J, l]^′), -((l (Subscript[J, l]))/((Subscript[δ, 1])ρ (Subscript[k, 0])))}, {((l (Subscript[J, l]) (Subscript[n, eff]))/((Subscript[δ, 1])ρ (Subscript[k, 0]))), i  (Subscript[J, l]^′)}, {((l (Subscript[J, l]) (Subscript[ϵ, 1]))/((Subscript[δ, 1])ρ (Subscript[k, 0]))), i (Subscript[n, eff])  (Subscript[J, l]^′)}, {i  (Subscript[ϵ, 1]) (Subscript[J, l]^′), -((l (Subscript[J, l]) (Subscript[n, eff]))/((Subscript[δ, 1])ρ (Subscript[k, 0])))}})).(({{(Subscript[D, 1])}, {(Subscript[D, 2])}})) 
Similarly for the cladding: 
(({{(Subscript[E, z]) }, {(Subscript[H, z])}, {(Subscript[E, ρ])}, {-(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))=(({{(Subscript[K, l]), 0}, {0, (Subscript[K, l])}, {((i q2 β (Subscript[K, l]^′))/(Subscript[Q, 2])), -((l (Subscript[k, 0]) (Subscript[K, l ]))/(ρ (Subscript[Q, 2])))}, {((l β (Subscript[K, l]))/(ρ (Subscript[Q, 2]))), ((i q2 (Subscript[k, 0]) (Subscript[K, l]^′))/(Subscript[Q, 2]))}, {((l (Subscript[k, 0]) (Subscript[K, l]) (Subscript[ϵ, 2]))/(ρ (Subscript[Q, 2]))), ((i q2 β (Subscript[K, l]^′))/(Subscript[Q, 2]))}, {((i q2 (Subscript[k, 0]) (Subscript[ϵ, 2]) (Subscript[K, l]^′))/(Subscript[Q, 2])), -((l β (Subscript[K, l]))/(ρ (Subscript[Q, 2])))}})).(({{(Subscript[c, 3])}, {(Subscript[c, 4])}})) ; and we change the coef:
(({{(Subscript[E, z]) }, {(Subscript[H, z])}, {(Subscript[E, ρ])}, {-(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))=(({{-(Subscript[K, l]) (Subscript[δ, 2]), 0}, {0, -(Subscript[K, l]) (Subscript[δ, 2])}, {i (Subscript[n, eff]) (Subscript[K, l]^′), -((l (Subscript[K, l]))/((Subscript[δ, 2])ρ (Subscript[k, 0])))}, {((l (Subscript[K, l]) (Subscript[n, eff]))/((Subscript[δ, 2])ρ (Subscript[k, 0]))), i  (Subscript[K, l]^′)}, {((l (Subscript[K, l]) (Subscript[ϵ, 2]))/((Subscript[δ, 2])ρ (Subscript[k, 0]))), i (Subscript[n, eff]) (Subscript[K, l]^′)}, {i  (Subscript[ϵ, 2]) (Subscript[K, l]^′), -((l (Subscript[K, l]) (Subscript[n, eff]))/((Subscript[δ, 2])ρ (Subscript[k, 0])))}})).(({{(Subscript[D, 3])}, {(Subscript[D, 4])}})) ; where (({{(Subscript[D, 3])}, {(Subscript[D, 4])}}))= -(1/(Subscript[δ, 2]))(({{(Subscript[c, 3])}, {(Subscript[c, 4])}}))

The only requirement of the coefficients is that they satisfy the eq.
 (({{(Subscript[J, l]), 0, -(Subscript[K, l]), 0}, {((i  q1 (Subscript[k, 0]) (Subscript[ϵ, 1]) (Subscript[J, l]^′))/(Subscript[Q, 1])), -(( l β (Subscript[J, l]))/(ρ (Subscript[Q, 1]))), -((i  q2 (Subscript[k, 0]) (Subscript[ϵ, 2]) (Subscript[K, l]^′))/(Subscript[Q, 2])), (( l β (Subscript[K, l]))/(ρ (Subscript[Q, 2])))}, {0, (Subscript[J, l]), 0, -(Subscript[K, l])}, {(( l β (Subscript[J, l]))/(ρ (Subscript[Q, 1]))), ((i  q1 μ (Subscript[k, 0]) (Subscript[J, l]^′))/(Subscript[Q, 1])), -(( l β (Subscript[K, l]))/(ρ (Subscript[Q, 2]))), -((i  q2 μ (Subscript[k, 0]) (Subscript[K, l]^′))/(Subscript[Q, 2]))}})).(({{c1}, {c2}, {c3}, {c4}}))==(({{0}, {0}, {0}, {0}}))
 or:
(({{(Subscript[J, l]), 0, -(Subscript[K, l]), 0}, {((i  q1 (Subscript[k, 0]) (Subscript[ϵ, 1]) (Subscript[J, l]^′))/(Subscript[Q, 1])), -(( l β (Subscript[J, l]))/(ρ (Subscript[Q, 1]))), -((i  q2 (Subscript[k, 0]) (Subscript[ϵ, 2]) (Subscript[K, l]^′))/(Subscript[Q, 2])), (( l β (Subscript[K, l]))/(ρ (Subscript[Q, 2])))}, {0, (Subscript[J, l]), 0, -(Subscript[K, l])}, {(( l β (Subscript[J, l]))/(ρ (Subscript[Q, 1]))), ((i  q1 μ (Subscript[k, 0]) (Subscript[J, l]^′))/(Subscript[Q, 1])), -(( l β (Subscript[K, l]))/(ρ (Subscript[Q, 2]))), -((i  q2 μ (Subscript[k, 0]) (Subscript[K, l]^′))/(Subscript[Q, 2]))}})).(({{(Subscript[δ, 1]), 0, 0, 0}, {0, (Subscript[δ, 1]), 0, 0}, {0, 0, -(Subscript[δ, 2]), 0}, {0, 0, 0, -(Subscript[δ, 2])}}))(({{D1}, {D2}, {D3}, {D4}}))==(({{0}, {0}, {0}, {0}}))

(({{(Subscript[J, l]) (Subscript[δ, 1]), 0, (Subscript[K, l]) (Subscript[δ, 2]), 0}, {i (Subscript[ϵ, 1]) (Subscript[J, l]^′), -((l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1]))), -i (Subscript[ϵ, 2]) (Subscript[K, l]^′), ((l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2])))}, {0, (Subscript[J, l]) (Subscript[δ, 1]), 0, (Subscript[K, l]) (Subscript[δ, 2])}, {((l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1]))), i (Subscript[J, l]^′), -((l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2]))), -i (Subscript[K, l]^′)}})).(({{D1}, {D2}, {D3}, {D4}}))==(({{0}, {0}, {0}, {0}})) 
 which now is reduced to only three unknowns since its determinant =0.

5.4.2 Discussion: TE and TM mode

Let's consider some special case

(({{(Subscript[J, l]) (Subscript[δ, 1]), 0, (Subscript[K, l]) (Subscript[δ, 2]), 0}, {I (Subscript[ϵ, 1]) (Subscript[J, l]^′), -((l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1]))), -I (Subscript[ϵ, 2]) (Subscript[K, l]^′), ((l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2])))}, {0, (Subscript[J, l]) (Subscript[δ, 1]), 0, (Subscript[K, l]) (Subscript[δ, 2])}, {((l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1]))), I (Subscript[J, l]^′), -((l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2]))), -I (Subscript[K, l]^′)}})) /. l-> 0

Let's look at the characteristic Eq. for l=0

((q2/(Subscript[k, 0])) (Subscript[K, l]) (Subscript[J, l]^′)+(q1/(Subscript[k, 0])) (Subscript[J, l]) (Subscript[K, l]^′)) ((q2/(Subscript[k, 0])) (Subscript[K, l]) (Subscript[ϵ, 1]) (Subscript[J, l]^′)+(q1/(Subscript[k, 0])) (Subscript[J, l]) (Subscript[ϵ, 2]) (Subscript[K, l]^′))==((Subscript[k, 0]^4)/((q1^2) (q2^2)))((β^2)/(Subscript[k, 0]^2))((((l^2)((Subscript[ϵ, 1])-(Subscript[ϵ, 2])))^2)/((Subscript[k, 0]^2) (ρ^2)))(Subscript[J, l]^2) (Subscript[K, l]^2)  /. l->0

(((q2 (Subscript[K, 0]) (Subscript[J, 0]^′))/(Subscript[k, 0]))+((q1 (Subscript[J, 0]) (Subscript[K, 0]^′))/(Subscript[k, 0]))) (((q2 (Subscript[K, 0]) (Subscript[ϵ, 1]) (Subscript[J, 0]^′))/(Subscript[k, 0]))+((q1 (Subscript[J, 0]) (Subscript[ϵ, 2]) (Subscript[K, 0]^′))/(Subscript[k, 0])))==0

Whenever we see a product like above, the Eq. can be satisfied by either: 
(((q2 (Subscript[K, 0]) (Subscript[J, 0]^′))/(Subscript[k, 0]))+((q1 (Subscript[J, 0]) (Subscript[K, 0]^′))/(Subscript[k, 0]))) ==0 or (((q2 (Subscript[K, 0]) (Subscript[ϵ, 1]) (Subscript[J, 0]^′))/(Subscript[k, 0]))+((q1 (Subscript[J, 0]) (Subscript[ϵ, 2]) (Subscript[K, 0]^′))/(Subscript[k, 0])))==0. We recall that this means each Eq. corresponds to some symmetry.  
(({{(Subscript[J, 0]) (Subscript[δ, 1]), 0, (Subscript[K, 0]) (Subscript[δ, 2]), 0}, {i (Subscript[ϵ, 1]) (Subscript[J, 0]^′), 0, -i (Subscript[ϵ, 2]) (Subscript[K, 0]^′), 0}, {0, (Subscript[J, 0]) (Subscript[δ, 1]), 0, (Subscript[K, 0]) (Subscript[δ, 2])}, {0, i (Subscript[J, 0]^′), 0, -i (Subscript[K, 0]^′)}})).(({{D1}, {D2}, {D3}, {D4}}))
			             =(({{D1 (Subscript[J, 0]) (Subscript[δ, 1])+D3 (Subscript[K, 0]) (Subscript[δ, 2])}, {i D1 (Subscript[ϵ, 1]) (Subscript[J, 0]^′)-i D3 (Subscript[ϵ, 2]) (Subscript[K, 0]^′)}, {D2 (Subscript[J, 0]) (Subscript[δ, 1])+D4 (Subscript[K, 0]) (Subscript[δ, 2])}, {i D2 (Subscript[J, 0]^′)-i D4 (Subscript[K, 0]^′)}}))==(({{0}, {0}, {0}, {0}}))
We can just reduce to 2 Eqs.:
(({{D1 (Subscript[J, 0]) (Subscript[δ, 1])+D3 (Subscript[K, 0]) (Subscript[δ, 2])}, {D2 (Subscript[J, 0]) (Subscript[δ, 1])+D4 (Subscript[K, 0]) (Subscript[δ, 2])}}))==(({{0}, {0}}))
Notice that D1 and D3 are related and D2 and D3 are related, but D1 and D2 are independent of each other. For example, we can choose D1=1; D3=-(Subscript[J, 0]) (Subscript[δ, 1])/(Subscript[K, 0]) (Subscript[δ, 2]) and let D2=0; or D2=1; D4=-(Subscript[J, 0]) (Subscript[δ, 1])/(Subscript[K, 0]) (Subscript[δ, 2]) and D1=D3=0. What does this mean?

(({{(Subscript[E, z]) }, {(Subscript[H, z])}, {(Subscript[E, ρ])}, {-(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))=(({{(Subscript[J, l]) (Subscript[δ, 1]), 0}, {0, (Subscript[J, l]) (Subscript[δ, 1])}, {i (Subscript[n, eff]) (Subscript[J, l]^′), -((l (Subscript[J, l ]))/((Subscript[δ, 1])ρ (Subscript[k, 0])))}, {((l (Subscript[J, l]) (Subscript[n, eff]))/((Subscript[δ, 1])ρ (Subscript[k, 0]))), i  (Subscript[J, l]^′)}, {((l (Subscript[J, l]) (Subscript[ϵ, 1]))/((Subscript[δ, 1])ρ (Subscript[k, 0]))), i (Subscript[n, eff])  (Subscript[J, l]^′)}, {i  (Subscript[ϵ, 1]) (Subscript[J, l]^′), -((l (Subscript[J, l]) (Subscript[n, eff]))/((Subscript[δ, 1])ρ (Subscript[k, 0])))}})).(({{1}, {0}}))=  (({{(Subscript[J, l]) (Subscript[δ, 1])}, {0}, {i (Subscript[n, eff]) (Subscript[J, l]^′)}, {((l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1])))}, {((l (Subscript[J, l]) (Subscript[ϵ, 1]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1])))}, {i (Subscript[ϵ, 1]) (Subscript[J, l]^′)}}))
Hz=0: this is TM mode: the magnetic field is purely transverse. Similarly, the other case give purely TE mode.

But if l is not equal 0, could we have TE or TM mode?  Suppose we choose D1=0

(({{(Subscript[J, l]) (Subscript[δ, 1]), 0, (Subscript[K, l]) (Subscript[δ, 2]), 0}, {I (Subscript[ϵ, 1]) (Subscript[J, l]^′), -((l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1]))), -I (Subscript[ϵ, 2]) (Subscript[K, l]^′), ((l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2])))}, {0, (Subscript[J, l]) (Subscript[δ, 1]), 0, (Subscript[K, l]) (Subscript[δ, 2])}, {((l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1]))), I (Subscript[J, l]^′), -((l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2]))), -I (Subscript[K, l]^′)}})) .(({{0}, {D2}, {D3}, {D4}}))==(({{0}, {0}, {0}, {0}}))

(({{D3 (Subscript[K, l]) (Subscript[δ, 2])}, {-((D2 l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1])))+((D4 l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2])))-I D3 (Subscript[ϵ, 2]) (Subscript[K, l]^′)}, {D2 (Subscript[J, l]) (Subscript[δ, 1])+D4 (Subscript[K, l]) (Subscript[δ, 2])}, {-((D3 l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2])))+I D2 (Subscript[J, l]^′)-I D4 (Subscript[K, l]^′)}}))==(({{0}, {0}, {0}, {0}})) /. D3->0

Solve[ ((D4 l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2])))-((D2 l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1])))==0 ,{D4}]
Solve[ D2 (Subscript[J, l]) (Subscript[δ, 1])+D4 (Subscript[K, l]) (Subscript[δ, 2])==0 ,{D4}]
Solve[ I I D2 (Subscript[J, l]^′)-I D4 (Subscript[K, l]^′)==0 ,{D4}]

Simplify[Det[(({{-(( l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1]))), (( l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2])))}, {(Subscript[J, l]) (Subscript[δ, 1]), (Subscript[K, l]) (Subscript[δ, 2])}}))]]==0

We see that three simultaneous equation must be satisfied, and they can't be all the time! Thus, unlike slab waveguide, fiber mode is much more complicated and cannot be classified in terms of some simple geometrical symmetry.
Thus, it is crucial to choose the coefficients that correspond to physically insightful solution. It is not useful to blindly follow numerical method and get numerical results without getting the "physics" of the problem.

5.4.2 HE and EH mode

Let's look at the field mode expression:
 (({{(Subscript[E, z]) }, {(Subscript[H, z])}, {(Subscript[E, ρ])}, {-(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))=(({{(Subscript[J, l]) (Subscript[δ, 1]), 0}, {0, (Subscript[J, l]) (Subscript[δ, 1])}, {i (Subscript[n, eff]) (Subscript[J, l]^′), -((l (Subscript[J, l ]))/((Subscript[δ, 1])ρ (Subscript[k, 0])))}, {((l (Subscript[J, l]) (Subscript[n, eff]))/((Subscript[δ, 1])ρ (Subscript[k, 0]))), i  (Subscript[J, l]^′)}, {((l (Subscript[J, l]) (Subscript[ϵ, 1]))/((Subscript[δ, 1])ρ (Subscript[k, 0]))), i (Subscript[n, eff])  (Subscript[J, l]^′)}, {i  (Subscript[ϵ, 1]) (Subscript[J, l]^′), -((l (Subscript[J, l]) (Subscript[n, eff]))/((Subscript[δ, 1])ρ (Subscript[k, 0])))}})).(({{(Subscript[D, 1])}, {(Subscript[D, 2])}})) 
 We recall ealier that we have a degeneracy. We can make any linear combination of (e^(i l ϕ)) and (e^(-i l ϕ)) to make any other solution. So it is natural to ask if we should choose coefficient in certain way that reflects some sort of "natural" solution.
The answer is yes. We can make any transformation: (({{(Subscript[D, 1])}, {(Subscript[D, 2])}}))=(({{(Subscript[s, 11]), (Subscript[s, 12])}, {(Subscript[s, 21]), (Subscript[s, 22])}}))(({{(Subscript[W, 1])}, {(Subscript[W, 2])}})) we like.
Let's consider the follow: (({{(Subscript[D, 1])}, {(Subscript[D, 2])}}))=(({{(1/(i (Subscript[n, eff]))), (1/(Subscript[ϵ, 1]))}, {-1, (1/(i (Subscript[n, eff])))}}))(({{(Subscript[W, 1])}, {(Subscript[W, 2])}})), which gives:
(({{(Subscript[E, z]) }, {(Subscript[H, z])}, {(Subscript[E, ρ])}, {-(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))=(({{-((i (Subscript[J, l]) (Subscript[δ, 1]))/(Subscript[n, eff])), (((Subscript[J, l]) (Subscript[δ, 1]))/(Subscript[ϵ, 1]))}, {-(Subscript[J, l]) (Subscript[δ, 1]), -((i (Subscript[J, l]) (Subscript[δ, 1]))/(Subscript[n, eff]))}, {((l (Subscript[J, l]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1])))+(Subscript[J, l]^′), ((i l (Subscript[J, l]))/(ρ (Subscript[k, 0]) (Subscript[n, eff]) (Subscript[δ, 1])))+((i (Subscript[n, eff]) (Subscript[J, l]^′))/(Subscript[ϵ, 1]))}, {-((i l (Subscript[J, l]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1])))-i (Subscript[J, l]^′), ((l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1]) (Subscript[ϵ, 1])))+((Subscript[J, l]^′)/(Subscript[n, eff]))}, {-((i l (Subscript[J, l]) (Subscript[ϵ, 1]))/(ρ (Subscript[k, 0]) (Subscript[n, eff]) (Subscript[δ, 1])))-i (Subscript[n, eff]) (Subscript[J, l]^′), ((l (Subscript[J, l]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1])))+(Subscript[J, l]^′)}, {((l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1])))+(((Subscript[ϵ, 1]) (Subscript[J, l]^′))/(Subscript[n, eff])), ((i l (Subscript[J, l]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1])))+i (Subscript[J, l]^′)}}))(({{(Subscript[W, 1])}, {(Subscript[W, 2])}})) . We recall further that: 
 (Subscript[J, l]^′)+ ((l (Subscript[J, l]))/x)=(Subscript[J, l-1])  ;  (Subscript[J, l]^′)- ((l (Subscript[J, l]))/x)=-(Subscript[J, l+1])  ;  (Subscript[K, l]^′)+ ((l (Subscript[K, l]))/x)=-(Subscript[K, l-1])  ;  (Subscript[K, l]^′)- ((l (Subscript[K, l]))/x)=-(Subscript[K, l+1])   
It becomes: 
(({{(Subscript[E, z]) }, {(Subscript[H, z])}, {(Subscript[E, ρ])}, {-(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))=(({{-((i (Subscript[J, l]) (Subscript[δ, 1]))/(Subscript[n, eff])), (((Subscript[J, l]) (Subscript[δ, 1]))/(Subscript[ϵ, 1]))}, {-(Subscript[J, l]) (Subscript[δ, 1]), -((i (Subscript[J, l]) (Subscript[δ, 1]))/(Subscript[n, eff]))}, {(Subscript[J, l-1]), ((i l (Subscript[J, l]))/(ρ (Subscript[k, 0]) (Subscript[n, eff]) (Subscript[δ, 1])))+((i (Subscript[n, eff]) (Subscript[J, l]^′))/(Subscript[ϵ, 1]))}, {-i (Subscript[J, l-1]), ((l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1]) (Subscript[ϵ, 1])))+((Subscript[J, l]^′)/(Subscript[n, eff]))}, {-((i l (Subscript[J, l]) (Subscript[ϵ, 1]))/(ρ (Subscript[k, 0]) (Subscript[n, eff]) (Subscript[δ, 1])))-i (Subscript[n, eff]) (Subscript[J, l]^′), (Subscript[J, l-1])}, {((l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1])))+(((Subscript[ϵ, 1]) (Subscript[J, l]^′))/(Subscript[n, eff])), i (Subscript[J, l-1])}}))(({{(Subscript[W, 1])}, {(Subscript[W, 2])}}))

Similarly for the cladding
(({{(Subscript[E, z]) }, {(Subscript[H, z])}, {(Subscript[E, ρ])}, {-(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))=(({{-(Subscript[K, l]) (Subscript[δ, 2]), 0}, {0, -(Subscript[K, l]) (Subscript[δ, 2])}, {i (Subscript[n, eff]) (Subscript[K, l]^′), -((l (Subscript[K, l]))/((Subscript[δ, 2])ρ (Subscript[k, 0])))}, {((l (Subscript[K, l]) (Subscript[n, eff]))/((Subscript[δ, 2])ρ (Subscript[k, 0]))), i  (Subscript[K, l]^′)}, {((l (Subscript[K, l]) (Subscript[ϵ, 2]))/((Subscript[δ, 2])ρ (Subscript[k, 0]))), i (Subscript[n, eff]) (Subscript[K, l]^′)}, {i  (Subscript[ϵ, 2]) (Subscript[K, l]^′), -((l (Subscript[K, l]) (Subscript[n, eff]))/((Subscript[δ, 2])ρ (Subscript[k, 0])))}})).(({{(1/(i (Subscript[n, eff]))), (1/(Subscript[ϵ, 2]))}, {-1, (1/(i (Subscript[n, eff])))}}))(({{(Subscript[W, 3])}, {(Subscript[W, 4])}}))
(({{(Subscript[E, z]) }, {(Subscript[H, z])}, {(Subscript[E, ρ])}, {-(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))=(({{((i (Subscript[K, l]) (Subscript[δ, 2]))/(Subscript[n, eff])), -(((Subscript[K, l]) (Subscript[δ, 2]))/(Subscript[ϵ, 2]))}, {(Subscript[K, l]) (Subscript[δ, 2]), ((i (Subscript[K, l]) (Subscript[δ, 2]))/(Subscript[n, eff]))}, {-(Subscript[K, l-1]), ((i l (Subscript[K, l]))/(ρ (Subscript[k, 0]) (Subscript[n, eff]) (Subscript[δ, 2])))+((i (Subscript[n, eff]) (Subscript[K, l]^′))/(Subscript[ϵ, 2]))}, {i (Subscript[K, l-1]), ((l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2]) (Subscript[ϵ, 2])))+((Subscript[K, l]^′)/(Subscript[n, eff]))}, {-((i l (Subscript[K, l]) (Subscript[ϵ, 2]))/(ρ (Subscript[k, 0]) (Subscript[n, eff]) (Subscript[δ, 2])))-i (Subscript[n, eff]) (Subscript[K, l]^′), -(Subscript[K, l-1])}, {((l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2])))+(((Subscript[ϵ, 2]) (Subscript[K, l]^′))/(Subscript[n, eff])), -i (Subscript[K, l-1])}})).(({{(Subscript[W, 3])}, {(Subscript[W, 4])}}))

(({{(Subscript[J, l]) (Subscript[δ, 1]), 0, (Subscript[K, l]) (Subscript[δ, 2]), 0}, {I (Subscript[ϵ, 1]) (Subscript[J, l]^′), -((l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1]))), -I (Subscript[ϵ, 2]) (Subscript[K, l]^′), ((l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2])))}, {0, (Subscript[J, l]) (Subscript[δ, 1]), 0, (Subscript[K, l]) (Subscript[δ, 2])}, {((l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1]))), I (Subscript[J, l]^′), -((l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2]))), -I (Subscript[K, l]^′)}})) .(({{(1/(I (Subscript[n, eff]))), (1/(Subscript[ϵ, 1])), 0, 0}, {-1, (1/(I (Subscript[n, eff]))), 0, 0}, {0, 0, (1/(I (Subscript[n, eff]))), (1/(Subscript[ϵ, 2]))}, {0, 0, -1, (1/(I (Subscript[n, eff])))}})).(({{W1}, {W2}, {W3}, {W4}}))==(({{0}, {0}, {0}, {0}}))

(({{-((I (Subscript[J, l]) (Subscript[δ, 1]))/(Subscript[n, eff])), (((Subscript[J, l]) (Subscript[δ, 1]))/(Subscript[ϵ, 1])), -((I (Subscript[K, l]) (Subscript[δ, 2]))/(Subscript[n, eff])), (((Subscript[K, l]) (Subscript[δ, 2]))/(Subscript[ϵ, 2]))}, {((l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1])))+(((Subscript[ϵ, 1]) (Subscript[J, l]^′))/(Subscript[n, eff])), I (Subscript[J, l-1]), -((l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2])))-(((Subscript[ϵ, 2]) (Subscript[K, l]^′))/(Subscript[n, eff])), I (Subscript[K, l-1])}, {-(Subscript[J, l]) (Subscript[δ, 1]), -((I (Subscript[J, l]) (Subscript[δ, 1]))/(Subscript[n, eff])), -(Subscript[K, l]) (Subscript[δ, 2]), -((I (Subscript[K, l]) (Subscript[δ, 2]))/(Subscript[n, eff]))}, {-I (Subscript[J, l-1]), ((l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1]) (Subscript[ϵ, 1])))+((Subscript[J, l]^′)/(Subscript[n, eff])), -I (Subscript[K, l-1]), -((l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2]) (Subscript[ϵ, 2])))-((Subscript[K, l]^′)/(Subscript[n, eff]))}})).(({{W1}, {W2}, {W3}, {W4}}))==(({{0}, {0}, {0}, {0}}))

We can have a Ez dominant mode, which we can call EH or the reverse and call HE mode. This is the equivalent to TE or TM mode in planar waveguide.

5.4.4 Summary

We see that by choosing appropriate combination, we have "simple" description of the field profile. The field above belongs to the HE or EH class: they are called HE or EH depending on whether the Hz or Ez component is dominant. Earlier on, we have seen TE or TM.
The modes are naturally grouped together by their propagation constants, which cluster together (sometime very closely, almost degenerate). Inspite of that fact that their polarization are different, their ultimate field profile are similar. When they can be treated as "nearly" degenerate, they can actually be transformed into modes that have linear polarization just like a paraxial beam. This approximation gives rise to the classification (Subscript[LP, l m]) (linear polarization). This approximation simplifies the problem greatly and is used very often. But remember that it is still only an approx.
More about them later.

5.5 LP approximation

Link to paraxial approx. Skip.

6 Illustration for step-index fiber

6.1 Review: The number of modes and mode search

Demo mode effective index calc (with λ)

6.2 Key concepts: single-mode vs. multimode, cut-off wavelength

6.2.1 Single mode vs. multi-mode - change of diameter discussion

6.2.2 Single mode vs. multi-mode - change of index discussion

6.2.3 Cut-off wavelength  discussion

Does a fiber or waveguide support all wavelengths? Not necessarily!
When we say large or small diameter, it means relatively to wavelength. Intuitively, given the same size fiber, it appears large to very small wavelength, and small to large wavelength. It means that a fiber may be multimode to small wavelength and may not even have anymode for large wavelength!

What do we see? ONLY one value l=1 has a solution of the characteristic function, which means only one mode. This is a single-mode fiber. 
What if the wavelength is NOT 1.55? longer or shorter?

Notice that there is a cut-off wavelength: beyond which the wavelength is just "too big" to go through a waveguide.

6.3 Calculation package

CharG[nindex_List,λ_,ρ_,n_,l_]:=Module[{ϵ1,ϵ2,δ1, δ2,ξ1, ξ2,k0, J,K,Jp,Kp},
		{ϵ1,ϵ2}=nindex^2 ; k0=2*Pi/λ;
    	δ1=Sqrt[Abs[ϵ1 -n^2]]; δ2=Sqrt[Abs[ϵ2-n^2] ];
		ξ1=δ1*ρ *k0 ; ξ2=δ2*ρ* k0 ;
		J=BesselJ[l,ξ1] ; K=BesselK[l, ξ2];
		Jp=0.5*(BesselJ[l-1,ξ1] -BesselJ[l+1,ξ1]);
		Kp=-0.5*(BesselK[l-1,ξ2] + BesselK[l+1,ξ2]) ;
     (δ1 J Kp+ δ2 K Jp) (δ1  ϵ2 J Kp+δ2 ϵ1 K Jp)
	-(( (n^2))/( (k0^2) (ρ^2) ((δ1)^2) ((δ2)^2)))(l^2) (((ϵ1-ϵ2))^2) K^2 J^2   ]    ;

Coeff[nindex_List,λ_,ρ_,n_,l_]:=Module[{ϵ1,ϵ2,q1, q2,ξ1, ξ2,k0, J,K,Jp,Kp},
		{ϵ1,ϵ2}=nindex^2 ; k0=2*Pi/λ;
    	q1=Sqrt[Abs[ϵ1 -n^2]]; q2=Sqrt[Abs[ϵ2-n^2] ];
		ξ1=q1*ρ *k0 ; ξ2=q2*ρ* k0 ;
		J=BesselJ[l,ξ1] ; K=BesselK[l, ξ2];
		Jp=0.5*(BesselJ[l-1,ξ1] -BesselJ[l+1,ξ1]);
		Kp=-0.5*(BesselK[l-1,ξ2] + BesselK[l+1,ξ2]) ;
		{ 1, ((I l J K n (ϵ1-ϵ2))/(ρ k0 q1 q2 (q2 K   Jp+q1 J Kp))),
			-((J q1)/(K q2)) , 
                             -((I l (J^2)  n (ϵ1-ϵ2))/(ρ k0 (q2^2) (q2 K Jp+q1 J Kp)))} 
 ]

FieldFiber[coor_List, nindex_List,λ_,rad_,n_,l_,coef_List]:=Module[{ρ, ϕ, ϵ1,ϵ2,q1, q2,ξ1, ξ2,k0, J,K,Jp,Kp, field},
		{ρ, ϕ}=coor ;
			{ϵ1,ϵ2}=nindex^2 ; k0=2*Pi/λ;
    	q1=Sqrt[Abs[ϵ1 -n^2]]; q2=Sqrt[Abs[ϵ2-n^2] ];
		ξ1=q1*ρ *k0 ; ξ2=q2*ρ* k0 ;
		J=BesselJ[l,ξ1] ; K=BesselK[l, ξ2];
		Jp=0.5*(BesselJ[l-1,ξ1] -BesselJ[l+1,ξ1]);
		Kp=-0.5*(BesselK[l-1,ξ2] + BesselK[l+1,ξ2]) ;
		field=If[ Abs[ρ]<=rad,
			 { J  q1
				, coef[[2]]*  J  q1
				, I n Jp -coef[[2]]*l *J/ξ1
				, -l *n J/ξ1 -coef[[2]] * I *Jp
				, -l *ϵ1  J/ξ1+coef[[2]]*I n Jp
				, I ϵ1 Jp-coef[[2]]*l*J*n/ξ1}
			,{-coef[[3]]*K q2
				,-coef[[4]]*K q2
				,coef[[3]]*I n Kp-coef[[4]]*l K/ξ2
				,-coef[[3]]*l n K/ξ2-coef[[4]]*I *Kp
				,coef[[3]]*l K ϵ2/ξ2+coef[[4]] I n Kp
				,coef[[3]] I ϵ2 Kp-coef[[4]]l K n/ξ2}] ;
		{I field * (E^(I ϕ l)),I field *(E^(-I ϕ l))}
		          ]

5.6.1 Set up notes

Let's choose D1=1
(({{(Subscript[J, l]) (Subscript[δ, 1]), 0, (Subscript[K, l]) (Subscript[δ, 2]), 0}, {i (Subscript[ϵ, 1]) (Subscript[J, l]^′), -((l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1]))), -i (Subscript[ϵ, 2]) (Subscript[K, l]^′), ((l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2])))}, {0, (Subscript[J, l]) (Subscript[δ, 1]), 0, (Subscript[K, l]) (Subscript[δ, 2])}, {((l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1]))), i (Subscript[J, l]^′), -((l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2]))), -i (Subscript[K, l]^′)}})).(({{1}, {D2}, {D3}, {D4}}))

(({{(Subscript[J, l]) (Subscript[δ, 1]), 0, (Subscript[K, l]) (Subscript[δ, 2]), 0}, {I (Subscript[ϵ, 1]) (Subscript[J, l]^′), -((l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1]))), -I (Subscript[ϵ, 2]) (Subscript[K, l]^′), ((l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2])))}, {0, (Subscript[J, l]) (Subscript[δ, 1]), 0, (Subscript[K, l]) (Subscript[δ, 2])}, {((l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1]))), I (Subscript[J, l]^′), -((l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2]))), -I (Subscript[K, l]^′)}})).(({{1}, {D2}, {D3}, {D4}}))

Solve[ (Subscript[J, l]) (Subscript[δ, 1])+D3 (Subscript[K, l]) (Subscript[δ, 2])==0 ,{D3}]

Simplify[((l (Subscript[J, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 1])))-((D3 l (Subscript[K, l]) (Subscript[n, eff]))/(ρ (Subscript[k, 0]) (Subscript[δ, 2])))+I D2 (Subscript[J, l]^′)-I D4 (Subscript[K, l]^′) /. D3->-(((Subscript[J, l]) (Subscript[δ, 1]))/((Subscript[K, l]) (Subscript[δ, 2])))]

((l (Subscript[J, l]) (Subscript[n, eff]) ((Subscript[δ, 1]^2)+(Subscript[δ, 2]^2)))/(ρ (Subscript[k, 0]) (Subscript[δ, 1]) (Subscript[δ, 2]^2)))+I (D2 (Subscript[J, l]^′)-D4 (Subscript[K, l]^′))/. (Subscript[δ, 1]^2)+(Subscript[δ, 2]^2)->(Subscript[ϵ, 1])-(Subscript[ϵ, 2])

Solve[{ D2 (Subscript[J, l]) (Subscript[δ, 1])+D4 (Subscript[K, l]) (Subscript[δ, 2])==0 ,
		((l (Subscript[J, l]) (Subscript[n, eff]) ((Subscript[ϵ, 1])-(Subscript[ϵ, 2])))/(ρ (Subscript[k, 0]) (Subscript[δ, 1]) (Subscript[δ, 2]^2)))+I (D2 (Subscript[J, l]^′)-D4 (Subscript[K, l]^′))==0},
		           {D2,D4}]

Core:
 (({{(Subscript[E, z]) }, {(Subscript[H, z])}, {(Subscript[E, ρ])}, {-(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))=(({{(Subscript[J, l]) (Subscript[δ, 1]), 0}, {0, (Subscript[J, l]) (Subscript[δ, 1])}, {i (Subscript[n, eff]) (Subscript[J, l]^′), -((l (Subscript[J, l]))/((Subscript[δ, 1])ρ (Subscript[k, 0])))}, {((l (Subscript[J, l]) (Subscript[n, eff]))/((Subscript[δ, 1])ρ (Subscript[k, 0]))), i  (Subscript[J, l]^′)}, {((l (Subscript[J, l]) (Subscript[ϵ, 1]))/((Subscript[δ, 1])ρ (Subscript[k, 0]))), i (Subscript[n, eff])  (Subscript[J, l]^′)}, {i  (Subscript[ϵ, 1]) (Subscript[J, l]^′), -((l (Subscript[J, l]) (Subscript[n, eff]))/((Subscript[δ, 1])ρ (Subscript[k, 0])))}})).(({{1}, {((i l (Subscript[(OverHat[J]), l]) (Subscript[(OverHat[K]), l]) (Subscript[n, eff]) ((Subscript[ϵ, 1])-(Subscript[ϵ, 2])))/(a (Subscript[k, 0]) (Subscript[δ, 1]) (Subscript[δ, 2]) ((Subscript[(OverHat[K]), l]) (Subscript[δ, 2]) (Subscript[(OverHat[J]), l]^′)+(Subscript[(OverHat[J]), l]) (Subscript[δ, 1]) (Subscript[(OverHat[K]), l]^′))))}})) 
 Cladding:
 (({{(Subscript[E, z]) }, {(Subscript[H, z])}, {(Subscript[E, ρ])}, {-(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))=(({{-(Subscript[K, l]) (Subscript[δ, 2]), 0}, {0, -(Subscript[K, l]) (Subscript[δ, 2])}, {i (Subscript[n, eff]) (Subscript[K, l]^′), -((l (Subscript[K, l]))/((Subscript[δ, 2])ρ (Subscript[k, 0])))}, {((l (Subscript[K, l]) (Subscript[n, eff]))/((Subscript[δ, 2])ρ (Subscript[k, 0]))), i  (Subscript[K, l]^′)}, {((l (Subscript[K, l]) (Subscript[ϵ, 2]))/((Subscript[δ, 2])ρ (Subscript[k, 0]))), i (Subscript[n, eff]) (Subscript[K, l]^′)}, {i  (Subscript[ϵ, 2]) (Subscript[K, l]^′), -((l (Subscript[K, l]) (Subscript[n, eff]))/((Subscript[δ, 2])ρ (Subscript[k, 0])))}})).(({{-(((Subscript[J, l]) (Subscript[δ, 1]))/((Subscript[K, l]) (Subscript[δ, 2])))}, {-((i l (Subscript[J, l]^2) (Subscript[n, eff]) ((Subscript[ϵ, 1])-(Subscript[ϵ, 2])))/(a (Subscript[k, 0]) (Subscript[δ, 2]^2) ((Subscript[K, l]) (Subscript[δ, 2]) (Subscript[J, l]^′)+(Subscript[J, l]) (Subscript[δ, 1]) (Subscript[K, l]^′))))}}))

5.6.2 Load routines

CharG[nindex_List,λ_,ρ_,n_,l_]:=Module[{ϵ1,ϵ2,δ1, δ2,ξ1, ξ2,k0, J,K,Jp,Kp},
		{ϵ1,ϵ2}=nindex^2 ; k0=2*Pi/λ;
    	δ1=Sqrt[Abs[ϵ1 -n^2]]; δ2=Sqrt[Abs[ϵ2-n^2] ];
		ξ1=δ1*ρ *k0 ; ξ2=δ2*ρ* k0 ;
		J=BesselJ[l,ξ1] ; K=BesselK[l, ξ2];
		Jp=0.5*(BesselJ[l-1,ξ1] -BesselJ[l+1,ξ1]);
		Kp=-0.5*(BesselK[l-1,ξ2] + BesselK[l+1,ξ2]) ;
     (δ1 J Kp+ δ2 K Jp) (δ1  ϵ2 J Kp+δ2 ϵ1 K Jp)
	-(( (n^2))/( (k0^2) (ρ^2) ((δ1)^2) ((δ2)^2)))(l^2) (((ϵ1-ϵ2))^2) K^2 J^2   ]    ;

Coeff[nindex_List,λ_,ρ_,n_,l_]:=Module[{ϵ1,ϵ2,δ1, δ2,ξ1, ξ2,k0, J,K,Jp,Kp},
		{ϵ1,ϵ2}=nindex^2 ; k0=2*Pi/λ;
    	δ1=Sqrt[Abs[ϵ1 -n^2]]; δ2=Sqrt[Abs[ϵ2-n^2] ];
		ξ1=δ1*ρ *k0 ; ξ2=δ2*ρ* k0 ;
		J=BesselJ[l,ξ1] ; K=BesselK[l, ξ2];
		Jp=0.5*(BesselJ[l-1,ξ1] -BesselJ[l+1,ξ1]);
		Kp=-0.5*(BesselK[l-1,ξ2] + BesselK[l+1,ξ2]) ;
		{ 1, ((I l J K n (ϵ1-ϵ2))/(ρ k0 δ1 δ2 (δ2  K   Jp+δ1 J Kp))),
			-((J δ1)/(K δ2)) , 
                             -((I l (J^2)  n (ϵ1-ϵ2))/(ρ k0 ((δ2)^2) (δ2 K Jp+δ1 J Kp)))} 
 ]

FieldFiber[coor_List, nindex_List,λ_,rad_,n_,l_,coef_List]:=Module[{ρ, ϕ, ϵ1,ϵ2,q1, q2,ξ1, ξ2,k0, J,K,Jp,Kp, field},
		{ρ, ϕ}=coor ;
			{ϵ1,ϵ2}=nindex^2 ; k0=2*Pi/λ;
    	q1=Sqrt[Abs[ϵ1 -n^2]]; q2=Sqrt[Abs[ϵ2-n^2] ];
		ξ1=q1*ρ *k0 ; ξ2=q2*ρ* k0 ;
		J=BesselJ[l,ξ1] ; K=BesselK[l, ξ2];
		Jp=0.5*(BesselJ[l-1,ξ1] -BesselJ[l+1,ξ1]);
		Kp=-0.5*(BesselK[l-1,ξ2] + BesselK[l+1,ξ2]) ;
		field=If[ Abs[ρ]<=rad,
			 { J  q1
				, coef[[2]]*  J  q1
				, I n Jp -coef[[2]]*l *J/ξ1
				, -l *n J/ξ1 -coef[[2]] * I *Jp
				, l *ϵ1  J/ξ1+coef[[2]]*I n Jp
				, I ϵ1 Jp-coef[[2]]*l*J*n/ξ1}
			,{-coef[[3]]*K q2
				,-coef[[4]]*K q2
				,coef[[3]]*I n Kp-coef[[4]]*l K/ξ2
				,-coef[[3]]*l n K/ξ2-coef[[4]]*I *Kp
                  ,coef[[3]]*l K ϵ2/ξ2+coef[[4]] I n Kp
	,coef[[3]] I ϵ2 Kp-coef[[4]]l K n/ξ2}]  
;
		{field * (E^(I ϕ l)),field *(E^(-I ϕ l))}
		          ]  ;
(*  Note: the field output six components are: (({{(Subscript[E, z]) }, {(Subscript[H, z])}, {(Subscript[E, ρ])}, {(Subscript[E, ϕ])}, {(Subscript[H, ρ])}, {(Subscript[H, ϕ])}}))    *)

5.6.4 Plot 1  l=1

Indexprof={1.48,1.45};l=1;λ=1.55; rad=5.5/2;
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
Plot[Log[Abs[CharG[Indexprof,λ, rad,x,l]]] , {x,Indexprof[[2]]+dn,Indexprof[[1]]-dn}
		,PlotStyle->{{RGBColor[1,0.,0.]},{RGBColor[1,0,0]},{RGBColor[0,0,1]}} ]

segment={1.472,1.474};l=1;
neff= x /. FindRoot[CharG[Indexprof,λ, rad,x,l] , {x,segment[[1]],segment[[2]]}
	          ] 
coef=Coeff[Indexprof,λ,rad,neff,l]

Plot E

ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,5*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,3]] ]^2
+5*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,4]] ]^2
		},
    {ρ , 0.01, 7},
    {ϕ, 0, 2 Pi }, PlotRange->All]

A concept is mode size. It is NOT a measurement of NOT some fundamental property, but a definition to describe roughly how large the spatial extent of the field is. The definition is:

5.6.5 Plot 2; l=0

Indexprof={1.48,1.45};l=0;λ=1.55; rad=5.5/2;
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
Plot[Log[Abs[CharG[Indexprof,λ, rad,x,l]]] , {x,Indexprof[[2]]+dn,Indexprof[[1]]-dn}
		,PlotStyle->{{RGBColor[1,0.,0.]},{RGBColor[1,0,0]},{RGBColor[0,0,1]}} 
,PlotRange-> {-20,0}]

segment={1.455,1.46};
		Plot[Log[Abs[CharG[Indexprof,λ, rad,x,l]]] , {x,segment[[1]],segment[[2]]}
		,PlotStyle->{{RGBColor[1,0.,0.]},{RGBColor[1,0,0]},{RGBColor[0,0,1]}} , AspectRatio->0.4 ]

segment={1.4575,1.4581};
neff= x /. FindRoot[CharG[Indexprof,λ, rad,x,l] , {x,segment[[1]],segment[[2]]}
	          ]
coef=Coeff[Indexprof,λ,rad,neff,l]

Plot E

grEρ=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,50*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,3]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Eρ|)^2 for l=1",DisplayFunction-> Identity];
grEϕ=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,50*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,4]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Eϕ|)^2 for l=1",DisplayFunction-> Identity] ;
grEz=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,50*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,1]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Ez|)^2 for l=1",DisplayFunction-> Identity] ;
Show[GraphicsRow[{grEρ,grEϕ,grEz}], DisplayFunction->$DisplayFunction, PlotLabel-> "This is (|E|)^2 for l=0"]

Plot H

grHρ=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,50*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,5]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Hρ|)^2 for l=1",DisplayFunction-> Identity];
grHϕ=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,50*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,6]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Hϕ|)^2 for l=1",DisplayFunction-> Identity] ;
grHz=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,50*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,2]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Hz|)^2 for l=1",DisplayFunction-> Identity] ;
Show[GraphicsRow[{grHρ,grHϕ,grHz}], DisplayFunction->$DisplayFunction, PlotLabel-> "This is (|H|)^2 for l=0"]

5.6.6 Plot 3; l=2

Indexprof={1.48,1.45};l=2;λ=1.55; rad=12./2;
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
Plot[Log[Abs[CharG[Indexprof,λ, rad,x,l]]] , {x,Indexprof[[2]]+dn,Indexprof[[1]]-dn}
		,PlotStyle->{{RGBColor[1,0.,0.]},{RGBColor[1,0,0]},{RGBColor[0,0,1]}} ]

segment={1.47,1.475};
neff= x /. FindRoot[CharG[Indexprof,λ, rad,x,l] , {x,segment[[1]],segment[[2]]}
	          ]
coef=Coeff[Indexprof,λ,rad,neff,l]

Plot E field

grEρ=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,5*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,3]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Eρ|)^2 for l=2",DisplayFunction-> Identity];
grEϕ=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,5*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,4]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Eϕ|)^2 for l=2",DisplayFunction-> Identity] ;
grEz=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,5*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,1]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Ez|)^2 for l=2",DisplayFunction-> Identity] ;
Show[GraphicsRow[{grEρ,grEϕ,grEz}], DisplayFunction->$DisplayFunction, PlotLabel-> "This is (|E|)^2 for l=2"]

Plot H field

grHρ=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,5*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,5]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Hρ|)^2 for l=2",DisplayFunction-> Identity];
grHϕ=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,5*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,6]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Hϕ|)^2 for l=2",DisplayFunction-> Identity] ;
grHz=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,5*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,2]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Hz|)^2 for l=2",DisplayFunction-> Identity] ;
Show[GraphicsRow[{grHρ,grHϕ,grHz}], DisplayFunction->$DisplayFunction, PlotLabel-> "This is (|H|)^2 for l=2"]

segment={1.455,1.460};
neff= x /. FindRoot[CharG[Indexprof,λ, rad,x,l] , {x,segment[[1]],segment[[2]]}
	          ]
coef=Coeff[Indexprof,λ,rad,neff,l]

Plot E field

grEρ=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,5*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,3]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Eρ|)^2 for l=2",DisplayFunction-> Identity, Mesh-> None];
grEϕ=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,5*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,4]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Eϕ|)^2 for l=2",DisplayFunction-> Identity, Mesh-> None] ;
grEz=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,5*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,1]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Ez|)^2 for l=2",DisplayFunction-> Identity, Mesh-> None] ;
Show[GraphicsRow[{grEρ,grEϕ,grEz}], DisplayFunction->$DisplayFunction, PlotLabel-> "This is (|E|)^2 for l=2"]

5.6.7 Plot 4;

Indexprof={1.48,1.41};l=3;λ=1.55; rad=5.5/2;
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
Plot[Log[Abs[CharG[Indexprof,λ, rad,x,l]]] , {x,Indexprof[[2]]+dn,Indexprof[[1]]-dn}
		,PlotStyle->{{RGBColor[1,0.,0.]},{RGBColor[1,0,0]},{RGBColor[0,0,1]}} ]

segment={1.42,1.44};
neff= x /. FindRoot[CharG[Indexprof,λ, rad,x,l] , {x,segment[[1]],segment[[2]]}
	          ]
coef=Coeff[Indexprof,λ,rad,neff,l]

Plot E field

grEρ=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,10*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,3]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Eρ|)^2 for l=2",DisplayFunction-> Identity, Mesh-> None];
grEϕ=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,10*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,4]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Eϕ|)^2 for l=2",DisplayFunction-> Identity, Mesh-> None] ;
grEz=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,50*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,1]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Ez|)^2 for l=2",DisplayFunction-> Identity, Mesh-> None] ;
Show[GraphicsRow[{grEρ,grEϕ,grEz}], DisplayFunction->$DisplayFunction, PlotLabel-> "This is (|E|)^2 for l=2"]

5.6.8 Plot 5

Indexprof={1.48,1.3};l=4;λ=1.55; rad=12./2;
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
Plot[Log[Abs[CharG[Indexprof,λ, rad,x,l]]] , {x,Indexprof[[2]]+dn,Indexprof[[1]]-dn}
		,PlotStyle->{{RGBColor[1,0.,0.]},{RGBColor[1,0,0]},{RGBColor[0,0,1]}} ]

Plot[Log[Abs[CharG[Indexprof,λ, rad,x,l]]] , {x,1.33,1.35}
		,PlotStyle->{{RGBColor[1,0.,0.]},{RGBColor[1,0,0]},{RGBColor[0,0,1]}} ]

segment={1.335,1.345};
neff= x /. FindRoot[CharG[Indexprof,λ, rad,x,l] , {x,segment[[1]],segment[[2]]}
, AccuracyGoal->20          ]
coef=Coeff[Indexprof,λ,rad,neff,l]

Plot E field

grEρ=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,10*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,3]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Eρ|)^2 for l=4",DisplayFunction-> Identity, Mesh-> None];
grEϕ=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,10*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,4]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Eϕ|)^2 for l=4",DisplayFunction-> Identity, Mesh-> None] ;
grEz=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,50*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,1]] ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Ez|)^2 for l=4",DisplayFunction-> Identity, Mesh-> None] ;
Show[GraphicsRow[{grEρ,grEϕ,grEz}], DisplayFunction->$DisplayFunction, PlotLabel-> "This is (|E|)^2 for l=4"]

grEρ=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,25*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,3]] ]
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Eρ|)^2 for l=4",DisplayFunction-> Identity, Mesh-> None];
grEϕ=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,25*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,4]] ]
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Eϕ|)^2 for l=4",DisplayFunction-> Identity, Mesh-> None] ;
grEz=ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,25*Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coef][[1,1]] ]
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, PlotLabel-> "This is (|Ez|)^2 for l=4",DisplayFunction-> Identity, Mesh-> None] ;
Show[GraphicsRow[{grEρ,grEϕ,grEz}], DisplayFunction->$DisplayFunction, PlotLabel-> "This is (|E|)^2 for l=4"]

6.4 Field profile illustration

Source code - Demo step index fiber - field illustration

Manipulate[
If[INIT== 1,
{off,on}={Button["OFF",BaseStyle->{16, Bold,White,FontFamily-> "Arial"}
,Background->Red]
,Button["ON",BaseStyle->{16, Bold,White,FontFamily-> "Arial"}
,Background->Hue[0.36,1,0.4]]} ;
CharG[nindex_List,λ_,ρ_,n_,l_]:=Module[{ϵ1,ϵ2,δ1, δ2,ξ1, ξ2,k0, J,K,Jp,Kp},
		{ϵ1,ϵ2}=nindex^2 ; k0=2*Pi/λ;
    	δ1=Sqrt[Abs[ϵ1 -n^2]]; δ2=Sqrt[Abs[ϵ2-n^2] ];
		ξ1=δ1*ρ *k0 ; ξ2=δ2*ρ* k0 ;
		J=BesselJ[l,ξ1] ; K=BesselK[l, ξ2];
		Jp=0.5*(BesselJ[l-1,ξ1] -BesselJ[l+1,ξ1]);
		Kp=-0.5*(BesselK[l-1,ξ2] + BesselK[l+1,ξ2]) ;
     (δ1 J Kp+ δ2 K Jp) (δ1  ϵ2 J Kp+δ2 ϵ1 K Jp)
	-(( (n^2))/( (k0^2) (ρ^2) ((δ1)^2) ((δ2)^2)))(l^2) (((ϵ1-ϵ2))^2) K^2 J^2   ] ; 

Coeff[nindex_List,λ_,ρ_,n_,l_]:=Module[{ϵ1,ϵ2,q1, q2,ξ1, ξ2,k0, J,K,Jp,Kp},
		{ϵ1,ϵ2}=nindex^2 ; k0=2*Pi/λ;
    	q1=Sqrt[Abs[ϵ1 -n^2]]; q2=Sqrt[Abs[ϵ2-n^2] ];
		ξ1=q1*ρ *k0 ; ξ2=q2*ρ* k0 ;
		J=BesselJ[l,ξ1] ; K=BesselK[l, ξ2];
		Jp=0.5*(BesselJ[l-1,ξ1] -BesselJ[l+1,ξ1]);
		Kp=-0.5*(BesselK[l-1,ξ2] + BesselK[l+1,ξ2]) ;
		{ 1, ((I l J K n (ϵ1-ϵ2))/(ρ k0 q1 q2 (q2 K   Jp+q1 J Kp))),
			-((J q1)/(K q2)) , 
                             -((I l (J^2)  n (ϵ1-ϵ2))/(ρ k0 (q2^2) (q2 K Jp+q1 J Kp)))} 
 ] ;

FieldFiber[coor_List, nindex_List,λ_,rad_,n_,l_,coef_List]:=Module[{ρ, ϕ, ϵ1,ϵ2,q1, q2,ξ1, ξ2,k0, J,K,Jp,Kp, field},
		{ρ, ϕ}=coor ;
			{ϵ1,ϵ2}=nindex^2 ; k0=2*Pi/λ;
    	q1=Sqrt[Abs[ϵ1 -n^2]]; q2=Sqrt[Abs[ϵ2-n^2] ];
		ξ1=q1*ρ *k0 ; ξ2=q2*ρ* k0 ;
		J=BesselJ[l,ξ1] ; K=BesselK[l, ξ2];
		Jp=0.5*(BesselJ[l-1,ξ1] -BesselJ[l+1,ξ1]);
		Kp=-0.5*(BesselK[l-1,ξ2] + BesselK[l+1,ξ2]) ;
		field=If[ Abs[ρ]<=rad,
			{ J  q1
				, coef[[2]]*  J  q1
				, I n Jp -coef[[2]]*l *J/ξ1
				, -l *n J/ξ1 -coef[[2]] * I *Jp
				, l *ϵ1  J/ξ1+coef[[2]]*I n Jp
				, I ϵ1 Jp-coef[[2]]*l*J*n/ξ1}
			,{-coef[[3]]*K q2
				,-coef[[4]]*K q2
				,coef[[3]]*I n Kp-coef[[4]]*l K/ξ2
				,-coef[[3]]*l n K/ξ2-coef[[4]]*I *Kp
                  ,coef[[3]]*l K ϵ2/ξ2+coef[[4]] I n Kp
	,coef[[3]] I ϵ2 Kp-coef[[4]]l K n/ξ2}]  
;
		I  {field * (E^(I ϕ l)),field *(E^(-I ϕ l))}
		          ];
ModeSearch[{on_,off_},{Indexprof_,λ_,rad_,l_},segment_]:=
Module[{neff,coefp,coefm}
,neff=SetPrecision[ x /. FindRoot[CharG[Indexprof,λ, rad,x,l] , {x,segment[[1]],segment[[2]]}
,AccuracyGoal-> 16], 8]  ; 

coefp=Coeff[Indexprof,λ,rad,neff,l] ; coefm=Coeff[Indexprof,λ,rad,neff,-l] ;

{neff,{coefp,coefm}}

];
SelectAppend[allmode_,{{l_,neff_},{coefp_,coefm_}}]:=Module[{u,am},
am=allmode;
u=Select[allmode,#[[1,1]]== l &];
If[Length[u]==0,AppendTo[am,{{l,neff},{coefp,coefm}}]
,If[Abs[u[[1]][[1,2]]-neff]≥ (10^(-6)),AppendTo[am,{{l,neff},{coefp,coefm}}], ]
] ;am
];

INIT=0 
,TextCell["Not ready"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]] 
;


Manipulate[plotready=0 ;
Column[{Button["Mode Search for Step Index Fiber"
,BaseStyle->{20, Bold,Black,FontFamily-> "Arial"}
,Background->Hue[0.5,0.7,1] ]

,If[RUN== 1,nclad=1.45;allmode={};
Manipulate[
If[λ≠ λG || ncore≠ ncoreG || rad≠ radG,allmode={}, ];
ncoreG=ncore ;
Indexprof={ncore,1.45};
λG=λ ; radG=rad;

dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
Row[{Plot[Log[Abs[CharG[Indexprof,λ, rad,x,l]]] , {x,Indexprof[[2]]+dn,Indexprof[[1]]-dn}
		,PlotStyle->{{RGBColor[1,0.,0.]},{RGBColor[1,0,0]},{RGBColor[0,0,1]}} 
, PlotRange->{{Min[{nplow,nphi}],Max[{nplow,nphi}]+0.001}
, All}
,Frame-> True, GridLines->Automatic
,FrameLabel->{"Effective mode index","Log CharFunction"}
,LabelStyle->{14, FontFamily-> "Arial"}
,ImageSize->{300,200} ]
,Spacer[25]
,
Manipulate[
If[TrueQ[SEARCH== on]
,plotready=0;
{neff,{coefp,coefm}}=ModeSearch[{on,off},{Indexprof,λ,rad,l},{nplow,nphi}];
allmode=SelectAppend[allmode,{{l,neff},{coefp,coefm}}];
allmode=Union[allmode] ;
{Indexprofp,λp,radp}={{ncoreG,1.45},λG,radG} ;
SEARCH=off 
,] ;
Manipulate[If[MDELETE==1,
If[Length[allmode]>0 && NumberQ[MSELECT] && MSELECT>0,allmode=Drop[allmode,{MSELECT}], ] ;
MDELETE=0, ];
Grid[Join[{{"l","neff","Deletion"}},Table[Append[allmode[[i,1]]
,If[MSELECT==i,
Button["X",BaseStyle->{Red,14,Bold,FontFamily-> "Arial"}],""]],{i,1,Length[allmode]}]]
,BaseStyle-> {16, FontFamily-> "TimesNewRoman"},Background->Hue[0.9,0.1,1]
,Frame-> All,FrameStyle->Directive[Black]]
,{MSELECT,Join[{"no selection"},Range[Length[allmode]]]}
,{{MDELETE,0,Button["Click to delete->",BaseStyle->{Red,16,FontFamily-> "Calibri"},Background->GrayLevel[0.9]]},{0,1},ControlPlacement->Bottom}]

,{{SEARCH,off},{on,off},ControlType-> Toggler}] }]
,Grid[{
{Control[{{ncore,1.48,Button["core index",BaseStyle->{16, Bold,Black,FontFamily-> "Calibri"},Background->Hue[0.3,1,1,1]]},1.46,1.52}]
,Control[{{rad,2.75,Button["core radius (μm)",BaseStyle->{16, Bold,Black,FontFamily-> "Calibri"},Background->Hue[0.45,1,1,1]]},2.5,4}]}
,{Control[{{l,0,Button["angular mode",BaseStyle->{16, Bold,Black,FontFamily-> "Calibri"},Background->Hue[0.15,1,1,1]]},0,5,1,ControlType-> SetterBar}]
,
Control[{{λ,1.55,Button["wavelength λ (μm)",BaseStyle->{16, Bold,Black,FontFamily-> "Calibri"},Background->Hue[0.036,0.5,1,1]]},0.5,3}]}}]
,Control[{{nplow,nclad+0.0005,"lower limit"},nclad+0.0005,ncore-0.001,0.000005,ControlPlacement-> Bottom}]
,Control[{{nphi,ncore-0.0005,"upper limit"},nclad+0.001,ncore-0.0005,0.000005,ControlPlacement-> Bottom}]
]
,TextCell["Not running"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]]

,Button["Field plot"
,BaseStyle->{20, Bold,Black,FontFamily-> "Arial"}
,Background->Hue[0.15,1,1] ]
,Manipulate[
If[TrueQ[FIELD==on],

Manipulate[If[MODESELECT==1,plotready=0 ;
If[Length[allmode]>0 && NumberQ[FIELDSELECT] && FIELDSELECT>0
,
{{lp,neffp},{coefp,coefm}}=allmode[[FIELDSELECT]];
{ρ,ϕ}= {0.1,0.05} ;
If[EH=="E",
grEρ=ParametricPlot3D[{ρ Cos[ϕ], ρ Sin[ϕ], Abs[ FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,lp,coefp][[1,3]]
+FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,-lp,coefm][[1,3]] ]^2}
,{ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All
, PlotLabel->StringJoin["This is Eρ^2 for l=±",ToString[lp]]
, Mesh-> None
,BoxRatios->{1,1,0.25},ImageSize-> {150,150}];
grEϕ=ParametricPlot3D[{ρ Cos[ϕ], ρ Sin[ϕ],Abs[ FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,lp,coefp][[1,4]]
+FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,-lp,coefm][[1,4]] ]^2},{ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All
, PlotLabel-> StringJoin["This is Eϕ^2 for l=±",ToString[lp]]
 , Mesh-> None
,BoxRatios->{1,1,0.25},ImageSize-> {150,150}] ;
grEz=ParametricPlot3D[{ρ Cos[ϕ], ρ Sin[ϕ],Abs[ FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,lp,coefp][[1,1]]
+FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,-lp,coefm][[1,1]] ]^2}
,{ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All
, PlotLabel->StringJoin["This is Ez^2 for l=±",ToString[lp]]
, Mesh-> None
,BoxRatios->{1,1,0.25},ImageSize-> {150,150}] ;
grET=ParametricPlot3D[{ρ Cos[ϕ], ρ Sin[ϕ],Re[ FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,lp,coefp][[1,3]] ]^2
+Re[FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,lp,coefp][[1,4]]]^2}
,{ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange-> All
, PlotLabel->StringJoin["This is (|ETransverse|)^2 for l=+",ToString[lp]], Mesh-> None
,BoxRatios->{1,1,0.25},ImageSize-> {150,150}] ;
plotready=1
,
grEρ=ParametricPlot3D[{ρ Cos[ϕ], ρ Sin[ϕ], Abs[ FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,lp,coefp][[1,5]]
+FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,-lp,coefm][[1,5]] ]^2}
,{ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All
, PlotLabel->StringJoin["This is Hρ^2 for l=±",ToString[lp]]
, Mesh-> None
,BoxRatios->{1,1,0.25},ImageSize-> {150,150}];
grEϕ=ParametricPlot3D[{ρ Cos[ϕ], ρ Sin[ϕ],Abs[ FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,lp,coefp][[1,6]]
+FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,-lp,coefm][[1,6]] ]^2},{ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All
, PlotLabel-> StringJoin["This is Hϕ^2 for l=±",ToString[lp]]
 , Mesh-> None
,BoxRatios->{1,1,0.25},ImageSize-> {150,150}] ;
grEz=ParametricPlot3D[{ρ Cos[ϕ], ρ Sin[ϕ],Abs[ FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,lp,coefp][[1,1]]
+FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,-lp,coefm][[1,1]] ]^2}
,{ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All
, PlotLabel->StringJoin["This is Hz^2 for l=±",ToString[lp]]
, Mesh-> None
,BoxRatios->{1,1,0.25},ImageSize-> {150,150}] ;
grET=ParametricPlot3D[{ρ Cos[ϕ], ρ Sin[ϕ],Re[ FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,lp,coefp][[1,5]] ]^2
+Re[FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,lp,coefp][[1,6]]]^2}
,{ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange-> All
, PlotLabel->StringJoin["This is (|HTransverse|)^2 for l=+",ToString[lp]], Mesh-> None
,BoxRatios->{1,1,0.25},ImageSize-> {150,150}] ;
plotready=1
]
, plotready=0] ;
MODESELECT=0, ];
Grid[{{Grid[Join[{{"l","neff","Selection"}},Table[Append[allmode[[i,1]]
,If[FIELDSELECT==i,
Button["X",BaseStyle->{Hue[0.35,1,0.5],14,Bold,FontFamily-> "Arial"}],""]],{i,1,Length[allmode]}]]
,BaseStyle-> {16, FontFamily-> "TimesNewRoman"},Background->Hue[0.9,0.1,1]
,Frame-> All,FrameStyle->Directive[Black]]

,

}
,{If[TrueQ[plotready== 1],
Show[GraphicsGrid[{{grEρ,grEϕ},{grEz,grET}}],  PlotLabel->StringJoin["E or H field for l=",ToString[l]]], ]
,If[TrueQ[plotready== 1],
ρL=Table[ρ,{ρ,0.3,1.2 radp,1.2 radp/12}] ;
ϕL=Table[ϕ,{ϕ,0,2 π, π/12}] ;
nρ=Length[ρL];mϕ=Length[ϕL];
FHp=Table[FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,lp,coefp],{ρ,0.3,1.2 radp,1.2 radp/12},{ϕ,0,2 π, π/12}] ;
FHm=Table[FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,-lp,coefm],{ρ,0.3,1.2 radp,1.2 radp/12},{ϕ,0,2 π, π/12}] ;
Manipulate[
If[EH=="E", 
Graphics[Table[{(* Hue[
0.7*(Abs[FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,lp,coefp][[1,3]] ]^2
+Abs[FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,lp,coefp][[1,4]]]^2 )
,1,1]*)
Hue[0,0.8,0.8]

,Arrow[{{
ρL[[i]] Cos[ϕL[[j]]]-s Re[ (E^(-I 2 π t)) FHp[[i,j]][[1,3]] ] Cos[ϕL[[j]]]
                +s Re[ (E^(-I 2 π t)) FHp[[i,j]][[1,4]] ] Sin[ϕL[[j]]]
, ρL[[i]] Sin[ϕL[[j]]]-s Re[ (E^(-I 2 π t)) FHp[[i,j]][[1,3]] ] Sin[ϕL[[j]]]
                  -s Re[ (E^(-I 2 π t)) FHp[[i,j]][[1,4]] ] Cos[ϕL[[j]]] }
,{ρL[[i]] Cos[ϕL[[j]]]+s Re[ (E^(-I 2 π t)) FHp[[i,j]][[1,3]] ] Cos[ϕL[[j]]]
                -s Re[ (E^(-I 2 π t)) FHp[[i,j]][[1,4]] ] Sin[ϕL[[j]]]
, ρL[[i]] Sin[ϕL[[j]]]+s Re[ (E^(-I 2 π t)) FHp[[i,j]][[1,3]] ] Sin[ϕL[[j]]]
                  +s Re[ (E^(-I 2 π t)) FHp[[i,j]][[1,4]] ] Cos[ϕL[[j]]] }
}]}  ,{i,1,nρ},{j,1,mϕ}] 
,PlotRange->{{-1.5 radp,1.5 radp},{-1.5 radp,1.5 radp}}
,ImageSize-> {200,200}]
,
Graphics[Table[{(* Hue[
0.7*(Abs[FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,lp,coefp][[1,3]] ]^2
+Abs[FieldFiber[{ρ,ϕ},Indexprofp,λp,radp,neffp,lp,coefp][[1,4]]]^2 )
,1,1]*)
Hue[0.65,1,1]

,Arrow[{{
ρL[[i]] Cos[ϕL[[j]]]-s Re[ (E^(-I 2 π t)) FHp[[i,j]][[1,5]] ] Cos[ϕL[[j]]]
                +s Re[ (E^(-I 2 π t)) FHp[[i,j]][[1,6]] ] Sin[ϕL[[j]]]
, ρL[[i]] Sin[ϕL[[j]]]-s Re[ (E^(-I 2 π t)) FHp[[i,j]][[1,5]] ] Sin[ϕL[[j]]]
                  -s Re[ (E^(-I 2 π t)) FHp[[i,j]][[1,6]] ] Cos[ϕL[[j]]] }
,{ρL[[i]] Cos[ϕL[[j]]]+s Re[ (E^(-I 2 π t)) FHp[[i,j]][[1,5]] ] Cos[ϕL[[j]]]
                -s Re[ (E^(-I 2 π t)) FHp[[i,j]][[1,6]] ] Sin[ϕL[[j]]]
, ρL[[i]] Sin[ϕL[[j]]]+s Re[ (E^(-I 2 π t)) FHp[[i,j]][[1,5]] ] Sin[ϕL[[j]]]
                  +s Re[ (E^(-I 2 π t)) FHp[[i,j]][[1,6]] ] Cos[ϕL[[j]]] }
}]}  ,{i,1,nρ},{j,1,mϕ}] 
,PlotRange->{{-1.5 radp,1.5 radp},{-1.5 radp,1.5 radp}}
,ImageSize-> {200,200}]
],
{ {t,0,Button["time",BaseStyle->{16, Bold,FontFamily-> "Calibri"}
,Background->Yellow]},0,1,0.02, DisplayAllSteps-> True}

,{{s,0.3,Button["vector scale",BaseStyle->{16, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.45,1,0.7]]},0.1,1.5}], ]}}]

,Row[{Control[{FIELDSELECT,Join[{"no selection"},Range[Length[allmode]]], ControlType-> SetterBar}],Spacer[10]
,Control[{{MODESELECT,0,Button["Click to select->",BaseStyle->{Hue[0.35,1,0.5],16,FontFamily-> "Calibri"},Background->GrayLevel[0.9]]},{0,1}}],Spacer[30]
,Control[{{EH,"E",Button["E or H field",BaseStyle->{White,16,FontFamily-> "Calibri"},Background->Hue[0.8,1,0.5]]},{"E","H"},ControlPlacement->Top}]
}]
]
 ,TextCell["Not running"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]]
,{{FIELD,off},{on,off},ControlType-> Toggler}]


}]
,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}}]
}]
]
,Row[{Button["Click once to initialize",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Red],Control[{{INIT,0,"->"},{0,1}}]
}]
]

Demo step index fiber - field illustration

7. Other types of fibers

7.1 Discussion

There are many ways to make fibers and waveguides, and in practice, there are other types of fibers.

7.2 Summary key points and issues in fiber design

1- Optical fibers are waveguides with cylindrical geometry although this is NOT necessary and some special-application fibers have elliptical cross section (polarization-preservation fibers). (Note that planar waveguide has rectangular cross section.)
2- A fiber can be designed to be single-mode or multimode, depending on the index profile and the size of the core. The larger the index difference and the larger the diameter, the more modes a fiber may have. Multi-mode fibers are designed with many times larger diameter than single-mode fibers. Each mode has an effective refractive index.
3- A fiber is designed for a certain wavelength range. For a given fiber, it can be multi-mode with shorter wavelengths, single-mode for certain wavelength range and may not have any mode for wavelength longer than certain value called cut-off wavelength (Subscript[λ, c])
4- The modes depend on the detailed index profile of the fiber, which is a design parameter so that specific properties can be optimized.
	- Mode size is a measure (defined, not fundamental) to describe roughly how large the field spatial profile is (it is roughly around the diameter of the core in many cases, but NOT all cases). [Note: Mode spatial profile is not important for the detection, which depends on only the total power.]
	- A very important property is dispersion. 
5- Dispersion includes: material dispersion, (modal) chromatic dispersion, intermodal dispersion (if multi-mode), and polarization dispersion. This topic is important and will be treated in details in subsequent chapters. It determines how much information that can be sent through a fiber (its bandwidth x reach)
6- Designs are done with numerical methods (e. g. software package) using basic principles of how the various properties vary vs. design parameters.

7.3 Graded index fiber

Graded index fiber has an index profile that is a function of radial distance, usually in the core only, although the grading can be anywhere, core and cladding.
Example:                             n[ρ]=(Subscript[n, 1])(Sqrt[1-2Δ ((((ρ/a)))^γ)])
which is for Δ<< 1
In fact, we can approximate:
                                               n[ρ]≈(Subscript[n, 1])(1-Δ ((((ρ/a)))^γ))    
The square root is to represent the dielectric constant:
                                          (Subscript[ϵ, 1])[ρ]= ((n[r])^2)=((Subscript[n, 1])^2) (1-2Δ ((((r/a)))^γ))=(Subscript[ϵ, 1;0]) (1-2Δ ((((r/a)))^γ))

a=1. ; Δ= 0.02  ; n1=1.5  ;
Plot[Table[ If[ Abs[ρ]≤ a, n1(Sqrt[1-2Δ (((((Abs[ρ])/a)))^γ)])  ,n1(Sqrt[1-2Δ])],
 {γ,1,4,1}],{ρ,-2,2}, PlotRange-> {1.45,1.5}]

A quadratic index profile:   n[ρ]=(Subscript[n, 1])(1-Δ ((((ρ/a)))^2)) acts like an infinite lens. (See figure above) 
A characteristic of graded index fiber is that the higher the order of a mode, the lower the index profile that it sees, hence although it has more delay in its path, it is compensated by traveling faster than the lower order modes. This leads to less inter-modal dispersion. Thus, for multi-mode fibers, graded index can be an advantage.
Generally, it is impractical to solve the wave equation analytically, it is done numerically.
Some special case like the quadratic profile has some reasonable approximations.

ρ=.  ;

DSolve[ ρ G'[ρ] + (ρ^2) G''[ρ] +(((k0^2)  ϵ1  (1-2Δ (ρ^2)) - (β^2))(ρ^2) -(l^2) )G[ρ ]== 0 ,
	                             G[ρ ], ρ ]

Notice that it has analytic solution expressed in terms of Hypergeometric function U and Laguerre polynomial.

It can be aproximated that the effective index of various modes (l, m) is:
                                    (Subscript[n, eff])[l,m]=(Subscript[n, 1])(Sqrt[((1-2(Sqrt[2 Δ]))/((Subscript[n, 1]) (Subscript[k, 0]) a))(2 m + l +1)])
where l is the azimuthal number from (e^(i l ϕ)) and m is the radial mode number.

7.4 Holey fibers or PCF (photonic crystal fiber) and PBG (photonic bandgap) fiber

7.4.1 Concept introduction

The general concept of dielectric waveguide is to have a high-index core relative to low index cladding. An innovative way of doing this is to “drill” or remove materials of a uniform medium unevenly to leave a core that is “denser” (hence  higher index) than the surrounding (hence lower index). Holey fiber is invented as shown below

Although the calculation is complicated, it can be shown that waveguiding occurs. The hexagon shape is the tightest 2-D periodic lattice packing of the “holes”, but other arrangement, like a square lattice would work also, although it is less well matched to the cylindrical geometry. The essence of the waveguiding mechanism is still of the fact that the center region index is slightly higher than the surrounding, similarly to the dielectric fibers discussed above.

It turns out that there was a very similar concept called photonic bandgap materials or photonic crystals. This concept is related to the problem of electrons in crystals.
In a periodic structure like a crystal, quantum mechanics shows that the electron energy and wavefunction have certain unique properties. In particular, the energy forms structure knows as energy bands with “gaps” between them. In other words, there are “forbidden” energy range in which the electron does not have any state. It means that at certain frequency (related to energy with the Planck’s constant: E=ℏω), the electron doesn’t have any propagating wave solution.

There is a similar situation with optics in periodic media.  For one-dimensional periodic thin films, there is a structure known as Bragg reflector. It is a periodic structure in which certain range of frequencies will be strongly reflected (no propagating solution exists) while others can transmit freely. These frequency ranges are similar to energy “band gap”, hence the expression “photonic bandgap”.

The concept is applicable to 2-D or 3-D periodic optical materials, and hence the term photonic crystals, which can be applied to make waveguide fibers. However, here, there is a subtle difference with the waveguiding mechanism in the holey fibers discussed above. The light is confined by a surrounding 2D lattice that is designed such that in the transverse dimension, its frequency is on the “band gap” and hence cannot escape (transmitted) through the wall.
Another way to think of it is like a mirror waveguide.
Nevertheless, all of these are also described as “microstructured fibers” a generic terms for various caterogies.

www.ipht-jena.de

7.4.2 Various examples of microstructured fibers

From spie.org

www.ipht-jena.de

Holey fibers are example of microstructure fibers. Above is a comparison of different fibers 
For the PBG fibers, there is no need for index guiding and the core can be hollow as illustrated below.

park.itc.u-tokyo.ac.jp

www.math.utwente.nl/aamp/ex_mathopt.html

7.4.3 Summary of advanced types of fiber

Summary:
1- Holey fibers and PBG fibers or PCF offer other approaches of making waveguides that can provide additional flexibilities in designing fiber properties
2- Various aspects of the microstructures: hole sizes, patterns, lattice configuration can be use as design parameters to obtain a diverse range of properties 
3- Some of the properties:	
	- A wide wavelength band
	- A wide range of mode size (from very small to very large)
	- Dispersion properties
	- Up to large NA, even for multi-mode fibers.
	- Asymmetry can be used to make polarization maintaining fiber or to design specific polarization property
	- A number of nonlinear optics applications
In short, they offer complementary approaches to obtain waveguide fibers beside the dielectric waveguide fibers.

8. Optical fiber and cables: manufacturing

Please see the below

http://www.youtube.com/watch?v=H1Yj-hg2UiU

http://www.youtube.com/watch?v=liKOYbgIC_c

http://www.youtube.com/watch?v=7tsF3mSpqX8

http://www.youtube.com/watch?v=Ncf4jRUQclI

http://www.youtube.com/watch?v=0RYiBqHgOl0

Remember cymatics:

http://www.youtube.com/watch?v=GtiSCBXbHAg

http://www.youtube.com/watch?v=CGiiSlMFFlI