Cylindrical geometry waveguide - Optical Fiber

1. Introduction and review

Suppose we now have a condition with cylindrical geometry like a fiber: we choose coordinate ρ, ϕ, z. What is the Helmholtz equation?

(((∂^2) E)/(∂ (z^2))) +(1/ρ)((∂ E)/(∂ ρ)) + (((∂^2) E)/(∂ (ρ^2))) + (1/(ρ^2))(((∂^2) E)/(∂ (ϕ^2))) +((ω^2)/(c^2))ϵ(ρ,ϕ, z) E = 0 
There are 2 types of optical waveguides with this geometry: fiber and circular waveguide. In fiber, the propagation direction is z, in circular waveguide, the propagation direction is ϕ. They are superficially similar, but in fact quite different in behavior.

2. Fiber geometry

2.0 Intro to fiber

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.

2.1 Approach to solve fiber mode

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
For fiber, z is the direction of propagation:
Let         E=(e^(i(β z - ω t)))  F[ρ,ϕ] 
Then:
(1/ρ)((∂ F)/(∂ ρ)) + (((∂^2) F)/(∂ (ρ^2))) + (1/(ρ^2))(((∂^2) F)/(∂ (ϕ^2)))  +(((ω^2)/(c^2))ϵ(ρ,ϕ)  -  (β^2) )F= 0  
Let's assume further angular geometry: E = (e^(i(β z - ω t))) F[ρ, ϕ]->(e^(i(β z - ω t))) (e^(i l ϕ)) F[ρ]
(1/ρ)((∂ F)/(∂ ρ)) + (((∂^2) F)/(∂ (ρ^2)))  +(((Subscript[k, 0])^2) ϵ(ρ)  -((l^2)/(ρ^2)) - (β^2) )F= 0 
What does the vector field look like? E, H  = (({{(Subscript[(E or H), ρ])}, {(Subscript[(E or H), ϕ])}, {(Subscript[(E or H), z])}}));

2.2 Boundary condition for fiber geometry

-(1/c)μ((∂ (Overscript[H, ->]))/(∂ t))= ∇×(Overscript[E, ->])

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

-(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, ρ])[ρ]}}))
and: ((-i ω)/c)ϵ  (Overscript[E, ->])  = ∇×(Overscript[H, ->])    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, ρ])[ρ]}}))

2.3 Reduction

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])[ρ]}})) 
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])[ρ]}}))    
and from:
((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])[ρ]}})) 
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])[ρ]}}))  
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])[ρ]}}))  
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])[ρ]}}))  ; 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])[ρ]}}))    
(({{(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])[ρ]}}))

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])[ρ]}}))  
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])[ρ]}})) 
 
(({{(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])[ρ]}})) 
 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])[ρ]}})) 
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])[ρ]}}))
So, we have solve for all components of E and H in terms of (Subscript[E, z]) and (Subscript[H, z]).

2.4 The use of S matrix

Can S-matrix formalism be used in this case for boundary condition? Why not? 
Let : (({{(Subscript[E, z ])}, {(∂_ρ) (Subscript[E, z])}}))= (({{(Subscript[ψ, A]), (Subscript[ψ, B])}, {(Subscript[ψ, A])', (Subscript[ψ, B])'}})).(({{(Subscript[c, A])}, {(Subscript[c, B])}})) ; (({{(Subscript[H, z ])}, {(∂_ρ) (Subscript[H, z])}}))= (({{(Subscript[ψ, A]), (Subscript[ψ, B])}, {(Subscript[ψ, A])', (Subscript[ψ, B])'}})). (({{(Subscript[d, A])}, {(Subscript[d, B])}})) . 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])}})) 
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])[ρ]}}))  
 (({{ϵ (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])}}))  
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])}}))  
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])}}))  
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, ϕ])}}))  ; (1/((ϵ μ ((Subscript[k, 0])^2)- (β^2) )))  = (1/Q) 
(({{(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])}}))

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])}}))   
(({{(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])}}))    
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])'}}))
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]))}}))]

3. Step-index fiber

3.1 Discussion

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
A common case is the step-index fiber, where  ϵ(ρ) is constant of ρ for each layer of the fiber. Although real fibers generally have more complex index of refraction profile, the case of step index fiber is very useful to illustrate the more general features common to all 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

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

3.1.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,100}
, 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]}} ]

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.

3.1.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 
(Subscript[S, l])= (((({{(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])'}}))
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 or 2])((Subscript[k, 0])^2) -(β^(2 ))]) ρ] in the core region, and (Subscript[K, l])[(Sqrt[(β^(2 ))-(Subscript[ϵ , 1 or 2])((Subscript[k, 0])^(2 ))]) ρ] in the clad. We further define: 
(Subscript[q, 1])=(Sqrt[(Subscript[ϵ , 1 or 2])((Subscript[k, 0])^2) -(β^(2 ))]) and (Subscript[q, 2])=(Sqrt[(β^(2 ))-(Subscript[ϵ , 1 or 2])((Subscript[k, 0])^2)]). 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);
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}}))

The determinant has to vanish for nontrivial solution to exist.

3.2 Characteristic equation (algebra)

3.2.1 Equation set up

FullSimplify[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]))}}))  ]]

3.2.2 Algebraic reduction

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

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

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

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

(q2 (Subscript[K, l]) (Subscript[J, l]^′)+q1 (Subscript[J, l]) (Subscript[K, l]^′)) (q2 (Subscript[K, l]) (Subscript[ϵ, 1]) (Subscript[J, l]^′)+q1 (Subscript[J, l]) (Subscript[ϵ, 2]) (Subscript[K, l]^′))==(1/((q1^2) (q2^2) (ρ^2) (Subscript[k, 0]^(2  ))))(l^2) (β^2) (Subscript[J, l]^2) (Subscript[k, 0]^4) (Subscript[K, l]^2) ((((Subscript[ϵ, 1])-(Subscript[ϵ, 2])))^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)

((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)

3.2.3 Numerical characteristic Eq. calculation

CharG[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]) ;
     (q1 J Kp+ q2 K Jp) (q1  ϵ2 J Kp+q2 ϵ1 K Jp)
	-(( (n^2))/((k0^2) (ρ^2) (q1^2) (q2^2)))(l^2) (((ϵ1-ϵ2))^2) K^2 J^2   ]    ;

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]}} ]

Indexprof={1.48,1.45};l=1;λ=1.55; rad=5.5/2;
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
Plot[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-> All]

l=0;
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-> All]

l=0;
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
Plot[Log[Abs[CharG[Indexprof,λ, rad,x,l]]] , {x,1.457,1.46}
		,PlotStyle->{{RGBColor[1,0.,0.]},{RGBColor[1,0,0]},{RGBColor[0,0,1]}} ]

l=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]}} ]

l=3;
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]}} ]

Indexprof={1.48,1.465};l=1;
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]}} ]

l=0;
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]}} ]

l=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]}} ]

Indexprof={1.48,1.465};l=1;λ=1.55; rad=5.5/2;
segment={1.472,1.474};
Plot[CharG[Indexprof,λ, rad,x,l] , {x,segment[[1]],segment[[2]]}
		,PlotStyle->{{RGBColor[1,0.,0.]},{RGBColor[1,0,0]},{RGBColor[0,0,1]}} ]

segment={1.472,1.474};l=1;
FindRoot[CharG[Indexprof,λ, rad,x,l] , {x,segment[[1]],segment[[2]]}
	          ]

3.2.4  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 ϕ))

3.3 Determine the mode profile

3.3.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.

3.3.2 Discussion: TE and TM mode (LP approx)

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.

3.3.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}}))

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.

3.4 LP approximation

Link to paraxial approx

4. Illustration of simple step-index fiber

4.1 Setup and calculation packages

4.1.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}]

(({{(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]^′))))}}))

4.1.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,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}] ;
		{field * (E^(I ϕ l)),field *(E^(-I ϕ l))}
		          ]

4.2 Fiber profile and mode structures

4.2.1 Effect of core diameter on the number of modes: single mode vs. multi-mode

Consider a very thin core fiber

Indexprof={1.49,1.48};l=1;λ=1.55; rad=5.5/2;
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
allp={} ;
For[l=0,l≤ 4, l++,
AppendTo[allp, 
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]}} ,
TicksStyle->Directive[Black,10]]
]
] ;
Show[GraphicsRow[allp]]

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.
Consider the same profile, but larger diameter:

Indexprof={1.49,1.48};l=1;λ=1.55; rad=50./2;
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
allp={} ;
For[l=0,l≤ 4, l++,
AppendTo[allp, 
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]}} ,
TicksStyle->Directive[Black,10]]
]
] ;
Show[GraphicsRow[allp]]

We see that for every l, there are NOT just one, but many radial modes. This is a multi-mode fiber. Let's see how this happens as we increase the diameter.

Indexprof={1.49,1.48};λ=1.55; 
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
l=1  ;
Manipulate[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]}} ,
TicksStyle->Directive[Black,10], PlotRange-> {-30,0}]
,{rad,5.5/2.,30./2}, AnimationRunning-> False]

Indexprof={1.49,1.48};λ=1.55; 
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
l=0  ;
Animate[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]}} ,
TicksStyle->Directive[Black,10], PlotRange-> {-30,0}]
,{rad,5.5/2.,30./2}, AnimationRunning-> False]

Indexprof={1.49,1.48};λ=1.55; 
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
l=2  ;
Manipulate[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]}} ,
TicksStyle->Directive[Black,10], PlotRange-> {-30,0}]
,{rad,5.5/2.,30./2}]

Indexprof={1.49,1.48};λ=1.55; 
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
l=3  ;
Manipulate[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]}} ,
TicksStyle->Directive[Black,10], PlotRange-> {-30,0}]
,{rad,5.5/2.,30./2}, AnimationRunning-> False]

Indexprof={1.49,1.48};λ=1.55; 
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
l=4  ;
Manipulate[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]}} ,
TicksStyle->Directive[Black,10], PlotRange-> {-30,0}]
,{rad,5.5/2.,30./2}, AnimationRunning-> False]

4.2.2 Effect of core-cladding index on the number of modes

Consider a very thin core fiber

Indexprof={1.49,1.48};l=1;λ=1.55; rad=5.5/2;
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
allp={} ;
For[l=0,l≤ 4, l++,
AppendTo[allp, 
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]}} ,
TicksStyle->Directive[Black,10]]
]
] ;
Show[GraphicsRow[allp]]

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.
Consider the same profile, but larger step index

Indexprof={1.53,1.48};l=1;λ=1.55; rad=5.5/2. ;
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
allp={} ;
For[l=0,l≤ 4, l++,
AppendTo[allp, 
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]}} ,
TicksStyle->Directive[Black,10]]
]
] ;
Show[GraphicsRow[allp]]

We see that for every l, there are NOT just one, but many radial modes. This is a multi-mode fiber. Let's see how this happens as we increase the step index.

λ=1.55; rad=5.5/2. ;
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
l=1  ;
Manipulate[Plot[Log[Abs[CharG[{n1,1.48},λ, rad,x,l]]] 
, {x,1.4805,n1-0.0005}
		,PlotStyle->{{RGBColor[1,0.,0.]},{RGBColor[1,0,0]},{RGBColor[0,0,1]}} ,
TicksStyle->Directive[Black,10], PlotRange-> {-30,0}]
,{n1,1.482,1.55}, AnimationRunning-> False]

λ=1.55; rad=5.5/2. ;
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
l=0  ;
Animate[Plot[Log[Abs[CharG[{n1,1.48},λ, rad,x,l]]] 
, {x,1.4805,n1-0.0005}
		,PlotStyle->{{RGBColor[1,0.,0.]},{RGBColor[1,0,0]},{RGBColor[0,0,1]}} ,
TicksStyle->Directive[Black,10], PlotRange-> {-30,0}]
,{n1,1.49,1.55}, AnimationRunning-> False]

λ=1.55; rad=5.5/2. ;
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
l=2  ;
Animate[Plot[Log[Abs[CharG[{n1,1.48},λ, rad,x,l]]] 
, {x,1.4805,n1-0.0005}
		,PlotStyle->{{RGBColor[1,0.,0.]},{RGBColor[1,0,0]},{RGBColor[0,0,1]}} ,
TicksStyle->Directive[Black,10], PlotRange-> {-30,0}]
,{n1,1.49,1.55}, AnimationRunning-> False]

4.2.3 Concept of cut-off wavelength

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!

Indexprof={1.49,1.48};l=1;λ=1.55; rad=5.5/2;
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
allp={} ;
For[l=0,l≤ 4, l++,
AppendTo[allp, 
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]}} ,
TicksStyle->Directive[Black,10]]
]
] ;
Show[GraphicsRow[allp]]

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?

Indexprof={1.49,1.48}; rad=5.5/2 ; 
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
l=1  ;
Manipulate[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]}} ,
TicksStyle->Directive[Black,10], PlotRange-> {-20,0}]
,{{λ,1.55},0.5,5}]

Indexprof={1.49,1.48}; rad=5.5/2 ; 
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
l=1  ;
Manipulate[Plot[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]}} ,
TicksStyle->Directive[Black,10], PlotRange->All]
,{λ,0.5,5}, AnimationRunning-> False]

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

5. Illustration of mode profiles

The following will illustrate a multimode step-index fiber. The wavelength for illustration is 1 μm.

5.1 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}] ;
		{field * (E^(I ϕ l)),field *(E^(-I ϕ l))}
		          ]

5.2 Plot mode search

Below shows the results.

Indexprof={1.49,1.45}; rad=5.5/2 ;  λ=1. ;
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
Manipulate[ 
NL=Min[nl,nh] ;NH=Max[nl,nh] ;
If[NL== NH,NH=NL+0.0005, ];
dn=Abs[NL-NH]/100;
Plot[Log[Abs[CharG[Indexprof,λ, rad,x,l]]] 
, {x,NL+dn,NH-dn}
,PlotStyle->{{RGBColor[1,0.,0.]},{RGBColor[1,0,0]},{RGBColor[0,0,1]}} ,
TicksStyle->Directive[Black,10], PlotRange-> {-NR,0}]
,{λ,1,10},{l,0,4,1},{{NR,40},0,60,5},{nl,Indexprof[[2]],Indexprof[[1]]},{{nh,Indexprof[[1]]},Indexprof[[2]],Indexprof[[1]]}]

5.3 Plot for l=1: fundamental l mode.

5.3.1 l=1, m=0 Intensity

segment={1.4852,1.4853};l=1;
neff= x /. FindRoot[CharG[Indexprof,λ, rad,x,l] , {x,segment[[1]],segment[[2]]}
	          ] 
coefp=Coeff[Indexprof,λ,rad,neff,l] ;
coefm=Coeff[Indexprof,λ,rad,neff,-l] ;

Plot intensity

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

Total intensity:

ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,Abs[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]]
+FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,-l,coefm][[1,3]] ]^2
+Abs[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] 
+FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,-l,coefm][[1,4]]  ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All
, PlotLabel->StringJoin["This is (|Eρ|)^2 +(|Eϕ|)^2 for l=",ToString[l]]
]

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:

Note that we can also plot total intensity of just l=1 or l=-1, without having to create 2 degenerate mode interference.

ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,Abs[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ]^2
+Abs[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]]    ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All
, PlotLabel->StringJoin["This is (|Eρ|)^2 +(|Eϕ|)^2 for l=",ToString[l]]
]

5.3.1 E field amplitude for m=0.

We can look at actual E field:

{ρ,ϕ}= {0.1,0.05} ;
Manipulate[GraphicsRow[ {ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ],3 Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ]
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, Mesh-> False, PlotPoints-> {8,15}
, PlotLabel->StringJoin["This is Eρ for l=",ToString[l] ] ]
, ParametricPlot3D[{ρ Cos[ϕ], ρ Sin[ϕ]
		,3 Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ]
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All,Mesh-> False, PlotPoints-> {8,15}
, PlotLabel->StringJoin["This is Eϕ for l=",ToString[l]] 
]}
] ,{t,0,1}]

To understand it, we need vector plot. It is sufficient to plot in 2D.

t=0. ; s=0.5 ; l=1 ;
Manipulate[Graphics[Table[{Hue[0.7,1,1.],Arrow[{{
ρ Cos[ϕ]-s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Cos[ϕ]
                +s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Sin[ϕ]
, ρ Sin[ϕ]-s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Sin[ϕ]
                  -s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Cos[ϕ] }
,{ρ Cos[ϕ]+s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Cos[ϕ]
                -s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Sin[ϕ]
, ρ Sin[ϕ]+s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Sin[ϕ]
                  +s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Cos[ϕ] }
}]}  ,{ρ,0.3,1.2 rad,1.2 rad/8},{ϕ,0,2 π, π/8}] ],{t,0,1},{{s,0.3},0.1,1.5}]

5.3.1 E field amplitude for m=1.

segment={1.4686,1.4687};l=1;
neff= x /. FindRoot[CharG[Indexprof,λ, rad,x,l] , {x,segment[[1]],segment[[2]]}
	          ] 
coefp=Coeff[Indexprof,λ,rad,neff,l] ;
coefm=Coeff[Indexprof,λ,rad,neff,-l] ;

We can look at actual E field:

Manipulate[GraphicsRow[ {ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ],3 Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ]
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, Mesh-> False, PlotPoints-> {8,15}
, PlotLabel->StringJoin["This is Eρ for l=",ToString[l] ] ]
, ParametricPlot3D[{ρ Cos[ϕ], ρ Sin[ϕ]
		,3 Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ]
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All,Mesh-> False, PlotPoints-> {8,15}
, PlotLabel->StringJoin["This is Eϕ for l=",ToString[l]] 
]}
] ,{t,0,1}]

To understand it, we need vector plot. It is sufficient to plot in 2D.

t=0. ; s=0.5 ; l=1 ;
Manipulate[Graphics[Table[{Hue[0.7,1,1.],Arrow[{{
ρ Cos[ϕ]-s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Cos[ϕ]
                +s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Sin[ϕ]
, ρ Sin[ϕ]-s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Sin[ϕ]
                  -s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Cos[ϕ] }
,{ρ Cos[ϕ]+s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Cos[ϕ]
                -s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Sin[ϕ]
, ρ Sin[ϕ]+s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Sin[ϕ]
                  +s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Cos[ϕ] }
}]}  ,{ρ,0.3,1.2 rad,1.2 rad/8},{ϕ,0,2 π, π/8}] ],{t,0,1},{{s,0.3},0.1,1.5}]

5.3.1 E field amplitude for m=2.

segment={1.4654,1.466};l=1;
neff= x /. FindRoot[CharG[Indexprof,λ, rad,x,l] , {x,segment[[1]],segment[[2]]}
	          ] 
coefp=Coeff[Indexprof,λ,rad,neff,l] ;
coefm=Coeff[Indexprof,λ,rad,neff,-l] ;

We can look at actual E field:

Manipulate[GraphicsRow[ {ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ],3 Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ]
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, Mesh-> False, PlotPoints-> {8,15}
, PlotLabel->StringJoin["This is Eρ for l=",ToString[l] ] ]
, ParametricPlot3D[{ρ Cos[ϕ], ρ Sin[ϕ]
		,3 Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ]
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All,Mesh-> False, PlotPoints-> {8,15}
, PlotLabel->StringJoin["This is Eϕ for l=",ToString[l]] 
]}
] ,{t,0,1}]

To understand it, we need vector plot. It is sufficient to plot in 2D.

t=0. ; s=0.5 ; l=1 ;
Manipulate[Graphics[Table[{Hue[0.7,1,1.],Arrow[{{
ρ Cos[ϕ]-s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Cos[ϕ]
                +s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Sin[ϕ]
, ρ Sin[ϕ]-s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Sin[ϕ]
                  -s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Cos[ϕ] }
,{ρ Cos[ϕ]+s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Cos[ϕ]
                -s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Sin[ϕ]
, ρ Sin[ϕ]+s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Sin[ϕ]
                  +s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Cos[ϕ] }
}]}  ,{ρ,0.3,1.2 rad,1.2 rad/8},{ϕ,0,2 π, π/8}] ],{t,0,1},{{s,0.3},0.1,1.5}]

5.4 Plot for l=0

5.4.1 l=0, m=0

segment={1.478,1.4785};l=0; λ=1. ;
neff= x /. FindRoot[CharG[Indexprof,λ, rad,x,l] , {x,segment[[1]],segment[[2]]}
	          ] 
coefp=Coeff[Indexprof,λ,rad,neff,l] ;
coefm=Coeff[Indexprof,λ,rad,neff,-l] ;

Plot intensity

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

Total intensity:

ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,Abs[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]]
+FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,-l,coefm][[1,3]] ]^2
+Abs[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] 
+FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,-l,coefm][[1,4]]  ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All
, PlotLabel->StringJoin["This is (|Eρ|)^2 +(|Eϕ|)^2 for l=",ToString[l]]
]

5.3.1 E field amplitude for m=0.

We can look at actual E field:

Manipulate[ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ],5 Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ]
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, Mesh-> False, PlotPoints-> {8,15}
, PlotLabel->StringJoin["This is Eρ for l=",ToString[l] ] 
, PlotRange-> {{-5.5,5.5},{-5.5,5.5},{-5,5}}
] ,{t,0,1}]

To understand it, we need vector plot. It is sufficient to plot in 2D.

t=0. ; s=0.5 ; l=0 ;
Manipulate[Graphics[Table[{Hue[0.7,1,1.],Arrow[{{
ρ Cos[ϕ]-s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Cos[ϕ]
                +s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Sin[ϕ]
, ρ Sin[ϕ]-s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Sin[ϕ]
                  -s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Cos[ϕ] }
,{ρ Cos[ϕ]+s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Cos[ϕ]
                -s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Sin[ϕ]
, ρ Sin[ϕ]+s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Sin[ϕ]
                  +s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Cos[ϕ] }
}]}  ,{ρ,0.3,1.2 rad,1.2 rad/8},{ϕ,0,2 π, π/8}] ],{t,0,1},{{s,0.3},0.1,1.5}]

5.4.1 l=0, m=1

segment={1.4779,1.478};l=0; λ=1. ;
neff= x /. FindRoot[CharG[Indexprof,λ, rad,x,l] , {x,segment[[1]],segment[[2]]}
	          ] 
coefp=Coeff[Indexprof,λ,rad,neff,l] ;
coefm=Coeff[Indexprof,λ,rad,neff,-l] ;

Plot intensity

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

Total intensity:

ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,Abs[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]]
+FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,-l,coefm][[1,3]] ]^2
+Abs[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] 
+FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,-l,coefm][[1,4]]  ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All
, PlotLabel->StringJoin["This is (|Eρ|)^2 +(|Eϕ|)^2 for l=",ToString[l]]
]

5.3.1 E field amplitude for m=1.

We can look at actual E field:

Manipulate[ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ],5 Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ]
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, Mesh-> False, PlotPoints-> {8,15}
, PlotLabel->StringJoin["This is Eρ for l=",ToString[l] ] 
, PlotRange-> {{-5.5,5.5},{-5.5,5.5},{-5,5}}
] ,{t,0,1}]

To understand it, we need vector plot. It is sufficient to plot in 2D.

t=0. ; s=0.5 ; l=0 ;
Manipulate[Graphics[Table[{Hue[0.7,1,1.],Arrow[{{
ρ Cos[ϕ]-s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Cos[ϕ]
                +s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Sin[ϕ]
, ρ Sin[ϕ]-s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Sin[ϕ]
                  -s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Cos[ϕ] }
,{ρ Cos[ϕ]+s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Cos[ϕ]
                -s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Sin[ϕ]
, ρ Sin[ϕ]+s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Sin[ϕ]
                  +s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Cos[ϕ] }
}]}  ,{ρ,0.3,1.2 rad,1.2 rad/8},{ϕ,0,2 π, π/8}] ],{t,0,1},{{s,0.3},0.1,1.5}]

5.4.1 l=0, m=2

segment={1.4526,1.4527};l=0; λ=1. ;
neff= x /. FindRoot[CharG[Indexprof,λ, rad,x,l] , {x,segment[[1]],segment[[2]]}
	          ] 
coefp=Coeff[Indexprof,λ,rad,neff,l] ;
coefm=Coeff[Indexprof,λ,rad,neff,-l] ;

Plot intensity

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

Total intensity:

ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,Abs[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]]
+FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,-l,coefm][[1,3]] ]^2
+Abs[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] 
+FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,-l,coefm][[1,4]]  ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All
, PlotLabel->StringJoin["This is (|Eρ|)^2 +(|Eϕ|)^2 for l=",ToString[l]]
]

5.3.1 E field amplitude for m=2.

We can look at actual E field:

Manipulate[ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ],5 Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ]
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, Mesh-> False, PlotPoints-> {8,15}
, PlotLabel->StringJoin["This is Eρ for l=",ToString[l] ] 
, PlotRange-> {{-5.5,5.5},{-5.5,5.5},{-5,5}}
] ,{t,0,1}]

To understand it, we need vector plot. It is sufficient to plot in 2D.

t=0. ; s=0.5 ; l=0 ;
Manipulate[Graphics[Table[{Hue[0.7,1,1.],Arrow[{{
ρ Cos[ϕ]-s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Cos[ϕ]
                +s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Sin[ϕ]
, ρ Sin[ϕ]-s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Sin[ϕ]
                  -s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Cos[ϕ] }
,{ρ Cos[ϕ]+s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Cos[ϕ]
                -s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Sin[ϕ]
, ρ Sin[ϕ]+s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Sin[ϕ]
                  +s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Cos[ϕ] }
}]}  ,{ρ,0.3,1.2 rad,1.2 rad/8},{ϕ,0,2 π, π/8}] ],{t,0,1},{{s,0.3},0.1,1.5}]

5.3 Plot for l=2: all modes

Indexprof={1.49,1.45}; rad=5.5/2 ;  λ=1. ;
dn=Abs[Indexprof[[1]]- Indexprof[[2]]]/100;
Manipulate[ 
NL=Min[nl,nh] ;NH=Max[nl,nh] ;
If[NL== NH,NH=NL+0.0005, ];
dn=Abs[NL-NH]/100;
Plot[Log[Abs[CharG[Indexprof,λ, rad,x,l]]] 
, {x,NL+dn,NH-dn}
,PlotStyle->{{RGBColor[1,0.,0.]},{RGBColor[1,0,0]},{RGBColor[0,0,1]}} ,
TicksStyle->Directive[Black,10], PlotRange-> {-NR,0}]
,{λ,1,3},{l,0,4,1},{{NR,40},0,60,5},{nl,Indexprof[[2]],Indexprof[[1]]},{{nh,Indexprof[[1]]},Indexprof[[2]],Indexprof[[1]]}]

5.3.1 l=2, m=0 Intensity

segment={1.47785,1.478};l=2;
neff= x /. FindRoot[CharG[Indexprof,λ, rad,x,l] , {x,segment[[1]],segment[[2]]}
	          ] 
coefp=Coeff[Indexprof,λ,rad,neff,l] ;
coefm=Coeff[Indexprof,λ,rad,neff,-l] ;

Plot intensity

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

Total intensity:

ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		, Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]]
-FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,-l,coefm][[1,3]]
  ]^2
+ Re[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]]
-FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,-l,coefm][[1,4]]  
  ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, Mesh-> False
, PlotLabel->StringJoin["This is (|Eρ|)^2 +(|Eϕ|)^2 for l=",ToString[l]]
]

5.3.1 E field amplitude for m=0.

We can look at actual E field:

Manipulate[GraphicsRow[ {ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ],3 Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ]
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, Mesh-> False, PlotPoints-> {8,15}
, PlotLabel->StringJoin["This is Eρ for l=",ToString[l] ] ]
, ParametricPlot3D[{ρ Cos[ϕ], ρ Sin[ϕ]
		,3 Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ]
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All,Mesh-> False, PlotPoints-> {8,15}
, PlotLabel->StringJoin["This is Eϕ for l=",ToString[l]] 
]}
] ,{t,0,1}]

To understand it, we need vector plot. It is sufficient to plot in 2D.

t=0. ; s=0.5 ; 
Manipulate[Graphics[Table[{Hue[0.7,1,1.],Arrow[{{
ρ Cos[ϕ]-s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Cos[ϕ]
                +s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Sin[ϕ]
, ρ Sin[ϕ]-s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Sin[ϕ]
                  -s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Cos[ϕ] }
,{ρ Cos[ϕ]+s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Cos[ϕ]
                -s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Sin[ϕ]
, ρ Sin[ϕ]+s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Sin[ϕ]
                  +s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Cos[ϕ] }
}]}  ,{ρ,0.3,1.2 rad,1.2 rad/6},{ϕ,0,2 π, π/12}] ],{t,0,1},{{s,0.3},0.1,1.5}]

5.3.1 E field amplitude for m=1.

segment={1.4574,1.4578};l=2;
neff= x /. FindRoot[CharG[Indexprof,λ, rad,x,l] , {x,segment[[1]],segment[[2]]}
	          ] 
coefp=Coeff[Indexprof,λ,rad,neff,l] ;
coefm=Coeff[Indexprof,λ,rad,neff,-l] ;

We can look at actual E field:

Manipulate[GraphicsRow[ {ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ],3 Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ]
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, Mesh-> False, PlotPoints-> {8,15}
, PlotLabel->StringJoin["This is Eρ for l=",ToString[l] ] ]
, ParametricPlot3D[{ρ Cos[ϕ], ρ Sin[ϕ]
		,3 Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ]
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All,Mesh-> False, PlotPoints-> {8,15}
, PlotLabel->StringJoin["This is Eϕ for l=",ToString[l]] 
]}
] ,{t,0,1}]

To understand it, we need vector plot. It is sufficient to plot in 2D.

t=0. ; s=0.5 ; 
Manipulate[Graphics[Table[{Hue[0.7,1,1.],Arrow[{{
ρ Cos[ϕ]-s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Cos[ϕ]
                +s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Sin[ϕ]
, ρ Sin[ϕ]-s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Sin[ϕ]
                  -s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Cos[ϕ] }
,{ρ Cos[ϕ]+s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Cos[ϕ]
                -s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Sin[ϕ]
, ρ Sin[ϕ]+s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Sin[ϕ]
                  +s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Cos[ϕ] }
}]}  ,{ρ,0.3,1.2 rad,1.2 rad/6},{ϕ,0,2 π, π/15}] ],{t,0,1},{{s,0.3},0.1,1.5}]

5.3.1 E field amplitude for m=2.

segment={1.4525,1.453};l=2;
neff= x /. FindRoot[CharG[Indexprof,λ, rad,x,l] , {x,segment[[1]],segment[[2]]}
	          ] 
coefp=Coeff[Indexprof,λ,rad,neff,l] ;
coefm=Coeff[Indexprof,λ,rad,neff,-l] ;

Plot intensity

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

Total intensity:

ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ]
		,Abs[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]]
+FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,-l,coefm][[1,3]] ]^2
+Abs[ FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] 
+FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,-l,coefm][[1,4]]  ]^2
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, Mesh-> False
, PlotLabel->StringJoin["This is (|Eρ|)^2 +(|Eϕ|)^2 for l=",ToString[l]]
]

We can look at actual E field:

Manipulate[GraphicsRow[ {ParametricPlot3D[
  {ρ Cos[ϕ], ρ Sin[ϕ],3 Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ]
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All, Mesh-> False, PlotPoints-> {8,15}
, PlotLabel->StringJoin["This is Eρ for l=",ToString[l] ] ]
, ParametricPlot3D[{ρ Cos[ϕ], ρ Sin[ϕ]
		,3 Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ]
		},
    {ρ , 0.01,7}, {ϕ, 0, 2 Pi}, PlotRange->All,Mesh-> False, PlotPoints-> {8,15}
, PlotLabel->StringJoin["This is Eϕ for l=",ToString[l]] 
]}
] ,{t,0,1}]

To understand it, we need vector plot. It is sufficient to plot in 2D.

t=0. ; s=0.5 ; l=2 ;
Manipulate[Graphics[Table[{Hue[0.7,1,1.],Arrow[{{
ρ Cos[ϕ]-s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Cos[ϕ]
                +s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Sin[ϕ]
, ρ Sin[ϕ]-s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Sin[ϕ]
                  -s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Cos[ϕ] }
,{ρ Cos[ϕ]+s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Cos[ϕ]
                -s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Sin[ϕ]
, ρ Sin[ϕ]+s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,3]] ] Sin[ϕ]
                  +s Re[ (E^(-I 2 π t)) FieldFiber[{ρ,ϕ},Indexprof,λ,rad,neff,l,coefp][[1,4]] ] Cos[ϕ] }
}]}  ,{ρ,0.2 rad,1 rad,0.8 rad/6},{ϕ,0,2 π, π/15}] ],{t,0,1},{{s,0.3},0.1,1.5}]