Approach to calculate normal modes by decomposing the dyadic Green’s function

Wenhai Yu; Wencheng Yue; Peijun Yao; Yonghua Lu; Wen Liu

doi:10.1364/OE.22.026712

1. Introduction

Modes are a fundamental ubiquitous concept in classical optics, for example laser physics. Modes are also popular in quantum optics because they simplify the quantum description of light. Electromagnetic processes in a perfect lossless isolated system could be described, in a mathematically precise way, in terms of normal modes, however, from a fundamental point of view, all practical systems are inevitably leaky. In much of quantum optics, for example the early days of the laser, it is assumed that the normal cavity modes remain and the damping is modeled by coupling these normal modes to the normal modes of a reservoir [1]. However, this is a good approximation only when the cavity quality factor Q is high. To deal with general leaky systems, an approach based on quasi-normal modes is frequently used, the quasi-normal modes provide an intuitively appealing description of the allowed electromagnetic vibrations in the leaky system. Fox and Li showed that, in the case of a Fabry-Perot-type cavity with perfectly reflecting mirrors, there exists a discrete set of quasi-normal modes for which the diffraction leakage from the cavity is small [2]. Ho et al. formulated the second quantization of a scalar field in an open cavity in terms of the quasi-normal modes, single-resonance domination of the density of states (DOS) and the spontaneous decay rate is then given a proper foundation [3]. Leung et al. showed that the quasi-normal modes of a leaky cavity form a complete set inside the cavity, the quasi-normal modes are also orthogonal under a modified definition of the inner product. The completeness and orthogonality hold even though the cavity is not a Hermitian system by itself [4, 5].

Alternatively, one might formulate the electromagnetic processes in terms of normal modes of the closed system formed by the leaky system and the rest of the universe. Here normal modes are defined in the usual sense: eigensolutions to the source-free wave equation with real eigenfrequencies. The boundary condition for the normal modes calculated by our method is that the fields decay to zero when |r| → ∞. One might argue that normal modes are not allowed for a photonic dielectric structure because of leaky properties, such as a dielectric cavity. In fact, the so-called “leakage” is relative to the local region of interest not for the whole space or the universe, for the universe, the system is still Hermitian [4, 6] and the eigenmodes could be normal modes [7, 8], yet the normal modes are continuous because of the infinity of universe. Much work of quantum optics is based on normal modes, for example, an exact quantum description of the field in a leaky system has been developed using the normal modes [7, 9–11]; Wubs et al. quantized the electromagnetic field by expanding the vector potential and its canonically conjugate field in normal modes, and investigated multipole interaction between atoms and their photonic environment [12]. Although much theoretical work in optics is described by normal modes, there is less work on how to calculate normal modes for a leaky system. Except for some typical structures, one generally use numerical methods with outgoing wave boundary conditions to calculate the eigenmodes, which are actually the quasi-normal modes [13], for example, frequency-domain methods (e.g. FDFD [14]) or time-domain method(e.g. FDTD) [15]. One of the main characteristics of quasi-normal modes is the exponential divergence in the far field as the imaginary part of the corresponding complex wave number is always less than zero for stable optical systems.

In this paper, we report a way to calculate the normal modes of an arbitrary photonic dielectric structure by decomposing Green’s function tensor. The essence of our method is to excite the corresponding structure with one or several dipoles and get the imaginary part of the excited field, the Purcell factors calculated through the normal modes match very well with the direct calculations by FDTD. By our method, the modes can be directly normalized without additional integration, the degeneracy could be removed by simply placing a few dipolar sources according to the symmetry of the calculated structure. As examples, calculations are provided for the normal modes of a planar photonic crystal cavity.

2. Derivation

Our derivation starts with the electric Green’s function tensor G(r, r′, ω) which is defined as follows,

[\nabla \times \nabla \times - \frac{ω^{2}}{c^{2}} ε (r)] G (r, r^{'}, ω) = \frac{ω^{2}}{c^{2}} 1 δ (r - r^{'}),

together with the Sommerfeld radiation condition, here, ε(r) is the relative permittivity, ω the angular frequency, c the light speed in vacuum and 1 denotes the unit dyad. We assume the system is non-lossy and non-dispersive, i.e., ε(r) is real. Although the operator

ℋ = \frac{1}{\sqrt{ε (r)}} \nabla \times \nabla \times \frac{1}{\sqrt{ε (r)}}

is not Hermitian for the cavity itself because of the leaky property, it is still Hermitian in terms of the whole universe [4]. Then due to the Hermitian property of operator ℋ in the universe, Green’s function tensor can be written in terms of eigenmodes of the system,

\begin{matrix} G (r, r^{'}, ω) = \sum_{n} \int_{0}^{\infty} d ν \frac{ω^{2}}{(ν - ω - i δ) (ν + ω + i δ)} \\ \times f_{n}^{T} (ν, r) {[f_{n}^{T} (ν, r^{'})]}^{*} - \sum_{λ} f_{λ}^{L} (r) {[f_{λ}^{L} (r^{'})]}^{*}, \end{matrix}

where δ is a positive infinitesimal that assures the causality,

f_{n}^{T} (ν, r)

represents the continuous transverse modes with real eigenfrequency ν,

f_{λ}^{L} (r)

the longitudinal modes with zero eigenfrequency, and the degenerate modes are characterized by n. The completeness of the eigenfunction set is

\sum_{n} \int d ν ε (r) f_{n}^{T} (ν, r) {[f_{n}^{T} (ν, r^{'})]}^{*} + \sum_{λ} ε (r) f_{λ}^{L} (r) {[f_{λ}^{L} (r^{'})]}^{*} = 1 δ (r - r^{'})

, the normalization condition is

\int d r ε (r) f_{n}^{T} (ν, r) f_{n}^{T} (ν^{'}, r) = δ (ν - ν^{'})

. Generally,

f_{n}^{T} (ν, r)

is complex, however, without loss of generality, here we assume it is real. If

f_{n}^{T} (ν, r)

is complex, the real normal modes can be constructed through

{f_{n}^{T} (ν, r) + {[f_{n}^{T} (ν, r)]}^{*}} / 2

and

{f_{n}^{T} (ν, r) - {[f_{n}^{T} (ν, r)]}^{*}} / 2 i

. It is worth mentioning that the choice of real modes will simplify the quantization procedure [12]. The Green’s function tensor for a photonic dielectric structure can also be expanded by discrete quasi-normal modes [5].

Extending the integral interval in Eq. (2) from (0, ∞) to (−∞, 0)∪(0, ∞), then

\begin{array}{r} G (r, r^{'}, ω) = \frac{1}{2} \sum_{n} (\int_{- \infty}^{0} + \int_{0}^{\infty}) d ν \frac{ω^{2}}{(ν - ω - i δ) (ν + ω + i δ)} \\ \times f_{n}^{T} (ν, r) {[f_{n}^{T} (ν, r^{'})]}^{*} - \sum_{λ} f_{λ}^{L} (r) {[f_{λ}^{L} (r^{'})]}^{*}, \end{array}

In view of the analytical properties of

f_{n}^{T} (ν, r)

in the space of frequency domain (ν = 0 is excluded), using the residue theorem in the upper or lower half space of frequency domain, we have

\begin{array}{l} G (r, r^{'}, ω) & = i \frac{π ω}{2} \sum_{n} f_{n}^{T} (ω, r) f_{n}^{T} (ω, r^{'}) - \sum_{λ} f_{λ}^{L} (r) {[f_{λ}^{L} (r^{'})]}^{*} \\ - \frac{1}{2} \sum_{n} \int_{C_{\inf} + C_{0}} d ν \frac{ω^{2} f_{n}^{T} (ν, r) {[f_{n}^{T} (ν, r^{'})]}^{*}}{(ν - ω - i δ) (ν + ω + i δ)}, \end{array}

where C_inf (C₀) denotes the integral path of infinite (infinitesimal) half cycle with center 0 either in upper or lower space of frequency domain. Noticing that the second and third terms on the right hand side of Eq. (4) are real (in fact, the second term is corresponding to longitudinal modes and has analytical expression [16]), and taking the imaginary part of the above equation, we have

Im {G (r, r^{'}, ω)} = \frac{π ω}{2} \sum_{n} f_{n}^{T} (ω, r) f_{n}^{T} (ω, r^{'}),

If several dipole sources with the real dipole moment are added into the system, then the imaginary part of the excited field is

Im {E (r, ω)} = \frac{π ω}{2 ε_{0}} \sum_{n, m} f_{n}^{T} (ω, r) f_{n}^{T} (ω, r_{m}) \cdot p_{m} (ω)

with p_m the real moment of the m-th dipole. Only one dipole is needed to obtain the normal mode if the system is single-mode or under single-mode approximation, and several dipoles, placed according to symmetry, are necessary to extract the normal modes if the system is degenerate. By our method, the normal mode can be directly normalized, for example, if the system is single-mode or under single-mode approximation, then the normalized mode can be obtained as

f^{T} (ω, r) = \frac{Im {E (r, ω)}}{\sqrt{π ω p_{0} (ω) \cdot Im {E (r_{0}, ω)} / 2 ε_{0}}},

where p₀(ω) and r₀ are the dipole moment and the position of dipole source respectively. For the degenerate case, we can place symmetrically several suitable dipole sources to excite the mode of interest and cancel the others, then the normalized normal mode can be calculated from Eq. (6) similar to the single-mode case. We emphasize that the mode obtained by Eq. (7) is not the scattered field of a dipole: (i) the scattered field of a dipole is complex, but our mode is real; (ii) the scattered field of a dipole cannot satisfy source-free wave equation, but our mode does, as it can be easily validated by taking the imaginary part on Eq. (1) and using the fact of E(r, ω) = G(r, r₀, ω) · p₀(ω)/ε₀. According to our derivation, it is not only a normal mode but also a normalized normal mode. The normal modes exist at any frequency with the leaky system, i.e., we can get normal modes for any specified frequency. As is known, there are discrete resonant normal modes for a closed lossless system, and the resonant normal mode is corresponding to the infinite peak of density of states, however, when the system becomes increasingly leaky, the resonant peak will broaden. In view of this point, we associate the resonant normal mode (or modes if there exists degeneracy for the resonance) in our work with the peak of DOS. Hereafter, the normal mode will refer to resonant normal mode unless it is expressly stated.

3. Example

As an example, we calculate the normal mode of photonic microcavity in a photonic crystal slab, the cavity illustrated in Fig. 1(a) is formed by removing three air holes in the center from a PC slab, the background is air, its thickness is 0.5a, the lattice constant a is 420nm, the radius of air holes is 0.275a, its relative dielectric constant is 12, the slab’s length is 7.5a and width 5a. The highlighted holes with the red circles are shifted a/40 outwards along x from the nominal lattice position. A single resonant mode exists at the eigenfrequency of 1.221 × 10¹⁵rad/s, which is non-degenerate according to the analysis based on group theory. The single mode will be excited as long as only one dipole is randomly placed except for some special position of vanishing electric fields. The modes can be excited by either harmonic or pulsed dipoles. In the FDTD simulation, the dipole source is located at the center of the PC slab and parallel to y axis, and we use a Gaussian pulse to excite the mode, the pulse length is 9 fs, the center frequency is 1.221 × 10¹⁵rad/s, the field’s distribution of E(r, ω) at 1.221 × 10¹⁵rad/s is obtained through Fourier integral where we integrate the field from the first time step, then the normalized normal mode is obtained through Eq. (7). This way is easier than quasi-mode method which integrates after the source field disappears, especially for low-Q modes, because it is hard to choose a proper start time. The distribution of y-component of normal mode $f_{y}^{T}$ is shown in Fig. 1(b).

Fig. 1. The sketch of photonic crystal slab cavity. The highlighted holes with the red circles are shifted outwards along x from the nominal lattice position. Right: the spatial distribution of y-component of the cavity’s normal mode in the plane z = 0.

Download Full Size | PDF

Through the previously obtained normal mode, the distribution of Purcell factor is calculated under the single-mode approximation. We first calculate the yy component of the imaginary part of Green’s function tensor through Eq. (5), i.e., Im{G_yy(r, r′ = r, ω)}, the Purcell factor corresponding y-polarized dipole is thus obtained by $F_{y} = Im {G_{y y} (r, r, ω)} / Im {G_{y y}^{hom} (r, r^{'} = r ω)}$ , where G^hom(r, r′, ω) indicates Green’s function tensor for homogeneous medium. The relative dielectric constant of homogeneous medium is 12 in our calculation. The results are shown in Fig. 2, the blue solid curve shows the distribution of Purcell factor along x axis, while the red circles indicate the Purcell factors calculated directly by FDTD. We can see that the Purcell factors calculated through the normal mode match very well with the direct calculations by FDTD. It’s noted that we need to perform a calculation for each frequency of interest, fortunately, usually only the normal mode at the resonance frequency is needed to characterize a cavity, for example, the mode volume actually corresponds to the resonant peak.

Fig. 2. Purcell factor corresponding to y-polarized dipole as a function of position x. Purcell factor indicated by the blue solid curve is obtained through the normal mode, while the red circles corresponds to the direct calculations by FDTD.

Download Full Size | PDF

To prove that the modes calculated by our method are not quasi-normal modes, the distributions of the normal mode and the quasi-normal mode are calculated. Considering that the quasi-normal mode for the resonant structure with a low quality factor will diverge in a small domain and thus bring us the convenience of numerical calculation, a silicon cuboid encompassed with air is chosen in our simulation, the cuboid is of length=500nm, width=200nm, and height=150nm and possesses a resonance at the frequency of 2.832 × 10¹⁵rad/s. The dipole source is located at the center of the cuboid and parallel to y axis. We use a Gaussian pulse to excite the mode, the pulse length is 4 fs, the center frequency is 2.832 × 10¹⁵rad/s. The distribution of the y-component of the normalized normal mode along x axis is plotted in Fig. 3(a), it can be seen that the field oscillates and approaches to zero when the position moves away from the cavity, in fact, the imaginary part of far field Im{E(r, ω)} ∝ sin(k₀r + ϕ₀)/r with k₀ = ω/c and constant phase ϕ₀, which is coincident with the fitting curve in Fig. 3(a). We note that our normal modes don’t satisfy the Sommerfeld radiation condition any longer although Green’s function tensor does. For comparison, the corresponding quasi-normal mode Ẽ(r) is also calculated by FDTD with the initial field filtered, the quasi-normal mode is not normalized and its y-component Ẽ_y versus x is shown in Fig. 3(b). As we can see from Fig. 3(b), the field of quasi-normal mode is already exponentially divergent within the same calculation region as normal mode.

Fig. 3. (a) The y-component of the normalized normal mode of the silicon cuboid as a function of position x. Blue solid line shows the normalized mode. Red solid line in the inset shows the fitting curve. The fitting function is 10.974/(x − 0.364), the shift on x could be neglected at very far field; (b) the y-component of the quasi-normal mode of the silicon cuboid as a function of position x

Download Full Size | PDF

Our method can easily remove the degeneracy. According to group theory, the above-stated structure belongs to D_2h and there are in total eight irreducible representations, i.e., A_g,B_1g,B_2g,B_3g,A_u,B_1u,B_2u,B_3u. For simplicity, we assume that the symmetry is even with respect to the plane z = 0, which is achieved by restricting all the dipoles in the symmetric plane, in other words, D_2h is reduced to C_2v, thus there are four possible normal modes for a given resonance frequency. The modes can then be extracted by placing four suitable dipole sources. In fact, we find it is only two-fold degenerate at the eigenfrequency of 1.305×10¹⁵rad/s, they are corresponding to B_3u, B_1g respectively. In our FDTD simulation, four dipole sources with Gaussian pulse excitation are placed at (±100, ±50, 0) nm and their dipole moments are parallel to y axis. The pulse length of the excitation source is 9 fs, the center frequency is 1.305×10¹⁵rad/s. When the orientations of the dipoles in the first and the third quadrants are opposite to that in the second and the fourth (the dipoles with opposite orientation are of phase difference π), the mode of B_3u is excited and the others are cancelled. Thus the field’s distribution of E(r, ω) at 1.305 × 10¹⁵rad/s is obtained through Fourier integral, and the normalized normal mode of B_3u can be calculated from Eq. (6) similar to the single-mode case, here we will use the fact that $f_{B_{3 u}}^{T} (ω, r_{m}) \cdot p_{m} (ω)$ is equal to each other for different m because of the symmetry. If the orientation of the dipoles in the first and the fourth quadrants are opposite to that in the second and the third, only the mode of B_1g is excited. The normal modes’ distributions of electric field are shown in Fig. 4. The degeneracy of D_2h point group can also be removed by imposing the symmetrical boundary condition [17], but it will be very hard for more complex point group, for example, C_3v. However it does not challenge our method, the degeneracy can be similarly removed by placing 6 dipoles with suitable real amplitude according to symmetry.

Fig. 4. The spatial distribution of x- and y-component of mode B_3u and B_1g in the plane z = 0, the resonance frequency is 1.305 × 10¹⁵rad/s. (a) $f_{x}^{T}$ for mode B_3u, (b) $f_{y}^{T}$ for mode B_3u, (c) $f_{x}^{T}$ for mode B_1g, (d) $f_{y}^{T}$ for mode B_1g.

Download Full Size | PDF

4. Conclusion

In summary, we have presented a method to calculate normal modes for arbitrary photonic dielectric structure. By our method, the modes can be directly normalized without additional integration. As an example, a single normal mode of planar photonic crystal cavity is calculated by FDTD. If the normal modes are degenerate, we can extract the different degenerate modes when several dipoles are placed according to symmetry. With the same dielectric structure, we demonstrate how to remove mode degeneracy according to group theory. Note that our calculation is only suitable to non-lossy and non-dispersive systems, and the extension to lossy materials will merit further research.

Acknowledgments

The authors would like to gratefully acknowledge the helpful discussions with S. Hughes of Department of Physics, Queen’s University, Canada. This work is supported by National Key Basic Research Program of China No. 2012CB922003, and National Natural Science Foundation of China No. 61177053, and Fundamental Research Funds for the Central Universities No. WK2470000002.

References and links

1. C. W. Gardiner and M. J. Collett, “Input and output in damped quantum systems: Quantum stochastic differential equations and the master equation,” Phys. Rev. A 31, 3761–3774 (1985). [CrossRef] [PubMed]

2. A. Fox and T. Li, “Resonant modes in a maser interferometer,” Bell Syst. Tech. J 40, 453–488 (1961). [CrossRef]

3. K. C. Ho, P. T. Leung, A. Maassen van den Brink, and K. Young, “Second quantization of open systems using quasinormal modes,” Phys. Rev. E 58, 2965–2978 (1998). [CrossRef]

4. P. T. Leung, S. Y. Liu, and K. Young, “Completeness and orthogonality of quasinormal modes in leaky optical cavities,” Phys. Rev. A 49, 3057–3067 (1994). [CrossRef] [PubMed]

5. K. M. Lee, P. T. Leung, and K. M. Pang, “Dyadic formulation of morphology-dependent resonances. i. completeness relation,” J. Opt. Soc. Am. B 16, 1409–1417 (1999). [CrossRef]

6. J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light (Princeton University Press, 2011).

7. R. Lang, M. O. Scully, and W. E. Lamb Jr, “Why is the laser line so narrow? a theory of single-quasimode laser operation,” Physical Review A 7, 1788 (1973). [CrossRef]

8. J. C. Penaforte and B. Baseia, “Quantum theory of a one-dimensional laser with output coupling: Linear approximation,” Phys. Rev. A 30, 1401–1406 (1984). [CrossRef]

9. R. J. Glauber and M. Lewenstein, “Quantum optics of dielectric media,” Phys. Rev. A 43, 467–491 (1991). [CrossRef] [PubMed]

10. S. M. Dutra and P. L. Knight, “Spontaneous emission in a planar fabry-pérot microcavity,” Phys. Rev. A 53, 3587–3605 (1996). [CrossRef] [PubMed]

11. K. Ujihara, “Quantum theory of a one-dimensional optical cavity with output coupling. field quantization,” Phys. Rev. A 12, 148–158 (1975). [CrossRef]

12. M. Wubs, L. G. Suttorp, and A. Lagendijk, “Multipole interaction between atoms and their photonic environment,” Phys. Rev. A 68, 013822 (2003). [CrossRef]

13. P. T. Kristensen and S. Hughes, “Modes and mode volumes of leaky optical cavities and plasmonic nanoresonators,” ACS Photonics 1, 2–10 (2014). [CrossRef]

14. C.-P. Yu and H.-C. Chang, “Yee-mesh-based finite difference eigenmode solver with pml absorbing boundary conditions for optical waveguides and photonic crystal fibers,” Optics Express 12, 6165–6177 (2004). [CrossRef] [PubMed]

15. A. Taflove, Computational Electrodynamics: The Finite - Difference Time - Domain Method (Artech House, Incorporated, 1995).

16. M. Wubs, L. G. Suttorp, and A. Lagendijk, “Multiple-scattering approach to interatomic interactions and super-radiance in inhomogeneous dielectrics,” Phys. Rev. A 70, 053823 (2004). [CrossRef]

17. Y. Tang, A. Mintairov, J. Merz, V. Tokranov, and S. Oktyabrsky, “Characterization of 2d-photonic crystal nanocavities by polarization-dependent photoluminescence,” in “Nanotechnology, 2005. 5th IEEE Conference on,” (IEEE, 2005), pp. 35–38.

Approach to calculate normal modes by decomposing the dyadic Green’s function

Abstract

1. Introduction

2. Derivation

3. Example

4. Conclusion

Acknowledgments

References and links

Cited By

Figures (4)

Equations (7)

Optics Express