Geometric variations in high index-contrast waveguides , coupled mode theory in curvilinear coordinates

Perturbation theory formulation of Maxwell’s equations gives a theoretically elegant and computationally efficient way of describing small imperfections and weak interactions in electro-magnetic systems. It is generally appreciated that due to the discontinuous field boundary conditions in the systems employing high dielectric contrast profiles standard perturbation formulations fail when applied to the problem of shifted material boundaries. In this paper we developed a novel coupled mode and perturbation theory formulations for treating generic non-uniform (varying along the direction of propagation) perturbations of a waveguide cross-section based on Hamiltonian formulation of Maxwell equations in curvilinear coordinates. We show that our formulation is accurate and rapidly converges to an exact result when used in a coupled mode theory framework even for the high index-contrast discontinuous dielectric profiles. Among others, our formulation allows for an efficient numerical evaluation of induced PMD due to a generic distortion of a waveguide profile, analysis of mode filters, mode converters and other optical elements such as strong Bragg gratings, tapers, bends etc., and arbitrary combinations of thereof. To our knowledge, this is the first time perturbation and coupled mode theories are developed to deal with arbitrary non-uniform profile variations in high index-contrast waveguides. c © 2002 Optical Society of America OCIS codes: (060.2310) Fiber optics; (060.2400) Fiber properties ; (060.2280) Fiber design and fabrication References and links 1. Steven G. Johnson, Mihai Ibanescu, M. Skorobogatiy, Ori Weisberg, Torkel D. Engeness, Marin Soljačić, Steven A. Jacobs, J. D. Joannopoulos, and Yoel Fink, “Low-loss asymptotically single-mode propagation in large-core OmniGuide fibers,” Opt. Express 9, 748 (2001), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-9-13-748. 2. Steven G. Johnson, Mihai Ibanescu, M. Skorobogatiy, Ori Weisberg, J. D. Joannopoulos, and Yoel Fink, “Perturbation theory for Maxwell’s equations with shifting material boundaries,” Phys. Rev. E 65, 66611 (2002). 3. M. Skorobogatiy, Mihai Ibanescu, Steven G. Johnson, Ori Weisberg, Torkel D. Engeness, Marin Soljačić, Steven A. Jacobs, and Yoel Fink, “Analysis of general geometric scaling perturbations in a transmitting waveguide: fundamental connection between polarization-mode dispersion and group-velocity dispersion,” Opt. Soc. Am. B 19, (2002). 4. M. Skorobogatiy, Steven A. Jacobs, Steven G. Johnson, and Yoel Fink, “Dielectric profile variations in high index-contrast waveguides, coupled mode theory and perturbation expansions,” to be published in J. Opt. Soc. Am. B, 2003. (C) 2002 OSA 21 October 2002 / Vol. 10, No. 21 / OPTICS EXPRESS 1227 #1640 $15.00 US Received August 30, 2002; Revised October 18, 2002 5. M. Lohmeyer, N. Bahlmann, and P. Hertel, “Geometry tolerance estimation for rectangular dielectric waveguide devices by means of perturbation theory,” Opt. Communications 163, pp. 86–94 (1999). 6. N. R. Hill,“Integral-equation perturbative approach to optical scattering from rough surfaces,” Phys. Rev. B 24, p. 7112 (1981). 7. D. Marcuse, Theory of dielectric optical waveguides (Academic Press, 2nd ed., 1991). 8. A. W. Snyder and J. D. Love, Optical waveguide theory (Chapman and Hall, London, 1983). 9. B. Z. Katsenelenbaum, L. Mercader del Ŕıo, M. Pereyaslavets, M. Sorolla Ayza, and M. Thumm, Theory of Nonuniform Waveguides: The Cross-Section Method (Inst. of Electrical Engineers, London, 1998). 10. L. Lewin, D. C. Chang, and E. F. Kuester, Electromagnetic waves and curved structures (IEE Press, Peter Peregrinus Ltd., Stevenage 1977). 11. F. Sporleder and H. G. Unger, Waveguide tapers transitions and couplers (IEE Press, Peter Peregrinus Ltd., Stevenage 1979). 12. H. Hung-Chia, Coupled mode theory as applied to microwave and optical transmission (VNU Science Press, Utrecht 1984). 13. Steven G. Johnson, P. Bienstman, M. A. Skorobogatiy, M. Ibanescu, E. Lidorikis, and J. D. Joannopoulos “The adiabatic theorem and a continuous coupled-mode theory for efficient taper transitions in photonic crystals,” to be published in Phys. Rev. E, 2002. 14. L. D. Landau and E. M. Lifshitz, Quantum mechanics (non-relativistic theory) (Butterworth Heinemann, 2000). 15. C. Vassallo, Optical waveguide concepts (Elsevier, Amsterdam, 1991). 16. R. Holland, “Finite-difference solution of Maxell’s equation in generalized nonorhogonal coordinates,” IEEE Trans. Nucl. Sci. 30, 4589 (1983). 17. J.P. Plumey, G. Granet, and J. Chandezon, “Differential covariant formalism for solving Maxwell’s equations in curvilinear coordinates: oblique scattering from lossy periodic surfaces,” IEEE Trans. on Antennas Propag. 43, 835 (1995). 18. E.J. Post, Formal Structure of Electromagnetics (Amsterdam: North-Holland, 1962). 19. F.L. Teixeira and W.C. Chew, “Analytical derivation of a conformal perfectly matched absorber for electromagnetic waves,” Microwave Opt. Technol. Lett. 17, 231 (1998). 20. P. Bienstman, software at http://camfr.sf.net.


Introduction
Standard perturbation and coupled mode theory formulations are known to fail or exhibit a very slow convergence [1][2][3][4][5][6] when applied to the analysis of geometrical variations in the structure of high index-contrast fibers.
In a uniform coupled mode theory framework (waveguide profile remains unchanged along the direction of propagation), eigenvalues of the matrix of coupling elements approximate the values of the propagation constants of a uniform waveguide of perturbed cross-section.When large enough number of modes are included coupled mode theory, in principle, should converge to an exact solution for perturbations of any strength.Perturbation theory is numerically more efficient method than coupled mode theory, but it is mostly applicable to the analysis of small perturbations.For stronger perturbations, higher order perturbation corrections must be included converging, in the limit of higher orders, to an exact solution.In a non-uniform (waveguide profile is changing along the direction of propagation) coupled mode and perturbation theories one propagates the modal coefficients along the length of a waveguide using a first order differential equation involving a matrix of coupling elements.Both uniform and non-uniform coupled mode and perturbation theory expansions rely on the knowledge of correct coupling elements.
Conventional approach to the evaluation of the coupling elements proceeds by expansion of the solution for the fields in a perturbed waveguide into the modes of an unperturbed system, then computes a correction to the Hamiltonian of a problem due to the perturbation in question and, finally, computes the required coupling elements.Unfortunately, this approach encounters difficulties when applied to the problem of find- ing perturbed electromaghetic modes in the waveguides with shifted high index-contrast dielectric boundaries.In particular, it was demonstrated (see [4], for example) that for a uniform geometric perturbation of a fiber profile with abrupt high index-contrast dielectric interfaces, expansion of the perturbed modes into an increasing number of the modes of an unperturbed system does not converge to a correct solution when standard form of the coupling elements [7,8] is used.Mathematical reasons of such a failure are still not completely understood but probably lie either in the incompletness of the basis of eigenmodes of an unpertubed waveguide in the domain of the eigenmodes of a perturbed waveguide or in the fact that the standard mode orthogonality condition (6) does not constitute a strict norm.We would like to point out that standard coupled mode theory can still be used even in the problem of finding the modes of a high index-contrast waveguide with sharp dielectric interfaces.One can calculate such modes by using as an expansion basis eigen modes of some "smooth" dielectic profile (empty metallic waveguide, for example).However, the convergence of such a method with respect to the number of basis modes is very slow (at most linear).Perturbation formulation withing this approach is also problematic, and even for small geometric variations of waveguide profile matrix of coupling element has to be recomputed anew.
Other methods developed to deal with shifting metallic boundaries and dielectric interfaces originate primarily from the works on metallic waveguides and microwave circuits [7,[9][10][11][12][13].There, however, Hermitian nature of the Maxwell's equations in the problem of radiation propagation along the waveguides is not emphasized, and consequent development of perturbation expansions is usually omitted.Moreover, dealing with non-uniform waveguides, these formulations usually employ an expansion basis of "instantaneous modes".Such modes have to be recalculated at each different waveguide cross-section leading to very computationally demanding propagation schemes.
In this paper we introduce a method of evaluating the coupling elements which is valid for any analytical geometrical waveguide profile variations and high index contrast using the eigen modes of an unperturbed waveguide as an expansion basis.To derive a correct form of the coupling elements we, first, define a convenient way of specifying geometrical perturbations using a technique of coordinate mapping, and then construct a novel expansion basis using spatially stretched modes of an unperturbed waveguide.Stretching is performed in such a way as to match the regions of the field discontinuities in the expansion modes with the positions of the perturbed dielectric interfaces.Thus defined expansion modes are designed to have continuous fields in the regions of the continuous dielectric of a perturbed index profile.By substituting such expansions into the Maxwell's equations, we later find the required expansion coefficients.It becomes more convenient to perform further algebraic manipulations in a coordinate system where stretched expansion modes become again unperturbed modes of an original waveguide.Thus, the final steps of evaluation of the coupling elements involve transforming and manipulation of Maxwell's equations in the perturbation matched curvilinear coordinates.
In further discussions we formulate geometrical waveguide profile variations in terms of an analytical mapping of an unperturbed dielectric profile onto a perturbed one.Given a perturbed dielectric profile (x , y , z ) in a Euclidian system of coordinates (x , y , z ) (where z is a direction of propagation) we define a mapping (x (x, y, z), y (x, y, z), z (x, y, z)) such that (x (x, y, z), y (x, y, z), z (x, y, z)) corresponds to an unperturbed profile in a curvilinear coordinate system associated with (x, y, z).We then perform a coordinate transformation from a Euclidian system of coordinates (x , y , z ) into a corresponding curvilinear coordinate system (x, y, z) by rewriting Maxwell's equations in such a coordinate system.Finally, because dielectric profile in a coordinate system (x, y, z) is that of an unperturbed waveguide, we can use the basis set of eigen modes of an unperturbed Fiber profile remains cylindrically symmetric, while the radii of the dielectric interfaces along the direction of propagation s become ρ i (1 + δ s L ). c) Scaling variation -sinusoidal Bragg gratings.Fiber profile remains cylindrically symmetric, while the radii of the dielectric interfaces along the direction of propagation s become ρ i (1+δSin(2π s Λ )).d) Non-concentric variation.Each dielectric interface stays cylindrically symmetric, while the center line is bent.
Hamiltonian in (x, y, z) coordinates to calculate coupling matrix elements due to the geometrical perturbation of a waveguide profile.
Our paper is organized as following.We first describe some typical geometrical variations of fiber profiles.Next, we discuss properties of general curvilinear coordinate transformations and relate them to a particular case of coordinate mappings describing geometrical variations of a fiber profile.We then formulate Maxwell's equations in a general curvilinear coordinate system.We apply this formulation to develop the coupled mode and perturbation theories using stretched eigen states of an unperturbed Hamiltonian as an expansion basis.We conclude with a detailed study of the convergence of the novel coupled mode and perturbation theories for a variety of uniform and non-uniform variations of the geometric profiles in high index-contrast waveguides.
In the following, we focus on the geometrical variation in fibers, although geometric variations in generic waveguides can be readily described by the same theory (see [3] for planar waveguides, for example).

Types of geometrical variations of fiber profiles
We start by considering some common geometrical variations of a fiber profile that can be either deliberately designed or arise during manufacturing as perturbations.
Let (x, y, z) correspond to Euclidian coordinate system.In cylindrical coordinates Then, for step index multilayer fibers, for example, the position of the i th dielectric interface of radius ρ i can be described by a set of points (1), ρ = ρ i ,θ ∈ (0, 2π) (see figure 1a)).Geometrical variations in multilayer fibers that will be considered in this paper for the purpose of demonstration of our method include high index contrast Bragg gratings, tapers and bends.
If all the radii of the dielectric interfaces are scaled from their original values ρ i to become ρ i (1 + f (s)) where f (s) is an analytic function of propagation distance s, then the following coordinate mapping will describe such a scaling variation of interface radii along the direction of propagation Here, ρ = ρ i , θ ∈ (0, 2π) define the points on the i th unperturbed dielectric interface, while the corresponding (x(ρ, θ, s), y(ρ, θ, s), z(ρ, θ, s)) describe the points on the perturbed interfaces.
Using description of the points on the perturbed dielectric interfaces in terms of the coordinate transformation (2) we introduce the following optical elements: 1) Dielectric tapers described by the monotonic variations in the radii of the dielectric layers in a cylindrically symmetric fiber (see figure 1b)).For the case of linear tapers f (s) = δ s L where L is a length of a taper and δ is a relative change in the taper dimensions from the beginning to the end.
2) Dielectric Bragg gratings of any strength, where the radii of the dielectric layers change according to some analytic periodic function along the direction of propagation in a cylindrically symmetric fiber (see figure 1c)).For the case of sinusoidal gratings, for example, f (s) = δsin(2π s Λ ), where δ defines the magnitude of the radial variation of the grating dimensions, while Λ corresponds to the pitch of the grating.
In the same manner as for the case of concentric perturbations 1) and 2), we define non-concentric perturbations through the respective coordinate transformation which for the case of bends, for example, is defined as following: 3) Assuming that the center of a fiber spans an arc of radius R (see figure 1d)), then the following coordinate mapping will describe the points on the perturbed dielectric interfaces due to a bend where s is a length along a circular arc.Here, ρ = ρ i , θ ∈ (0, 2π) define the points on the i th unperturbed dielectric interface, while the corresponding (x(ρ, θ, s), y(ρ, θ, s), z(ρ, θ, s)) of (3) describe the points on the perturbed bent interfaces.

Coupled mode theory for Maxwell's equations
In the following, we introduce a Hamiltonian formulation of Maxwell's equations in terms of transverse field coordinates to address radiation propagation in generic nonuniform waveguides.Hamiltonian formulation in regular Euclidian coordinates in terms of the transverse fields is known and detailed in [1,3].Hamiltonian formulation in curvilinear perturbation matched coordinates for the case of uniform fibers of arbitrary cross-sections is detailed in [3,4], where perturbation induced deterministic P MD is considered in length.
Assuming the standard time dependence of the electro-magnetic fields F(x, y, z, t) = F(x, y, z) exp −iωt, (F denotes electric or magnetic field vector) and introducing transverse and longitudinal components of the fields F = F t + F z with respect to the propagation direction z, Maxwell's equations can be written in terms of the transverse field components as where operators Â and B are and Dirac notation is introduced |ψ = E t H t .In this form operators on the left and on the right of equation ( 4) are Hermitian.
In the case of a uniform waveguide profile operator Â is going to be a function of the transverse coordinates (x, y) only, and we will refer to such an operator as Â0 .In that case, fields -solutions of ( 4) will possess an extra symmetry F(x, y, z, t) = F(x, y)expiβz − iωt.An orthogonality condition between the modes β and β could then be taken in the form (see [3]) where pure real β's correspond to the guided modes, pure imaginary β's correspond to the evanescent modes and δ β,β stands for the Kronecker δ.Propagation problem (4) will then reduce to a generalized Hermitian eigen value problem In the current paper we assume a closed non-uniform waveguide with a linear isotropic lossless and non-magnetic media so that only the pure guided and evanescent waves are present in the field expansions.Strictly speaking, this formulation avoids the analysis of radiation loss.However, current framework still allows study of radiation propagation in leaky systems (Bragg Fibers and Photonic Crystal (PC) waveguides in example below) in certain regimes.Particular methods and their justifications are beyond the scope of this paper and they will only be mentioned in passing.For Bragg Fibers and PC waveguides of infinite number of confining layers, for example, radiation leakage in the band gap is zero with modes in the band gap being pure guided and integrable.If scattering to other guided modes in the band gap dominates scattering to the continuum of radiation modes (as frequently is the case in the perturbed multimoded air core PC of Bragg Fibers) one can use current formulation without modification.In the case of leaky modes with small radiation loss (originating from the finite number of the confining layers in the PC or Bragg fibers, for example) one can still use the modes of an ideal lossless system to approximate the modal scattering for the propagation lengthes on which the radiation losses of the corresponding leaky modes are negligible.In general, introducing absorbing boundary conditions together with conducting boundary conditions outside the absorber recasts the problem of dealing with non-integrable leaky modes of an open waveguide onto the problem of dealing with lossy (due to absorption and inosotropy of absorbing dielectric layer) but integrable modes of a closed waveguide.Modification of the current method is possible to include absorption loss and inosotropy of corresponding dielectric materials thus, in principle, allowing the treatment of leaky modes in non-uniform waveguides.
In the following, we analyze variations of a waveguide profile that are non-uniform along the direction of propagation.Uniform perturbations of a waveguide profile has been analyzed in our previous works [2][3][4].If a variation such as a continuous change in the sizes of transverse dimensions of a structure is introduced into a system, it will modify the operator Â making it, generally, dependent on all the coordinates (x, y, z).We denote correction to an original uniform operator Â0 by ∆ Â(z).Then, propagation equation ( 4) can be rewritten as where |ψ is a corresponding non-uniform solution.One solves equation ( 8) by expanding a solution |ψ into a linear combination of some basis functions.If variation corresponding to ∆ Â(z) is small, then the natural choice of such basis functions would be uniform solutions |ψ 0 β of an unperturbed operator Â0 , where C i (z) form a vector of expansion coefficients C. Substitution of ( 9) into ( 8) and further use of the orthogonality condition ( 6) lead to the following set of coupled linear differential equations with non-uniform coefficients where normalization matrix B has entries B i,j = ψ 0 β * i | B|ψ 0 βj , and matrix of coupling elements ∆A has entries ∆A i,j = ψ 0 Moreover, if perturbation of a fiber profile is small scattering problem (10) allows a perturbative analysis [14].Thus, assuming that all incoming radiation is initially in the guided mode n, then C n (z) ≈ 1 and scattered power into the guided mode j can be found from population coefficients Standard methods of evaluation of the coupling matrix elements due to the geometrical distortions of the dielectric interfaces [7,8,15] employ eigen modes of an unperturbed waveguide while shifting the dielectric boundaries to extract additional contributions from Maxwell's equations to the original Hamiltonian.While intuitive, these methods are known to fail when dealing with high index-contrast sharp interfaces (see detailed discussion in [1][2][3][4]).
In the following we derive the form of the coupling matrix ∆A of (10) using transformation of the Maxwell's equations into the curvilinear coordinates that match the geometrical variation in question.This way geometrical variations are interpreted as a change in the geometrical metric of space, while the shapes of dielectric interfaces stay intact allowing the use of the eigen modes of an unperturbed waveguide in a corresponding curvilinear coordinate system.

Curvilinear coordinate systems
Following [16,17], we first introduce general properties of the curvilinear coordinate transformations.Let (x 1 , x 2 , x 3 ) be the coordinates in a Euclidian coordinate system.We introduce an analytical mapping into a new coordinate system with coordinates (q 1 , q 2 , q 3 ) as (x 1 (q 1 , q 2 , q 3 ), x 2 (q 1 , q 2 , q 3 ), x 3 (q 1 , q 2 , q 3 )).A new coordinate system can be characterized by its covariant basis vectors a i defined in the original Euclidian system as Note, if all the coordinates except q j are kept fixed, then a j is tangential to the set of points described in Euclidian coordinate system by a curve (x 1 (q i = const, q j , q k = const), x 2 (q i = const, q j , q k = const), x 3 (q i = const, q j , q k = const)).Now, define reciprocal (contravariant) vector ai as where metric g ij is defined as and g = det(g ij ).From the definition of reciprocal vector a i , it is perpendicular to the surface of constant q i described by (x 1 (q i = const, q j , q k ), x 2 (q i = const, q j , q k ), x 3 (q i = const, q j , q k )) where (j = k = i).Vectors a i and their reciprocal a i satisfy the following orthogonality conditions where g ij is an inverse of the metric g ij .In general, a vector may be represented by its covariant components E = e i a i or by its contravariant components E = e i a i .These components might have unusual dimensions because the underlying vectors a i and a i are not properly normalized in a Euclidian coordinate system.Components having the usual dimensions are defined by where i i , i i are unitary vectors.Normalized covariant and contravariant components are connected by We demonstrate the properties of curvilinear coordinate transformations on the example of a linear taper variation of a fiber profile (2) where f (s) = δ s L .To ease further presentation we will use a more common notation for Euclidian coordinates (x, y, z) instead of (x 1 , x 2 , x 3 ) and curvilinear coordinates (ρ, θ, s) instead of (q 1 , q 2 , q 3 ).
Coordinate mapping due to a linear taper variation of a fiber profile (see figure 1b)) expressed in the polar coordinates (ρ, θ, s) is Covariant and contravariant basis vectors corresponding to a curvilinear coordinate transformation (16) are Thus, coordinate transformation corresponding to a linear taper variation is not an orthogonal transformation, and consequently the covariant and contravariant components behave differently.Covariant field component E ρ i ρ is always perpendicular to the tapered dielectric interface, while a covariant field component E s i s is not parallel to the tapered interface, however, when δ → 0, E s i s remains close to being parallel to the elliptical interface up to the first order accuracy in δ.In a similar fashion, a contravariant field component E s i s is always parallel to the elliptical dielectric interface, while a contravariant field component E ρ i ρ is not perpendicular to the elliptical interface, however, when δ → 0, E ρ i ρ remains close to being perpendicular to the elliptical interface up to the first order accuracy in δ.
In the examples above and in the forthcoming presentation we concentrate on the analysis of the geometrical fiber profile perturbations that can be described by a general coordinate transformation of the form x(ρ, θ, s), y(ρ, θ, s), z(ρ, θ, s), (18) where (x, y, z) are the coordinates in a Euclidian coordinate system.In an unperturbed cylindrically symmetric fiber, the i th dielectric interface can be characterized by its radius ρ i and a corresponding set of points (x = ρCos(θ), y = ρSin(θ), z = s), where ρ = ρ i , θ ∈ (0, 2π) for any s.Given the mapping (18), perturbed i th dielectric interface while changing along the direction of propagation s can still be characterized by its unperturbed radius ρ i and for every fixed s is represented by a set of points (x(ρ, θ, s), y(ρ, θ, s), z(ρ, θ, s)), ρ = ρ i , θ ∈ (0, 2π) (frequently, it is convenient to assume z = s).As follows from the discussion of this section vectors forming local covariant and contravariant coordinate systems associated with coordinates (ρ, θ, s) will have the following properties on the dielectric interfaces.Contravariant basis vectors: i ρ is strictly perpendicular to the perturbed dielectric interface, i s is almost parallel to the perturbed dielectric interface more so the less is the perturbation of a fiber profile.Covariant basis vectors: i ρ is almost perpendicular to the perturbed dielectric interface more so the less is the perturbation of a fiber profile, i s is strictly parallel to the perturbed dielectric interface.
To summarize, geometric perturbations of a fiber profile can be characterized by a small parameter (δ in the examples above) that defines a curvilinear coordinate mapping of an unperturbed profile expressed in coordinates (ρ, θ, s) onto a perturbed profile expressed in Euclidian coordinates (x(ρ, θ, s), y(ρ, θ, s), z(ρ, θ, s)).This mapping defines a curvilinear and, generally, a nonorthogonal coordinate system of coordinates (ρ, θ, s).The smaller is the profile perturbation the closer is the corresponding curvilinear coordinate system to an original orthogonal system.

Maxwell's equations in curvilinear coordinates
We now present a formulation of Maxwell's equations in general curvilinear coordinates.A well known form of Maxwell's equations in curvilinear coordinates can be found in a variety of references [16][17][18][19].These expressions are compactly expressed in terms of the normalized covariant and contravariant components of the fields and in the absence of free electric currents they are (q 1 , q 2 , q 3 ) where e ijk is a Levi-Civita symbol.For further derivations it is most convenient to write Maxwell's equations in terms of the normalized covariant components only.This will allow us later in the paper to express the coupling elements only in terms of the integrals over the fields avoiding the integrals over the field derivatives.For a non-uniform fiber, propagation constant is not a conserved property, and the only symmetry possessed by a solution is Assuming non-magnetic materials µ = 1 we express Maxwell's equations (19) in terms of the normalized covariant coordinates only Note that by substituting s components of the electric and magnetic fields ( 22) into (21) we can form the non-uniform Hamiltonian formulation of the Maxwell's equations in the curvilinear coordinates in terms of the differential equation relating the covariant components of transverse fields only.

Coupled mode theory for Maxwell's equations in curvilinear coordinates
In the previous sections we introduced geometrical variations of waveguide cross-section in terms of the analytical mapping of an unperturbed onto a perturbed waveguide profile.Next, we construct an expansion basis using the eigen fields of an unperturbed fiber by spatially stretching them in such a way as to match the regions of discontinuity of their field components with the position of the perturbed dielectric interfaces.Finally, we find expansion coefficients by satisfying Maxwell's equations (21,22).In the following, we first define an expansion basis and then demonstrate how perturbation theory and a coupled mode theory can be formulated in such a basis.

Expansion basis
We showed previously [1,3,4] that eigen modes of an unperturbed waveguide can be described in the Hermitian Hamiltonian formulation in terms of the transverse fields satisfying an orthogonality condition (6).In the following we will concentrate on the case of fibers, while analogous derivations can be applied to any type of waveguides.
For cylindrically symmetrical fibers the natural choice of coordinates are the polar coordinates.Condition of orthogonality ( 6) for the eigen fields of an unperturbed fiber can be expanded as where in polar coordinates Jacobian of transformation is J 0 (ρ) = ρ.Now, consider a waveguide with a geometry slightly different from an original one.Let (x, y, z) to define a Euclidian coordinate system associated with a perturbed waveguide, and (ρ, θ, s) be a coordinate system corresponding to an unperturbed, cylindrically symmetric fiber, where s is a length of a curve transversed by the center of a waveguide cross-section.Define a coordinate transformation relating such two coordinate systems as (x(ρ, θ, s), y(ρ, θ, s), z(ρ, θ, s)) .Using solutions for the fields of an unperturbed fiber expressed in the coordinates (ρ, θ, s) we form an expansion basis in a Euclidian coordinate system (x, y, z) as following

y, z).
(27) There are several important properties that basis vector fields (27) possess.First, if coordinate transformation is orthogonal ( i ρ = i ρ , i θ = i θ ) then one can show that any two basis fields of ( 27) are orthogonal in a sense of the orthogonality condition (6,26) Moreover, in the case of non-orthogonal transformations, basis fields (27) will be almost orthogonal with an amount of non-orthogonality proportional to the strength of perturbation.
Second, basis fields (27) are continuous in the regions of continuous dielectric of a perturbed fiber profile.In other words, regions of discontinuity in the field components of the basis fields (27), by construction, coincide with the positions of the perturbed dielectric interfaces.
Finally, due to the properties of the generalized coordinate transformations discussed in the previous section i ρ will be strictly normal to the perturbed dielectric interfaces, while i θ will be either strictly tangential (orthogonal coordinate transformation) or almost tangential to the perturbed dielectric interfaces, deviating from strict tangentiality by an amount proportional to the strength of perturbation.As the unperturbed components of the fields E 0 ρ , H 0 ρ and E 0 θ , H 0 θ satisfy the correct boundary conditions on the unperturbed dielectric interfaces the field components ( 27) of ( 27) basis functions will either exactly satisfy the proper boundary conditions on the perturbed interfaces (orthogonal coordinate transformation) or will satisfy them approximately deviating only by an amount proportional to the strength of perturbation.

Coupled mode theory
In the previous section we derived the form (21,22) of Maxwell's equations in the curvilinear coordinates (ρ, θ, s).While seemingly complicated these equations involve an unperturbed, cylindrically symmetric dielectric profile (ρ, θ, s).We look for a solution of Maxwell's equations (21,22) in terms of the basis fields (27) which in (ρ, θ, s) coordinate system are essentially the eigen fields of an unperturbed fiber.We look for a solution of Maxwell's equations (21,22) in curvilinear coordinates (ρ, θ, s) in terms of an expansion into the basis fields (27) entering with corresponding varying along the direction of propagation coefficients C β m (s).Thus, in the covariant coordinates By substitution of expansion ( 28) into ( 21) and manipulation of the resultant equations (see Appendix) we arrive to the following set of equations 6), M is a matrix of coupling elements given by ) where integration is performed over an unperturbed fiber profile and non-zero elements of 6 × 6 matrix M are  Equations (28) present a system of first order linear coupled differential equations with respect to a vector of expansion coefficients C β m (s).Boundary conditions such as a modal content of an incoming and an outgoing signal defines a boundary value problem that can be further solved numerically.
Note that the matrix of coupling elements M is Hermitian and, thus . In (30) we can further perform an integration of matrix M over θ, M exp i(m − m)θdθ as such an integration is decoupled from the integration over the fields which are the functions of ρ only.
To summarize, given a geometric perturbation of a fiber profile, we characterize such a perturbation in terms of a coordinate transformation (x(ρ, θ, s), y(ρ, θ, s), z(ρ, θ, s)) where an original cylindrically symmetric fiber profile in coordinates (ρ, θ, s) is mapped onto a perturbed fiber profile in Euclidian coordinates (x, y, z).Such a coordinate transformation defines a metric g ij .Solving propagation equation governed by the matrix of coupling elements (29,30,31) and applying the proper boundary conditions one finds the transmission and reflection of the incident radiation into the modes of the outgoing and incoming waveguides to any order of accuracy.
Moreover, (29) allows a perturbative interpretation.As a metric of a slightly perturbed coordinate system in only slightly different from the metric of an unperturbed coordinate system that will naturally introduce a small parameter for small geometrical perturbations of fiber profile.Exactly the same perturbative expansion for the values of the modal coefficients due to scattering as in (11) holds in our case, where instead of operator ∆A one would use an operator M − β 0 , where β 0 is a matrix of propagation constants of unperturbed modes.

Study of convergence of the coupled mode theory
We test our theory on a set of geometrical variations in the all-dielectric high indexcontrast waveguides.First we are going to consider long and short tapers and a limiting case of an abrupt junction between two waveguides.Next, we will consider strong Bragg gratings.Finally, bends will be addressed.
An underlying waveguide used in our study is a single core high index-contrast waveguide of radius R core = 1.525a, n core = 3, n clad = 1 surrounded by a metallic jacket of R metal = 4a.Frequency of operation is fixed and equal ω = 0.2 2πc a , propagation constant is measured in units 2π a , where a defines the scale and can be chosen at will.There are altogether 10 -m = 0, 16 -m = ±1, 12 -m = ±2 and 8 -m = ±3 guided modes in this fiber.Metallic boundary conditions due to a waveguide outer jacket eliminate the radiation continuum leaving a discrete spectrum of guided (propagation constant is pure real) and evanescent modes (propagation constant is pure imaginary).We will refer further to this waveguide as a "basic" waveguide.
First, we consider a linear taper variation of a waveguide profile where all the interfaces (including an outside metal jacket) are scaled by the same amount along the direction of propagation (see figure 1b)).A corresponding coordinate mapping will be that of (2) with f (s) = δ s L , where L is a length of a taper, and δ is its strength.Coupling elements can be calculated from (30) and with a selection rule ∆m = 0 they are For a strong taper variation of δ = 100% and a boundary condition of all the incoming radiation coming in the fundamental m = 1 mode we solve (29) for the power in the transmitted and reflected modes propagating in the incoming and outgoing waveguides.We assume propagation from a larger waveguide into a smaller "basic" waveguide.Expansion basis consisted of 8 -m = 1 guided and 42 -m = 1 evanescent modes of an outgoing "basic" waveguide.On figure 2a) we plot in dotted lines a transmitted power T 1 in the fundamental m = 1 mode, along with the transmitted powers T 2 , T 3 in the second and third m = 1 parasitic modes as a function of the taper length L, and calculated by our coupled mode theory .For comparison, results of an asymptotically exact transfer matrix based CAMFR code [20] are presented on the same figure in circles exhibiting an excellent agreement with our results.When taper length goes to zero we retrieve the case of an abrupt junction between different waveguides.In the case when taper length is large enough we are in the regime of a standard taper with most of the transmitted power staying in the same mode.
On figure 2b), we present convergence of the transmitted and reflected coefficients for a linear taper of length L = 10a.Solid lines correspond to the relative errors in the transmission coefficients while dotted lines correspond to the relative errors in the reflected coefficients.Note a very robust convergence of our coupled mode theory with respect to the number of modes in the expansion basis.Errors in the coefficients exhibit faster than quadratic convergence.For comparison, standard coupled mode theory dealing with perturbed high index contrast discontinuous interfaces either does not converges at all or exhibit a very slow linear convergence of errors [2,4].We also observe that in all the regimes of short and long taper lengths our coupled mode theory gives the correct values of the main transmission and reflection coefficients to a percent error with less then 40 basis functions.
Next, we consider a strong index-contrast sinusoidal Bragg grating variation of a waveguide profile where all the interfaces (including an outside metal jacket) are scaled by the same amount along the direction of propagation (see figure 1c)).A corresponding coordinate mapping will be that of (2) with f (s) = δSin(2π s Λ ), where Λ is a period of the grating, and δ is its strength.We assume a starting profile of a "basic" waveguide at s = 0. Coupling elements can be calculated exactly from (30) and with a selection rule ∆m = 0 they are radiation coming in the fundamental m = 1 mode we solve (29) for the power in the transmitted and reflected modes propagating in the incoming and outgoing waveguides.We choose the Bragg grating period Λ to match a beat length between the fundamental (transmission coefficient T 1 ) and a first excited mode (transmission coefficient T 2 ), so in the limit of long gratings all the power will be transferred to the excited mode.Expansion basis consisted of 8 -m = 1 guided and 42 -m = 1 evanescent modes of an incoming "basic" waveguide.On figure 3a) we plot in crosses the transmitted powers T 2 , T 3 , T 4 in the second, third and forth m = 1 modes for the grating lengths [ Λ 2 , 3Λ] in the Λ 2 increments, as calculated by our coupled mode theory.In this geometry the incoming and outgoing waveguides are the same.For comparison, results of an asymptotically exact transfer matrix based CAMFR code [20] are presented on the same figure in circles exhibiting an excellent agreement with our results.When grating length is increased the power transfer to the first excited mode is monotonically increased as expected.On figure 3b), we present convergence of the transmitted and reflected coefficients for a sinusoidal taper of length L = Λ 2 .Solid lines correspond to the relative errors in the transmission coefficients.Note a greater than linear convergence of our coupled mode theory with respect to the number of modes in the expansion basis.We also observe that our coupled mode theory gives the correct values of the main transmission and reflection coefficients to a percent error with less then 20 basis functions.
Finally, we derive an expression for the coupling elements due to a uniform bend (see figure 1d)).A corresponding coordinate mapping will be that of (3), where R is the radius of a bend.Coupling elements can be calculated exactly from (30) and with a selection rule ∆m = 1 they are These matrix elements can be showed to correspond exactly to the well known coupling matrix elements for the coupling of the guided modes in uniform bends of metallic waveguides (see [9], for example).This correspondence is expected as in the case of the uniform bends our method is identical to the "method of cross-sections" as detailed in ( [9]).

Conclusion
In this work, we presented a novel form of the coupled mode and perturbation theories to treat general geometric variations of a waveguide profile with an arbitrary dielectric index contrast.This formulation, unlike much previous work, holds equally well for high index contrasts.We demonstrated that for general scaling variations (including tapers and Bragg gratings) of high index-contrast profiles our methods exhibit a faster than linear convergence with the number of expansion modes.To our knowledge, this is the first time coupled mode and perturbation theories are developed to deal with arbitrary profile variations along the direction of propagation in the high index-contrast waveguides.

Appendix: Derivation of coupling elements
We denote a product of two vector components F 1 , F 2 with respect to a general operator Ô as F 1 | ÔF 2 = (F † 1 ÔF 2 )J 0 (ρ)dρdθ ( † defines a transpose conjugate).First, we substitute expansion (28) into (21) multiplying left and right sides of the resultant equations by the respective factors of g 0ρρ and g 0θθ .We then multiply the modified equations (21) on the left by the proper components of the unperturbed fields and integrate over the fiber cross-section to exploit an orthogonality condition satisfied by the unperturbed solutions (26).Modified equations ( 21 34) into it, while B is a normalization matrix (6).
Evaluation of the coupling elements using expression (33) involves integrals over the derivatives of the unperturbed field components.Using the properties of the modes of an unperturbed structure (24) and integrating by parts we can substantially simplify evaluation of the coupling elements in the form of expression (33).In fact, evaluation of coupling elements can be reduced to the evaluation of integrals over the products of the field components sandwiched between a 6 × 6 matrix of analytical functions M given by (30), totally eliminating the derivatives of the field components.

Fig. 1 .
Fig. 1. a) Dielectric profile of a cylindrically symmetric fiber.Concentric dielectric interfaces are characterized by their radii ρ i .b) Scaling variation -linear tapers.Fiber profile remains cylindrically symmetric, while the radii of the dielectric interfaces along the direction of propagation s become ρ i (1 + δ s L ). c) Scaling variation -sinusoidal Bragg gratings.Fiber profile remains cylindrically symmetric, while the radii of the dielectric interfaces along the direction of propagation s become ρ i (1+δSin(2π s Λ )).d) Non-concentric variation.Each dielectric interface stays cylindrically symmetric, while the center line is bent.

Fig. 2 .
Fig. 2. a) Transmitted power T 1 in the fundamental m = 1 mode, along with the transmitted powers T 2 , T 3 in the second and third m = 1 parasitic modes as a function of the taper length L, calculated by our coupled mode theory .Results of an asymptotically exact transfer matrix based CAMFR code are presented in circles.b) Convergence of the errors in the transmitted and reflected coefficients for a taper length of L = 10a as a function of the number of expansion modes.Solid lines correspond to the relative errors in the transmission coefficients while dotted lines correspond to the relative errors in the reflected coefficients.Calculated by our coupled mode theory, errors in the transmission and reflection coefficients exhibit faster than a quadratic convergence.

Fig. 3 .
Fig. 3. a) Transmitted powers T 2 , T 3 , T 4 in the second third and forth m = 1 modes for the grating lengths [ Λ 2, 3Λ] in the Λ 2 increments are plotted in crosses, calculated by our coupled mode theory.In this geometry the incoming and outgoing waveguides are the same.Results of an asymptotically exact transfer matrix based CAMFR code are presented in circles.When grating length is increased the power transfer to the first excited mode is monotonically increased as expected.b) Convergence of the errors in the transmitted and reflected coefficients for a grating of L = Λ 2 as a function of the number of expansion modes.Solid lines correspond to the relative errors in the transmission coefficients.Calculated by our coupled mode theory, errors in the transmission and reflection coefficients exhibit faster than a linear convergence.