Optical mode solving for complex waveguides using a finite cloud method

D. R. Burke; T. J. Smy

doi:10.1364/OE.20.017783

1. Introduction

With the recent development of both optical fibers with intricate geometric structures [1] and the integration of a wide variety of geometrically complex optical waveguides into optoelectronic circuits there is a need to accurately model the propagation constants and modes of these devices. Determination of the guided modes for optical transmission provides a description of the field intensities, as well as the effective index of refraction (n_eff) and group velocity of the signals. The complex geometries of these structures precludes analytical approaches and mandates the use of numerical techniques. Typical vectorial field mode solvers use either; the Finite Difference Method (FDM) and a regular rectangular mesh [2], or the Finite Element Method (FEM) and a triangular mesh [3]. The FEM has previously shown itself to be the most practical for complicated curvilinear geometries present in structured fibers.

This paper proposes as an alternative to the above methods the use of the Finite Cloud Method (FCM) – a meshless method – using an unconstrained point distribution, originally described in [4]. The FCM is well suited to complicated curvilinear geometries and problems requiring flexibility in nodal distributions. The method will be briefly described below and has previously been adapted to solve inhomogeneous partial differential equations (PDE), such as the heat transfer equation [5] and stress-strain in solids [6].

A significant challenge in the numerical modeling of complex optical waveguides is the wide variety of possible structures formed out of multiple materials. Previous work solving inhomogeneous PDEs using FCM has been presented using the coupling of material regions through interface equations [5]. This paper presents the implementation of two methods for modeling the multiple material regions present in the optical waveguide structures. The hetero-junction is either modeled as an abrupt junction (step-index) or using a quasi-homogeneous method with a graded junction (graded-index). Both solution methods are compared with solutions from commercial solvers COMSOL MultiPhysics and Rsoft FemSIM [7, 8].

This initial work implements a full vectorial coupled H-field solution, and is applied to isotropic graded-index or step-index structures which are invariant in the z-direction. For this initial paper Dirichlet boundary conditions are used, thus restricting the solutions to guided modes with purely real propagation constants. However, extension to anisotropic media and leaky boundary conditions is expected to be possible. Although not appropriate for this type of problem, the method is readily adaptable to one, two and three dimensional solutions for a variety of other differential equations.

2. Finite cloud method

The Finite Cloud Method (FCM) is a meshless approach to solving partial differential equations. Instead of a geometrically structured region as in Finite Element or a regular grid as in Finite Difference, the FCM uses a simple distribution of points or nodes. These nodes are present throughout the entire computational domain and represent the unknown field values. Each node has a corresponding list of nearby nodes referred to as its finite cloud.

These clouds of nodes are then used to obtain an approximate form of the function to be solved for. The nodes can be given an irregular distribution and can vary from high density, for regions of high solution gradients or intricate small geometries; to low density regions of low solution gradients and simple geometries. The ease with which nodes can be distributed and the lack of overhead due to not creating a formal “mesh” is quite attractive and would be particularly useful in the automated meshing of complex geometric structures. The method is quite suitable, as well, to adaptable meshing as additional nodes can be easily added in any desired areas.

The FCM attempts to solve for the values of a discretized function u(x), with an approximate solution u^a(x). The form of the approximation for a 1D case is

u {(x)}^{a} = \sum_{I = 1}^{N P} ζ (x, x_{I}) φ (x_{K} - x_{I}) u_{I}

with a summation over I of the NP nodes present in the cloud, ζ(x, x_I) is a correction function and φ(x_K − x_I) an integration kernel centred on the cloud at point x_K. The integration kernel is a weighting function with delta distribution type properties, such as a gaussian.

The correction function must be calculated for each node after they have all been placed throughout the computational domain Ω. The unknown function is reproduced using a quadratic basis,

P^{T} (s) = [1, s, s^{2}], m = 3

of size m × 1. The correction function is created using

ζ (x, x_{I}) = P^{T} (x_{I}) C (x)

with P^T being the quadratic basis, and C(x) the correction function coefficients. Given that the unknown function is approximated with this basis function, the correction function must be able to reproduce the given quadratic function. Replacing the function u(x) in Eq. (1) with P(x) from Eq. (2) we have,

p_{i} (x) = \sum_{I = 1}^{N P} ζ (x, x_{I}) φ (x_{K} - x_{I}) p_{i} (x_{I}), i = 1, \dots, m .

This is the consistency condition which is enforced for every cloud to determine the correction coefficients for each corresponding node. The consistency condition process involves the creation and decomposition of a moment matrix, M, with

M C (x) = P (x) .

A unique moment matrix must be created for each node. The moment matrix must be created and decomposed for each particular cloud to determine the correction function coefficients relating the neighbouring nodes to the corresponding cloud’s node.

The interpolation function for the I’th node can be written as,

N_{I} (x) = P^{T} (x) M^{- 1} P (x_{I}) φ (x_{K} - x_{I})

where N_I(x) is the interpolation vector for a single node and also the I’th row of N(x), the interpolation matrix for all points in the domain. The derivatives of the interpolation matrix can be easily determined using the derivatives of P^T (x) in Eq. (6) as this is the only function in the equation dependent on the variable x. The interpolation matrix derivatives are represented as N_x and N_xx for the first and second derivatives in x respectively.

Each row I of the matrix N_xx represents the derivative operator d²/dx² acting on the function u(x) at the point x = x_I. The entire matrix is thus operating on a vector U approximating a continuous function u(x) at the nodal positions and can be used to solve discretized approximations of partial differential equations. A PDE for the unknown function u can therefore be solved by direct use of the derivatives of the shape array. For example the two dimensional PDE over a simple single material domain,

σ (\frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}}) = ρ (x, y),

becomes the equation,

σ (N_{xx} U + N_{yy} U) = R

where σ is a constant material parameter and R is the value of ρ at each node.

For a simple domain consisting of a single material with smoothly varying material parameters the use of the operator matrices is also straightforward with the rows of the N matrices weighted by the value of the material parameter of the appropriate node.

Of course boundary nodes must be defined and appropriate modifications made to N_xx to incorporate these into the equation. For example when a fixed or Dirichlet boundary is applied to node I all the elements on the I’th row of the matrix formed by N_xx + N_yy would be set to zero except for the diagonal element which should be set to one. Similarly Neumann boundary conditions would use the shape functions derived from N_x and N_y on the corresponding row.

When multiple materials are present in the domain with “hard” boundaries two approaches are possible [5]. The first approach is to use multiple coupled domains. Such a formulation defines interface nodes which are shared between the material domains. Coupling equations (which may involve the surface normal of the interface) are then defined for each interface node and introduced into the solution matrices. In effect, the problem is defined as a coupled set of single material domains, where the coupling is through the boundary conditions. The method has the advantage of creating an ideal “hard” interface and we refer to it as the “step-index” method for this paper.

The second approach is to treat a multi-material model as a single material. At the boundaries of the differing materials the parameters used for a cloud are calculated from an average of two materials present. As we shall see for a vectorial PDE the spatial derivatives of the material parameters may also be required. For this paper we refer to this method as the “graded-index” solution and we will require knowledge of the permittivity throughout the domain as well as its derivatives.

3. Optical mode solving

Optical propagation along a waveguide is described by Maxwell’s equations [9] and is a 2D eigenvalue problem in (x, y) with propagation in the z direction and having a free space propagation constant given by k₀ = 2π/λ₀ with free space wavelength λ₀. The determined eigen-vectors H_x and H_y correspond to the eigenvalue, β, which is the guided propagation constant which can be used to find the effective refractive index of the mode and material, n_eff = β/k₀.

Two non-coupled equations used for homogeneous isotropic regions are as follows [2]

\begin{array}{l} \frac{\partial^{2} H_{x}}{\partial x^{2}} + \frac{\partial^{2} H_{x}}{\partial y^{2}} + (ε_{r} k_{0}^{2} - β^{2}) H_{x} = 0 \\ \frac{\partial^{2} H_{y}}{\partial x^{2}} + \frac{\partial^{2} H_{y}}{\partial y^{2}} + (ε_{r} k_{0}^{2} - β^{2}) H_{y} = 0 \end{array}

with ε_r the relative permittivity of the homogeneous material.

When representing the waveguide as coupled regions of constant permittivity (step-index method) the boundary conditions for the H-fields on a junction between two regions of differing refractive index, regions a and b, must also be satisfied. These conditions also provide a coupling between the H_x and H_y fields. Ensuring the continuity of the H_z field across a junction gives [2]

{\frac{\partial H_{x}}{\partial_{x}} |}_{a} + {\frac{\partial H_{y}}{\partial_{y}} |}_{a} = {\frac{\partial H_{x}}{\partial_{x}} |}_{b} + {\frac{\partial H_{y}}{\partial_{y}} |}_{b}

and as well, ensuring the continuity of the E_z field yields [2]

{\frac{1}{ε_{a}} \frac{\partial H_{x}}{\partial_{y}} |}_{a} - {\frac{1}{ε_{a}} \frac{\partial H_{y}}{\partial_{x}} |}_{a} = {\frac{1}{ε_{b}} \frac{\partial H_{x}}{\partial_{y}} |}_{b} - {\frac{1}{ε_{b}} \frac{\partial H_{y}}{\partial_{x}} |}_{b}

3.1. Eigenmodes and eigenvalues

To solve for the eigenmodes and eigenvalues one must enforce the above conditions, ensuring that each matrix line includes the eigen-vector and -value on the right hand side of the equation. First the 2D homogeneous equation is rearranged to isolate for the eigenvalues and vectors,

\begin{array}{l} \frac{\partial^{2} H_{x}}{\partial x^{2}} + \frac{\partial^{2} H_{x}}{\partial y^{2}} + ε_{r} k_{0}^{2} H_{x} = β^{2} H_{x} \\ \frac{\partial^{2} H_{y}}{\partial x^{2}} + \frac{\partial^{2} H_{y}}{\partial y^{2}} + ε_{r} k_{0}^{2} H_{y} = β^{2} H_{y} \end{array}

An example of the implementation of the FCM in matrix form on this homogeneous equation is

[\begin{matrix} N_{x x} + N_{y y} + I ε_{r} k_{0}^{2} & 0 \\ 0 & N_{x x} + N_{y y} + I ε_{r} k_{0}^{2} \end{matrix}] [\begin{matrix} H_{x} \\ H_{y} \end{matrix}] = β^{2} [\begin{matrix} H_{x} \\ H_{y} \end{matrix}]

with a constant material permittivity, ε_r, and I being the identity matrix. When multiple materials are present the appropriate ε_r for a particular node must be used in forming these equations.

Boundary conditions, which couple the x and y fields, must also be satisfied for the nodes which lie along a material heterojunction. The matrix lines for these nodes must also include the eigenvalues and vectors which are to be solved in this system of equations. Two methods of modeling heterojunctions have been implemented. The first method, the step-index method, calculates a separate homogeneous cloud for both sides of the junction (clouds “a” and “b”). Each cloud approximates the region on either side of the heterojunction and each is given a weighting factor, which sum to unity, to ensure that the addition of both eigenvalue equations sum to a single eigenvalue. There are also two heterojunction conditions, (10) and (11), which must be enforced. The first condition, Eq. (10) is included in the H_x line of the matrix and can be summed along with the weighted eigenvalue equations as,

\begin{array}{l} M_{a x} (\nabla_{a}^{2} H_{x} + k_{0}^{2} ε_{a} H_{x}) + M_{b x} (\nabla_{b}^{2} H_{x} + k_{0}^{2} ε_{b} H_{x}) + \\ ({\frac{\partial H_{x}}{\partial x} |}_{a} + {\frac{\partial H_{y}}{\partial y} |}_{a}) - ({\frac{\partial H_{x}}{\partial_{x}} |}_{b} + {\frac{\partial H_{y}}{\partial_{y}} |}_{b}) = β^{2} H_{x} \end{array}

and the second, Eq. (11) is included in the H_y line of the matrix as

\begin{array}{l} M_{a y} (\nabla_{a}^{2} H_{y} + k_{0}^{2} ε_{a} H_{y}) + M_{b y} (\nabla_{b}^{2} H_{y} + k_{0}^{2} ε_{b} H_{y}) + \\ ({\frac{1}{ε_{a}} \frac{\partial H_{x}}{\partial_{y}} |}_{a} - {\frac{1}{ε_{a}} \frac{\partial H_{y}}{\partial_{x}} |}_{a}) - ({\frac{1}{ε_{b}} \frac{\partial H_{x}}{\partial_{y}} |}_{b} - {\frac{1}{ε_{b}} \frac{\partial H_{y}}{\partial_{x}} |}_{b}) = β^{2} H_{y} . \end{array}

The cloud weighting factor M_nm is constructed so that the the interface conditions (10) and (11) are satisfied and also that M_ax + M_bx = 1.0 and M_ay + M_by = 1.0,

\begin{array}{l} M_{n x} = [(\frac{ε_{n}}{ε_{a} + ε_{b}}) {cos}^{2} ϕ + (\frac{1}{2} {sin}^{2} ϕ)] \\ M_{n y} = [(\frac{ε_{n}}{ε_{a} + ε_{b}}) {sin}^{2} ϕ + (\frac{1}{2} {cos}^{2} ϕ)] \end{array}

with ϕ being the angle normal to the plane of the interface measured counter-clockwise from the positive x-axis. The weighting factor, seen in (16), comes about due to the inverse weighting of the fields perpendicular to the interface normal seen in (11), and ensures that the values on the right hand side of the equation sum to unity giving exclusively β²H on the right hand side.

For the second method, the graded-index method, the solution is treated as homogeneous with each finite cloud having its own value of permittivity and corresponding spatial derivatives. The permittivity and its derivative being taken at the centre of the cloud. The eigenvalue equation solved for this solution is [10],

\begin{array}{l} \frac{\partial^{2} H_{x}}{\partial x^{2}} + \frac{\partial^{2} H_{x}}{\partial y^{2}} - \frac{1}{ε} \frac{\partial ε}{\partial y} \frac{\partial H_{x}}{\partial y} + \frac{1}{ε} \frac{\partial ε}{\partial y} \frac{\partial H_{y}}{\partial x} + ε_{r} k_{0}^{2} H_{x} = β^{2} H_{x} \\ \frac{\partial^{2} H_{y}}{\partial x^{2}} + \frac{\partial^{2} H_{y}}{\partial y^{2}} - \frac{1}{ε} \frac{\partial ε}{\partial x} \frac{\partial H_{y}}{\partial x} + \frac{1}{ε} \frac{\partial ε}{\partial x} \frac{\partial H_{x}}{\partial y} + ε_{r} k_{0}^{2} H_{y} = β^{2} H_{y} \end{array}

which can be transformed into an eigenvalue equation similar to Eq. (13) by direct application of N and its derivatives.

3.2. Implementation

An optical mode solving routine based on the equations above has been created using the FCM, as described in [4, 5]. The implementation is in the “C” computer language and initially creates a node distribution from geometric information. This node distribution is then used in conjunction with interface normal information to create the shape matrix N and its derivatives such as N_x, N_xx as well as the material parameter vector ε, the permittivity at each node, and its spatial derivatives ε_x and ε_y. All of this information is then used to form an eigenvalue problem similar to (17). The eigenvalues and vectors are then found using either Matlab [11] or Octave [12].

The boundary conditions used in this paper are Dirichlet fixed values giving a perfect magnetic and electric conductor. This condition restricts solutions to solely real valued n_eff and fully guided modes. For the solution to include complex values from leaky modes a more complicated leaky boundary condition could be used.

3.3. Symmetry

In cases where there is a physical symmetry in the waveguide, the domain of the simulation can be reduced to save on computational resources with the use of symmetry boundary conditions. These conditions are the Perfect Electric Conductor (PEC),

H_{⊥} = 0

and the Perfect Magnetic Conductor (PMC),

H_{| |} = 0

Solution accuracy for a variety of fibers will be compared for a full solution, half symmetry and quarter symmetry. The time required to create the domain, matrices and solve the eigenvalue equation, as well as the domain size are also compared for the three symmetry situations. In this paper the use of symmetry was investigated by the simple creation of a reduced geometric structure and then the imposition of appropriate boundary conditions, however, other more mathematical approaches would also be possible [13, 14].

4. Results

The following section details initial results of using the FCM to solve for the guided modes of several microstructured optical waveguides. In all cases the figures show the node placements, physical parameters and the entire computational domain of the structure. Dirichlet boundary conditions are used, holding the fields to be zero at the domain edges.

4.1. Optical waveguides

For an initial investigation of the applicability of the FCM to optical mode solving two simple waveguides were modeled: 1) a ridge waveguide typical of a Silica on Silicon integrated optical platform and 2) a step-index optical fiber.

4.1.1. Ridge waveguide

As an initial test, a simple ridge waveguide, with parameters as shown in Fig. 1(a), was analyzed using the FCM and compared with the solutions from commercial modelling software COMSOL MultiPhysics [7] and Rsoft FemSIM [8]. A similar grid size and computation window were used for all of the solutions, width dx ≈ 0.05 μsm. For the FCM four material regions were defined. Three layers producing the basic layer structure consisting of a bottom Silicon layer (n=3.34), a thin 0.2 μm layer of Silica (n=3.44) and a top layer of Air (n=1). Within this structure was placed a ridge of dimensions 1.1 μm × 2 μm (n=3.44). The nodes were distributed uniformly except at the material interfaces where two extra rows of nodes were added on either side of the interface (see Fig. 1(a)). For the graded-index method the permittivity was smoothed over this region of nodes.

Fig. 1 Parameters and dimensions for (a) a ridge waveguide (b) a solid core step-index fiber.

Download Full Size | PDF

The effective index of refraction for the first six modes of the ridge are shown in Table 1, with the eigenvectors from the FCM shown in Fig. 2. Comparison of the different results shows a high degree of agreement for all of the methods. Numerical investigation determined that the results of the FCM for n_eff were quite robust with regards to node density and distribution and the specification of the gradient used for the graded-index method.

Fig. 2 The first six modes of the ridge waveguide.

Download Full Size | PDF

Table 1. Comparison of effective index of refraction for the first six modes of the ridge waveguide. The s-FCM and g-FCM being the step-index FCM and the graded-index FCM, compared with results from Rsoft FemSIM and COMSOL MultiPhysics.

View Table | View all tables in this article

4.1.2. Solid core waveguide

A second test, that of a fiber optic cable with a solid core was performed and also compared with the two commercial solvers. The parameters and computational window of the solid core waveguide are shown in Fig. 1(b). This problem involves a very low contrast waveguide and has circular symmetry, as opposed to the Manhattan layout of the ridge waveguide. Therefore with this structure a radial distribution of nodes was used with a constant density throughout the region. At the interface between the core and the cladding additional nodes were inserted as for the ridge example.

Results for the first six modes of the solid core fiber are shown in Table 2 which show a high agreement with the commercial solvers on the same structure, the corresponding eigenvectors are seen in Fig. 3.

Fig. 3 The first six modes of the solid core waveguide.

Download Full Size | PDF

Table 2. Comparison of effective index of refraction for the first six modes of the solid core waveguide. The s-FCM and g-FCM being the step-index FCM and the graded-index FCM, compared with results from Rsoft FemSIM and COMSOL MultiPhysics.

View Table | View all tables in this article

The eigenvalue solutions for the step-index fiber have also been compared for various levels of symmetry, the entire domain, half symmetry and quarter symmetry. Results for all three situations and both graded and step methods are compared in Table 3. Values for the number of solved nodes and time for complete solution, including domain creation, matrix creation, and sparse eigenvector solving are also shown in the table. As can be seen all configurations compare well with the commercial simulations and it was found that as with the ridge the results were robust with respect to changes in the node distribution. For example an orthogonal node distribution with a radially defined interface produced very similar numbers.

Table 3. Comparison of effective index of refraction for the first six modes of the step fiber waveguide. The s-FCM and g-FCM being the step-index FCM and the graded-index FCM, compared with half and quarter symmetry using PEC and PMC boundaries. ^*Simulation time using a dual core iMac at 2.4GHz.

View Table | View all tables in this article

4.2. Two structured waveguides

The previous two examples were relatively standard optical waveguides which are based on total internal reflection due to an optical difference between core and cladding. In this section we will present the analysis of two structured optical fibers. The first is an idealized air core fiber suitable for applications where non-linearities are to be kept to a minimum and is based on Bragg diffraction. The second is a structured fiber based on Photonic Crystal principles.

4.2.1. Bragg diffraction air core

The Bragg diffraction based waveguide, Fig. 4(a), used in this section is a simplification of the structure examined in [3] and is an idealized air core waveguide. The guide consists of a silica outer region at a radius of 19 μm and an inner air core of 10 μm. The air core is surrounded by three annular shaped air regions having thickness t_annular = 2.3 μm each separated by a thin silica ring of thickness t_ring = 0.2 μm. The index for the silica is calculated to be n₁ = 1.213567 and the index of air is n₀ = 1, with free space wavelength λ₀ = 1.06 μm.

Fig. 4 Parameters and dimensions for (a) a Bragg diffraction air core fiber and (b) a photonic crystal fiber with 6 circular air holes.

Download Full Size | PDF

Radial meshing, with single row at the interface, is used and the s-FCM results were obtained using a full model and approximate spacing of d_rad ≈ 0.1μm. The graded method results required a slightly finer mesh at the interfaces as it was found that a single node in the ring thickness of 0.2μm was not sufficient to capture the sharp index change from air-silica-air. Thus, a spacing of d_rad ≈ 0.066 μm was used for the g-FCM results and quarter symmetry was also needed to reduce the computational window size resulting from the added nodes. The first six modes, Fig. 5, are compared in Table 4 for the two methods and two commercial solvers all in very good agreement.

Fig. 5 The first six modes of the Bragg diffraction air core structure.

Download Full Size | PDF

Table 4. Comparison of effective index of refraction for the first six modes of the Bragg diffraction air core structure. The s-FCM and g-FCM being the step-index FCM and the graded-index FCM, compared with results from Rsoft FemSIM and COMSOL Multi-Physics.

View Table | View all tables in this article

4.2.2. Photonic crystal fibers

Finally, a complex microstructured air hole waveguide is also examined using the FCM and compared with the commercial solvers. The structure is a photonic crystal fiber with 6 circular air holes with parameters and dimensions seen in Fig. 4(b), as published in [3]. Radial meshing was used for all seven discs; one primary disc and the 6 inset air holes, with triple row interface used for boundaries.

The fundamental quasi-degenerate modes of the structure are compared for the two methods and results shown in Table 5, and the modes shown in Fig. 6. Again, the two FCM solutions give highly accurate results in agreement with the compared methods.

Fig. 6 The first six modes of the circular air hole photonic crystal fiber.

Download Full Size | PDF

Table 5. Comparison of effective index of refraction for the first six modes of the air hole waveguide. The s-FCM and g-FCM being the step-index FCM and the graded-index FCM, compared with results from Rsoft FemSIM and COMSOL MultiPhysics.

View Table | View all tables in this article

As with the step-fiber example, the air hole fiber has a physical symmetry which can be used to reduce the domain size and computational load. The three domain sizes, full, half, and quarter, are compared in Table 6 along with the full solution time and number of nodes in the domain.

Table 6. Comparison of effective index of refraction for the first six modes of the air hole waveguide. The s-FCM and g-FCM being the air-hole FCM and the graded-index FCM, compared with half and quarter symmetry using PEC and PMC boundaries.

View Table | View all tables in this article

A further test of the photonic crystal fiber is to investigate the dispersion properties of the fiber in the region of operation around its operating wavelength. The dispersion parameter is based on the second derivative of the effective index with respect to the wavelength and is calculated as [3]

Dispersion = - \frac{λ}{c} \frac{\partial^{2}}{\partial λ^{2}} [ℜ (n_{eff})]

For simplicity in this comparison the index of refraction for the materials is kept constant with a varying wavelength. A more rigorous fiber characterization would include the variation of the index of refraction as a function of the wavelength. The dispersion of the fundamental mode and the third mode from Fig. 6 for the photonic crystal fiber is compared using the FCM against the results from COMSOL MultiPhysics and seen in Fig. 7. The results show that the FCM is in very good agreement with the dispersion found using a commercial solver.

Fig. 7 Dispersion parameter comparison between COMSOL MultiPhysics and the FCM for modes 1 and 3 of the photonic crystal fiber.

Download Full Size | PDF

4.3. Solution interpolation

A beneficial feature of the Finite Cloud Method is that, as mentioned previously, during the creation of the shape matrix all of the derivatives at the nodes can be trivially calculated and stored. This allows the final solution points and the derivative matrix to be used along with a Taylor series expansion to interpolate for any given point in the domain. The solution is found using a quadratic basis, hence a second order Taylor series expansion as shown in (21) is used for the interpolation.

f (x) = f (a) + \frac{f^{'} (a)}{1!} (x - a) + \frac{f^{″} (a)}{2!} {(x - a)}^{2}

A point is first chosen for the location at which one wishes to determine the unknown field. This unknown field is found using the closest node to the point along with its field solution, and first and second order derivatives. This is demonstrated in Fig. 8 examining two second order modes of the ridge and air-hole structures with the fields sampled spatially at four times the node density. Results from the interpolation show a smooth and continuous function throughout the domain.

Fig. 8 Mode amplitudes for two structures. (a) Field magnitude of a 2nd order mode for ridge (b) H_x plotted for a 2nd order mode for the air-hole fiber.

Download Full Size | PDF

5. Conclusion

In this paper a meshless method of solving vectorial field modes for micro-structured optical waveguides has been presented, the Finite Cloud Method (Readers interested in the code are encouraged to contact the authors). The FCM has been used to solve the propagation constants and fields for four waveguide structures, and compared with numerical solutions from commercial solvers COMSOL MultiPhysics and Rsoft FemSIM. The dispersion parameter for the photonic crystal fiber was also compared with a commercial solver.

Both the step index and graded index methods for dealing with materially inhomogeneous problems produced very accurate results. The step method has the benefit of forcing hard boundaries between the regions. As a drawback, however, the added constraints on the interface nodes can add numerical problems or difficulties due to enforcing the interface conditions. The graded index enforces a slightly softer boundary to material interfaces, which can depend on the node density, but gives a more robust solution as additional matching conditions are not required. There exists as well the possibility of a hybrid method which could involve a variable index throughout the domain as well as hard material transitions. This would make use of both methods playing to their strengths and advantages.

Further work would include the implementation of a leaky boundary condition such as a TBC or PML allowing for the solution of leaky modes with complex effective indices of refraction. Fundamentally, the meshless method is very flexible and straightforward in its implementation; enabling adaptive node mapping techniques which may prove to be valuable in solving the modes of finely structured fibers and the incorporation of complex and anisotropic material permittivities present in fibers used for sensors.

References and links

1. P. S. Russell, “Photonic-crystal fibers,” J. Lightwave Technol. 24, 4729–4749 (2006). [CrossRef]

2. P. Lusse, P. Stuwe, J. Schule, and H.-G. Unger, “Analysis of vectorial mode fields in optical waveguides by a new finite difference method,” J. Lightwave Technol. 12, 487–494 (1994). [CrossRef]

3. H. Uranus and H. Hoekstra, “Modelling of microstructured waveguides using a finite-element-based vectorial mode solver with transparent boundary conditions,” Opt. Express 12, 2795–2809 (2004). [CrossRef] [PubMed]

4. N. Aluru and G. Li, “Finite cloud method: a true meshless technique based on a fixed reproducing kernel approximation,” Int. J. Numer. Methods Eng. 50, 2373–2410 (2001). [CrossRef]

5. D. Burke, S. Moslemi-Tabrizi, and T. Smy, “Simulation of inhomogeneous models using the finite cloud method,” Materialwiss. Werkstofftech. 41, 336–340 (2010). [CrossRef]

6. J.-S. Chen, S. Yoon, H.-P. Wang, and W. K. Liu, “An improved reproducing kernel particle method for nearly incompressible finite elasticity,” Comput. Methods Appl. Mech. Eng. 181, 117–145 (2000). [CrossRef]

7. COMSOL Multiphysics, Version 4.1 Comsol Inc. (2011), http://www.comsol.com.

8. Rsoft FemSim, Version 3.3 Rsoft Inc. (2011), http://www.rsoftdesign.com/products.php?sub=Component+Design$\&$itm=FemSIM.

9. J. D. Jackson, Classical Electrodynamics (Wiley, 1998).

10. K. Ramm, P. Lusse, and H.-G. Unger, “Multigrid eigenvalue solver for mode calculation of planar optical waveguides,” IEEE Photon. Technol. Lett. 9, 967–969 (1997). [CrossRef]

11. MathWorks (2011), http://www.mathworks.com/products/matlab/.

12. GNU Octave (2011), http://www.gnu.org/s/octave/.

13. P. McIsaac, “Symmetry-induced modal characteristics of uniform waveguides-i: Summary of results,” IEEE Trans. Microwave Theory Tech. 23(5), 421–429 (1975). [CrossRef]

14. J.M. Fini, “Improved symmetry analysis of many-moded microstructure optical fibers,” J. Opt. Soc. Am. B 21(8), 1431–1436 (2004). [CrossRef]

Mode	s-FCM	g-FCM	Rsoft	COMSOL

1	3.388803	3.387818	3.388718	3.388699
2	3.387870	3.387500	3.387854	3.387872
3	3.333842	3.335003	3.333806	3.336465
4	3.330638	3.330539	3.330954	3.334248
5	3.330209	3.329923	3.330079	3.332082
6	3.329599	3.327569	3.329560	3.331420

Mode	s-FCM	g-FCM	Rsoft	COMSOL

1	3.413099	3.413097	3.413093	3.413093
2	3.413097	3.413097	3.413093	3.413093
3	3.410550	3.410549	3.410527	3.410531
4	3.410525	3.410524	3.410509	3.410512
5	3.410525	3.410520	3.410509	3.410511
6	3.410504	3.410499	3.410494	3.410495

Mode	s-FCM	s-FCM $\frac{1}{2}$	s-FCM $\frac{1}{4}$	g-FCM	g-FCM $\frac{1}{2}$	g-FCM $\frac{1}{4}$

1	3.413099	3.413108	3.413107	3.413097	3.413098	3.413097
2	3.413097	3.413097	3.413096	3.413097	3.413097	3.413097
3	3.410550	3.410580	3.410576	3.410549	3.410550	3.410549
4	3.410525	3.410540	3.410539	3.410524	3.410525	3.410523
5	3.410525	3.410525	3.410525	3.410520	3.410525	3.410521
6	3.410504	3.410501	3.410500	3.410499	3.410503	3.410499

# of Nodes	18445	9495	4840	18445	9495	4840

Time^* (s)	38	15	7	29	13	7

Mode	s-FCM	g-FCM	Rsoft	COMSOL

1	0.999111	0.998704	0.999107	0.999109
2	0.999111	0.999370	0.999107	0.999109
3	0.997792	0.997788	0.997787	0.997788
4	0.997750	0.997495	0.997742	0.997744
5	0.997750	0.997498	0.997742	0.997744
6	0.997712	0.997268	0.997694	0.997704

Mode	s-FCM	g-FCM	Rsoft	COMSOL

1	1.445402	1.445411	1.445394	1.445395
2	1.445400	1.445409	1.445394	1.445395
3	1.438605	1.438635	1.438581	1.438582
4	1.438384	1.438409	1.438355	1.438358
5	1.438466	1.438494	1.438442	1.438440
6	1.438466	1.438490	1.438440	1.438440

Optical mode solving for complex waveguides using a finite cloud method

Abstract

1. Introduction

2. Finite cloud method

3. Optical mode solving

3.1. Eigenmodes and eigenvalues

3.2. Implementation

3.3. Symmetry

4. Results

4.1. Optical waveguides

4.1.1. Ridge waveguide

4.1.2. Solid core waveguide

4.2. Two structured waveguides

4.2.1. Bragg diffraction air core

4.2.2. Photonic crystal fibers

4.3. Solution interpolation

5. Conclusion

References and links

Cited By

Figures (8)

Tables (6)

Equations (21)

Optics Express

Mode	s-FCM	s-FCM $\frac{1}{2}$	s-FCM $\frac{1}{4}$	g-FCM	g-FCM $\frac{1}{2}$	g-FCM $\frac{1}{4}$
1	1.445402	1.445404	1.445404	1.445411	1.445402	1.445412
2	1.445400	1.445378	1.445378	1.445409	1.445391	1.445411
3	1.438605	1.438569	1.438569	1.438635	1.438610	1.438633
4	1.438384	1.438366	1.438365	1.438409	1.438363	1.438417
5	1.438466	1.438442	1.438441	1.438494	1.438471	1.438498
6	1.438466	1.438441	1.438440	1.438490	1.438445	1.438491

# of Nodes	21181	10625	5454	21181	10625	5454

Time (s)	40	16	8	36	16	8