## Abstract

A parameterization of the dispersive conductivity of highly-doped graphene has been developed and is presented for use in finite-difference time-domain simulation of near infrared graphene-based photonic and plasmonic devices. The parameterization is based on fitting a Padé approximant to the conductivity arising from interband electronic transitions. The resulting parameterization provides an accurate spectral representation of the conductivity in the wavelength range 1.3 – 2.3*μ*m which is important for near infrared graphene plasmonics. Finite-difference time-domain simulations of straight graphene plasmonic waveguides of infinite and finite width are presented.

© 2012 Optical Society of America

## 1. Introduction

Graphene is an exciting new material that is receiving significant attention in the research community as a building block for both flexible and inherently novel devices. Graphene is a two-dimensional monoatomic layer of carbon arranged in a honeycomb lattice [1,2]. It has a number of interesting features of interest from both a fundamental physics standpoint as well as an engineering one. The physics of graphene is novel in that it represents one of the only examples of a truly two-dimensional solid that can be completely isolated from a substrate or supporting structure. Due to graphene’s particular crystal structure, the electrons in graphene lose their effective mass and behave like Dirac fermions [3] with velocities approaching 0.01*c*. Graphene is a zero gap semiconductor, and its electronic band structure is conical about its K points.

Already the interesting photonic properties of graphene have been explored for a range of applications. Demonstrations of graphene-based photodetectors capable of operating frequencies as high as 40 GHz at a wavelength of 1.55*μ*m have been carried out, and the authors estimate that the peak bandwidth of a graphene photodetector could be as large as 500 GHz [4]. Subsequent work demonstrated error free detection at 10 Gbps [5]. Fabrizio *et al.* have explored the effect of a graphene layer on a photonic metamaterial which resulted in resonance shifts and transmission enhancements when compared to devices without the graphene [6]. Theoretical proposals have been published for modulators operating at 10*μ*m whose maximum speed is limited by a 0.1 ps carrier relaxation time [7] and graphene-based parallel plate waveguides designed to operate at terahertz frequencies [8]. Blake *et al.* demonstrated liquid crystal switching devices using graphene as a transparent electrode [9]. The nonlinear saturable absorption properties of graphene have been intensely studied recently in demonstrations of passively mode-locked lasers emitting around 1.5 *μ*m with pulse durations shorter than 1 ps [10–13].

Recently Vakil and Engheta discussed metamaterials and transformation optics in which the surface plasmon polariton (SPP) properties of highly doped graphene play a key role [14]. Using graphene for SPP devices instead of noble metals such as gold or silver is advantageous because the material properties of graphene are tunable via external electric or magnetic field, and it supports SPPs with longer propagation lengths. Metamaterials, plasmonics and integrated photonics deal with electromagnetic behavior associated with micrometer and nanometer geometrical features. Typically finding the electromagnetic fields in these subwavelength geometries requires numerical approaches to solve Maxwell’s equations. One of the most popular numerical techniques is the finite difference time domain (FDTD) method [15, 16] which is an explicit discretization of Maxwell’s equations in space and time. FDTD is attractive due to its generality, ease of implementation, linear scaling in execution time with problem size, straight forward parallelizability, favorable speedup with parallelization and ability to handle dispersive and nonlinear materials. It is often the default method when modeling large three-dimensional irregular geometries.

In the frequency domain the relationship between the electric field and current density for isotropic materials is given by *J⃗*(*ω*) = *σ*(*ω*)*E⃗*(*ω*). Because FDTD solves Maxwell’s equations in the time domain, incorporating dispersive materials requires special handling. One of the most efficient approaches uses auxilary differential equations (ADE) [16–18] to model the dispersive material properties in the time domain. The ADE method requires that the dispersive conductivity be represented as a ratio of polynomials in *ω*. Popular forms include the Drude [19], Debye [17], Lorentz [18] and Cole [20] models. While the intraband relaxation contribution to graphene’s conductivity can be represented using the Drude model alone, the inter-band transition component results in a dispersive conductivity function that cannot be directly incorporated into FDTD using the ADE method. In this work a Padé approximant [16,21–23] is used to represent the interband term in a form that can be incorporated into FDTD. The general recipe is described for performing the fit, fitting parameters are given and modeling results of SPPs in graphene are presented.

Recent work has modeled an infinite graphene sheet using FDTD [24]; however, only the intraband (Drude) conductivity was used. While neglecting the interband term may provide good results for long wavelengths, it must be included when considering optical behavior with optical energies near the chemical potential. This work presents a form for the interband conductivity term that is especially important for near infrared SPP modeling in graphene. FDTD is an extremely powerful tool for nanophotonic modeling [25–28], and the results presented here will further enable discovery and design of photonic devices that exploit the exciting properties of graphene.

## 2. Graphene conductivity

Graphene has a frequency-dependent complex-valued conductivity which is shown in Fig. 1 for values of chemical potential ranging from 0.1 to 0.6 eV. The chemical potential is a function of applied static electric field, applied magnetic field and chemical doping [24, 29–31]. Experimental measurements of room temperature, undoped graphene have confirmed that in the frequency regime dominated by interband transitions, the real part of the optical conductivity is approximately constant with a value of *e*^{2}/4*h*̄ [32–34].

For photon energies below two times the chemical potential, the imaginary part of the optical conductivity is negative opening the possibility for transverse magnetic SPP waves [14,31,35]. From the time-harmonic Maxwell equations, the relationship between conductivity and electric permittivity is given by

where*ε*(

_{r}*ω*) is the real part of the electric permittivity and −

*σ*(

*ω*)/

*ω*is the imaginary part [36]. However, if the conductivity

*σ*(

*ω*) =

*σ*(

_{r}*ω*) +

*iσ*(

_{i}*ω*) is complex, then the electric permittivity is more clearly written by performing the complex arithmetic as in

*σ*(

_{i}*ω*) will result in a negative real part of the electric permittivity when

*ε*(

_{r}*ω*) +

*σ*(

_{i}*ω*)/

*ω*< 0. In this case the graphene layer behaves like a thin metal capable of supporting SPP waves.

In the absence of an applied electric or magnetic field, the conductivity is associated with electronic intraband relaxation and interband transitions according to *σ* = *σ _{intra}* +

*σ*. Following the notation of [24, 31], an expression for the conductivity as a function of frequency (

_{inter}*ω*), chemical potential (

*μ*), carrier scattering rate (Γ) and temperature (

_{c}*T*) is given by

*f*(

_{d}*ε*) = 1/(

*e*

^{(ε−μc)/kBT}+ 1) is the Fermi-Dirac distribution function and

*k*is the Boltzmann constant. The first integral describes the intraband carrier relaxation contribution, and the second integral describes the interband carrier transition contribution. The first integral can be evaluated resulting in

_{B}*k*≪ |

_{B}T*μ*|,

_{c}*h*̄

*ω*, the second integral may be approximated by

*μ*= 0.6 eV, so the condition

_{c}*k*≪ |

_{B}T*μ*|,

_{c}*h*̄

*ω*is satisfied even at room temperature

*T*= 300 K. A useful property of graphene is the tunability of its chemical potential with an external DC electric or magnetic field. However, in this work chemically doped graphene is considered where the relationship between the carrier density

*n*and chemical potential is given by

*v*has a value around 10

_{F}^{6}m/s. In [31, 35, 37], values of the chemical potential less than or equal to 0.2 eV were discussed in the context of practical graphene-based electromagnetic devices. These small values of the chemical potential are consistent with low to moderate levels of chemical doping in graphene. Because the imaginary part of the conductivity is negative for optical excitations less than two times the chemical potential, these SPP devices are expected to function only for wavelengths longer than 3

*μ*m. This work is motivated by integration of graphene into standard fiber optic components operating near 1.55

*μ*m. Therefore, in order for SPP waves to exist in these wavelength ranges, the chemical potential of the graphene will need to be at least 0.5 eV. Zhao

*et al.*reported one example of how such a large doping level may be achieved in graphene [38]. Graphene flakes were submerged in solutions of 14 molar concentration of H

_{2}SO

_{4}and the doping concentration was monitored

*in situ*using Raman spectroscopy. They achieved carrier concentrations of 0.2 − 0.4 × 10

^{13}cm

^{−2}in bilayer graphene which is the predicted carrier concentration necessary for near infrared graphene plasmonics [35].

Figure 1(a) displays the real and imaginary parts of the conductivity of graphene normalized to *σ*_{0} for different values of the chemical potential (*μ _{c}*) ranging from 0.1 eV to 0.6 eV. The imaginary part of the conductivity passes from positive to negative when

*h*̄

*ω*< 2

*μ*. As

_{c}*μ*is increased to values exceeding 0.5 eV, the imaginary part of

_{c}*σ*becomes negative for wavelengths longer than 1.5

*μ*m.

## 3. Padé fit to interband conductivity term

A Padé approximant fits a rational polynomial of the form

*jω*for

*ω*preserves the required conjugate symmetry (

*σ*(

*ω*) =

*σ*

^{*}(−

*ω*)) in the Padé fit assuming the

*a*and

_{i}*b*are real. Polynomials of order

_{i}*M*=

*N*= 2 in the numerator and denominator lead to straight forward and efficient FDTD implementation. Use of higher order polynomials did not result in significantly improved fits. The fitting procedure consists of moving the denominator polynomial to the right side and then moving all terms containing

*ω*back to the left as shown below.

*a*and

_{i}*b*. In order to enforce purely real expansion coefficients Eq. (9) must be separated into its real and imaginary parts according to

_{i}*γ*(

*ω*) = Re[

*σ*(

_{inter}*ω*)] and

*η*(

*ω*) = Im[

*σ*(

_{inter}*ω*)], the resulting matrix equation is

*ω*

_{1},

*ω*

_{2}and

*ω*

_{3}are three distinct frequency values at which

*σ*(

_{inter}*ω*) is sampled. Figure 2 displays the real and imaginary parts of the resulting Padé fit to the interband conductivity term associated with

*μ*= 0.6 eV. The fit is very good for wavelengths longer than 1.0

_{c}*μ*m. Table 1 displays the coefficients used in the fit.

Previous work has determined that the dispersion relation for a transverse magnetic (TM) SPP wave on an infinite graphene sheet is given by *β* = *k*_{0}[1 + 2/*η*_{0}*σ*(*ω*)]^{1/2} where *k*_{0} is the free space wave number, and *η*_{0} is the impedance of free space [14, 31]. Figure 3 displays the dispersion relation of a TM SPP wave using the exact formulas (Eqs. (4) and (5)) for *σ* = *σ _{intra}* +

*σ*as well as the Padé approximation to

_{inter}*σ*. For comparison, the dispersion relation using only the intraband term

_{inter}*σ*=

*σ*is shown to illustrate the importance of including the interband term in the near infrared region. The shaded yellow region shows the frequency region where a significantly improved fit is obtained. The yellow shaded region also corresponds to the targeted frequency region for near infrared graphene plasmonics which makes this fitting procedure highly effective for this wavelength region. As expected, for shorter wavelengths above the yellow shaded region, there is less agreement between the frequency dispersion using the Padé fit and the exact formulas.

_{intra}## 4. Compact FDTD

To demonstrate the quality of the Padé fit to the interband conductivity, the compact FDTD (cFDTD) method is used to simulate an infinite sheet. cFDTD is a fully vectorial numerical method for solving Maxwell’s equations [39–45]. For geometries uniform along the *z* direction, the electric and magnetic fields can be written as

*i*th vector component of the electric or magnetic field,

*z*is the propagation direction, and

*β*is the propagation constant. If Eq. (12) is inserted into Maxwell’s curl equations in cartesian coordinates, derivatives with respect to

*z*can be evaluated explicitly and replaced with

*jβ*. The FDTD method can then be used to solve for the two-dimensional field quantities ${f}_{i}^{\beta}(x,y)$. In this case, the problem still incorporates six vector components, but the computational domain requires only the two-dimensional cross section of the waveguide. Examples of the FDTD update equations for

*E*and

_{y}*H*are shown below.

_{x}*J*(

_{y}*ω*) =

*σ*(

*ω*)

*E*(

_{y}*ω*) to the time domain by substituting

*jω*→

*∂/∂t*results in an auxilary differential equation relating

*J*back to

_{y}*E*. The approach allows for the analysis of waveguide structures one user-specified

_{y}*β*-value at a time. The method is labeled “compact” because it takes advantage of geometry uniformity along the propagation direction requiring only a two-dimensional computational grid that represents the cross-section of the graphene geometry. The computational domain is illustrated in Fig. 4(a). Because the graphene sheet is assumed to be infinite in both the

*x*and

*z*directions, derivatives with respect to

*x*can be set to zero. This implies a SPP wave uniform along

*x*and propagating along

*z*. The resulting computational domain is only one-dimensional.

The thin atomic monolayer of graphene is represented as a single grid point along the *y* direction. It should be noted that the conductivity expressions presented in Section 2 correspond to a surface conductivity associated with the atomically thin graphene sheet. The conductivity appearing in the three-dimensional Maxwell’s equations is a volume conductivity with units of Siemens per meter [14, 24]. In order to incorporate an effective volume conductivity, *σ*/Δ*y* is used in the simulation.

The spatial discretization was no larger than Δ*y* = 2.5nm, and a time step no larger than 0.8 times the Courant stability limit is used. The spatial domain included 200 discretization points along *y*. The boundary is truncated with 15 layers of perfectly matched layers (PMLs). TM SPP waves are excited on the graphene sheet using a broadband initial condition. The initial excitation consists of a one time step pulse assigned to one or more field components at low symmetry points in close proximity to the graphene sheet. A time sequence on the order of a few hundred thousand time steps is recorded. The discrete Fourier transform is taken, and Padé interpolation is applied to extract the propagation frequency [23, 46]. This process is repeated for all *β* values of interest. A comparison between the FDTD calculated dispersion relation and the exact result is shown in Fig. 4(c). The agreement is good with a maximum error no larger than 2%. Figure 5 displays the dominant electric field components of the TM SPP wave. Both fields are shown at the same moment in time. The tight field confinement to the graphene layer is clear.

The cFDTD method is capable of calculating the propagation length, and a typical calculated value is 10*λ* or (20*μ*m) which is consistent with [35]. These results generally agree with previous theoretical work which incorporate electron-phonon losses via the relaxation time approximation. However, the relaxation time can be influenced by a number of experimental parameters including doping mechanism and temperature and can vary depending on sample preparation method. However, one of the strengths of the presented fitting procedure is that it can incorporate more realistic experimental measurements of the material properties of graphene in the near infrared region.

## 5. Finite width graphene plasmonic waveguide

To conclude this work, results are presented for a graphene SPP waveguide of finite width. The SPP modes of a finite width structure do not have closed form solutions which makes a numerical calculation necessary. A complete modal analysis of the SPP modes of a finite width graphene waveguide is beyond the scope of the current work. However, an example of the dispersion and electric fields of a TM-like mode is shown in Fig. 6. The FDTD simulation parameters are the same as those in the previous section with the exception of using a two-dimensional grid with 300 grid points along *x* and a Δ*x* no larger than 10 nm. As the width of the graphene is reduced, the dispersion relation is shifted to higher frequencies consistent with well known waveguide theory. The field profiles display the same localized field behavior and symmetry along the *y* direction as the fields for the infinite graphene sheet (Fig. 5). However, the fields along the *x* direction are no longer uniform and display higher order mode behavior.

## 6. Conclusion

This work presents a method for introducing the graphene interband conductivity term into FDTD simulations. The fitting procedure is based on Padé interpolation which leads to straight forward integration with FDTD via the ADE method. The approach is illustrated for highly doped graphene which supports SPP waves in the near infrared. This work is important for modeling and design of near infrared graphene plasmonic devices. FDTD modeling results for relatively simple structures are presented to demonstrate the method; however, much more complicated SPP waveguiding structures may be modeled in future studies such as multilayer structures or graphene layers with curved or non-planar shapes. Furthermore, these results are not restricted to waveguiding geometries; three-dimensional structures such as curved waveguides or graphene-based SPP microcavities can also be investigated with this method.

## Acknowledgments

This research is supported in part by Funds for Advancing Research from the College of Science and Technology at Central Michigan University.

## References and links

**1. **K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, “Electric field effect in atomically thin carbon films,” Science **306**, 666–669 (2004). [CrossRef] [PubMed]

**2. **A. K. Geim, “Graphene: status and propects,” Science **324**, 1530–1534 (2009). [CrossRef] [PubMed]

**3. **A. K. Geim and K. S. Novoselov, “The rise of graphene,” Nat. Mater. **6**, 183–191 (2007). [CrossRef] [PubMed]

**4. **F. Xia, T. Mueller, Y. Lin, A. Valdes-Garcia, and P. Avouris, “Ultrafast graphene photodetector,” Nat. Nanotechnol. **4**, 839–843 (2009). [CrossRef] [PubMed]

**5. **T. Mueller, F. Xia, and P. Avouris, “Graphene photodetectors for high-speed optical communications,” Nat. Photonics **4**, 297–301 (2010). [CrossRef]

**6. **E. D. Fabrizio, A. E. Nikolaenko, and N. I. Zheludev, “Graphene in a photonic metamaterial,” Opt. Express **18**, 8359–8353 (2010).

**7. **D. R. Andersen, “Graphene-based long-wave infrared tm surface plasmon modulator,” J. Opt. Soc. Am. B **27**, 818–823 (2010). [CrossRef]

**8. **G. W. Hanson, “Quasi-transverse electromagnetic modes supported by a graphene parallel-plate waveguide,” J. Appl. Phys. **104**, 084314 (2008). [CrossRef]

**9. **P. Blake, P. D. Brimicombe, R. R. Nair, T. J. Booth, D. Jiang, F. Schedin, L. A. Ponomorenko, S. V. Morozov, H. F. Gleeson, E. W. Hill, A. K. Geim, and K. S. Novoselov, “Graphene-based liquid crystal device,” Nano Lett. **8**, 1704–1708 (2008). [CrossRef] [PubMed]

**10. **H. Zhang, D. Y. Tang, L. M. Zhao, Q. L. Bao, and K. P. Loh, “Large energy mode locking of an erbium-doped fiber laser with atomic layer graphene,” Opt. Express **17**, 17630–17635 (2009). [CrossRef] [PubMed]

**11. **Z. Sun, T. Hasan, F. Torrisi, D. Popa, G. Privitera, F. Wang, F. Bonaccorso, D. M. Basko, and A. C. Ferrari, “Graphene mode-locked ultrafast laser,” ACS Nano **4**, 803–810 (2010). [CrossRef] [PubMed]

**12. **Y.-W. Song, S.-Y. Jang, W.-S. Han, and M.-K. Bae, “Graphene mode-lockers for fiber lasers functioned with evanescent field interaction,” Appl. Phys. Lett. **96**, 051122 (2010). [CrossRef]

**13. **W. D. Tan, C. Y. Su, R. J. Knize, G. Q. Zie, L. J. Li, and D. Y. Tang, “Mode locking of ceramic Nd:yttrium aluminum garnet with graphene as a saturable absorber,” Appl. Phys. Lett. **96**, 031106 (2010). [CrossRef]

**14. **A. Vakil and N. Engheta, “Transformation optics using graphene,” Science **332**, 1291–1294 (2011). [CrossRef] [PubMed]

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

**16. **A. Taflove and S. C. Hagness, *Computational Electrodynamics* (Artech House, Massachusetts, 2000).

**17. **M. Okoniewski, M. Mrozowski, and M. A. Stuchly, “Simple treatment of multi-term dispersion in fdtd,” IEEE Microw. Guid. Wave Lett. **7**, 121–123 (1997). [CrossRef]

**18. **F. Hao and P. Nordlander, “Efficient dielectric function for FDTD simulation of the optical properties of silver and gold nanoparticles,” Chem. Phys. Lett. **446**, 115–118 (2007). [CrossRef]

**19. **A. Vial, A.-S. Grimault, D. Macías, D. Barchiesi, and M. L. de la Chapelle, “Improved analytical fit of gold dispersion: Application to the modeling of extinction spectra with a finite-difference time-domain method,” Phys. Rev. B **71**, 085416 (2005). [CrossRef]

**20. **I. T. Rekanos and T. G. Papadopoulos, “FDTD modeling of wave propagation in Cole–Cole media with multiple relaxation times,” IEEE Antennas Wireless Propag. Lett. **9**, 67–69 (2010). [CrossRef]

**21. **G. A. Baker and P. Graves-Morris, *Padé Approximants* (Cambridge University Press, New York, 1996). [CrossRef] [PubMed]

**22. **S. Dey and R. Mittra, “Efficient computation of resonant frequencies and quality factors of cavities via a combination of the finite-difference time-domain technique and the Padé approximation,” IEEE Microw. Guid. Wave Lett. **8**, 415–417 (1998). [CrossRef]

**23. **A. Mock and J. D. O’Brien, “Direct extraction of large quality factors and resonant frequencies from Padé interpolated resonance spectra,” Opt. Quantum Electron. **40**, 1187–1192 (2008). [CrossRef]

**24. **G. D. Bouzianas, N. V. Kantartzis, C. S. Antonopoulos, and T. D. Tsiboukis, “Optimal modeling of innite graphene sheets via a class of generalized FDTD schemes,” IEEE Trans. Magn. **48**, 379–382 (2012). [CrossRef]

**25. **H. Wang, Y. Wu, B. Lassiter, C. L. Nehl, J. H. Hafner, P. Nordlander, and H. J. Halas, “Symmetry breaking in individual plasmonic nanoparticles,” Proc. Natl. Acad. Sci. USA **103**, 10856–10860 (2006). [CrossRef] [PubMed]

**26. **M. T. Hill, Y.-S. Oei, B. Smalbrugge, Y. Zhu, T. de Vries, P. J. van Veldhoven, F. W. M. van Otten, T. J. Eijkemans, J. P. Turkiewicz, H. de Waardt, E. J. Geluk, S.-H. Kwon, Y.-H. Lee, R. Notzel, and M. K. Smit, “Lasing in metal-coated nanocavities,” Nat. Photonics **1**, 589–594 (2007). [CrossRef]

**27. **I. Ahmed, E. H. Khoo, O. Kurniawan, and E. P. Li, “Modeling and simulation of active plasmonics with the FDTD method by using solid state and Lorentz–Drude dispersive model,” J. Opt. Soc. Am. B **28**, 352–359 (2011). [CrossRef]

**28. **A. Mock, “Modal analysis of nanoplasmonic multilayer spherical resonators,” IEEE Photonics J. **3**, 765–776 (2011). [CrossRef]

**29. **K. Ziegler, “Minimal conductivity of graphene: Nonuniversal values from the kubo formula,” Phys. Rev. B **75**, 233407 (2007). [CrossRef]

**30. **L. A. Falkovsky and S. S. Pershoguba, “Optical far-infrared properties of a graphene monolayer and multilayer,” Phys. Rev. B **76**, 153410 (2007). [CrossRef]

**31. **G. W. Hanson, “Dyadic greens functions and guided surface waves for a surface conductivity model of graphene,” J. Appl. Phys. **103**, 064302 (2008). [CrossRef]

**32. **F. Mak, M. Y. Sfeir, Y. Wu, C. H. Lui, J. A. Misewich, and T. F. Heinz, “Measurement of the optical conductivity of graphene,” Phys. Rev. Lett. **101**, 196405 (2008). [CrossRef] [PubMed]

**33. **R. R. Nair, P. Blake, A. N. Grigorenko, K. S. Novoselov, T. J. Booth, T. Stauber, N. M. R. Peres, and A. K. Geim, “Fine structure constant defines visual transparency of graphene,” Science **320**, 1308 (2008). [CrossRef] [PubMed]

**34. **J. M. Dawlaty, S. Shivaraman, J. Strait, P. George, M. Chandrashenkar, F. Rana, M. G. Spencer, D. Veksler, and Y. Chen, “Measurement of the optical absorption spectra of epitaxial graphene from terahertz to visible,” Appl. Phys. Lett. **93**, 131905 (2008). [CrossRef]

**35. **M. Jablan, H. Buljan, and M. Soljačić, “Plasmonics in graphene at infrared frequencies,” Phys. Rev. B **80**, 245435 (2009). [CrossRef]

**36. **D. K. Cheng, *Field and Wave Electromagnetics* (Addison Wesley, New York, 1992).

**37. **T. Stauber, N. M. Peres, and A. K. Geim, “Optical conductivity of graphene in the visible region of the spectrum,” Phys. Rev. B **78**, 085432 (2008). [CrossRef]

**38. **W. Zhao, P. Tan, J. Zhang, and J. Liu, “Charge transfer and optical phonon mixing in few-layer graphene chemically doped with sulfuric acid,” Phys. Rev. B **82**, 245423 (2010). [CrossRef]

**39. **A. Asi and L. Shafai, “Dispersion analysis of anisotropic inhomogeneous waveguides using compact 2D-FDTD,” Electron. Lett. **28**, 1451–1452 (1992). [CrossRef]

**40. **S. Xiao, R. Vahldieck, and H. Jin, “Full-wave analysis of guided wave structures using a novel 2-D FDTD,” IEEE Microw. Guid. Wave Lett. **2**, 165–167 (1992). [CrossRef]

**41. **S. Xiao and R. Vahldieck, “An efficient 2-D FDTD algorithm using real variables [guided wavestructure analysis],” IEEE Microw. Guid. Wave Lett. **3**, 127–129 (1993). [CrossRef]

**42. **Y. Chen and R. Mittra, “A highly efficient finite-difference time domain algorithm for analyzing axisymmetric waveguides,” Microwave Opt. Technol. Lett. **15**, 201–203 (1997). [CrossRef]

**43. **M. Qiu, “Analysis of guided modes in photonic crystal fibers using the finite-difference time-domain method,” Microwave Opt. Technol. Lett. **30**, 327–330 (2001). [CrossRef]

**44. **A. Mock and J. D. O’Brien, “Dependence of silicon-on-insulator waveguide loss on lower oxide cladding thickness,” in Integrated Photonics and Nanophotonics Research and Applications Topical Meeting (Optical Society of America, Boston, MA, USA, 2008),p. IWG4.

**45. **A. Mock and P. Trader, “Photonic crystal fiber analysis using cylindrical FDTD with Bloch boundary conditions,” PIERS Online **6**, 783–787 (2010). [CrossRef]

**46. **W. Kuang, W. J. Kim, A. Mock, and J. D. O’Brien, “Propagation loss of line-defect photonic crystal slab waveguides,” IEEE J. Sel. Top. Quantum Electron. **12**, 1183–1195 (2006). [CrossRef]