Supermodes in multiple coupled photonic crystal waveguides

We analyze the supermodes in multiple coupled photonic crystal waveguides for long-wavelengths. In the tight-binding limit we obtain analytic results that agree with fully numerical calculations. We find that when the field flips sign after a single photonic crystal period, and there is an odd number of periods between adjacent waveguides, the supermode order is reversed, compared to that in conventional coupled waveguides, generalizing earlier results obtained for two coupled waveguides. © 2006 Optical Society of America OCIS codes: (130.2790) Guided waves; (230.7370) Waveguides; (050.1960) Diffraction theory References and links 1. See, e.g., Y. Sugimoto, Y. Tanaka, N. Ikeda, Y. Nakamura, K. Asakawa, and K. Inoue, “Low propagation loss of 0.76 dB/mm in GaAs-based single-line-defect two-dimensional photonic crystal slab waveguides up to 1 cm in length” Opt. Express 12, 1090-1096 (2004), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-6-1090 2. A. Martinez, F. Cuesta, and J. Marti, “Ultrashort 2-D photonic crystal directional couplers,” IEEE Photonics Technol. Lett. 15, 694-696 (2003). 3. C.M. de Sterke, L.C. Botten, A.A. Asatryan, T.P. White, and R.C. McPhedran, “Modes of coupled photonic crystal waveguides,” Opt. Lett. 29, 1384-1386 (2004). 4. P. Yeh, Optical Waves in Layered Media (Wiley, Hoboken, 1988). 5. L. D. Landau and E.M. Lifshitz, Quantum Mechanics (Non-Relativistic Theory), 3rd Ed. (Pergamon, Oxford, 1977). 6. L.C. Botten, N.A. Nicorovici, R.C. McPhedran, C.M. de Sterke, and A.A. Asatryan, “Photonic bandstructure calculations using scattering matrices,” Phys. Rev. E 64, 046603:1-18 (2001). 7. L. C. Botten, N. A. Nicorovici, A. A. Asatryan, R. C. McPhedran, C. M. de Sterke, and P. A. Robinson, “Electromagnetic scattering and propagation through grating stacks of metallic and dielectric cylinders for photonic crystal calculation. Part 1: Formulation”, J. Opt. Soc. Am. 17, 2165-76 (2000). 8. M. Born and E. Wolf, Principles of Optics 6th Edition (Pergamon, Oxford, 1980), p. 66. 9. L. C. Botten, T. P. White, A. A. Asatryan, T. N. Langtry, C. M. de Sterke and R. C. McPhedran, “Bloch mode scattering matrix methods for modelling extended photonic crystal structures. Part I: Theory”, Phys. Rev. E 70, 056606, 1-13 (2004). 10. D. Felbacq, A. Moreau, and Rafik Smaâli, “Goos-Hänchen effect in the gaps of photonic crystals,” Opt. Lett. 28, 1633-1635 (2003). 11. C.M. de Sterke, “Superstructure gratings in the tight-binding approximation,” Phys. Rev. E 57 3502-3509 (1998). 12. M. Bayindir, B. Temelkuran, and E. Özbay, “Tight-binding description of the coupled defect modes in threedimensional photonic crystals,” Phys. Rev. Lett. 84, 2140 (2000). 13. S. Mookherjea and A. Yariv, “Optical pulse propagation in the tight-binding approximation,” Opt. Express 9, 91-96 (2001). 14. N.W. Ashcroft and N.D. Mermin Solid state physics, (Saunders College, Philadelphia, 1976). #9886 $15.00 USD Received 6 December 2005; revised 28 December 2005; accepted 3 January 2006 (C) 2006 OSA 9 January 2006 / Vol. 14, No. 1 / OPTICS EXPRESS 387


Introduction
One of the most studied applications of photonic crystals (PCs) are photonic crystal waveguides (PCWs).These are formed by creating a line defect in a PC and by operating at frequencies inside a bandgap.At such frequencies, the light cannot propagate in the photonic crystal and is therefore forced to remain in the defect, which thus acts as a waveguide.Though the losses of PCWs were initially very high, values below 1 dB/cm were reported in recent work [1].The key difference compared to conventional waveguides is that confinement is provided by Bragg reflection in the PC, rather than by total internal reflection.
Coupled waveguides, waveguides that run parallel to each other and are sufficiently close so that their evanescent tails overlap, have applications such as directional couplers.They are particularly interesting in a PC geometry since the beat length between the supermodes of the structure, the size of which determines the length of the device, can be made very short [2].These structures were shown to have some properties that differ from conventional waveguides [3].For example, the fundamental mode of two coupled PCWs can be odd, whereas the fundamental mode of a pair of coupled conventional waveguide is always even [4].The fundamental mode was found to be odd when the electric field in the barrier between the waveguides changes sign after a period and when the number of periods in this barrier is odd [3].Such a situation can never occur in conventional planar structures since in the barriers the field decays monotonically.Nevertheless, it was also shown that the oscillation theorem [5] applies in coupled PCWs: consecutive modes have one additional electric field node, and thus the odd fundamental mode has one fewer node than the second mode [3].Briefly, the oscillation theorem is consistent with the fundamental mode being odd, because we are interested only in bound modes, that is, modes that grow or decay exponentially in the PC cladding.In contrast, the oscillation theorem applies to all modes, both bound and unbound [3].
Here we generalize the previously developed theory to an arbitrary number of identical coupled PCWs and consider this geometry at long-wavelengths where the analogy to conventional waveguides is strongest [3].We then concentrate on the limit in which the guides are sufficiently spaced so that their mutual coupling is weak.In this tight-binding limit, analytic results can be obtained that are the PC equivalent of well-known expressions in the literature [4].We develop this theory in Section 2 and present the results, both in terms of the propagation constants and in terms of the associated fields, in Section 3. Section 4 contains a discussion of our results and concludes.

Theory
We consider the two-dimensional structure in Fig. 1 and find its modes by considering the multiple reflection and transmission of fields in a systematic, recursive manner.We regard the bulk photonic crystal as a sequence of grating layers, of period d, arranged in square lattice.With β denoting the propagation constant of a mode, the fields satisfy the quasiperiodicity condition V (r r r + de e e x ) = V (r r r) exp(iβ d), with e e e x denoting the unit vector pointing in the waveguide direction.In each guide j, it is natural to express the fields in an expansion of upward (+) and downward (−) propagating plane waves, the direction of which are related by the grating equation, i.e. [6] where Inside the photonic crystal, we expand the field in its Bloch mode basis.Representing the fields in the gaps between the layers of the bulk crystal by plane wave expansions of the form (1), with coefficients, T , we apply Bloch's theorem to derive a transfer matrix equation T f f f = μ f f f .Here, the transfer matrix T may be calculated using the reflection and transmission scattering matrices of a single grating layer (c.f.Eq. ( 6) below).The set of Bloch modes that arise from the solution of the eigenvalue problem may be partitioned into sets of forward and backward propagating states, with the forward modes associated with a matrix of eigenvalues Λ Λ Λ = diag μ j (with |μ j | = 1 for propagating modes and |μ i | < 1 for evanescent modes) and eigenvectors For a square lattice, the corresponding set of backward modes have Bloch factors Λ Λ Λ −1 and eigenvectors The two key quantities to emerge from this eigenproblem are thus the matrix of Bloch factors Λ Λ Λ, which determines the propagation of modes in the bulk crystal, and the scattering matrix R R R ∞ = F F F + (F F F − ) −1 which characterizes the reflection of a field incident on a semi-infinite bulk crystal from free space.Here, F F F ∓ are matrices whose columns are the eigenvector components ff f ∓ , with the matrices F F F − and F F F −1 − respectively performing changes of basis from plane waves to Bloch modes and vice versa [6].
From the PC Bloch mode quantities, we may then write down expressions for the reflection and transmission scattering matrices of each of the l-layer barriers resembling the corresponding scalar forms for a Fabry-Perot etalon.In Equations ( 2) and (3), quantifies the propagation of the Bloch modes within the l-layer barriers.
The geometry of Fig. 1 shows that the structure can be regarded as a composite of (m − 1) barriers that are sandwiched between two semi-infinite PCs.The presence of each of the guides of width h is taken into account, by "padding, " relative to the field origins located in the center of the guides, the relevant scattering matrices to allow for propagation to and from the phase origin.The fields in the uppermost and lowermost guides thus satisfy with the leftmost expressions derived from the end-mirror conditions and with the rightmost equations derived from the propagation of plane waves through the composite (m − 1) barrier structure.In Equations ( 4) and ( 5), R R R ∞ = P P P 1/2 R R R ∞ P P P 1/2 , with P P P = diag[exp(iχ p h)] representing, for each plane wave order, the phase accumulation when traversing a waveguide of width h.
While the scattering matrices, R R R n and T T T n , for the n = m − 1 "padded" barriers can be derived by recursion [7] from those for a single padded barrier, R R R 1 = P P P 1/2 R R R l P P P 1/2 and T T T 1 = P P P 1/2 T T T l P P P 1/2 , we will follow an alternative route that involves the calculation of the Bloch modes of an infinite structure composed of these padded barrier sections.This leads us to the dispersion equation in a form that is suitable for computational use, and which is also amenable to an asymptotic analysis in the long wavelength and tight binding limits.To do so, we introduce the transfer matrix T for crossing a padded layer.Proceeding as before, we may solve the eigenvalue problem for this transfer matrix, deduce the corresponding Bloch factors Λ Λ Λ and eigenvector matrices F F F ∓ , and compute the reflection and transmission matrices R R R n and T T T n in a form analogous to Equations ( 2) and (3).We next introduce the transfer matrix to cross the n = m − 1 padded barriers, namely and note that Because of the up-down symmetry of the layers, and of the entire structure, the transfer matrix T satisfies the symmetry relationship T −1 = U T U , where U is the matrix that reverses the partition order.Accordingly, some straightforward manipulation reveals that and since this must have the form of the transfer matrix equation ( 7), we deduce that the modes must be either symmetric (σ = 1) or antisymmetric (σ = −1) according to g g g ± m = σ g g g ∓ 1 .The dispersion equation is then derived by substituting the end mirror constraints, and the symmetry relations g g g ± m = σ g g g ∓ 1 into the transfer matrix equation (7).This leads to in which g g g = g g g + 1 .Each of the two equations ( 9) are equivalent and reduce to subject to T T T n being invertible.The modes of the structure may now be found by solving det M M M(k, β ) = 0 for the propagation constant β , after which the mode may be reconstructed from the null vector g g g.
At long wavelengths, all fields can be well approximated by expansions comprising only a single propagating plane wave term, thus reducing all matrix quantities to scalars and facilitating an asymptotic analysis in the tight binding limit.The transfer matrix then becomes a 2 × 2 matrix with unit determinant, for which simple formulae exist [8] for the elements of its n th power, namely in which u n (t) = sin[(n + 1)φ ]/ sinφ denotes the Chebyshev polynomial of the second kind of degree n in the variable t = cos φ , and Furthermore, R l and T l are given respectively by the scalar versions of Equations ( 2) and (3), in which we have replaced the matrix R R R ∞ by its specular (0, 0) th order component ρ, and the matrix Q Q Q by the single, dominant Bloch factor μ for the bulk crystal, a number which is real and has a magnitude |μ| < 1 (since we are operating in a band gap).
With this simplification, we now consider the tight binding limit in which the number of barriers layers l is sufficiently large that ξ = μ l is small.Guided by the scalar form of the dispersion equation for a single waveguide, Pρ = ±1, in which both P = exp(iχh) and ρ have unit magnitude, the latter because we are operating in a band gap, we seek a solution of the form where v is a small argument that represents the effect on the phase upon reflection off a waveguide-cladding interface due to presence of all other waveguides.While we may expand ν(ξ ) in a power series and seek a self-consistent perturbation solution to arbitrary order, we need to work only to first order with ν(ξ ) = αξ to obtain the results that we are after.In what follows, we work with the form Pρ = exp[iν(ξ )] corresponding to the fundamental symmetric mode for each of the single waveguides.The ansatz Pρ = − exp[iν(ξ )] gives essentially the same result and it is therefore not considered further here.
Then, returning to dispersion equations ( 9), we form the corresponding scalar form that is appropriate in the long wavelength limit.We then apply the tight binding approximation and expand each of the factors and variables in Eq. ( 16) in a power series in ξ : Substituting these asymptotic forms into Eq.( 16) and exploiting the recurrence relation that is satisfied by the {u n (t)}, namely u n (t) = 2tu n−1 (t) − u n−2 (t), we derive the tight binding form of the dispersion relation which, correct to first order, is The simultaneous requirement that T T T n is invertible, translates in the long wavelength limit to T n = 0 which, physically, requires the guides be not completely isolated from one another.To first order, in the tight binding limit, this requires that u n−1 (t) = 0. Subject to this constraint we solve (18) to deduce the solutions t s = cos ϑ s , with ϑ s = sπ/(m + 1), for s = 1, 2,...,m, with odd values of s corresponding to the symmetric, or even, modes (σ = 1) and with even values of s corresponding to the antisymmetric, or odd, modes (σ = −1).Here, the integer s enumerates the different supermodes of the structure.
The reduced form of the dispersion equation of the supermodes follows from the argument form of Eq. ( 15), namely χh + arg(ρ) = αξ into which we substitute the possible values of α = α s = 2ℑ(ρ) cosϑ s .Thus, where β s is the propagation constant of mode s and j is an integer.Since the corresponding form for a single waveguide is k 2 − β 2 0 h + arg[ρ(β 0 )] = jπ, we deduce from a simple power series expansion of the left hand side of Eq. ( 19) that Since this indicates that the right-hand side should be evaluated at β s , this leads to an implicit equation that needs to be solved self consistently.However, in the spirit of the tight-binding approximation we evaluate it at β 0 leading to an explicit expression for the β s .
The interpretation of Eq. ( 20) is facilitated by recognizing that the term ∂ argρ(β 0 )/∂ β corresponds to the lateral beam displacement δ x due to the Goos-Hänchen shift [4].The displacement can be associated with a barrier penetration t = −χ 0 δ x/(2β 0 ); note that in the special case of a planar waveguide the parameter t corresponds to the 1/e decay length of the evanescent field in the cladding.With this, the magnitude of the denominator of Eq. ( 20) can be written as (h + 2t)/χ 0 ≡ h eff /χ 0 , where the effective width h eff , a standard parameter in the theory of planar waveguides, can be considered to be the width of the mode that includes the effect of the evanescent tails.Equation (20) can thus be simplified as Having solved the propagation problem and determined the modes of the coupled waveguide system, it remains only to calculate the fields in each of the waveguides which, in turn, lets us reconstruct the field everywhere.For waveguide j, the plane wave field (1) originates from multiple reflections and transmissions through the structure and we write by a consideration of how these fields originate.Solving these, we may derive expressions for g g g ± j in terms of g g g − 1 and g g g + m , from which we calculate the fields in each of the guides, using the null vectors of the dispersion equation (10) for both the symmetric (σ = 1) and antisymmetric (σ = −1) driving fields g g g + 1 = σ g g g − m = g g g.Of particular interest in the long wavelength limit is the total field which is represented by the amplitude of the specular order plane wave Then, for each mode s, we assign g + 1 = σ g − m = 1 according to its symmetry, and perform an asymptotic analysis of the total field (23) in the tight binding limit (|ξ | 1), leading to a result that is consistent with that for conventional waveguides [4].

Results
Here we consider the tight-binding results and compare these with exact results.We first consider the modal propagation constants and after that we then consider the associated eigenfunctions.The results given are for a two-dimensional geometry in which cylinders of radius a and refractive index n = 3 are arranged in a square lattice with period d, where a = 0.3d.We consider W1 waveguides, i.e. waveguides that are formed by the removal of one row of inclusions, so that h = d, and solutions in which the electric field points in the direction parallel to the inclusions.At the normalized frequency d/λ = 0.32, the single mode of a single, isolated waveguide has a propagation constant β 0 d = 1.16223.The dominant Bloch function here has μ = −0.468476,and ρ = exp(−0.52224πi).Since μ < 0, this is the situation in which the fundamental mode of two coupled waveguides is odd when As a first set of results, we show in Figures 2, β 2 s − β 2 0 with increasing thickness of the separating layer l for m = 5 (Fig. 2(a)) and for m = 6 (Fig. 2(b)).The propagation constants are calculated using both an independent numerical method [9], and the full matrix form of the formulation outlined here (10), involving a combination of multipole and Bloch mode methods [6] with 7 plane wave orders, and which for our purposes here can be considered to be exact.Note that the pattern of supermode propagation constants is roughly the same for all l and that, as expected from the μ l factor in Eq. ( 20), the β s values approach β 0 with increasing l.The color of the dots indicates whether the super mode is even (red dots) or odd (blue dots).Note that for m = 5 the fundamental mode, the mode with the highest propagation constant, is always even, whereas for m = 6 the fundamental mode is alternatingly even (when l is even), or odd (when l is odd).We return to this below.A practical issue to be considered when presenting quantitative results for the propagation constants β s of the supermodes is that the μ l in (20) causes them all to converge exponentially to β 0 with large l, thus making the verification of the key asymptotic form (20) difficult.Accordingly, we scale out the exponential term by plotting values of in the case of the tight binding approximation of Equations ( 20) and ( 21).Note that we have included |μ| l , rather than μ l , so that the multiplying factor, and therefore the order in which the modes appear, does not change sign with l.The results are depicted by the dots in Figures 3, with m = 5 in Fig. 3(a) and m = 6 in Fig. 3(b).The dots again indicate exact numerical results, with the red dots corresponding to even supermodes and blue dots corresponding to odd supermodes.Note that the color patterns are the same as in Figures 2. According to Equations ( 20) and ( 21), and plotted in this way, the normalized propagation constants should not depend on l and should be given by t = (−1) l cos ϑ s where ϑ s = sπ/(m + 1) (10).The horizontal lines in Figures 3 show that Eq. ( 20) is an excellent approximation to the propagation constants of the supermodes, particularly when l ≥ 6, since the dots converge to the the results given by Eq. ( 25).In both figures, the dots converge to the asymptotic lines.The quality of the agreement is not surprising since |μ| l < |μ| 6 ≈ 0.01, indicating that ξ 1 as required for the tight-binding approximation to apply.However, even for l = 3, for which |μ| 3 ≈ 0.1, the tight-binding result (20) is fairly accurate, though the need for higher order corrections can be discerned.Fig. 3(b) is similar to Fig. 3(a), but is for a set of m = 6 coupled waveguides and similar conclusions can be drawn as for Fig. 3(a).Note also that the colors of the dots in Figures 3 are consistent with the earlier observation that modes associated with cos ϑ s for odd s are even while those for even s are odd.
We now consider the symmetry of the modes more closely and observe from Figures 2 that for m = 5 the fundamental mode is always even, whereas for m = 6 the fundamental mode alternates between even and odd symmetry according to whether l is odd or even.We find that this applies to all even and odd m, respectively, and is consistent with the result for m = 2 [3].We now consider the eigenmodes in more detail.The field amplitudes in each of the guides j are given by Eq. ( 24), which are consistent with results for conventional waveguides.Figures 4  and 5 show the supermodes for m = 5 and m = 6, respectively, both for l = 6 (Figures 4(a) and 5(a)), and for l = 7 (Figures 4(b) and 5(b)).For l = 6 (Figures 4(a) and 5(a)) the order of the modes is identical to those for conventional structures: in the fundamental supermodes the fields in the waveguides are mutually in phase, whereas in the highest order supermode the fields in adjacent waveguides are out of phase.The key difference with conventional structures occurs when l is odd (and μ < 0, the case we consider here).In this case the order of the supermodes is reversed with respect to that for conventional coupled waveguides: for example, in the highest order supermode the fields in all the waveguides are in phase, whereas in conventional waveguides the field of the fundamental supermode has this property.This is  3).Finally, though we do not show this explicitly here, we have confirmed that the oscillation theorem [5] holds in all cases: the fundamental supermode has the smallest number of nodes, with consecutive higher order modes, having one extra node, also consistent with the earlier results for m = 2.

Discussion and Conclusions
The use of the tight-binding approximation in optics is not new [4,11,12,13].However, the derivation to obtain the desired results typically involves the overlap integral of the field in a single, isolated well or waveguide, with the perturbation associated with neighboring wells or waveguides [4,11,12,13].This type of approach ultimately originates from condensed matter physics, where the tight-binding method was developed to calculate band structures [14].In contrast, here we develop the theory purely in terms of the reflection coefficients within the guiding structure, and it is for this reason that the final results look quite different.It is true that all tight-binding results have a term of the type μ l , which is associated with exponential decay of the field between adjacent wells or waveguides.However, here we are interested in  21) is consistent with the perturbative results for planar waveguides, where this parameter often appears.The origin of the term ℑ(ρ) in Eq. ( 21) can be understood using a simple perturbation treatment which, for simplicity, we outline for a two waveguide structure.In each of the m = 2 guides, we represent the fields by plane wave expansions (1) and write down relationships analogous to Equations ( 4) and ( 5) in which, in scalar form, R∞ = ρP, R1 = R l P, and T1 = T l P. In the tight binding limit, in which we neglect terms of O(ξ 2 ) or higher, we approximate R l ≈ ρ and T l ≈ (1 − ρ 2 )ξ .Appropriate to the case of a symmetric single guide mode, we construct the following total field with v j = g − j + g + j by adding the pairs of Equations ( 4) and ( 5) for g ± j , while for an antisymmetric mode, we would construct a quantity proportional to the field derivative, v j = g − j − g + j .For the symmetric case, then, with v j = g − j + g + j , we have which may be further simplified by respectively relating the driving terms g − 1 and g + 2 to the total field according to g − 1 = ρP/(1 + ρP)v 1 ≈ ρPv 1 /2 and similarly g + 2 ≈ ρPv 2 /2 for the symmetric single mode case, for which Pρ ≈ 1.The final form of the dispersion relation is then derived from the coupled system of equations in which the single guide dispersion equations v j = ρPv j for j = 1, 2 are perturbed by the presence of the neighboring guide through the term −iℑ(ρ)ξ v l for l = 2, 1.It is thus clear that the factor ℑ(ρ) arises directly from the transmission coefficient for a sequence of l layers, and we observe that if ρ = ±1 (as for an electric or magnetic mirror) the perturbation vanishes, in keeping with the guides being completely isolated.We note, in closing, that in the tight binding limit the solution of the dispersion equation may be recast in terms of an algebraic eigenvalue problem derived from ( 27) and (28).The approach generalizes naturally to an arbitrary number of guides, leading to the eigenvalue problem for a tridiagonal matrix (corresponding to nearest neighbor interactions), similar to that derived for a system of conventional waveguides using a coupled mode approach [4].
In conclusion, we have analyzed the supermodes of coupled PCW.We have done so for long wavelengths and in the tight-binding limit, for which simple analytic results can be obtained.They are similar to the modes of coupled conventional waveguides, except that the order of the supermodes is reversed when μ < 0 and l is odd.As a consequence of this, we find that the fundamental supermode is then odd when m is even, consistent with earlier findings [3].

Fig. 1 .
Fig. 1.Schematic of the geometry that we are considering, consisting of m coupled PCWs, separated by m − 1 barriers, all sandwiched between two semi-infinite PCs.Each barrier consists of l PC periods.

Fig. 2 .
Fig. 2. Exact numerical results for β 2 s − β 2 0 with increasing layer thickness l for (a) m = 5, and (b) m = 6 identical waveguides.Parameter values are as given in the text, with the red and blue dots respectively showing values for even and odd symmetric supermodes.

FFig. 3 .
Fig. 3. Comparison of exact numerical results for a set of (a) m = 5, and (b) m = 6 identical waveguides with parameters given in the text.The horizontal lines give the limiting values of F(β s ) in the tight binding approximation, cos[sπ/(m + 1)], whereas the dots give exact results as discussed in the text.The red dots are for even modes while the blue dots indicate odd modes.
illustrated in Figures 4(b) and 5(b).This shows that when m is even, the fundamental mode is odd (Figures 5(b)), consistent with the earlier result for m = 2 [3], and with the results in Figures (2) and (

Fig. 4 .Fig. 5 .
Fig. 4. Fields of the supermodes for m = 5 and (a) l = 6 and (b) l = 7, with the associated propagation constants β s d indicated above the figures.In each case, the fundamental mode is shown on the right, with consecutive higher order modes shown towards the left.The guides are oriented parallel to the x-axis.
The presence of h eff in the denominator of (