**3. FDTD and FDM algorithms for the MHOF analysis**

Optical light is electromagnetic wave in nature, and thus its propagation properties are governed by the laws of electrodynamics which are collectively known as Maxwell's equations. It is known that guiding of light signals in the MHOF with photonic bandgap structures relies on constructive interference effect due to the periodic arrangement of identical air holes. On the other hand, when the air holes become random in size, location, or both, the phenomenon of total internal reflection takes place and the index difference between the core and cladding provides light confinement and hence guidance of light along the fiber. Although the propagation characteristics of complicated structures like arbitrary MHOFs cannot be calculated easily using analytical methods, there are ways to solve electromagnetic problems numerically.

In this section, two numerical techniques of finite-difference time-domain (FDTD) and finite difference method (FDM) are addressed considering extension to the analysis of holey optical fibers with arbitrary air-hole distributions. Each of these techniques has certain advantages. Using the FDTD method, the continuous electromagnetic field in a finite volume of space is sampled at distinct points in a space lattice and at equally spaced sampling points in time. The sampled data at the points are used for numerical calculations of allowed modes, without generating spurious mode solutions, in a given waveguide. Despite being an effective technique for calculation of propagation constants of guided modes, the FDTD method is not well suited for the evaluation of individual mode field distributions. This is because the source is an impulse function in the time domain covering an infinite spectrum, thus field-distribution solutions are superposition of all possible modes. To alleviate this problem with propagation constants available from FDTD, individual mode field distributions are obtained using the FDM, which can quickly and conveniently provide individual mode field solutions.

The FDTD has gained considerable popularity in recent years, because this method provides robust solutions, based on Maxwell's equations [8], and can readily accommodate complexvalued material properties. An arbitrary material object can be approximated by building up unit cells for which field component positions are disposed with the desired values of permittivity and permeability. Once the geometry of the object is specified in the numerical simulation region, source condition is modelled somewhere in the region. Initially, it is assumed that all fields within the calculation domain are identically zero. Then, an incident wave is enforced to enter the numerical calculation region.

28 Optical Communication

2

(1)

2

*x y dxdy*


**E**

**E**

where **E**(,) *x y* is the electric field on a transverse plane. For the optical waveguide of Figure 1(a), the results of the effective area are 1.7925 m2, 2.0624 m2, and 2.2148 m2 at the

Optical light is electromagnetic wave in nature, and thus its propagation properties are governed by the laws of electrodynamics which are collectively known as Maxwell's equations. It is known that guiding of light signals in the MHOF with photonic bandgap structures relies on constructive interference effect due to the periodic arrangement of identical air holes. On the other hand, when the air holes become random in size, location, or both, the phenomenon of total internal reflection takes place and the index difference between the core and cladding provides light confinement and hence guidance of light along the fiber. Although the propagation characteristics of complicated structures like arbitrary MHOFs cannot be calculated easily using analytical methods, there are ways to

In this section, two numerical techniques of finite-difference time-domain (FDTD) and finite difference method (FDM) are addressed considering extension to the analysis of holey optical fibers with arbitrary air-hole distributions. Each of these techniques has certain advantages. Using the FDTD method, the continuous electromagnetic field in a finite volume of space is sampled at distinct points in a space lattice and at equally spaced sampling points in time. The sampled data at the points are used for numerical calculations of allowed modes, without generating spurious mode solutions, in a given waveguide. Despite being an effective technique for calculation of propagation constants of guided modes, the FDTD method is not well suited for the evaluation of individual mode field distributions. This is because the source is an impulse function in the time domain covering an infinite spectrum, thus field-distribution solutions are superposition of all possible modes. To alleviate this problem with propagation constants available from FDTD, individual mode field distributions are obtained using the FDM, which can quickly and

The FDTD has gained considerable popularity in recent years, because this method provides robust solutions, based on Maxwell's equations [8], and can readily accommodate complexvalued material properties. An arbitrary material object can be approximated by building up unit cells for which field component positions are disposed with the desired values of permittivity and permeability. Once the geometry of the object is specified in the numerical simulation region, source condition is modelled somewhere in the region. Initially, it is

 

 

*eff*

**3. FDTD and FDM algorithms for the MHOF analysis** 

*A*

wavelengths of 0.8 μm, 1.3 μm, and 1.55 μm, respectively.

solve electromagnetic problems numerically.

conveniently provide individual mode field solutions.

(,)

4

*x y dxdy*

Using the MKS system of units, let us first consider Maxwell's curl equations expressed as:

$$
\nabla \times \mathbf{E} = -\mu \frac{d\mathbf{H}}{dt} \tag{2}
$$

$$
\nabla \times \mathbf{H} = \varepsilon \frac{d\mathbf{E}}{dt} \tag{3}
$$

where is the electrical permittivity constant in F/m and is the magnetic permeability constant in H/m. Expanding the curl expressions and equating the like components, the system of six coupled partial differential equations are formed for the FDTD analysis of electromagnetic wave interactions with general three-dimensional objects. It should be noted that the electric and magnetic field components (Ex, Ey, Ez, Hx, Hy, and Hz) are interrelated. That is, Maxwell's equations do not directly yield electric and magnetic field values, but rather relate the rate of change between electric and magnetic field values.

Adopting central finite difference approximation for space and time derivatives with accuracy to the second order, the following approximations as representative examples in a three dimensional (3D) FDTD formulation can be developed:

$$E\_x^{n+1}(i+\frac{1}{2},j,k) = E\_x^n(i+\frac{1}{2},j,k) + \frac{\Delta t}{\varepsilon\_0 c\_r (i+\frac{1}{2},j,k)} \left[ \left[ \frac{H\_z^{n+\frac{1}{2}}(i+\frac{1}{2},j+\frac{1}{2},k) - H\_z^{n+\frac{1}{2}}(i+\frac{1}{2},j-\frac{1}{2},k)}{\Delta y} \right] \right] \tag{4}$$

$$-\left[ \frac{H\_y^{n+\frac{1}{2}}(i+\frac{1}{2},j,k+\frac{1}{2}) - H\_y^{n+\frac{1}{2}}(i+\frac{1}{2},j,k-\frac{1}{2})}{\Delta z} \right] \tag{5}$$

$$\begin{aligned} H\_y^{n+\frac{1}{2}}(i+\frac{1}{2},j,k+\frac{1}{2}) &= H\_y^{n-\frac{1}{2}}(i+\frac{1}{2},j,k+\frac{1}{2}) + \frac{\Delta t}{\mu} \left[ \left[ \frac{E\_z^n(i+1,j,k+\frac{1}{2}) - E\_z^n(i,j,k+\frac{1}{2})}{\Delta x} \right] \right] \\ &- \left[ \frac{E\_x^n(i+\frac{1}{2},j,k+1) - E\_x^n(i+\frac{1}{2},j,k)}{\Delta z} \right] \end{aligned} \tag{5}$$

$$H\_z^{n+\frac{M}{2}}(i+\frac{1}{2},j+\frac{1}{2},k) = H\_z^{n-\frac{M}{2}}(i+\frac{1}{2},j+\frac{1}{2},k) + \frac{\Delta t}{\mu} \left[ \left[ \frac{E\_x^n(i+\frac{1}{2},j+1,k) - E\_x^n(i+\frac{1}{2},j,k)}{\Delta y} \right] \right. \tag{6}$$

$$-\left[ \frac{E\_y^n(i+1,j+\frac{1}{2},k) - E\_y^n(i,j+\frac{1}{2},k)}{\Delta x} \right] \tag{7}$$

where i, j, k, and n are integers for x, y, z, and t, respectively, as the space and time increments [9].

Since optical fibers such as MHOFs generally have no variations along the direction of propagation and variations of material properties are limited to the transverse directions as shown in Figure 3, the 3D FDTD formulation can be simplified to the compact two dimensional (2D) FDTD algorithm [10]. By using phasor notation with the axial propagation constant (), the first-order partial derivatives with respect to z are replaced with -jβ, because the z-dependence of fields is as exp(-jβz). And two adjacent fields required for the first-order derivatives in the discretized space region can be represented by a field at the mid point between them. Based on these two facts, the following formulation as an example is obtained:

$$\mathbf{H}\_{\mathbf{y}}^{n+\frac{\chi'}{2}}(\mathbf{i}+\frac{1}{2},\mathbf{j}) = \mathbf{H}\_{\mathbf{y}}^{n-\frac{\chi'}{2}}(\mathbf{i}+\frac{1}{2},\mathbf{j}) + \frac{\Delta t}{\mu} \left\{ \left[ \frac{\mathbf{E}\_{z}^{n}(\mathbf{i}+\mathbf{1},\mathbf{j}) - \mathbf{E}\_{z}^{n}(\mathbf{i},\mathbf{j})}{\Delta \mathbf{x}} \right] + \mathbf{j}\boldsymbol{\beta} \cdot \mathbf{E}\_{x}^{n}(\mathbf{i}+\frac{1}{2},\mathbf{j}) \right\} \tag{7}$$

The resulting 2D algorithm takes advantage of significant reduction in the required computer memory allocation and running time. Thus, for computer-calculation of arbitrary waveguides that are uniform along the direction of wave propagation, only modelling of the cross sections of waveguides is sufficient.

Along with this efficient algorithm, infinite media in the 2D space for an arbitrary electromagnetic object need to be modelled carefully, because computer memory is limited in the calculation region even with advanced current technology. In order to model regions extending to infinity, a perfectly matched layer (PML) as a highly effective absorbing boundary condition (ABC) is designed at the outer lattice boundary of a calculation domain. Ideally, the absorbing medium is only as thick as a few lattice cells, highly absorbing, reflectionless to all impinging waves, and effective over the full range of operating wavelengths.

**Figure 3.** Schematic of a cross section for an MHOF with one layer of air holes in a hexagonal arrangement

Similarly to the development of the FDTD algorithm, the FDM formulation can be derived from the coupled Maxwell's equations [8]. For continuous waves in linear and isotropic media, combining Eqs (2) and (3) results in the following vectorial wave equation:

30 Optical Communication

is obtained:

arrangement

1 1

cross sections of waveguides is sufficient.

Air

Since optical fibers such as MHOFs generally have no variations along the direction of propagation and variations of material properties are limited to the transverse directions as shown in Figure 3, the 3D FDTD formulation can be simplified to the compact two dimensional (2D) FDTD algorithm [10]. By using phasor notation with the axial propagation constant (), the first-order partial derivatives with respect to z are replaced with -jβ, because the z-dependence of fields is as exp(-jβz). And two adjacent fields required for the first-order derivatives in the discretized space region can be represented by a field at the mid point between them. Based on these two facts, the following formulation as an example

> 2 2 1 1 1 2 2 2 ( 1, ) ( , ) ( ,) ( ,) ( ,)

The resulting 2D algorithm takes advantage of significant reduction in the required computer memory allocation and running time. Thus, for computer-calculation of arbitrary waveguides that are uniform along the direction of wave propagation, only modelling of the

Along with this efficient algorithm, infinite media in the 2D space for an arbitrary electromagnetic object need to be modelled carefully, because computer memory is limited in the calculation region even with advanced current technology. In order to model regions extending to infinity, a perfectly matched layer (PML) as a highly effective absorbing boundary condition (ABC) is designed at the outer lattice boundary of a calculation domain. Ideally, the absorbing medium is only as thick as a few lattice cells, highly absorbing, reflectionless to all

*x*

r

**x** 

**y** 

**z** 

(7)

*n n n n z z <sup>n</sup> y y x <sup>t</sup> E i j E ij H i jH i j j Ei j*

impinging waves, and effective over the full range of operating wavelengths.

**Figure 3.** Schematic of a cross section for an MHOF with one layer of air holes in a hexagonal

**d**

Glass

$$
\nabla \times \nabla \times \mathbf{E} - \nu^2 k\_0^2 \mathbf{E} = 0 \tag{8}
$$

where n is the refractive index and k0 is the propagation constant in free space. Many waveguiding devices, like optical fibers, can be viewed as z-invariant, or piecewise zinvariant structures. For those structures, the refractive index n(x,y,z) varies slowly along the propagation direction z, which is valid for most photonic guided-wave devices. By using the vector identity of <sup>2</sup> , Eq (8) can be written as

$$
\nabla^2 \mathbf{E} + n^2 k\_0^2 \mathbf{E} = \nabla \left( \nabla \cdot \mathbf{E} \right) \tag{9}
$$

Also with the reasonable assumption of negligible time dependency along the z-axis, the FDM formulation as in Eq (9) can be implemented by replacing spatial derivatives with finite difference approximations. Here, it is noted that the transverse component of (9) is

$$
\nabla^2 \mathbf{E}\_t + n^2 k\_0^2 \mathbf{E}\_t = \nabla\_t \left( \nabla\_t \cdot \mathbf{E}\_t + \frac{\partial E\_z}{\partial z} \right) \tag{10}
$$

where the subscript "t" stands for the transverse components. Since the longitudinal component may be readily obtained by application of the following zero divergence (Gauss's law) constraint:

$$\nabla \cdot \left( n^2 \mathbf{E} \right) = \mathbf{0} \,, \tag{11}$$

the transverse components are sufficient to describe the full-vectorial natures of the electromagnetic field in an optical waveguide.

For the initial investigation of guidance properties of MHOFs, the optical fiber shown in Figure 3 is computer-analyzed. Generally, the MHOF geometry can be described with two parameters, pitch length () and diameter (d), as indicated in Figure 3. Here, the pitch length is the distance between centers of two closest air holes with the cylindrical shape. For the MHOF of Figure 3, each small air hole has a diameter of 1.4 μm, constituting a hexagon with = 1.7 μm. The glass portion surrounding the six air holes of the yellow regions has a refractive index of 1.45. The outer radius (r) of the holey fiber is assumed to be 10 μm. Also, the outside region of the MHOF is air.

Once the cross section of a holey fiber is defined in a proper calculation domain, the FDTD simulation can be undertaken with several specified parameters, such as τ in defining a Gaussian source, t for stable simulation, the total number (ntot) of time steps for sampling data in the time domain, and reasonable values of . Here, in order to avoid numerical divergence and ensure stability of the FDTD algorithm, an appropriate Δt needs to be selected to satisfy the following stability condition:

$$
\Delta t \le \frac{1}{c\_{\text{M}} \left[ \frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} + \frac{1}{\Delta z^2} \right]^{\frac{1}{\Delta x}}} \tag{12}
$$

where cM is the maximum wave phase velocity within a given numerical model. Summarizing the mechanism behind the FDTD analysis, the computer simulation proceeds by the following steps:


Figure 4 illustrates the characteristic curves, which are obtained from the FDTD calculation, for the first three lower-order modes of the MHOF defined in Figure 3.

**Figure 4.** Effective refractive index versus wavelength for the first three modes in an MHOF with one layer of air holes

The red curve with the star symbols plots the normalized propagation constant for the first mode versus wavelength, while the blue and green curves show the normalized propagation constants for the second and the third modes, respectively. The results indicate that the HMOF with a single hexagonal air-hole cladding layer supports multimode guiding.

### **4. Effects of cladding structure on optical guidance properties**

32 Optical Communication

by the following steps:

layer of air holes

1 2 (12)

M 222

where cM is the maximum wave phase velocity within a given numerical model. Summarizing the mechanism behind the FDTD analysis, the computer simulation proceeds

Figure 4 illustrates the characteristic curves, which are obtained from the FDTD calculation,

**Figure 4.** Effective refractive index versus wavelength for the first three modes in an MHOF with one

The red curve with the star symbols plots the normalized propagation constant for the first mode versus wavelength, while the blue and green curves show the normalized propagation constants for the second and the third modes, respectively. The results indicate

*t*

 Choose appropriate parameter values (τ, t, ntot, and ) Data sampling of a field component in the time domain

 Take the Fourier transform of the time data Get spectral data of a field component

Collect and mode frequency data

Pick mode frequencies associated with the value

Make a plot of mode index versus wavelength

*c*

for the first three lower-order modes of the MHOF defined in Figure 3.

1

111

*xyz*

Based on the FDTD and the FDM techniques, important optical guidance properties of the fundamental mode, such as normalized propagation constant, chromatic dispersion, field distributions, and effective area, have been investigated for the HMOF with different structural parameter values. In this section, how the cladding structure of the HMOF affects optical guidance properties is first analyzed. Then in the next section, the investigation is extended to find out the effects of the core in the holey fiber design on optical propagation properties.

For the case when the MHOF has = 2.0 μm, d = 1.2 μm, and r = 12.0 μm with three hexagonal cladding layers of air holes as in Figure 5, the effective refractive index characteristics are computer-calculated by the FDTD analysis approach. Assuming that the core index, which is the same as the surrounding medium of the small air holes, is 1.45 and the geometry of the MHOF is uniform along the z-axis, Figure 6 shows the results for the effective refractive index of the fundamental mode obtained from the FDTD method without accounting for the material dispersion (Dmat). In this analysis, a Gaussian pulse source, which excites the MHOF at one end, is used with the parameter value of τ = 0.25x10- 6. The expression describing the source is given as

$$S\_G = \exp\left[-\frac{\left(\chi^2 + \chi^2\right)}{2\pi^2}\right] \tag{13}$$

where x and y are perpendicular coordinates in the transverse plane on a cross section of the holey fiber. In Figure 6, the 3D plot of the source is also shown inset.

In order to see the effect of the thickness of the dielectric layer surrounding the air holes, the outer radius (r) of the holey fiber is reduced from 12 μm to 9 μm. As shown in Figure 6, no significant change on the propagation constant occurs. The propagation constants for the 12 μm and 9-μm radius MHOFs are almost the same except for the longest wavelength. The reason is that the field of the fundamental mode is nearly entirely confined to a region of 9 μm radius, which is about 6 times the operating wavelength of = 1.55 μm. It should be also noted that both of the MHOFs allow only a single mode over the wavelength range from 0.9 μm to 1.8 μm.

Now, with the same parameter values ( = 2.0 μm and d = 1.2 μm) as in Figure 5, the MHOF is numerically analyzed by changing the number of hexagonal air-hole cladding layers. Since it is demonstrated that increasing the outer radius (r) about 6 times larger than the operating wavelength doesn't affect the propagation characteristics significantly, the dielectric glass material surrounding the air holes has been assumed to have the enough thickness. Here, however, it is assumed that the high-index material of the MHOF is pure fused silica (100 mol % SiO2) with refractive index n1, while the low-index region is air. That

**Figure 5.** Transverse structure of the MHOF (The air holes in the fifth layer from the core are linked by the dashed line, forming a hexagonal ring.)

is, in order to account for the material dispersion effect of the glass [7], the wavelength dependence of the glass refractive index is directly incorporated in the effective index calculation through the following widely-utilized Sellmeier's equation:

$$m\_1(\mathcal{A}) = \left[1 + \sum\_{j=1}^3 \frac{A\_j \mathcal{A}^2}{\mathcal{A}^2 - \mathcal{A}\_j^2}\right]^{1/2} \tag{14}$$

where material constants, Aj and j, for the pure fused silica are given as (A1 = 0.6961663, 1 = 0.0684043), (A2 = 0.4079426, 2 = 0.1162414), and (A3 = 0.8974794, 3 = 9.8961610).

Continuing numerical analysis for the optical guidance properties of the step-index rod fiber case, chromatic dispersion is calculated using the FDM analysis and compared to the exact analytical results. With the effective refractive index results in section 2, the chromatic dispersion (Dch) can be calculated by using the following relationship:

$$D\_{\rm ch} = -\frac{\lambda}{c} \frac{d^2 \overline{\beta}}{d\lambda^2} \tag{15}$$

where c is the speed of light in free space [6]. The comparison of the chromatic dispersion results is shown in Figure 7(a) for the dielectric rod waveguide. Here, the blue curve is from the analytical solution, while the red curve is obtained for the FDM analysis. It is noticed

the dashed line, forming a hexagonal ring.)

**Figure 5.** Transverse structure of the MHOF (The air holes in the fifth layer from the core are linked by

**s b**

is, in order to account for the material dispersion effect of the glass [7], the wavelength dependence of the glass refractive index is directly incorporated in the effective index

> 1 2 2 1 () 1 *<sup>j</sup>*

 

where material constants, Aj and j, for the pure fused silica are given as (A1 = 0.6961663, 1 =

Continuing numerical analysis for the optical guidance properties of the step-index rod fiber case, chromatic dispersion is calculated using the FDM analysis and compared to the exact analytical results. With the effective refractive index results in section 2, the chromatic

> ch 2 *<sup>d</sup> <sup>D</sup>*

where c is the speed of light in free space [6]. The comparison of the chromatic dispersion results is shown in Figure 7(a) for the dielectric rod waveguide. Here, the blue curve is from the analytical solution, while the red curve is obtained for the FDM analysis. It is noticed

*c d* 

1 2 <sup>2</sup> <sup>3</sup>

(14)

(15)

*j j*

2

*A*

 

calculation through the following widely-utilized Sellmeier's equation:

*n*

dispersion (Dch) can be calculated by using the following relationship:

0.0684043), (A2 = 0.4079426, 2 = 0.1162414), and (A3 = 0.8974794, 3 = 9.8961610).

**Figure 6.** Normalized propagation constants of 3-layer MHOFs with different radii (Inset: Gaussian source for the FDTD analysis)

that the two results agree well. This provides additive confidence that the FDM analysis can accurately predict chromatic dispersion in holey fibers.

**Figure 7.** Dispersion results: (a) for the dielectric rod waveguide of Figure 1(a), (b) for the MHOFs with = 2.0 μm, d = 1.2 μm, and different numbers of air-hole cladding layers

Again for the MHOF shown in Figure 6, initially, the number of air-hole layers is chosen to be 5. The background glass material (n1) is assumed to be pure fused silica with enough extension. The blue curve in Figure 7(b) shows the chromatic dispersion as a function of

wavelength for the fundamental mode in the MHOF with five hexagonal cladding layers of air holes. The cyan (light-blue) curve illustrates the waveguide dispersion for the same holey fiber without material dispersion included. Compared with the chromatic dispersion result, the waveguide dispersion varies rather gradually from about 78.985 ps/nmkm to 59.537 ps/nmkm over the wavelength range between 1.0 μm and 1.7 μm. It is noticed that the fivelayer MHOF provides positive chromatic and waveguide dispersions over the same wavelength range.

When the number of hexagonal air-hole cladding layers is changed from 5 to 3, it is noticed that the chromatic dispersion curves for the both cases are very close at the shorter wavelength range but deviation occurs from 1.2 μm. However, with the increased number of air-hole layers to 7, the Dch dispersion changes very close to the case of the five-layer MHOF. And adding more air-hole cladding layers like nine or eleven layers does not significantly change the chromatic dispersion result from the five-layer holey fiber in the operation wavelength of interest for general optical communications. This behaviour is reasonably attributed to the fact that when the operating wavelength is shorter, electromagnetic fields are more confined to the core region and only the core region has a major impact on the optical guidance properties. By comparison, when the operating wavelength is longer, fields spread more to the cladding region and the index profile of the cladding region has more influence on the effective refractive index.

Incorporating the material dispersion for the MHOF with five layers of air holes, the field solution of the fundamental mode is obtained by using the finite difference technique. For instance, Figure 8 illustrates the normalized field pattern for the Ex electric component of the MHOF at the operating wavelength of 1.55 μm. Figure 8(a) demonstrates the top view with the colorbar scale in the right side and Figure 8(b) does the 3D view. As expected, it is observed that most of the energy is confined within the core region. Also for the fundamental mode of the MHOF with three, seven or nine cladding layers of air holes, about the same field shape is maintained at the same wavelength.

**Figure 8.** Normalized field pattern for the Ex electric component of the fundamental mode at = 1.55 μm in the MHOF with five layers of air holes and the parameter values of = 2.0 μm and d = 1.2 μm: (a) top view with the colorbar scale; (b) 3D view

Furthermore, similarly to the results of optical guidance properties for the normalised propagation constant, chromatic dispersion, and field distribution, those for the effective core area also tend to change very little when the number of hexagonal air-hole cladding layers becomes larger than five. In other words, adding more air-hole layers to the cladding beyond a certain point will not significantly affect the propagation constant, dispersion, and, in fact, all guidance properties of the MHOF waveguide.
