## Introduction

Many kinds of optical fiber can be described, or in the case of graded index fibers,

approximated by, a series of concentric layers of loss-less dielectric. When the index

contrast in the structure is small, it is common to use the scalar wave equation to

obtain the linearly-polarized modes (LP modes). However, when the index contrast

becomes larger, the LP approximation becomes inaccurate, and a full vector analysis

is necessary [1]. Within any one of the concentric layers, the field components are all

linear combinations of Ordinary and Modified Bessel functions, and the solution is a

matter of matching the tangential field components of adjacent layers, leading to the

solution of a linear system. There are some choices to make about the best way to

arrange the linear system (see Ref. [2] for a complete treatment) but probably the best

choice for a general, quick algorithm is the method of Yeh and Lindgren [3]. The

method has intuitive appeal in its similarity to the transfer matrix method used for the

analysis of planar waveguides [4]. Both methods use a matrix to express the fields at

one side of a layer given their values at the other side, and the whole structure is then

characterized by cascading the layers by matrix multiplication. For vector fields in

fibers, there is no natural separation into TE and TM polarization, so one must include

the two tangential components for both electric and magnetic field simultaneously. In

total there are four transverse components at the layer boundaries, so the transfer

matrix in fiber is a 4×4 matrix, instead of 2×2.

The method used by the Optiwave fiber mode solvers is different from the previous

works in three ways. The 4×4 matrix method for vector modes in fibers is

implemented with these improvements to make the algorithm faster in execution time,

and more reliable in the finding of modes, particularly in cases where the modes are

almost degenerate. The first improvement involves a reformulation of the basic

equations to make a real-valued numerical implementation. Since it is loss-less

modes on a fiber of loss-less material that is being sought, there is no reason to

include complex numbers in the formulation. Real-valued implementation uses less

computer resources. Second, the complete transfer across one layer boundary

requires a matrix inverse. It is better to find this inverse analytically, rather than rely

on numerical inversion at each stage. The layer matrix can be decomposed into two

4×4 matrices, and the new matrices are of a form where it is easy to find the inverse

by inspection, thereby finding the analytic solution to the original inversion problem.

Third, it is usual to construct the dispersion function (the zeros of which are located at

the modal indices) as the determinant of the 4×4 matrix system. While this is

theoretically correct, it is not the most intuitive prescription. Worse, the determinant is not the most convenient construction for locating the zeros, particularly in the case

where two zeros are very close to each other. Optiwave uses an eigen-value analysis

that splits the dispersion function into two functions. This splitting resolves the almost

degenerate pairs and helps the simple-minded computer algorithms to find the zeros

of the dispersion function (i.e. the modes) more reliably.

## Real-valued formulation

The time independent Maxwell curl equations with a positive time convention ( *e ^{jωt }*) are

and the divergence equations are

The factor j does not suit the current purpose, so make the substitution

where is the free space impedance, so that h and E have the same

units. Then

where . Eliminating the h gives the electric wave equation

## Debye Potential

For finding electromagnetic fields in a cylindrical geometry it is convenient to use

cylindrical coordinates r, θ and z, and Debye potentials parallel to the axis of rotation, [2]

It is supposed that the permittivity is constant in the region in which the above

equations apply. The fibers are defined as a series of concentric layers of constant

dielectric, so for the fiber mode solver, the region is the annulus between the layer

boundaries. Equation 6 and Equation 7 are applied in a piecemeal way. Strictly

speaking, the ψ, φ and ε should have subscripts to indicate to which layer the

solution applies, but they are dropped here to simplify the representation. The

complete solution for the multilayer fiber will be constructed by using a separate pair

of functions for each layer, and by matching the tangential field components at the

layer boundaries

The divergence of the curl of any vector field is identically zero. Therefore the choice

of (6) and (7) to represent the electric and magnetic fields means that the divergence

equations (2) will automatically be satisfied. In the remainder of this section, we show

that if the two potentials ψ and φ are solutions to the scalar Helmholtz equation, then

(6) will be a solution to the Maxwell wave equation (5), at least within any given layer.

In subsequent sections, the particular solution for the mode will be found by observing

the boundary conditions imposed by physical considerations on **E** and **h**.

Suppose ψ is a continuous function of position that satisfies the Helmholtz equation

in 3 dimensions

Consider the following equation, true for any ψ satisfying (9) in a region of constant

permittivity.

Associate the first two terms together and apply the vector identity

The right hand side of (11) is a gradient of a scalar function, so the curl of the left hand

side must be zero. take the curl of (11)

and define **E** by the first term of (6). With **E** defined this way, the Maxwell electric

wave equation (5) follows, and therefore this **E** is a possible solution for electric field.

Solutions defined from the Ψ function solely (that is to say, solutions with φ = 0), will

have electric fields of the form

from which we can see there is no longitudinal (z) component to the electric field, i.e.

these are transverse electric fields. The magnetic field associated with these

transverse electric fields is constructed from the first Maxwell curl equation in (4):

which is the second term of (7).

The same sequence (7) – (11) can be applied to the φ function to give transverse

magnetic fields, and the electric components of these are given as the second term

of (6). A linear superposition of transverse magnetic and transverse electric fields will

give solutions that are neither transverse magnetic nor transverse electric. In fact,

most modes of the fiber are of this hybrid kind. In the hybrid case, equations (6) and

(7) are used as written, and the relative value of the ψ and φ is now important, since

it is a specific linear combination that matches the boundary conditions at the layer

boundaries.

## Separation of Variables

The potentials themselves are solutions of the scalar Helmholtz equation, and the

particular solution is found by observing the boundary conditions imposed by physical

considerations on **E** and **h**. The potentials are supposed to form modes, so a solution

where the variables are separated is appropriate

and a similar expression applies for the other Debye potential, φ . The Helmholtz

equation (8) is expanded in cylindrical co-ordinates, and then (15) is substituted.For

regions where ε is constant, the radial functions follow the Ordinary or Modified

Bessel equation, and so the solutions are linear combinations of Bessel functions of

integer order v. In any given layer,

where . In layers where the propagation constant squared is

larger than k^{2}ε, the Bessel functions J and Y are replaced by the Modified Bessel

functions I and K, respectively.,

where . Equations (15), (16), and (17) are substituted in (6) and (7)

to find the electromagnetic fields. It is the tangential components θ and z that are

needed explicitly, since it is the tangential components of the electric and magnetic

fields that should match at layer boundaries. These field components are related to

the coefficients A, B, C, and D by a 4×4 matrix

where *n = β/k *is the modal index, and the common factor exp[*j(vθ – βz)*] is suppressed. In layers where the propagation constant squared is larger than *k ^{2}ε*, the Bessel functions J and Y are replaced by the Modified Bessel functions I and K, respectively, and the u is replaced by w.

A similar equation to (18) applies in adjacent layers, with different constants A, B, C,

D. Given the constants in one layer, the constants in the adjacent layer can be found

by solving the linear system created by the field matching condition. The f found by

the two matrices should be the same field vector at the boundary between layers. The

difference between this formulation and that in Ref. [2] and [3] is that this formulation

uses real numbers only.

## Solution of the linear system

The order of the components in the vector f is chosen to help factor the matrix in (9)

into two matrices as **Ff** = **Dc** where **c **= [A, B, C, D]* ^{T}* ,

The advantage of this factoring is that it is easy to find the expressions for **F**^{-1} and **D**^{-1}. **F**^{-1} is found by two row reductions and **D**^{-1} by inverting the 2×2 matrices in the

block diagonal. These matrices are required to find the coefficients in the next layer.

If **c**_{m }are the coefficients in the *m*-th layer, then the coefficients in the next layer is

found from the matrix product,

The use of the factored matrix means four matrix multiplications instead of two are

required to cross the boundary, but it avoids having to find the matrix inverse

numerically.

## Dispersion equation

The matrix to describe the whole multilayer can be written as one 4×4 matrix, the

product of all layers, one layer of which is shown in (21). The dispersion equation is

usually constructed by setting the determinant of this 4×4 matrix product to 0, but this

is not a good choice from the numerical point of view. In the limit of low index contrast,

the modal indices can be very close to each other, making them hard to find by

numerical methods, which tend to skip over closely spaced zeros. A more convenient

condition comes from using rectangular matrices at the first and last layers, which are

special cases. In the first layer, which contains r = 0, coefficients B_{0} and D_{0} must

be zero, so the first matrix, **D**_{0}, is 4 x 2 consisting of the first and third columns of the

**D** matrix in (20). Similarly, in the last layer, it is the coefficients A_{M} and C_{M} that must

go to zero in (17), since those contain solutions I_{V} with unphysical growth. The last

matrix B^{-1}_{m} will be of size 2 x 4, and consist of the first and third rows of **D**^{-1}. The

product of all these matrices will be of size 2×2, which we call S(n). (The dependence

on n is written explicitly since it is the modal index to be found). A mode is obtained

at a modal index n when a non-trivial solution is found to,

The condition det(**S**) = 0 is sufficient for a non-trivial solution. However this condition

can be inconvenient for the case of certain modes of the low index contrast fiber. For

example, consider a step index fiber with radius 2 μm, having a core index of 1.46

and a cladding of 1.45. For order zero ( v = 0 ), the determinant of **S** will vary with

the modal index n as shown in Figure 1.

**Figure 1: The determinant of the system matrix vs. the modal index. Step index fiber with radius 2 μm, core index of 1.46, and a cladding of 1.45. The optical wavelength is 0.5** **μm**,** and the order v = 0.**

It is difficult to find the zero for curves of this type. The position of the zero is easily

passed over in the initial search for zero crossings. The shape of the curve in the

region of the zero is very nearly quadratic, so the best solution is to split the

determinant function into a product of two functions. In that case, the zeros can be

found by looking for zero crossings of the two functions separately. Suitable functions

can be found by considering the vector [A_{0}, C_{0}], as an eigenvector of S with

eigenvalue, equal to 0. The eigenvalues of S(n) are ,

Where T = S_{11} + S_{22}. The determinant of a matrix is equal to the product of its

eigenvalues, so the product det(S(n)) = λ_{1}(*n*)λ_{2}(*n*) is a factoring convenient for

the zero finding task. A mode exists when either λ_{1}(*n*)* = 0 *or λ_{2}(*n*) = 0. The

advantage of searching for zeros of the two functions (23) is that the functions λ_{1}(*n*)

and λ_{2}(*n*) have zeros that are widely spaced apart, and are more easily found than

the zeros of the single determinant function.

**Figure 2: Plot of the eigenvalues and determinant vs modal index.**

**Figure 3: Expanded view of the zero crossing region**

## LP Modes

Linearly Polarized modes can be found with the same style of approach as above, by

a transfer matrix analysis of multilayer fibers.In the LP analysis, the initial excitation is

assumed to have a linear polarization, and so a modal analysis supporting only one

field component is used. It can be shown this is a good approximation for fibers with

low refractive index contrast, as shown below.

Consider again the wave equation of (5), but this time use the vector identity (10)

directly on the electric wave equation instead of introducing the Debye potential. The

result is

From the Maxwell divergence equation (2)

In the case of low index contrast waveguides, such as are found in most optical fibers,

the last term in (26) is negligible. The last term is the only one that contributes to

coupling among the field components of the electric field, so if this term is neglected,

the model will be a polarization preserving model. If the excitation was a linear

polarization parallel to the X axis, the model could be further simplified by considering

only the X component of (26). This shows that for the LP model, the governing

equation for the field is the Helmholtz equation (8) applied to the field component E_{x}.

The solution for LP modes is constructed in a similar way as for the vector modes.

The E_{x} is constructed in each layer as a linear combnation of Bessel functions (16)

and (17), except this time only two coefficients, *A _{1 }*and

*B*

_{1}, are required, since only

one function needs to be constructed. The physical consideration at the layer

boundaries is the continuity of E

_{x}and its derivative. The two conditions are related

to the two coefficients A and B for each layer by 2×2 matrices. The modes are found

by setting up a calculation that assumes

*A*

_{1}= 1 and

*B*for the first (inner

_{1 }= 0most) layer, and then calculating the coefficients in each subsequent layer by matrix

manipulations. The last layer has the condition that the coefficient (

*A*) of the Bessel

_{M }function I

_{V}must be zero, since this Bessel function is not bounded at infinity. This

last coefficient is a function of the order

*v*and the modal index n, so the modes are

found by finding the zeros of

*A*

_{M}(n).

## References

[1] K. Okamoto, Fundamentals of Optical Waveguides, (Academic Press, San Diego, 2000).

[2] C. Tsao, Optical Fibre Waveguide Analysis, (Oxford University Press, 1992), Part III

[3] C. Yeh and G. Lindgren, “Computing the propagation characteristics of radially stratified fibers: an efficient method”, Applied Optics 16(2) p483-493 (1977)

[4] J. Chilwell and I. Hodgkinson, “Thin-films field-transfer matrix theory of planar multilayer waveguides and reflection from prism-loaded waveguide”, Journal of the Optical Society of America A, 1 p742-753 (1984)