1. Introduction

The Boundary Element Method (BEM) has been shown to be accurate and rigorous in solving two-dimensional optical waveguide problems [1]. Guan et al. [2,3] has applied BEM to the study of photonic crystal fiber (PCF), which is a single-material fiber with air holes in the cladding running through its length. The nature of BEM is such that it requires the integration along the boundary only. This implies that meshing is extremely easy for two-dimensional problems. Furthermore, good accuracy is achieved with extremely small number of unknowns. Another distinct advantage of BEM is that external absorbing boundary conditions (ABC) are no longer needed for open boundary problems. However, a non-standard eigenvalue problem is encountered in the BEM. A direct search method is required. It needs to find the distribution of the determinant of the coefficient matrix in order to obtain the wave number k. During the direct search process, the integrals have to be re-evaluated. Consequently, the scheme is very time-consuming.

We propose the use of the multiple reciprocity boundary element method (MRBEM) to analyze PCF. In MRBEM, the Helmholtz equation is interpreted as a Poisson equation with a solution dependent body force term. The method can effectively separate k from the integration by using a series of fundamental solutions of the Laplace equation instead of the usual fundamental solution of the Helmholtz equation. Consequently, the efficiency in root searching is significantly improved.

2. MRBEM formulation

The PCF is taken to be uniform in the propagation direction. Thus the Helmholtz equation to be solved in each homogenous region is of the following form,

where ${\nabla }_t^2$ is the Laplacian operator in transverse plane.

where k₀ and β are free space wave number and propagation constant respectively.

In the direct BEM, by applying the Green’s second identity, the field at the observation point r′ can be expressed by the boundary integral as:

where r is the field point, n is the normal outward direction to the boundary and u*₀ is the zero order fundamental solution of the two dimensional Helmholtz equation given by,

Here $H_0^{\left[2\right]}$ is zero order second kind Hankel function with the unknown k and |r′-r| (distance between source point and field point) as variables. Here k appears in the argument of the Hankel function.

The value of c(r′) is 1 when the field point r is inside the domain V, 0 when r is outside the domain and 1/2 when r is on the boundary, which is smooth in the PCF problem.

During the root-searching process, different values of k are substituted into the Green’s function and the boundary integrals involved are evaluated again and again, making the search very time consuming.

In the MRBEM, the Helmholtz equation is first interpreted as a Poisson equation with a solution dependent body force term given as follows:

where b₀ (x)=k ² u(x). Thus Eq. (2) becomes [4],

where

$u_0^{*}(r\prime ,r)=-\frac{1}{2\pi }\ln \left(\mid r-r\prime \mid \right).$

To evaluate the domain integral on the right hand side of Eq. (4), the domain may be discretized into cells in direct boundary element method. However, the attractive ‘boundary only’ feature of the BEM, as mentioned earlier, is lost.

The MRBEM provides a method for transforming the domain integral into boundary integrals. It introduces a sequence of higher order fundamental solutions as well as a sequence of body force Laplacians by the recurrence formulae:

where j is the order of fundamental solutions which has the value of 0, 1 and so on [4]. With this, the domain integral on the right hand side of Eq. (4) becomes,

Fortunately, as m increase, the last term in Eq. (5) converges to zero. Thus when a sufficiently large m is used to truncate the number of terms, Eq. (4) becomes, [5]

where u* _j is jth order fundamental solution given by,

where

$s_j=\{\begin{array}{cc}0& \left(j=0\right)\\ {\displaystyle \sum _{n=1}^j}\frac{1}{n}& \left(j\ge 1\right)\end{array}$

Equation (6) can be discretized to give,

It is noteworthy that matrices [H _j] and [G_j ] (j=0,1,…m) are evaluated once only because the fundamental solution (Eq. (7)) is independent of k.

Two coefficient matrices can be further expressed as follows:

Equation (8) can be simplified as,

Equation (9) is the eigenvalue equation for each homogeneous region, which is exactly of the same form as in the direct BEM. Here the matrix H and G are both real matrix instead of complex matrix in the case of BEM. This is because the fundamental solutions of Laplace equation are of real form instead of the complex fundamental solution for Helmholtz equation. By comparing the final matrix of BEM and MRBEM, one can find that MRBEM is equivalent to BEM but without the imaginary part. By doing so, the time taken to evaluate the integrations is greatly reduced, however, the imaginary part of the root will not be found.

By applying the boundary conditions of continuity of transverse magnetic field components (H_x and H_y ) and longitudinal field components (E_z and H_z ), an eigenvalue equation of the following form can be obtained which is the same as in direct BEM [2].

where x is the unknown field values at the boundaries, and A is the coefficient matrix with unknown propagation constant.

To find non-trivial solutions of x, different values of k are substituted into A. The distribution of the determinant of A is plotted against k. A rapid decrease in the value of the determinant indicates the presence of a particular value of k, which gives Eq. (10) a non-trivial solution. Subsequently, the field components at all the boundary segments are found by substituting k back into Eq. (10). The field components at any point of the domain can be obtained by Eq. (6) in a straightforward way.

3. Simulation results

Because direct BEM has been shown to be very accurate [1][2,3][6], the validity of MRBEM is verified by comparing its results with those obtained by the direct BEM for a 3-ring holey fiber. Figure 1 shows the meshing of the air holes boundaries. The air holes are segmented into 32, 16 and 8 segments respectively from the inner to the outer ring.

 

Fig. 1. 3-ring holey fiber with radius=0.4µm and Λ=2.3µm

Download Full Size | PDF

In MRBEM, it is critical to set the number of terms m correctly for different regions. Results obtained with different values of m at wavelength λ=1.55µm are tabulated in Table 1.

Table 1. Root searching results with different values of m

View Table | View all tables in this article

Notice that the value of m for silica region is higher than the one for air region. This is because the maximum distance r between elements in silica region is much larger than that in air region, hence more terms are needed in order to make u* _j in Eq. (7) converges to very small values. From Table 1, it is found the results converge with the increase of m.

The results obtained using MRBEM and BEM for the same structure with different meshes at λ=1.55µm are compared in table 2. The results by BEM are in complex form and its accuracy has been verified with multipole method [7]. The number of segments used to discretize the holes at the 1^st, 2^nd and 3^rd ring are represented by n₁ , n₂ and n₃ respectively. From Table 2, the results of the MRBEM agree well with the ones of the direct BEM. The number of terms m used for silica and air region are 20 and 7 respectively in the MRBEM.

Table 2. Comparison between the MRBEM and direct BEM

View Table | View all tables in this article

The root searching figures in the real domain by MRBEM and direct BEM are shown in Fig. 2. For BEM and MRBEM, n₁ , n₂ and n₃ are 32, 16, and 8 respectively. In the MRBEM, m is 20 and 7 for silica and air region respectively. The root is present at the local minimum of the determinant of matrix A. The insets in Fig. 2 are the surface plots of x component of the magnetic field by the corresponding methods.

 

Fig. 2. Root searching figure: (a) direct BEM (b) MRBEM

Download Full Size | PDF

With the obtained field profile, the effective mode areas (A_eff ) can be calculated using

The values of A_eff calculated based on the mode profiles obtained with the MRBEM and direct BEM in Fig. 2 are 11.2115 µm² and 11.1686 µm² respectively.

A two-ring birefringent PCF shown in Fig. 3 is also studied using the direct BEM and MRBEM. The radii of the small and big air holes are 0.4µm and 0.8 µm respectively.

 

Fig. 3. Birefringent PCF with Λ=2.3µm

Download Full Size | PDF

The birefringence at λ=1.55µm obtained by BEM and MRBEM are 7.4×10^-4 and 7.3×10^-4 respectively. Figure 4 shows the vector field plots of the x and y polarized modes.

 

Fig. 4. Vector plots of the fields as calculated using MRBEM: (a) $H_x^{11}$ (b) $H_y^{11}$

Download Full Size | PDF

Moreover, the performance of the MRBEM in terms of simulation speed is evaluated against direct BEM as shown in Fig. 5. The number of terms m in Eq. (5) used to express the integral in silica and air region are 15 and 5 respectively. It is obvious that the MRBEM outperforms the direct BEM in terms of simulation speed. It is also noteworthy that the speed of MRBEM can be further increased by putting a remainder term when the series starts converging instead of using very large number of terms [5].

 

Fig. 5. Comparison of computation time

Download Full Size | PDF

4. Conclusion

The MRBEM is found to be an efficient and accurate numerical tool for finding the eigenvalues of the Helmholtz equation. By using the fundamental solutions of the Laplace equation, the unknown wave number k is separated from the integrand. The integrations are evaluated once only, which significantly improves the efficiency of the algorithm. Although only a real propagation constant is considered and MRBEM cannot be used to predicate the leakage loss due to a finite number of rings, it is effective in studying other PCF properties such as birefringence and dispersion as well as effective mode area. In conclusion, when confinement loss is not a design concern, MRBEM is a much faster approach to analyze the properties (e.g., dispersion, birefringence and nonlinearity) of PCFs.