## Abstract

Coupled mode theory for waveguide arrays is extended to next-nearest neighbor interactions using propagation equations. Both lateral diffraction and propagation of Floquet-Bloch waves are altered respectively by extra coupling and non-orthogonality between isolated waveguide modes. The analytical formula describing the distortions of the diffraction relation is validated by direct numerical simulation for weakly coupled InP and GaAs shallow ridge waveguides and for strongly coupled Si-SiO_{2} buried strip waveguides. The impact of extended coupled mode theory on propagation and diffraction design in waveguide arrays is discussed with reference to available experimental work.

©2010 Optical Society of America

## 1. Introduction

Light propagation in arrays of weakly coupled waveguides has been studied extensively for fundamental reasons as well as for applications to optical processing devices. Light injected in a single waveguide couples to more and more waveguides as it propagates under the form of Floquet-Bloch (FB) waves. Specific features of these waves have been studied in the linear (discrete diffraction, beam-like behavior, non-divergent beams…) and nonlinear (self-focusing, discrete solitons…) regimes [1,2]. They have also been exploited to demonstrate diffraction engineering [3], Bloch oscillations [4,5], quantum mechanical decay [6], Anderson localization [7,8], Rabi oscillations [9], and optical routing schemes [10–12]. Arrays have then become an important tool of optics and as such require a robust though simple assessment basis, which can treat efficiently in the linear and nonlinear regimes the complex patterns already built or envisioned in the future [14–16].

Waveguide arrays are series of identical equally-spaced (period = S) confining optical structures. Propagation in a single isolated structure is permitted for discrete optical indices corresponding to the guided modes. Coupling between N structures splits those into sets of N sublevels indexed by K_{x}∈[-π/S, π/S] corresponding to “supermodes” encompassing all structures. For an infinite array bands are formed and allow for the free propagation of FB waves across the whole space [3]. Most experimental situations involve mainly the upper band generated around the fundamental level of the isolated structure, all the more when this structure is designed to be monomode. The band shape n(K_{x}) known as the diffraction relation controls all the features of FB waves: direction of rays, divergence of beams [3], refraction and reflection at interfaces [17]... Describing n(K_{x}) is then the key objective of any model.

This basis is at present provided by a long-established approximation, the coupled mode theory (CMT), which describes “weak” coupling situations by assuming (1) the slowly varying envelope scheme and (2) the orthogonality of the modes of the quasi-isolated confining structures even for finite S. CMT predicts a cosine-shaped relation $n\left({K}_{x}\right)={n}_{irw}+2C/{\beta}_{vacuum}\mathrm{cos}\left({K}_{x}S\right)$ where n is the sublevel effective index, n_{irw} is the effective index of the fundamental mode of the isolated waveguide, C is the coupling coefficient between neighboring waveguides, and β_{vacuum}=2π/λ. In addition simple equations describe the propagation of the electric field reduced to the amplitudes of the FB wave in each confining structure. If refined theories give a more accurate description of systems such as pairs of asymmetric waveguides [18,19] they become too cumbersome for the much more complex case of the waveguide array so that CMT is used in most studies because of its very simplicity and in spite of its drastic assumptions.

Using the slowly-varying envelope scheme assumes that C is much smaller than the propagation constant of the upper guided mode β_{irw}. This is satisfied in most systems since C<1mm^{−1} whereas β_{irw}~10 µm^{−1}. On the other hand the orthogonality assumption and deviations from the pure CMT calculations have not been assessed explicitly in waveguide arrays, even though they have an important impact on the behavior of FB waves through deformations of the diffraction relation. Furthermore, the design of functional devices requires as short lengths of active waveguides as possible, i.e. tends to maximize the coupling strength. This trend pushes CMT to the limits of its range of validity, which calls for the check whether CMT holds within a required accuracy for a given waveguide design and for the determination of the necessary corrections to CMT when this is no longer the case.

Only a few papers have tried to assess the validity of CMT in waveguides arrays. Eyges and Wintersteiner [20] have compared sets of supermode levels for arrays of N≤6 circular waveguides obtained by CMT and numerical calculation and they reported large discrepancies (~30%) at very high couplings. At C~0.1mm^{−1} Kaplan and Ruschin [21] found only small discrepancies between output intensities calculated by CMT and observed in experiments. In a recent and stimulating paper Cooper and Mookherjea [22] have renewed the field by developing a simple analytical model which takes non-orthogonality and further-neighbor effects into account to find out the supermodes. With a modal approach they obtain accurate values for the full CMT matrix from numerical data through an inverse-problem calculation and call this approach numerically-assisted CMT.

In this paper, using an analytical model similar to theirs but limited to first and second neighbors, we adopt a direct approach based on propagation equations as is usual in CMT rather than a modal approach. All the parameters – C and the correction factors which express non-orthogonality and second-neighbor effects – can be determined knowing only the mode of the isolated waveguide. This directly yields the corrections to the diffraction relation, which straightforwardly reflects the impact of the parameters. Our extension of CMT to the “moderate” coupling regime is successfully tested against supermode calculations in widely used GaAs, InP, and Si-based shallow-ridge or rectangular waveguide arrays. Scaling the various results demonstrates very general rules which give a clear insight into the corrections to CMT and can help to optimize waveguide arrays. The approach also delimits the “weak” coupling regime unambiguously, which makes it possible to design structures for which CMT applies within a tolerated correction level. We finally discuss the impact of the corrections on experiments and potential devices.

## 2. Analytical modeling: a first-order extension of the coupled mode theory

The derivation of the extended CMT model is given in detail in Appendix, and is based on the following assumptions:

- • propagation takes place along Z; confining structures and mode shapes lie in the XY plane.
- • ε varies slowly on the wavelength scale, except at some particular interfaces where the variations are abrupt and can be taken into account using appropriate boundary conditions. As a result one can substitute the Helmholtz equation to the Maxwell equations.
- • considering a single polarization we search solutions in electric field E which are combinations $E\left(X,Y,Z\right)={\displaystyle \sum _{m}{a}_{m}\left(Z\right)\text{\hspace{0.17em}}{M}_{m}\left(X-mS,Y\right)}\text{\hspace{0.17em}}{e}^{i{\beta}_{irw}Z}$ of the modes of the isolated waveguides M
_{m}(X,Y). The dominant part of the spatial phase variation along Z is accounted for by β_{irw}and the remaining part requires to solve propagation equations for the a_{m}(Z) which contain all the distinctive characteristics of the supermodes. - • under the slowly-varying envelope assumption we neglect second derivatives of a
_{m}in Z with respect to first derivatives. - • isolated waveguides are monomode and lossless and only one of the main polarizations is considered. Then only the first mode m=1 need to be retained in the sum for E(X,Y,Z).
- • all terms involving more-than-second-neighbor confining structures are neglected throughout the calculation. Second-neighbor terms are eliminated only at the end of the calculation in order to assess their perturbation order relatively to the usual CMT equations.

_{x}=K

_{x}S and k

_{z}=(β-β

_{irw})/C propagate across the infinite array according to the diffraction relation

Our model thus yields as desired an extension of the CMT with next-order terms of clear-cut significance in an extended diffraction relation: extra coupling and non-orthogonality between isolated waveguide modes alter respectively lateral diffraction and propagation and arise respectively in the numerator and the denominator. Long range couplings to successive neighbors and the corresponding harmonics of the diffraction relation can be derived exactly using Wannier functions, i.e. a basis of orthogonal functions localized at each site and built from the exact Bloch functions (see e.g [23].). We have checked that, to first order in η and using the diffraction relation given by Eq. (5), Eqs. (3) are equivalent to the set of equations obtained with Wannier functions (see Appendix). Actually the amplitudes a_{m} are not expected to depend strongly on the shape of the localized base functions. It may be stressed that all the parameters entering the diffraction relation can be obtained here from the mode of the isolated waveguide, owing to integrations over either the whole transverse plane or the ridge perturbation section so that the more involved computation of the supermodes is not necessary, in contrast to the Wannier function approach [23] or to Ref [22]. It may also be noted that such a perturbative extension of CMT remained hidden in the modal approach [22] where the mode propagation constants are input data for an inverse problem, although the connection between the numerically-calculated modes and CMT relies upon the same integrals as in our propagative approach.

We now validate our extended CMT model on selected examples by reference to direct numerical solutions of the Maxwell equations. More precisely we compare the predicted diffraction relation with the set of sublevels of the supermodes obtained for a finite number of waveguides in the array, which is actually a sample of this diffraction relation.

## 3. Numerical simulation

#### 3.1 The waveguide arrays used for validation

The validation is performed mainly on prototypal systems, shallow-ridge structures built in the GaAs/GaAlAs and InP/InGaAs systems, widely used in published works on waveguide arrays [24,15] and operated at the telecom wavelength λ=1.55µm; Si/SiO_{2} buried-strip waveguides will also be evaluated in Section 4.1. We shall vary only the ridge width L_{r} and the array period S, which tune respectively the individual ridge waveguide and its coupling to its neighbors. All other parameters will be kept at common literature values (Fig. 1
).

Three reference optical structures are of interest: the lower planar waveguide (lpw) corresponding to the absence of ridges, the upper planar waveguide (upw) corresponding to complete coverage by merged ridges (S=L_{r}), and the isolated ridge waveguide (irw) corresponding to the presence of a single ridge (S=∞). Their effective indices (n_{cladding} < n_{lpw} < n_{irw} < n_{upw} < n_{core}), together with mode maps, are calculated by direct numerical solution of the Maxwell equations with the finite-element method (FEM) implemented by the COMSOL^{®} package. We checked that results do not depend significantly on boundary conditions. Due to the monomode structure chosen for the planar waveguides, n_{lpw} and n_{upw} take single values. n_{irw} can take several values depending on the lateral confinement i.e. on L_{r}. (Fig. 2
).Results for both systems and both polarizations are similar. Indeed empirically scaling both coordinates merges all curves (Fig. 2(c)). The monomode character is then identical in all cases and is essentially governed by the evanescence length (horizontal scale) and confinement factor (vertical scale) of the fundamental mode of the isolated waveguides (see Eqs. (A2) and A3). Figure 2(c) thus validates the one-mode model for the waveguides of Fig. 1. As expected, they remain monomode as long as n_{irw} remains below the middle of [n_{lpw}, n_{upw}].

We will consider in the following typical cases, i.e. L_{r}=1.5, 3, and 4µm for InP and L_{r}=4µm for GaAs. All the results for these structures, as well as results for Si/SiO_{2} buried-strip waveguides will be merged in Section 4 within a conclusive scaled plot which will emphasize the broad applicability of extended CMT.

#### 3.2 The isolated ridge waveguide

Results for all shallow-ridge structures are similar and will be exemplified by the case of the InP-based structure with L_{r}=3µm for the horizontal polarization. Index reference values are n_{lpw}=3.22974, n_{irw}=3.231632, n_{upw}=3.234355. Figure 3
shows the horizontal section of the fundamental isolated-ridge mode. For all cases calculated the mode tails outside the ridge can be quite well fitted by an exponential decay curve with a decay length L_{decay} very near to the estimation of the lateral evanescence length ${L}_{evan}=\sqrt{{n}_{iw}^{2}-{n}_{lp}^{2}}/{\beta}_{vacuum}$: L_{decay}=2.20µm while L_{evan}=2.23µm for the horizontal polarization in Fig. 3. Based on the map of the isolated ridge mode all parameters of the model (C, η, ξ, and ζ) can be determined by numerical integration (Fig. 4
).As expected by considering only the mode tails, their variation for S above 3µm follows exponential laws with a characteristic length of L_{decay} or 2L_{decay}, depending on the order of the parameter. More precisely η can be quite well fitted by the law $\eta \left(S\right)=\left(1+S/{L}_{decay}\right)\text{\hspace{0.17em}}{e}^{-S/{L}_{decay}}$. Deviations from these laws at high S can be attributed to numerical errors and at low S by the erroneous description of the mode map. However they apply in the useful range of S.

#### 3.3 Supermodes of finite arrays

We now introduce the waveguide array by considering a set of N identical ridges with a center-to-center separation S. As stated above this system now exhibits N sublevels around the center level n_{irw} corresponding to supermodes${a}_{m}\text{\hspace{0.17em}}{e}^{i{k}_{z}z}$ [14]. Their maps are quite well reproduced by superposition of isolated-ridge modes (see Fig. 5
) with weights given by the CMT ${a}_{m}=\mathrm{sin}\left(m\frac{p+1}{N+1}\pi +\left(p+1\right)\frac{\pi}{2}\right)$ where p is the order of the supermode (p=0 for the fundamental mode). No clear evidence of deviation from mere mode superposition has been observed even for strongly-coupled arrays. Polarization hybridization remains small: for instance the TE mode in dense InP-based arrays with L_{r}=3µm, S=4µm has only |E_{y}|/|E_{x}|<0.05.

The set of sublevels n(p) indexed by p are represented in Fig. 6(a) and corresponds to a discrete sample of the diffraction relation (Eq. (5). Data show that an excellent sampling is indeed obtained whatever the number of waveguides (see an example in Fig. 6(b)). We have mostly used arrays with N=7 as a compromise: at higher numbers more modes and hence more data are obtained at the expense of simplicity of numerical calculation, while 7 points correctly describe the diffraction relation. These results are a strong indication that the perturbative treatment CMT relies upon can be used to derive adequate diffraction relations.

## 4. Comparison of analytical modeling to numerical simulation and discussion

#### 4.1 Test of the extended CMT model

Within CMT n(p) is expected to follow a cosine law:

_{r}=NL

_{r}encompassing all merged ridges with ${n}_{r}{}^{2}\left(p\right)\approx {\beta}_{upw}{}^{2}-{\left(p\pi /L{\text{'}}_{r}\right)}^{2}$.

The fit of effective index data by the extended CMT model is two-fold: C controls the height of the diffraction relation i.e. the width of the effective index band while η, ξ, and ζ control its shape. The FEM-calculated width is very well fitted by the model over 3 orders of magnitude in all cases considered, provided that all the values of C obtained from the model are increased by a common factor <1.5 (Fig. 4). This small discrepancy which is actually smaller in other cases (e.g. 1.2 for L_{r}=4µm) can be traced back to the limitations of a perturbative treatment which makes use of isolated waveguide modes (see Appendix). The shape of the band and especially its change at smaller S is fairly well reproduced by the model (Fig. 6(c), dashed lines) and a nearly perfect fit can be obtained (Fig. 6(c), solid lines) by slightly decreasing η. These small discrepancies can be visualized in Fig. 4 through the small shift between directly calculated and FEM-derived values. A similar fit is obtained for all InP or GaAs shallow-ridge waveguides, using smaller corrections of C and η.

Si-SiO_{2} rectangular waveguides at λ=1.55 µm are also considered here for comparison with the work of Cooper and Mookherjea [22], using 200 nm × 500 nm Si core strips (n=3.48) embedded in SiO_{2} cladding medium (n=1.45). In such a system the index contrast is much larger than in shallow ridge waveguides and the coupling can get much stronger, thus making the validity of CMT and even extended CMT questionable. Figure 7(a)
shows the fan-out diagrams of our results – quite similar to those of Ref [22]. – and Fig. 7(b) a comparison with the prediction of extended CMT. In this case the fit is quite good and does not require any adjustment. Extended CMT then applies also well to those strongly-coupled systems.

Finally the excellent overall fit on many various systems validates our extended CMT model and in addition the pertinence of the parameters which describe the strength of the coupling.

#### 4.2 The limit of “weak” coupling – towards a general criterion

We first discuss the case of shallow ridge waveguides. We note that C/β is always <10^{−3}, which confirms the validity of the slowly-varying envelope assumption. On the other hand the correction imposed by the non-orthogonality of individual ridge modes (mostly η, which overrides ξ and ζ) becomes negligible (<1%) only at large values of S i.e. at very low and indeed useless coupling constants (<0.01mm^{−1}). Hence this correction should be considered in most practical cases. In this framework the limit of “weak” coupling below which CMT can be used is defined as the deviation to the CMT diffraction relation that can be tolerated. This maximal tolerated deviation corresponds to a maximum value of η, ξ, and ζ, and in turn for a given system to a maximum coupling strength C determined by curves as those shown in Fig. 4, which are obtained using only the fundamental mode map of the isolated ridge waveguide.

A still simpler though less accurate way may be suggested. Since shallow-ridge structures are all similar, reducing curves to dimensionless coordinates should merge them into a single one. The relevant scale for S is L_{evan}. With the reduced abscissa s=S/L_{evan}, η(s) and ζ(s) curves – and to a lesser extent the noisier ξ(s) curve – for various L_{r} values and for InP as well as GaAs systems do coincide (Fig. 8
). For the C scale, we note that within the CMT the sublevel band of width Δn=4C/β_{vacuum} must remain within the 2D guidance band [n_{lpw},n_{upw}] (see Fig. 2a,b). We then propose to use as a reference the coupling constant C_{max} defined by C_{max}=(β_{vacuum}/4).min(n_{upw}-n_{irw},n_{irw}-n_{lpw}). Figure 8 shows that C/C_{max}(s) curves for all systems also coincide very well. These data can be correctly approximated by the formulas derived from our extended CMT model for narrow ridges (Eqs. (A10-)A13): C/C_{max} ~8e^{-s}, η ~(1+s)e^{-s}, and ξ ~ζ ~2e^{-s}.

Within this “general” rule, prediction becomes easy. For a given structure one calculates the indices of the isolated ridge and planar waveguides and deduces L_{evan} and C_{max}. For a desired maximal correction, the “universal” η(s) yields the minimal s=S/L_{evan} and hence the minimal S. Finally this minimal s yields the maximal C/C_{max} and hence the maximal C. For given planar structure (layer stack) and deviation from CMT, if C has to be maximized, the best choice for L_{r} is to maximize C_{max} by bringing n_{irw} near the middle of [n_{upw}, n_{lpw}]. For the planar structures considered here, C_{max} is limited to ~2.3 and 3.9 mm^{−1} for InP in respectively the horizontal and vertical polarizations, and 0.6 mm^{−1} for GaAs in the horizontal polarization.

Data obtained for Si-SiO_{2} waveguides in vertical (TE type) polarization which involve much shorter evanescence lengths (<0.2 µm) and much stronger confinement factors fit also quite well in this scheme (Fig. 8). Such an overall agreement is somewhat surprising considering the polarization hybridization reported [22]. This effect – which is nevertheless unexpected and has to be taken with caution due to the strong singularities occurring at the corners of the strips – does not seem to prevent satisfactory predictions of interactions between neighboring strips, at least as long as S is not too small. Thus Fig. 8 emphasizes the “universal” character of the above rules which may still be used with rather strongly coupled systems where C reaches 2 µm^{−1} and C/β about 0.2.

#### 4.3 Distortion of the diffraction relation at strong couplings – impact on experiments

On this basis we now describe the distortions of the diffraction relation with respect to its CMT cosine shape (Fig. 9a
) with decreasing s at the moderately small values of s where the model has been validated – at lower values the extended model is no longer valid as attested by the apparent divergence of its diffraction relation. The most distorted part of the curve lies at high k_{x} where a very strong change of negative k_{z} takes place. At low k_{x} the change appears only as a moderate drop of the whole curve; indeed if curves are normalized to take the CMT value of 2 at k_{x}=0 distortion is barely visible (Fig. 9b).

Hence in most experiments involving only low k_{x} the main effect will be a reduction of the coupling constant C to an effective one C_{eff}. Many literature studies have been performed in the GaAs system [14]. For instance a value of C=0.75mm^{−1} is reported for the planar structure we used above with L_{r}=4µm and S=8µm, and is applied to diffraction management studies [3]. Our assessment of this structure on the lines described above leads to C=0.58mm^{−1} but s=1.4 which gives C_{eff}=0.37mm^{−1}. In our own study of “channels” performed on InP-based arrays involving the planar structure we used above with L_{r}=1.5µm and S=5µm, the published value of 1.5mm^{−1} was obtained by fits of n(p) but without taking distortion into account. The correct value is C=1.32mm^{−1} but s=1.24 so C_{eff}=0.81mm^{−1}. Those examples show the importance of corrections in such moderately coupled arrays. However in both cases since only low k_{x} were mainly investigated qualitative agreement was obtained.

If on the other hand large k_{x} are used, the asymmetry of the diffraction curve can lead to important deviations from the CMT behavior. This happens for instance in refraction experiments involving oblique interfaces or following direct injection in high-k_{x} states [25,26]. The distortion can lead to considerable errors on the direction of FB wave packets (beams) which is controlled by the derivative of the diffraction relation and on their divergence which is controlled by its second derivative, if one considers only the low-k_{x} value C_{eff}. In Ref [3]. the authors observe a strong deviation from CMT at k_{x} above π/2 that they attribute to the contribution of the second band but that could be explained as well by mere high-coupling distortion. Finally in arrays which involve a patterning of the coupling constant the contrast between neighboring high-C and low-C zones is reduced by the distortion. In the case of channels [15] – high-C stripes surrounded by large low-C zones – the C contrast 2.83 reduces to a C_{eff} contrast of 2.58, still comfortably large enough to insure superguiding.

## 5. Conclusion and perspectives

In conclusion the present perturbative treatment of arrayed isolated waveguides when they are increasingly coupled to their neighbors yields a first-order extension of CMT. Such a treatment based on a propagative approach gives a simple analytical expression of the propagation constant, which naturally gets rid of the high-order polynomial dependencies involved in the modal approaches. The associated diffraction relation exhibits the distinct consequences of coupling to next-nearest neighbors on lateral diffraction and of non-orthogonality between nearest neighbors on propagation, respectively through two transfer integrals and one overlap integral. The distortion of the shape of the diffraction relation with respect to the CMT shape at increasingly stronger coupling situations is thus nicely accounted for. Quantitative agreement with fully numerical simulation is obtained even nearby the strong coupling regime for various waveguide structures, requiring only slight adjustments of the ab initio parameters. In a broad range of weakly and moderately coupled systems the distortion of the diffraction relation can be assessed using the non-orthogonality overlap integral as an indicator, so that the designer can confidently avoid excessive distortions. Alternatively he can take advantage of fairly simple empirical formulas to improve propagation and diffraction management in the design of arrayed devices.

## Appendix – Analytical derivations

In the above discussion numerical resolution of the Maxwell equations for shallow ridge waveguide arrays is compared to more straightforward results derived from approximate models. As shown on the left hand side of Fig. A1 the array is considered as the periodical juxtaposition of planar waveguide portions.

Since the dielectric constant is piecewise constant, the Maxwell equations can be reduced to the Helmholtz equation:

The solutions are superpositions of planar modes${g}_{m}^{R,E}\left(Y\right)\text{\hspace{0.17em}}$, m=1 … (taken as real) in the “ridge” or “etch” (R,E) regions, which have to be connected at the region boundaries using continuity equations. In the following we assume that only the fundamental planar modes with propagation constants ${\beta}_{1}^{R,E}$ need to be considered, either in TE or in TM polarization (${\beta}_{1}^{R}$ and ${\beta}_{1}^{E}$ are denoted β_{upw} and β_{lpw} in the text). This is quite legitimate for shallow ridge waveguides in which the core layer is separated from the etched mesas by a thick cladding layer, and provided the mesas exhibit good verticality in order to minimize polarization mode mixing. Then the confined mode of the isolated waveguide and the modes of the array – the super-modes – are built using two leftward and rightward propagating waves, with vertical profiles given by the two fundamental planar modes in either region. The horizontal profile of the fundamental (symmetric) mode of an isolated waveguide can thus be cast into the form derived for three-layer planar waveguides [27]:

where *w* is the waveguide width, and q=1/L_{evan} obeys:

*r _{R/E}* being an averaged value of the ratio of the dielectric profiles, either equal or very close to 1 depending on the polarization.

On the right-hand side of Fig. A1, the array is decomposed into the usual perturbation scheme of an isolated waveguide plus a dielectric perturbation $\sum _{n}\delta {\epsilon}_{n}\left(X,Y\right)$ brought about by the mesas of all the other waveguides. In this approach the Helmholtz equation is solved using a linear superposition of the fundamental modes M_{1} of all the isolated waveguides:

The Helmholtz equation is not strictly equivalent to the Maxwell equations because of the abrupt changes of the dielectric constant around the perturbation by the mesas, and the approach is valid only for small perturbations. Another way to envision its limitations is to recognize that expansion of a supermode over the isolated waveguide modes Eq. A4 cannot satisfy all the boundary conditions at the mesa edges because the evanescent tails of modes m’≠m do not at mesa m; this would require a deeper discussion of the polarization properties of the modes, which is beyond the scope of the present analysis. Using the definitions of C and of the overlaps by Eqs. (1) and (2) and approximating up to the second neighbor, the perturbation approach [27] yields:

Since the mode has exponential tails, $\u300801\u3009/\u300800\u3009$,$\u3008101\u3009/\u3008001\u3009$,$\u3008\overline{1}01\u3009/\u3008001\u3009$, and$\u3008002\u3009/\u3008001\u3009$ are first order while$\u300802\u3009/\u300800\u3009$ is second order with respect to evanescent coupling. At 0th order one finds the CMT result. At 1st order, denoting z=CZ and using definitions of η, ξ, and ζ by Eqs. (4), one finds:

The diffraction relation of FB waves ${\text{e}}^{i{k}_{x}x+i{k}_{z}z}$ then becomes:

To first order in η, the propagation equation Eq. A6 is identical to what can be derived from the approach based on Wannier functions [23] when the diffraction relation obeys Eq. A7.

Now, using Eq. A3 and the profile of Eq. A2, we calculate the overlap integrals and the extended CMT parameters C, η, ξ, and ζ, and for instance:

At low ridge widths, formulas reduce to simple expressions:

where ${\epsilon}_{1}^{R}={\beta}_{1}^{R}{}^{2}/{\beta}_{vacuum}{}^{2}$, ${\epsilon}_{1}^{E}={\beta}_{1}^{E}{}^{2}/{\beta}_{vacuum}{}^{2}$ (in the main text n_{upw}=$\sqrt{{\epsilon}_{1}^{R}}$ and n_{lpw}=$\sqrt{{\epsilon}_{1}^{E}}$). Taking the planar waveguide in the etched zone as an unperturbed system, a perturbative treatment of the Helmholtz equation for the planar waveguide in the ridge zone straightforwardly gives to first perturbation order:

so that assuming r_{R/E}=1:

The above treatment can also be applied to the case of buried strips [22] though in this case a modal method through the periodical slit and strip arrangement might be more effective. Equivalent roles can be played between the above “ridge” and “etched” regions and, respectively, the “strip” and “aperture” regions. A difficulty arises in the aperture regions which are uniform and do not support localized modes. This difficulty can be circumvented by symmetrically adding new strips to the system in the aperture regions, far from the original strip array in the ±Y directions. Such strips give rise to extra confined modes, the highest of which can be both degenerate and weakly coupled with the planar mode for the strip regions, and hence serve as a planar mode for the aperture regions. As a perturbative approach shows, the new system will exhibit similar strip array modes as the original one near the planar mode for the strip region, due to weak evanescent coupling of the latter with the extra modes. Eqs. A10 and A13 are then also valid estimates for the fundamental band in buried strip arrays as long as higher order planar modes can be neglected.

## Acknowledgements

These results are within the scope of C’nano IdF; C’nano IdF is a CNRS, CEA, MESR and Région Ile-de-France Nanosciences Competence Center.

## References and links

**1. **D. N. Christodoulides, F. Lederer, and Y. Silberberg, “Discretizing light behaviour in linear and nonlinear waveguide lattices,” Nature **424**(6950), 817–823 (2003). [CrossRef] [PubMed]

**2. **J. Fleischer, G. Bartal, O. Cohen, T. Schwartz, O. Manela, B. Freedman, M. Segev, H. Buljan, and N. Efremidis, “Spatial photonics in nonlinear waveguide arrays,” Opt. Express **13**(6), 1780–1796 (2005), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-13-06-1780. [CrossRef] [PubMed]

**3. **H. S. Eisenberg, Y. Silberberg, R. Morandotti, and J. S. Aitchison, “Diffraction management,” Phys. Rev. Lett. **85**(9), 1863–1866 (2000). [CrossRef] [PubMed]

**4. **T. Pertsch, P. Dannberg, W. Elflein, A. Brauer, and F. Lederer, “Optical Bloch oscillations in temperature tuned waveguide arrays,” Phys. Rev. Lett. **83**(23), 4752–4755 (1999). [CrossRef]

**5. **R. Morandotti, U. Peschel, J. S. Aitchison, H. S. Eisenberg, and Y. Silberberg, “Experimental observation of linear and nonlinear optics Bloch oscillations,” Phys. Rev. Lett. **83**(23), 4756–4759 (1999). [CrossRef]

**6. **F. Dreisow, A. Szameit, M. Heinrich, T. Pertsch, S. Nolte, A. Tünnermann, and S. Longhi, “Decay control via discrete-to-continuum coupling modulation in an optical waveguide system,” Phys. Rev. Lett. **101**(14), 143602 (2008). [CrossRef] [PubMed]

**7. **T. Schwartz, G. Bartal, S. Fishman, and M. Segev, “Transport and Anderson localization in disordered two-dimensional photonic lattices,” Nature **446**(7131), 52–55 (2007). [CrossRef] [PubMed]

**8. **Y. Lahini, A. Avidan, F. Pozzi, M. Sorel, R. Morandotti, D. N. Christodoulides, and Y. Silberberg, “Anderson localization and nonlinearity in one-dimensional disordered photonic lattices,” Phys. Rev. Lett. **100**(1), 013906 (2008). [CrossRef] [PubMed]

**9. **K. G. Makris, D. N. Christodoulides, O. Peleg, M. Segev, and D. Kip, “Optical transitions and Rabi oscillations in waveguide arrays,” Opt. Express **16**(14), 10309–10314 (2008), http://www.opticsexpress.org/abstract.cfm?URI=oe-16-14-10309. [CrossRef] [PubMed]

**10. **D. N. Christodoulides and E. D. Eugenieva, “Blocking and routing discrete solitons in two-dimensional networks of nonlinear waveguide arrays,” Phys. Rev. Lett. **87**(23), 233901 (2001). [CrossRef] [PubMed]

**11. **A. Fratalocchi, G. Assanto, K. A. Brzdakiewicz, and M. A. Karpierz, “All-optical switching and beam steering in tunable waveguide arrays,” Appl. Phys. Lett. **86**(5), 051112 (2005). [CrossRef]

**12. **R. A. Vicencio, M. I. Molina, and Y. S. Kivshar, “Switching of discrete optical solitons in engineered waveguide arrays,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **70**(2), 026602 (2004). [CrossRef] [PubMed]

**13. **A. L. Jones, “Coupling of optical fibers and scattering in fibers,” J. Opt. Soc. Am. **55**(3), 261–269 (1965). [CrossRef]

**14. **F. Lederer, G. I. Stegeman, D. N. Christodoulides, G. Assanto, M. Segev, and Y. Silberberg, “Discrete solitons in optics,” Phys. Rep. **463**(1-3), 1–126 (2008) (and references therein). [CrossRef]

**15. **N. Belabas, S. Bouchoule, I. Sagnes, J. A. Levenson, C. Minot, and J. M. Moison, “Confining light flow in weakly coupled waveguide arrays by structuring the coupling constant: towards discrete diffractive optics,” Opt. Express **17**(5), 3148–3156 (2009), http://www.opticsexpress.org/abstract.cfm?URI=oe-17-5-3148. [CrossRef] [PubMed]

**16. **J. M. Moison, N. Belabas, C. Minot, and J. A. Levenson, “Discrete photonics in waveguide arrays,” Opt. Lett. **34**(16), 2462–2464 (2009), http://www.opticsinfobase.org/ol/abstract.cfm?URI=ol-34-16-2462. [CrossRef] [PubMed]

**17. **A. Szameit, H. Trompeter, M. Heinrich, F. Dreisow, U. Peschel, T. Pertsch, S. Nolte, F. Lederer, and A. Tünnermann, “Fresnel’s laws in discrete optical media,” N. J. Phys. **10**(10), 103020 (2008). [CrossRef]

**18. **A. Hardy and W. Streifer, “Coupled mode theory of parallel waveguides,” J. Lightwave Technol. **3**(5), 1135–1146 (1985). [CrossRef]

**19. **W. P. Huang, “Coupled-mode theory for optical waveguides: an overview,” J. Opt. Soc. Am. A **11**(3), 963–983 (1994). [CrossRef]

**20. **L. Eyges and P. Wintersteiner, “Modes of an array of dielectric waveguides,” J. Opt. Soc. Am. **71**, 1351–1360 (1981), http://www.opticsinfobase.org/abstract.cfm?URI=josa-71-11-1351.

**21. **A. Kaplan and S. Ruschin, “Characterization and performance evaluation of coupled multiwaveguide arrays,” J. Lightwave Technol. **17**(10), 1884–1889 (1999), http://jlt.osa.org/abstract.cfm?URI=JLT-17-10-1884. [CrossRef]

**22. **M. L. Cooper and S. Mookherjea, “Numerically-assisted coupled-mode theory for silicon waveguide couplers and arrayed waveguides,” Opt. Express **17**(3), 1583–1599 (2009), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-17-03-1583. [CrossRef] [PubMed]

**23. **G. L. Alfimov, P. G. Kevrekidis, V. V. Konotop, and M. Salerno, “Wannier functions analysis of the nonlinear Schrödinger equation with a periodic potential,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **66**(4), 1046608 (2002). [CrossRef]

**24. **H. S. Eisenberg, Y. Silberberg, R. Morandotti, A. Boyd, and J. S. Aitchison, “Discrete spatial optical solitons in waveguide arrays,” Phys. Rev. Lett. **81**(16), 3383–3386 (1998). [CrossRef]

**25. **T. Pertsch, T. Zentgraf, U. Perchel, A. Brauer, and F. Lederer, “Anomalous refraction diffraction in discrete optical systems,” Phys. Rev. Lett. **88**(9), 093901 (2002). [CrossRef] [PubMed]

**26. **D. Mandelik, H. S. Eisenberg, Y. Silberberg, R. Morandotti, and J. S. Aitchison, “Band-gap structure of waveguide arrays and excitation of Floquet-Bloch solitons,” Phys. Rev. Lett. **83**, 4752–4754 (1999).

**27. **B. E. A. Saleh, and M. C. Teich, “Fundamentals of Photonics,” Wiley-Interscience, 2007.