1. Introduction

Beside ordinary fibers that guide by total internal reflection, fibers with solid claddings employing the photonic band gap effect (PBG) allow the efficient guidance of light in low refractive index (RI) cores mostly in the fundamental mode [1]. The modal properties of such fibers mainly depend on the used materials and geometrical parameters, and have led to numerous applications in, for instance, high-extinction spectral filtering [2], mode-locked lasing [3], light amplification [4], high power transmission [5] and pulse dispersion engineering [6,7].

One important class of PBG fibers are composed of hexagonal arrays of longitudinal invariant cylindrical glass strands, which are made from a high RI material and located inside a low RI matrix, commonly referred as all-solid PBG (APBG) fibers (schematic shown in Fig. 1(a)). A core section providing transverse light confinement is formed by omitting one or more strands in the center of the structure, which creates one or more states in the PBG and is known as the defect core. Examples of APBG fibers include arrays of GeO₂ doped silica strands in pure silica [8,9] or combinations of various soft glasses or even metal nanowires [10–15]. Moreover very high RI-contrast APBG fibers consisting of a tellurite or chalcogenide glass inside silica or other soft glasses have been implemented, showing, for instance, extremely high extinction ratios of more than 60 dB/cm [16,17]. Pivotal parameters characterizing all mentioned APBG fibers, namely the strand diameters d and pitches Λ (center-to-center distance between the strands) are of the order of several micrometers, which are dimensions that can be easily accessed with current fabrication technology such as fiber drawing [1] or pressure-assisted melt filling [2,18].

 

Fig. 1 (a) Schematic of an all-solid photonic band gap fiber with a hexagonal array of high index rods in a low index matrix and a central defect core (dark green: high refractive index cylindrical strands, bright green: host material). (b) Flow chart indicating how the binary function is obtained. The individual steps are (i) approximating the hexagonal unit cell by a circle, (ii) deriving the dispersion functions of the band edge modes (the two modes are indicated by the purple and dark yellow line) and (iii) plotting the binary function in the parameter domain of interest.

Download Full Size | PDF

A key feature of APBG fibers is that the modes guided in the HI strands interact, i.e. couple together inside the low RI material and form supermodes (SMs). These SMs can phase-match to a central defect mode in case their effective mode indices n_eff are identical, opening a loss channel for the electromagnetic energy to laterally tunnel into the cladding. The phase-matching leads to strong transmission dips, which are commonly referred as resonances. It is important to note that even outside these resonances the guided mode is intrinsically leaky (i.e., continuously sheds energy out of the domain of the defect), since the strand array is finite, i.e. consists of a finite number of strand rings, which imposes core-cladding reflection coefficient below unity. This may suggest high propagation loss, whereas the simulations have revealed that the PBG effect is so strong that using four to six strand rings is sufficient to provided losses down to the dB/km level [19].

The prerequisite to obtain low propagation loss thus defines the key challenge in the field of APBG fibers, which is to identify that parameter space where no cladding SMs exist, i.e., to find the domains of zero density of cladding states (zDOS). These domains can be used for guiding light using the PBG effect with effective indices below the RI of core material [1, 8]. Identifying the PBG regions requires numerical methods, typically relying on using Floquet-Bloch modes [20] and in the case of cylindrical strands calculating a large number of modes at the high symmetric points of the 1st Brillouin zone [13]. These scans can be extensive and time-consuming particular for high index contrast structures and at large normalized frequencies, potentially leading to long simulation times and convergence issues which can make detailed parameter scans unfeasible within realistic time scales.

One very successful approximate approach to identify the PBG regions was introduced by T. Birks et al. [21] and relies on analyzing only the dispersions of the extremal modes of the cladding SMs (referred in the following as PBG edge modes). All other SMs are not considered, since for low loss guidance it is only of relevance to find the PBG domains. This method requires numerical solving two independent transcendental equations which give rise to the n_eff dispersions of the edge modes and is favorable to full numerical simulations due to substantially reduced calculation effort and memory consumption. This original model (which will be referred here as Birks model) relies on a scalar field approach and has been recently extended to include the full vectorial nature of the modes [22], which, e.g., is important if cladding structures such as arrays of metallic nanowires are analyzed [13].

One drawback of the Birks model is the prerequisite to numerically solve transcendental equations, which imposes using root search methods. Moreover, information about the DOS cannot be calculated for one desired combination of wavelength and n_eff; it is rather required to find all the corresponding edge modes which demands searching starting values, which are typically obtained above the cladding RI. Using these values the modal dispersions of the edge modes have to be tracked accurately into the domains below the cladding RI. Due to the complex behavior of the edge modes in this effective index domain, fine step sizes in the modal tracking have to be used not to lose the dispersions, overall increasing simulation time. As a consequence the strategy of tracking the edge mode dispersions fundamentally does not allow for parameter sweeps of, e.g., pitch and/or strand diameter, and thus it is hard if not impossible to optimize an APBG structure for a given wavelength or finding PBGs at low effective index values.

Here we introduce a new approach for calculating the dispersions of the PBG edge modes relying on the Birks model but using a root-search free ansatz. The idea relies on plotting one binary function as function of any relevant parameter, allowing for the direct identification of the PBG domains, which is a fundamentally different approach than solving the transcendental dispersion equations. It is possible to directly target the effective index values of interest using the binary function, thus allowing for extensive parameter sweeps at a fixed value of n_eff, which is important for finding PBGs at effective mode indices far below the cladding index – a domain which is again hard to access with root-search approaches due to the complex modal behavior in this domain. Our method does not require any kind of discretization, which is especially important when addressing systems with large refractive index contrast at high normalized frequencies.

2. Method

The main idea behind the Birks model is to determine the dispersions of the PBG edge modes, which can be used to identify the PBG regions, i.e. the domains being relevant for guiding light in a defect core using the PBG effect. The general motivation of this very successful ansatz can be found within molecular physics, precisely in the formation of diatomic molecules: Here, the most-bonding state, i.e. the state with the largest electron probability density in the center between the atoms, is created when the (total) electronic wave function Ψ has a maximum exactly in the center of the molecule, imposing the von-Neumann boundary condition (BC) dΨ/ds = 0 [23], with in the molecular science case s being the coordinate along the connection line between the centers of the two atoms. In contrast, the most-antibonding state is characterized by a node Ψ = 0 at half center-to-center distance between the atoms (Dirichlet BC) [23]. Any other molecular state formed by the two atoms lies within these two extremal cases, and the domain of zDOS is outside this regime.

This straightforward concept has been adapted by Birks et al. and applied in an electrodynamical manner to the PBG edge modes of an infinite array of high RI dielectric cylinders embedded in a low RI matrix (Fig. 1(b)) [21]. Assuming an interaction of only the nearest neighboring strands, guided cylinder modes of the same LP mode order can couple together via their fields in the low RI cladding. Similar to the case of diatomic molecules, the bonding state is characterized by a maximum of the electric field E exactly at half the distance between the centers of the cylinders, giving the above mentioned von-Neumann BC (dE/ds = 0 at s = Λ/2 with s being the spatial coordinate along the connection line between the centers of the two cylindrical strands). This corresponds to the state with the highest possible value of propagation constant, i.e. effective mode index, and is referred in the following as top band edge mode. Conversely, the Dirichlet BC, E = 0 at s = Λ/2, imposes a mode with the lowest possible propagation constant, and is referred as bottom band edge mode. The effective index of any other mode resulting from a coupling between neighboring cylinders lies in-between these two extremal cases, and thus the PBG domains, i.e. the domains of zDOS, can be found by analyzing the band edge modes only.

The original Birks model as introduced by Birks el al. includes two major assumptions: First, the hexagonal unit cell of the hexagonal lattice is approximated by a single circular unit cell, which is a valid assumption for small ratios of strand diameter and pitch, d/Λ. This approximation greatly simplifies the calculations, as the dimension of the initial two-dimensional unit cell is reduced to a system which has no azimuthal dependence, i.e. it is sufficient to calculate the edge modes only along one direction, which results in a one-dimensional mathematical problem. The outer radius of the circular unit cell b is found by imposing that the area fraction of high RI material is identical to that of hexagonal geometry, leading to b = Λ(√3/(2π))^1/2.

The second assumption is that the modal fields are scalar (which is a valid assumption if the longitudinal electromagnetic field components can be neglected), imposing that the Birks model is applicable for moderate RI contrasts and low RI in general [24], overall reducing the number of required BCs from six to three for both top and bottom dispersion functions. Using two BCs at the core/cladding interface (continuity of tangential field components) and either a Dirichlet or a von-Neumann BC at s = b, the dispersion functions for the top and bottom PBG modes, f_T and f_B, respectively, are obtained,

The roots of Eqs. (1) give in fact the solutions to the problem discussed here, i.e. the effective indices of top and bottom edge modes for a given configuration (combination of wavelength, strand and cladding index, radius and pitch), whereas the actual value of the dispersion functions are irrelevant. The particular forms of Eqs. (1) and Eqs. (2) use the Bessel functions J_l and Y_l, which result in real-valued dispersion functions if effective indices below the cladding index are of interest. To obtain real-valued dispersion functions for n_eff > n_host, J_l and Y_l have to be replaced by their modified counterparts (see Appendix).

Both dispersion functions cross the zero line, i.e. change their sign, at defined effective indices (example calculation illustrates this in Fig. 2(a)), whereas the crossing point of top edge band has a higher effective mode index than that of the bottom edge mode as expected. This particular behavior is in fact the key feature of our presented calculation approach: Within the PBG domain (bright yellow regions of Fig. 2(a)), both dispersion functions have equal signs, whereas f_T and f_B possess opposite signs within the regions of non-zero DOS. As a consequence, a PBG domain can be identified by the sign of the product of the two dispersion functions. As the actual value of the dispersion function is irrelevant and to indicate the PBG regimes by a zero, we define the normalized dispersion binary function

 

Fig. 2 Dependence of the dispersion functions of top and bottom band edge modes (purple: bottom band edge function f_B, green: top band edge function f_T) on effective mode index. In both diagrams, the blue curves correspond to the binary function which depends on the sign of the product of top and bottom band edge functions. The yellow areas refer to the PBG domains, the gray regions to situation of SM formation in the cladding. The dashed black lines are a guide-to-the-eye to indicate the zero line. (a) is calculated for an example configuration (n_strand = 1.48716, n_host = 1.458, d = 2µm, Λ = 3 µm, λ = 1um, l = 0) showing close-up view of a region with one selected SM domain. (b) is calculated for a second example configuration (n_strand = 1.48716, n_host = 1.458, d/Λ = 0.41, kΛ = 150, l = 0) showing a series of SM guiding region.

Download Full Size | PDF

As a result of this particular definition, the binary function produces a zero in case of zDOS or unity in all other cases, thus allowing for a clear and simple identification of PBG regions (Fig. 2(b)). The idea of our method, i.e., calculating a binary function by Eq. (4), is fundamentally different from a root search ansatz typically used for analyzing PBG fibers, since it does not require tracking the dispersions of the edge mode; it is, however, sufficient to plot this binary function for any desired parameter combinations and check whether f_bin indicates zero values, i.e., PBG domains. A direct identification of PBG regions is also important from the simulation point of view, as the numerical difficulty to access region close to the host RI does not need to be addressed, which solves one issue pointed out in the original work of Birks et al. [21].

3. Results

To verify the validity of our approach we simulate the dispersions of the PBG edge modes using our binary function (Eq. (4)) and using a root-finding algorithm with parameters identical to those of the original work of Birks et al. (Fig. 3) [21]. In detail, we have calculated the effective index dispersion of the top and bottom band edge modes as function of normalized frequency for the following parameters: n_strand = 1.48716, n_host = 1.458, d/Λ = 0.41.

 

Fig. 3 Dispersion of the effective indices of the top and bottom band edge modes as function of normalized frequency. The plots have been obtained using the binary function Eq. (5), light yellow: PBG regions, gray: domains of SM formation in the cladding). To show the validity of our approach the parameters used in the original work of Birks et al. have been used here (n_strand = 1.48716, n_host = 1.458, d/Λ = 0.41) [21]. The dispersion of the edge modes have also been calculated using a standard root finding algorithm (included as points (bright blue circles: top band edge modes, dark yellow circles: bottom band edge modes)). The horizontal dashed black lines separate the regions of total-internal reflection guidance (positive ordinate numbers) from the domain of photonic band gap guidance (negative ordinate numbers). (a) Dispersions of the edge modes of lowest Bessel function order (l = 0). (b) Full dispersion map after overlaying the first seven dispersion maps (l = 0…7).

Download Full Size | PDF

The analysis starts by considering a fixed value of Bessel function order l = 0 (Fig. 3(a)): the effective index of top and bottom edge modes are in close proximity for n_eff values above the cladding index, which is a result of the comparably weak interaction of the modes in the neighboring cylinders. This changes dramatically for effective index values below the cladding index (below the horizontal dashed lines in Fig. 3(a)), which is the relevant regime for PBG guidance, as the effective indices of defect core modes are located within this domain: Here, strong modal interaction between the modes of neighboring strands leads to a substantial modification of the dispersions of the band edge modes compared to the isolated modes or the modes with indices above the cladding index. The regions of non-zero DOS are located in-between the top and bottom edge modes (gray regions in Fig. 3(a)), whereas the DOS is zero outside those regions (bright yellow regions in Fig. 3(a)). It is important to note that regions of non-zero DOS broadens significantly towards smaller values of n_eff, strongly reducing the available parameter space for PBG guidance.

The dispersions of the band edge modes calculated by the root finding ansatz (points in Fig. 3(a)) is entirely reproduced by our binary function with no single deviation from it (colored regions in Fig. 3(a)). Thus our method allows the direct determination of PBG regions for any parameter combination (in Fig. 3(a) kΛ and (β - kn_host)Λ) without calculating the corresponding top and bottom edge modes. We plot the dispersion map using the binary function from −0.4 < (β - kn_host)Λ < 0.15 and 5 < kΛ < 160 without conducting any sort of mode tracking. The underlying mechanism for the perfect agreement between our binary function approach and the root finding ansatz relates to the fact that for a fixed value of effective index (i.e., (β - kn_host)Λ = const., horizontal level in Fig. 3(a)), the dispersion of the top edge mode never crosses the dispersion of the bottom edge mode, i.e., the sequence of f_T and f_B is for a given value of n_eff always top/bottom/top/bottom etc. The physical reason for this is the following: For a fixed Bessel function order (modes labeled by LP_lx), the adjacent top and bottom edge modes can couple and thus anti-cross, which results in the fact that their dispersions are always separated. This is different for modes with different value of l, as such modes are orthogonal and thus the modal dispersions are allowed to cross.

A full PBG map is obtained by overlaying all the dispersion maps which include any kind of edge mode dispersions (Fig. 3(b)), leading to

Due to the different evolutions of the edge mode dispersions particularly at large normalized frequencies, the areas of PBG guidance reduce substantially after overlaying the first seven dispersion maps (light yellow areas in Fig. 3(b)). Again a 100% correspondence between the results obtained using the root-finding algorithm and our binary function approach is obtained (Fig. 3(b)).

One important feature which is evident from Fig. 3(b) is that PBG guidance is difficult, if not impossible, to achieve at small values of effective indices (i.e., n_eff far below the zero line). These regions are of interest for various important applications requiring guidance of light in a water or even a hollow air core. These regions, however, are typically difficult to access from the simulation perspective using root-searching, since the dispersions of the edge modes have to be tracked down to these regimes starting at high effective index values (presumably above the cladding RI), making large parameter sweeps of combinations of pitch, radius and RIs extremely time consuming, as one has to calculate almost the entire PBG map even though only a small effective index region is of interest. Here we believe our approach can provide a straightforward and time-effective solution, since we can directly determine whether the DOS is zero or not by calculating the value of the binary function for a given configuration.

Thus our approach allows designing and optimizing the cladding properties for application in various effective mode index regimes via fixing the value of n_eff. It is well known that the effective index of the guided core mode is only slightly below the index of the core material (e.g., horizontal black dashed line in Fig. 3(b)) in case the defect core diameter is sufficiently large, i.e. the mode having most of its fields inside the defect region. Thus for designing a PBG fiber operating in a predefined wavelength and effective index domain it is thus sufficient in many practical cases to assume an effective index slightly below the RI of the core material.

To demonstrate the design capability of our approach on the just mentioned example (presented in Fig. 3(b)), we determine the regions of zero and non-zero DOS as function of pitch and strand diameter by plotting the binary function using the material RIs associated with Fig. 3. We consider a wavelength of λ = 1 µm and assume an effective index n_eff = 1.456, which is slightly below the cladding RI (n_host = 1.458). This parameter sweep results in a PBG map, in which the domains of vanishing DOS are indicated by zeros (light yellow areas in Fig. 4), whereas the regions of SM formation are represented by finite values (gray areas in Fig. 4).

 

Fig. 4 Representation of photonic band gap maps as function pitch and strand diameter, calculated using the binary function for a fixed effective index (n_eff = 1.4561) and the material parameters defined in Fig. 3 (λ = 1 µm). In both diagrams the light yellow (gray) domains refer to PBG regions (SM formation). The white regions correspond to the ratio of diameter to pitch which is unrealistic from the fiber fabrication perspective. (a) Example of an overlay of the first two lowest order maps (l = 0 and l = 1 (plotted semitransparent)). (b) Complete photonic band gap map, obtained by overlaying all relevant maps 0≤ l ≤10. The inset shows the guided mode in the central defect region with the parameters indicated by the dark yellow dot in (b) (Λ = 7.33 µm, d = 4.7277µm).

Download Full Size | PDF

The results show that various domains of PBG guidance can be identified, which are mainly determined by the two maps associated with the first two lowest Bessel function orders (l = 0 and l = 1, yellow areas in Fig. 4(a)). This observation of the first two Bessel being most relevant for defining the domain of PBG has been observed in several of our calculation and is in correspondence with literature [17]. Overlaying the first ten PBG-maps (maps with varying d and Λ, l_max = 10) then reveals the regions relevant for PBG guidance (Fig. 4(b)), which are located for this example in the vicinity of d/Λ ≈0.54 (example of a propagating defect mode is shown in the inset of Fig. 4(b) (corresponding to the configuration indicated by the dark yellow dot in Fig. 4(b)). Calculating the entire PBG map including all relevant Bessel function order is straightforward and takes with our current self-written non-optimized Python code only a few hours, whereas numerical simulations using the Finite Element method results in a calculation time of several days.

Our calculation method is also particularly important to understand the susceptibility of the individual PBG to structure variation and imperfections: For instance analyzing the PBG domain at around d = 5 µm and Λ = 8 µm shows that this band gap is less sensitive to pitch variation than to variations in strand diameter. Assuming that the fiber cladding is designed with the parameters corresponding to those exactly in the middle of the PBG domain (d = 4.68 µm, Λ = 8.33 µm), a 20% variation of the pitch is still acceptable for PBG guidance (along the vertical direction), whereas the maximum tolerable deviation in diameter is only 10%. Thus our method allows evaluating PBG-fiber design in terms of robustness against structural imperfections and deviations from the ideal geometry.

4. Conclusion

In this contribution we have introduced a new and straightforward calculation method to find the regions of zero density of states for photonic band gap guiding fibers, consisting of hexagonal arrays of high refractive index cylindrical strands in a low index cladding. The method relies on plotting a single binary function resulting from combining the dispersion functions of the edge modes introduced by Birks et al. [21]. Regions of zero density of states are identified by a single binary value for any combination of geometric parameters, materials and wavelength requiring neither a root-finding algorithm nor tracking the dispersion of edge modes. By applying our model to well-established results, we show that our approach perfectly reproduces the results from Birks et al. [21] and is extremely valuable in cases extensive parameter scans are considered, particularly for large strand diameters, which intrinsically have a large number of guided modes. One interesting observation we found is that the photonic band gap domains are mainly determined by the modes with the two lowest order Bessel functions (l = 0 and l = 1). Due to its straightforwardness, our method shows a substantially reduced simulation effort compared to full numerical simulations and allows the evaluation of photonic band gap structures with respect to structural and material inaccuracies.

The idea of finding the domains of zero density of states using a single binary function is not restricted to the Birks model and should also work for more sophisticated dispersion functions of edge modes, as the range of validity of our method stands and falls with the validity range of the model used for deriving the dispersion functions. In a next step we plan to extend the presented scheme to more sophisticated dispersion functions (e.g., those taking into account the vector nature of the fields [13,22]), which include fewer assumptions and thus make the presented approach applicable to more experimentally relevant configurations.