Efficient numerical method for analyzing optical bistability in photonic crystal microcavities

Lijun Yuan; Ya Yan Lu

doi:10.1364/OE.21.011952

1. Introduction

Optical bistability in dielectric structures with the Kerr nonlinearity has been extensively studied since 1980s [1]. It has applications in all-optical devices such as optical switches, but these applications are often limited by the weak nonlinear response of conventional dielectric materials, thus require high operating power. The development of photonic crystals (PhCs) has made the applications of optical bistability feasible, because nonlinear effects can be dramatically enhanced by PhC microcavities [2–6].

Except for simple one-dimensional problems [7–12], not many numerical methods are available for simulating nonlinear optical phenomena, including optical bistability, in two- or three-dimensional structures. Most researchers use a time domain method, such as the finite-difference time-domain (FDTD) method [13], since the nonlinear terms are relatively easy to handle in the time domain. However, FDTD requires a small grid size to resolve large index-contrast and a small time step to ensure numerical stability. Furthermore, the implementation of realistic material dispersion and the truncation of unbounded domains (especially when the structure is inhomogeneous at infinity) can be rather complicated. For two-dimensional (2D) structures, the nonlinear interactions between light and Kerr media can be modeled by a nonlinear Helmholtz equation. In general, the nonlinear Helmholtz equation can be discretized by a standard numerical method such as the finite difference or finite element method, and solved by an iterative method, such as Newton’s method or the nested iteration scheme [14, 15]. However, this approach is not very popular, since the iterative methods often converge very slowly or fail to converge for interesting problems with strong nonlinear effects and unstable solutions. Meanwhile, more efficient but approximate frequency-domain methods have been developed for problems with special geometric properties. For example, the multi-scattering theory was used to analyze optical bistability in finite-size 2D photonic crystals with microcavities [16]. Effective discrete equations were derived by expanding the wave field in localized modes to analyze 2D PhC waveguides with nonlinear defects [17]. A perturbation theory and a coupled-mode theory have also been developed to analyze nonlinear PhC devices [6, 18]. However, all these methods are based on some analytic assumptions. For example, in the multi-scattering theory, the wave field in the nonlinear cylinders is assumed to be a constant. This assumption becomes invalid when the wavelength is comparable to the radii of the cylinders.

In this paper, we develop a rigorous and efficient numerical method for analyzing nonlinear optical phenomena enhanced by microcavities. The method is presented for a PhC waveguide-cavity structure with the Kerr nonlinearity, where the nonlinear effect is concentrated in a microcavity. Our approach is to reduce the nonlinear problem to a small region around the microcavity. The linear problem in the exterior domain is solved first (and only once) to establish a boundary condition around the nonlinear region. The small nonlinear problem is then solved by an iterative method. To take advantage of the many identical unit cells in the PhC waveguide-cavity structure, we use the Dirichlet-to-Neumann (DtN) map method [19] to solve the linear problem. For a linear unit cell, the DtN map is an operator that maps the wave field to its normal derivative on the cell boundary, and it is approximated by a small matrix. Only a few different DtN maps are needed, since most of the unit cells are identical. This paper is organized as follows. In section 2, the problem is formulated and linear results are presented. In section 3, the procedure for establishing the boundary condition for the small nonlinear problem is described. In section 4, we present a numerical method for solving the nonlinear Helmholtz equation based on Newton’s iterative scheme and a pseudospectral method. Numerical results for the nonlinear problem are given in section 5.

2. Problem formulation and linear results

We consider a 2D PhC waveguide-cavity structure as previously studied in [6, 18] and shown in Fig. 1. It is formed by introducing a defect cylinder and two semi-infinite waveguides in an otherwise perfectly periodic and infinite PhC. The background PhC consists of a square lattice of parallel and infinitely long circular dielectric cylinders surrounded by a homogenous medium with a lower refractive index. The lattice constant and the radii of the cylinders are L and a, respectively. The dielectric constants of the cylinders and the surrounding medium are ε₁ and ε₂, respectively. A waveguide is formed by reducing the radii of one array of cylinders to a₁ (a₁ < a) in the PhC. A microcavity is created by increasing the radius of one cylinder to a₂ (a₂ > a), it is placed in the waveguide core with three regular cylinders on each side. In a Cartesian coordinate system where the waveguide axis and the cylinder axes are parallel to the x- and z-axis, respectively, the governing equation for polarized electromagnetic waves in a Kerr medium is [20]

ρ \frac{\partial}{\partial x} (\frac{1}{ρ} \frac{\partial u}{\partial x}) + ρ \frac{\partial}{\partial y} (\frac{1}{ρ} \frac{\partial u}{\partial y}) + k_{0}^{2} [ε (x, y) + \frac{3}{4} χ^{(3)} (x, y) {| u |}^{2}] u = 0,

where u is the z-component of the electric or magnetic field for the transverse electric (TE) or transverse magnetic (TM) polarization, respectively, ε(x, y) is the dielectric function, k₀ = ω/c is the free space wavenumber, ω is the angular frequency, c is the speed of light in vacuum, χ⁽³⁾(x, y) is an element of the third order nonlinear susceptibility and it is only non-zero in nonlinear media, and ρ = 1 or ρ = ε for the TE or TM polarization, respectively.

Fig. 1 A PhC waveguide-cavity system with a microcavity at the center and two semi-infinite waveguides. The truncated domain S (for m = 5) is the rectangle enclosed by the red lines.

Download Full Size | PDF

The PhC waveguide is periodic in x with period L. If the medium is linear, the wave field in the waveguide is a superposition of Bloch modes. A Bloch mode is a solution of the Helmholtz equation given as

u (x, y) = ϕ (x, y) e^{i β x},

where ϕ(x, y) is periodic in x with period L and β is a constant. The Bloch modes can be found by solving an eigenvalue problem. We follow the approach developed in [21], where the eigenvalue problem (for a given frequency) is linear and the eigenvalue is related to β. In the following, we assume that the incident wave is the fundamental propagating Bloch mode (denoted as ϕ₀) with amplitude A₀ and it comes from x = −∞. For the PhC waveguide-cavity system shown in Fig. 1, where ε₁ = 12.25, ε₂ = 2.25, a = 0.25L, a₁ = a/3, and a₂ = 0.33L, we calculate the linear transmission spectrum using the DtN-map method formulated in [19]. The results are shown in Fig. 2, where T represents the relative power (normalized by the power of the incident wave) carried by the fundamental Bloch mode propagating towards x = +∞. It can be observed that almost 100% percent of the incident power is transmitted at frequency ωL/(2πc) = 0.28815. Figure 3 shows the magnitude of the total wave field at frequency ωL/(2πc) = 0.28815. We notice that the wave field is highly localized around the microcavity. Therefore, it is reasonable to consider the nonlinear effect only in the PhC microcavity even though all dielectric cylinders are nonlinear in principle.

Fig. 2 Transmission spectrum of the PhC waveguide-cavity system.

Download Full Size | PDF

Fig. 3 Magnitude of the wave field excited by an incident Bloch mode with amplitude A₀ = 1 at frequency ωL/(2πc) = 0.28815.

Download Full Size | PDF

As a linear structure, the PhC waveguide-cavity system shown in Fig. 1 has a leaky cavity mode. It is a special solution of the linear Helmholtz equation satisfying outgoing radiation conditions in the waveguides, but it only exists for a complex frequency. The real part of the complex frequency corresponds to the peak frequency in the transmission spectrum shown in Fig. 2, and the imaginary part gives the damping rate of the field amplitude. The wave field of the leaky cavity mode looks like the wave field at the peak frequency, but its magnitude actually increases as the horizontal distance to the microcavity increases and blows up at infinity. We calculate the leaky cavity mode by the DtN-map method presented in [22], and obtain the complex frequency ω_cL/(2πc) = 0.28815 − 0.00058i (the time dependence is assumed to be e^−iωt). Figure 4(a) shows the wave field of the leaky cavity mode to the left of the cavity center. Notice that the wave field increases in the waveguide core as the distance to the microcavity is increased. As a comparison, we show the wave field excited by an incident Bloch mode at the peak frequency ωL/(2πc) = 0.28815 in Fig. 4(b).

Fig. 4 (a): Magnitude of the leaky cavity mode at complex frequency ω_cL/(2πc) = 0.28815−0.00058i. (b): Magnitude of the wave field excited by an incident Bloch mode at frequency ωL/(2πc) = 0.28815. The fields are normalized for comparison.

Download Full Size | PDF

3. Domain reduction process

In the previous section, it is shown that the electromagnetic field excited by an incident Bloch mode is particularly strong in the microcavity. Since the third order nonlinear effect is proportional to the field intensity, it is reasonable to consider the nonlinear effect only in the central cylinder of the microcavity. The original problem can be considered as two coupling problems: a nonlinear problem in the central cylinder and a linear problem outside the cylinder. Our approach is to establish a boundary condition on the boundary of the central cylinder. Let us denote the cross section of the central cylinder by Ω_c. Our boundary condition has the following form

\partial_{r}^{+} u = 𝒟 u + A_{0} f on \partial Ω_{c},

where ∂Ω_c denotes the boundary of Ω_c, r is the radial coordinate with the center of Ω_c as the origin,

\partial_{r}^{+}

denotes one-side (from the exterior of Ω_c) partial derivative with respect to r, f is a function defined on ∂Ω_c, 𝒟 is an operator that acts on functions defined on ∂Ω_c, and A₀ is the amplitude of the incident Bloch mode. In a discrete approximation, ∂Ω_c is sampled by M points, f is approximated by a column vector of length M, and 𝒟 is approximated by an M × M matrix. In the following, we present a numerical method for computing f and 𝒟 by solving the linear problem outside Ω_c.

To solve the linear problem, the infinite PhC waveguide-cavity structure has to be truncated. For frequencies in the bandgap of the background PhC, the wave field decays exponentially to zero as |y| → ∞. Therefore, the wave field sufficiently far away form the waveguide core in the y direction can be accurately approximated by zero. The computation domain S is the rectangular region enclosed by the red lines in Fig. 1. It covers seven unit cells in the horizontal direction and 2m + 1 unit cells in the vertical direction. That is

S = {(x, y) | 0 < x < 7 L, - (m + 0.5) L < y < (m + 0.5) L},

where m is a positive integer (m = 5 in Fig. 1). Zero Dirichlet boundary conditions are imposed on the top and bottom edges of S, i.e.

u = 0, y = \pm (m + 0.5) L .

On the two vertical edges of S, the boundary conditions are

\frac{\partial u}{\partial x} = ℒ^{-} u + A_{0} (ℒ^{+} - ℒ^{-}) ϕ_{0}, x = 0,

\frac{\partial u}{\partial x} = ℒ^{+} u, x = 7 L,

where ϕ₀ is the fundamental propagating Bloch mode of the PhC waveguide, ℒ⁻ and ℒ⁺ are operators acting on functions of y. These operators are approximated by matrices when y is discretized. Equations (5) and (6) are derived by expanding the field in the waveguides (x < 0 or x > 7L) in a complete set of Bloch modes. A detailed derivation of these boundary conditions are given in [19].

The linear problem is governed by the linear Helmholtz equation in domain Ω_e (the exterior of Ω_c in S) with boundary conditions (4), (5), (6) and an arbitrary Dirichlet boundary condition on ∂Ω_c, i.e. u on ∂Ω_c is assumed to be given and arbitrary. The DtN-map method developed in [19, 23] can be used to solve the linear problem efficiently. In this method, the domain S is divided into unit cells as shown in Fig. 1. The unit cells are

Ω_{j k} = {(x, y) | x_{j - 1} < x < x_{j}, y_{k - 1} < y < y_{k}}

for j = 1, 2,...,7 and k = 1, 2,...,2m + 1, where x_j = jL and y_k = kL − (m + 0.5)L. There are only two different kinds of unit cells: the regular unit cell containing a cylinder of radius a and the defect unit cell Ω_4,m+1 containing the nonlinear cylinder of radius a₂. For the regular unit cell Ω_jk with (j, k) ≠ (4, m + 1), we can find the DtN map Λ satisfying

[\begin{matrix} \partial_{x} v_{j - 1, k} \\ \partial_{x} v_{j k} \\ \partial_{y} h_{j, k - 1} \\ \partial_{y} h_{j k} \end{matrix}] = Λ [\begin{matrix} v_{j - 1, k} \\ v_{j k} \\ h_{j, k - 1} \\ h_{j k} \end{matrix}],

where v_jk = u(x_j, y) for y_k₋₁ < y < y_k, h_jk = u(x, y_k) for x_j₋₁ < x < x_j, ∂_xv_jk = ∂_xu(x_j, y) for y_k₋₁ < y < y_k, etc. The DtN map of a unit cell is an operator that maps u to its normal derivative on the cell boundary. In a unit cell containing a circular cylinder, the wave field can be expanded in cylindrical waves as

u (x, y) = \sum_{n = - \infty}^{\infty} c_{n} ψ_{n} (r) e^{i n θ},

where (r, θ) are the local polar coordinates with the origin at the center of the cylinder, and the function ψ_n(r) is related to the Bessel functions. The matrix approximation of Λ can be constructed by truncating expansion (8) to finite terms and evaluating both u and its normal derivative on the boundary of Ω_jk[24, 25]. Notice that the DtN maps for all regular unit cells in Ω_e are identical, thus we only need to calculate Λ once. Beside the regular unit cells, the domain Ω_e contains Ω̃_4,m+1 which is the part of the defect unit cell Ω_4,m+1 with the nonlinear cylinder Ω_c removed. For Ω̃_4,m+1 we can find a DtN map Λ̃ such that

[\begin{matrix} \partial_{x} v_{3, m + 1} \\ \partial_{x} v_{4, m + 1} \\ \partial_{y} h_{4 m} \\ \partial_{y} h_{4, m + 1} \\ \partial_{r} u |_{\partial Ω_{c}} \end{matrix}] = \tilde{Λ} [\begin{matrix} v_{3, m + 1} \\ v_{4, m + 1} \\ h_{4 m} \\ h_{4, m + 1} \\ u |_{\partial Ω_{c}} \end{matrix}],

where u and its normal derivative on ∂Ω_c are also involved. A numerical approximation of Λ̃ can also be constructed by a cylindrical wave expansion in Ω̃_4,m+1. Using N points on each edge of the unit cell and M points on ∂Ω_c, the DtN maps Λ and Λ̃ are approximated by (4N) × (4N) and (4N + M) × (4N + M) matrices, respectively.

With the DtN maps Λ and Λ̃, an equation can be set up for each edge of the unit cells. First, we write Λ and Λ̃ in block forms as

Λ = [\begin{matrix} Λ_{11} & Λ_{12} & Λ_{13} & Λ_{14} \\ Λ_{21} & Λ_{22} & Λ_{23} & Λ_{24} \\ Λ_{31} & Λ_{32} & Λ_{33} & Λ_{34} \\ Λ_{41} & Λ_{42} & Λ_{43} & Λ_{44} \end{matrix}], \tilde{Λ} = [\begin{matrix} {\tilde{Λ}}_{11} & {\tilde{Λ}}_{12} & {\tilde{Λ}}_{13} & {\tilde{Λ}}_{14} & {\tilde{Λ}}_{15} \\ {\tilde{Λ}}_{21} & {\tilde{Λ}}_{22} & {\tilde{Λ}}_{23} & {\tilde{Λ}}_{24} & {\tilde{Λ}}_{25} \\ {\tilde{Λ}}_{31} & {\tilde{Λ}}_{32} & {\tilde{Λ}}_{33} & {\tilde{Λ}}_{34} & {\tilde{Λ}}_{35} \\ {\tilde{Λ}}_{41} & {\tilde{Λ}}_{42} & {\tilde{Λ}}_{43} & {\tilde{Λ}}_{44} & {\tilde{Λ}}_{45} \\ {\tilde{Λ}}_{51} & {\tilde{Λ}}_{52} & {\tilde{Λ}}_{53} & {\tilde{Λ}}_{54} & {\tilde{Λ}}_{55} \end{matrix}],

where Λ_jk and Λ̃_jk for j, k = 1, 2, 3, 4 are N × N matrices, Λ̃_j₅ for j = 1, 2, 3, 4 are N × M matrices, Λ̃₅_k for k = 1, 2, 3, 4 are M × N matrices, and Λ̃₅₅ is an M × M matrix. The equation for an edge is established by matching the normal derivative of u using the DtN maps of two neighboring unit cells or boundary conditions (5) and (6). For example, if a vertical edge (related to field v_jk) is not on the boundary of S or the defect unit cell, the x-derivative ∂_xv_jk can be evaluated from the DtN maps of the two unit cells Ω_jk and Ω_j_+1,_k. This gives rises to

\begin{matrix} Λ_{21} v_{j - 1, k} + (Λ_{22} - Λ_{11}) v_{j k} = Λ_{23} h_{j, k - 1} + Λ_{24} h_{j k} \\ - Λ_{12} v_{j + 1, k} - Λ_{13} h_{j + 1, k - 1} - Λ_{14} h_{j + 1, k} = 0. \end{matrix}

If the edge is on the left boundary of S (related to field v₀_k), we have

\sum_{n = 1, n \neq k}^{2 m + 1} ℒ_{k n}^{-} v_{0 n} + (Λ_{11} - ℒ_{k k}^{-}) v_{0 k} - Λ_{12} v_{1 k} - Λ_{13} h_{1, k - 1} - Λ_{14} h_{1 k} = - A_{0} g_{n},

where

ℒ_{k n}^{-}

are blocks of the matrix approximation of ℒ⁻, and g = [g₁, g₂,..., g₂_m₊₁]^T is the vector approximation of (ℒ⁺ − ℒ⁻)ϕ₀(0, y). For an edge on the boundary of the defect unit cell, we establish the equation using Λ and Λ̃. For example, for the edge with the corresponding field v_3,_m₊₁, the equation is

\begin{array}{l} Λ_{21} v_{2, m + 1} + (Λ_{22} - {\tilde{Λ}}_{11}) v_{3, m + 1} + Λ_{23} h_{3 m} + Λ_{24} h_{3, m + 1} - {\tilde{Λ}}_{12} v_{4, m + 1} - {\tilde{Λ}}_{13} h_{4 m} \\ - {\tilde{Λ}}_{14} h_{4, m + 1} = {\tilde{Λ}}_{15} u |_{\partial Ω_{c}}, \end{array}

where u|_∂_Ω_c is assumed to be given.

Assembling all equations like (11), (12) and (13) together, and using the zero boundary condition (4), we obtain the following linear system

𝒜 U = ℬ u |_{Ω_{c}} + A_{0} G,

where U is a vector for wave field u on the edges of all unit cells in domain S or on the vertical boundaries of S, G is a vector related to the right hand side of Eq. (12), and ℬ is the matrix related to the right hand side of Eq. (13). The sizes of matrices 𝒜 and ℬ are (30mN + 8N) × (30mN + 8N) and (30mN + 8N) × M, respectively, and the coefficient matrix 𝒜 is sparse. Notice that v_3,_m₊₁, v_4,m+1, h₄_m, and h_4,m+1 (related to the edges of the defect unit cell Ω_4,m+1) are components of vector U. Solving Eq. (14), we obtain the following relation:

[\begin{matrix} v_{3, m + 1} \\ v_{4, m + 1} \\ h_{4 m} \\ h_{4, m + 1} \end{matrix}] = 𝒞 u |_{\partial Ω_{c}} + A_{0} f_{0},

where 𝒞 is a (4N) × M matrix, f₀ is a vector of length 4N. Using Eq. (15) and the last row of Eq. (9), we can eliminate v_3,_m₊₁, v_4,m+1, h₄_m, h_4,m+1, and obtain the boundary condition (3). In fact, the matrix 𝒟 and vector f are given by

𝒟 = [\begin{matrix} {\tilde{Λ}}_{51} & {\tilde{Λ}}_{52} & {\tilde{Λ}}_{53} & {\tilde{Λ}}_{54} \end{matrix}] 𝒞 + {\tilde{Λ}}_{55}, f = [\begin{matrix} {\tilde{Λ}}_{51} & {\tilde{Λ}}_{52} & {\tilde{Λ}}_{53} & {\tilde{Λ}}_{54} \end{matrix}] f_{0} .

Notice that only the vector f is related to the incident wave.

4. Nonlinear solver

With the boundary condition (3), we can solve the nonlinear Helmholtz equation (1) in the non-linear cylinder Ω_c. For that purpose, it is natural to use the polar coordinate system. Equation (1) becomes

\frac{\partial^{2} u}{\partial r^{2}} = \frac{1}{r} \frac{\partial u}{\partial r} + \frac{1}{r^{2}} \frac{\partial^{2} u}{\partial θ^{2}} + k_{0}^{2} (ε_{1} + \frac{3}{4} χ^{(3)} {| u |}^{2}) u = 0

for r ∈ (0, a₂) and θ ∈ [0, 2π]. Newton’s method is widely used to solve nonlinear problems. For Eq. (17), it gives

\begin{array}{l} (\frac{\partial^{2}}{\partial r^{2}} + \frac{1}{r} \frac{\partial}{\partial r} + \frac{1}{r^{2}} \frac{\partial^{2}}{\partial θ^{2}} + k_{0}^{2} ε_{1}) u^{(n)} + \frac{3}{4} k_{0}^{2} χ^{(3)} [2 {| u^{(n - 1)} |}^{2} u^{(n)} + {(u^{(n - 1)})}^{2} {\bar{u}}^{(n)}] \\ = \frac{3}{3} k_{0}^{2} χ^{(3)} {| u^{(n - 1)} |}^{2} u^{(n - 1)}, \end{array}

where u⁽ⁿ⁻¹⁾ is a previous iteration of the solution, u⁽ⁿ⁾ is the current iteration to be determined, and ū⁽ⁿ⁾ is the complex conjugate of u⁽ⁿ⁾. To start the iterative process, we need an initial guess u⁽⁰⁾.

To solve the inhomogeneous Helmholtz equation (18), we use a mixed Fourier-Chebyshev pseudospectral method [26, 27]. More precisely, the Chebyshev collocation method is used to discretize the radial direction and the Fourier pseudospectral method is used to discretize the angular direction. To avoid the singularity at r = 0, we follow the approach presented in [26], and extend the radial variable r to (−a₂, a₂) by

u (r, θ) = u (- r, \tilde{θ}) for - a_{2} < r < 0 and \tilde{θ} = (θ + π) mod (2 π),

so that Eq. (18) is valid for r ∈ (−a₂, a₂) and θ ∈ [0, 2π]. We discretize r and θ by

\begin{array}{l} r_{j} = a_{2} cos (j π / q) for 0 \leq j \leq q, \\ θ_{k} = (2 k - 1) π / M for 1 \leq k \leq M, \end{array}

where q is an odd integer and M is an even integer (the number of discretization points on ∂Ω_c). Notice that r₀ = a₂, r_q = −a₂, and r = 0 is avoided since q is an odd integer. Denoting u(r_j, θ_k) by u_jk, then Eq. (19) implies that for (q + 1)/2 ≤ j ≤ q

u_{j k} = {\begin{array}{l} u_{q - j, k + M / 2}, for 1 \leq k \leq M / 2, \\ u_{q - j, k - M / 2}, for 1 + M / 2 \leq k \leq M, \end{array}

Therefore, although the computation domain is doubled by the extension of r to negative values, the total number of unknowns still corresponds to the original domain. In the pseudospectral method, the partial derivative with respect to r is approximated as

\frac{\partial}{\partial r} [\begin{matrix} u_{0 k} \\ u_{1 k} \\ ⋮ \\ u_{q k} \end{matrix}] \approx 𝒲 [\begin{matrix} u_{0 k} \\ u_{1 k} \\ ⋮ \\ u_{q k} \end{matrix}] = [\begin{matrix} w_{00} & {\tilde{w}}_{0} & w_{0 q} \\ {\hat{w}}_{0} & \hat{𝒲} & {\hat{w}}_{q} \\ w_{q 0} & {\tilde{w}}_{q} & w_{q q} \end{matrix}] [\begin{matrix} u_{0 k} \\ u_{1 k} \\ ⋮ \\ u_{q k} \end{matrix}],

where 𝒲 is the Chebyshev differentiation matrix [26], 𝒲̂ is the middle (q − 1) × (q − 1) block of 𝒲, w̃₀ and w̃_q are row vectors of length q − 1, and ŵ₀ and ŵ_q are column vectors of length q − 1. The second order partial derivative with respect to r can be approximated by the matrix 𝒯 = 𝒲². Similarly, the second order partial derivative with respect to θ is approximated as

\frac{\partial^{2}}{\partial θ^{2}} [\begin{matrix} u_{j 1} \\ u_{j 2} \\ ⋮ \\ u_{j M} \end{matrix}] \approx ℱ [\begin{matrix} u_{j 1} \\ u_{j 2} \\ ⋮ \\ u_{j M} \end{matrix}],

where ℱ is the M × M second order Fourier differentiation matrix [26].

The pseudospectral method applies the matrix approximations for ∂_r, $\partial_{r}^{2}$ and $\partial_{θ}^{2}$ to Eq. (18), and assumes that the approximate equation remains valid at all interior points (r_j, θ_k) for 1 ≤ j ≤ q − 1 and 1 ≤ k ≤ M. With a further application of the symmetry relation (20), Eq. (18) is discretized as

ℋ_{1} u^{(n)} + 2 γ diag {{| u^{(n - 1)} |}^{2}} u^{(n)} + γ diag {{[u^{(n - 1)}]}^{2}} {\bar{u}}^{(n)} + ℋ_{2} u_{0}^{(n)} = γ diag {{| u^{(n - 1)} |}^{2}} u^{(n)},

where

u^{(n)} = {[u_{11}^{(n)}, u_{12}^{(n)}, \dots, u_{1 M}^{(n)}, u_{21}^{(n)}, u_{22}^{(n)}, \dots, u_{2 M}^{(2)}, \dots, u_{(q - 1) / 2, 1}^{(n)}, \dots, u_{(q - 1) / 2, M}^{(n)}]}^{T}

is a column vector for u⁽ⁿ⁾ at all interior points, and

ℋ_{1} = \hat{ℐ} \otimes I + (ℛ \hat{𝒲}) \otimes I + ℛ^{2} \otimes ℱ + k_{0}^{2} ε_{1}, ℋ_{2} = {\hat{w}}_{0} \otimes I + {\hat{w}}_{q} \otimes \hat{I} .

In the above, ⊗ is the Kronecker product, 𝒯̂ is the middle (q−1)×(q−1) block of matrix 𝒯, I is the M × M identity matrix, ℛ = diag{1/r₁, 1/r₂,..., 1/r_(q−1)/2},

u_{0}^{(n)} = {[u_{01}^{(n)}, u_{02}^{(n)}, \dots, u_{0 M}^{(n)}]}^{T}

is the vector for u⁽ⁿ⁾ on ∂Ω_c, Î is the M × M matrix given by

\hat{I} = [\begin{matrix} 0 & I_{2} \\ I_{2} & 0 \end{matrix}],

where I₂ is the (M/2) × (M/2) identity matrix. On the boundary of Ω_c, i.e. at r = a₂, the partial derivative of u with respect to r can be evaluated using the first row of Eq. (21):

\partial_{r}^{-} u_{0}^{(n)} = {\tilde{w}}_{0} \otimes I u^{(n)} + (w_{00} + w_{0 q} \hat{I}) u_{0}^{(n)},

where the superscript “−” indicates the one-side derivative taken from the interior of Ω_c. Meanwhile, boundary condition (3) gives

\partial_{r}^{+} u_{0}^{(n)} = 𝒟 u_{0}^{(n)} + A_{0} f,

where the one-side derivative is taken from the exterior of Ω_c. The continuity condition of ρ⁻¹∂_ru gives us another relation between u and u₀ :

{\tilde{w}}_{0} \times I u^{(n)} + (w_{00} + w_{0 q} \hat{I} - σ 𝒟) u_{0}^{(n)} = A_{0} σ f,

where σ = 1 or σ = ε₁/ε₂ for TE or TM polarization, respectively. If u⁽ⁿ⁻¹⁾ is given, we can solve Eqs. (23) and (26) to find u⁽ⁿ⁾ and

u_{0}^{(n)}

. Since the complex conjugate of u⁽ⁿ⁾ appears in Eq. (23), it is necessary to work with the real and imaginary parts of u⁽ⁿ⁾. On the other hand, since the structure is symmetric about the x-axis, we can reduce the size of the linear system by one half. As a result, we only need to solve a (qM/2) × (qM/2) linear system in each iteration.

5. Numerical results: nonlinear case

In this section, we present numerical results for the PhC waveguide-cavity system shown in Fig. 1, where the cylinder Ω_c at the center of the structure has a third order nonlinearity with χ⁽³⁾ = 2 × 10⁻¹²m²/V², while all other cylinders and the surrounding medium are linear. The lattice constant of the background PhC is assumed to be L = 1μm. The radii of the cylinders and the dielectric constants are given in section 2. The structure exhibits the optical bistability phenomenon. In Fig. 5, we show the normalized output power as a function of the normalized input power for ωL/(2πc) = 0.286397. The input and output power are related to the incoming and outgoing propagating Bloch modes in the left and right semi-infinite PhC waveguides respectively, and they are proportional to the squared modulus of the coefficients of the Bloch mode. The reference power P₀ corresponds to an incident Bloch mode with the amplitude $A_{0} = \sqrt{30} \times 10^{4} \approx 54772$ . Since the structure is infinite and invariant in z, it is only possible to consider the power of the Bloch mode over a given interval of z. Based on the field profiles of the Bloch mode, we found that P₀ over 100μm in z is approximately 0.45mW. From Fig. 5, it is clear that the nonlinear Helmholtz equation has three solutions for certain range of the incident power. Presumably, the solutions corresponding to low and high output power are stable, while the solution corresponding to the medium output power is unstable. Their stability cannot be analyzed using Eq. (1). Instead, it should be analyzed using equations derived from the original time-domain nonlinear Maxwell’s equations through a linearization process around these solutions. The three solutions marked by A, B and C in Fig. 5 are obtained with an incident Bloch mode with amplitude $A_{0} = \sqrt{180} \times 10^{4} \approx 1.3416 \times 10^{5}$ or power P_in = 6P₀ (P_in over 100μm in z is approximately 2.7mW). The corresponding field patterns are shown in Fig. 6.

Fig. 5 Normalized output power as a function of the nomalized incident power at frequency ωL/(2πc) = 0.286397.

Download Full Size | PDF

Fig. 6 Wave field patterns (magnitude of u) corresponding to the three solutions marked as A, B and C in Fig. 5.

Download Full Size | PDF

The results in Fig. 5 and Fig. 6 are obtained with m = 10, N = 11, M = 44, and q = 101, where m is the number of remaining unit cells in each side of the waveguide in the truncated domain, N is the number of points for discretizing each edge of the square unit cells, M is the number of points for discretizing the boundary of the nonlinear cylinder Ω_c, and q is the number of point for discretizing the radial variable r in Ω_c. On a personal computer with a 2.4GHz CPU, the boundary condition (3) can be constructed within 8 seconds. For the nonlinear problem in Ω_c, the linear system for each iteration involves a 2222 × 2222 coefficient matrix, and it can be solved in less than 1 second.

The lower branch of the solution curve in Fig. 5 can be easily obtained with a zero initial guess u⁽⁰⁾ = 0. Typically, Newton’s method converges in less than 20 iterations. It can be calculated more efficiently by a continuation scheme, where the solution obtained for one value of P_in is used as the initial guess for a slightly larger P_in. For sufficiently large P_in, the non-linear problem has only one solution. In that case, Newton’s method again converges for the zero initial guess. The upper branch of the solution curve is then calculated by a continuation scheme that slightly decreases P_in in each step. Since the middle branch is unstable, it cannot be calculated by time-domain methods, such as the nonlinear FDTD [6, 18], but it can still be calculated in the frequency-domain. To find the solution marked as B in Fig. 5, we start Newton’s iteration where u⁽⁰⁾ is the product of solution C with a number s ∈ (0, 1). For a properly chosen s, the iterations converge to a solution (i.e. solution B) which is different from A or C. Once the solution B is obtained, we use the continuation scheme to calculate the left and right parts of the middle branch where P_in is decreased and increased in each step, respectively.

6. Conclusions

In the previous sections, we developed an efficient frequency-domain numerical method to analyze a PhC waveguide-cavity system which exhibits the optical bistability phenomenon. The basic idea is to assume the nonlinear effect is only important in a small region around the microcavity, remove the linear part by calculating a boundary condition for the small nonlinear region, and solve the small nonlinear problem iteratively. Since the nonlinear problem is solved in a small region, the number of unknowns involved in each iteration is significantly reduced, therefore, convergence is much easier to achieve. Although we have only worked on the particular two-dimensional structure involving a PhC microcavity connected by two semi-infinite PhC waveguides, the basic idea should be applicable to any problem where the nonlinear effect is dominated in a small region.

For the PhC structure considered in this paper, we make use of the so-called Dirichlet-to-Neumann (DtN) maps of the unit cells to calculate Bloch modes of the PhC waveguides as in [21], truncate domains with PhC waveguides extending to infinity as in [19], set up linear system of equations on edges of the unit cells only, and finally replace the linear part by a boundary condition. The nonlinear problem is solved iteratively by Newton’s method and the discretization is based on a combined Fourier-Chebyshev pseudospectral method. The optical bistability of the structure appears in a strong nonlinear regime where multiple solutions exist. To find the relation between the input and output power, we have also used a continuation scheme to improve the efficiency. The DtN map and the pseudospectral methods are used to take advantage of the the special geometry of the concerned structure. For a general problem with a microcavity enhancing the nonlinearity, a more general numerical method, such as the finite element method, may be more appropriate.

Acknowledgments

This research was partially supported by the National Natural Science Foundation of China (project No. 11201508), the Scientific Research Foundation of Chongqing Technology and Business University (project No. 20125605), and by City University of Hong Kong (project No. 7002864).

References and links

1. H. Gibbs, Optical Bistability: Controlling Light with Light (Academic, 1985).

2. C. M. Bowden and A. M. Zheltikov, “Nonlinear optics of photonic crystals,” J. Opt. Soc. Am. B 19,2046–2048 (2002) [CrossRef] .

3. R. E. Slusher and B. J. Eggleton, Nonlinear Photonic Crystals (Springer-Verlag, Berlin, 2003).

4. K. J. Vahala, “Optical microcavities,” Nature 424,839–846 (2003) [CrossRef] [PubMed] .

5. M. Soljacic and J. D. Joannopoulos, “Enhancement of nonlinear effects using photonic crystals,” Nat. Mater. 3,211–219 (2004) [CrossRef] [PubMed] .

6. J. Bravo-Abad, A. Rodriguez, P. Bermel, S. G. Johnson, J. D. Joannopoulos, and M. Soljacic, “Enhance nonlinear optics in photonic-crystal microcavities,” Opt. Express 15,16161–16176 (2007) [CrossRef] [PubMed] .

7. G. S. Agarwal and S. D. Gupta, “Effect of nonlinear boundary conditions on nonlinear phenomena in optical resonators,” Opt. Lett. 12,829–831 (1987) [CrossRef] [PubMed] .

8. J. Danckaert, K. Fobelets, I. Veretennicoff, G. Vitrant, and R. Reinisch, “Dispersive optical bistability in stratified structures,” Phys. Rev. B 44,8214–8225 (1991) [CrossRef] .

9. M. Midrio, “Shooting technique for the computation of the plane-wave reflection and transmission through one-dimensional nonlinear inhomogeneous dielectric structres,” J. Opt. Soc. Am. B 18,1866–1871 (2001) [CrossRef] .

10. A. Suryanto, E. van Groesen, M. Hammer, and H. J. W. M. Hoekstra, “A finite element scheme to study the nonlinear optical response of a finite grating without and with defect,” Opt. Quant. Electron. 35,313–332 (2003) [CrossRef] .

11. A. Suryanto, E. van Groesen, and M. Hammer, “Finite element analysis of optical bistability in one-dimensional nonlinear photonic band gap structures with defect,” J. Nonlinear Opt. Phy. Mater. 12,187–204 (2003) [CrossRef] .

12. P. K. Kwan and Y. Y. Lu, “Computing optical bistability in one-dimensional nonlinear structures,” Opt. Commun. 238,169–175 (2004) [CrossRef] .

13. A. Talflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech, 2000).

14. G. Baruch, G. Fibich, and S. Tsynkov, “A high-order numerical method for the nonlinenar Helmholtz equation in multi-dimensional layered media,” J. Comput. Phys. 228,3789–3815 (2009) [CrossRef] .

15. Z. Xu and G. Bao, “A numerical scheme for nonlinear Helmholtz equations with strong nonlinear optical effects,” J. Opt. Soc. Am. A 27,2347–2353 (2010) [CrossRef] .

16. E. Centeno and D. Felbacq, “Optical bistability in finite-size nonlinear bidimensional photonic crystals doped by a microcavity,” Phys. Rev. B 62,R7683–R7686 (2000) [CrossRef] .

17. S. F. Mingaleev and Y. S. Kivshar, “Nonlinear transmission and light localization in photonic-crystal waveguides,” J. Opt. Soc. Am. B 19,2241–2249 (2002) [CrossRef] .

18. J. Bravo-Abad, S. Fan, S. G. Johnson, J. D. Joannopoulos, and M. Soljacic, “Modeling nonlinear optical phenomena in nanophotonics,” J. Lightw. Technol. 25,2539–2546 (2007) [CrossRef] .

19. Z. Hu and Y. Y. Lu, “Efficient analysis of photonic crystal devices by Dirichlet-to-Neumann maps,” Opt. Express 16,17383–17399 (2008) [CrossRef] [PubMed] .

20. R. W. Boyd, Nonlinear Optics (Academic, 1992).

21. Y. Huang, Y. Y. Lu, and S. Li, “Analyzing photonic crystal waveguides by Dirichlet-to-Neumann maps,” J. Opt. Soc. Am. B 24,2860–2867 (2007) [CrossRef] .

22. S. Li and Y. Y. Lu, “Efficient method for analyzing leaky cavities in two-dimensional photonic crystals,” J. Opt. Soc. Am. B 26,2427–2433 (2009) [CrossRef] .

23. Z. Hu and Y. Y. Lu, “A simple boundary condition for terminating photonic crystal waveguides,” J. Opt. Soc. Am. B 29,1356–1360 (2012) [CrossRef] .

24. Y. Huang and Y. Y. Lu, “Scattering from periodic arrays of cylinders by Dirichlet-to-Neumann maps,” J. Lightw. Technol. 24,3448–3453 (2006) [CrossRef] .

25. J. Yuan and Y. Y. Lu, “Photonic bandgap calculations using Dirichlet-to-Neumann maps,” J. Opt. Soc. Am. A 23,3217–3222 (2006) [CrossRef] .

26. L. N. Trefethen, Spectral Methods in MATLAB (Society for Industrial and Applied Mathematics, 2000) [CrossRef] .

27. L. Yuan and Y. Y. Lu, “Analyzing second harmonic generation from arrays of cylinders using the Dirichlet-to-Neumann maps,” J. Opt. Soc. Am. B 26,587–594 (2009) [CrossRef] .

Efficient numerical method for analyzing optical bistability in photonic crystal microcavities

Abstract

1. Introduction

2. Problem formulation and linear results

3. Domain reduction process

4. Nonlinear solver

5. Numerical results: nonlinear case

6. Conclusions

Acknowledgments

References and links

Cited By

Figures (6)

Equations (32)

Optics Express