Frequency-domain formulation of photonic crystals using sources and gain

We present a formulation to analyze photonic periodic structures from viewpoints of sources and gain. The approach is based on a generalized eigenvalue problem and mode expansions of sources which sustain optical fields with phase boundary conditions. Using this scheme, we calculate power spectra, dispersion relations, and quality factors of Bloch modes in one-dimensional periodic structures consisting of dielectrics or metals. We also compare the results calculated from this scheme with those from the complex-frequency method. The outcomes of these two approaches generally agree well and only deviate slightly in the regime of low quality factors. © 2013 Optical Society of America OCIS codes: (230.5298) Photonic crystals; (260.2110) Electromagnetic optics; (260.2030) Dispersion; (260.5740) Resonance; (000.4430) Numerical approximation and analysis. References and links 1. L. Brillouin, Wave Propagation in Periodic Structures (Dover, 1953). 2. E. Yablonovitch, “Inhibited spontaneous emission in solid-state physics and electronics,” Phys. Rev. Lett. 58, 2059–2062 (1987). 3. S. John, “Strong localization of photons in certain disordered dielectric superlattices,” Phys. Rev. Lett. 58, 2486– 2489 (1987). 4. J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light (Princeton University Press, 2008), 2nd ed. 5. P. Russell, “Photonic crystal fibers,” Science 299, 358–362 (2003). 6. O. Painter, R. K. Lee, A. Scherer, A. Yariv, J. D. O’Brien, P. D. Dapkus, and I. Kim, “Two-dimensional photonic band-gap defect mode laser,” Science 284, 1819–1821 (1999). 7. S. Noda, M. Yokoyama, M. Imada, A. Chutinan, and M. Mochizuki, “Polarization mode control of twodimensional photonic crystal laser by unit cell structure design,” Science 293, 1123–1125 (2001). 8. H. G. Park, S. H. Kim, S. H. Kwon, Y. G. Ju, J. K. Yang, J. H. Baek, S. B. Kim, and Y. H. Lee, “Electrically driven single-cell photonic crystal laser,” Science 305, 1444–1447 (2004). 9. F. Lemarchand, A. Sentenac, and H. Giovannini, “Increasing the angular tolerance of resonant grating filters with doubly periodic structures,” Opt. Lett. 23, 1149–1151 (1998). 10. H. Raether, Surface Plasmons on Smooth and Rough Surfaces and on Gratings (Springer-Verlag, 1988). 11. M. I. Stockman, “Electromagnetic theory of SERS,” in Surface-Enhanced Raman Scattering, Topics in Applied Physics, K. Kneipp, M. Moskovits, and H. Kneipp, Eds. (Springer-Verlag, 2006). 47–65. 12. S. I. Bozhevolnyi and V. M. Shalaev, “Nanophotonics with surface plasmons–part I,” Photonics Spectra, Jan. Issue, 58–66 (2006). 13. V. M. Shalaev and S. I. Bozhevolnyi, “Nanophotonics with surface plasmons–part II,” Photonics Spectra, Feb. Issue, 66–73 (2006). 14. W. Cai, R. Sainidou, J. Xu, A. Polman, and F. J. G. de Abajo, “Efficient generation of propagating plasmons by electron beams,” Nano Lett. 9, 1176–1181 (2009). 15. A. P. Hibbins, B. R. Evans, and J. R. Sambles, “Experimental verification of designer surface plasmons,” Science 308, 670–672 (2005). #179977 $15.00 USD Received 15 Nov 2012; revised 11 Jan 2013; accepted 11 Jan 2013; published 17 Jan 2013 (C) 2013 OSA 28 January 2013 / Vol. 21, No. 2 / OPTICS EXPRESS 1972 16. C. W. Cheng, M. N. Abbas, M. H. Shih, and Y. C. Chang, “Characterization of the surface plasmon polariton band gap in an Ag/SiO2/Ag T-shaped periodical structure,” Opt. Express 19, 23698–23705 (2011). 17. W. L. Barnes, T. W. Preist, S. C. Kitson, and J. R. Sambles, “Physical origin of photonic energy gaps in the propagation of surface plasmons on gratings,” Phys. Rev. B 54, 6227–6244 (1996). 18. A. Kocabas, S. S. Senlik, and A. Aydinli, “Plasmonic band gap cavities on biharmonic gratings,” Phys. Rev. B 77, 195130 (2008). 19. S. I. Bozhevolnyi, J. Erland, K. Leosson, P. M. W. Skovgaard, and J. M. Hvam, “Waveguiding in surface plasmon polariton band gap structures,” Phys. Rev. Lett. 86, 3008–3011 (2001). 20. F. Wu, D. Han, X. Hu, X. Liu, and J. Zi, “Complete surface plasmon-polariton band gap and gap-governed waveguiding, bending and splitting,” J. Phys.: Condens. Matter 21, 185010 (2009). 21. V. Kuzmiak and A. A. Maradudin, “Photonic band structures of oneand two-dimensional periodic systems with metallic components in the presence of dissipation,” Phys. Rev. B 55, 7427–7444 (1997). 22. S. Guo and S. Albin, “Simple plane wave implementation for photonic crystal calculations,” Opt. Express 11, 167–175 (2003). 23. C. C. Chang, R. L. Chern, C. C. Chang, and R. R. Hwang, “Interfacial operator approach to computing modes of surface plasmon polaritons for periodic structures,” Phys. Rev. B 72, 205112 (2005). 24. R. L. Chern, C. C. Chang, and C. C. Chang, “Analysis of surface plasmon modes and band structures for plasmonic crystals in one and two dimensions,” Phys. Rev. E 73, 036605 (2006). 25. E. Moreno, D. Erni, and C. Hafner, “Band structure computations of metallic photonic crystals with the multiple multipole method,” Phys. Rev. B 65, 155120 (2002). 26. A. G. Hanif, T. Arima, and T. Uno, “Finite-difference frequency-domain algorithm for band-diagram calculation of 2-D photonic crystals composed of Debye-type dispersive materials,” IEEE Antennas Wireless Propag. Lett. 11, 41–44 (2012). 27. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Arctech House, 2005), 3rd ed. 28. P. G. Petropoulos, “Phase error control for FD-TD methods of second and fourth order accuracy,” IEEE Trans. Antennas Propag. 42, 859–862 (1994). 29. Y. Cao, Z. Hou, and Y. Liu, “Finite difference time domain method for band-structure calculations of twodimensional phononic crystals,” Solid State Commun. 132, 539–543 (2004). 30. Y. Tanaka, Y. Tomoyasu, and S. Tamura, “Band structure of acoustic waves in phononic lattices: Two-dimensional composites with large acoustic mismatch,” Phys. Rev. B 62, 7387–7392 (2000). 31. T. Søndergaard and B. Tromborg, “General theory for spontaneous emission in active dielectric microstructures: example of a fiber amplifier,” Phys. Rev. A 64, 033812 (2001). 32. Y. Chen, T. R. Nielsen, N. Gregersen, P. Lodahl, and J. Mørk, “Finite-element modeling of spontaneous emission of a quantum emitter at nanoscale proximity to plasmonic waveguides,” Phys. Rev. B 81, 125431 (2010). 33. V. Siahpoush, T. Søndergaard, and J. Jung, “Green’s function approach to investigate the excitation of surface plasmon polaritons in a nanometer-thin metal film,” Phys. Rev. B 85, 075305 (2012). 34. A. G. Vlasov and O. P. Skliarov, “An electromagnetic boundary value problem for a radiating dielectric cylinder with reflectors at both ends,” Radio. Eng. Electron. Phys. 22, 17–23 (1977). 35. E. I. Smotrova and A. I. Nosich, “Mathematical study of the two-dimensional lasing problem for the whisperinggallery modes in a circular dielectric microcavity,” Opt. Quantum Electron. 36, 213–221 (2004). 36. E. I. Smotrova, A. I. Nosich, T. M. Benson, and P. Sewell, “Cold-cavity thresholds of microdisks with uniform and nonuniform gain: quasi-3-d modeling with accurate 2-d analysis,” IEEE J. Sel. Top. Quantum Electron. 11, 1135–1142 (2005). 37. A. I. Nosich, E. I. Smotrova, S. V. Boriskina, R. M. Benson, and P. Sewell, “Trends in microdisk laser research and linear optical modelling,” Opt. Quantum Electron. 39, 1253–1272 (2007). 38. S. W. Chang, “Full frequency-domain approach to reciprocal microlasers and nanolasers-perspective from Lorentz reciprocity,” Opt. Express 19, 21116–21134 (2011). 39. S. W. Chang, “Confinement factors and modal volumes of microand nanocavities invariant to integration regions,” IEEE J. Sel. Top. Quantum. Electron. 18, 1771–1780 (2012). 40. E. Merzbacher, Quantum Mechanics (Wiley and Sons, New York, 1998), 3rd ed. 41. H. A. Lorentz, “The theorem of poynting concerning the energy in the electromagnetic field and two general propositions concerning the propagation of light,” Verh. K. Akad. Wet. Amsterdam, Afd. Natuurkd. 4, 176–187 (1896). 42. J. R. Carson, “A generalization of reciprocal theorem,” Bell System Technical Journal 3, 393–399 (1924). 43. C. A. Balanis, Advanced Engineering Electromagnetics (Wiley and Sons, 1989). 44. F. Olyslager, “Properties of and generalized full-wave transmission line models for hybrid (bi)(an)isotropic waveguides,” IEEE Trans. Microw. Theory Techn. 44, 2064–2075 (1996). 45. D. Pissoort and F. Olyslager, “Study of eigenmodes in periodic waveguides using the lorentz reciprocity theorem,” IEEE Trans. Microw. Theory Techn. 52, 542–553 (2004). 46. A. D. Yaghjian, “Bidirectionality of reciprocal, lossy or lossless, uniform or periodic waveguides,” IEEE Microw. Wireless Compon. Lett. 17, 480–482 (2007). #179977 $15.00 USD Received 15 Nov 2012; revised 11 Jan 2013; accepted 11 Jan 2013; published 17 Jan 2013 (C) 2013 OSA 28 January 2013 / Vol. 21, No. 2 / OPTICS EXPRESS 1973


Introduction
Photonic periodic structures [1][2][3][4] have been widely utilized in optical fibers [5], laser cavities [6][7][8], and gratings [9] due to the effect of photonic bands.For these structures composed of generic dielectrics at optical frequencies, the material dispersions are often mild.As a result, the time-harmonic Maxwell's equations in a source-free periodic structure can be written in the form of a generalized eigenvalue (GE) problem, in which the frequency ω b of the Bloch mode E(r) is solved: εr (r + R m ) = εr (r); where εr (r) is the relative permittivity tensor of the structure, which is approximately dispersionless in the frequency range of interest; c is the speed of light in vacuum; R m (m = 1 − 3) are primitive vectors of the unit cell; k is the wave vector corresponding to the crystal momentum; and F(r) is the Bloch periodic part of E(r).In Eq. (1a), the factor (ω b /c) 2 plays the role of eigenvalues, and the convention of exp(−iω b t) is adopted.The real part Re[ω b ] indicates the dispersion curve while the ratio −Re[ω b ]/(2Im[ω b ]) is the quality (Q) factor.Equation (1a) is often converted into an approximate matrix form and solved numerically.The condition of minor frequency dispersions is not always valid.This fact becomes especially true if dispersive metals at optical frequencies are incorporated into periodic structures.Since surface plasmon polaritons (SPPs) may be excited in the presence of metals [10][11][12][13], the effect of dispersions on wave propagations has to be taken into account a priori.Under these circumstances, the permittivity tensor, now denoted as ε(r, ω), does depend on the frequency ω, and solving the SPP band diagrams with frequency dispersions is necessary in the first place [14][15][16][17][18][19].The frequency dispersion, nevertheless, turns the GE problem in Eq. (1a) into a nonlinear one.A few authors have proposed specific techniques [20][21][22][23][24][25][26] for the nonlinear GE problem in the presence of dispersions.In addition to complicated numerical analyses, one concern is that ε(r, ω) needs to be extended to the complex frequency (complex ω) because the imaginary part of the solution ω b is present due to the loss.The generalization often requires some preexisting theories such as the Drude (-Lorentz) model and is not always straightforward.
The finite-difference-time-domain (FDTD) method can model metallic components by evolving Maxwell's equations in the time domain [27,28].However, the approach has loworder spatial and temporal accuracies.In addition, to obtain ω b and field distributions of Bloch modes, numerical outcomes corresponding to different locations of the excitation source usually have to be examined because they may critically depend on positions of the source [29,30].
To overcome these shortcomings, we generalize the idea of sources, gain, and eigenfunction expansions [31][32][33][34][35][36][37][38] to fields, sources, and Bloch modes in photonic crystals (PhCs).At a given real frequency ω, the GE problem in Eq. ( 1a) is transformed into a modified one in order to construct a basis set for expansions of the optical field and sustaining sources, in analogy to the formulation of cavities in Ref. [38].The presented formulation is able to take two factors which are not often addressed simultaneously under a single frame into account.First, the resonance as a characteristic of optical modes is manifest.The Bloch mode in this formalism can be regarded as or at least well approximated by the obtained optical field when ω hits the so-called resonance frequency.Second, the fact that sources and gain only exist in specific active regions is embedded in the formulation explicitly.This feature enables us to obtain the threshold gains and confinement factors of Bloch modes with geometries of unit cells and active regions as well as the distribution of absorption automatically taken into account [39].In fact, the presented formalism is more like a resonance calculation based on the variation of complex permittivity rather than the direct perturbation of the resonance frequency.With these points, we apply the formulation to a few one-dimensional (1D) periodic structures including the nondispersive and lossless dielectric/dielectric, nondispersive but lossy metal/dielectric, as well as dispersive and lossy metal/dielectric bilayers.Band diagrams, whitesource power spectra, and Q factors of Bloch modes are calculated based on this approach and compared to those from the dispersive complex-ω method similar to Eq. (1a).The results from these two approaches in general agree with each other and only deviate slightly as Q factors are low.The minor difference can be attributed to the fact that this scheme is constructed from viewpoints of active photonic devices, while the complex-ω method is based on passive ones.
The rest of this paper is organized as follows.The GE formulation, band diagrams, power spectra, and Q factors are described in section 2. The biorthogonality relation of modes in PhCs based on this formulation is also provided in the appendix.Calculations of 1D periodic bilayers using this approach are presented and compared with those of the complex-ω method in section 3. The conclusion is given in section 4.

Formulation
The unit cell Ω uc of a typical PhC and its surface S uc are shown in Fig. 1.It contains an active region Ω a , solely inside which the current source sustaining the field exists.The remaining part of the unit cell is denoted as Ω b .Our goal is to obtain the field E k (r) generated by a given source J s,k (r) only present in Ω a , which additionally satisfies the Bloch phase boundary condition (BC) of the wave vector k in the first Brillouin zone (BZ) at a frequency ω: The invariance of εr (r, ω) under translations of primitive vectors in Eq. (1b) ensures that E k (r) satisfies the same phase BC.One procedure to solve Eq. (2a) is to expand E k (r) and J s,k (r) (at least partially) with some basis sets of fields {f nk (r, ω)} and sources {j s,nk (r, ω)}, respectively, under the constraint of the phase BC (n is the label of the basis).For simplicity, we demand that j s,nk (r, ω) generates f nk (r, ω) so that the two fields also satisfy the wave equation in Eq. (2a).In this way, the field E k (r) and source J s,k (r) can be approximately expanded with the same set of expansion coefficients {c nk (ω)} as A sensible construction of j s,nk (r, ω) is to set it proportional to f nk (r, ω) in Ω a but force it to vanish elsewhere [38], namely, where U(r) is also a periodic function in repeated unit cells and is unity in Ω a but zero in Ω b ; and Δε r,nk (ω) is a proportional factor.The function U(r) assures that J s,k (r) is only present in Ω a through the expansion of the basis set {j s,nk (r, ω)}.The substitution of the ansatz in Eq. (4a) into the wave equation of f nk (r, ω) leads to a GE problem as follows: In Eq. (5a), the factor (ω/c) 2 Δε r,nk (ω) at the right-hand side plays the role of eigenvalues and derives its dependence on ω through solving the GE problem.If the phase BC is imposed on f nk (r, ω), it is applied to j s,nk (r, ω) in Eq. (4a) automatically.This condition ensures that J s,k (r) also satisfies the phase BC.In addition, the phase BC in Eq. (5b) and Bloch theorem enables us to express f nk (r, ω) as where ψ ψ ψ nk (r) is the Bloch periodic part of f nk (r, ω), Similarly, we write the source j s,nk (r, ω) in the Bloch form as where ς ς ς s,nk (r, ω) is the Bloch periodic part of the source j s,nk (r, ω).
We note that the basis set {j s,nk (r, ω)} may not exhaust all the possibilities of the current source J s,k (r).The basis current density in {j s,nk (r, ω)} could be approximately transverse in Ω a [∇ • j s,nk (r, ω) ≈ 0, r ∈ Ω a ], for example, when Ω a is composed of a homogeneous and isotropic medium [38].Under such circumstances, the charge density may not be restored with the expansion based solely on this set.Another complementary set of sources can be introduced to make the expansion more complete.Here, we do not run into this issue further because characteristics of Bloch modes can be often grasped well with the set {j s,nk (r, ω)}.
Extracting the amplitude c nk (ω) of a basis source j s,nk (r, ω) in J s,k (r) confined in Ω a may be necessary in more detailed calculations such as the Purcell effect and spontaneous emission coupling factors [see Eq. (3a) and (3b)].Rather than the orthogonality relation for reciprocal cavities [38], we have to use the biorthogonality relation for the Bloch periodic parts ψ ψ ψ nk (r, ω) and ς ς ς s,nk (r, ω) for reciprocal PhCs.This biorthogonality condition also serves as a tool to uniquely determine the expansion coefficient c nk (ω).The derivation is presented in the appendix.
Following similar procedures to those for reciprocal cavities in Ref. [38], we obtain the white-source power spectrum P nk (ω), (resonance) frequency ω nk , and quality factor Q nk of mode (n, k).We first consider a source J s,k (r) that is proportional to j s,nk (r, ω) and demand the spatial integration of |J s,k (r)| 2 in Ω a to be frequency-independent (white source):  where a(ω) is a frequency-dependent amplitude; J is the spatially averaged strength of the source J s,k (r) and is frequency-independent; and V a is the volume of Ω a .The expression of |a(ω)| 2 can be obtained through the white-source condition in Eq. (8b).After the substitution of this expressions for |a(ω)| 2 in the generated power, the white-source power spectrum P nk (ω) is connected to the parameter Δε nk (ω) as follows: The (resonance) frequency ω nk is the one which minimizes the absolute value |ωΔε r,nk (ω)| of the denominator that appears inside the imaginary part of Eq. ( 9).The dispersion curve of ω nk versus k is taken as the band diagram of band n in this formulation.This curve would be compared with that obtained from the complex-ω method.With the identification of ω nk , the field f nk (r, ω nk ) is considered as (at least approximately) the Bloch mode of band n at k. Also, if one sets ω = ω nk in Eq. (5a), the imaginary part Im[Δε r,nk (ω nk )] corresponds to the amount of threshold gain required for the lasing of mode (n, k) in the active region Ω a .The effects of geometries of unit cells and active regions as well as the distribution of loss are incorporated in Im[Δε r,nk (ω nk )] automatically.
The magnitude P nk (ω nk ) at ω = ω nk is usually close to the peak value if the lineshape around ω nk is sufficiently narrow.Under these circumstances, we can approximate P nk (ω nk ) around ω nk with a Lorentzian.The ratio between the resonance frequency ω nk and full width at half maximum of the Lorentzian is the Q factor of mode (n, k).Its expression can be derived from the linear expansion of ωΔε r,nk (ω) around ω = ω nk [38] and is written as

Calculations and discussions
We utilize the 1D periodic bilayer structure in Fig. 2 to demonstrate the aforementioned formulation.Layer 1 (2) is the active (complement) region Ω a (Ω b ) with a width d 1 (d 2 ).The unit cell has a size a = d 1 + d 2 .Various combinations of materials will be considered.
The modes in this PhC can be analytically solved with transcendental equations.The derivation is briefly presented in section 3.1.The wave vector k = k x x + k can be decomposed into the component k x ∈ [−π/a, π/a] (the first BZ) perpendicular to layer planes and in-plane part k = k y ŷ + k z ẑ which is unlimited in magnitude due to infinitesimal in-plane periods.We calculate the parameters Δε r,nk (ω) of a few low-order bands in the first BZ for k x ∈ [0, π/a].In general, the in-plane wave vector k is set to the null vector except for some cases.The whitesource power spectra P nk (ω) of these modes will be depicted on the k x a/π − ωa/(2πc) plane for easy visualizations.The resonance frequencies ω nk obtained with minimum estimations of |ωΔε r,nk (ω)| are also plotted as a function of k x a/π on the same plane.Both ω nk and Q nk will be compared with those calculated from the complex-ω method.

Analytical solution of 1D periodic bilayer structure
For simplicity, we only take the wave propagating along the y direction into account and set k z = 0 since the structure is rotationally invariant in layer planes.In this way, the modes can be classified as either transverse-magnetic (TM) or transverse-electric (TE) to the y direction.We rewrite the mode label n = (s, l) as a two-tuple, where s indicates TM/TE, and l is the mode index within each category.Consistent with Eq. (5a), the GE problem is conveniently cast into Helmholtz equations for other cartesian field components in different regions: where φ nk (x), aside to a factor exp(ik y y), is proportional to the z component of the magnetic (electric) field for the TM (TE) mode; and ε a (ω) [ε b (ω)] is the relative permittivity in Ω a (Ω b ).
The field φ nk (x) can be separately solved in Ω a and Ω b as where A, B, C, and D are amplitudes to be determined.With Eq. ( 12), we demand the continuity of φ nk (x) at x = 0 and impose the phase BC at x = −d 1 and d 2 with a phase factor exp(ik x a).
In addition, y components of the electric field (TM case) and magnetic field (TE case) also satisfy the same continuity condition and phase BC.These constraints lead to the following characteristic equation in the matrix form: The determinant of of the matrix at the right-hand side of Eq. (13a) has to vanish so that nontrivial solutions of A, B, C, and D exist.This condition leads to a transcendental equation analogous to that of the Kronig-Penney model for electrons in a periodic rectangular potential [40]: For given ω, k x , and k y , parameters Δε r,nk (ω) are solved self-consistently using Eq. ( 14).These parameters are then processed to extract the information of interest.

Nondispersive and lossless structure of dielectric/dielectric bilayers
For the nondispersive and lossless dielectric/dielectric structure, we set ε a (ω) = 1 and ε b (ω) = 12.25.The two widths d 1 and d 2 are both equal to a/2.We do not specify the size a of the unit cell here because the nondispersive Maxwell's equations bring about scale-invariant band diagrams when depicted in the k x a/π − ωa/(2πc) plane.Also, since there is no dispersion, the complex-ω method in Eq. (1a) can be simplified into a transcendental equation similar to that in Eq. ( 14).The differences between two transcendental equations lie in that the parameter Δε r,nk (ω) is no longer present in the equation corresponding to the complex-ω method, and as a result of lossless dielectrics, it is real frequencies of Bloch modes that are solved.
Without the absorption, we project that the counterparts ω nk calculated from the proposed GE problem are identical to those obtained from the complex-ω method.Since no source is required to sustain the field in absence of the loss, we expect Δε r,nk (ω nk ) = 0, namely, the magnitudes |ωΔε r,nk (ω)| have minima equal to zero at resonances.In addition, the parameter Δε r,nk (ω) is real when ω is off resonances because the power P nk (ω) in Eq. ( 9) does not need to balance any loss.
Figure 3(a) shows |ωΔε r,n0 (ω)| at the BZ center as a function of the normalized frequency ωa/(2πc) for TE modes with l = 2 and 3.The resonance frequencies ω (TE,2),0 and ω (TE,3),0 are positions of the respective minima.For simplicity, at each frequency ω, we only show the smallest among magnitudes |ωΔε r,n0 (ω)| that are taken into account.The yellow lines in Fig. 3 mark adjacencies of neighboring modes in the frequency domain, at which the behavior of |ωΔε r,n0 (ω)| abruptly changes due to mode switching.As expected, the magnitude |ωΔε r,n0 (ω)| approaches zero as the frequency comes close the resonance frequency.We further carry out the above procedure to the two TE branches at k x ∈ [0, π/a] and k y = 0 and show the corresponding curves of |ωΔε r,nk (ω)| in Fig. 3(b).For clear illustrations, these curves are depicted in a 3D manner so that magnitudes |ωΔε r,nk (ω)| increase as the curves bend into the figure.The loci of minima (zeros) on each curve are projected onto the k x a/π − ωa/(2πc) plane.These loci correspond to the band diagrams calculated from the proposed GE problem.The counterparts from the complex-ω method are also shown for comparisons.The two sets of band diagrams are identical, as anticipated.Additional band diagrams using the two approaches at k y = π/(2a) are shown in Fig. 4. The two lowest-order TE and TM modes are taken into ac-  4. Comparisons between band diagrams of the two lowest-order TE and TM modes using the proposed formulation (solid lines) and complex-ω method (red circles) at k y = π/(2a) in the periodic structure of nondispersive and lossless dielectric/dielectric bilayers.count in this case.Again, the two sets of dispersion curves from the two schemes agree with each other.
In lossless structures, we can alternatively imagine that the parameter Δε r,nk (ω) acquires an infinitesimal imaginary part at the resonance.In this way, the power spectrum P nk (ω) in Eq. ( 9) is proportional to a delta function centered at ω nk , which is the sign of a persistent mode with an infinitely long lifetime.Accordingly, the Q factor Q nk in Eq. ( 10) also approaches infinity.

Nondispersive but lossy structure of dielectric/metal bilayers
For the nondispersive but lossy structure of dielectric/metal bilayers, we set ε a (ω) = 1 and ε b (ω) = −140 + 48i.The widths d 1 and d 2 are 0.99 a and 0.01 a, respectively.For this nondispersive structure, it is also possible to calculate band diagrams and Q factors with the complexω method [23] and the proposed GE problem using the transcendental equation in Eq. ( 14).
However, the mode frequencies from the complex-ω are no longer real.Under these circumstances, we take real parts of these complex frequencies as the mode frequencies.
The band diagrams of the four lowest-order TE modes calculated through the complex-ω method and the GE problem are shown in Fig. 5(a).The propagation constant k y is set to zero.From the inset, the two approaches lead to slightly but distinguishably different dispersion curves for the lowest-order TE mode (l = 1) near the BZ center around 0 ≤ k x a/π 0.3.Except for this special case, other dispersion curves calculated from two approaches nearly coincide with each other.From Fig. 5(b), the Q factors obtained through the two approaches are almost identical in most cases except for the lowest-order mode in the same range of k x a/π.Reflecting these differences, the modes with high (low) Q factors in Fig. 5(b) always correspond to narrow (wide) power spectra in Fig. 5(c).
A quick observation reveals that the agreements (disagreements) tend to take place when Q factors are high (low).In Table 1, the comparison of resonance frequencies and Q factors at the BZ center calculated from the two methods makes the correlation between Q factors and deviations more transparent.These observations indicate a connection between loss and frequency deviations.In fact, the origin of deviations may be inferred from the parameter Δε nk (ω nk ).At the resonance, the parameter Δε nk (ω nk ) is nearly imaginary, as can be observed from Table 1.In this case, Δε nk (ω nk ) reflects the amount of gain to be inserted into Ω a so that the metal loss Table 1.Comparisons between resonance frequencies as well as Q factors, and a list of Δε r,n0 (ω n0 ) for the four lowest TE modes at the BZ center in Fig. 5.

Complex-ω method [Eq. (1a)]
GE problem [Eq.(5a)] l ω n0 a/(2πc) is compensated, and the mode (n, k) can start to self oscillate (lase) at the real frequency ω nk .On the other hand, a huge gain may also shift the resonance significantly, which just reflects the deviated resonance frequencies obtained from the complex-ω method (passive structures) and GE problem (active structures) for the lowest TE mode.Therefore, the low-Q (high-loss) modes tend to bring about observable deviations between the two approaches while the high-Q (low-loss) modes do not.

Dispersive and lossy structure of dielectric/metal bilayers
For the periodic structure of dispersive and lossy dielectric/metal bilayers, we set ε a (ω) = 1 and adopt the Drude model ε b (ω) = 1 − ω 2 p /(ω 2 + iγω) for the dispersive metal, where the plasma frequency ω p is 10 15 rad•s −1 , and the damping term γ is set to 0.01ω p .The width d 1 and d 2 are set to 0.9 a and 0.1 a, respectively, with a = 1 μm.
Figure 6(a) and (b) shows the band diagrams and power spectra P nk (ω) of the three lowestorder TE modes, respectively.The calculations of the complex mode frequencies are implemented with the transfer-matrix method [21].The mode frequencies ω nk corresponding to the GE problem are still obtained from the transcendental equation in Eq. (14).Fairly good agreements for all the modes are present between calculations based on these two approaches due to decently high Q factors, as can be inferred from the narrow lineshapes of P nk (ω) in Fig. 6(b).
We further consider TM modes at k y = 0.9π/a because SPP waves propagating along layer planes are of TM characteristics.The band diagrams of the three lowest-order TM modes are shown in Fig. 7(a).The x component of the electric-field magnitude at the BZ center and boundary are also illustrated.From the profile of the field magnitude, the lowest-order TM modes exhibit maxima at the metal-dielectric interface, which is the feature of generic surface waves.On the other hand, although the magnitude profiles of the second TM branch (l = 2) still exhibit surface-like behaviors, the flat field profile and relatively high strength of the x electric field in the metal layer indicates ω ω p / √ 2 so that |ε b (ω)| 1.In a single dielectric/metal interface, such a surface wave should not exist, but the double interfaces and periodic structure help sustain these surface-wave-like Bloch modes.For this branch, the portion of the field strength in lossy metal does not change significantly as k x increases from the center to boundary.As a result, the lineshapes of corresponding power spectra P nk (ω) in Fig. 7(b) do not vary much, indicating that the Q factor is quite flat across the first BZ.In contrast, the third TM branch (l = 3) is not surface-wave-like in dielectric layers.The mode profiles there indicate that these modes resemble the guided mode of a metal/dielectric/metal waveguide.For this branch, the portion of the field strength in metal layers decreases significantly as k x varies across the first BZ.Accordingly, the corresponding power spectra P nk (ω) in Fig. 7(b) becomes sharper as k x comes closer to the BZ boundary, indicating a high Q factor of the branch there.

Conclusion
We have presented a formulation in the form of a GE problem to calculate power spectra, band diagrams, and Q factors of periodic structures.The approach is based on viewpoints of active photonic devices and can be applied to PhCs with arbitrary frequency dispersions and absorption.The calculations indicate that in most cases, outcomes from this formalism agree well with those from other approaches such as the complex-ω method.Slight deviations may occur in the low-Q regime.However, with decently high-Q factors, the formulation leads to reliable outcomes but is free from the complexity brought by frequency dispersions.In this way, the features of frequency dispersions can be easily grasped.

Appendix: Biorthogonality between Bloch modes
For a given wave vector k, the reciprocity theorem in Lorentz or Rayleigh-Carson form [41][42][43] is inapplicable among the corresponding basis fields because sources are present in each unit cell and does not vanish at infinity.Under these circumstances, we appeal to an alternative reciprocity theorem for Bloch periodic parts.With Faraday's law, we first define the magnetic field g nk (r, ω) from f nk (r, ω) as follows: where ϕ ϕ ϕ nk (r, ω) is the Bloch periodic part of the magnetic field g nk (r, ω).In terms of the Bloch periodic parts ψ ψ ψ nk (r, ω), ϕ ϕ ϕ nk (r, ω), and ς ς ς s,nk (r, ω), Maxwell's equations are rewritten as In reciprocal PhCs, the relative permittivity tensor is symmetric when represented in real and orthonormal coordinate bases, namely, εr (r, ω) = εT r (r, ω).With this property, we apply an analogous derivation of the general reciprocity theorem (integral form in a unit cell Ω uc ) [43] to two sets of Bloch periodic parts with labels (n, k) and (n, k ) and obtain In Eq. ( 17), the surface integral on the first line always vanishes because the integrand repeats on opposite sides of S uc while the corresponding outward normal vectors of differential surface areas (da) are antiparallel.If we further choose k = −k, the volume integral on the second line also becomes zero.These two conditions convert Eq. ( 17) into an analogy to the Rayleigh-Carson reciprocity for fields ψ ψ ψ nk (r, ω), ς ς ς s,nk (r, ω), ψ ψ ψ n −k (r, ω), and ς ς ς s,n −k (r, ω): From the second line of Eq. ( 18), if Δε r,n −k (ω) = Δε r,nk (ω), the volume integral of the dot product ψ ψ ψ nk (r, ω) • ψ ψ ψ n −k (r, ω) in Ω a has to be zero.The integral can be nonvanishing if the two parameters Δε r,n −k (ω) and Δε r,nk (ω) are identical.These characteristics motivate us to adopt the following biorthogonality relation after proper mutual orthogonalizations between Bloch periodic parts with wave vectors k and −k: where δ nn is the Kroneckers delta; Λ nk (ω) = Λ n−k (ω) is the complex normalization constant for f nk (r, ω) and f n−k (r, ω); and Θ nk (ω) = Θ n−k (ω) is the counterpart for j s,nk (r, ω) and j s,n−k (r, ω).Thus, in Eq. (3b), if the amplitude c nk (ω) of j s,nk (r, ω) in J s,k (r) is required, we should extract it by carrying out the integral of [exp(−ik • r)J s,k (r)] • ς ς ς s,n−k (r, ω) rather than [exp(−ik • r)J s,k (r)] • ς ς ς s,nk (r, ω) in Ω a because the Bloch periodic part ς ς ς s,nk (r, ω) would not necessarily be orthogonal to other counterparts with identical k, namely, Equation (20) uniquely determines the expansion coefficient c nk (ω) once a source J s,k (r) which is only present in Ω a and satisfies the phase BC is given.Still, if equation (19a) and (19b) needs to be utilized properly, the parameters from the set {Δε r,nk (ω)} have to repeat in {Δε r,n −k (ω)} and vice versa.It also means that the corresponding degenerate subsets in {Δε r,nk (ω)} and {Δε r,n −k (ω)} must have identical degeneracies.In this way, the biorthogonalization leading to Eq. (19a) and (19b) among various nondegenerate and degenerate sets can proceed without ambiguities.Otherwise, some degrees of freedom (amplitudes of basis sources) in the source J s,k could be left behind.Without going into more details, we assert that this property is expected in reciprocal (guiding) structures, in which the bidirectionality is present [44][45][46].Further understanding on this point shall be pursued in the future study.

Fig. 1 .
Fig. 1.The schematic diagram of the unit cell Ω uc in a generic periodic structure.The active region and its complement are denoted as Ω a and Ω b , respectively.

Fig. 2 .
Fig. 2. The unit cell of a 1D periodic bilayer structure.The widths of Ω a , Ω b , and the whole unit cell are d 1 , d 2 and a = d 1 + d 2 , respectively.

#(Fig. 3 .
Fig. 3. (a) The behaviors of |ωΔε r,n0 (ω)| as a function of the normalized frequency for TE modes with l = 2 and 3.The periodic structure is composed of nondispersive and lossless dielectric/dielectric bilayers.(b) The band diagrams from the loci of minima of |ωΔε r,nk (ω)| at k y = 0 (sold lines) and counterparts from the complex-ω method (red circles).For clear illustrations, the magnitudes |ωΔε r,nk (ω)| increase as the curves bend into the figure.
Fig. 4. Comparisons between band diagrams of the two lowest-order TE and TM modes using the proposed formulation (solid lines) and complex-ω method (red circles) at k y = π/(2a) in the periodic structure of nondispersive and lossless dielectric/dielectric bilayers.

#(Fig. 5 .
Fig. 5. (a) Band diagrams, (b) Q factors Q nk versus k x a/π, and (c) quasi-3D views of P nk (ω) as a function of k x a/π.The periodic structure is composed of nondispersive but lossy dielectric/metal bilayers.The propagation constant k y is set to zero.

Fig. 6 .
Fig. 6.(a) The band diagrams and (b) quasi-3D views of P nk (ω) of the three lowest TE modes in the periodic structure of dispersive and lossy dielectric/metal bilayers.The propagation constant k y is set to zero.

Fig. 7 .
Fig. 7. (a) The band diagrams for the three lowest-order TM modes and (b) lateral 3D views of P nk (ω) for TM modes with l = 2 and 3.The periodic structure is composed of dispersive and lossy dielectric/metal bilayers.The x component of the electric-field magnitude at the BZ center and boundary are also shown in (a).The propagation constant k y is set to 0.9π/a.