Frequency-domain formulation of photonic crystals using sources and gain

Po-Jui Chiang; Shu-Wei Chang

doi:10.1364/OE.21.001972

1. Introduction

Photonic periodic structures [1–4] have been widely utilized in optical fibers [5], laser cavities [6–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:

\nabla \times \nabla \times E (r) = {(\frac{ω_{b}}{c})}^{2} \bar{\bar{ε}} (r) E (r), E (r) = e^{i k \cdot r} F (r),

{\bar{\bar{ε}}}_{r} (r + R_{m}) = {\bar{\bar{ε}}}_{r} (r); F (r + R_{m}) = F (r), m = 1 - 3,

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)² plays the role of eigenvalues, and the convention of exp(−iω_bt) 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–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–19]. The frequency dispersion, nevertheless, turns the GE problem in Eq. (1a) into a nonlinear one. A few authors have proposed specific techniques [20–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 low-order 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–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, white-source 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.

2. 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 ω:

\nabla \times \nabla \times E_{k} (r) - {(\frac{ω}{c})}^{2} {\bar{\bar{ε}}}_{r} (r, ω) E_{k} (r) = i ω μ_{0} J_{s, k} (r) = 0, \forall r \notin Ω_{a},

J_{s, k} (r + R_{m}) = e^{i k \cdot R_{m}} J_{s, k} (r), m = 1 - 3.

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_n_k(r, ω)} and sources {j_s,_n_k(r, ω)}, respectively, under the constraint of the phase BC (n is the label of the basis). For simplicity, we demand that j_s,_n_k(r, ω) generates f_n_k(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_n_k(ω)} as

E_{k} (r) \approx \sum_{n} c_{n k} (ω) f_{n k} (r, ω),

J_{s, k} (r) \approx \sum_{n} c_{n k} (ω) j_{s, n k} (r, ω) .

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.

Download Full Size | PDF

A sensible construction of j_s,nk(r, ω) is to set it proportional to f_n_k(r, ω) in Ω_a but force it to vanish elsewhere [38], namely,

j_{s, n k} (r, ω) = - i ω ε_{0} Δ ε_{r, n k} (ω) U (r) f_{n k} (r, ω),

U (r + R_{m}) = U (r), m = 1 - 3,

where U(r) is also a periodic function in repeated unit cells and is unity in Ω_a but zero in Ω_b; and Δε_r,_n_k(ω) 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,_n_k(r, ω)}. The substitution of the ansatz in Eq. (4a) into the wave equation of f_n_k(r, ω) leads to a GE problem as follows:

\nabla \times \nabla \times f_{n k} (r, ω) - {(\frac{ω}{c})}^{2} {\bar{\bar{ε}}}_{r} (r, ω) f_{n k} (r, ω) = {(\frac{ω}{c})}^{2} Δ ε_{r, n k} (ω) U (r) f_{n k} (r, ω),

f_{n k} (r + R_{m}, ω) = e^{i k \cdot R_{m}} f_{n k} (r, ω), m = 1 - 3.

In Eq. (5a), the factor (ω/c)²Δε_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_n_k(r, ω), it is applied to j_s,_n_k(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_n_k(r, ω) as

f_{n k} (r, ω) = e^{i k \cdot r} ψ_{n k} (r, ω), ψ_{n k} (r + R_{m}, ω) = ψ_{n k} (r, ω), m = 1 - 3,

where ψ_n_k(r) is the Bloch periodic part of f_n_k(r, ω), Similarly, we write the source j_s,_n_k(r, ω) in the Bloch form as

j_{s, n k} (r, ω) = e^{i k \cdot r} ς_{s, n k} (r, ω), ς_{s, n k} (r + R_{m}, ω) = ς_{s, n k} (r, ω), m = 1 - 3,

ς_{s, n k} (r, ω) = - i ω ε_{0} Δ ε_{r, n k} (ω) U (r) ψ_{n k} (r, ω),

where ς_s,_n_k(r, ω) is the Bloch periodic part of the source j_s,_n_k(r, ω).

We note that the basis set {j_s,_n_k(r, ω)} may not exhaust all the possibilities of the current source J_s,_k(r). The basis current density in {j_s,_n_k(r, ω)} could be approximately transverse in Ω_a [∇ · j_s,_n_k(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,_n_k(r, ω)}.

Extracting the amplitude c_n_k(ω) of a basis source j_s,_n_k(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 ψ_n_k(r, ω) and ς_s,_n_k(r, ω) for reciprocal PhCs. This biorthogonality condition also serves as a tool to uniquely determine the expansion coefficient c_n_k(ω). 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_n_k(ω), (resonance) frequency ω_n_k, and quality factor Q_n_k of mode (n, k). We first consider a source J_s,_k(r) that is proportional to j_s,_n_k(r, ω) and demand the spatial integration of |J_s,_k(r)|² in Ω_a to be frequency-independent (white source):

J_{s, k} (r) = a (ω) j_{s, n k} (r, ω),

𝒥^{2} V_{a} = \int_{Ω_{a}} d r {| J_{s, k} (r) |}^{2} = {| a {(ω)}^{2} ε_{0}^{2} ω Δ ε_{r, n k} (ω) |}^{2} \int_{Ω_{a}} d r {| f_{n k} (r, ω) |}^{2},

where a(ω) is a frequency-dependent amplitude; 𝒥 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(ω)|² can be obtained through the white-source condition in Eq. (8b). After the substitution of this expressions for |a(ω)|² in the generated power, the white-source power spectrum P_n_k(ω) is connected to the parameter Δε_n_k(ω) as follows:

\begin{array}{l} P_{n k} (ω) & = - \frac{1}{2} \int_{Ω_{a}} d r Re [J_{s, k}^{*} (r) \cdot E_{k} (r)] = - \frac{{| a (ω) |}^{2}}{2} \int_{Ω_{a}} d r Re [j_{s, n k}^{*} (r, ω) \cdot f_{n k} (r, ω)] \\ = \frac{𝒥^{2} V_{a}}{2 ε_{0}} Im [\frac{1}{ω Δ ε_{r, n k} (ω)}] . \end{array}

The (resonance) frequency ω_n_k is the one which minimizes the absolute value |ωΔε_r,_n_k(ω)| of the denominator that appears inside the imaginary part of Eq. (9). The dispersion curve of ω_n_k 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 ω_n_k, the field f_n_k(r, ω_n_k) is considered as (at least approximately) the Bloch mode of band n at k. Also, if one sets ω = ω_n_k in Eq. (5a), the imaginary part Im[Δε_r,nk(ω_n_k)] 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,_n_k(ω_n_k)] automatically.

The magnitude P_n_k(ω_n_k) at ω = ω_n_k is usually close to the peak value if the lineshape around ω_n_k is sufficiently narrow. Under these circumstances, we can approximate P_n_k(ω_n_k) around ω_n_k with a Lorentzian. The ratio between the resonance frequency ω_n_k 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,_n_k(ω) around ω = ω_n_k[38] and is written as

Q_{n k} = \frac{i}{2 Δ ε_{r, n k} (ω_{n k})} {\frac{\partial [ω Δ ε_{r, n k} (ω)]}{\partial ω} |}_{ω = ω_{n k}} .

3. 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₁ (d₂). The unit cell has a size a = d₁ + d₂. Various combinations of materials will be considered.

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

Download Full Size | PDF

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_xx̂ + 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_zẑ 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 white-source power spectra P_n_k(ω) of these modes will be depicted on the k_xa/π − ωa/(2πc) plane for easy visualizations. The resonance frequencies ω_n_k obtained with minimum estimations of |ωΔε_r,_n_k(ω)| are also plotted as a function of k_xa/π on the same plane. Both ω_n_k and Q_n_k will be compared with those calculated from the complex-ω method.

3.1. 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:

\frac{\partial^{2} ϕ_{n k} (x)}{\partial x^{2}} + {{(\frac{ω}{c})}^{2} [ε_{a} (ω) + Δ ε_{r, n k} (ω)] - k_{y}^{2}} ϕ_{n k} (x) = 0, x \in Ω_{a},

\frac{\partial^{2} ϕ_{n k} (x)}{\partial x^{2}} + [{(\frac{ω}{c})}^{2} ε_{b} (ω) - k_{y}^{2}] ϕ_{n k} (x) = 0, x \in Ω_{b},

where ϕ_n_k(x), aside to a factor exp(ik_yy), 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 ϕ_n_k(x) can be separately solved in Ω_a and Ω_b as

ϕ_{n k} (x) = {\begin{array}{l} A e^{i M x} + B e^{- i M x}, & x \in Ω_{a}, \\ C e^{i N x} + D e^{- i N x}, & x \in Ω_{b}, \end{array} {\begin{array}{l} M^{2} = {(ω / c)}^{2} [ε_{a} (ω) + Δ ε_{r, n k} (ω)] - k_{y}^{2}, \\ N^{2} = {(ω / c)}^{2} ε_{b} (ω) - k_{y}^{2}, \end{array}

where A, B, C, and D are amplitudes to be determined. With Eq. (12), we demand the continuity of ϕ_n_k(x) at x = 0 and impose the phase BC at x = −d₁ and d₂ with a phase factor exp(ik_xa). 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:

(\begin{array}{l} 0 \\ 0 \\ 0 \\ 0 \end{array}) = (\begin{matrix} 1 & 1 & - 1 & - 1 \\ \frac{M}{N} & - \frac{M}{N} & - τ & τ \\ e^{i (k_{x} a - M d_{1})} & e^{i (k_{x} a + M d_{1})} & - e^{i N d_{2}} & - e^{- i N d_{2}} \\ \frac{M}{N} e^{i (k_{x} a - M d_{1})} & - \frac{M}{N} e^{i (k_{x} a + M d_{1})} & - τ e^{i N d_{2}} & τ e^{- i N d_{2}} \end{matrix}) (\begin{matrix} A \\ B \\ C \\ D \end{matrix}),

τ = {\begin{matrix} [ε_{a} (ω) + Δ ε_{r, n k} (ω)] / ε_{b} (ω) & (TM), \\ 1 & (TE) . \end{matrix}

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]:

cos (k_{x} a) = cos (M d_{1}) cos (N d_{2}) + \frac{1}{2} (\frac{M}{τ N} + \frac{τ N}{M}) sin (M d_{1}) sin (N d_{2}) .

For given ω, k_x, and k_y, parameters Δε_r,_n_k(ω) are solved self-consistently using Eq. (14). These parameters are then processed to extract the information of interest.

3.2. 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₁ and d₂ 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_xa/π − ω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,_n_k(ω) 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 ω_n_k 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,_n_k(ω_n_k) = 0, namely, the magnitudes |ωΔε_r,_n_k(ω)| have minima equal to zero at resonances. In addition, the parameter Δε_r,_n_k(ω) is real when ω is off resonances because the power P_n_k(ω) in Eq. (9) does not need to balance any loss.

Figure 3(a) shows |ωΔε_r,_n₀(ω)| 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),₀ and ω_(TE,3),₀ are positions of the respective minima. For simplicity, at each frequency ω, we only show the smallest among magnitudes |ωΔε_r,_n₀(ω)| 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,_n₀(ω)| abruptly changes due to mode switching. As expected, the magnitude |ωΔε_r,_n₀(ω)| 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,_n_k(ω)| in Fig. 3(b). For clear illustrations, these curves are depicted in a 3D manner so that magnitudes |ωΔε_r,_n_k(ω)| increase as the curves bend into the figure. The loci of minima (zeros) on each curve are projected onto the k_xa/π − ω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 account in this case. Again, the two sets of dispersion curves from the two schemes agree with each other.

Fig. 3 (a) The behaviors of |ωΔε_r,_n₀(ω)| 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,_n_k(ω)| at k_y = 0 (sold lines) and counterparts from the complex-ω method (red circles). For clear illustrations, the magnitudes |ωΔε_r,_n_k(ω)| increase as the curves bend into the figure.

Download Full Size | PDF

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.

Download Full Size | PDF

In lossless structures, we can alternatively imagine that the parameter Δε_r,_n_k(ω) acquires an infinitesimal imaginary part at the resonance. In this way, the power spectrum P_n_k(ω) in Eq. (9) is proportional to a delta function centered at ω_n_k, which is the sign of a persistent mode with an infinitely long lifetime. Accordingly, the Q factor Q_n_k in Eq. (10) also approaches infinity.

3.3. 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₁ and d₂ 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_xa/π ≲ 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_xa/π. 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).

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

Download Full Size | PDF

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 Δε_n_k(ω_n_k). At the resonance, the parameter Δε_n_k(ω_n_k) is nearly imaginary, as can be observed from Table 1. In this case, Δε_n_k(ω_n_k) reflects the amount of gain to be inserted into Ω_a so that the metal loss is compensated, and the mode (n, k) can start to self oscillate (lase) at the real frequency ω_n_k. 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.

Table 1. Comparisons between resonance frequencies as well as Q factors, and a list of Δε_r,_n₀(ω_n₀) for the four lowest TE modes at the BZ center in Fig. 5.

View Table

3.4. 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 - ω_{p}^{2} / (ω^{2} + i γ ω)$ for the dispersive metal, where the plasma frequency ω_p is 10¹⁵ rad·s⁻¹, and the damping term γ is set to 0.01ω_p. The width d₁ and d₂ 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_n_k(ω) of the three lowest-order TE modes, respectively. The calculations of the complex mode frequencies are implemented with the transfer-matrix method [21]. The mode frequencies ω_n_k 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_n_k(ω) in Fig. 6(b).

Fig. 6 (a) The band diagrams and (b) quasi-3D views of P_n_k(ω) 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.

Download Full Size | PDF

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} / \sqrt{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 BZ center to boundary. As a result, the lineshapes of corresponding power spectra P_n_k(ω) 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_n_k(ω) in Fig. 7(b) becomes sharper as k_x comes closer to the BZ boundary, indicating a high Q factor of the branch there.

Fig. 7 (a) The band diagrams for the three lowest-order TM modes and (b) lateral 3D views of P_n_k(ω) 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.

Download Full Size | PDF

4. 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–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_n_k(r, ω) from f_n_k(r, ω) as follows:

g_{n k} (r, ω) = \frac{1}{i ω μ_{0}} \nabla \times f_{n k} (r, ω) \equiv e^{i k \cdot r} φ_{n k} (r, ω),

φ_{n k} (r + R_{m}, ω) = φ_{n k} (r), m = 1 - 3,

where φ_n_k(r, ω) is the Bloch periodic part of the magnetic field g_n_k(r, ω). In terms of the Bloch periodic parts ψ_n_k(r, ω), φ_n_k(r, ω), and ς_s,_n_k(r, ω), Maxwell’s equations are rewritten as

\nabla \times ψ_{n k} (r, ω) + i k \times ψ_{n k} (r, ω) = i ω μ_{0} φ_{n k} (r, ω),

\nabla \times φ_{n k} (r, ω) + i k \times φ_{n k} (r, ω) = i ω ε_{0} {\bar{\bar{ε}}}_{r} (r, ω) ψ_{n k} (r, ω) + ς_{s, n k} (r, ω) .

In reciprocal PhCs, the relative permittivity tensor is symmetric when represented in real and orthonormal coordinate bases, namely, ε̳_r(r, ω) = ε̳_r (r, ω_T). 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

\begin{array}{l} \oint_{S_{uc}} d a \cdot [ψ_{n k} (r, ω) \times φ_{n^{'} k^{'}} (r, ω) - ψ_{n^{'} k^{'}} (r, ω) \times φ_{n k} (r, ω)] \\ = - \int_{Ω_{a}} d r i (k + k^{'}) \cdot [ψ_{n k} (r, ω) \times φ_{n^{'} k^{'}} (r, ω) - ψ_{n^{'} k^{'}} (r, ω) \times φ_{n k} (r, ω)] \\ - \int_{Ω_{uc}} d r [ψ_{n k} (r, ω) \cdot ς_{s, n^{'} k^{'}} (r, ω) - ψ_{n^{'} k^{'}} (r, ω) \cdot ς_{s, n k} (r, ω)] . \end{array}

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 ψ_n_k(r, ω), ς_s,_n_k(r, ω), ψ_n_′−_k(r, ω), and ς_s,_n_′−_k(r, ω):

\begin{array}{l} 0 = \int_{Ω_{uc}} d r [ψ_{n k} (r, ω) \cdot ς_{s, n^{'} - k} (r, ω) - ψ_{n^{'} - k} (r, ω) \cdot ς_{s, n k} (r, ω)] \\ = - i ω ε_{0} [Δ ε_{r, n^{'} - k} (ω) - Δ ε_{r, n k} (ω)] \int_{Ω_{a}} d r ψ_{n k} (r, ω) \cdot ψ_{n^{'} - k} (r, ω) . \end{array}

From the second line of Eq. (18), if Δε_r,_n_′−_k(ω) ≠ Δε_r,_n_k(ω), the volume integral of the dot product ψ_n_k(r, ω)· ψ_n_′−_k(r, ω) in Ω_a has to be zero. The integral can be nonvanishing if the two parameters Δε_r,_n_′−_k(ω) and Δε_r,_n_k(ω) 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:

\int_{Ω_{a}} d r ψ_{n k} (r, ω) \cdot ψ_{n^{'} - k} (r, ω) = \int_{Ω_{a}} d r f_{n k} (r, ω) \cdot f_{n^{'} - k} (r, ω) \equiv δ_{n n^{'}} Λ_{n k} (ω),

\int_{Ω_{a}} d r ς_{s, n k} (r, ω) \cdot ς_{s, n^{'} - k} (r, ω) = \int_{Ω_{a}} d r j_{s, n k} (r, ω) \cdot j_{s, n^{'} - k} (r, ω) \equiv δ_{n n^{'}} Θ_{n k} (r, ω),

Θ_{n k} (ω) = - {[ω ε_{0} Δ ε_{r, n k} (ω)]}^{2} Λ_{n k} (ω),

where δ_nn_′ is the Kroneckers delta; Λ_n_k(ω) = Λ_n₋_k(ω) is the complex normalization constant for f_n_k(r, ω) and f_n₋_k(r, ω); and Θ_n_k(ω) = Θ_n₋_k(ω) is the counterpart for j_s,_n_k(r, ω) and j_s,_n₋_k(r, ω). Thus, in Eq. (3b), if the amplitude c_n_k(ω) of j_s,_n_k(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,_n_k(r, ω) in Ω_a because the Bloch periodic part ς_s,_n_k(r, ω) would not necessarily be orthogonal to other counterparts with identical k, namely,

\begin{array}{l} c_{n k} (ω) & = \frac{1}{Θ_{n k} (ω)} \int_{Ω_{a}} d r e^{- i k \cdot r} J_{s, k} (r) \cdot ς_{s, n - k} (r, ω) \\ = \frac{1}{Θ_{n k} (ω)} \int_{Ω_{a}} d r J_{s, k} (r) \cdot j_{s, n - k} (r, ω) . \end{array}

Equation (20) uniquely determines the expansion coefficient c_n_k(ω) 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,_n_k(ω)} have to repeat in {Δε_r,_n_′−_k(ω)} and vice versa. It also means that the corresponding degenerate subsets in {Δε_r,_n_k(ω)} 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–46]. Further understanding on this point shall be pursued in the future study.

Acknowledgments

This work is sponsored by Research Center for Applied Sciences, Academia Sinica, Taiwan, and National Science Council, Taiwan, under grant number NSC100-2112-M-001-002-MY2.

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). [CrossRef] [PubMed]

3. S. John, “Strong localization of photons in certain disordered dielectric superlattices,” Phys. Rev. Lett. 58, 2486–2489 (1987). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

7. S. Noda, M. Yokoyama, M. Imada, A. Chutinan, and M. Mochizuki, “Polarization mode control of two-dimensional photonic crystal laser by unit cell structure design,” Science 293, 1123–1125 (2001). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef]

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. [CrossRef]

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). [CrossRef] [PubMed]

15. A. P. Hibbins, B. R. Evans, and J. R. Sambles, “Experimental verification of designer surface plasmons,” Science 308, 670–672 (2005). [CrossRef] [PubMed]

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/SiO₂/Ag T-shaped periodical structure,” Opt. Express 19, 23698–23705 (2011). [CrossRef] [PubMed]

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). [CrossRef]

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). [CrossRef] [PubMed]

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 one- and two-dimensional periodic systems with metallic components in the presence of dissipation,” Phys. Rev. B 55, 7427–7444 (1997). [CrossRef]

22. S. Guo and S. Albin, “Simple plane wave implementation for photonic crystal calculations,” Opt. Express 11, 167–175 (2003). [CrossRef] [PubMed]

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). [CrossRef]

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). [CrossRef]

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). [CrossRef]

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). [CrossRef]

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). [CrossRef]

29. Y. Cao, Z. Hou, and Y. Liu, “Finite difference time domain method for band-structure calculations of two-dimensional phononic crystals,” Solid State Commun. 132, 539–543 (2004). [CrossRef]

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). [CrossRef]

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). [CrossRef]

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). [CrossRef]

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). [CrossRef]

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 whispering-gallery modes in a circular dielectric microcavity,” Opt. Quantum Electron. 36, 213–221 (2004). [CrossRef]

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). [CrossRef]

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). [CrossRef]

38. S. W. Chang, “Full frequency-domain approach to reciprocal microlasers and nanolasers-perspective from Lorentz reciprocity,” Opt. Express 19, 21116–21134 (2011). [CrossRef] [PubMed]

39. S. W. Chang, “Confinement factors and modal volumes of micro- and nanocavities invariant to integration regions,” IEEE J. Sel. Top. Quantum. Electron. 18, 1771–1780 (2012). [CrossRef]

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). [CrossRef]

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). [CrossRef]

46. A. D. Yaghjian, “Bidirectionality of reciprocal, lossy or lossless, uniform or periodic waveguides,” IEEE Microw. Wireless Compon. Lett. 17, 480–482 (2007). [CrossRef]

	Complex-ω method [Eq. (1a)]		GE problem [Eq. (5a)]
l	ω_n₀a/(2πc)	Q_n₀	ω_n₀a/(2πc)	Q_n₀	Δε_r,_n₀(ω_n₀)
1	0.33863	1.99319	0.31900	1.66022	−0.03136 − 0.20758i
2	1.00044	3522.60	1.00044	3522.60	(0.048423 − 2.8i) × 10⁻⁴
3	1.46613	47.614	1.46619	47.595	(−0.05718 − 1.9774i) × 10⁻²
4	2.00310	1179.21	2.00310	1179.21	(−0.00979 − 8.5i) × 10⁻⁴

Frequency-domain formulation of photonic crystals using sources and gain

Abstract

1. Introduction

2. Formulation

3. Calculations and discussions

3.1. Analytical solution of 1D periodic bilayer structure

3.2. Nondispersive and lossless structure of dielectric/dielectric bilayers

3.3. Nondispersive but lossy structure of dielectric/metal bilayers

3.4. Dispersive and lossy structure of dielectric/metal bilayers

4. Conclusion

Appendix: Biorthogonality between Bloch modes

Acknowledgments

References and links

Cited By

Figures (7)

Tables (1)

Equations (33)

Optics Express