## Abstract

This article investigates the evolution of the distribution of propagating modes in a graded-index multimode fibre with random imperfections. These perturbations may be microscopic random bends and ellipticity, or an index of refraction fluctuations, introduced during manufacture. For large-diameter fibres (with respect to a typical electromagnetic wavelength), the number of guided modes becomes large, and so we can regard the discrete set as a continuum in which the mode number is treated as a continuous variable. The theory is applied to the propagation of guided optical modes in round fibres with a parabolic refractive index profile. This continuous model offers the great advantage of elucidating the relevant physical parameters that participate in the mode mixing dynamics, and leads to a numerically tractable problem. Numerical results are presented for the specific problem of random micro-bends, demonstrating that the theoretically predicted behaviour is achieved.

## 1. Introduction

Fibre optics are now ubiquitous in their application to the telecommunications world for the rapid transfer of vast amounts of data and are the backbone of the Internet. Such cables are electromagnetic waveguides that transport optical pulses at wavelengths usually in the near infrared portions of the spectrum. In practice, optical fibres are of small diameters (typically from a few micrometres up to 1 mm), highly flexible and constituted of highly transparent dielectric materials. They are composed of basically three layers as illustrated in figure 1, an inner glass or polymer core that can be of constant or graded refractive index *n*, a glass or polymer cladding layer of constant refractive index *n*_{c} where *n*_{c}<*n*_{0} (*n*_{0} is the maximum index on the axis) for guiding of waves, and an outer layer of acrylate employed for protection of the fibre that plays no role in the effective optical transmission properties of the waveguide. As the field strength of guided waves is extremely small in the protective jacket it is often reasonable, for mathematical simplicity, to assume that the cladding occupies the unbounded domain exterior to the core. An exact description of the optical field inside the fibre may be constructed by writing down the solution of Maxwell's equations as a sum of eigenmodes, including both guided modes (those that have fields which decay exponentially in the cladding) and radiation modes that describe the non-guided field (Snyder & Love 1983). Each of these guided modes carries a portion of the optical power with a different group velocity, and multimode fibres permit many hundreds or even tens of thousands of such modes to propagate. These fibres are finding increasing application in short- and medium-distance networks for communications in buildings, aircraft, trains, cars, etc. Their attractiveness stems from their low cost and from the ease by which they may be spliced together (owing to their large core size) and excited using inexpensive sources, couplers, splitters and connectors.

In practice, random variations of the optical and geometrical properties of fibres from the ideal model are impossible to avoid despite the advanced and sophisticated methods for manufacturing graded-index multimode fibres. These small imperfections may include index of refraction fluctuations intrinsic to the material constituting the fibre core or extrinsic variations such as microscopic random bends and ellipticity caused by external stress, diameter variations and other fibre core defects (Oates & Young 1989; Heimrath *et al.* 1992). Perturbations having a sufficiently high spatial frequency influence the signal propagation as a result of mode coupling; their cumulative effects may become important after a long propagation length so that they need to be taken into account when calculating the power attenuation, the signal distortion and the bandwidth of the fibre. This is particularly relevant to plastic optical fibres for which experimental and theoretical results indicate that the modes are highly coupled (Jiang *et al.* 1997; Garito *et al.* 1998; Golowich *et al.* 2003).

Many aspects of mode coupling, radiation loss and intermodal dispersion are well known and have been documented since the early 1970s (Marcuse 1969*a*,*b*, 1972; Rowe & Young 1972). The most common approach is to derive and solve the coupled power equations describing the evolution of the average optical power carried by the propagating modes; an excellent summary of the theory can be found in the well-known textbook of Marcuse (1991). This standard procedure requires the diagonalization of an *N*×*N* matrix for the *N* guided modes. Thus, when *N* becomes too large, a direct real-time algebraic treatment of the coupled system is ruled out owing to the huge computational overhead. Plastic optical fibres, for instance, having large core diameters of, say, 1 mm typically, can support several hundred thousand modes when operating at near-infrared wavelengths.

In some cases, nearest neighbour coupling predominates over other power transfer mechanisms and, within appropriate limits, the coupling process can be ideally described in terms of a diffusion equation in which the mode number is treated as a continuous variable. This idea originated in the mid-1970s as a method to understand and model the mode coupling process and to allow quantitative predictions to be made. It was, in particular, developed to deal with the specific problem of random bends (Gloge 1973; Marcuse 1973; Olshansky 1975). However, to our knowledge, no attempt has been carried out to quantify the error made in this continuum approximation and, curiously enough, no progress has been made since this period to further develop the original approach; thus, until recently the same diffusion theory has always been the starting point for evaluating mode conversion in step-index multimode fibres (Zubia *et al.* 2003; Djordjevich & Savović 2004).

In recent papers, the authors offered a new contribution to the diffusion approach for TE mode propagation in randomly perturbed slab waveguides (Perrey-Debain & Abrahams 2005, 2007). In the present work, the method is extended to round fibres with a parabolic graded-index profile. The theory presented here, while relying on the following approximations and assumptions, is, we believe, applicable to a wide range of multimode fibres supporting a sufficiently high number of modes. Firstly, the problem is simplified by assuming that the modes of the parabolic index fibre are the same as the modes of an infinitely extended square-law medium. The fibre boundary is included in the description only insofar as requiring that modes interacting with it suffer high losses so that an effective cut-off exists; at cut-off, we demand that the modes do not carry power. Secondly, the total index change within the guiding core region is only a few hundredths, so the propagating modes can be considered essentially as transverse electromagnetic. Finally, the coupling between the modes is taken to be weak and caused by *random imperfections* or small random variations of the refraction index along the fibre axis; thus, changes in the power distribution between modes take place over distances that are very long compared with the wavelength of light.

The structure of the paper is as follows. After reducing the physical model to a mathematical problem in §2, we present the standard coupled mode theory in §3. After a statistical treatment taking into account the random properties of the distribution of the imperfections, the coupled power equation system1 that describes the evolution of the average power carried by the propagating modes is given at the end of the section. In §4, we present the continuum analogue of the discrete coupled system. Using a regular expansion with respect to a small parameter *ϵ*∼*N*^{−1/2} (*N* is the number of modes), it is shown that the coupled power equations can be approximated as a diffusion equation with an error of order (*ϵ*). We show that the problem is self-adjoint and, consequently, using standard orthogonality arguments, the theory gives access to the explicit form for the first-order correction terms; by this means, the error of approximation can easily be reduced to (*ϵ*^{2}). In §5, we discuss the practical implementation of our approach, and we present some numerical results for the problem of random micro-bends, validating the theoretically predicted behaviour. Finally, in §6, we summarize our results and we offer suggestions for future work.

## 2. Problem statement

We aim at studying the propagation of a monochromatic electromagnetic field with frequency *ω*, , in a circular dielectric waveguide of core radius *a* whose graded-index profile *n* is affected by small random perturbations; to be precise, we perturb *n*^{2} by a small amount −*δn*^{2}. The electric field is governed by the time-harmonic Maxwell vector wave equation for a source-free isotropic medium,(2.1)where ** ∇** is the classical nabla operator; is the position vector in Cartesian coordinates;

*Z*is the guide axis coordinate; and the subscript ⊥ refers to the transversal coordinate on the waveguide cross section, . Further, is the vacuum wavenumber, where

*ϵ*

_{0}and

*μ*

_{0}denote the usual free-space dielectric and permeability constants. The refractive index distribution has the parabolic profile(2.2)in the waveguide region and in the infinite cladding, . The constant

*Δ*is the usual profile height parameter and

*n*

_{0}is the maximum index at the centre of the waveguide. We introduce the waveguide parameter(2.3)which is large for multimode fibres. We will consider weakly guiding waveguides (

*Δ*≪1) as typical multimode fibres operate in this limit, and in fact

*Δ*≈0.01. By expressing the transverse electric field as a power series in

*Δ*and treating

*V*and

*Δ*as independent parameters, the vectorial wave equation (2.1) reduces to the scalar equation (Snyder & Love 1983)(2.4)where the superscript 0 refers to the lowest order solution in

*Δ*. The lowest order term for the transverse magnetic field is given by . The total electromagnetic field can be recovered, with an approximation error (

*Δ*), as(2.5)As the perturbations to the homogeneous fibre are very small, it is reasonable to consider the forward scattering approximation to the reduced wave equation (2.4). We therefore seek a solution of the form , where is the transverse polarization vector. This is substituted in (2.4) and is neglected so that, under the appropriate scaling(2.6)the problem can be conveniently formulated as the following Schrödinger-type equation:(2.7)where

*H*is the

*z*-independent Hamiltonian(2.8)The potential

*v*stands for the quadratic well of finite depth , where and

*v*(

**)=**

*x**V*otherwise. We introduce the small parameter

*ϵ*which characterizes the relative amplitude of the perturbation

*δn*

^{2}and we define the normalized perturbation shape function as so that . The

*z*-dependent perturbed Hamiltonian can then be shown to be(2.9)

## 3. Evolution equations

One method for solving (2.7) is to expand the field *Ψ* in the *z*-dependent eigenmode basis of the perturbed operator *H*+*δH*. This is a very tedious task as *δH*(** x**,

*z*) has generally a complicated shape and is a random function of the normalized axis coordinate

*z*. To circumvent this, we consider sufficiently small perturbations and use standard perturbation theory (Landau & Lifshitz 1977; Snyder & Love 1983; Skorobogatiy

*et al.*2003). In this case,

*Ψ*is expanded in the eigenmode basis of the unperturbed operator

*H*as(3.1)The first summation extends over all the discrete spectrum of the guided modes satisfying the eigenvalue problem (note that we identify the mode index

**(in bold italic) with the multi-index (**

*n**ν*,

*μ*), and use them interchangeably)(3.2)Guided modes propagate without attenuation and their propagation constants lie in the range 0<

*β*

_{n}<

*β*

_{c}, where

*β*

_{c}=

*V*/2 is the cut-off wavenumber. They have an evanescent behaviour in the cladding as(3.3)so they have finite transverse energy, . For a very large number of modes (

*V*≫1), the largest permitted eigenvalue below cut-off can be shown to be approximately given by the upper bound, . The integral in (3.1) extends over the continuum spectrum of

*H*corresponding to the radiation modes and the discrete summation extends over their azimuthal dependence (index

*l*). These modes are composed of a combination of inward and outward travelling waves and consequently they do not satisfy the radiation condition at infinity. Nevertheless, since they still behave like as , they can be conveniently normalized so that they satisfy the orthogonality relation (Sammut 1982)(3.4)where the asterisk symbol refers to the complex conjugate; denotes the Kronecker symbol; and is the Dirac delta function. To be consistent with the forward scattering approximation introduced earlier, we may assume that the integration domain in (3.1) can be restricted to a narrow range of propagation constants, . Using orthogonality properties for both guided and radiation modes, equation (2.7) is transformed into the system of ordinary differential equations(3.5)where the coupling coefficients are given by(3.6)and similar formulae hold for the radiation coupling coefficients . It is now well established that, under some additional assumptions, system (3.5) can be averaged over an ensemble of

*N*

_{w}similar waveguide realizations. More precisely, if the guided modes are weakly coupled over a distance which is large compared with the correlation length of the random process

*δH*then the ensemble average power(3.7)carried by mode

**, can be shown to satisfy the master equation system (Marcuse 1991)(3.8)where the transition probability matrix coefficients are given from the spectral density of evaluated at the wavenumber spacing , i.e.(3.9)where the operator stands for the usual Fourier transform. The last equality is a consequence of the Wiener–Khintchine theorem for sufficiently well-behaved stationary random processes. The content of (3.9) is that the averaged transfer of power between modes**

*n***and**

*n***is given by the power spectral density of the process, , evaluated at the spatial frequency**

*m**β*

_{m}−

*β*

_{n}. Therefore, the condition for appreciable coupling between the two modes is that the normalized perturbation must both induce a substantial overlap integral in (3.6) and have spatial frequency support at

*β*

_{m}−

*β*

_{n}. Note that a similar derivation can be found in the context of quantum mechanics (Garnier 1999) and in acoustics (Dozier & Tappert 1978). The positive quantities

*α*

_{n}take into account the averaged coupling between mode

**and the modes of the continuum of the radiation modes. A rigorous analysis of the radiation loss is a difficult task as it requires both an accurate description of the perturbation in the vicinity of the core–cladding interface as well as the mathematical form for the guided and radiation modes close to cut-off. Nevertheless, some simplification will be considered later for fibres supporting a sufficiently large number of guided modes.**

*n*Owing to the symmetry of the transition probability matrix, the solution of (3.8) is explicitly given by(3.10)with initial conditions(3.11)Here, the initial wavefield is given by the illumination of the laser beam at the input of the fibre (launch conditions). Column vectors *U*^{(i)} are the eigenmodes of the real symmetric system(3.12)The transition probability operator ** W** may be written explicitly as(3.13)The modal attenuation matrix

**contains the attenuation factors, , and the diagonal matrix, contains the real eigenvalues in descending order. Owing to the special structure of (3.12), the Gerschgorin discs(3.14)with radius , lie in the left part of the complex plane and all eigenvalues are negative with . Above a sufficiently large distance**

*Γ**z*>

*L*

_{s}, say, the power distribution becomes independent of the initial condition and reaches the steady-state distribution(3.15)The steady-state mode coupling length

*L*

_{s}strongly depends upon the launch condition coefficients . In most practical situations, it is expected that , where (Marcuse 1991). Note that neglecting the radiation loss yields the adiabatic case,

*λ*

_{1}=0, corresponding to the long-distance solution(3.16)This ideal scenario is obviously unrealistic as losses from highest order modes are unavoidable (though mode coupling can be tailored to minimize losses; Marcuse 1976). The purpose of this paper is to evaluate the set of eigenmodes

*U*^{(i)}corresponding to the lowest |

*λ*

_{i}|'s when the number of guided modes is so large that a direct algebraic treatment is out of the question from a computational point of view. Fortunately, the density of guided modes allows the alternative continuum approach to be successful, as we shall now show.

## 4. The continuous model

### (a) Simplification when *V*≫1

Though eigenfunctions *Ψ*_{n} of (3.2) can be represented exactly by Whittaker's functions (Snyder & Love 1983), a scalar eigenvalue problem still needs to be solved. Moreover, these functions are not easy to manipulate. In order to simplify the calculation, the lowest order modes can be treated as solutions of the two-dimensional harmonic oscillator since these modes essentially see an infinite parabolic well. These ideal modes are given by the Cartesian product , , where(4.1)and *H*_{ν} denotes the usual Hermite polynomial. These functions are normalized such that the orthogonality property holds(4.2)The eigenvalues for these appropriate eigensolutions are given by . We make the important note that writing as a Cartesian product in any way does not restrict the form of the waveguide cross section, and the present analysis could be equally carried out in polar coordinates using Laguerre–Gauss functions.

In the quadratic well, , we assume that the perturbation admits the separable form(4.3)where the *Χ*_{j}'s are real-valued zero-mean, stationary and ergodic processes with respect to the normalized waveguide axis coordinate *z* and the 's describe the shape function perturbation corresponding to the *j*th process. Without loss of generality, we may assume that these shape functions are sufficiently regular so that they can be fairly approximated by their truncated Taylor series(4.4)The new coordinate system is introduced above to ensure that the amplitude of the perturbation is independent of the waveguide parameter *V* and remains of order unity everywhere in the waveguide. For the ideal mode (4.1), the evaluation of the coupling coefficients can be simplified by extending the perturbation over the whole waveguide cross section to get(4.5)The lemma 4.1 established previously by the authors (Perrey-Debain & Abrahams 2006) allows an analytical treatment of the overlap integrals (4.5).

*Given positive integers* *, then the following integration formula holds:*(4.6)*where* *and*(4.7)*and* *is the set of indices* *associated with the nested sum*(4.8)

By common convention, the products above take the value unity when the lower limit exceeds the upper, and the notation indicates that *σ* takes the value of *ξ* when it is an integer, or else the sum is zero. Note that the parameter *ν* has been written as an argument of *F* and *G* because we will soon generalize it to take non-integer values. Several other quantities will soon be defined which will also use this convention. Now, let us introduce the separation vector defined as(4.9)The lemma 4.1 shows that the coupling coefficient may be written as(4.10)where(4.11)with and . Thus, the transition probability matrix from mode ** n** to mode

**+**

*n***has the analytical form(4.12)in which**

*s**w*

_{s}is the polynomial series defined over the entire

*continuous*plane ,(4.13)and stands for the correlation spectrum(4.14)Equation (4.12) is the most important result of the paper. It reveals the underlying structure of the transition probability matrix which is perfectly suited for a continuum analysis approach. Expression (4.12) holds under the assumption that the modes of the parabolic fibre are essentially the same as the modes of an infinitely extended square-law medium. This is a reasonable approximation for highly multimode fibres, i.e.

*V*≫1, as the finiteness of the fibre will have little effect on most of the guided modes. Only the highest order guided modes close to the cut-off condition carry non-negligible power at the fibre boundary, and therefore suffer from very high losses. To simplify the analysis, we will assume that whenever , which means that the highest order modes above cut-off carry no power(4.15)whereas modes below cut-off are assumed to be conveniently described by the harmonic oscillator basis (4.1). The cut-off condition (4.15) was originally introduced by Marcuse (1973) and was recently found to be in agreement with experimental measurements carried out by Golowich

*et al.*(2003). By virtue of (4.7), the coupling between guided modes is

*local*, i.e.(4.16)where and are the maximum exponent in the polynomial expansion (4.4). Thus, the cut-off condition (4.15) needs to be imposed only in a narrow strip above the cut-off line to ensure the transfer of modal energy between the modes below and above cut-off. Finally, since there is no coupling with negative indices, we have necessarily(4.17)In fact, this condition is automatically satisfied in (4.12). Consider, for instance, a ‘virtual’ mode

**located in the no-coupling region as in figure 2. The point is such that , which is equivalent to . Thus, by definition and (4.17) holds.**

*m*We consider an infinitely differentiable function of the continuous variable , say, which interpolates the value of the discrete eigenmode , i.e. for all corresponding to the guided modes in the fibre. Note, however, that, in the cut-off strip, we have simply . The interpolating Lagrangian polynomial is one such function, for instance, which satisfies the above properties. We can now exploit the analytical form (4.12) to determine the continuous analogue of the eigensystem (3.12). Direct application of Taylor's theorem for infinitely differentiable functions yields the following lemma 4.2.

*Let* *be the infinitely differentiable function introduced above. Then the transition probability operator may be written formally as*(4.18)*where* *stands for the directional derivative* . *For non-pathological functions, that is, functions which have a regular Taylor series expansion, the residual term* *is given by the infinite series*(4.19)

By virtue of (4.16) and (4.17), the finite sum in (3.13) can be extended over the whole set of integers ^{2} and, using the identity (4.12), we can write(4.20)Expanding the (infinite) Taylor series of in the neighbourhood of the mid-point gives(4.21)Since *w*_{s} is a polynomial series, we can expand the r.h.s. term of (4.20) in its infinite Taylor series around the point ** n**. This completes the proof. ▪

The reader's attention is drawn to the fact that the leading term on the r.h.s. of (4.18) is just the diffusion operator acting on . We comment on this further below.

### (b) Asymptotic series expansion

We are now interested in the eigenmode solutions of (3.12) when the number of guided modes tends to infinity. Thus, we introduce the small parameter *ϵ*=1/*V*, and we write every term in (3.12) as a series of powers of *ϵ*. Equating coefficients at each order then yields a hierarchy of equations. Let us first observe that, by introducing the normalized continuous variables , where , the function *w*_{s} admits the regular series expansion(4.22)where(4.23)are all polynomial series with respect to the variables *u* and *v*. The leading order term is found to be(4.24)where(4.25)and *C*_{pη} is the number of summations in the nested sum (4.8), i.e.(4.26)if and otherwise. Thus, terms in the summation (4.24) are non-zero provided *p* and *p*′ (respectively *q* and *q*′) are of the same parity, and therefore the square root term never leads to non-integer powers. Hence, (4.22) suggests that we seek a solution in the form of a simple regular perturbation expansion(4.27)and(4.28)Similarly, we introduce a regular function *γ* such that *γ*(** n**)=

*α*

_{n}for all modes below cut-off

*ν*+

*μ*<

*β*

_{c}. Furthermore, we may assume that

*γ*admits the regular expansion(4.29)By introducing these series into the eigensystem (3.12), and using the continuous form for the transition probability matrix given in lemma 4.2, we obtain (at leading and next order)(4.30)(4.31)where and

^{α}is the symmetric matrix(4.32)which can be interpreted as a diffusion tensor controlling the average transfer of modal power at the ‘mode number’ . At higher orders, the equations are more complicated since they involve higher derivatives in (4.19).

The two equations (4.30) and (4.31) must be solved over the triangular domain with edges , and . The leading order solution must satisfy the ‘cut-off’ boundary conditions(4.33)Conditions on emerge naturally after realizing that there is no transfer of modal energy from negative indices, which implies that(4.34)where is the unit vector normal to the boundary. This condition is in fact the continuum analogue of (4.17) and is automatically satisfied for any sufficiently regular solution. With the conditions (4.33) and (4.34), the leading order problem forms a self-adjoint boundary-value problem on *L*^{2}(*Ω*) and, consequently, the following orthogonality property holds:(4.35)where ‖.‖ stands for the usual *L*^{2}(*Ω*) norm. Identity (4.35) is valid provided the set of eigenvalues {*λ*_{i,0}} are all distinct. To avoid unnecessary complications, we will assume in the following that this condition holds. Note that the orthogonality of the eigenfunctions is associated with the orthogonality of the eigenvectors(4.36)

### (c) Leading order solution and first-order correction

Let us now define vectors , as the discrete versions of their continuous counterparts: for all guided modes. Rewriting the perturbation expansion (4.27) in its vectorial form gives(4.37)According to (4.36), the norm of the leading order solution ‖*U*_{i,0}‖ must be chosen such that . This can be easily shown to be satisfied by simply taking . The first-order correction is explicitly obtained by expanding *U*_{i,1} in the orthogonal basis *U*_{i,0}, i.e.(4.38)Substituting (4.38) into (4.31) and using orthogonality properties yields(4.39)and(4.40)The diagonal correction terms *a*_{ii} stem from the discrete normalization (4.36). To first order, this is equivalent to(4.41)It can be shown (see appendix A) that this condition is satisfied if(4.42)This completes evaluation of the first two terms of in (4.37).

## 5. Some applications

In this section, we wish to simulate typical deformations observed in optical fibres as illustrated in figure 3. In the core region (), the unperturbed waveguide has the normalized index profile of the parabolic type: . The perturbed waveguide has the following index profile:(5.1)where (*δx*, *δy*) are the coordinates of the origin of the perturbed fibre; and *δa* and *δb* refer to the small ellipticity with respect to the horizontal and vertical axes, respectively. All these quantities are random functions of the guide axial coordinate *z*. To follow the notation introduced earlier in §4*a*, we put . Away from the core–cladding interface, straightforward calculation shows that, to first order,(5.2)By assuming that all the *Χ*_{j}'s are independent processes, the leading order diffusion tensor becomes(5.3)and the first-order tensor is given by(5.4)

### (a) Micro-bending

There is no analytical solution to the eigenvalue problem defined in (4.30) and numerical methods have to be used. Nevertheless, in the particular case of random eccentricity for which the random processes *Χ*_{1}(*z*) and *Χ*_{2}(*z*) are identical, the eigenvalue problem has the relatively simple form (we ignore modal attenuation)(5.5)where for the sake of clarity we set . The associated original discrete system (3.12) can be written explicitly as(5.6)In this example, mode coupling occurs only between adjacent modes and so the representation of the power-coupling process in terms of a diffusion equation is clearly revealed since (5.6) is nothing else but the finite-difference discretization of (5.5). Equation (5.5) has the remarkable property of being separable with respect to the new variables *θ*=*v*/*u* and *ζ*=*u*+*v* (see appendix B). It turns out that the unique set of regular eigenfunctions is , where(5.7)and(5.8)where the function is the Legendre polynomial of degree *k*′; is the Bessel function of the first kind of order 2*k*′+1; and is the *k*th zero of . The correspondence is given by ordering the associated eigenvalues in descending order(5.9)The first-order diffusion tensor ^{1} is simply the identity (see (5.4)), which yields(5.10)The proof for the second equality is rather tedious and is omitted for brevity. Unfortunately, there is no analytical form for the correction coefficients given by formula (4.39) and numerical integration has to be performed. In table 1, we display the values of the first three eigenvalues calculated from the original discrete system (5.6) and from the first-order approximation . The number of digits of accuracy given by the continuous model is in agreement with the expected *ϵ*^{2} law. This is clearly confirmed in figure 4 where the quadratic errors (in per cent) for the leading order solutions, and for the first-order solutions are plotted against the waveguide parameter *V*. Note that the first-order correction vectors are computed with the first 30 terms in the infinite sum (4.38). Figure 5 shows the graph with the first three eigenfunctions (5.7) for the micro-bending problem and the corresponding discrete modal distribution from the eigenvalue problem (5.6). Note that the second eigenmode is antisymmetric with respect to the line *u*=*v* and therefore is unlikely to be excited under standard launch conditions (i.e. the laser spot is centred on the waveguide axis and *Ψ*(** x**, 0) is only a function of |

**|).**

*x*### (b) General case: the Rayleigh–Ritz method

In a more general situation such as that given by the diffusion tensor (5.3), eigenfunctions of (4.30) can be numerically recovered using the Rayleigh–Ritz method. For the sake of notational simplicity, we define the symmetric bilinear form(5.11)for any sufficiently regular functions *ϕ* and *ϕ*′ satisfying the Dirichlet condition (4.33). The positiveness of is a direct consequence of the positiveness of the quantity on *Ω*. Indeed, it is then clear that . Furthermore, we have the following identity:(5.12)Consequently, for any non-trivial solution and we can define the energy functional . (We note in passing that, if the cut-off condition is neglected, *ϕ*=const. is an eigensolution corresponding to the adiabatic case (3.16).) Solutions of (4.30) are stationary points of the energy functional subject to the normalization constraint . The set of functions forms an orthogonal system in *L*^{2}(*Ω*) satisfying the cut-off boundary condition at *∂Ω*_{0} and can therefore serve as a natural basis for an approximate solution , that is, we consider the truncated generalized Fourier expansion(5.13)Following standard variational techniques, stationary points are reached at the approximate eigenvalue where the expansion coefficients satisfy the matrix eigenvalue problem(5.14)The symmetry of the bilinear form implies that the basis set is orthogonal in *L*^{2}(*Ω*). Furthermore, the true solution is recovered by taking the limit . Similarly, and, from the Rayleigh–Ritz principle, approximate eigenvalues are upper bounds for the true eigenvalues of the infinite-dimensional problem. The convergence of the method depends upon the properties of the perturbation such as its shape and its power spectrum. In most cases of practical interest, the micro-bending is the predominant perturbation and the truncated series (5.13) is expected to converge very rapidly with *K*. It is found that the first eigenfunctions (which are the quantities of interest) can be computed with excellent accuracy at very low computational cost for a matrix size that does not exceed *K*=100. For examples computed in the one-dimensional case, we refer the interested reader to a recent article (Perrey-Debain & Abrahams 2007).

## 6. Conclusion

In this paper, we have analysed the evolution of the modal power distribution of the electromagnetic wave field as it propagates along a randomly deformed multimode fibre with a quadratic index profile. We showed that, for fibres supporting a sufficiently large number of guided modes, the mode coupling mechanism can be ideally described as a continuous diffusion equation. Even for a moderate number of modes, the regular expansion method allows us to obtain very accurate solutions when first-order correction terms are taken into account. This continuous model offers the great advantage of elucidating the relevant physical parameters that participate in the mode mixing dynamics and leads to a numerically tractable and efficient problem. Work is ongoing to apply the theory to study the effect of typical imperfections introduced during the manufacturing process. One direction of particular interest is to adapt the time-dependent theory to the diffusive regime. The ultimate goal is to identify relationships between (i) the random perturbations, (ii) the associated diffusion tensor, and (iii) the characteristic bandwidth as well as the loss in the fibre. By studying these connections, one can hope to find ways of tailoring the perturbations to further reduce the loss in fibres while still maintaining the high bandwidth.

## Acknowledgements

The research described herein was carried as part of the activities of the Smith Institute Knowledge Transfer Network. The authors are most grateful to the UK Engineering and Physical Sciences Research Council (EPSRC) and Photon Design (Oxford, UK) for the award of a postdoctoral grant for E.P.-D., and to Dr David Allwright and Dr John Ockendon (OCIAM, University of Oxford, UK), and Dr Dominic Gallagher and Dr Tom Felici (Photon Design) for their substantial input to the project.

## Footnotes

↵Commonly referred to as master equations by the mathematical physics community.

- Received July 20, 2007.
- Accepted January 2, 2008.

- © 2008 The Royal Society