A vector boundary matching technique for efficient and accurate determination of photonic bandgaps in photonic bandgap fibers

Liang Dong

doi:10.1364/OE.19.012582

1. Introduction

All-solid photonic bandgap fibers have demonstrated some very promising novel applications recently [1,2]. One area of very interesting developments is in the use of such fibers for distributed spectral filtering. This has proven to be very useful for making efficient high power lasers at wavelengths where conventional fiber lasers have found impossible to operate due to strong gain competitions [2]. Light is guided in a lattice defect within the bandgaps of the cladding lattice in these fibers [3–6]. It is, therefore, very important to be able to determine the bandgaps of the cladding lattice in a fast and accurate fashion for the designs of next generations of fibers with increasing complexities. Conventional electronic bandgaps are found using a propagation model based on Bloch wave theorem, since electrons are moving in the plane of the lattice with certain degree of freedom in the propagation direction. In photonic bandgap fibers, the guided light propagates in a well-defined direction along the axis of the 2D cladding lattice. The spatial modes supported in the 2D cladding lattice determine the photonic bandgaps. The problem is an eigenvalue problem and needs a spatial mode solver. Typically, the modes in an infinite 2D lattice of high index nodes are found for the determination of photonic bandgaps in all-solid photonic bandgap fibers.

Plane-wave decomposition technique has been proven to be a very useful tool for solving this infinite lattice problem [7], with infinite lattice being already implies in its algorithm and implementation. It is, however, very inefficient computationally for the lattices of concern in photonic bandgap fibers. The small physical size of the high index nodes and their steep refractive index boundaries require a broad spatial spectral coverage including very high spatial frequencies (see Fig. 1 for an illustration of a lattice). Recently a new technique using a curvilinear space has been proposed to ease this problem somewhat [8]. Additionally, photonic bandgaps of lattices with nodes of complex refractive index profiles have not been demonstrated so far. There are increasing needs for considering complex refractive index distributions of the nodes. All practical fibers have smoothed-out refractive index distributions, instead of step index profiles. More advanced designs also require the optimizations of refractive index profiles of the nodes. Some recent work on photonic lattice of rings has also pointed in this direction [7].

Fig. 1 Illustration of a photonic lattice.

Download Full Size | PDF

Mode solvers based on multipole algorithm are very accurate and highly efficient for solving discrete circular boundary problems [9,10]. This is largely due to the fact that the Fourier-Bessel expansion series represents very well the scattered fields around a circular scatter. The algorithm has been very successful in determining leakage losses of modes in finite cladding structures with circular scatters [9,10]. The plane-wave decomposition technique cannot find leakage loss due to the implied infinity in the cladding lattice. In theory, multipole algorithm with a very large number of nodes in the cladding can be used to approximate the case of infinite cladding. In practice, this is very costly, because the computational efficiency goes down very quickly with an increased number of nodes. A propagation model based on Bloch theorem and lattice sum has been proposed to find photonic bandgaps in photonic bandgap fibers [11]. It is, however, fundamentally ill-suited for this eigenvalue problem.

Recently, a very interesting approach has been proposed for a quick determination of photonic bandgaps in a 2D photonic lattice [12]. This technique is based on recognizing that only the modes at the edges of the photonic bandgaps need to be determined to know where the bands are. Furthermore, the authors recognized that the edge modes have some well-defined characteristics, namely the transverse fields of the modes determining the upper edge of bands of modes are likely in-phase among neighboring nodes and the fields of the modes determining lower edge of bands of modes are likely anti-phase among the neighboring nodes. The original implementation of this idea used a circular boundary for each unit cell, a very rough approximation. This can introduce significant errors especially in cases of strong couplings among neighboring nodes.

In this work, a boundary matching technique is proposed and demonstrated to take account of the mode symmetry and the real geometry of the lattice for a much improved accuracy while using the efficient Fourier-Bessel expansion series. The much improved efficiency, comparing to the conventional plane wave decomposition techniques, comes for two reasons. The first is that only the modes at the band edges are identified and sought. The second is the use of the Fourier-Bessel expansion for efficient and accurate representation of the field around the circular scatters. Furthermore, it is recognized that, even though, modes that have in-phase transverse field distributions among any two neighboring nodes would have the highest modal index among a specific band of modes, these modes may not actually exist due to the specific geometric constraints of modes and lattice symmetries. The same is true for the modes that have anti-phase transverse field distributions among any two neighboring nodes. These modes would have the lowest modal index among the specific band of modes, but they may not actually exist due to the geometric constrains. Furthermore, these modes may have significant different mode field distributions over neighboring nodes while satisfying in-phase or anti-phase conditions, an aspect that has not been considered previously. These additional aspects need to be considered for accurate determination of bandgaps. We have implemented a technique which, in addition to taking into considerations of these new aspects, also enables determination of bandgaps having nodes with any arbitrary refractive index profiles. It incorporates a transfer matrix technique to take account of the refractive index profile of the nodes. Furthermore, this new technique is implemented in a vector form to enable accurate determinations of photonic bandgaps with more than just few percent of index contrasts, another potentially interesting area of future developments using more exotic glasses.

2. Fields around a circular node with a complex refractive index profile

The complex refractive index of each circular node is represented by N concentric rings (not necessarily equally spaced) of uniform indexes. Using the Fourier-Bessel series for field representations, component of the electric and magnetic fields along the propagation direction z in layer p can be written as [9]

e_{z}^{p} = \sum_{m = - \infty}^{+ \infty} [a_{m}^{p} J_{m} (k_{p} r_{}) + b_{m}^{p} H_{m}^{1} (k_{p} r_{})] e^{i m θ}

H_{z}^{p} = \sum_{m = - \infty}^{+ \infty} [C_{m}^{p} J_{m} (k_{p} r_{}) + D_{m}^{p} H_{m}^{1} (k_{p} r_{})] e^{i m θ}

where a, b, C, D are field coefficients to be determined for each layer; J_m is the Bessel function of the first kind; H¹ _m is the Hankel function of the first kind, r is radial coordination centered at the node center in the transverse plane; and θ is the corresponding azimuthal coordination. The first term of each field representation in Eq. (1) and Eq. (2) can be viewed as incoming wave originated from outside of the node and the second term as outgoing scattered wave from the node as described in [9]. We also have

k_{p}^{2} = k^{2} n_{p}^{2} - β^{2}

where k is vacuum wave vector; n_p is refractive index of layer p; and β is the propagation constant of the mode. It is worth noting that we have allowed β>kn_p in Eq. (3), resulting in imaginary k_pr for J_m and H¹ _m in Eq. (1) and (2) in this case. The fields in Eq. (1) and (2) are still valid in this case considering the transformation K_m(-iz) = (π/2)i^{m + 1}H¹ _m(z) and I_m(-iz) = i^˗mJ_m(z) where I_m and K_m are modified Bessel function of the first and second kind respectively. H in Eq. (2) is a normalized magnetic field component for easier representations. It is connected to regular magnetic field h by

H = \sqrt{\frac{μ_{0}}{ε_{0}}} h

Other field components can be determined from the z field components [13]. Boundary conditions require the continuity of all tangential field components, i.e. z and θ components, at each ring boundary. This produces four equations to be satisfied at each boundary by z and θ components of electric and magnetic fields. This gives us the following relationship deriving all field coefficients of one layer from an adjacent layer.

[\begin{matrix} a_{m}^{p + 1} \\ b_{m}^{p + 1} \\ C_{m}^{p + 1} \\ D_{m}^{p + 1} \end{matrix}] = M_{m}^{p} [\begin{matrix} a_{m}^{p} \\ b_{m}^{p} \\ C_{m}^{p} \\ D_{m}^{p} \end{matrix}]

The innermost circle is denoted by p = 0 and the outmost region is denoted by p = N in our denominations. The details of the transfer matrix are given in the appendix. The first thing to note is that there is no coupling between each set of m-field coefficients. This is due to the circular boundary between layers which does not couple between the m terms. The second thing to note is that there is generally some coupling between electric and magnetic fields at each boundary. It is then easy to see that we have the following relationship between fields outside the node and fields at the center of the node

[\begin{matrix} a_{m}^{N} \\ b_{m}^{N} \\ C_{m}^{N} \\ D_{m}^{N} \end{matrix}] = M_{m} [\begin{matrix} a_{m}^{0} \\ b_{m}^{0} \\ C_{m}^{0} \\ D_{m}^{0} \end{matrix}]

where

M_{m} = \prod_{p = 0}^{N - 1} M_{m}^{p}

Bearing in mind that scattering fields are only generated at the discontinuous boundaries, there is no outgoing wave over the innermost circle, i.e. b_m = D_m = 0 for p = 0. It can now be understood that all the four field coefficients for the outmost region, i.e. a_m, b_m, C_m and D_m for p = N, can be determined by the only two non-zero components a_m and C_m of the innermost circle for each m. Since our following analysis is going to mostly involve the field in the outmost region, we will write them down.

e_{z}^{N} = \sum_{m = - \infty}^{+ \infty} [a_{m}^{N} J_{m} (k_{p} r_{}) + b_{m}^{N} H_{m}^{1} (k_{p} r_{})] e^{i m θ}

H_{z}^{N} = \sum_{m = - \infty}^{+ \infty} [C_{m}^{N} J_{m} (k_{p} r_{}) + D_{m}^{N} H_{m}^{1} (k_{p} r_{})] e^{i m θ}

It needs to be emphasized that there are only two independent coefficients in these field representations in Eq. (8) and Eq. (9). Other field components can be derived from these. The transverse field can then be derived for the region around the node

e_{x}^{N} = e_{r}^{N} {\cos (θ) - e}_{θ}^{N} {\sin (θ)}_{}

e_{y}^{N} = e_{r}^{N} {\sin (θ) + e}_{θ}^{N} {\cos (θ)}_{}

Corresponding magnetic fields can be similarly written.

To implement the boundary matching technique for the edge modes, symmetries of the mode field distributions at the upper and lower edges have to be determined first for each band of modes. This can be done using either a mode solver or alternatively by a mode symmetry analysis. It is worth noting mode symmetry analysis can be complicated, especially considering that fields over adjacent nodes can be dominated by different m-terms in the Fourier-Bessel expansion series and the specific geometric constraints of the lattice. A Multipole mode solver is heavily relied on in this work to be sure of the mode symmetry. It needs to be noted that these edge mode symmetries need to be determined only once for each band of modes for each lattice arrangement. This is because these edge mode symmetries are only dependent on the dominant azimuthal mode number and lattice symmetry. Therefore, as it will be discussed in more details later, each band needs to be determined separately as the mode symmetry is expected to be dependent on the dominant azimuthal mode number.

3. Vector boundary matching technique for the first band of modes

In this section, the upper and lower edge modes of the first band of modes, which are derivatives of the fundamental mode of the node, are used to illustrate the boundary matching technique. Although hexagonal lattice is used throughout this work, this technique can be similarly implemented for any other lattice arrangements.

Symmetries of the edge modes of the first band of modes is characterized by the dominance of m = 1 terms in e_z and h_z fields. Unlike in the scalar mode approximations, other m terms are non-zero, but smaller in comparison. The z components are dominated by two anti-phase symmetric lobes over the node (see Fig. 2 and Fig. 3 ). This gives a transverse field over the node resembles the quasi-Gaussian mode with which we are familiar. The orientations of each z field components need to be such so that the transverse fields are in-phase between adjacent nodes for modes determining the upper edge of this band of modes. This can be easily achieved for this mode as shown in Fig. 2.

Fig. 2 Amplitudes of e_z (a), h_z (b), and e_t (c) for a lattice with node V = 0.9 and d/Λ = 0.19. This is a mode located at upper edge of the first band.

Download Full Size | PDF

Fig. 3 Amplitudes of e_z (a), h_z (b), and e_t (c) for a lattice with node V = 0.8 and d/Λ = 0.17. This is a mode located at lower edge of the first band.

Download Full Size | PDF

The mode in Fig. 2 is determined by a multipole mode solver with a finite number of nodes. The field amplitudes of the e_z, h_z and transverse electric field e_t are given in Fig. 2 for a hexagonal lattice with d/Λ = 0.19 where d is node diameter and Λ is node spacing. The node normalized frequency V is 0.9 in this case, where normalized frequency V is defined as for conventional fibers. In this case, this upper edge mode has an identical distribution over each node. The hexagonal boundary of each unit cell is apparent (see the white cell boundary in Fig. 2). Each corner of this unit cell is located at the center of the equilateral triangle formed by the three adjacent nodes. The transverse field is in-phase between any adjacent nodes. The hexagonal boundary of the center unit cell is illustrated in the e_t plot along with radial lines from center of coordination of this cell. In this case, there are 24 radial lines equally spaced by π/12 degree angle. At the hexagonal cell boundary, we should have ∂\e_t\/∂r = 0 and ∂\h_t\/∂r = 0, due to the symmetry of this mode. This allows us to write a total of 48 boundary equations in this case, two at each intersection point for each radial line with the cell boundary. Bearing in mind that the field parameters around the node can be derived from a_m and C_m at the center of the node for each m from Eq. (6). We can assemble the following cell boundary matching equations.

M [\begin{matrix} \begin{matrix} a_{- M}^{0} \\ ... \\ a_{M}^{0} \end{matrix} \\ \begin{matrix} b_{- M}^{0} \\ ... \\ b_{M}^{0} \end{matrix} \end{matrix}] = 0

Each two rows of M are from ∂\e_t\/∂r = 0 and ∂\h_t\/∂r = 0 respectively at each intersection points, making a total 48 rows. We used ∂\e_x\/∂r = 0 and ∂\h_x\/∂r = 0 in our implementations. This is adequate to constrain all the field coefficients, as e_y and h_y are similarly constrained. The m is chosen from –m_max to m_max, making a total of 2(2m_max + 1) variables. It is worth noting that we are looking for modes with a minimum, i.e. trench, along the cell boundary. ∂\e_x\/∂r = 0 and ∂\h_x\/∂r = 0 are valid even though most of the radial lines are not crossing the cell boundary at a right angle.

For maximum numerical stability, M needs to be close to a rectangular matrix. We have chosen m_max = 12 in our case, making a total of 50 variables in this case, giving M a dimension of 48 × 50. Typically, if a band from the equivalent scalar LP_rs mode needs to be resolved, m_max needs to be chosen so that it is slightly larger than r + 1 to allow accurate resolution of the desired band. The number of radial lines can be chosen accordingly to allow efficient computation as well as accurate resolution of all desired bands. It is not necessary to have equally spaced radial lines. It is also worth noting that the boundary conditions in Eq. (12) lead to coupling among various m terms. This allows the field distributions arising from summation of all m-terms in Eq. (8) and Eq. (9) to produce the desired non-circular field distributions. Search for singular value of M provides solutions for mode propagation constant β from eigenvalue and the associated field coefficients from the corresponding eigenvalue vector, which allow calculation of all fields around a node.

For the lower edge of the first band, a simple symmetry analysis shows that it is impossible to find a field distribution which allows anti-phase between any adjacent nodes in a hexagonal lattice while maintaining a similar power distribution at each node. The field amplitudes given by multipole mode solver for a lower edge mode are given in Fig. 3. It can be seen that there are nodes with zero fields so that the remaining nodes can have anti-phase field between any adjacent nodes with non-zero fields. The unit cell boundary in this case is triangular (see Fig. 3) and its three corners are at the centers of the three nodes with zero field. Due to the field symmetries and the anti-phase nature of field between adjacent nodes, we have e_t = 0 and h_t = 0 on the triangular unit cell boundary of the transverse electric field. We again used 24 radial lines in our case with m_max = 12 to set up a matrix M with a dimension of 48 × 50 to find this band edge.

Since we are expecting the unit cell for edge modes is only dependent on the dominant azimuthal mode number and lattice arrangement, we can then apply the hexagonal and triangular unit cell obtained here for upper and lower band edges respectively for all first band of modes for hexagonal lattices.

4. Determination of higher order bandgaps

The mode symmetry around each node in an infinite lattice of nodes is largely governed by the dominant m terms in the field summations, which is equivalent to the azimuthal mode number r in scalar mode denomination LP_rs, and lattice symmetry. The edge modes of the bands of higher order modes which are dominated by m = 1 term, equivalent to the scalar LP_0s modes, are, therefore, expected to have similar symmetry as the edge modes of the first band for the same lattice which was studied in the last section for a hexagonal lattice. From previous works [12], bandgaps of interests are largely determined by scalar LP_0s and LP_1s modes in a typical photonic bandgap fibers, equivalent to modes dominated by m = 1 and m = 2 terms respectively. Modes with higher azimuthal mode number, i.e. s>1 (equivalent to field being dominated by terms with m>2), tend to be more confined to the nodes. Related bands of these modes tend to be narrow and do not contribute as much to the bandgaps.

The amplitudes of e_z, h_z and e_t of a mode on the upper edge of the second band of modes are given in Fig. 4 for a hexagonal lattice. This mode symmetry is complicated by the geometry constrains of the fields dominated by m = 2 terms. The resulting transverse fields have alternating field distribution from one node to an adjacent node. This makes precise determination of the boundary cell very difficult for this case. We resort to use a similar hexagonal unit cell as that for the upper edge of the first band of modes in this case. The hexagonal unit cell is illustrated in Fig. 4. On this unit cell boundary, we have ∂\e_t\/∂r = 0 and ∂\h_t\/∂r = 0. A total number of 24 equally spaced radial lines (not illustrated) similar to that in last section is used to generate 48 boundary equations in this study.

Fig. 4 Amplitudes of e_z (a), h_z (b), and e_t (c) for a lattice with node V = 2.6 and d/Λ = 0.46. This is a mode located at upper edge of the second band.

Download Full Size | PDF

The amplitudes of e_z, h_z and e_t of a mode on the lower edge of the second band of modes are given in Fig. 5 for a hexagonal lattice of nodes. In this case, the field symmetry is straight forward, resulting in a hexagonal unit cell (see Fig. 5). Transverse fields should be zero over the boundary of this unit cell due to the anti-phase symmetries. A total number of 24 equally spaced radial lines (not illustrated) similar to that in last section is used to generate 48 boundary equations in this study. The mode symmetry of the analysis in this section equally applied to all other modes dominated by m = 2 terms in the field summations as it has been pointed out earlier in the discussion of modes dominated by m = 1 terms.

Fig. 5 Amplitudes of e_z (a), h_z (b), and e_t (c) for a lattice with node V = 2.65 and d/Λ = 0.47. This is a mode located at lower edge of the second band.

Download Full Size | PDF

5. Impact of unit cell boundaries of the edge modes

Since a great deal of efforts has been spent so far in determining the unit cell of the edge modes, some discussion of the impact of the choice of unit cell on the precision in determination of the bandgap edges is warranted. A triangular and hexagonal unit cell have been determined in the last section for the lower edge modes for modes dominated by m = 1 and m = 2 terms respectively in the z field summations. In order to determine the impact of the choice of unit cell on the band edges, the lower band edges of modes dominated by all m-terms for kΛ<180 are determined using both triangular and hexagonal unit cells for the fiber simulated in [12]. The results are shown in Fig. 6 , where the edges determined by the triangular unit cell are given in black dots and those by hexagonal unit cell red circles.

Fig. 6 Lower band edges determined using the triangular unit cell (black dots) and hexagonal unit cell (red circles) for the fiber in [12] (d/Λ = 0.41, node index n₀ = 1.48716 and background index n_b = 1.458.).

Download Full Size | PDF

The two unit cells provide very close results for the upper part of the edges, i.e. regime of low coupling among adjacent nodes. This part of the band edge happens to be an area of strong interests. This is where most relevant bandgaps are for typical photonic bandgap fibers. There are strong deviations especially for lower part of the edges where there is strong coupling among adjacent nodes.

6. Full bandgap simulations

The bandgap diagram for the fiber studied in [12] is simulated as a demonstration of the new technique. This fiber is well studied by the plane wave decomposition technique as well as by the approximation method used in [12]. This provides an easy point of comparison. The hexagonal unit cell described earlier is exclusively used for all the upper edges. The lower edges for band dominated by m = 1 terms are determined by the triangular unit cell and all the other lower edges are determined by the hexagonal unit. The results are plotted in Fig. 7 . The white area is the bandgaps where defect modes are guided while the light blue shaded area is where bands of modes can be guided in the infinite lattice. The hexagonal unit cell used for the upper band edges sometime does find additional in-band modes. Since the modes on upper edge always have the highest modal index, it is easy to recognize the edge modes.

Fig. 7 Bandgap diagram determined by the vector boundary matching technique the fiber in [12] (d/Λ = 0.41, node index n₀ = 1.48716 and background index n_b = 1.458.). The upper band edges are determined by the hexagonal unit cell (black dots) and lower edges are determined by a combination of the triangular and hexagonal unit cells (red circles).

Download Full Size | PDF

A comparison with the bandgap diagram in [12] simulated with plane wave technique for the same fiber indicates a good agreement and an improved accuracy over the approximate method used in [12]. A precise comparison is difficult due to the complexity of the bandgap diagram. The bandgap diagram from the plane wave technique is a mode density plot and edges do not always show up clearly. The location of the bandgaps is, however, reproduced with excellent accuracy.

7. Bandgap of lattice with nodes of graded index

In order to illustrate the ability to calculate bandgaps for lattice with nodes of arbitrary refractive index profiles, hexagonal lattices with nodes of step index and parabolic index profile are simulated. Both nodes have peak refractive index difference of 0.003. In case of the parabolic profile, the refractive index of the nodes is simulated with N = 50 layers of rings of uniform indexes. The nodes with the parabolic index also have strong practical interest, since most fabricated all-solid photonic bandgap fibers are made with nodes of similar parabolic index profile originally intended for graded index multimode fibers. The results are shown in Fig. 8(a) and 8(b) respectively for lattices with the step index nodes and parabolic index nodes. The horizontal axis is V/π, where V is the normalized frequency of the node. This makes sense, because the modes in the lattice are largely determined by the modes in the nodes. The normalization factor π is meant to show that the bandgaps are roughly spaced by ΔV = π/2. The vertical axis is plotted as Δn = n-n_b, where n is modal index and n_b is the index of the background glass. Waveguide strength for guided defect modes is directly related to -Δn of the upper band edge defining the relevant bandgap and can, therefore, be estimated directly. The horizontal axis in the plot for the parabolic index nodes in Fig. 8(b) is rescaled so that the lower edge of the second band crosses Δn = n-n_b = 0 at the same V value as that in Fig. 8(a) for the step index nodes. Hexagonal unit cell is used for both upper and lower band edges in this case. One notable difference in bandgap diagram for a hexagonal lattice of nodes with the parabolic index profile in Fig. 8(b) is that the bands of modes dominated by m>2 terms have moved closer at the upper part of the diagram towards nearest band of modes dominated by m = 1 and m = 2 terms. The bandgaps are in general shallower for the nodes with parabolic index profiles.

Fig. 8 Bandgap diagrams for hexagonal lattices with (a) nodes of a step index profile and (b) nodes of a parabolic index profile. Peak index difference is 0.003, d/Λ = 0.4, n_b = 1.45, M = 6, and 12 equally spaced radial lines. Upper band edges are marked by black dots and lower band edges red circles.

Download Full Size | PDF

8. Summary

A vector boundary matching algorithm has been formulated and demonstrated for efficient and accurate determination of photonic bandgaps for photonic bandgap fibers with circular nodes. The algorithm also incorporates a transfer matrix method for dealing with nodes with arbitrary refractive index profile. The much improved accuracy, comparing to the earlier works, comes from the use of more accurate cell boundary for the implementation of the boundary matching. These cell boundaries need to be determined for each band modes and each lattice configurations, in addition to whether it is upper or lower band edge. Some examples are given to illustrate the technique. This new technique is a useful tool for exploring more advanced designs in the next generation of photonic bandgap fibers.

9. Appendix

The elements for the transfer matrix can be obtained by first deriving other field terms from the z field components for layer p [13].

e_{θ}^{p} = \sum_{m = - \infty}^{\infty} {- \frac{{m β}_{}^{}}{{r k}_{p}^{2}} [a_{m}^{p} J_{m} (k_{p} r_{}) + b_{m}^{p} H_{m}^{1} (k_{p} r_{})] - \frac{{i k}_{}}{k_{p}} [C_{m}^{p} \frac{J_{m - 1} (k_{p} r_{}) - J_{m + 1} (k_{p} r_{})}{2_{}} + D_{m}^{p} \frac{H_{m}^{1} (k_{p} r_{}) - H_{m + 1}^{1} (k_{p} r_{})}{2_{}}]} e^{i m θ}

H_{θ}^{p} = \sum_{m = - \infty}^{\infty} {- \frac{{m β}_{}^{}}{{r k}_{p}^{2}} [C_{m}^{p} J_{m} (k_{p} r_{}) + D_{m}^{p} H_{m}^{1} (k_{p} r_{})] + \frac{{i k n}_{p}^{2}}{k_{p}} [a_{m}^{p} \frac{J_{m - 1} (k_{p} r_{}) - J_{m + 1} (k_{p} r_{})}{2_{}} + b_{m}^{p} \frac{H_{m}^{1} (k_{p} r_{}) - H_{m + 1}^{1} (k_{p} r_{})}{2_{}}]} e^{i m θ}

e_z, H_z, e_θ and H_θ are continuous at each boundary p, resulting in four equations for adjacent layers p and p + 1 for each m. The transfer matrix can then be obtained. The following matrix elements are for a specific m; but subscript m is omitted for simplicity.

M_{11}^{p} = \frac{1}{M_{0}} {J_{m} (k_{p} ρ_{p}) [H_{m - 1}^{1} (k_{p + 1} ρ_{p}) - H_{m + 1}^{1} (k_{p + 1} ρ_{p})] - \frac{n_{p}^{2} k_{p + 1}}{n_{p + 1}^{2} k_{p}} H_{m}^{1} (k_{p + 1} ρ_{p}) [J_{m - 1} (k_{p} ρ_{p}) - J_{m + 1} (k_{p} ρ_{p})]}

M_{12}^{p} = \frac{1}{M_{0}} {H_{m}^{1} (k_{p} ρ_{p}) [H_{m - 1}^{1} (k_{p + 1} ρ_{p}) - H_{m + 1}^{1} (k_{p + 1} ρ_{p})] - \frac{n_{p}^{2} k_{p + 1}}{n_{p + 1}^{2} k_{p}} H_{m}^{1} (k_{p + 1} ρ_{p}) [H_{m - 1}^{1} (k_{p} ρ_{p}) - H_{m + 1}^{1} (k_{p} ρ_{p})]}

M_{13}^{p} = - \frac{1_{}}{M_{0}} \frac{{i 2 m β k}_{p + 1}}{ρ_{p} n_{p + 1}^{2} k_{}} (\frac{1_{}}{k_{p}^{2}} - \frac{1_{}}{k_{p + 1}^{2}}) H_{m}^{1} (k_{p + 1} ρ_{p}) J_{m} (k_{p} ρ_{p})

M_{14}^{p} = - \frac{1_{}}{M_{0}} \frac{{i 2 m β k}_{p + 1}}{ρ_{p} n_{p + 1}^{2} k_{}} (\frac{1_{}}{k_{p}^{2}} - \frac{1_{}}{k_{p + 1}^{2}}) H_{m}^{1} (k_{p + 1} ρ_{p}) H_{m}^{1} (k_{p} ρ_{p})

M_{21}^{p} = - \frac{1_{}}{M_{0}} {J_{m} (k_{p} ρ_{p}) [J_{m - 1} (k_{p + 1} ρ_{p}) - J_{m + 1} (k_{p + 1} ρ_{p})] - \frac{n_{p}^{2} k_{p + 1}}{n_{p + 1}^{2} k_{p}} J_{m} (k_{p + 1} ρ_{p}) [J_{m - 1} (k_{p} ρ_{p}) - J_{m + 1} (k_{p} ρ_{p})]}

M_{22}^{p} = - \frac{1_{}}{M_{0}} {H_{m}^{1} (k_{p} ρ_{p}) [J_{m - 1} (k_{p + 1} ρ_{p}) - J_{m + 1} (k_{p + 1} ρ_{p})] - \frac{n_{p}^{2} k_{p + 1}}{n_{p + 1}^{2} k_{p}} J_{m} (k_{p + 1} ρ_{p}) [H_{m - 1}^{1} (k_{p} ρ_{p}) - H_{m + 1}^{1} (k_{p} ρ_{p})]}

M_{23}^{p} = \frac{1_{}}{M_{0}} \frac{{i 2 m β k}_{p + 1}}{ρ_{p} n_{p + 1}^{2} k_{}} (\frac{1_{}}{k_{p}^{2}} - \frac{1_{}}{k_{p + 1}^{2}}) J_{m} (k_{p + 1} ρ_{p}) J_{m} (k_{p} ρ_{p})

M_{24}^{p} = \frac{1_{}}{M_{0}} \frac{{i 2 m β k}_{p + 1}}{ρ_{p} n_{p + 1}^{2} k_{}} (\frac{1_{}}{k_{p}^{2}} - \frac{1_{}}{k_{p + 1}^{2}}) J_{m} (k_{p + 1} ρ_{p}) H_{m}^{1} (k_{p} ρ_{p})

M_{31}^{p} = \frac{1_{}}{M_{0}} \frac{{i 2 m β k}_{p + 1}}{ρ_{p} k_{}} (\frac{1_{}}{k_{p}^{2}} - \frac{1_{}}{k_{p + 1}^{2}}) H_{m}^{1} (k_{p + 1} ρ_{p}) J_{m} (k_{p} ρ_{p})

M_{32}^{p} = \frac{1_{}}{M_{0}} \frac{{i 2 m β k}_{p + 1}}{ρ_{p} k_{}} (\frac{1_{}}{k_{p}^{2}} - \frac{1_{}}{k_{p + 1}^{2}}) H_{m}^{1} (k_{p + 1} ρ_{p}) H_{m}^{1} (k_{p} ρ_{p})

M_{33}^{p} = \frac{1_{}}{M_{0}} {J_{m} (k_{p} ρ_{p}) [H_{m - 1}^{1} (k_{p + 1} ρ_{p}) - H_{m + 1}^{1} (k_{p + 1} ρ_{p})] - \frac{k_{p + 1}}{k_{p}} H_{m}^{1} (k_{p + 1} ρ_{p}) [J_{m - 1} (k_{p} ρ_{p}) - J_{m + 1} (k_{p} ρ_{p})]}

M_{34}^{p} = \frac{1_{}}{M_{0}} {H_{m}^{1} (k_{p} ρ_{p}) [H_{m - 1}^{1} (k_{p + 1} ρ_{p}) - H_{m + 1}^{1} (k_{p + 1} ρ_{p})] - \frac{k_{p + 1}}{k_{p}} H_{m}^{1} (k_{p + 1} ρ_{p}) [H_{m - 1}^{1} (k_{p} ρ_{p}) - H_{m + 1}^{1} (k_{p} ρ_{p})]}

M_{41}^{p} = - \frac{1_{}}{M_{0}} \frac{{i 2 m β k}_{p + 1}}{ρ_{p} k_{}} (\frac{1_{}}{k_{p}^{2}} - \frac{1_{}}{k_{p + 1}^{2}}) J_{m} (k_{p + 1} ρ_{p}) J_{m} (k_{p} ρ_{p})

M_{42}^{p} = - \frac{1_{}}{M_{0}} \frac{{i 2 m β k}_{p + 1}}{ρ_{p} k_{}} (\frac{1_{}}{k_{p}^{2}} - \frac{1_{}}{k_{p + 1}^{2}}) J_{m} (k_{p + 1} ρ_{p}) H_{m}^{1} (k_{p} ρ_{p})

M_{43}^{p} = - \frac{1_{}}{M_{0}} {J_{m} (k_{p} ρ_{p}) [J_{m - 1} (k_{p + 1} ρ_{p}) - J_{m + 1} (k_{p + 1} ρ_{p})] - \frac{k_{p + 1}}{k_{p}} J_{m} (k_{p + 1} ρ_{p}) [J_{m - 1} (k_{p} ρ_{p}) - J_{m + 1} (k_{p} ρ_{p})]}

M_{44}^{p} = - \frac{1_{}}{M_{0}} {H_{m}^{1} (k_{p} ρ_{p}) [J_{m - 1} (k_{p + 1} ρ_{p}) - J_{m + 1} (k_{p + 1} ρ_{p})] - \frac{k_{p + 1}}{k_{p}} J_{m} (k_{p + 1} ρ_{p}) [H_{m - 1}^{1} (k_{p} ρ_{p}) - H_{m + 1}^{1} (k_{p} ρ_{p})]}

where

M_{0} = J_{m} (k_{p + 1} ρ_{p}) [H_{m - 1}^{1} (k_{p + 1} ρ_{p}) - H_{m + 1}^{1} (k_{p + 1} ρ_{p})] - H_{m}^{1} (k_{p + 1} ρ_{p}) [J_{m - 1} (k_{p + 1} ρ_{p}) - J_{m - 1} (k_{p + 1} ρ_{p})]

The e_θ and H_θ for p = N describe the azimuthal fields around the node and can be derived similarly to these in Eq. (14). Radial field components e_r and H_r for p = N can be obtained through the following equations.

e_{r}^{N} = \sum_{m = - \infty}^{\infty} {- \frac{{m k}_{}^{}}{{r k}_{b}^{2}} [C_{m}^{N} J_{m} (k_{b} r_{}) + D_{m}^{N} H_{m}^{1} (k_{b} r_{})] + \frac{{i β}_{}}{k_{b}} [a_{m}^{N} \frac{J_{m - 1} (k_{b} r_{}) - J_{m + 1} (k_{b} r_{})}{2_{}} + b_{m}^{N} \frac{H_{m}^{1} (k_{b} r_{}) - H_{m + 1}^{1} (k_{b} r_{})}{2_{}}]} e^{i m θ}

H_{r}^{N} = \sum_{m = - \infty}^{\infty} {\frac{{m k n}_{b}^{2}}{{r k}_{b}^{2}} [a_{m}^{N} J_{m} (k_{b} r_{}) + b_{m}^{N} H_{m}^{1} (k_{b} r_{})] + \frac{{i β}_{}^{}}{k_{b}} [C_{m}^{N} \frac{J_{m - 1} (k_{b} r_{}) - J_{m + 1} (k_{b} r_{})}{2_{}} + D_{m}^{N} \frac{H_{m}^{1} (k_{b} r_{}) - H_{m + 1}^{1} (k_{b} r_{})}{2_{}}]} e^{i m θ}

The transverse fields around the node can then be obtained using Eq. (10) and Eq. (11). The boundary conditions at each intersection of the radial lines with the unit cell boundary can be obtained so that matrix M for the linear equations in Eq. (12) can be set up.

References and links

1. C. B. Olausson, C. I. Falk, J. K. Lyngsø, B. B. Jensen, K. T. Therkildsen, J. W. Thomsen, K. P. Hansen, A. Bjarklev, and J. Broeng, “Amplification and ASE suppression in a polarization-maintaining ytterbium-doped all-solid photonic bandgap fibre,” Opt. Express 16(18), 13657–13662 (2008). [CrossRef] [PubMed]

2. A. Shirakawa, C. B. Olausson, H. Maruyama, K. Ueda, J. K. Lyngsø, and J. Broeng, “High power ytterbium fiber lasers at extremely long wavelengths by photonic bandgap fiber technology,” Opt. Fiber Technol. 16(6), 449–457 (2010). [CrossRef]

3. F. Luan, A. K. George, T. D. Hedley, G. J. Pearce, D. M. Bird, J. C. Knight, and P. St. J. Russell, “All-solid photonic bandgap fiber,” Opt. Lett. 29(20), 2369–2371 (2004). [CrossRef] [PubMed]

4. A. Argyros, T. A. Birks, S. G. Leon-Saval, C. M. B. Cordeiro, F. Luan, and P. St. J. Russell, “Photonic bandgap with an index step of one percent,” Opt. Express 13(1), 309–314 (2005). [CrossRef] [PubMed]

5. A. Argyros, T. A. Birks, S. G. Leon-Saval, C. M. B. Cordeiro, and P. St J Russell, “Guidance properties of low-contrast photonic bandgap fibres,” Opt. Express 13(7), 2503–2511 (2005). [CrossRef] [PubMed]

6. G. Bouwmans, L. Bigot, Y. Quiquempois, F. Lopez, L. Provino, and M. Douay, “Fabrication and characterization of an all-solid 2D photonic bandgap fiber with a low-loss region (< 20 dB/km) around 1550 nm,” Opt. Express 13(21), 8452–8459 (2005). [CrossRef] [PubMed]

7. J. M. Stone, G. J. Pearce, F. Luan, T. A. Birks, J. C. Knight, A. K. George, and D. M. Bird, “An improved photonic bandgap fiber based on an array of rings,” Opt. Express 14(13), 6291–6296 (2006). [CrossRef] [PubMed]

8. G. J. Pearce, T. D. Hedley, and D. M. Bird, “Adaptive curvilinear coordinates in a plane-wave solution of Maxwell’s equations in photonic crystals,” Phys. Rev. B 71(19), 195108 (2005). [CrossRef]

9. T. P. White, B. T. Kuhlmey, R. C. McPhedran, D. Maystre, G. Renversez, C. M. de Sterke, and L. C. Botten, “Multipole method for microstructured optical fibers, I. formulation,” J. Opt. Soc. Am. B 19(10), 2322–2330 (2002). [CrossRef]

10. B. T. Kuhlmey, T. P. White, G. Renversez, D. Maystre, L. C. Botten, C. M. de Sterke, and R. C. McPhedran, “Multipole method for microstructured optical fibers, II. implementation and results,” J. Opt. Soc. Am. B 19(10), 2331–2340 (2002). [CrossRef]

11. T. White, R. McPhedran, L. Botten, G. Smith, and C. M. de Sterke, “Calculations of air-guided modes in photonic crystal fibers using the multipole method,” Opt. Express 9(13), 721–732 (2001). [CrossRef] [PubMed]

12. T. A. Birks, G. J. Pearce, and D. M. Bird, “Approximate band structure calculation for photonic bandgap fibres,” Opt. Express 14(20), 9483–9490 (2006). [CrossRef] [PubMed]

13. A. W. Snyder, and J. D. Love, Optical Waveguide Theory (Chapman and Hall, 1983).

A vector boundary matching technique for efficient and accurate determination of photonic bandgaps in photonic bandgap fibers

Abstract

1. Introduction

2. Fields around a circular node with a complex refractive index profile

3. Vector boundary matching technique for the first band of modes

4. Determination of higher order bandgaps

5. Impact of unit cell boundaries of the edge modes

6. Full bandgap simulations

7. Bandgap of lattice with nodes of graded index

8. Summary

9. Appendix

References and links

Cited By

Figures (8)

Equations (33)

Optics Express