Full frequency-domain approach to reciprocal microlasers and nanolasers – perspective from Lorentz reciprocity

We develop a frequency-domain formulation in the form of generalized eigenvalue problems for reciprocal microlasers and nanolasers. While the goal is to explore the resonance properties of dispersive cavities, the starting point of our approach is the mode expansion of arbitrary current sources inside the active regions of lasers. Due to the Lorentz reciprocity, a mode orthogonality relation is present and serves as the basis to distinguish various cavity modes. This scheme can also incorporate the asymmetric Fano lineshape into the emission spectra of cavities. We show how to obtain the important parameters of laser cavities based on this formulation. The proposed approach could be an alternative to other computation schemes such as the finite-difference-time-domain method for reciprocal cavities. © 2011 Optical Society of America OCIS codes: (140.3430) Laser theory; (140.3945) Microcavities. References and links 1. 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 metallic-coated nanocavities,” Nat. Photonics 1, 589–594 (2007). 2. M. T. Hill, M. Marell, E. S. P. Leong, B. Smalbrugge, Y. Zhu, M. Sun, P. J. van Veldhoven, E. J. Geluk, F. Karouta, Y. S. Oei, R. Nötzel, C. Z. Ning, and M. K. Smit, “Lasing in metal-insulator-metal sub-wavelength plasmonic waveguides,” Opt. Express 17, 11107–11112 (2009). 3. M. P. Nezhad, A. Simic, O. Bondarenko, B. Slutsky, A. Mizrahi, L. Feng, V. Lomakin, and Y. Fainman, “Roomtemperature subwavelength metallo-dielectric lasers,” Nat. Photonics 4, 395–399 (2010). 4. K. Yu, A. Lakhani, and M. C. Wu, “Subwavelength metal-optic semiconductor nanopatch lasers,” Opt. Express 18, 8790–8799 (2010). 5. S. Kita, S. Hachuda, K. Nozaki, and T. Baba, “Nanoslot laser,” Appl. Phys. Lett. 97, 161108 (2010). 6. C. Y. Lu, S. W. Chang, S. L. Chuang, T. D. Germann, and D. Bimberg, “Metal-cavity surface-emitting microlaser at room temperature,” Appl. Phys. Lett. 96, 251101 (2010). 7. M. Nomura, Y. Ota, N. Kumagai, S. Iwamoto, and Y. Arakawa, “Zero-cell photonic crystal nanocavity laser with quantum dot gain,” Appl. Phys. Lett. 97, 191108 (2010). 8. B. Ellis, M. A. Mayer, G. Shambat, T. Sarmiento, J. Harris, E. E. Haller, and J. Vuckovic, “Ultralow-threshold electrically pumped quantum-dot photonic-crystal nanocavity laser,” Nat. Photonics 5, 297–300 (2011). 9. M. W. Kim and P. C. Ku, “Semiconductor nanoring lasers,” Appl. Phys. Lett. 98, 201105 (2011). 10. C. Y. Lu, S. L. Chuang, A. Mutig, and D. Bimberg, “Metal-cavity surface-emitting microlaser with hybrid metaldbr reflectors,” Opt. Lett. 36, 2447–2449 (2011). 11. K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. 14, 302–307 (1966). 12. A. Taflove and M. E. Brodwin, “Numerical solution of steady-state electromagnetic scattering problems using the time-dependent Maxwell’s equations,” IEEE Trans. Microwave Theory Tech. 23, 623–630 (1975). #152167 $15.00 USD Received 1 Aug 2011; revised 18 Sep 2011; accepted 27 Sep 2011; published 10 Oct 2011 (C) 2011 OSA 24 October 2011 / Vol. 19, No. 22 / OPTICS EXPRESS 21116 13. A. Taflove, “Application of the finite-difference time-domain method to sinusoidal steady state electromagnetic penetration problems,” IEEE Trans. Electromagn. Compat. 22, 191–202 (1980). 14. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. (Artech House, 2005). 15. 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). 16. S. H. Chang and A. Taflove, “Finite-difference time-domain model of lasing action in a four-level two-electron atomic system,” Opt. Express 12, 3827–3833 (2004). 17. Y. Huang and S. T. Ho, “Computational model of solid-state, molecular, or atomic media for FDTD simulation based on a multi-level multi-electron system governed by Pauli exclusion and Fermi-Dirac thermalization with application to semiconductor photonics,” Opt. Express 14, 3569–3587 (2006). 18. S. V. Zhukovsky and D. N. Chigrin, “Numerical modelling of lasing in microstructures,” Phys. Stat. Solidi B 244, 3515–3527 (2007). 19. S. V. Zhukovsky, D. N. Chigrin, A. V. Lavrinenko, and J. Kroha, “Switchable lasing in multimode microcavities,” Phys. Rev. Lett. 99, 073902 (2007). 20. S. V. Zhukovsky, D. N. Chigrin, and J. Kroha, “Bistability and mode interaction in microlasers,” Phys. Rev. A 79, 033803 (2009). 21. J. Jin, The Finite Element Method in Electromagnetics (Wiley and Sons, 2002). 22. K. Busch, M. König, and J. Niegemann, “Discontinuous Galerkin methods in nanophotonics,” Laser Photonics Rev., doi: 10.1002/lpor.201000045 (2011). 23. S. M. Spillane, T. J. Kippenberg, K. J. Vahala, K. W. Goh, E. Wilcut, and H. J. Kimble, “Ultrahigh-Q toroidal microresonators for cavity quantum electrodynamics,” Phys. Rev. A 71, 013817 (2005). 24. J. Wiersig, “Boundary element method for resonances in dielectric microcavities,” J. Opt. A: Pure Appl. Opt. 5, 53–60 (2003). 25. J. Wiersig, “Hexagonal dielectric resonators and microcrystal lasers,” Phys. Rev. A 67, 023807 (2003). 26. S. Y. Lee, M. S. Kurdoglyan, S. Rim, and C. M. Kim, “Resonance patterns in a stadium-shaped microcavity,” Phys. Rev. A 70, 023809 (2004). 27. M. S. Kurdoglyan, S. Y. Lee, S. Rim, and C. M. Kim, “Unidirectional lasing from a microcavity with a rounded isosceles triangle shape,” Opt. Lett. 29, 2758–2760 (2004). 28. T. Nobis and M. Grundmann, “Low-order optical whispering-gallery modes in hexagonal nanocavities,” Phys. Rev. A 72, 063806 (2005). 29. J. Wiersig, “Formation of long-lived, scarlike modes near avoided resonance crossings in optical microcavities,” Phys. Rev. Lett. 97, 253901 (2006). 30. 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). 31. B. Klein, L. F. Register, M. Grupen, and K. Hess, “Numerical simulation of vertical cavity surface emitting lasers,” Opt. Express 2, 163–168 (1998). 32. B. Klein, L. F. Register, K. Hess, D. G. Deppe, and Q. Deng, “Self-consistent Green’s function approach to the analysis of dielectrically apertured vertical-cavity surface-emitting lasers,” Appl. Phys. Lett. 73, 3324–3326 (1998). 33. 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). 34. 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). 35. E. I. Smotrova, A. I. Nosich, T. M. Benson, and P. Sewell, “Optical coupling of whispering-gallery modes of two identical microdisks and its effect on photonic molecule lasing,” IEEE J. Sel. Top. Quantum Electron. 12, 78–85 (2006). 36. L. D. Landau, E. M. Lifshitz, and L. P. Pitaevskii, Electrodynamics of Coninuous Media, 2nd ed. (ButterworthHeinemann, 1984). 37. C. A. Balanis, Advanced Engineering Electromagnetics, 1st ed. (Wiley and Sons, 1989). 38. 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). 39. N. Qureshi, H. Schmidt, and A. R. Hawkins, “Cavity enhancement of the magneto-optic Kerr effect for optical studies of magnetic nanostructures,” Appl. Phys. Lett. 85, 431–433 (2004). 40. N. Qureshi, S. Wang, M. A. Lowther, A. R. Hawkins, S. Kwon, A. Liddle, J. Bokor, and H. Schmidt, “Cavityenhanced magnetooptical observation of magnetization reversal in individual single-domain nanomagnets,” Nano. Lett. 5, 1413–1417 (2005). 41. Y. Arakawa and H. Sakaki, “Multidimensional quantum well laser and temperature dependence of its threshold current,” Appl. Phys. Lett. 40, 939–941 (1982). 42. H. Aoki, “Novel Landau level laser in the quantum Hall regime,” Appl. Phys. Lett. 48, 559–560 (1986). #152167 $15.00 USD Received 1 Aug 2011; revised 18 Sep 2011; accepted 27 Sep 2011; published 10 Oct 2011 (C) 2011 OSA 24 October 2011 / Vol. 19, No. 22 / OPTICS EXPRESS 21117 43. G. Scalari, C. Walther, L. Sirigu, M. L. Sadowski, H. Beere, D. Ritchie, N. Hoyler, M. Giovannini, and J. Faist, “Strong confinement in terahertz intersubband lasers by intense magnetic fields,” Phys. Rev. B 76, 115305 (2007). 44. A. Wade, G. Fedorov, D. Smirnov, S. Kumar, B. S. Williams, Q. Hu, and J. L. Reno, “Magnetic-field-assisted terahertz quantum cascade laser operating up to 225 K,” Nat. Photonics 3, 41–45 (2009). 45. G. Scalari, D. Turčinková, J. Lloyd-Hughes, M. I. Amanti, M. Fischer, M. Beck, and J. Faist, “Magnetically assisted quantum cascade laser emitting from 740 GHz to 1.4 THz,” Appl. Phys. Lett. 97, 081110 (2010). 46. S. W. Chang, C. Y. Lu, S. L. Chuang, T. D. Germann, U. W. Pohl, and D. Bimberg, “Theory of metal-cavity surface-emitting microlasers and comparison with experiment,” IEEE. J. Sel. Top. Quantum. Electron. (to be published). 47. A. Yariv, Optical Electronics in Modern Communications, 5th ed. (Oxford, 1996). 48. L. A. Coldren and S. W. Corzine, Diode Lasers and Photonic Integrated Circuits, 1st ed. (Wiley and Sons, 1995). 49. S. W. Chang, T. R. Lin, and S. L. Chuang, “Theory of plasmonic Fabry-Perot nanolasers,” Opt. Express 18, 15039–15053 (2010).


Introduction
There is significant progress in microlasers and nanolasers [1][2][3][4][5][6][7][8][9][10] as the top-down fabrication technology advances.The geometries and sizes of ultrasmall laser cavities can be well controlled, which makes the fine tuning of mode characteristics such as lasing wavelengths more accurate and flexible.Nevertheless, it often takes much effort and time to implement active photonic devices in the subwavelength regime.To save cost and time, the preliminary cavity calculation for performance estimations is often preferred before device fabrications.After characterization measurements, if the optimization of device structures is necessary, the more serious modeling has to be carried out.With these demands, a physical and efficient formulation for the modeling of microcavities and nanocavities is indispensable.
Two common approaches to the resonance modes of microcavities and nanocavities are the finite-difference-time-domain (FDTD) method [11][12][13][14] and complex eigenfrequency method [15].Both have their own advantages and disadvantages.The merits of the FDTD method, for instances, are the easy implementations on the Yee lattice [11], simultaneous spectral searches for resonance modes with broadband excitations, and the lower memory usage than those of other computation schemes.On the other hand, some drawbacks of the FDTD method may bring inconveniences into the cavity modeling.These shortcomings include (1) time-consuming computations for modes with high quality (Q) factors due to the stability issue from the size of time steps [12]; (2) poor mesh adaptivity to arbitrary cavity structures; (3) excitation-dependent outcomes of the mode patterns; and (4) few convenient ways to incorporate arbitrary frequency dispersions of materials into time-domain calculations.Albeit the drawbacks, the FDTD method has been applied in the modeling of active photonic devices [16][17][18][19][20].The scheme is particularly useful when the dynamic properties of lasers, which are not easily conceivable from other approaches, are of interest.
The complex eigenfrequency method resolves some of the issues in the FDTD method and is represented as a generalized eigenvalue problem from the source-free Maxwell's equations: where c is the speed of light in vacuum; ω cav and E cav (r) are the complex eigenfrequency and electric field of the cavity mode; and εr (r, ω cav ) is the relative permittivity tensor.The implementation of this generalized eigenvalue problem with the finite-element method (FEM) [21][22][23] removes the problem of mesh adaptivity in issue (2).In addition to the FEM, the eigenvalue problem can be also implemented using the integral-equation method [24][25][26][27][28][29].For issue (3), as an eigenvalue-type approach without sources, the cavity structure itself determines the mode profiles, and multimode excitations from arbitrary sources are avoided.
The complex eigenfrequency method is not flawless, however.For dispersive cavities, the eigenfrequency ω cav also comes into play in the relative permittivity tensor εr (r, ω cav ).Taking the dispersion into account implies that the eigenfrequency ω cav should be obtained iteratively and self-consistently, which is the penalty for issue (4).In addition, the relative permittivity tensor εr (r, ω) needs to be extended to the complex frequency ω.It is not always clear how this generalization is made theoretically or empirically if the experimental data at real frequencies are the only reliable sources of dispersions.Also, even though the Q factor is easily calculated from the real and imaginary parts of ω cav [Q = −Re[ω cav ]/2Im[ω cav ], where Im[ω cav ] < 0 for the time dependence exp(−iω cav t)], the complex eigenfrequency results in a nonphysical divergent far field due to the complex vacuum wave vector k 0 = ω cav /c (Im[k 0 ] < 0) and the outgoing-wave boundary condition [| exp(ik 0 r)| r→∞ = exp(−Im[k 0 ]r)| r→∞ → ∞, assuming the cavity is surrounded by vacuum].All these issues limit the applicability of the eigenfrequency method to the modeling of dispersive cavities.
Excitation sources are the origin of many critical issues mentioned above.At the level of classical electrodynamics, as long as the system loss, including the absorption and radiation ones, outnumbers the system gain (even slightly), the physically meaningful fields must be generated by some sorts of sources.For instance, a plane wave in the source-free region is actually generated by a current sheet somewhere else.The nonphysical divergent far field in the eigenfrequency method is also due to the absence of sources.In fact, the dispersion effect is easier to manage if carefully-constructed current sources sinusoidally oscillating at real frequencies are properly embedded into the cavity.
We present a full frequency-domain formulation for reciprocal microlasers and nanolasers and demonstrate how to obtain important parameters of laser cavities.The central part of this approach is a generalized eigenvalue problem which resembles that of self-supporting lasing modes [15,[30][31][32][33][34][35] (a set of discrete real frequencies and corresponding threshold gains are searched in this case).However, unlike the frequency-domain schemes including the eigenfrequency method and approach of self-supporting modes, our motivation is the mode expansion of arbitrary sinusoidal current sources inside the active regions of lasers.The proposed formulation is more advantageous than the conventional schemes in the following respects: (1) no divergent far fields when compared with the eigenfrequency method because the source is present; (2) no multimode excitations when compared with the FDTD method because this approach is a generalized eigenvalue problem inherently; (3) the ability to directly include arbitrary frequency dispersions of materials when compared with all the schemes mentioned above because the frequency ω is real and given in the first place; (4) the straightforward mode expansion due to a natural orthogonality relation brought by the Lorentz reciprocity; and (5) the capability of modeling spectral properties of modes, including the asymmetric Fano lineshape of the emission spectra, which is nontrivial in the eigenfrequency method and approach of selfsupporting modes.These advantages make the formulation an alternative to other computation schemes such as the FDTD method in the cavity modeling.
The rest of the paper is organized as follows.We first briefly review the reciprocity theorem and introduce the concept of reciprocal cavities and lasers (section 2).The construction of the formulation is then discussed in detail (section 3 and appendix).We also demonstrate how to obtain the important parameters of microlaser and nanolaser cavities such as resonance frequencies, quality factors, threshold gains, and spontaneous emission coupling factors (section 4).The application of this formulation to a one-dimensional (1D) Fabry-Perot (FP) cavity is illustrated as an example (section 5).Additional remarks on the implementation and variation of the formulation are addressed (section 6), and a conclusion is given at the end (section 7).

Reciprocity theorem and reciprocal cavities/lasers
In this section, we briefly introduce the concept of (Lorentz) reciprocity and reciprocal cavities/lasers.The Lorentz reciprocity plays an important role in the construction of the formula- tion.Many of the results in the following sections need modifying if the Lorentz reciprocity does not hold.We will address the applicability of the formulation.
The schematic diagram of a microlaser or nanolaser cavity is depicted in Fig. 1.The active region of the laser (might be unconnected) is denoted as Ω a and is embedded in the cavity region, which could be an open structure and need not have physical boundaries from the surrounding.S a is the surface of Ω a (or union of the surfaces from unconnected parts).We consider a cavity which is characterized by the relative permittivity tensor εr (r, ω) but a relative permeability of unity.We also construct a region Ω which contains Ω a but does not necessarily cover the cavity region.The surface of Ω is denoted as S.
For two current sources J s,1 (r) and J s,2 (r) confined in Ω a but vanishing elsewhere, the reciprocity theorem states that if responses of the material system to any electromagnetic fields are linear, and the relative permittivity tensor εr (r, ω) in the cartesian basis is symmetric inside Ω: the fields generated by the respective sources satisfy the following integral identity [36,37]: where E 1 (r) and H 1 (r) are the electric and magnetic fields generated by J s,1 (r) alone; and E 2 (r) and H 2 (r) are the counterparts generated by J s,2 (r).Note that the symmetric form in Eq. ( 2) is not restricted to the cartesian basis.As long as Eq. ( 2) holds, εr (r, ω) is symmetric in any real orthogonal local bases inside Ω.If εr (r, ω) is nonsymmetric, an additional volume integral over Ω with an integrand proportional to E 1 (r) • [ εr (r, ω) − εT r (r, ω)]E 2 (r), where the superscript "T" means transpose, is present at the left-hand side of Eq. (3).
The surface S belongs to a region Ω containing Ω a and therefore encloses Ω a .Let us assume that far away from the cavity is the isotropic free space or absorptive media/structures.In this way, the surface integral at the right hand-side of Eq. ( 3) vanishes as a result of extending the surface S to infinity, at which the two cross products nearly cancel each other due to the plane-wave approximation of outgoing waves in the far-field regime, or the fields just turn exponentially small due to the absorption present at infinity.If other outer spaces can lead to a vanishing surface integral in Eq. (3), they can also be considered.Under such circumstances, the general reciprocity theorem in Eq. (3) turns into the Lorentz reciprocity theorem [38]: The Lorentz reciprocity theorem is more restrictive than the general reciprocity theorem because the relative permittivity tensor has to be symmetric everywhere (Ω is the full space now).A reciprocal cavity is a cavity to which the Lorentz reciprocity in Eq. ( 4) is applicable for two arbitrary current sources confined in Ω a .Practical cavities which exhibit a symmetric permittivity tensor in Eq. (2) everywhere, for example, dielectric spheres, cleaved ridge wavegeuides, and microdisks (assuming that outside the substrate is the free space), all belong to this type.For lasers, the nonlinearity such as the optical feedback from the stimulated emission is impermissible in the reciprocity theorem.Therefore, we adopt a loosened definition for reciprocal lasers, by which the effective permittivity tensor dressed by the nonlinearity remains symmetric, and the integral form in Eq. ( 4) stays valid for two weak current sources (dropping the perturbation to gain dynamics) inside Ω a .
The Lorentz reciprocity in Eq. ( 4) is critical to the formulation because it introduces a natural orthogonality relation between modes and provides a way to distinguish them as well as extract their magnitudes from an arbitrary source distribution in Ω a (see section 3.2).Its failure leads to the inapplicability of the proposed formulation to some cavity and laser systems.If the (effective) permittivity tensor εr (r, ω) inside or outside the cavity becomes nonsymmetric and breaks the Lorentz reciprocity, the formulation should not be applied, though some qualitative estimations can be still made when the asymmetry is small.The nonsymmetric permittivity tensor εr (r, ω) can take place when the time-reversal symmetry is broken.Such examples include cavities with the magneto-optic effect [39,40] or an external magnetic field [41][42][43][44][45].Still, most of the microlasers and nanolasers nowadays remain reciprocal during operation, namely, the effective permittivity tensors remain symmetric.Thus, this approach can be useful in a wide range of active photonic devices in the subwavelength regime.

Formulation of reciprocal cavities
In the following subsections, we will describe the generalized eigenvalue problem, mode orthogonality relation, and dyadic Green's function for the sources localized in Ω a .We note that unless fixed by the real device layout, the active region should have a higher (identical) symmetry group than (to) that of the cavity structure (if there is any) so that the symmetry is preserved.

Generalized eigenvalue problem
In the frequency domain, the Maxwell's equations in the presence of sources are where E(r) and H(r) are the electric and magnetic fields; ε 0 and μ 0 are the vacuum permittivity and permeability, respectively; and J s (r) is the current source.Both the effects of absorption (cold cavity, but inter-state dipole absorption might be excluded) and gain (warm cavity) can be incorporated into εr (r, ω), depending on the operation condition of the cavity.Also, in Eq. ( 5a) and (5b), the real frequency ω is given and can be continuously varied.It is different from the so-called resonance frequencies of cavity modes, which are a set of discrete real frequencies to be sought with specific criteria.
After eliminating H(r) in Eq. ( 5a) and (5b), we obtain the wave equation for E(r) We now consider the current sources J s (r) that are only present in Ω a , namely, n Re[ r,n ( )] shifts the resonance frequency ω n to ω, while the imaginary part Im[Δε r,n (ω)] compensates the loss and converts P n (ω ) into a delta function centered at ω.
The constraint in Eq. ( 7) is often physical because the spontaneous emission dipole moments which trigger the lasing field usually coexist with gain in the active region only.Our goal is to find a set {j s,n (r, ω)} of current sources, which is present only in Ω a and labeled by index n, for the mode expansion of J s (r) in Eq. (7).Denote the set of electric fields generated by {j s,n (r, ω)} as {f n (r, ω)}.We want to prevent the multi-mode excitations from carelesslyconstructed sources.A solution to this issue is to excite the electric field with a source which is proportional to the electric field itself in Ω a but vanishes elsewhere, namely, a self-duplicated vector field in Ω a .Thus, we assign the following ansatz to j s,n (r, ω): where U(r) is the indicator function for Ω a ; and Δε r,n (ω) is a complex parameter.With Eq. (8a), the wave equation of f n (r, ω) is transformed into a generalized eigenvalue problem: where (ω/c) 2 Δε r,n (ω) acts as the eigenvalue; and mode quantization indicated by index n is justified from the source confinement in Ω a , cavity structure, and outgoing-wave boundary condition.Note that f n (r, ω) is not divergent in the far-field zone because it is effectively generated by a source j s,n (r, ω).Once f n (r, ω) is obtained, the corresponding magnetic field g n (r, ω) is derived from Faraday's law in Eq. (5a): The parameter Δε r,n (ω) acquires its frequency dependence when Eq. ( 9) is solved for different ω's.Its physical interpretation becomes clear if we rewrite Eq. ( 9) as where Ī is the identity tensor.Equation ( 11) resembles the source-free wave equation of which the solution f n (r, ω) oscillates at the real frequency ω, namely, a self-supporting mode.As indicated in Fig. 2, the real part Re[Δε r,n (ω)] can be regarded as the average permittivity variation in Ω a required to shift the original resonance frequency ω n of mode n to the given frequency ω, while the imaginary part Im[Δε r,n (ω)] is related to the necessary gain which compensates the loss and converts the emission spectrum P n (ω ) (ω is a dummy variable for the frequency) into the delta-function spectrum centered at ω [for details about ω n and P n (ω ), see section 4.1].The self-supporting mode is a bonus of this formulation though the original motivation is aimed at current sources.We also note that Im[Δε r,n (ω)] has to be negative [in the convention exp(−iωt)] because it represents gain.The magnitude |Im[Δε r,n (ω)]| is smaller if the more gain is incorporated into εr (r, ω), but Im[Δε r,n (ω)] never turns positive under physical circumstances.If too much gain is present initially, the stimulated emission would clamp it so that the steady-state gain still corresponds to a negative imaginary part Im[Δε r,n (ω)].
On the other hand, some current sources are not expandable with the set {j s,n (r, ω)}.This incompleteness is more evident for a homogeneous and isotropic active region characterized by a scalar permittivity ε r,a (ω).In this case, we take the divergence in Eq. ( 11) for r ∈ Ω a : Since Δε r,n (ω) = −ε r,a (ω) [otherwise, the fact that ∇ × ∇ × f n (r, ω) ∝ ∇ × g n (r, ω) = 0 in Ω a may lead to discontinuous tangential magnetic fields across S a , and thus fictitious surface currents], the divergence of the field f n (r, ω) vanishes in Ω a .Therefore, the corresponding current source j s,n (r, ω) is also divergenceless (solenoidal) in Ω a : From Eq. ( 13), any current sources confined in Ω a that contribute to volume charge densities oscillating at ω (nonzero divergence) cannot be fully expanded by the set {j s,n (r, ω)} (however, j s,n (r, ω) may result in the surface charge density on S a ).Although real-space charge densities oscillating around the resonance frequencies of cavity modes are uncommon in typical lasers, these charge densities, once induced, affect both the near-field and far-field profiles and should be taken into account.For general reciprocal permittivity tensors in Ω a , the situation becomes less transparent because the set {j s,n (r, ω)} may not be divergenceless.However, from the lesson of homogeneous and isotropic active region, it is probable that {j s,n (r, ω)} does not span all the source configurations.Therefore, we introduce an analogous set {i s,m (r, ω)} of current sources, where m is the mode index, to complement {j s,n (r, ω)}.Similar to j s,n (r, ω), we require i s,m (r, ω) to be confined in Ω a only.The corresponding sets of electric and magnetic fields generated by {i s,m (r, ω)} are denoted as {u m (r, ω)} and {w m (r, ω)}, respectively, and they also have to satisfy the outgoing-wave boundary condition.We summarize the construction of these additional sets in the appendix.
With the sets {j s,n (r, ω)} and {i s,m (r, ω)}, we expand an arbitrary current source J s (r) confined in Ω a and oscillating at ω as follows: where c n and d m are the expansion coefficients.From the superposition principle of linear systems, the corresponding electric field E(r) and magnetic field H(r) are expressed as Our next goal is the extractions of expansion coefficients c n and d m .This step is straightforward if various modes are orthogonal to each other via a certain form of inner product.We will show that the Lorentz reciprocity provides such a handy relation.

Mode orthogonality
From the Lorentz reciprocity in Eq. ( 4), we are now in a position to show how a natural orthogonality relation between various modes can be derived.With the sets {j s,n (r, ω)} and {f n (r, ω)} constructed in section 3.1, we make the following assignments: and substitute them into the Lorentz reciprocity theorem in Eq. ( 4): In Eq. ( 16), if Δε r,n (ω) = Δε r,n (ω), the volume integral of the dot product f n (r, ω) • f n (r, ω) must vanish.On the other hand, if a degeneracy exists such that Δε r,n (ω) = Δε r,n (ω), we can still utilize the volume integral of the dot product as the inner-product rule to orthogonalize the modes.Thus, we obtain a natural orthogonality relation for the set {f n (r, ω)} as follows: where δ n n is the Kronecker's delta; and Λ n (ω) is the complex normalization constant of f n (r, ω).We also define an analogous orthogonality relation to Eq. ( 17) for the set {j s,n (r, ω)}: where Θ n (ω) is the normalization constant of j s,n (r, ω).
For the set {i s,m (r, ω)}, we also demand a similar orthogonality relation based on the same inner-product rule: where Ξ m (ω) is the normalization constant of i s,m (r, ω).In addition, two current sources, one from {j s,n (r, ω)} and the other from {i s,m (r, ω)}, always have to be orthogonal to each other: The conditions in Eqs.(18a), (19), and (20) necessitate only one single orthogonality relation for any basis vector functions from {j s,n (r, ω)} and {i s,m (r, ω)}.This property, however, stems from a specifically-constructed set {i s,m (r, ω)}, which is summarized in the appendix.
With the orthogonality relations in Eqs.(18a), (19), and (20), we can extract the expansion coefficients c n and d m in Eq. (14a), (14b), and (14c).The expressions of these coefficients can help construct the dyadic Green's function, which is useful in the calculations of spontaneous emission coupling factors.

Dyadic Green's function
For an arbitrary current source J s (r) confined in Ω a , we need to find its mode expansion coefficients c n and d m in Eq. (14a) in order to reconstruct the electric field E(r) and magnetic field H(r) from Eq. (14b) and (14c), respectively.We first dot product both sides of Eq. (14a) with the current source j s,n (r, ω) and integrate it over Ω a .With the orthogonality relations in Eqs.(18a), (19), and (20), only the term corresponding to c n at the right-hand side of Eq. (14a) remains, and we obtain the expression of c n as follows: With the same procedure but adopting i s,m (r, ω) rather than j s,n (r, ω), we derive the analogous expression of d m : For the applications in optics, the spatial profile and far-field pattern of the electric field E(r) are important.With the expressions of expansion coefficients in Eq. ( 21a) and (21b), we substitute them into the expansion series of E(r) in Eq. (14b).After renaming the dummy index n (m ) into n (m) and interchanging the variable r with r , we then link the electric field E(r) to the current source J s (r) through the dyadic Green function Ḡee (r, r , ω): The matrix multiplications j T s,n (r , ω)J s (r ) and i T s,m (r , ω)J s (r ) in the tensor operation of Ḡee (r, r , ω) on J s (r ) represent the dot products j s,n (r , ω) • J s (r ) and i s,m (r , ω) • J s (r ), respectively.On the other hand, not every current source can be substituted into Eq.(22a).The dyadic Green's function Ḡee (r, r , ω) in Eq. (22b) is not applicable to the sources which spread outside the active region Ω a .
We will utilize the expression of the dyadic Green's function Ḡee (r, r , ω) in Eq. (22b) to calculate the spontaneous emission coupling factor in section 4.2.

Characteristics of reciprocal cavities
In section 3, we constructed the sets of current sources {j s,n (r, ω)}, electric fields {f n (r, ω)}, and magnetic fields {g n (r, ω)} at a given real frequency ω.To further investigate the spectral properties of these modes, we need to link the calculations at different ω's together.Through this connection, we can then define the cavity resonance.
With the field f l (r, ω) of a nondegenerate mode l at ω, the formal way to associate it with its counterpart f l (r, ω + Δω) at ω + Δω, where Δω is a small frequency difference, is to track how the field profile evolves from ω to ω + Δω and make a one-to-one link from {f n (r, ω)} to {f n (r, ω + Δω)}.An easier approach is the connection through eigenvalue Δε r,l (ω), assuming a continuous and smooth frequency variation.On the other hand, in degenerate cases due to the cavity symmetry, once a field profile f l (r, ω) in the set {f l (r, ω)|Δε r,l (ω) = Δε r,l (ω)} of degenerate modes at ω (excluding the accidental degeneracy) is chosen, its counterpart f l (r, ω + Δω) has to satisfy a condition similar to the orthogonality relation in Eq. ( 17): Unlike Eq. ( 17), however, Eq. ( 23) originates from the symmetry viewpoint even though two equations coincide with each other at Δω = 0. Note that in degenerate cases, the continuity and smoothness of Δε r,l (ω) may only provide the group-to-group rather than one-to-one correspondence of modes and are insufficient to uniquely identify a particular mode.
In the following derivations, we will assume that the proper one-to-one connection has been established for mode l at different frequencies.In this way, we can identify one of these frequencies as its resonance frequency and obtain the mode information from there.

Resonance frequency, lineshape, quality factor, and threshold gain
With the identification of mode l at each frequency, we look into the frequency dependence of the power generated by a current source proportional to j s,l (r, ω).In this way, the forms of the current source J s (r) and electric field E(r) are where a(ω) is the source strength.For a fair comparison between the responses at different frequencies, we demand a frequency-independent volume integral of where J is a real constant; and V a is the volume of the active region.The constraint in Eq. ( 25) is the white-noise condition [46] for mode l so that the frequency-dependent strength a(ω) of J s (r) does not interfere with the intrinsic nature of mode l on the power spectrum.In this way, the expressions of the square magnitude |a(ω)| 2 and power P l (ω) become Equation (26b) indicates that the lineshape of the white-noise power is determined by the term Im{[ωΔε r,l (ω)] −1 }.To understand its behavior, we define a frequency ω l such that the absolute value |ω l Δε r,l (ω l )| is the minimum.As shown in Fig. 3(a), if we plot the locus of η(ω) ≡ ωΔε r,l (ω) parameterized by ω on the complex η plane, the differential change δ η due to a small frequency variation δ ω around ω ≈ ω l , when viewed as a two-dimensional (2D) vector, has to be perpendicular to that of η(ω l ).If these two vectors were not perpendicular, |η(ω)| would have further reduced its magnitude by either moving forward or backward on the curve.In terms of complex numbers, the illustration in Fig. 3(a where Δω l is a parameter that must be real.The presence of imaginary number i in Eq. ( 27) indicates a ±π/2 phase change (depending on the sign of δ ω), namely, the 2D vectors of δ η and η(ω l ) on the complex η plane are perpendicular.From Eq. ( 27) and δ η η (ω l )δ ω, where η (ω) is the frequency derivative of η(ω), we can express the parameter Δω l as With Eq. ( 27), we approximate η(ω) near ω l with the linear expansion in ω − ω l and substitute it into the power spectrum P l (ω) in Eq. (26b): (29b) Equation (29b) indicates that P l (ω) has a Fano lineshape near ω l (weighted sum of the Lorentzian centered at ω l and its Hilbert transformation).As shown in Fig. 3(b), this lineshape is asymmetric with respect to ω l .In contrast to the Lorentzian, which has a maximum at ω l , the maximum of Fano lineshape is blueshifted (Re[ω l Δε r,l (ω l )] > 0) or redshifted (Re[ω l Δε r,l (ω l )] < 0).Although the asymmetry on the power spectrum is derived based on the white-noise source, the phenomenon is generic to most frequency-dependent sources.
If |Re[ω l Δε r,l (ω l )]| |Im[ω l Δε r,l (ω l )]| (valid for most cavities with high Q factors), the physical interpretations of ω l and Δω l are more evident.In this case, the frequency ω l is nearly the peak frequency of P l (ω).Therefore, we may view ω l as the resonance frequency of mode l and generalize it to the cases of asymmetric lineshapes.The electric field f l (r, ω l ) at ω l is then the field profile of mode l.Correspondingly, we can identify the parameter Δω l as the full-width-at-half-maximum (FWHM) linewidth of the approximate Lorentzian.The (cold-or warm-cavity) quality factor Q l of mode l is defined as the ratio between ω l and Δω l : Recall that Δε r,l (ω l ) is the permittivity variation required for the self-supporting mode l at ω l .With a homogeneous and isotropic active region [ εr (r, ω)| r∈Ω a = ε a (ω) Ī] and the cold-cavity condition [inter-state dipole absorption excluded from ε a (ω)], the threshold gain g th,l is

Spontaneous emission coupling factor and Purcell effect
The spontaneous emission coupling factor β l is defined as the ratio of the spontaneous emission power into mode l over the total spontaneous emission power.Classically, various spontaneous emission events are generated by uncorrelated dipoles.Therefore, the spontaneous emission current source J sp (r, ω) per unit square-root energy [46] at ω satisfies the following relation: where α, α = x, y, z; D cv (r, ω) is the position-dependent emitter density per unit energy due to transitions between states (bands) c and v, which includes the information of state occupations and vanishes for r / ∈ Ω a ; j sp,cv (ω) is the microscopic source responsible for the spontaneous emission; . . .means the ensemble average; and overline means averages or expectations of other types such as crystal symmetries.For dipole transitions, the source j sp,cv (ω) is where qd cv (ω) is the dipole moment; and a factor of two in Eq. ( 33) is due to the consistent usages between the harmonic sum and phasor notation ( With Eq. ( 32), the spontaneous emission rate r sp (ω) per unit energy is expressed as where E(r, ω) is the electric field generated by J sp (r, ω); and the factor of 1/(hω) is to convert the power into rate.We then substitute the expansion of the dyadic Green's function Ḡee (r, r, ω) in Eq. (22b) into Eq.( 34) and obtain where r sp,n (ω) is the spontaneous emission rate per unit energy into mode n in {f n (r, ω)}; and rsp,m (ω) is the counterpart into mode m in {u m (r, ω)}.
The total spontaneous emission rate R sp , spontaneous emission rate R sp,n into mode n in {f n (r, ω)}, and the counterpart Rsp,m into mode m in {u m (r, ω)}, are the integrations of r sp (ω), r sp,n (ω), and rsp,m (ω) over photon energy hω, respectively: and the spontaneous emission coupling factor β l into mode l in {f n (r, ω)} is written as For a single dipole [see Eq. ( 33)] due to a two-level system (the occupied excited state |c , empty ground state |v , and transition frequency ω cv ) and located at position r s ∈ Ω a , namely, the spontaneous emission rate R sp,l into mode l in Eq. (36b) becomes Equation ( 39) is indicative of the Purcell effect on a single dipole.The spatial enhancement comes from the modal strength f l (r s , ω cv ) at position r s , and spectral enhancement originates from the lineshape effect from 1/[ωΔε r,l (ω cv )] if ω cv is sufficiently close to ω l .To compare it with the counterpart W sp,l from Fermi's Golden rule (Lorentzian density of states of mode l): where Êph,l (r s ) is the single-photon field of mode l, we rewrite R sp,l in Eq. ( 39) with a few simplifications.First, for ω cv ≈ ω l , the fastest-varying part with respect to ω cv − ω l is the factor 1/[ω cv Δε r,l (ω cv )].Except for this factor, we replace ω cv with ω l in other parts of R sp,l .Second, if we set Λ l (ω l ) to a positive real number and assume that the phase of f l (r, ω l ) does not vary rapidly in Ω a (except for ±π jumps near nodes or nodal surfaces), f l (r, ω l ) is close to a real vector field in Ω a , namely, f * l (r, ω l ) ≈ f l (r, ω l ) for r ∈ Ω a .Third, we assume that the asymmetry of the lineshape is minor and only keep the symmetric Lorentzian.With these simplifications, the form of R sp,l is reduced to that of W sp,l in Eq. (40) with the following connection between the single-photon field Êl (r s ) and mode profile f l (r s , ω l ):

One-dimensional Fabry-Perot cavity
To see if the proposed formulation leads to reasonable outcomes, we apply it to a onedimensional Fabry-Perot cavity, of which many analytical or approximate properties are known.  of mode n can be estimated from the nature of a FP resonator and represented as [47][48][49] where m n is the number of standing waves of mode n in the cavity; v g is the group velocity of the wave in the cavity and is identical to the phase velocity (c/ √ ε c ) in this case; and Γ z,n = (L a /L c )[1 ± sinc(ω n L a /c)] is the longitudinal confinement factor [+ (−) for even (odd) modes; and sinc(x) = sin(x)/x].The quantization of ω n in Eq. (43a) is due to an integral number of standing waves in the cavity region while Q n in Eq. (43b) is estimated from the fractional loss of the energy due to the power leakage at the two outputs in a round-trip period [47].The threshold gain g th,n in Eq. (43c) is obtained from the round-trip balance condition, taking into account the effects of L a = L c and the standing-wave pattern, as described by the expression of the longitudinal confinement factor Γ z,n [48].In Table 1, we show the theoretical estimations of hω n , Q n , and g th,n and their FP counterparts from Eq. (43a) to (43c).The spectral parameters hω n , Q n , and g th,n obtained from these two different approaches agree well.Therefore, we believe that the proposed approach has grasped the essential points of the cavity modeling.Figure 5(a) shows the square field magnitude of the even mode at a resonance photon energy of about 1.42 eV (n = 3 in Table 1).In addition to the standing-wave pattern of typical FP modes, the effect of the permittivity variation Δε r,3 (ω) can be observed in the active region (|z| < 1 μm).The growing envelop of the mode profile toward both ends of the active region indicates that the amplification in the active region compensates the radiation loss at two outputs of the cavity.Figure 5(b) shows the locus of hη(ω) = hωΔε r,3 (ω) on the complex hη plane.The permittivity variation Δε r,3 (ω 3 ) at resonance is about 0.00256-0.286i.Although the real part of the permittivity variation is much lower than the imaginary part in magnitude, this small but positive number indicates that the white-noise lineshape of this FP mode is a little bit asymmetric, and its peak photon energy is slightly blueshifted from hω 3 = 1.42 eV.
To make the phenomenon of asymmetry more significant, we lower the relative cavity permittivity to one quarter of its original value (ε c = 3.0625) and show the lineshape of an even mode l (solid black) in Fig. 6(a) (hω l = 1.281 eV).For comparison, we also show the lineshape of an even mode l (dashed red) with a close resonance photon energy when ε c = 12.25 (hω l = 1.278 eV).The square magnitudes of the mode profiles corresponding to the two cases are shown in Fig. 6(b).The lower relative cavity permittivity makes the radiation loss larger.Therefore, the quality factor of mode l is lower than that of mode l (Q l = 21.52 versus Q l = 96.34).The larger radiation loss of mode l also leads to the more rapid field growth in the active region, as shown in Fig. 6(b).Usually, the asymmetry on the lineshape is more significant in the more lossy or leaky cavities.From Fig. 6(a), the asymmetry on the lineshape of mode l is indeed more prominent than that of l .This behavior is also reflected in the more significant real part of the permittivity variation for mode l at its resonance frequency [Δε r,l (ω l ) = −0.0733− 0.359i versus Δε r,l (ω l ) = −0.00190− 0.311i].In addition, a plateaulike feature takes place between photon energies 1.4 to 1.45 eV, which even goes beyond the applicability of Fano lineshape in Eq. (29b).

Additional remarks
A practical issue is how to numerically implement the computation of the mode profiles.We focus on the set {f n (r, ω)} and show a generic computation domain based on the FEM in Fig. 7. Since we require the outgoing-wave boundary condition for the modes, perfect-matched layers (PMLs) are inserted in the inner sides of the computation domain to avoid unnecessary reflections.With fields significantly attenuated inside PMLs, the boundary conditions corresponding to the perfect electric conductor (PEC) or perfect magnetic conductor (PMC) can be imposed at the outer boundaries of the computation domain.The generalized eigenvalue problem is then implemented in an analogous numerical scheme to the mode calculations of cavities with the outer PEC or PMC boundaries.The only difference is that it is the permittivity variation Δε r,n (ω) rather than the eigenfrequency that is to be sought.The construction should be applicable to the modeling of most laser cavities as long as the Lorentz reciprocity holds.The computation time is also critical in scientific computations.The proposed formulation is based on a generalized eigenvalue problem solved at each sampling frequency.One can speed up the calculations with (1) efficient solvers for generalized eigenvalue problems, and (2) fewer frequency samplings.While (1) depends on details of solver implementations, (2) can be often improved with suitable estimation schemes.For example, the search for the resonance frequency ω l of mode l is equivalent to the minimization of the goal function |ωΔε r,l (ω)| 2 with respect to frequency ω.The convergence to the target frequency can be effective and robust using proper schemes such as quasi-Newton methods or conjugate gradient method.
Another issue is the ansatz of the current source.We may alternatively write the source as j s,n (r, ω) = −iωε 0 γ n U(r) T (r, ω)f n (r, ω), where γ n is related to the eigenvalue; and T (r, ω) is a symmetric tensor in the cartesian basis.The orthogonality relations take the form of volume integrals of f T n (r, ω) T (r, ω)f n (r, ω) and j T s,n (r, ω) T −1 (r, ω)j s,n (r, ω) in Ω a .As an example, we can use T (r, ω) = εr (r, ω) and directly set i s,m (r, ω) ∝ U(r) ε−1 r (r, ω)∇Φ m (r, ω) [Φ m (r, ω) is a scalar function vanishing at r ∈ S a ] without the projection process (see appendix).However, in this case, the physical interpretation of γ n related to the eigenvalue is less clear, and we should not pursue the related formulation here.

Conclusion
We have presented a frequency-domain approach to reciprocal microlasers and nanolasers.The motivation of the approach is the mode expansion of arbitrary current sources in the active region.The proposed formulation can take the frequency dispersion of materials into account and avoid the effect of multimode excitation.We have shown how to obtain the important spectral parameters of cavity modes using this approach, and demonstrated that the asymmetric Fano lineshape comes into play naturally in this formulation.With the constructed dyadic Green's function, we also derive the spontaneous emission coupling factor and Purcell effect on a single dipole.This approach may be adopted as an alternative to other computation schemes such as the FDTD method.

Fig. 1 .
Fig. 1.The schematic diagram of the laser cavity.The active region is denoted as Ω a , and Ω is an arbitrary region which contains Ω a .S a and S are the surfaces corresponding to Ω a and Ω, respectively.

Fig. 2 .
Fig.2.The effect of Δε r,n (ω) on the emission spectrum P n (ω ).The real part Re[Δε r,n (ω)] shifts the resonance frequency ω n to ω, while the imaginary part Im[Δε r,n (ω)] compensates the loss and converts P n (ω ) into a delta function centered at ω.

Fig. 3 .
Fig. 3. (a)The locus of η(ω) = ωΔε r,l (ω) parameterized by ω on the complex η plane.At resonance frequency ω l , η(ω l ) is closest to the origin of the complex plane.(b) The comparison between the Lorentzian and Fano lineshapes.The Fano lineshape is asymmetric with respect to ω l , and the peak is shifted from that of the Lorentzian.

Fig. 7 .
Fig. 7.The generic computation domain for the numerical implementation of the formulation.PMLs are inserted in the inner sides of the computation domain while the outer boundaries of the computation domain are set to PECs or PMCs.

Table 1 .
The Comparison of hω n , Q n , and g th,n between the Theoretical and FP Estimations