aa r X i v : . [ phy s i c s . op ti c s ] S e p Perturbation Theory for Plasmonic Modulation and Sensing
Aaswath Raman ∗ Department of Applied Physics, Stanford University, Stanford, CA, 94305 USA
Shanhui Fan † Department of Electrical Engineering, Stanford University, Stanford, CA, 94305 USA (Dated: August 21, 2018)We develop a general perturbation theory to treat small parameter changes indispersive plasmonic nanostructures and metamaterials. We specifically apply itto dielectric refractive index, and metallic plasma frequency modulation in metal-dielectric nanostructures. As a numerical demonstration, we verify the theory’s accu-racy against direct calculations, for a system of plasmonic rods in air where the metalis defined by a two-pole fit of silver’s dielectric function. We also discuss new opticalbehavior related to plasma frequency modulation in such systems. Our approachprovides new physical insight for the design of plasmonic devices for biochemicalsensing and optical modulation, and future active metamaterial applications.
I. INTRODUCTION
To design active optical devices such as sensors, switches and modulators, one needsto calculate how a small change in refractive index affects the device response function.For active devices based on dielectric structures described by a frequency-independentdielectric constant distribution ε ( r ), the effect of an index change can be understood interms of a frequency shift ω of the eigenmodes of the system, which is given by first-orderperturbation theory as: ω = − ω R d r ∆ ε ( r ) | E ( r ) | R d r ε ( r ) | E ( r ) | (1)The numerator in Eq. (1) only has contributions from the perturbed regions as describedby ∆ ε ( r ). The shift in the eigenfrequency thus depends on the overlap of the modal electricfield energy with the perturbed region.In recent years, there has been substantial interest in using plasmonics for active devices.While surface-plasmon sensors are already prominent in biochemical sensing applications, remarkable improvements in device performance have been achieved using plasmonic nanos-tructures and metamaterials. Active plasmonic devices have also been implemented formodulation and switching.
Given these developments, it is essential to develop a generalphotonic perturbation theory for plasmonic nanophotonic structures to firmly ground theanalysis and characterization of next-generation devices.Eq. (1), however, is not applicable for plasmonic systems. For example, ε ( r ) can benegative in a metal system, and hence directly applying Eq. (1) could lead to a prediction ofinfinite sensitivity, which is unphysical. Moreover, in plasmonic systems, recent experimentshave successfully varied the plasma frequency of the metal as well, introducing a newdegree of freedom that requires formal theoretical treatment.In this paper, we build upon a formulation recently developed to model the photonicband structure of dispersive material systems, to construct a perturbation theory thatcan predict modal frequency shifts due to changes in the dielectric constants of dispersivesystems involving both metals and dielectrics. We also show that a similar perturbationtheory can be developed when the plasma frequency of the dispersive metal is modulated,and highlight new physical behavior related to such modulation. Our theory thus allows usto treat variations in key parameters of both metal and dielectric components of plasmonicnanostructures.The paper is organized as follows. In Section II we review the photonic band theorypreviously developed to describe dispersive metamaterials and plasmonic nanostructures.We develop a general perturbation theory in Section III to describe the effect of smallparameter changes on the optical modes of dispersive nanostructures. We consider thespecific cases of refractive-index variation in Section IV, and plasma frequency variation inSection V, and demonstrate the accuracy of this perturbation theory numerically for bothcases by comparing to direct calculations. Finally, we conclude in Section VI. II. REVIEW OF DISPERSIVE METAMATERIAL BAND THEORY
As the basis for the developments of this paper, we first briefly review the theory recentlydeveloped to model the photonic modes of dispersive material systems. This theory isapplicable for a nanophotonic structure containing a dispersive material with a dielectricfunction ε ( ω ) = ε ∞ + ε ∞ N X n =1 ω p,n ω ,n − ω + iω Γ n . (2)To describe such a dispersive material, for the n -th pole in the dielectric function, oneintroduces a polarization field P n and a polarization velocity field V n , satisfying ∂ P n ∂t = V n (3) ∂ V n ∂t = ω p,n ε ∞ E − ω ,n P n − Γ n V n (4)These auxiliary fields are then coupled to Maxwell’s equations through: ∂ H ∂t = − µ ∇ × E (5) ∂ E ∂t = 1 ε ∞ ( ∇ × H − N X n =1 V n ) (6)For steady state, with fields varying as exp( iωt ), Eqs. (3)-(6) become iω H = − µ ∇ × E (7) iω E = 1 ε ∞ ( ∇ × H − N X n =1 V n ) (8) iω P n = V n (9) iω V n = ω p,n ε ∞ E − ω ,n P n − Γ n V n (10)and thus define an eigenvalue problem for ω . Eqs. (7)-(10) also define a total energy density W = 14 ( ε ∞ | E | + µ | H | ) + N X n =1 ε ∞ ω p,n ( ω ,n | P n | + | V n | ) (11)The spatial integral of Eq. (11), which represents the total energy of the system, is conservedwhen Γ n = 0 for all poles. Further, defining x = ( H , E , P , V , · · · , P N , V N ) T , whichrepresents a multi-component vector field varying over the whole space, we can write Eqs.(7)-(10) as ω Ax = Bx , (12)where A = diag (cid:0) µ , ε ∞ , · · · , ω ,N /ω p,N ε ∞ , /ω p,N ε ∞ (cid:1) and B = i ∇× · · · − i ∇× · · · i ... . . . ...0 0 · · · − i ω ,N ω p,N ε ∞ − i · · · i ω ,N ω p,N ε ∞ i Γ N ω p,N ε ∞ . (13)Finally, defining y = √ Ax , we can re-write this as an eigenvalue equation for ω : ω y = (cid:16) √ A (cid:17) − B (cid:16) √ A (cid:17) − y . (14)For the lossless case, Γ n = 0, Eq. (14) becomes a Hermitian eigenvalue equation whichresults in an orthogonality condition: Z d r "
14 ( ε ∞ E ∗ a · E b + µ H ∗ a · H b ) + N X n =1 ε ∞ ω p,n ( ω ,n P ∗ a ,n · P b ,n + V ∗ a ,n · V b ,n ) = δ ab (15) III. PERTURBATION THEORY
In this section, we develop a general perturbation theory based on the generalized eigen-value equation for the photonic bands of dispersive nanostructures, Eq. (12). We start fromthe unperturbed system ω A x = B x (16)In the presence of a perturbation, the system matrices become A = A + A , and B = B + B , and as a result, we have( ω + ω )( A + A )( x + x ) = ( B + B )( x + x ) (17)Using Eq. (16) and keeping only first-order terms in Eq. (17), we have ω A x + ω A x + ω A x = B x + B x (18)Eq. (16), in its most general form, describes a lossy system and cannot be written as aHermitian eigenvalue problem. Thus, to calculate ω we also need to determine the lefteigenvector z that satisfies ω z A = z B (19)Multiplying z through Eq. (18) and solving for ω , we find: ω = z B x − ω z A x z A x (20)Eq. (20) is the main result of this paper. In the following sections, we apply this equationto two relevant examples of perturbations in plasmonic systems. IV. APPLICATION I: DIELECTRIC REFRACTIVE-INDEX MODULATION
In this section we consider the specific case of a small change (∆ ε ( r )) in the dielectricconstant of a dielectric region, in a nanostructure consisting of both metal and dielectricregions. The metal region is assumed to be unperturbed. In this case, the perturbationtakes the form A = diag(0 , ∆ ε ( r ) , · · · , , , (21)while B = 0. We now determine the change in modal frequency for the cases when themetal in the metal-dielectric nanostructure is lossless and lossy. A. Lossless case, Γ n = 0 For the lossless case Eq. (12) both A and B are Hermitian. In this case, from Eq. (19),we have z = x † . Thus, Eq. (20) reduces to ω = − ω x † A x x † A x = − ω R d r ∆ ε ( r ) | E ( r ) | R d r W . (22)For the lossless dispersive system, we thus obtain a result that has the same form as Eq. (1)that is now appropriate for a system with dispersion, provided that we consider the total energy density in the system including contributions from the auxiliary mechanical fields.The expression for total energy, Eq. (11), has contributions from multiple Lorentz polesand is a multi-pole extension of the energy density expression previously derived by takingelectric polarization into account explicitly. For the lossless case this reduces to the usualexpression for energy density in metals W = 14 (cid:20) d ( ωεε ∞ ) dω (cid:21) | E | + µ | H | . (23) B. Lossy case, Γ n = 0 For the lossy case, the matrix B in Eq. (16) is no longer Hermitian and z = x † . Thus,Eq. (20) in this case reduces to: ω = − ω z A x z A x (24)While no explicit expression analogous to Eq. (1) can be written for the lossy case, Eq. (24)still allows one to calculate the frequency shift due to a dielectric refractive index change inthe presence of a lossy metal; an important ability in realistic plasmonic sensing schemes.Eq. (24) represents the technically correct way to do perturbation theory, where one needsto determine both the left and right eigenvectors of the general eigenvalue problem. Moreoverthe denominator in Eq. (24) cannot be interpreted as an energy integral. Empirically, onthe other hand, we will show numerically that in fact z ≈ x † , even for metals with realisticloss parameters, and thus the denominator of Eq. (24), z A x ≈ Z d r W (25)where W is the energy density of the mode for the lossy system as defined in Eq. (11),which includes contributions from the mechanical auxiliary fields. We note that in a lossysystem, when multiple poles are involved, there is no simple relation such as Eq. (23) thatcan be used to describe the total energy. Instead, the definition of Eq. (11), which explicitlytakes into account contributions from the auxiliary mechanical fields, must be used. C. Numerical Example
We now numerically verify the accuracy of the perturbation theory results presented.Motivated by a recent experiment, we consider as our model system a two-dimensionalperiodic array of square plasmonic rods in air that are uniform along the third z direction,and numerically solve its eigenmodes. The system has periodicity a = 130nm and the rodhas a side length of s = 0 . a = 58 .
1. Metal Dielectric Function Fit
The plasmonic metal’s dielectric function is defined by fitting silver’s tabulatedpermittivity with a Lorentz-Drude model, consisting of a Lorentz pole and a Drude
0 0.5 1 -10-505 0510
Frequency ( ωa /2 πc ) R e a l P e r m i tt i v i t y I m a g i n a r y P e r m i tt i v i t y ω ω p,1 ω p,2 R e a l P e r m i tt i v i t y
300 500 700 900-40-200
Lorentz-DrudeDrudeJohnson & Christy ε real > 0 (a)(b) Wavelength (nm)
FIG. 1: (a) Comparing tabulated data for the real part of silver’s permittivity against the Drudemodel and Lorentz-Drude fit used in this paper, at optical wavelengths. (b) The imaginaryand real parts of the Lorentz-Drude fit of silver’s permittivity over normalized frequencies for a = 130nm. pole [Eq. (2)]. In Fig. 1(a) we see that using the two poles is more accurate than usingthe Drude model alone for modeling silver’s dielectric function at optical frequencies. TheLorentz pole in this fit is defined by ω , = 0 . ω p, = 0 . = 0 . ω , = 0, ω p, = 0 . = 0 . πc/a . A plot of the real and imaginary parts of this fit over the πa k x F r e q u e n c y ( ω a / π c ) πa H y H x E y E x (a) ω ω c (b) Figs. 3 & 5
Wavevector k x Fig. 6
FIG. 2: Computed band structures for a square lattice of 2D plasmonic rods ( s = 0 . a ) in airfor the (a) TM and (b) TE polarization between ( k x = 0 , k y = 0) and ( k x = π/a, k y = 0) points.The metal is describe by a Lorentz-Drude fit of silver with the presence of the Lorentz pole at ω , highlighted. The cutoff frequency for TM modes is identified at ω c , and frequency regions whereRe[ ε ] > relevant normalized frequency range is presented in Fig. 1(b). In the formalism of SectionII, the Drude pole has zero resonant frequency, and hence only requires the V -field as itsauxiliary field. Thus, in our system, we describe the effects of dispersion in terms of threeauxiliary fields: P and V for the Lorentz pole, and V for the Drude pole.
2. Theory Verification
We analyze the TM and TE modes (with their electric and magnetic fields respectivelypolarized along the z -dimension) of this system in Fig. 2. These band structures arecalculated by implementing a finite-difference spatial discretization of the fields with a Yee πa Wavevector k x πa Lossless case (a) F r e q u e n c y ( ω a / π c ) UnperturbedPerturbedDirect πa Wavevector k x Lossy case (b) πaπa πa FIG. 3: Comparing the perturbation theory prediction (’Perturbed’) and direct solution (’Direct’)of ω in the region identified in Fig. 2(b), for dielectric ∆ ε ( r ) = 0 .
02. The (a) lossless metal(Γ n = 0) and (b) lossy metal cases shows excellent agreement for both the flat surface mode andnon-surface mode. grid, and then solving the eigenvalue equation using the Arnoldi method. In our finitedifference implementation, we have truncated the finite-difference grid appropriately at themetal-air interface to ensure that boundary conditions for the tangential field componentsare satisfied at these interfaces.
For the TM case, we find that a stop-band exists belowa cutoff frequency ω c = 0 . ω ≈ . ω , in theregion where Re[ ε ] >
0, highlighted yellow in Fig. 1(b) and Fig. 2.To verify the accuracy of our perturbation theory approach we highlight a region of the0
Unperturbed ∆ ε = 0 Perturbed ∆ ε = 1e-3 PMLPML as k Δ ω = -9.9e-05 R e fl e c t i v i t y Frequency ( ωa /2 πc ) FIG. 4: Reflectivity spectrum of a finite 2D square lattice of square plasmonic rods ( s = 0 . a ) inair, consisting of 50 layers. The results are obtained with a full-field FDFD simulation (shown inthe inset schematic) and show a reflectivity dip corresponding to the lowest-frequency propagatingmode in the system (the mode at ω c in Fig. 2(a)). Altering the air region by ∆ ε = 1e − ω = − . − ω = − . −
5, well.
TE band structure featuring both a surface and non-surface mode in Fig. 3. We alter thedielectric constant of the dielectric region by ∆ ε = 0 .
02, and calculate the resulting shift ineigenfrequency, using perturbation theory for both lossy and lossless cases. For the losslesssystem, we set Γ n = 0, and use Eq. (22). For the lossy system, we can use the exactperturbation theory result of Eq. (24), as well as the approximation in terms of energydensity in Eq. (25). The results from these two forms of perturbation theory are nearlyidentical to each other. These results, from both forms of the perturbation theory, are then1compared to the band structure obtained by directly solving Eq. (12) for the perturbedsystem. The results from the perturbation theory show excellent agreement with resultsfrom the direct calculation for both the lossless and lossy systems.
3. Practical Application
As an illustration of the practical significance of this result, we analyze a potentialrefractive-index sensing scheme by calculating the reflection/transmission spectrum of 50-layers of the plasmonic rod system considered above, using a full-field 2D finite-differencefrequency-domain (FDFD) simulation. The plasmonic metal, corresponding as before tothe two-pole fit of silver’s dielectric function, is assumed to be lossless for the purposes ofthis illustration. For the TM case, below the cutoff frequency ω c [Fig. 1(a)], no propagatingmodes are supported. Thus, the structure is strongly reflecting. The lowest-frequency dipin the reflection spectrum [Fig. 4] corresponds to the first mode supported by the system,identified previously in Fig. 2(a) at ω c . We can then examine the shift in this reflectivitydip when the dielectric surrounding the rods is changed by ∆ ε (due to, for example, theintroduction of a biochemical analyte). We observe a shift of the dip by ∆ ω = − . − ε = 1e −
3. Using Eq. (22) we calculate the shift of the lowest-frequency mode at k = 0, for the same ∆ ε , to be ω = − . −
5, which matches well with the shift observedin the full-field simulation. Thus, the perturbation theory with respect to the eigenmodes ofthe system can be used to predict shifts in features of the same system’s transmission andreflection spectra.
V. APPLICATION II: METALLIC PLASMA FREQUENCY MODULATION
In this section we demonstrate that we can also treat small changes in the plasma fre-quency of a metal present in a plasmonic nanostructure. This is of interest given a recentexperimental work with ITO where an applied electric field induced a change in carrierconcentration thereby changing its plasma frequency and behavior at optical frequencies. Inother recent work, the infrared plasmonic resonance of a metallic nanostructure (split-ringresonators) immersed in an electrolyte was affected by applying a voltage that altered thestructure’s average electron bulk density. N -th pole experiences a plasma frequency mod-ulation. Then, we note that to first order the effect of ∆ ω p,N affects our matrix elements by ω p,N → ω p,N + 2 ω p,N ∆ ω p,N . With this established we now separately consider the losslessand lossy cases. A. Lossless case, Γ n = 0 For the lossless case the following perturbing matrices, A and B are applicable: A = diag , , · · · , − ω ,N ω p,N ε ∞ (cid:18) ω p,N ω p,N (cid:19) , − ω p,N ε ∞ (cid:18) ω p,N ω p,N (cid:19)! (26) B = · · · · · · i ω ,N ω p,N ε ∞ (cid:16) ω p,N ω p,N (cid:17) · · · − i ω ,N ω p,N ε ∞ (cid:16) ω p,N ω p,N (cid:17) . (27)Since the eigenvalue equation Eq. (12) is Hermitian for the lossless case, Eq. (20) becomes: ω = x † B x − ω x † A x x † A x (28)which reduces to the following expression in terms of the constituent fields (as determinedby solving Eq. (12)): ω = ω R d r W Z d r ω p,N ( r ) ω p,N ω − ω ,N ω p,N ! | P N ( r ) | (29)The frequency shift ω for a given mode is thus directly proportional to the strengthof the modal mechanical field in the N -th pole that is undergoing modulation, and theproximity of the unperturbed modal frequency ω to ω ,N . Of particular interest is that,assuming ∆ ω p,N > ω is positive or negative depending on whether the unperturbedmodal frequency is above or below the Lorentz pole’s resonance frequency ω ,N . We verifythis behavior numerically in Section C below.3 B. Lossy case, Γ n = 0 As noted in the previous section, for the lossy case, Eq. (12) is no longer Hermitian and z = x † . Thus, Eq. (20) in this case reduces to: ω = z B x − ω z A x z A x (30)We note here that B is slightly altered in the lossy case: B = · · · · · · i ω ,N ω p,N ε ∞ (cid:16) ω p,N ω p,N (cid:17) · · · − i ω ,N ω p,N ε ∞ (cid:16) ω p,N ω p,N (cid:17) − i Γ N ω p,N ε ∞ (cid:16) ω p,N ω p,N (cid:17) . (31)As in Sec. IV.B, we note that empirically at optical frequencies z ≈ x † , and the denom-inator of Eq. (30) ∼ R d r W . C. Numerical Example
We consider the same example nanostructure of plasmonic rods in air used in the previoussection, this time to numerically verify the accuracy of the perturbation theory for plasmafrequency modulation in the plasmonic rod. Specifically, we focus again on the region ofthe TE band structure highlighted in Fig. 2(b) featuring both a surface and non-surfacemode. As before, in addition to the original structure’s TE bands, we directly calculatethe bands of the perturbed structure, where the metal rod’s Lorentz pole is modulated as∆ ω p, = 0 . ω p, . We now use the perturbation theory calculations to predict the effect of∆ ω p, and compare it to the direct result. We note excellent agreement in Fig. 5 betweenthe perturbation theory and direct approaches for both lossless and lossy cases.Moreover, we note that in both lossy and lossless cases, the surface plasmon mode expe-riences a substantially greater ω compared to the lower frequency non-surface mode. Thisis consistent with the observation from Raman and Fan that surface modes have strongermechanical field intensities, and the form of Eqs. (29) and (30), which predict greater ω for modes with stronger mechanical fields. This behavior is of practical interest for activeplasmonic device applications.4 Wavevector k x Lossless case (a) F r e q u e n c y ( ω a / π c ) UnperturbedPerturbedDirect
Wavevector k x Lossy case (b) πaπa πa πaπa πa FIG. 5: Comparing the perturbation theory prediction (’Perturbed’) and direct solution (’Direct’)of ω in the same zoomed-in region identified in Fig. 2(b), for ∆ ω p, = 0 . ω p, . The (a) lossless(Γ n = 0) and (b) lossy cases shows excellent agreement for both the flat surface mode and non-surface mode. Note the substantially greater ω for the flat plasmon mode which has a strongermechanical component to its eigenmode, and that ω < Finally, as observed in Eq. (29), for a positive shift of the plasma frequency ∆ ω p, > ω depending on whether a mode’s ω is greater or less than ω , . In Fig. 5 the modes presented are below ω , and ω <
0, asexpected for ∆ ω p, >
0. In Fig. 6, we examine modes in the region above ω , highlightedin Fig. 2(b), and note that ω > πa Wavevector k x πa Lossless case (a) F r e q u e n c y ( ω a / π c ) UnperturbedPerturbedDirect Wavevector k x Lossy case (b) πaπa
FIG. 6: Comparing the perturbation theory prediction (’Perturbed’) and direct solution (’Direct’)of ω in the region above ω , identified in Fig. 2(b), for ∆ ω p, = 0 . ω p, . The (a) lossless (Γ n = 0)and (b) lossy cases show ω > ω < VI. CONCLUSION
We have developed and verified a perturbation theory that provides considerable insightinto the effect of small variations in both the dielectric and metallic components of plasmonicand dispersive metamaterial nanostructures. Using a previously developed photonic bandformalism for general dispersive nanostructures described by an arbitrary number of Lorentzpoles, we are able to accurately predict the effect of modulation in the dielectric refractiveindex and metallic plasma frequency of a metal-dielectric plasmonic system. We specifi-cally highlight new behaviors related to plasma frequency modulation that are of potentialsignificance for future research in active plasmonics.6Given increasing interest in such dispersive systems for sensing and modulation appli-cations, and active devices more generally, the theory developed above will indeed be ofconsiderable utility in future designs and analyses.
Acknowledgments
This work was supported by the Center for Advanced Molecular Photovoltaics (CAMP)(Award No KUSC1-015-21), made by King Abdullah University of Science and Technology(KAUST), and the Interconnect Focus Center, funded under the Focus Center ResearchProgram (FCRP), a Semiconductor Research Corporation entity. ∗ Electronic address: [email protected] † Electronic address: [email protected] E. Chow, A. Grot, L. W. Mirkarimi, M. Sigalas, and G. Girolami, Opt. Lett. , 1093 (2004). N. Mortensen, S. Xiao, and J. Pedersen, Microfluid. Nanofluid. , 117 (2008). I. White and X. Fan, Optics Express , 1020 (2008). J. Robinson, L. Chen, and M. Lipson, Optics Express , 4296 (2008). F. Dell’Olio and V. Passaro, Optics Express , 4977 (2007). J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade,
Photonic Crystals: Moldingthe Flow of Light (Princeton University Press, 2008), 2nd ed. J. Homola, Analytical and Bioanalytical Chemistry , 528 (2003). J. Homola,
Surface plasmon resonance based sensors (Springer Series on Chemical Sensors andBiosensors, 2006). A. Kabashin, P. Evans, and S. Pastkovsky, Nat. Mat. , 867 (2009). J. Rosenberg, R. V. Shenoi, T. E. Vandervelde, S. Krishna, and O. Painter, Appl. Phys. Lett. , 161101 (2010). C. Alleyne, A. Kirk, and R. McPhedran, Optics Express , 8163 (2007). J. A. Dionne, K. Diest, L. A. Sweatlock, and H. A. Atwater, Nano Letters , 897 (2009). W. Cai, J. S. White, and M. L. Brongersma, Nano Letters , 4403 (2009). S. Fan, Nat. Photon. , 76 (2010). R. A. Pala, K. T. Shimizu, N. A. Melosh, and M. L. Brongersma, Nano Letters , 1506 (2008). D. Pacifici, H. J. Lezec, and H. A. Atwater, Nat. Photon. , 402 (2007). L. Shao, M. Ruther, S. Linden, J. Weissmueller, and M. Wegener, in
Photonic Metamaterialsand Plasmonics (Optical Society of America, 2010), p. MMC2. K. Diest, Ph.D. thesis, Caltech, CA, USA (2010). U. Guler and R. Turan, Opt. Express , 17322 (2010). A. Raman and S. Fan, Phys. Rev. Lett. , 087401 (2010). A. Taflove and S. C. Hagness,
Computational Electrodynamics: The Finite-Difference Time-Domain Method, Third Edition (Artech House Publishers, 2005). R. M. Joseph, S. C. Hagness, and A. Taflove, Opt. Lett. , 1412 (1991). R. Loudon, Journal of Physics A: General Physics , 233 (1970). R. Ruppin, Physics Letters A , 309 (2002). P. B. Johnson and R. W. Christy, Phys. Rev. B , 4370 (1972). V. P. Drachev, U. K. Chettiar, A. V. Kildishev, H.-K. Yuan, W. Cai, and V. M. Shalaev, Opt.Express , 1186 (2008). S. Fan, P. R. Villeneuve, and J. D. Joannopoulos, Phys. Rev. B , 11245 (1996). A. Ferrario, M. Yan, and M. Qiu, Submitted (2010). K. C. Huang, P. Bienstman, J. D. Joannopoulos, K. A. Nelson, and S. Fan, Phys. Rev. Lett. , 196402 (2003). G. Veronis, R. W. Dutton, and S. Fan, J. Appl. Phys.97