1. Introduction

Waveguide arrays in optical fibers have been implemented mostly as a periodic arrangement of dielectric strands along the fiber axis, resulting in the formation of photonic bands as well as band gaps. By omitting one or more strands in the center of the array a lattice defect (core), surrounded by a periodic microstructured cladding, is introduced. If the effective refractive index of the defect mode coincides with the effective indices of a pass band (i.e. eigenmode of the microstructured cladding), light can tunnel transversely thus leading to a strongly attenuated core mode. However, light can be trapped inside the core at wavelengths at which the cladding exhibits no eigenmodes (band gap) and therefore can propagate along the fiber. Applications of these special kinds of PCFs have emerged in various fields such as (tunable) spectral filtering [1–3], sensing [4], mode-field scaling [5] and fiber lasers [6, 7].

A particular new direction in the field of fiber arrays are metallic longitudinal nanowires [8,9] inside optical fibers (see Fig. 1(a) [10–13], leading to phenomena such as propagating spiraling planar surface plasmon polaritons (SPPs), ultra long-range SPPs [14] or plasmonic molecules [15]. In a periodic lattice of submicron-sized wires (diameter d, pitch Λ) hundreds of SPPs are coupled simultaneously, leading to superplasmonic modes inside the PCF cladding [16, 17]. Three classes of modes can be found in the current unit cell (Fig. 1(b): (i) volume (bulk) modes within the metal wire with extremely high loss (> 100dB/μm); (ii) surface (plasmonic) modes on the wire/silica interface; (iii) volume modes which are mainly located in the silica (grey shaded area in Fig. 1(b). The lowest order volume mode in the dielectric (having the highest effective mode index) is called fundamental space-filling mode (FSM) and separates the volume modes from the region of potential band gap guidance.

 

Fig. 1: (a) Schematic of a metallic nanowire array with central defect (core) in an optical fiber. (b) Cross section of the array and the corresponding real-space unit cell (lattice vectors v⃗₁ and v⃗₂, hole diameter d and pitch Λ)

Download Full Size | PDF

Simulations of the band structure of metallic nanowire arrays usually require sophisticated calculation tools such as finite element method (FEM) mode solvers. The FEM calculations yield all eigenmodes of the photonic structure and its effective refractive indices. These results can be used to calculate the density of states (DOS) in the microstructured cladding. These calculations are quite time-consuming due to various reasons such as fine wavelength discretization (required for an accurate DOS plot), a sufficiently resolved unit cell in reciprocal space and the large number of eigenmodes which need to be calculated. However, the exact knowledge of the DOS value is not required to understand and interpret the light guiding properties for the majority of practical cases. The most relevant information is the spectral position of the edges of the plasmonic bands. Band gap guidance in a lattice defect is only possible in regions where the defect mode is not phase matched to any of the cladding modes. Therefore the key to understanding the light guidance properties is to distinguish between regions of zero and nonzero DOS of the microstructured cladding.

A simple and straightforward approach to calculate the photonic band edges has been introduced by Birks and is based on approximating the hexagonal unit cell of the lattice by a circle (Fig. 2) [18]. This enables the separation of variables, and reduces the problem to finding effective indices of the eigenmodes under consideration of two particular boundary conditions on the edge of the unit cell. This model works quite well and has been applied in modeling all-solid photonic band gap fibers (all-dielectric structures). However, it fails in the case of (plasmonic) surface modes, as it fully neglects the vectorial nature of the electromagnetic fields. This is especially important in the case of cylindrical SPPs, since they show all six electromagnetic field components being non-zero (except the lowest order mode).

 

Fig. 2: Approximation of the hexagonal unit cell (left) by a circular unit cell (right). ε₁ and ε₂ are the permittivities of the metal wire and the cladding, respectively. The red arrows indicate the transverse unit vectors r̂ and θ̂ of the cylindrical coordinate system.

Download Full Size | PDF

Recently, a full-vectorial boundary matching technique was presented by Dong [19]. It is based on a boundary condition of the electric and magnetic fields on a finite number of points along the boundary of the unit cell. This method includes the vectorial properties of the electromagnetic fields and is therefore especially accurate for large ratios d/Λ since it takes into account the correct shape of the unit cell. However, the large number of boundary points which is needed for numerical stability and the inevitable need for numerical mode-solvers in order to find the correct symmetry relations for each set of eigenmodes makes this technique less attractive (especially in terms of demand on simulation time and hardware capacity) for application as practical device design tool.

We report on a straightforward-to-implement and full-vectorial method allowing the approximation of all relevant modal features (plasmonic band edges, the FSM and the modal dispersion of an isolated SPP) of an infinite hexagonal array of metal nanowires without the necessity of using sophisticated mode solvers. The model fully accounts for polarization effects, which are particularly important for metallic nanowire arrays. We achieve a significant reduction of calculation time by a factor of 100 compared to FEM calculations. For all considered cases, our model matches the numeric calculations very well and changing between the different boundary conditions (corresponding to different modes) only requires a simple redefinition of a few entries in a 4 × 4 matrix prior to calculation. The model can also be applied to dielectric structures by merely changing the respective permittivities accordingly.

2. Method

The periodic cladding of a metal-filled PCF, exhibiting a lattice constant (pitch) Λ, is characterized by a hexagonal unit cell and the lattice vectors v⃗₁ and v⃗₂ containing a cylindrical wire (radius a = d/2, permittivity ε₁) embedded in a background material (here silica, permittivity ε₂) (Fig. 2). The circular approximation of the hexagonal unit cell enables a separation of the variables r, θ and z and leads to a mathematical treatment comparable to classical step-index fibers with a finite cladding. This approximation can be interpreted as an effective dimension reduction of the two-dimensional Brillouin zone to a quasi-one-dimensional zone. We chose the radius $R={\left(\sqrt{3}/(2\pi )\right)}^{1/2}\mathrm{\Lambda }$ of the unit cell which corresponds to the equivalent area of the hexagonal unit cell. This choice is considered to be more accurate than the less frequently used definition of the radius R = Λ/2 according to [18, 20, 21]. The structure has then two relevant boundaries: the core-cladding boundary at r = a and the circular unit cell boundary at r = R (right hand side of Fig. 2). The top band edge of a plasmonic band corresponds to a fully bonding (in-phase) linear combination of all surface plasmons in the lattice. The bottom band edge corresponds to a complete antibonding of all surface plasmons, exhibiting a node between adjacent wires [18, 19, 22]. Thus, the one-dimensional wave function is either zero (bottom band edge) or extremal (top band edge) at the unit cell boundary, corresponding to a Dirichlet and Neumann boundary condition, respectively. In the current example of plasmonic nanowire arrays, the electromagnetic fields of the surface plasmons are quasi-transverse magnetic polarized (quasi-TM), having a large fraction of their transverse electromagnetic energy in the azimuthal magnetic and radial electric components. Therefore we anticipate TM-like boundary conditions to be appropriate for this problem, resulting in particular constraints of the radial electric and azimuthal magnetic field components at the radius r = R as shown in Table 1. Using the continuity of the tangential fields at the core-cladding boundary, a 4 × 4 matrix

Table 1:. Boundary conditions corresponding to the top and bottom band edge of a plasmonic band.

View Table | View all tables in this article

Calculating all relevant modes of the lattice by simply inserting the respective terms of Table 2 is the unique feature of this matrix. This can be seen in Eqs. (1) and (2), as the boundary conditions for the unit cell interface at r = R only appear in the parameters η_e and η_h. Therefore a full calculation of all relevant eigenmodes only requires a small change of the matrix itself without the necessity of implementing an entirely different matrix. The individual terms for open boundary conditions (i.e. single metallic wire in an infinite cladding), the top and bottom band edge as well as the FSM are given in Table 2. The FSM requires the boundary conditions E_z(R) = 0 and H_z(R) = 0 with an azimuthal mode order m = 1 [21]. The resulting parameters η_e and η_h (see Table 2) can be calculated straightforwardly from Eqs. (5) and (6).

Table 2:. Definition of the boundary-related parameters η_e and η_h. To address a desired mode it is simply required to use the entries of the respective row of this table in Eqs. (1) and (2).

View Table | View all tables in this article

The advantage of our model is that it provides a preselection of one particular kind of mode (e.g. FSM or one plasmonic mode order), whereas a FEM calculates all modes of the system with almost no option of mode selection. Therefore our model allows investigating the dispersion of the different types of modes independently, which is obviously key to understand the underlying physics of this kind of plasmonic system.

3. Results

We directly compare the results of the presented approximation model with those of a commercial FEM [24] (the details of the exact simulations are given in appendix C). As example system we examine a hexagonal array of gold nanowires embedded in a silica matrix (dielectric functions of silica and gold given in [25] and [26]) and calculate the DOS in the cladding as function of wavelength and effective mode index. The structures under investigation have a moderate ratio d/Λ = 0.20 (d = 800nm, Λ = 4μm) and d/Λ = 0.40 (d = 800nm, Λ = 2μm).

The FEM calculations in Figs. 3(a) and 3(c) show that the regions of non-zero DOS above the silica line correspond solely to plasmonic bands, formed in close proximity to the dispersion of isolated surface plasmons (labeled with their azimuthal mode order). Below the silica line, regions of non-zero DOS dominate (see colorbar in Figs. 3(a) and 3(c), which is a result of the large number of volume modes. However, there also exist areas of vanishing DOS even below the silica line, providing light guidance in the presence of a lattice defect, i.e. an omitted nanowire (core in Fig. 1(a) and are therefore called band gaps. For many applications the zero DOS regions between the silica line and the upper boundary of the volume modes, i.e. the FSM, is most relevant since the effective index of the fundamental core mode is located within this regime (green-shaded regions in Fig. 3(a)–3(d), different band gaps are are labeled (A)–(D)). When the effective index of the defect core mode crosses the plasmonic bands, they phase-match and light is able to transversely tunnel away from the core via the metal nanowire array towards the unstructured cladding region. Optical guidance inside the core is therefore achieved in the regions complementary to the plasmonic bands, the band gap regions.

 

Fig. 3: Comparison of the DOS calculated using a finite element method (left-hand column) and the presented semi-analytic model (right-hand column) for a wire diameter of d = 800nm. White and green regions correspond to zero DOS. The approximation on the right-hand side is composed from the refractive index of silica (blue line), the FSM (black dashed line), the plasmonic bands (yellow shaded regions) as well as the effective indices of the corresponding isolated surface plasmons (black lines). The vertical black dotted line indicates the wavelength (λ = 1000nm) of the fields presented in Fig. 4. (a, b) Λ = 4μm; (c, d) Λ = 2μm.

Download Full Size | PDF

In accordance with the exact simulations, the semi-analytical model (Figs. 3(b) and 3(d), using Eqs. (1) and (2) and the definitions of Table 2) predicts that the top and bottom band edges are located in close vicinity to the effective index of the corresponding isolated surface plasmon dispersion, and that plasmonic bands are formed by a coupling of SPPs of the same order. Below the silica line the isolated surface plasmon modes approach their cutoffs, leading to a strong extension of the modal fields into the silica. Inside the lattice this results in an enhanced interaction of modes with adjacent unit cells, leading to the band formation with increasing bandwidth towards longer wavelengths. For d/Λ = 0.20, the width and depth of every plasmonic band gap (A)–(D) in the region of interest is accurately represented by our model (Figs. 3(a) and 3(b); even for larger ratios such as d/Λ = 0.40, the band gaps (A), (C) and (D) of our model and the FEM spectrally overlap (Figs. 3(c) and 3(d). We observe a difference between our model and the FEM calculations at the long-wavelength side of bandgap (B) which we attribute to two reasons: (i) the relatively small unit cell compared to the actual wavelength, which makes the assumption of a circular unit cell more critical. (ii) higher order modes are generally less confined to the surface of the wire, making them more sensitive to the circular cell approximation. Apart from this, the important parts of the band gap regions (green shading) are very accurately reproduced, with a simulation time reduced by a factor of 100 (from several hours to a few minutes for a full map calculation). The agreement is always better for the bottom band edge, which presumably results from the circular unit cell approximation being less crucial for the corresponding boundary conditions (E_r = 0, H_θ = 0) than for the top band. The spectral positions of the plasmonic bands can be adjusted by the ratio λ/d, whereas the width is predominantly determined by the ratio λ/Λ. The agreement improves if material loss is included. Surface plasmons on an isolated metal wire exhibit a cutoff for m ≤ 2, corresponding to a vanishing imaginary part of their effective index [9, 14]. This behavior is missing in the case of an array of wires. The losses of the SPPs in the array beyond the cutoff wavelength of the corresponding isolated SPPs are in the magnitude of 10dB/mm.

A direct comparison (as example we choose λ = 1000nm, vertical dashed lines in Figs. 3(c) and 3(d)) reveals that the electric field distributions of the FEM simulations are extremely well reproduced by the model for both top and bottom band edge (Fig. 4). The only significant difference is observed for the top band edge of the quadrupole mode, resulting from a symmetry mismatch of plasmon mode and hexagonal unit cell. However, this difference does not have any severe impact on the accuracy of our model. It is noteworthy that for the bottom band edge the radial electric field is not perfectly zero for all positions on the boundary of the hexagonal unit cell. The only significant difference is observed for the top band edge of the quadrupole mode, resulting from a symmetry mismatch of plasmon mode and hexagonal unit cell. However, we believe that this difference does not have a crucial impact on the accuracy of our model. (In order to help reproducing our results, the numerical values of the corresponding effective mode indices can be found in Table 3 in appendix A.)

 

Fig. 4: Comparison of the radial electric field distributions in the real-space unit cells of model and FEM calculations (The top scale bar refers to the amplitude of the electric field). Each row refers to different plasmonic mode order. Left-hand side: bottom band edge. Right-hand side: top band edge.

Download Full Size | PDF

4. Conclusion

We presented a new full-vectorial semi-analytical model to determine plasmonic bands and band gaps in a hexagonal lattice of metal nanowires. It is based on an approximation of the hexagonal unit cell by a circle and introducing particular boundary conditions for the top and bottom band edges. It was shown that the results of our model match those of finite element method simulations, which was additionally confirmed by a direct electromagnetic field comparison. A reduction of simulation time by a factor of 100 was achieved. Moreover, this model allows calculating all relevant modes of the array (e.g. FSM), therefore providing an efficient tool for designing novel fiber devices.

Appendix A: Numerical values

Table 3:. Numerical results of the band edges shown in Fig. 3(d) at a wavelength λ = 1000 nm (vertical black dotted line).

View Table | View all tables in this article

Appendix B: Derivation of the 4 × 4 matrix M̂

In order to determine the eigenmodes of the metal wire in silica it is necessary to solve the full-vectorial Helmholtz equations [27]

(3)$$\mathrm{\Delta }\overrightarrow{E}+k_0^2{\epsilon }_j\overrightarrow{E}=0$$

(4)$$\mathrm{\Delta }\overrightarrow{H}+k_0^2{\epsilon }_j\overrightarrow{H}=0$$

where we assumed a nonmagnetic material and ε_j is the permittivity of the wire (j = 1) and cladding (j = 2). E⃗ and H⃗ are the electric and magnetic field vectors, respectively. Equations (3) and (4) can be rewritten as a modified Bessel differential equation by taking advantage of the separation of the variables and using cylindrical coordinates (separation of longitudinal and transverse field components). A general solution of this equation is a linear combination of the modified Bessel functions I_m(κ_jr) and K_m(κ_jr).

(5)$$E_z^j=[A_j{I}_m({\kappa }_{j}r)+B_j{K}_m({\kappa }_{j}r)]p_m(\theta )$$

(6)$$H_z^j=[C_j{I}_m({\kappa }_{j}r)+D_j{K}_m({\kappa }_{j}r)]q_m(\theta )$$

Here, ${\kappa }_j=+k_0{(n_{\text{eff}}^2-{\epsilon }_j)}^{1/2}$ is the (modified) transverse wave vector and p_m(θ) and q_m(θ) describe the azimuthal dependence of the fields.

(7)$$p_m(\theta )=p_m=\{\begin{array}{cc}\text{cos}(m\theta )& \text{even}\hspace{0.17em}\text{modes}\\ \text{sin}(m\theta )& \text{odd}\hspace{0.17em}\text{modes}\end{array}$$

(8)$$q_m(\theta )=q_m=\{\begin{array}{cc}-\text{sin}(m\theta )& \text{even}\hspace{0.17em}\text{modes}\\ \text{cos}(m\theta )& \text{odd}\hspace{0.17em}\text{modes}\end{array}$$

The parameters A_j, B_j, C_j and D_j are needed to determine the modal fields (i.e. the electromagnetic field components). Depending on the arrangement of the transverse wave vector κ_j, the modified Bessel functions I_m(z) and K_m(z) represent two independent solutions of the Helmholtz Eqs. (3) and (4). In case of real arguments, these functions are, to some extend comparable to exponential functions with a non-oscillating behavior. If their argument is, as in the case discussed here, complex they reveal “mixed” characteristic by showing features of oscillating as well as exponential-type functions. They are directly linked to the “regular” Bessel functions J_m(z) and Y_m(z) [28].

The transverse electric and magnetic field components can be expressed in terms of E_z and H_z:

(9)$$E_r^{(j)}=-i\frac{k_0}{{\kappa }_j^2}\left(n_{\text{eff}}\frac{\partial E_z^{(j)}}{\partial r}+\frac{1}{r}\sqrt{\frac{{\mu }_0}{{\epsilon }_0}}\frac{\partial H_z^{(j)}}{\partial \theta }\right)$$

(10)$$E_{\theta }^{(j)}=-i\frac{k_0}{{\kappa }_j^2}\left(\frac{n_{\text{eff}}}{r}\frac{\partial E_z^{(j)}}{\partial \theta }-\sqrt{\frac{{\mu }_0}{{\epsilon }_0}}\frac{\partial H_z^{(j)}}{\partial r}\right)$$

(11)$$H_r^{(j)}=-i\frac{k_0}{{\kappa }_j^2}\left(n_{\text{eff}}\frac{\partial H_z^{(j)}}{\partial r}-\frac{{\epsilon }_j}{r}\sqrt{\frac{{\epsilon }_0}{{\mu }_0}}\frac{\partial E_z^{(j)}}{\partial \theta }\right)$$

(12)$$H_{\theta }^{(j)}=-i\frac{k_0}{{\kappa }_j^2}\left(\frac{n_{\text{eff}}}{r}\frac{\partial H_z^{(j)}}{\partial \theta }+{\epsilon }_j\sqrt{\frac{{\epsilon }_0}{{\mu }_0}}\frac{\partial E_z^{(j)}}{\partial r}\right).$$

The tangential fields inside the wire (j = 1) and the cladding (j = 2) can then be expressed by the matrix equation

(13)$$\left(\begin{array}{c}{E}_z^j\\ E_{\theta }^j\\ H_z^j\\ H_{\theta }^j\end{array}\right)=\left(\begin{array}{cccc}{a}_j(r)\cdot p_m& a_j^{\prime }(r)\cdot p_m& 0& 0\\ -b_j(r)\cdot q_m& -b_j^{\prime }(r)\cdot q_m& c_j(r)\cdot q_m& c_j^{\prime }(r)\cdot q_m\\ 0& 0& a_j(r)\cdot q_m& a_j^{\prime }(r)\cdot q_m\\ d_j(r)\cdot p_m& d_j^{\prime }(r)\cdot p_m& b_j(r)\cdot p_m& b_j^{\prime }(r)\cdot p_m\end{array}\right)\hspace{0.17em}\hspace{0.17em}\left(\begin{array}{c}{A}_j\\ B_j\\ C_j\\ D_j\end{array}\right),$$

where the matrix elements are given in Table 4. The constant Z₀ = (μ₀/ε₀)^1/2 represents the wave-impedance of free space. Inside the wire the matrix elements a′₁, b′₁, d′₁ and d′₁ need to be zero because the function K_m(κ₁r) exhibits a singularity at r = 0.

Table 4:. Matrix elements of Eq. (13) in the wire (j = 1) and the cladding (j = 2).

View Table | View all tables in this article

In order to set up a 4 × 4 transfer matrix it is necessary to resolve the underdetermination of the system by substituting the matrix elements in Eq. (13) for j = 2 which are multiplied by the parameters A₂ and C₂. This is achieved by introducing special conditions at the unit cell boundary. In the general case of an open boundary (infinite cladding) these parameters are assumed to be zero. This is, however, not correct for the present case where A₂ and C₂ are non-zero. Physically, this emergence of incoming waves is caused by the outgoing waves of all other wires in the periodic lattice. In order to substitute A₂ and C₂ we use the boundary conditions at the approximated circular unit-cell boundary r = R. By applying the boundary conditions E_r(r = R) = 0 and H_θ (r = R) = 0 for the bottom band edge or ∂E_r/∂r|_r=R = 0 and ∂H_θ/∂r|_r=R = 0 for the top band edge, we are able to express the factors A₂ and C₂ in terms of B₂ and D₂.

(14)$$\frac{A_2}{B_2}={\eta }_e=\{\begin{array}{ll}-\frac{K_m^{\prime }({\kappa }_{2}R)}{I_m^{\prime }({\kappa }_{2}R)}\hfill & \text{bottom}\hspace{0.17em}\text{band}\hspace{0.17em}\text{edge}\hfill \\ -\frac{K_m^{″}({\kappa }_{2}R)}{I_m^{″}({\kappa }_{2}R)}\hfill & \text{top}\hspace{0.17em}\text{band}\hspace{0.17em}\text{edge}\hfill \end{array}$$

(15)$$\frac{C_2}{D_2}={\eta }_h=\{\begin{array}{ll}-\frac{K_m({\kappa }_{2}R)}{I_m({\kappa }_{2}R)}\hfill & \text{bottom}\hspace{0.17em}\text{band}\hspace{0.17em}\text{edge}\hfill \\ -\frac{K_m({\kappa }_{2}R)/{\kappa }_{2}R-K_m^{\prime }({\kappa }_{2}R)}{I_m({\kappa }_{2}R)/{\kappa }_{2}R-I_m^{\prime }({\kappa }_{2}R)}\hfill & \text{top}\hspace{0.17em}\text{band}\hspace{0.17em}\text{edge}\hfill \end{array}$$

The primes denote the derivatives of the Bessel functions and can be expressed as

(16)$$I_m^{\prime }(z)=\frac{1}{2}(I_{m-1}(z)+I_{m+1}(z))$$

(17)$$K_m^{\prime }(z)=-\frac{1}{2}(K_{m-1}(z)+K_{m+1}(z))$$

(18)$$I_m^{″}(z)=\frac{1}{4}(I_{m-2}(z)+2I_m(z)+I_{m+2}(z))$$

(19)$$K_m^{″}(z)=\frac{1}{4}(K_{m-2}(z)+2K_m(z)+K_{m+2}(z)),$$

where we have used the recurrence formulas [28]. The functions η_e and η_h describe the weighting of I_m(κ₂r) and K_m(κ₂r) (see Eqs. (5) and (6)) which is needed to fulfill the electromagnetic field conditions at the unit cell boundary.

Inserting equations (14) and (15) into (13) yields the modified matrix for the tangential fields outside the wire and using the continuity of the tangential fields (13) at r = a, it is possible to set up a 4 × 4 matrix equation system

(20)$$\left(\begin{array}{cccc}{a}_1(a)& -[a_2^{\prime }(a)+{\eta }_e{a}_2(a)]& 0& 0\\ -b_1(a)& [b_2^{\prime }(a)+{\eta }_e{b}_2(a)]& c_1(a)& -[c_2^{\prime }(a)+{\eta }_h{c}_2(a)]\\ 0& 0& a_1(a)& -[a_2^{\prime }(a)+{\eta }_h{a}_2(a)]\\ d_1(a)& -[d_2^{\prime }(a)+{\eta }_e{d}_2(a)]& b_1(a)& -[b_2^{\prime }(a)+{\eta }_h{b}_2(a)]\end{array}\right)\hspace{0.17em}\hspace{0.17em}\left(\begin{array}{c}{A}_1\\ B_2\\ C_1\\ D_2\end{array}\right)=0.$$

Substitution of the matrix elements given in Table 4 yields the final matrix (1).

When the effective refractive index of an eigenmode is found, the parameters A₁ to D₁ and A₂ to D₂ have to be determined in order to calculate the electric and magnetic fields. It is crucial to know that only 7 of these factors can be calculated, which then depend on the one remaining. This factor can be used to scale the power of the eigenmode. Using the matrix (1) as well as equations (14) and (15) we are able to derive the amplitudes based on A₁:

(21)$$\left(\begin{array}{c}{A}_1\\ B_1\\ C_1\\ D_1\end{array}\right)=\left(\begin{array}{c}1\\ 0\\ m_0\hspace{0.17em}{m}_{34}{m}_{12}\\ 0\end{array}\right)\cdot A_1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left(\begin{array}{c}{A}_2\\ B_2\\ C_2\\ D_2\end{array}\right)=\left(\begin{array}{c}-{\eta }_e\hspace{0.17em}\hspace{0.17em}{m}_{11}/m_{12}\\ -m_{11}/m_{12}\\ -{\eta }_h\hspace{0.17em}\hspace{0.17em}{m}_0\hspace{0.17em}\hspace{0.17em}{m}_{33}/m_{12}\\ -m_0\hspace{0.17em}\hspace{0.17em}{m}_{33}/m_{12}\end{array}\right)\cdot A_1$$

Here, m_ij are the elements of matrix (1) and m₀ is given by

(22)$$m_0=\frac{m_{11}{m}_{22}-m_{12}{m}_{21}}{m_{23}{m}_{34}-m_{24}{m}_{33}}$$

Appendix C: Finite element simulations

The calculations for the hexagonal unit cell were performed using a commercial FEM [24]. In order to simulate an infinite lattice of metal wires, a hexagonal unit cell containing a circular wire (see Fig. 1(b) was constructed and periodic boundary conditions were applied to the unit cell boundary.

The real-space lattice vectors of the hexagonal array are given by

(23)$${\overrightarrow{v}}_1=\mathrm{\Lambda }\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\overrightarrow{v}}_2=\mathrm{\Lambda }\left(\begin{array}{c}1/2\\ \sqrt{3}/2\\ 0\end{array}\right).$$

The reciprocal lattice vectors w⃗₁ and w⃗₂ can be calculated from (23) yielding

(24)$${\overrightarrow{w}}_1=\frac{2\pi }{\mathrm{\Lambda }}\left(\begin{array}{c}1\\ -1/\sqrt{3}\\ 0\end{array}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\overrightarrow{w}}_2=\frac{2\pi }{\mathrm{\Lambda }}\left(\begin{array}{c}0\\ 2/\sqrt{3}\\ 0\end{array}\right).$$

The first Brillouin zone, which is shown in Fig. 5, can be constructed by using the reciprocal lattice vectors. Due to symmetry reasons, it is only necessary to investigate one twelfth (irreducible wedge) of the first Brillouin zone. The Floquet-Bloch wave vector mesh was generated by superposing the vectors

(25)$${\overrightarrow{u}}_1=\frac{2\pi }{\mathrm{\Lambda }}\left(\begin{array}{c}1/2\\ -1/(2\sqrt{3})\\ 0\end{array}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\overrightarrow{u}}_2=\frac{2\pi }{\mathrm{\Lambda }}\left(\begin{array}{c}1/6\\ 1/(2\sqrt{3})\\ 0\end{array}\right).$$

The eigenmodes of the structure was then calculated for a fixed wavelength at every point on the mesh

(26)$${\overrightarrow{k}}_F=j\cdot ({\overrightarrow{u}}_1+l\cdot {\overrightarrow{u}}_2)$$

where the parameters j and l take values from 0 to 1, yielding a mesh of 1 + N · (N − 1) nodes. We have chosen N = 9 for our simulations which results in a mesh of 73 distinct points in the irreducible Brillouin zone (see right hand side of Fig. 5). Repeating this procedure for every wavelength it was possible to determine the band diagrams for the whole spectrum from 550nm to 1150nm.

 

Fig. 5: Hexagonal unit cell of the reciprocal lattice. Γ, M and K indicate the points of symmetry. The red grid lines illustrate the mesh of Bloch wave vectors in the irreducible wedge of the first Brillouin zone we have used in the FEM simulations. The red numbers correspond to the symmetry-induced weighting factors of the mesh points.

Download Full Size | PDF

In order to visualize the results in one single graph, we calculated the (normalized) density of states [29, 30]

(27)$$\text{DOS}(\lambda ,n_{\text{eff}})=\sum _k{w}_k\sum _i\delta (n-n_{\text{eff}}^{i,k}),$$

where w_k corresponds to a weighting factors for each mesh point and has to satisfy ∑_k w_k = 1. An equal weighting of all mesh points in the irreducible wedge of the Brillouin zone would not describe the symmetry of the hexagonal array correctly [31, 32]. Thus, it is necessary to use a particular relative weighting of each mesh point which is shown as red numbers on the right hand side of Fig. 5. The δ-function in Eq. (27) was approximated by a triangular function with a finite width Δn_eff (FWHM). For the calculations of the d/Λ = 0.20 and d/Λ = 0.40 nanowire array we used Δn_eff = 0.001 and Δn_eff = 0.002, respectively.

Acknowledgments

This work was supported by federal state of Thuringia in the framework of the European Regional Development Fund (ERDF) project SiFAS-P FKZ: B714-07035.

References and links

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

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

3. B. Kuhlmey, B. Eggleton, and D. Wu, “Fluid-filled solid-core photonic bandgap fibers,” J. Lightwave Technol. 27, 1617–1630 (2009). [CrossRef]  

4. W. Y. W. Yuan, G. Town, and O. Bang, “Refractive index sensing in an all-solid twin-core photonic bandgap fiber,” IEEE Sens. J. 10, 1192–1199 (2010). [CrossRef]  

5. S. Saitoh, K. Saitoh, M. Kashiwagi, S. Matsuo, and L. Dong, “Design optimization of large-mode-area all-solid photonic bandgap fibers for high-power laser applications,” J. Lightwave Technol. 32, 440–449 (2013). [CrossRef]  

6. A. Isomäki and O. G. Okhotnikov, “Femtosecond soliton mode-locked laser based on ytterbium-doped photonic bandgap fiber,” Opt. Express 14, 9238–9243 (2006). [CrossRef]   [PubMed]  

7. A. Wang, A. K. George, and J. C. Knight, “Three-level neodymium fiber laser incorporating photonic bandgap fiber,” Opt. Lett. 31, 1388–1390 (2006). [CrossRef]   [PubMed]  

8. C. C. Pfeiffer, E. E. Economou, and K. K. Ngai, “Surface polaritons in a circularly cylindrical interface: surface plasmons,” Phys. Rev. B 10, 3038–3051 (1974). [CrossRef]  

9. L. Novotny and C. Hafner, “Light propagation in a cylindrical waveguide with a complex, metallic, dielectric function,” Phys. Rev. E 50, 4094–4106 (1994). [CrossRef]  

10. H. W. Lee, M. A. Schmidt, R. F. Russell, N. Y. Joly, H. K. Tyagi, P. Uebel, and P. S. J. Russell, “Pressure-assisted melt-filling and optical characterization of Au nano-wires in microstructured fibers,” Opt. Express 19, 12180–12189 (2011). [CrossRef]   [PubMed]  

11. H. W. Lee, M. A. Schmidt, H. K. Tyagi, L. P. Sempere, and P. S. J. Russell, “Polarization-dependent coupling to plasmon modes on submicron gold wire in photonic crystal fiber,” Appl. Phys. Lett. 93, 111102 (2008). [CrossRef]  

12. H. K. Tyagi, H. W. Lee, P. Uebel, M. a. Schmidt, N. Joly, M. Scharrer, and P. S. J. Russell, “Plasmon resonances on gold nanowires directly drawn in a step-index fiber,” Opt. Lett. 35, 2573–2575 (2010). [CrossRef]   [PubMed]  

13. A. Nagasaki, K. Saitoh, and M. Koshiba, “Polarization characteristics of photonic crystal fibers selectively filled with metal wires into cladding air holes,” Opt. Express 19, 3799–3808 (2011). [CrossRef]   [PubMed]  

14. M. A. Schmidt and P. S. J. Russell, “Long-range spiralling surface plasmon modes on metallic nanowires,” Opt. Express 16, 13617–13623 (2008). [CrossRef]   [PubMed]  

15. H. W. Lee, M. A. Schmidt, and P. S. J. Russell, “Excitation of a nanowire molecule in gold-filled photonic crystal fiber,” Opt. Lett. 37, 2946–2948 (2012). [CrossRef]   [PubMed]  

16. M. A. Schmidt, L. Prill Sempere, H. K. Tyagi, C. G. Poulton, P. S. J. Russell, and L. N. P. Sempere, “Waveguiding and plasmon resonances in two-dimensional photonic lattices of gold and silver nanowires,” Phys. Rev. B 77, 033417 (2008). [CrossRef]  

17. C. G. Poulton, M. A. Schmidt, G. J. Pearce, G. Kakarantzas, and P. S. J. Russell, “Numerical study of guided modes in arrays of metallic nanowires,” Opt. Lett. 32, 1647–1649 (2007). [CrossRef]   [PubMed]  

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

19. L. Dong, “A vector boundary matching technique for efficient and accurate determination of photonic bandgaps in photonic bandgap fibers,” Opt. Express 19, 12582–12593 (2011). [CrossRef]   [PubMed]  

20. Y. Li, Y. Yao, M. Hu, L. Chai, and C. Wang, “Improved fully vectorial effective index method for photonic crystal fibers: evaluation and enhancement,” Appl. Optics 47, 399–406 (2008). [CrossRef]  

21. M. Midrio, M. Singh, and C. Someda, “The space filling mode of holey fibers: an analytical vectorial solution,” J. Lightwave Technol. 18, 1031–1037 (2000). [CrossRef]  

22. P. W. Atkins and R. S. Friedman, Molecular quantum mechanics (Oxford University Press, 1997).

23. 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. Lightwave Technol. 19, 2322–2330 (2002).

24. “www.comsol.com,”.

25. I. H. Malitson, “Interspecimen comparison of the refractive index of fused silica,” J. Opt. Soc. Am. 55, 1205–1208 (1965). [CrossRef]  

26. P. G. Etchegoin, E. C. Le Ru, and M. Meyer, “Erratum: An analytic model for the optical properties of gold,” J. Chem. Phys. 127, 189901 (2007). [CrossRef]  

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

28. F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, NIST Handbook of Mathematical Functions (Cambridge University Press, 2010).

29. J. Pottage, D. Bird, T. Hedley, J. Knight, T. Birks, P. Russell, and P. Roberts, “Robust photonic band gaps for hollow core guidance in PCF made from high index glass,” Opt. Express 11, 2854–2861 (2003). [CrossRef]   [PubMed]  

30. J. C. Flanagan, R. Amezcua, F. Poletti, J. R. Hayes, N. G. R. Broderick, and D. J. Richardson, “The effect of periodicity on the defect modes of large mode area microstructured fibers,” Opt. Express 16, 18631–18645 (2008). [CrossRef]  

31. L. Kleinman, “Error in the tetrahedron integration scheme,” Phys. Rev. B 28, 1139–1141 (1983). [CrossRef]  

32. K. Busch and S. John, “Photonic band gap formation in certain self-organizing systems,” Phys. Rev. E 58, 3896–3908 (1998). [CrossRef]