1 Introduction

In conventional optical fibers, the light confinement is achieved through total internal reflection and photons propagate mainly in the high index center core. A completely different confinement mechanism, Bragg reflection, provides an alternative way of guiding photons. This possibility was first pointed out by Yeh et. al. in Ref. [1], where the concept of a Bragg fiber was introduced.

Since Bragg fibers and conventional optical fibers utilize different guiding mechanisms, it is not surprising that Bragg fibers offer many possibilities that are difficult to achieve in conventional fibers. One such example is the possibility of guiding light in air, which has attracted much recent interest [2]–[6].

A Bragg fiber can also be designed to support a single guided mode without azimuthal dependence (i.e. the TE or TM mode). In contrast with the fundamental mode in conventional fibers, which is always doubly degenerate, these guided Bragg fiber modes are truly single mode. Consequently, many undesirable polarization dependent effects can be completely eliminated in Bragg fibers [6].

One variant of the original Bragg fiber, which has a uniform low-index core, is the so-called coaxial Bragg fiber [6]. Coaxial Bragg fibers also guide light through Bragg reflection, except now that a smaller high-index column is embedded in the center of the low-index core. Fig. 1 gives a schematic comparison of the two fibers. For ease of comparison, in this paper we choose the same cladding structure for both waveguides: layer 1 has an index of refraction n ₁=4.6 and a thickness d ₁=0.333a, whereas layer 2 has an index of refraction n ₂=1.6 and a thickness d ₂=0.667a. Here, a=d ₁+d ₂ is the unit length of periodicity of the multilayered structure. The cores of the two fibers are different: whereas the Bragg fiber has only a uniform core with an index of refraction n_o =1.0 (i.e. an air-core Bragg fiber) and a radius r_o , the coaxial fiber core contains an additional inner core with an index of refraction n_i and a radius r_i . For simplicity we choose to vary only three parameters in this paper: n_i , r_i , and r_o .

For the calculations presented here, we have also filled the coaxial waveguiding region with air. In practice, in order to provide structural support it may be replaced by some dielectric material with a low index of refraction without greatly affecting the main observation of this paper.

Many interesting features can potentially be introduced by the presence of a high-index inner core, but in this paper we are mainly interested in its effect on the dispersion relation of the fiber and whether or not it can aid us in our search for a frequency window of single mode transmission.

2 Theoretical analysis

2.1 Asymptotic calculations

In Ref. [7] we developed a general matrix theory for cylindrically symmetric waveguides with Bragg claddings. In this theory, an arbitrary number of inner dielectric layers are treated exactly and the outside cladding layers are approximated in the asymptotic limit. Approximate mode equations are derived by matching the boundary conditions at the interface and dispersion relations can be found by numerically solving these mode equations. For more information about the asymptotic calculations, the reader is referred to Ref. [7].

Using the asymptotic formalism, we can find confined modes for both the Bragg and the coaxial fibers in the bandgap. With n_i =4.6, r_o =1.4a, and the cladding parameters given in Fig. 1, We plot the dispersion relations in Fig. 2 for different inner core radius r_i .

The main observation here is that as the inner core radius r_i is increased from r_i =0 (air-core Bragg fiber) to r_i =0.267a, the TM band undergoes a large downward shift in frequency whereas the m=1 band moves very little, creating a large single-mode transmission window for the TM mode in the lower half of the bandgap.

Before we begin our investigation of the band shifts in the coaxial fiber, it is instructive to first take a look at the field profiles along the radial direction inside the two fibers. In Fig. 3 we plot the $\overrightarrow{E}$ field profiles for both the TM and the m=1 modes and for later convenience we choose modes with the same propagation wavevector β. Note the drastic reduction in field amplitude for the m=1 mode inside the inner core of the coaxial fiber (Compare Fig. 3(B) with Fig. 3(D)).

2.2 Perturbation Analysis

The disparate behaviours of the two bands under the influence of a high-index inner core can be explained by perturbation theory. Under the perturbation approach, we treat the inner core as a perturbation to an air-core Bragg fiber structure and solve the wave equation for correction to the unperturbed eigenmode frequency.

For a fixed propagation wavevector β, let the unperturbed modes of an air-core Bragg fiber be:

${\overrightarrow{E}}_m={\overrightarrow{\epsilon }}_m$ (r, θ)exp[i(βz-ω_mt)], m=0, 1, 2, …,

which satisfy the unperturbed eigenmode equation:

These modes form a complete orthogonal set and obey the orthonormal relation [8]:

We now consider the effect of a dielectric perturbation Δε(r, θ) induced by the high-index inner core. Let δ${\overrightarrow{E}}_m$ and δω_m be the changes in the mode functions and eigen frequencies, respectively. The eigenmode equation now takes the form:

In a perturbation analysis, we normally neglect the second-order terms such as Δεδ${\overrightarrow{E}}_m$ when we expand the right-hand side of Eq. (3). In this case, however, we realize that not all second-order terms are small. For example, as we noted earlier, the correction in mode function δ${\overrightarrow{E}}_m$ may be large in the inner core region of the coaxial fiber (Compare Fig. 3(B) with Fig. 3(D)). Also, the fractional change in ε(r, θ) may also be significant owing to the large index of refraction of the inner core (For n_i =4.6,Δε=$n_i^2$ -1=20.16≫1.0). So the only quantity here that can be made small is δω_m , provided the size of the inner core is small enough (See Fig. 2).

Assuming that the second-order term (δω_m )² can be neglected, expand Eq. (3) and use Eq. (1), then Eq. (3) simplifies to:

To solve this we expand δ${\overrightarrow{E}}_m$ in terms of the unperturbed mode functions:

where the a_mn ’s are constants. Substituting Eq. (5) for δ${\overrightarrow{E}}_m$ in Eq. (4) and using Eq. (1), we obtain:

If we now scalar-multiply by ${\stackrel{\to *}{E}}_m$ and integrate over all space, we observe that the expression on the left vanishes because of the orthogonal property (2). Thus we obtain the equation:

Since δω_m is a constant, it can be pulled outside the integral:

The above expression can now be used to explain the different effects that a high-index inner core has on the TM and the m=1 bands. First note that Δε is a windowing function that is nonzero only inside the inner core, hence for the mode function (${\overrightarrow{E}}_m$ +δ${\overrightarrow{E}}_m$ ) in the top integral, only the portion that falls inside the inner core will contribute.

For the TM band, we note from Fig. 3(C) that the TM mode has a large E_z component inside the inner core, which would make a significant contribution to the top integral in Eq. (8). This explains why the TM band makes a large downward shift in frequency even with the insertion of a relatively small inner core. In the case of the m=1 band, we refer back to Fig. 3(D), which shows a very small $\overrightarrow{E}$ field inside the inner core. Hence only a small δω is to be expected, verified by Fig. 2 where the m=1 band changes very little for the given range of r_i ’s.

But then why is it that the TM mode can have a large $\overrightarrow{E}$ field inside the inner core and the m=1 mode can’t? The Appendix provides us with a clue. For the TM mode, we see from the Appendix that it has two $\overrightarrow{E}$ field component in the inner core: E_z and E_r , which are proportional to the zeroth order Bessel’s function J ₀(k_cor) and first order Bessel’s function J ₁(k_cor) respectively, where k_co is the transverse wavevector inside the core. Assuming the core is small and knowing that J ₀(0)=1 and J ₁(0)=0, we see that E_z would necessarily be large and E_r would necessarily be small.

The situation is different for the m=1 mode. From the Appendix, we observe that all three $\overrightarrow{E}$ field components, E_z, E_r , and E θ, are present in the inner core. The E_z component can be neglected because it is proportional to J ₁(k_cor), whereas E_r and Eθ each contains a J ₀(k_cor) term and are equal in amplitude at r=0. For a high-index core with n_i =4.6 and an air region outside, boundary matching imposes a large jump discontinuity at the boundary for E_r . In fact, at the boundary $E_r^{\mathit{\text{air}}}$ =$n_i^2$ $E_r^{\mathit{\text{core}}}$ =21.16$E_r^{\mathit{\text{core}}}$ . This explains why E_r is necessarily small in a high-index center core after the mode function has been properly normalized. Eθ is also small in the core because it’s been tied down with E_r at r=0, (i.e., E_θ (0)=E_r (0) in amplitude). As a final remark, we note the same boundary conditions do not have the same “clamping” effect on the TM mode, where the predominant component is E_z , which is continuous across the core boundary.

To verify the validity of our perturbation analysis, we lower the index of refraction of the inner core and repeat the same asymptotic calculations. The resulting band diagrams for n_i =1.45 and r_o =1.4a are plotted in Fig. 4, and the field profiles in Fig. 5. According to our above analysis, the clamping effect imposed by the boundary conditions on the E_r component of the m=1 mode should not be nearly as prominent as before owing to the greatly reduced index of refraction in the core. As a result, the downward shifts of the TM band and the m=1 band should be much smaller than before and of the same order of magnitude. This is exactly what we observe in Fig. 4, although the downward shift of the TM band is still larger because the inner core still provides some clamping even with the reduced index of refraction.

Having found a single-mode transmission window for the TM band in Fig. 2(D), we now explore the possibility of including a zero-dispersion point in this window. If we reduce the outer core radius r_o initially to push both the TM and the m=1 bands higher up in the bandgap for the air-core Bragg fiber, then when we insert the inner core, we know the m=1 band will pretty much stay high up in the bandgap while the TM band will come down all the way to the bottom. In this way we can extend the single-mode frequency window and hopefully it will also include a zero-dispersion point. This is exactly what we did and the resulting band diagrams for n_i =4.6 and r_o =0.867a are plotted in Fig. 6. Note the arrow that points to the zero-dispersion point inside the single-mode window in Fig. 6(D).

2.3 FDTD simulations

Thus far we have used the asymptotic algorithm to study the dispersion properties of the air-core and coaxial Bragg fibers. Now we use 2D finite difference time domain (FDTD) simulations to verify the validity of the asymptotic approach. The 2D FDTD algorithm we used in this paper follows that in Ref. [9]. The key observation is that for any mode propagating in a dielectric structure uniform in the z direction, the z dependence of the fields is simply exp(iβz). For example, we can write the electric field as $\overrightarrow{E}$ (x, y, t)exp(iβz) and the magnetic field as $\overrightarrow{H}$ (x, y, t)exp(iβz). As a result, the 3D Maxwell equations are reduced to 2D ones:

where ∇_⊥ is defined as $\frac{\partial }{\partial x}{\overrightarrow{e}}_x+\frac{\partial }{\partial y}{\overrightarrow{e}}_y$. We can transform these 2D differential equations into a set of finite difference equations following the standard finite difference time domain procedure [10]. At the boundaries of the computational domain, we use the perfectly matched layer boundary condition (PML) [11, 12] to absorb all the outgoing radiation. For the 2D FDTD equations and other details about the implementation of the algorithm, the reader should consult Ref. [9].

We choose to study the coaxial fiber with reduced r_o whose band diagram is given in Fig. 6(D). Specifically the fiber parameters are given as follows: n_i =4.6, r_i =0.267a, n_o =1, r_o =0.867a, n ₁=4.6, d ₁=0.333a, n ₂=1.6 and d ₂=0.667a, where all the parameters are defined in Fig. 1. A total of 5 cladding pairs are used, which should provide good mode confinement based on our experiences.

To verify the dispersion curve in Fig. 6(D) with FDTD simulations, we first give the fiber an initial field and let it evolve in time in accordance with the Maxwell’s equations, with the value of a chosen field component recorded at each time step. We then proceed to perform a fast Fourier transform (FFT) on the time evolution data, and each confined mode would be represented by a peak in the FFT plot in frequency spectrum.

Both the asymptotic and the FDTD results are shown in Fig. 7, where the asymptotic curve is simply copied over from Fig. 6(D). We see that there is a good agreement between the two curves.

In Fig. 8, we show the E_z field distribution of a TM mode, filtered out by FDTD simulation. The frequency and propagation constant of the mode are respectively ω=0.2238(2πc/a) and β=0.2(2π/a). Fig. 8 clearly shows that the guided mode has no θ dependence, just as what we would expect from a TM mode function.

3 Conclusion

In this paper we have studied the dispersion properties of both the air-core and coaxial Bragg fibers. The effect of an inner core on fiber dispersion relation has also been characterized by a perturbation analysis. This analysis leads to the discovery of a large single-mode frequency window in the coaxial fiber for the TM mode, which also contains a zero-dispersion point. The validity of our analysis has been verified by FDTD simulations.

A Appendix

Let’s take the z axis as the direction of propagation, then every field component has the form

(11)$$\psi (r,\theta ,z,t)=\psi (r,\theta )e^{i\left(\beta z-\omega t\right)},$$

where ψ can be E_z, E_r , E_θ, H_z, H_r, H_θ. ω is the angular frequency and β is the propagation constant.

From waveguide theory we know that the transverse field components can be expressed in terms of E_z and H_z . In particular, E_r and E_θ can be written as [1]:

(12)$$E_r=\frac{i\beta }{{\omega }^2\mu \in -{\beta }^2}\left(\frac{\partial }{\partial r}{E}_z+\frac{\omega \mu }{\beta }\frac{\partial }{r\partial \theta }{H}_z\right),$$

(13)$$E_{\theta }=\frac{i\beta }{{\omega }^2\mu \in -{\beta }^2}\left(\frac{\partial }{r\partial \theta }{E}_z-\frac{\omega \mu }{\beta }\frac{\partial }{\partial r}{H}_z\right).$$

Inside the core,

TM mode: E_r is the only transverse $\overrightarrow{E}$ field component

(14)$$E_z=A{J}_0\left(k_{\mathit{co}}r\right),$$

where k_co is the transverse wavevector inside the core,

(15)$$k_{\mathit{co}}=\sqrt{\frac{n^2{\omega }^2}{c^2}-{\beta }^2}.$$

(16)$$E_r=\frac{i\beta }{\frac{n^2{\omega }^2}{c^2}-{\beta }^2}\left(\frac{\partial }{\partial r}{E}_z\right)$$
$$=\frac{i\beta }{\frac{n^2{\omega }^2}{c^2}-{\beta }^2}\frac{d}{\mathit{dr}}\left[A{J}_0\left(k_{\mathit{co}}r\right)\right]$$
$$=\frac{i\beta k_{\mathit{co}}A}{\frac{n^2{\omega }^2}{c^2}-{\beta }^2}{J\prime }_0\left(k_{\mathit{co}}r\right)$$
$$=\frac{-i\beta k_{\mathit{co}}A}{\frac{n^2{\omega }^2}{c^2}-{\beta }^2}{J}_1\left(k_{\mathit{co}}r\right).$$

l=1 mode:

(17)$$E_z=A{J}_1\left(k_{\mathit{co}}r\right)\mathit{cos}\theta ,$$

(18)$$H_z={\mathit{CJ}}_1\left(k_{\mathit{co}}r\right)\mathit{sin}\theta ,$$

(19)$$E_r=\frac{i\beta }{\frac{n^2{\omega }^2}{c^2}-{\beta }^2}\left(\frac{\partial }{\partial r}{E}_z+\frac{\omega \mu }{\beta }\frac{\partial }{r\partial \theta }{H}_z\right)$$
$$=\frac{i\beta }{\frac{n^2{\omega }^2}{c^2}-{\beta }^2}\left\{\frac{\partial }{\partial r}\left[{\mathit{AJ}}_1\left(k_{\mathit{co}}r\right)\mathit{cos}\theta \right]+\frac{\omega \mu }{r\beta }\frac{\partial }{\partial \theta }\left[{\mathit{CJ}}_1\left(k_{\mathit{co}}r\right)\mathit{sin}\theta \right]\right\}$$
$$=\frac{i\beta \mathit{cos}\theta }{\frac{n^2{\omega }^2}{c^2}-{\beta }^2}\{{\mathit{Ak}}_{\mathit{co}}{J\prime }_1\left(k_{\mathit{co}}r\right)+\frac{\mathrm{\omega C}}{r\beta }{\mathit{J}}_1(k_{\mathit{co}}r)\}$$
$$=\frac{i\beta \mathit{cos}\theta }{\frac{n^2{\omega }^2}{c^2}-{\beta }^2}\{{\mathit{Ak}}_{\mathit{co}}{J}_0\left(k_{\mathit{co}}r\right)+(\frac{\mathrm{\omega C}}{r\beta }{-\frac{{\mathit{Ak}}_{\mathit{co}}}{k_{\mathit{co}}r})\mathit{J}}_1\left(k_{\mathit{co}}r\right)\}$$
$$\to \frac{i\beta \mathit{cos}\theta }{\frac{n^2{\omega }^2}{c^2}-{\beta }^2}\left\{A+\frac{\omega C}{\beta }\right\}\frac{k_{\mathit{co}}}{2}\mathit{as}r\to 0,$$

where to go from the third equality to the fourth equality we’ve used the fact that ${J\prime }_1\left(k_{\mathit{co}}r\right)=J_0\left(k_{\mathit{co}}r\right)-\frac{1}{k_{\mathit{co}}r}{J}_1\left(k_{\mathit{co}}r\right)$. And to reach the last equality we have used the fact that as $J_1\left(k_{\mathit{co}}r\right)\to \frac{k_{\mathit{co}}r}{2}$. By the same token,

(20)$$E_{\theta }=\frac{i\beta }{\frac{n^2{\omega }^2}{c^2}-{\beta }^2}\left(\frac{\partial }{r\partial \theta }{E}_z-\frac{\omega \mu }{\beta }\frac{\partial }{\partial r}{H}_z\right)$$
$$=\frac{-i\beta \mathit{sin}\theta }{\frac{n^2{\omega }^2}{c^2}-{\beta }^2}\left\{\frac{A}{r}{\mathit{J}}_1\left(k_{\mathit{co}}r\right)+\frac{\mathit{\omega C}}{\beta }\frac{\partial }{\partial r}\left[{\mathit{J}}_1\left(k_{\mathit{co}}r\right)\right]\right\}$$
$$=\frac{-i\beta \mathit{sin}\theta }{\frac{n^2{\omega }^2}{c^2}-{\beta }^2}\{\frac{A}{r}{\mathit{J}}_1\left(k_{\mathit{co}}r\right)+\frac{\omega k_{\mathit{co}}C}{\beta }{[\mathit{J}}_0(k_{\mathit{co}}r\left)-\frac{1}{k_{\mathit{co}}r}{\mathit{J}}_1\left(k_{\mathit{co}}r\right)]\right\}$$
$$\frac{-i\beta \mathit{sin}\theta }{\frac{n^2{\omega }^2}{c^2}-{\beta }^2}\{\frac{\omega k_{\mathit{co}}C}{\beta }{\mathit{J}}_0\left(k_{\mathit{co}}r\right)\left[\frac{A}{r}-\frac{\omega k_{\mathit{co}}C}{\beta k_{\mathit{co}}r}\right]{\mathit{J}}_1\left(k_{\mathit{co}}r\right)\}$$
$$\to \frac{-i\beta \mathit{sin}\theta }{\frac{n^2{\omega }^2}{c^2}-{\beta }^2}\left\{A+\frac{\omega C}{\beta }\right\}\frac{k_{\mathit{co}}}{2}\mathit{as}r\to 0.$$

 

Fig. 1. Schematics of air-core and coaxial Bragg waveguide cross sections. (A) Air-core Bragg waveguide (B) Coaxial Bragg waveguide. Cladding layer 1 has an index of refraction n ₁=4.6 and a thickness d ₁=0.333a (not shown), whereas layer 2 an index of refraction n ₂=1.6 and a thickness d ₂=0.667a (not shown). Here, a=d ₁+d ₂ is the unit length of periodicity of the cladding. The white region in the core is assumed to be filled with air, with an index of refraction n_o =1. Three parameters are subject to change in this paper: n_i , the index of the refraction of the inner core in (B); r_i , the radius of the inner core in (B); and r_o , the outer core radius in both (A) and (B). Note that an air-core Bragg fiber can be thought as a coaxial Bragg fiber with r_i =0.

Download Full Size | PDF

 

Fig. 2. Band diagrams for n_i =4.6 and r_o =1.4a. r_i is varied as follows: (A) r_i =0 (air-core Bragg fiber), (B) r_i =0.133a, (C) r_i =0.2a, and (D) r_i =0.267a. Note in (D) the TM mode exhibits single-mode behaviour in the lower half of the bandgap, which does not include the zero-dispersion point identified by the arrow.

Download Full Size | PDF

 

Fig. 3. $\overrightarrow{E}$ field profiles for n_i =4.6 and r_o =1.4a. The air-core Bragg fiber has a r_i =0 whereas for the coaxial Bragg fiber r_i =0.133a. The propagation constant β=0.05(2π/a) for all four modes and the eigen frequencies ω are found to be: (A) ω=0.2430(2πc/a), (B) ω=0.2040(2πc/a), (C) ω=0.1943(2πc/a), (D) ω=0.2018(2πc/a). Note in (D) the large jump discontinuity at the core boundary for E_r .

Download Full Size | PDF

 

Fig. 4. Band diagrams for n_i =1.45 and r_o =1.4a. r_i is varied as follows: (A) r_i =0 (air-core Bragg fiber), (B) r_i =0.133a, (C) r_i =0.267a, and (D) r_i =0.4a. Note the downward shifts of the two bands are small compared to those of Fig. 2. Also note that the shifts are of the same order of magnitude for the two bands.

Download Full Size | PDF

 

Fig. 5. $\overrightarrow{E}$ field profiles for n_i =1.45, r_i =0.4a, and r_o =1.4a. The propagation constant β=0.05(2π/a) for both plots and the eigen frequencies ω are found to be: (A) ω=0.2267(2πc/a), (B) ω=0.1966(2πc/a). Note that in (B) the jump discontinuity for E_r at the core boundary is not as sharp as that of Fig. 3

Download Full Size | PDF

 

Fig. 6. Band diagrams for n_i =4.6 and r_o =0.867a. r_i is varied as follows: (A) r_i =0 (air-core Bragg fiber), (B) r_i =0.133a, (C) r_i =0.2a, and (D) r_i =0.267a. Note that in (D) the TM mode exhibits single-mode behaviour almost all the way through the bandgap, with the zero-dispersion point (identified by the arrow) falling inside the single-mode frequency range.

Download Full Size | PDF

 

Fig. 7. Coaxial fiber band diagram for n_i =4.6, r_i =0.267a, and r_o =0.867a. The asymptotic curves are copied from Fig. 6(D) and the FDTD structure is defined with 5 cladding pairs.

Download Full Size | PDF

 

Fig. 8. The E_z field distribution of a TM mode based on FDTD simulation. The parameters of the Bragg fiber are given in the caption of Fig. 1 with n_i =4.6, r_i =0.267a, and r_o =0.867a. The Bragg cladding consists of 5 cladding pairs and the whole fiber is immersed in air. The frequency and propagation constant of the mode are ω=0.2238(2πc/a) and β=0.2(2π/a).

Download Full Size | PDF

References and links

1. P. Yeh, A. Yariv, and E. Marom, “Theory of Bragg fiber,” J. Opt. Soc. Amer. 68, 1196–1201 (1978). [CrossRef]  

2. Y. Fink, D. J. Ripin, S. Fan, C. Chen, J. D. Joannopoulos, and E. L. Thomas, “Guiding optical light in air using an all-dielectric structure,” J. Lightwave Technol. 17, 2039–2041 (1999). [CrossRef]  

3. M. Miyagi, A. Hongo, Y. Aizawa, and S. Kawakami, “Fabrication of germanium-coated nickel hollow waveguides for infrared transmission,” Appl. Phys. Lett. 43, 430–432 (1983). [CrossRef]  

4. N. Croitoru, J. Dror, and I. Gannot, “Characterization of hollow fibers for the transmission of infrared radiation,” Appl. Opt. 29, 1805–1809 (1990). [CrossRef]   [PubMed]  

5. R. F. Cregan, B. J. Mangan, J. C. Knight, T. A. Birks, P. St. J. Russell, P. J. Roberts, and D. C. Allan, “Single-mode photonic band gap guidance of light in air,” Science 285, 1537–1539 (1999). [CrossRef]   [PubMed]  

6. M. Ibanescu, Y. Fink, S. Fan, E. L. Thomas, and J. D. Joannopoulos, “An all-dielectric coaxial waveguide,” Science 289, 415–419 (2000). [CrossRef]   [PubMed]  

7. Y. Xu, G. Ouyang, R. Lee, and A. Yariv, “Asymptotic matrix theory of Bragg fibers,” (submitted to J. Lightwave Technol.).

8. Roy J. Glauber and M. Lewenstein, “Quantum optics of dielectric media,” Phy. Rev. A 43, 467–491 (1991). [CrossRef]  

9. F. Zepparelli, P. Mezzanotte, F. Alimenti, L. Roselli, R. Sorrentino, G. Tartarini, and P. Bassi, “Rigorous analysis of 3D optical and optoelectronic devices by the compact-2D-FDTD method.” Opt. and Quantum Electron. 31, 827–841 (1999). [CrossRef]  

10. K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. AP-14, 302–307 (1966).

11. J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Computat. Phys. 114, 185–200 (1994). [CrossRef]  

12. S. D. Gedney, “An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD lattices,” IEEE Trans. Antennas Propag. 44, 1630–1639 (1996) [CrossRef]