Numerical simulation of optical coupling and light propagation in coupled optical resonators with size disorder

Applied Numerical Mathematics 57 (2007) 475–485 www.elsevier.com/locate/apnum

Numerical simulation of optical coupling and light propagation in coupled optical resonators with size disorder Shaozhong Deng Department of Mathematics and Statistics, University of North Carolina at Charlotte, Charlotte, NC 28223-0001, USA Available online 24 July 2006

Abstract We use a discontinuous spectral element method for the study of optical coupling via whispering gallery modes between two size mismatched dielectric microresonators, the building block of coupled resonator optical waveguides, and light propagation through chains of size mismatched dielectric microresonators. It is shown that such optical coupling does not depend monotonically on the size dispersion, and on the other hand, such light propagation is highly sensitive to the size dispersion in the coupled resonator optical waveguides. © 2006 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Whispering gallery modes; Coupled resonator optical waveguides; Discontinuous spectral element methods

1. Introduction In order to control the velocity of light and provide routing and switching functionality on a microscopic scale, there has been much active recent research in coupled optical resonators, and devices such as Coupled Resonator Optical Waveguides (CROWs) have been proposed [2,16–18,29,31]. Such waveguides can provide the possibilities of controlling light propagation speed and manipulating light propagation path in a very different way from that in the conventional optical fiber waveguides, and thereby have potential applications in telecommunication such as optical buffers and delay lines, nonlinear elements, and sensors. Theoretically, it is well known that in CROW devices, waveguiding is achieved by optical coupling of evanescent whispering gallery modes (WGMs) in individual optical resonators, and the mechanism of such optical coupling can be well understood qualitatively by the tight binding method [16,29,31]. However, the formalism is not convenient for real physical CROW structures because, for example, it does not account for different resonator sizes. Moreover, from a designer’s point of view, more quantitative study of optical coupling and light propagation in CROW structures is needed to address various issues such as the sensitivity of the coupling strength and the propagation efficiency to the inter-resonator gap size, and the resonator size dispersion of those structures due to the quality of their manufacturing, etc. With the aim of analyzing realistic CROW structures and also considering that recent research on CROW devices is mostly experimental, we have developed a numerical approach to study optical coupling and light propagation in CROW devices of size mismatched resonators. E-mail address: [email protected] (S. Deng). 0168-9274/$30.00 © 2006 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2006.07.001

476

S. Deng / Applied Numerical Mathematics 57 (2007) 475–485

Previous numerical simulations of optical coupling and light propagation in CROW structures have been done by assuming that all coupled resonators are identical, and various interesting findings have been obtained [5,7]. However, the performance of real physical CROW structures will be dependent on the size disorder of resonators which could have a profound, but not well studied, effect on optical coupling and light propagation. The optical coupling and transport phenomena in size mismatched microspheres has been recently investigated experimentally [2]. To provide valuable insight for the feasibility of, and critical guideline to the size disorder tolerance in, designing and building general practical CROW devices based on coupled microresonators, in this paper we shall specifically study optical coupling and light propagation in coupled microcylinder optical waveguides with slight size disorder. 2. Formulation In CROW systems, optical coupling occurs through a “photon hopping” transport between adjacent resonators at the frequencies of WGMs in individual resonators. Generally speaking, WGMs are resonance modes traveling in media of circular symmetric structures. The WGMs were first studied by Lord Rayleigh trying to understand the acoustic waves clinging to the dome of St. Paul’s Cathedral [19] and were then investigated by James Wait for a dielectric rod [26] including the explicit formulation which is reviewed below. 2.1. Cylinder fields of electromagnetic WGMs We shall consider a circular dielectric cylinder of radius a and infinite length with dielectric constant 1 and magnetic permeability μ1 , which is embedded in an infinite homogeneous medium of material parameters 2 and μ2 . As is well known for problems of this kind, the fields may be derived entirely from the axial components of the electric and magnetic fields. With respect to a cylindrical coordinate system (r, θ, z), for a time factor e−iωt , these axial components are solutions of Helmholtz equations (∇ 2 + k12 )ψ = 0 for r < a, and (∇ 2 + k22 )ψ = 0 for r > a, √ √ where ψ is either Ez or Hz . And k1 = ω 1 μ1 and k2 = ω 2 μ2 are the propagation constants inside and outside the cylinder, respectively. The fields must be finite at the center and, consequently, the wave functions within the cylinder will be represented by Bessel functions of the first kind Jn . Outside the cylinder, Hankel functions of the first kind (1) Hn ensure the Sommerfeld outgoing conditions at the infinity. Therefore, solutions are of the type  r < a, Jn (λ1 r), ψ = einθ+ihz−iωt (1) (1) Hn (λ2 r), r > a, where h is the axial propagation constant, λ21 = k12 − h2 , and λ22 = k22 − h2 . Representing the field components Ez and Hz in terms of solutions in (1), and then applying Maxwell’s equations, we can have the components of the magnetic field H = (Hr , Hθ , Hz ) and the electric field E = (Er , Eθ , Ez ) by the following equations [22,26]  ∞  2  s s nk s ih  an Hr = Gn (λr) + bn Gn (λr) Fn , λ μωλ2 r n=−∞  ∞   ik 2  nh ans Gn (λr) − bns 2 Gn (λr) Fn , Hθs = μωλ λ r n=−∞ Hzs = Ers

=

∞ 

bns Gn (λr)Fn ,

n=−∞ ∞   n=−∞ ∞ 

Eθs = − Ezs =



ih μωn ans Gn (λr) − bns 2 Gn (λr) λ

λ r

Fn ,

  nh iμω  ans 2 Gn (λr) + bns Gn (λr) Fn , λ λ r n=−∞

∞  n=−∞

ans Gn (λr)Fn ,

S. Deng / Applied Numerical Mathematics 57 (2007) 475–485

477

(1)

where Fn = einθ+ihz−iωt , and the function Gn ≡ Jn for r < a and Hn for r > a. Prime denotes differentiation with respect to the argument λr. Also, for r < a, μ = μ1 , k = k1 , λ = λ1 , and for r > a, μ = μ2 , k = k2 , λ = λ2 . Here, the superscript “s” is either “i” or “e”, representing “interior” and “exterior”, respectively. The coefficients of the expansions ani , bni , ane , bne , and the axial propagation constant h are determined by boundary conditions. At the cylindrical boundary r = a, the tangential components of the fields are continuous, yielding the following relations among the coefficients for each mode number n ik22 ik12  nh nh  Jn (u)ani − 2 Jn (u)bni = H (1) (v)ane − 2 Hn(1) (v)bne , μ1 ωu μ2 ωv n u v

(2)

Jn (u)bni = Hn(1) (v)bne , iμ1 ω  nh iμ2 ω (1) nh Jn (u)bni = 2 Hn(1) (v)ane + Hn (v)bne , Jn (u)ani + 2 u v u v Jn (u)ani = Hn(1) (v)ane ,

(3) (4) (5)

where u = λ1 a, v = λ2 a. Eqs. (2)–(5) constitute a homogeneous system of linear equations for the four coefficients ani , bni , ane , and bne . For a nontrivial solution, the determinant of the corresponding 4 × 4 coefficient matrix should vanish, resulting in an algebraic equation for the axial propagation constant h [26]     2    (1) k22 Hn(1) (v) k1 Jn (u) μ1 Jn (u) μ2 Hn (v) 1 2 2 2 1 − − =n h − , (6) u Jn (u) v Hn(1) (v) μ1 u Jn (u) μ2 v Hn(1) (v) v 2 u2 whose roots are the allowed values of the propagation constant h and determine the characteristics or natural modes of propagation. Note that for a given mode number n, Eq. (6) does not have a unique solution. As stated earlier, electromagnetic WGMs are resonance modes that can be excited in a resonator, and would circulate around the inside of the resonator’s boundary. However, the phase matching condition dictates that only at certain √ wavelengths does such resonance occur. In our case, that is nλ ≈ 2πa r μr , where λ is the azimuthal wavelength, r is the relative permittivity ( = r 0 ) and μr the relative permeability (μ = μr μ0 ) of the cylinder. In other words, the electromagnetic WGMs are represented by solutions of Eq. (6) when n is of the order of λ1 a [19]. Note that the mode number n is also the number of maxima of the field intensity in the azimuthal direction and is called the azimuthal number of the WGM. 2.2. Maxwell’s equations for coupled microcylinders In order to investigate optical coupling and light propagation via WGMs in microcylinder CROW devices, we shall turn to Maxwell’s equations. For a WGM with the axial propagation constant h, the magnetic field H = (Hx , Hy , Hz ) and the electric field E = (Ex , Ey , Ez ) in a rectangular coordinate system (x, y, z) may be expressed as H(x, y, z, t) = H(x, y, t)eihz ,

E(x, y, z, t) = E(x, y, t)eihz .

Substituting them into the three-dimensional Maxwell’s equations ∂E ∂H = −∇ × E,  = ∇ × H, μ ∂t ∂t we obtain the following hyperbolic system of equations in the metric system of units ∂Q ∂Q ∂Q + A(, μ) + B(, μ) = S, ∂t ∂x ∂y where Q = (μH, E)T and S = (ihEy , −ihEx , 0, −ihHy , ihHx , 0)T , and ⎛ ⎛ ⎞ 0 0 0 0 0 0 0 0 0 0 0 −1/ ⎟ ⎜ 0 ⎜0 ⎜ ⎜ ⎟ 0 0 0 1/ 0 ⎟ ⎜ 0 ⎜0 B(, μ) = ⎜ A(, μ) = ⎜ ⎟, 0 0 0 0 0 ⎟ ⎜ 0 ⎜0 ⎝ ⎝ ⎠ 0 0 0 1/μ 0 0 0 1/μ 0 −1/μ 0 0 0 0

0 0 0 0 0 0 0 0 −1/ 0 −1/μ 0 0 0 0 0 0 0

(7)

0 0 0 0 0 0

⎞ 1/ 0 ⎟ ⎟ 0 ⎟ ⎟. 0 ⎟ ⎠ 0 0

478

S. Deng / Applied Numerical Mathematics 57 (2007) 475–485

2.3. Discontinuous spectral element method (DSEM) for time domain Maxwell’s equations The finite difference time domain (FDTD) methods have been used extensively in computational electromagnetics [24,32]. The algorithms enforce Gauss’s Law, handle material interface naturally, and are easy to code. However, they have large phase errors, work well only for simple geometries or Cartesian grids, and cannot implement any local refinement. Also, a large number of grid points are needed for FDTD methods to have numerical solutions of high-order accuracy. For instance, if we use a second-order finite difference method to solve a wave equation, more than 30 grid points per wavelength may be needed to have a numerical solution with 1% error [33]. For CROW applications in optical buffers and delay lines, the propagation speed, the arrival time and the phase of signals are critical parameters to have. Thus, a high-order numerical technique is crucial for any numerical simulations so that we can gauge the influence of the resonator size disorder to the optical coupling strength and light propagation efficiency. It has been well known [4,10] that spectral element methods, in which approximations are based on highorder orthogonal polynomials such as Legendre or Chebyshev polynomials, have the least numerical phase error among all numerical methods including the widely-used Yee’s scheme [32] and other FDTD methods [24]. Therefore, in our simulation we will employ a high-order DSEM for solving Maxwell’s equations that govern optical coupling and light propagation via WGMs in microcylinder CROW structures. In contrast to those previous numerical studies by finite difference methods [11] or the method of matrix analysis [18] for related devices, the DSEM will provide higher-order accuracy and therefore the expected numerical simulation results will provide more valuable insight for designing and building general practical CROW structures. However, considering that all elements of this scheme have been well known and thoroughly documented elsewhere [12–14,20,21,25,27,28], the detail discussion of the scheme shall be omitted. To solve Maxwell’s equations (7) in a general two-dimensional domain Ω by the DSEM, the physical domain is first divided into non-overlapping quadrilateral and/or triangular physical elements. Each physical element is then mapped onto a reference element, either the square R0 = [−1, 1] × [−1, 1] or the triangle T0 = {(x, y) | 0  x, y; x + y  1}, by an isoparametric transformation x = χ(ξ ), where x = (x, y) and ξ = (ξ, η) are the coordinates in the physical element and the reference element, respectively. Then on the reference element, the solution is approximated by a linear combination of the basis functions of some approximation space, and the residual of the approximation is required to be orthogonal to the approximation space locally within each element, yielding a system of ordinary differential equations for the time-dependent expansion coefficients, which can be solved by explicit schemes. 2.4. Orthogonal bases on rectangles and triangles An important issue in the implementation of a DSEM is the choice of the approximation basis functions. For numerical stability and accuracy concerns when we solve the resulting system of ordinary differential equations, the local mass matrices should be made as simple as possible. Because the order of the basis functions may have to be taken large to obtain desired accuracy, in the framework of spectral element methods, it is critical to have local mass matrices that are easy to invert and well-conditioned. The ideal case will be that the local mass matrices are diagonal, which can be achieved by employing orthogonal basis. In our simulation, we use both quadrilateral and triangular elements to partition the underlying domain. For quadrilateral elements, the approximation space on the standard reference square R0 is chosen as PM,M (R0 ) = PM × PM , where PM represents the space of polynomials of degree M or less. An orthogonal nodal basis for PM,M (R0 ) is the set of tensor products of the Lagrange interpolation polynomials βmn (ξ, η) = φm (ξ )φn (η),

0  m, n  M,

(8)

where φi (ξ ) =

M  j =0, j =i

ξ − τj τi − τj

with the nodal points τi , i = 0, 1, . . . , M, being the Gauss points in the interval [−1, 1]. For triangular elements, the approximation space on the standard reference triangle T0 is frequently chosen as PM (T0 ) = span{ξ m ηn , 0  m, n, m + n  M}.

(9)

S. Deng / Applied Numerical Mathematics 57 (2007) 475–485

479

An orthogonal polynomial basis for PM (T0 ) can then be obtained by using transformations between triangles and squares and warped tensor product grids within triangles [8]   2ξ − 1 (1 − η)m Pn2m+1,0 (2η − 1), 0  m, n, m + n  M, gmn (ξ, η) = 2m Pm0,0 (10) 1−η α,β

where Pn (x) represents the nth order Jacobi polynomials on [−1, 1] which are orthogonal polynomials under Jacobi weight w(x) = (1 − x)α (1 + x)β , i.e., 1

α,β

(1 − x)α (1 + x)β Pl

(x)Pmα,β (x) dx = δlm .

−1

However, the orthogonal polynomial basis on the triangle T0 is not a nodal basis, which makes its implementation much more inconvenient than that of a nodal basis. Additionally, for the above orthogonal polynomial basis a warped tensor product grid has to be employed for the accurate approximation of integrals, but that grid is over-sampled since it requires twice as many grid points as there are degrees of freedom in the polynomial expansion. Therefore, we use the following orthogonal non-polynomial nodal basis on the triangle T0 [6]   2ξ ψmn (ξ, η) = φm − 1 φn (2η − 1), 0  m, n  M. (11) 1−η As implied in Eq. (11), this basis is constructed by the standard tensor product of the Lagrange interpolation polynomials on the interval [−1, 1] to form a basis on the square R0 , which is then transformed to a basis on the triangle T0 by the same “collapsing” mapping as that used in [8] to construct the orthogonal polynomial basis (10). Unlike the orthogonal polynomial basis (10), however, this non-polynomial basis is a nodal basis over the following collocation points on the triangle T0   (1 + τi )(1 − τj ) 1 + τj , , 0  i, j  M, (ξi , ηj ) = 4 2 namely, ψmn (ξi , ηj ) = δmi δnj , 0  i, j  M. Also, this non-polynomial nodal basis on the triangle T0 is orthogonal in the Legendre inner product defined by  ωm ωn (1 − τn ) δmi δnj . (ψmn , ψij ) = ψmn (ξ, η)ψij (ξ, η) dξ dη = 8 T0

Additionally, the approximation space characterized by the above orthogonal non-polynomial nodal basis is found to contain the space PM (T0 ). More precisely, let GM (T0 ) denote the finite element space on the triangle T0 spanned by the orthogonal non-polynomial nodal basis, i.e.,   GM (T0 ) = span ψmn (ξ, η), 0  m, n  M . Then PM (T0 ) ⊂ GM (T0 ). And actually a very interesting representation of the finite element space GM (T0 ) is given by [6]   ξt ξ i ηj , 0  i, j, i + j  M; , 1  s  M, s  t  M . (1 − η)s For simplicity of programming, we take the same order for both the approximation space G(T0 ) for triangles and the approximation space PM,M (R0 ) for rectangles. Mortars can be used to couple element faces along which the orders differ though [14]. 3. Numerical results and discussions In this section, we shall apply the DSEM to study optical coupling and light propagation via WGMs in microcylinder CROW structures with slight size disorder. But, first we should point out that the exponential convergence of the DSEM with respect to the expansion order M has been previously demonstrated in [7], and therefore will not be included in this paper.

480

S. Deng / Applied Numerical Mathematics 57 (2007) 475–485

Fig. 1. A typical computational mesh for three coupled microcylinders with slight size disorder.

Fig. 2. (a) The mapping between the reference square R0 and a curved quadrilateral element; (b) The mapping between the reference triangle T0 and a curvilinear triangle ABC.

3.1. Mappings for curved elements As mentioned earlier, both quadrilateral and triangular elements are used in our simulation. Fig. 1 shows a typical computational mesh for three coupled microcylinders with slight size disorder. As described in Section 2.3, each physical element in the subdivided domain shall be mapped individually onto the reference square R0 or the reference triangle T0 by an isoparametric transformation, which in general can be arrived by the blending function method [23]. Next we only give the transformation required to map those curvilinear elements shown in Fig. 1 onto the corresponding reference element R0 or T0 . For more detailed description of the blending function method as well as more examples the readers may consult [23]. First, we will need a mapping from the reference square R0 onto a quadrilateral element with one curved boundary. Suppose that the side BC is curved as shown in Fig. 2(a), and it is parametrizable by (χ(η), π(η)) such that χ(−1) = x2 , χ(1) = x3 , π(−1) = y2 , and π(1) = y3 . Then we can write 1 1 1+ξ x = (1 − ξ )(1 − η)x1 + (1 − ξ )(1 + η)x4 + χ(η), 4 4 2 1 1+ξ 1 π(η). y = (1 − ξ )(1 − η)y1 + (1 − ξ )(1 + η)y4 + 4 4 2 Each triangular element in Fig. 1 has also one curved boundary so we need a transformation that can map the reference triangle T0 onto a triangle with one curved boundary. Suppose that the curved side AB is parametrizable as shown in Fig. 2(b) by (χ(ξ ), π(ξ )) such that χ(0) = x1 , χ(1) = x2 , π(0) = y1 , and π(1) = y2 . Then we can write 1 − ξ − η  , x = (1 − ξ − η)x1 + ξ x2 + ηx3 + χ(ξ ) − (1 − ξ )x1 − ξ x2 1−ξ 1 − ξ − η  y = (1 − ξ − η)y1 + ξy2 + ηy3 + π(ξ ) − (1 − ξ )y1 − ξy2 . 1−ξ 3.2. Boundary and initial conditions Appropriate boundary conditions shall be posted at the artificial numerical boundaries of the computational domain to prevent outgoing waves from reflecting off the boundaries. Otherwise, the reflections may falsify the computational results and cause instabilities. In practical simulations, absorbing boundary conditions (ABC) such as perfectly

S. Deng / Applied Numerical Mathematics 57 (2007) 475–485

481

matched layer (PML) boundary conditions [1,3] are widely-used. Since most of the electromagnetic energy of a WGM is confined inside the cylinder and fields decay rapidly away from the cylindrical boundary, a simple matched layer (ML) technique [30] shall be sufficient. To investigate optical coupling and light propagation in CROW devices, we assume that initially there exists a WGM in the most left resonator and no fields exist inside all other cylinders. More specifically, the exact values of a WGM in the most left cylinder are taken as the initial conditions in the entire computational domain except for the inside of all other cylinders, where a zero field is initialized. In addition, to assure that the initial field satisfies the interface condition on the cylinder surfaces, in our numerical test we also assume that over the cylinder surfaces exist surface currents of the form Jm (x, t) = J0m (x)e−αt ,

Je (x, t) = J0e (x)e−αt ,

(12)

where the constant α > 0 is chosen so that the surface currents become negligible after a short period of time, and J0m and J0e are calculated from the initial fields E(x, 0) and H(x, 0) as     J0e (x) = n × H+ (x, 0) − H− (x, 0) . J0m (x) = n × E+ (x, 0) − E− (x, 0) , For such boundary currents, the numerical normal flux for the dielectric surfaces in the DSEM can be written as [7,15]   − + E+n×H)+ −Je n × (Y E−n×H) Y+(Y − Y −Y+Y + Jm − +Y + − (F · n) = , − + +J + m −n × (ZH+n×E) Z+(ZH−n×E) + Z −Z+Z + Je − +Z +

for the − side of the surface and   − − E+n×H)+ −Je n × (Y E−n×H) Y+(Y + Y −Y+Y + Jm − +Y + + (F · n) = , − + +J − m −n × (ZH+n×E) Z+(ZH−n×E) − Z −Z+Z + Je − +Z +

± ± for the + side of the surface, where n is the outward unit normal to the surface,  and Z and Y are the local impedance ± ± ± ± and the local admittance, respectively, and are defined as Z = 1/Y = μ / .

3.3. Optical coupling between two size mismatched cylinders Successful optical coupling via WGMs between two identical microcylinders has been numerically demonstrated [7]. For example, Fig. 3 illustrates the generation of a clockwise WGM of azimuthal number 8 in the right cylinder due to optical coupling and phase matching when a WGM of the same azimuthal number initially exists in the left cylinder. However, in making real physical CROW devices, we will have to consider many factors such as inter-resonator separation and different resonator sizes. It has been found that as the separation between two resonators increases, the optical coupling strength always appears to decrease. This monotonic dependence can easily be understood by the fact that the spatial overlap of the evanescent WGMs of the two resonators will reduce as the inter-resonator separation increases. To assess the effects of the size disorder on the optical coupling, we consider a model of two microcylinders either in contact or separated with a gap of 5% of the radius of the left cylinder. We then gradually increase the size of the other cylinder. Fig. 4 shows the computed coupling coefficient κ as a function of L/λ, the ratio of the perimeter L of the right cylinder to the azimuthal wavelength λ of the WGM in the left cylinder. The coupling coefficient κ is defined as the percentage of energy coupled into the right cylinder from the left cylinder after 20 rotations of the WGM in the left cylinder, where the electromagnetic energy stored in an electromagnetic field in a volume V is defined as    2 E μH 2 dv. (13) + P= 2 2 V

Now it is clear that the dependence of optical coupling on size disorder is not monotonic. For instance, for the case the two cylinders are separated, when the size of the right cylinder is 2% larger than that of the left one (i.e., L/λ = 8.16), around 40% of the total energy is transmitted from the left cylinder to the right indicating strong coupling. Then, we increase the size disorder of the right cylinder to 6% (i.e., L/λ = 8.48), and this time only about 10% of the total energy is transmitted. However, if we further increase the size of the right cylinder to 10% larger than that of the left

482

S. Deng / Applied Numerical Mathematics 57 (2007) 475–485

Fig. 3. Optical energy transport via WGMs between two identical microcylinders in contact. The initial state of the system is represented by a counterclockwise WGM of azimuthal number 8 in the left cylinder. The four sequential snapshots at the times t = 2, 6, 8 and 10 illustrate the generation of a clockwise WGM of the same azimuthal number in the right cylinder due to resonant optical coupling and phase matching. The expansion order M = 10 is used in this simulation.

Fig. 4. Percentage of energy coupled into the right cylinder from the left cylinder as a function of L/λ, the ratio of the perimeter L of the right cylinder to the azimuthal wavelength λ of the WGM in the left cylinder. The two cylinders are either in contact or separated by w = 0.05r1 , and the energy is measured at t = 20. The expansion order M = 10 is used in this simulation.

cylinder (i.e., L/λ = 8.8), the energy transport actually improves dramatically to the 2% case, which again has around 40% of the total energy coupled into the right cylinder. This optical coupling and transport phenomena appears to be, but actually is not, surprising because it can be appreciated in terms of the ratio of the perimeter of the right cylinder to the azimuthal wavelength of the WGM in the left cylinder. Once the perimeter of the right cylinder can fit a whole number of azimuthal waves of the WGM, then the strongest optical coupling and transport occurs.

S. Deng / Applied Numerical Mathematics 57 (2007) 475–485

483

3.4. Light propagation via WGMs in microcylinder CROWs with slight size disorder Next we shall simulate light propagation in CROW structures with slight size disorder in order to provide some guidelines to the size disorder tolerance in building real physical CROW devices. Since only “large” systems, say those with 10 or more than 10 resonators, are regarded as the new type of optical waveguide, we shall consider chains of cylindrical dielectric microresonators with random but slight size disorder. Again, successful light propagation via WGMs has been experimentally observed [9], and numerically demonstrated in perfect CROW devices [5]. For example, Fig. 5(a) shows the temporal history of the energy distribution in individual microcylinders during the light propagation through a 10 cylinder long perfectly ordered CROW structure. As indicated, as much as 85% of the electromagnetic energy in the system can be transmitted all the way through the chain to the most right cylinder at around t = 4200 fs. Also, as indicated by those equally spaced energy peaks in

Fig. 5. Temporal history of energy distribution in individual microcylinders during the light propagation via WGMs through a 10 cylinder long CROW structure. The radius of the most left cylinder is fixed at r1 = 2.283 µm, but the sizes of the other cylinders vary randomly between r1 and 1.02r1 . The initial WGM in the most left cylinder has azimuthal number 10 and time frequency 100 THz. The energies are normalized to the electromagnetic energy of the initial WGM in the most left cylinder. The expansion order M = 8 is used in this simulation. (a) There is no size disorder; (b) The maximal random size disorder is 0.5%; (c) The maximal random size disorder is 1%; and (d) The maximal random size disorder is 2%.

484

S. Deng / Applied Numerical Mathematics 57 (2007) 475–485

individual cylinders before t = 4200 fs, the light appears to travel through the chain at a uniform rate. (Our system does not account for input/output coupling at the both ends of the waveguide. And that is for this reason that, the most right cylinder can accumulate as much as 85% but not 100% of the total electromagnetic energy of the system.) Next we consider different disordered chains with random size dispersion varying in 0–2% range, that is, the cylinder sizes vary from r1 to 1.02r1 , where r1 again represents the radius of the most left cylinder. Figs. 5(b), (c), and (d) again show the temporal energy history during the light propagation through these systems. For the case (b) where the maximal random size dispersion is only 0.5%, optical transport through the chain is observed, but as illustrated in Fig. 5(b), the propagation does not exhibit the same uniform pattern as in otherwise perfectly ordered chains, and the overall optical transport efficiency appears to be much lower. For example, at around t = 5700 fs, only about 15% of the total electromagnetic energy in the system is accumulated in the right end of the system. For the system with the maximal size dispersion being 1%, as shown in Fig. 5(c), although a small part of the electromagnetic energy is transported through the chain, most of the energy seems to bounce forth and back between the first two microcylinders. And finally, for the case (d) whose maximal size dispersion is 2%, as suggested in Fig. 5(d), the optical transport through the chain looks stronger than the 1% size dispersion case. In contrast to perfectly ordered CROW structures, however, the steady, uniform propagation phenomena is again not observed here. It should be reminded that our simulations employ only a single WGM with a very small azimuthal number and our systems do not account for input/output coupling. The actual performance of real physical CROW structures may also depend on many other factors, including the actual distribution of the size dispersions through the chains. For instance, successful optical transport through a disordered chain of polystyrene microspheres with 1% size dispersion is actually observed in the lab [2]. Nevertheless, the above results clearly indicate the high sensitivity of optical transport through a CROW to its size disorder. 4. Conclusions In this paper, we have applied a DSEM for the study of optical coupling and light propagation via WGMs through microcylinder CROW structures with slight size disorders. We have demonstrated the optical coupling strength via WGMs between two size mismatched microcylinders does not depend monotonically on the size dispersion due to phase matching. We have also demonstrated the long-range propagation effects in slightly disordered chains of dielectric microcylinders, which indicates the feasibility of developing and building optical waveguides based on using high quality resonators even with some potential size disorder. On the other hand, the optical transport via WGMs through such chains is highly sensitive to their internal size dispersion due to the evanescent nature of the WGMs. Acknowledgement This work was supported, in part, by funds provided by the University of North Carolina at Charlotte. References [1] S. Abarbanel, D. Gottlieb, On the construction and analysis of absorbing layers in CEM, Appl. Numer. Math. 27 (1998) 331–340. [2] V.N. Astratov, J.P. Franchak, S.P. Ashili, Optical coupling and transport phenomena in chains of spherical dielectric microresonators with size disorder, Appl. Phys. Lett. 85 (2004) 5508–5510. [3] J. Berenger, A perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys. 114 (1994) 185–200. [4] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods in Fluid Dynamics, Springer-Verlag, New York, 1987. [5] S. Deng, W. Cai, Numerical study of light propagation via whispering gallery modes in microcylinder coupled resonator optical waveguides, Opt. Express 12 (2004) 6468–6480. [6] S. Deng, W. Cai, Analysis and application of an orthogonal nodal basis on triangles for discontinuous spectral element methods, Appl. Numer. Anal. Comp. Math. 2 (2005) 326–345. [7] S. Deng, W. Cai, Discontinuous spectral element method modeling of optical coupling by whispering gallery modes between microcylinders, J. Opt. Soc. Am. A 22 (2005) 952–960. [8] M. Dubiner, Spectral methods on triangles and other domains, J. Sci. Comput. 6 (1991) 345–390. [9] H. Furukawa, K. Tenjimbayashi, Light propagation in periodic microcavities, Appl. Phys. Lett. 80 (2002) 192–194. [10] D. Gottlieb, S. Orszag, Numerical Analysis of Spectral Methods: Theory and Applications, Society for Industrial and Applied Mathematics, Philadelphia, PA, 1977.

S. Deng / Applied Numerical Mathematics 57 (2007) 475–485

485

[11] S.C. Hagness, D. Rafizadeh, S.T. Ho, A. Taflove, FDTD microcavity simulations: Design and experimental realization of waveguide-coupled single-mode ring and whispering-gallery-mode disk resonators, J. Lightwave Technol. 15 (1997) 2154–2165. [12] J.S. Hesthaven, T. Warburton, High order nodal methods on unstructured grids I, Time domain solutions of Maxwell’s equations, J. Comput. Phys. 181 (2002) 186–221. [13] D. Komatitsch, R. Martin, J. Tromp, M.A. Taylor, B.A. Wingate, Wave propagation in 2-D elastic media using a spectral element method with triangles and quadrangles, Comput. Acoustics 9 (2001) 703–718. [14] D.A. Kopriva, S.L. Woodruff, M.Y. Hussaini, Computation of electromagnetic scattering with a non-conforming discontinuous spectral element method, Int. J. Numer. Meth. Engng. 53 (2002) 105–122. [15] A.H. Mohammadian, V. Shankar, W.F. Hall, Computation of electromagnetic scattering and radiation using a time domain finite volume discretization procedure, Comput. Phys. Comm. 68 (1991) 175–196. [16] S. Mookherjea, A. Yariv, Coupled resonator optical waveguides, IEEE J. Sel. Topics in Q. El. 8 (2002) 448–456. [17] S. Mookherjea, A. Yariv, Kerr-stabilized super-resonant modes in coupled-resonator optical waveguides, Phys. Rev. E 66 (2002) 046610. [18] J. Poon, J. Scheuer, S. Mookherjea, G.T. Paloczi, Y. Huang, A. Yariv, Matrix analysis of microring coupled-resonator optical waveguides, Opt. Express 12 (2004) 90–103. [19] L. Rayleigh, Further applications of Bessel functions of high order to the whispering gallery and allied problems, Phil. Mag. 27 (1914) 100–104. [20] S.J. Sherwin, G.E. Karniadakis, A new triangular and tetrahedral basis for high-order (hp) finite element methods, Int. J. Numer. Meth. Engng. 38 (1995) 3775–3802. [21] S.J. Sherwin, G.E. Karniadakis, A triangular spectral element method: Applications to the incompressible Navier–Stokes equations, Comput. Meth. Appl. Mech. Engrg. 123 (1995) 189–229. [22] J.A. Stratton, Electromagnetic Theory, McGraw-Hill, New York, 1941. [23] B. Szabo, I. Babuska, Finite Element Analysis, Wiley, New York, 1991. [24] A. Taflove, S.C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, Artech House, Boston, MA, 2000. [25] M.A. Taylor, B.A. Wingate, A generalized diagonal mass matrix spectral element method for non-quadrilateral elements, Appl. Numer. Math. 33 (2000) 259–265. [26] J.R. Wait, Electromagnetic whispering gallery modes in a dielectric rod, Radio Sci. 2 (1967) 1005–1017. [27] T. Warburton, Application of the discontinuous Galerkin method to Maxwell’s equations using unstructured polymorphic hp-finite elements, in: B. Cockburn, G. Karniadakis, C.-W. Shu (Eds.), Discontinuous Galerkin Methods: Theory, Computation and Applications, Springer, New York, 2000, p. 451. [28] B. Wingate, J.P. Boyd, Triangular spectral element methods for geophysical fluid dynamics applications, in: A.V. Ilin, L.R. Scott (Eds.), Proceedings of the 3rd International Conference on Spectral and Higher Order Methods, Houston, TX, 1990, Houston J. Mathematics, Houston, 1996, pp. 305–314. [29] Y. Xu, R.K. Lee, A. Yariv, Propagation and second-harmonic generation of electromagnetic waves in a coupled-resonator optical waveguide, J. Opt. Soc. Am. B 17 (2000) 387–400. [30] B. Yang, D. Gottlieb, J.S. Hesthaven, Spectral simulation of electromagnetic wave scattering, J. Comput. Phys. 34 (1997) 216–230. [31] A. Yariv, Y. Xu, R.K. Lee, A. Scherer, Coupled-resonator optical waveguide: A proposal and analysis, Opt. Express 24 (1999) 711–713. [32] K.S. Yee, Numerical solution of initial boundary value problems in solving Maxwell’s equations in isotropic media, IEEE Trans. Antennas Propag. P-14 (1966) 302–307. [33] A. Yefet, P.G. Petropoulos, A non-dissipative staggered fourth-order accurate explicit finite difference scheme for the time-domain Maxwell’s equations, ICASE Technical Report, No. 99-30, NASA/CR-1999-209514, 1999.