Efficient numerical method for analyzing optical bistability in photonic crystal microcavities

Nonlinear optical effects can be enhanced by photonic crystal microcavities and be used to develop practical ultra-compact optical devices with low power requirements. The finite-difference time-domain method is the standard numerical method for simulating nonlinear optical devices, but it has limitations in terms of accuracy and efficiency. In this paper, a rigorous and efficient frequency-domain numerical method is developed for analyzing nonlinear optical devices where the nonlinear effect is concentrated in the microcavities. The method replaces the linear problem outside the microcavities by a rigorous and numerically computed boundary condition, then solves the nonlinear problem iteratively in a small region around the microcavities. Convergence of the iterative method is much easier to achieve since the size of the problem is significantly reduced. The method is presented for a specific two-dimensional photonic crystal waveguide-cavity system with a Kerr nonlinearity, using numerical methods that can take advantage of the geometric features of the structure. The method is able to calculate multiple solutions exhibiting the optical bistability phenomenon in the strongly nonlinear regime. © 2013 Optical Society of America OCIS codes: (000.4430) Numerical approximation and analysis; (190.1450) Bistability; (050.5298) Photonic crystals; (050.1755) Computational electromagnetic methods. References and links 1. H. Gibbs, Optical Bistability: Controlling Light with Light (Academic, New York, 1985). 2. C. M. Bowden and A. M. Zheltikov, “Nonlinear optics of photonic crystals,” J. Opt. Soc. Am. B 19, 2046–2048 (2002). 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). 5. M. Soljacic and J. D. Joannopoulos, “Enhancement of nonlinear effects using photonic crystals,” Nat. Mater. 3, 211–219 (2004). 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). 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). 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). 9. M. Midrio, “Shooting technique for the computation of the plane-wave reflection and transmission through onedimensional nonlinear inhomogeneous dielectric structres,” J. Opt. Soc. Am. B 18, 1866–1871 (2001). 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). 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). 12. P. K. Kwan and Y. Y. Lu, “Computing optical bistability in one-dimensional nonlinear structures,” Opt. Commun. 238, 169–175 (2004). 13. A. Talflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech, Norwood, MA, 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). 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). 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). 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). 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). 19. Z. Hu and Y. Y. Lu, “Efficient analysis of photonic crystal devices by Dirichlet-to-Neumann maps,” Opt. Express 16, 17383–17399 (2008). 20. R. W. Boyd, Nonlinear Optics (Academic, San Diego, CA, 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). 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). 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). 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). 25. J. Yuan and Y. Y. Lu, “Photonic bandgap calculations using Dirichlet-to-Neumann maps,” J. Opt. Soc. Am. A 23, 3217–3222 (2006). 26. L. N. Trefethen, Spectral Methods in MATLAB (Society for Industrial and Applied Mathematics, 2000). 27. L. Yuan and Y. Y. Lu, “Analyzing second harmonic generation from arrays of cylinders using the Dirichlet-toNeumann maps,” J. Opt. Soc. Am. B 26, 587–594 (2009).


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][3][4][5][6].
Except for simple one-dimensional problems [7][8][9][10][11][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 finitedifference 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 indexcontrast 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 coupledmode 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 waveguidecavity 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 waveguidecavity 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.

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 ε 1 and ε 2 , respectively.A waveguide is formed by reducing the radii of one array of cylinders to a 1 (a 1 < a) in the PhC.A microcavity is created by increasing the radius of one cylinder to a 2 (a 2 > 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 xand z-axis, respectively, the governing equation for polarized electromagnetic waves in a Kerr medium is [20] 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 0 = ω/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.The PhC waveguide is periodic in x with period L. If the medium is linear, the wave field in y x Fig. 1.A PhC waveguide-cavity system with a microcavity at the center and two semiinfinite waveguides.The truncated domain S (for m = 5) is the rectangle enclosed by the red lines.the waveguide is a superposition of Bloch modes.A Bloch mode is a solution of the Helmholtz equation given as 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 φ 0 ) with amplitude A 0 and it comes from x = −∞.For the PhC waveguide-cavity system shown in Fig. 1, where ε 1 = 12.25, ε 2 = 2.25, a = 0.25L, a 1 = a/3, and a 2 = 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.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 ω c L/(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).

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 de-note the cross section of the central cylinder by Ω c .Our boundary condition has the following form where ∂ Ω c denotes the boundary of Ω c , r is the radial coordinate with the center of Ω c as the origin, ∂ + r denotes one-side (from the exterior of Ω c ) partial derivative with respect to r, f is a function defined on ∂ Ω c , D is an operator that acts on functions defined on ∂ Ω c , and A 0 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 D is approximated by an M × M matrix.In the following, we present a numerical method for computing f and D 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 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.
On the two vertical edges of S, the boundary conditions are where φ 0 is the fundamental propagating Bloch mode of the PhC waveguide, L − and L + 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 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 2 .For the regular unit cell Ω jk with ( j, k) = (4, m + 1), we can find the DtN map Λ satisfying where v jk = u(x j , y) for y k−1 < y < y k , h jk = u(x, y k ) for x j−1 < x < x j , ∂ x v jk = ∂ x u(x j , y) for y k−1 < 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 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 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 where Λ jk and Λ jk for j, k = 1, 2, 3, 4 are N × N matrices, Λ j5 for j = 1, 2, 3, 4 are N × M matrices, Λ5k for k = 1, 2, 3, 4 are M × N matrices, and Λ55 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 ∂ x v jk can be evaluated from the DtN maps of the two unit cells Ω jk and Ω j+1,k .This gives rises to If the edge is on the left boundary of S (related to field v 0k ), we have where L − kn are blocks of the matrix approximation of L − , and g = [g 1 , g 2 , . . ., g 2m+1 ] T is the vector approximation of (L + − L − )φ 0 (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+1 , the equation is 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 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 B is the matrix related to the right hand side of Eq. ( 13).The sizes of matrices A and B are (30mN + 8N) × (30mN + 8N) and (30mN + 8N) × M, respectively, and the coefficient matrix A is sparse.Notice that v 3,m+1 , v 4,m+1 , h 4m , 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: where C is a (4N) × M matrix, f 0 is a vector of length 4N.Using Eq. ( 15) and the last row of Eq. ( 9), we can eliminate v 3,m+1 , v 4,m+1 , h 4m , h 4,m+1 , and obtain the boundary condition (3).In fact, the matrix D and vector f are given by Notice that only the vector f is related to the incident wave.

Nonlinear solver
With the boundary condition (3), we can solve the nonlinear Helmholtz equation ( 1) in the nonlinear cylinder Ω c .For that purpose, it is natural to use the polar coordinate system.Equation (1) becomes for r ∈ (0, a 2 ) and θ ∈ [0, 2π].Newton's method is widely used to solve nonlinear problems.For Eq. ( 17), it gives where u (n−1) is a previous iteration of the solution, u (n) is the current iteration to be determined, and ū(n) is the complex conjugate of u (n) .To start the iterative process, we need an initial guess u (0) .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 2 , a 2 ) by u(r, θ ) = u(−r, θ ) for − a 2 < r < 0 and θ = (θ + π)mod(2π), (19) so that Eq. ( 18) is valid for r ∈ (−a 2 , a 2 ) and θ ∈ [0, 2π].We discretize r and θ by r j = a 2 cos( jπ/q) for 0 ≤ j ≤ q, where q is an odd integer and M is an even integer (the number of discretization points on ∂ Ω c ).Notice that r 0 = a 2 , r q = −a 2 , 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 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 where W is the Chebyshev differentiation matrix [26], Ŵ is the middle (q − 1) × (q − 1) block of W , w w w 0 and w w w q are row vectors of length q − 1, and ŵ w w 0 and ŵ w w q are column vectors of length q − 1.The second order partial derivative with respect to r can be approximated by the matrix T = W 2 .Similarly, the second order partial derivative with respect to θ is approximated as where F is the M × M second order Fourier differentiation matrix [26].
The pseudospectral method applies the matrix approximations for ∂ r , ∂ 2 r and ∂ 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 where is a column vector for u (n) at all interior points, and In the above, ⊗ is the Kronecker product, T is the middle (q−1)×(q−1) block of matrix T , I is the M × M identity matrix, R = diag{1/r 1 , 1/r 2 , . . ., 1/r (q−1)/2 }, u u u where I 2 is the (M/2) × (M/2) identity matrix.On the boundary of Ω c , i.e. at r = a 2 , the partial derivative of u with respect to r can be evaluated using the first row of Eq. ( 21): where the superscript "−" indicates the one-side derivative taken from the interior of Ω c .Meanwhile, boundary condition (3) gives where the one-side derivative is taken from the exterior of Ω c .The continuity condition of ρ −1 ∂ r u gives us another relation between u u u and u u u 0 : where σ = 1 or σ = ε 1 /ε 2 for TE or TM polarization, respectively.If u u u (n−1) is given, we can solve Eqs. ( 23) and ( 26) to find u u u (n) and u u u (n) 0 .Since the complex conjugate of u u u (n) appears in Eq. ( 23), it is necessary to work with the real and imaginary parts of u u u (n) .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.

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 χ (3) = 2 × 10 −12 m 2 /V 2 , 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 normal-ized 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 0 corresponds to an incident Bloch mode with the amplitude A 0 = √ 30 × 10 4 ≈ 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 0 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 = √ 180 × 10 4 ≈ 1.3416 × 10 5 or power P in = 6P 0 (P in over 100 µm in z is approximately 2.7mW).The corresponding field patterns are shown in Fig. 6.
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 u u (0) = 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 nonlinear 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 u u (0) 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.

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.

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