Comprehensive analytical model to characterize randomness in optical waveguides

Junhe Zhou; Philippe Gallion

doi:10.1364/OE.24.006825

1. Introduction

Optical waveguides, including the planar slab optical waveguides, rib-ridge waveguide, optical fibers etc. are the fundamental structures in the optical communication systems. The waveguides are categorized as single mode optical waveguides and multimode optical waveguides. Single mode waveguides have received much attention during the past decades, because of the lower modal dispersion, while multimode waveguides have been focused on recently due to the emerging mode division multiplexing technique.

Either single mode waveguides or multimode waveguides suffer from random refractive index variations. During the fabrication of the optical waveguides, the imperfect surface caused by etching resolution limit [1,2 ] and the irregularity of the materials inside the waveguides, will lead to random index variations either at the edges or within the cores of the waveguides. Such random index variations will introduce additional mode loss and mode coupling [1,2 ]. As indicated by the literatures [1,2 ], functional devices based on optical waveguides will have a degraded performance due to the effect.

The impact of the random index variations inside optical waveguides can be studied by the coupled mode theory (CMT) [3]. Marcuse [4–6 ] derived a set of coupled equations for the power evolution of the modes based on the CMT. His pioneering works has opened a gate for the analysis of the random index variations inside optical waveguides and has inspired abundant subsequent works [1,2,7–11 ].

Despite the ingenious works [1,2,4–12 ] led by Marcuse, there are a few unanswered questions. First of all, there has been no comprehensive theoretical framework, which brings difficulties to study the problem in a deeper sense. Although Marcuse has brought the coupled power equation, he did not provide the exact coupled equations for the higher order statistics, such as the power fluctuations and power correlation coefficients. In his book [12], Marcuse derived an approximate coupled power fluctuation equation and solved it by assuming that the coupling between the modes does not contribute to the correlation between them (i. e. the mode powers are independent variables), and the correlation maintains as 0. As pointed out by Marcuse himself [12], this is only an approximate assumption and is obviously not true in a lossless two-mode waveguide, where the powers of the two modes are perfectly correlated. Therefore, the power deviation calculated by the equation in [12] becomes inaccurate when the mode coupling is strong and the propagation distance is long. Furthermore, the correlation between the modes are always assumed to be zero and have not been studied in Marcuse' book [12]. The lack of study on these important statistical parameters brings difficulties during the evaluation of the system performance, such as the power sensitivity and overall system capacity. For example, a batch of multimode waveguides are fabricated by the same process, and therefore they should have similar edge roughness. However, due to the random nature, the mode coupling may induce different output mode power distributions for different waveguides when the inputs are the same. A more significant example is from the emerging mode division multiplexing technology. When the multimode waveguides are used for mode division multiplexing, random mode coupling evaluation will be essential to characterize the optical multiple input multiple output channel [13–16 ]. In propagation links, the optical waveguide index might be impacted by the temperature and strain variations [17]. For example, the temperature sensitivity of SiO₂ is 1.6 × 10⁻⁵°C [17], while temperature sensitivity of Ge, the widely used material doped in the optical communication fibers, is 2 × 10⁻³°C [17]. These values will be significant for a fiber link with hundreds of kilometers' length. The actual index temperature sensitivity for the optical fibers or waveguides depends on the material composition of the optical fibers/waveguides. Just like the case of free space optical communication links [18], the random changes of index from time to time will also bring power variations for each mode at the output. In free space communication links, such random power fluctuation is measured by the scintillation index, and it scales with the propagation distance [18] (because it scales with the Rytov variance and the Rytov variance scales with the propagation distance). Therefore, it can be expected that in the multimode fiber/waveguide transmission systems, the varying output power at each mode will bring more stringent sensitivity (defined as minimum input power required to produce the desired signal to noise ratio) requirement for the detectors at the receiver side. The degradation of the sensitivity can be measured by the standard deviation of the mode power. Furthermore, the correlations between different modes due to the random coupling will impact final channel covariance matrix, which is related to the channel capacity [16]. Such critical information can only be calculated from the higher order statistics such as the variance of the mode power, or the covariance between the powers of different modes.

To tackle these unresolved problems, we have conducted the following works: 1) A set of stochastic differential equations (SDEs) is proposed to describe the wave propagation and coupling inside the optical waveguides with random index variations; 2) A theoretical framework is established to find the ordinary differential equations (ODEs) of the higher order statistics of the modal amplitudes; 3) ODEs for the second order statistics (the optical power) and the fourth order statistics (the power fluctuations and the power correlation coefficients) are derived. It is found that the ODEs for the modal power evolution coincide with the coupled power model proposed by Marcuse. Analytical solutions to these ODEs are provided; 4) Detailed Mant-Carlo simulations are conducted based on the SDEs. The results are compared with the ones predicted by the ODEs and excellent agreement is achieved.

2. Mathematical model

2.1 Coupled mode theory inside optical waveguides

Without loss of generality, we start with the scalar Helmholtz equation [3–6 ] is to derive the coupled mode equation,

\nabla^{2} φ + k^{2} φ = 0

where k stands for the wave number in the waveguides. The theory will be extended to the vector case [3] afterwards.

Assuming that there is a perturbation in the refractive index, we have

\begin{array}{l} \nabla^{2} φ + δ k^{2} φ + 〈 k^{2} 〉 φ = 0 \\ δ k^{2} = k^{2} - 〈 k^{2} 〉 \end{array}

where < > stands for the average value. The field inside optical waveguide can be expanded by its orthogonal modes (including the leaky modes) of the unperturbed waveguides as [3]

φ = \sum a_{m} ϕ_{m} e^{- j β_{m} z}

where, a_m stands for the amplitude, β_m the propagation constant, and ϕ_m the modal field of the modes. Here, Σ not only stands for the summation for the guided modes, but also the integration for the leaky modes which have continuous varying propagation constants. These modes should fulfill

\nabla_{⊥}^{2} ϕ_{m} + (〈 k^{2} 〉 - β_{m}^{2}) φ = 0

Substituting Eq. (3-4) into Eq. (2), and neglecting the second order derivative of a_m with respect to z, one has [3]

- \sum ϕ_{m} 2 j β_{m} \frac{\partial a_{m} (z)}{\partial z} e^{- j β_{m} z} + δ k^{2} \sum a_{m'} (z) ϕ_{m'} e^{- j β_{m} z} = 0

Due to the orthogonality between the modes, it can be derived that [3]

\frac{d a_{m} (z)}{d z} = \sum C_{m m'} (z) e^{j (β_{m} - β_{m'}) z} a_{m'} (z)

where

C_{m m'} (z) = \frac{1}{2 j β_{m}} \frac{\iint_{x y} ϕ_{m}^{*} ϕ_{m'} δ k^{2} d S}{\iint_{x y} {| ϕ_{m} |}^{2} d S}

Assuming the modes has been normalized such as [3–6 ]

\iint_{x y} ϕ_{m} ϕ_{n}^{*} d S = \frac{2 ω μ}{β_{m}} δ_{m n}

one has

\begin{array}{l} C_{m m'} (z) = \frac{k_{0}^{2}}{j 4 ω μ} \iint_{x y} ϕ_{m}^{*} ϕ_{m'} δ n^{2} d S \\ δ k^{2} = δ n^{2} k_{0}^{2} \\ δ n^{2} = n^{2} - 〈 n^{2} 〉 \end{array}

where k ₀ is the wave number in the free space.

2.2 The Stochastic differential equation model

Assuming δn ² can be decomposed into two random independent processes in the x-y plane and the z direction, we have

\begin{array}{l} 〈 δ n^{2} (x, y, z) δ n^{2} (x', y', z') 〉 = R_{x y} (x - x', y - y') η (z - z') \\ 〈 C_{m m'} (z) C_{m m'}^{*} (z') 〉 = κ_{m m'} η (z - z') \end{array}

where

κ_{m m'} = \frac{k_{0}^{4}}{16 ω^{2} μ^{2}} \iint_{x y} \iint_{x' y'} ϕ_{m}^{*} (x, y) ϕ_{m} (x', y') ϕ_{m'} (x, y) ϕ_{m'}^{*} (x', y') R (x - x', y - y') d S d S'

In the special case that the randomness only occurs in the z direction while it is deterministic on the x-y plane (e.g the rough side wall), the coupling coefficient C_mm'(z) can be modeled as a constant multiplies a random process f(z) in the z direction.

\begin{array}{l} C_{m m'} (z) = ς_{m m'} f (z) \\ 〈 f (z) f (z') 〉 = η (z - z') \end{array}

where

ς_{m m'}

is determined by the position of the rough surface and the modal distributions of mode m and mode m', as indicated by Eq. (9). For example, when the rough surface is at the position of x = a,

δ n^{2} = δ (x - a) f (z)

it can be calculated as

ς_{m m'} = \frac{k_{0}^{2}}{j 4 ω μ} \int_{- \infty}^{\infty} ϕ_{m}^{*} (a, y) ϕ_{m'} (a, y) d y

With the introduction of the above random variables, Eq. (6) becomes a set of SDEs. The coupled mode equation should be regarded as Stratonovich sense, which uses the value in the middle of the propagation step during the numerical integration. It should be converted into Ito sense in order to take the average on both sides of the equation and let the random variable to be averaged to zero. The conversion is accomplished as follows

\begin{array}{l} d a_{m} (z) = a_{m} (z + d z) - a_{m} (z) \\ = \sum_{m'} a_{m'} (z) \int_{z}^{z + d z} C_{m m'} e^{j (β_{m} - β_{m'}) z'} d z' + \sum_{m'} \frac{d a_{m'} (z)}{2} \int_{z}^{z + d z} C_{m m'} e^{j (β_{m'} - β_{m}) z'} d z' \\ = \sum_{m'} a_{m'} (z) \int_{z}^{z + d z} C_{m m'} e^{j (β_{m} - β_{m'}) z'} d z' + \frac{1}{2} (\sum_{m'} \int_{z}^{z + d z} C_{m m'} e^{j (β_{m} - β_{m'}) z'} d z' \int_{z}^{z + d z} C_{m m'}^{*} e^{- j (β_{m} - β_{m'}) z'} d z') a_{m} (z) + O (d z^{3 / 2}) \end{array}

The detailed derivation of Eq. (14) is provided in the appendix. Equation (14) is quite different from the usual coupled mode equation which has been adopted by Marcuse [4–6 ]. It indicates that the changing ratios of the amplitudes are not only related to the amplitudes of the modes, but also the changing ratios themselves. The second term on the right hand side (RHS) does not vanish because the expectation of this term will not be zero. Equation (14) is the SDE which serves as the theoretical basis of the following analysis.

2.3 Energy conservation

Firstly, the validity of Eq. (14) can be examined by deriving the energy conservation law. The total power of the optical field is the summation of the power of each mode. We denote the total power as P _t

〈 P_{t} 〉 = 〈 a^{H} a 〉

where a is the vector whose elements are a_i, and H denote Hermitian transpose. Since the coupling coefficients C_mm _' forms an anti-Hermitian matrix, it can be verified from Eq. (14) that

\frac{d {〈 P 〉}_{t}}{d z} = 0

So the total power P _t is a constant along the z axis. However, the mode coupling induced energy transfer will take place between the modes and the power fluctuations will continue to grow. It should be noted that the power conversion law is only valid when the leaky modes are taken into account. If only the guided modes are considered, the total power will be attenuated due to the coupling to the leaky modes.

2.4 Coupled power equation

The power of each mode is [4]

P_{m} (z) = 〈 a_{m} (z) a_{m}^{*} (z) 〉

Its differential is

d (a_{m} (z) a_{m}^{*} (z)) = d (a_{m} (z)) a_{m}^{*} (z) + a_{m} (z) d (a_{m}^{*} (z)) + d a_{m} (z) d a_{m}^{*} (z)

It can be noticed that the differential from the SDE is different from the one from the ODE. The second order product of the differentials is included in Eq. (18), because it will not vanish after averaging. It can be assumed that a_m and δn² are independent to each other, which is reasonable since Ito sense SDE has been adopted. Therefore, substituting Eq. (14) into Eq. (18), and using the property that the average of the product of δn² and a_m is 0, one has

\begin{array}{l} d P_{m} (z) = \sum P_{m'} (z) 〈 \int_{z}^{z + d z} C_{m m'} e^{j (β_{m} - β_{m'}) z'} d z' \int_{z}^{z + d z} C_{m m'}^{*} e^{- j (β_{m} - β_{m'}) z''} d z'' 〉 \\ - \sum P_{m} (z) 〈 \int_{z}^{z + d z} C_{m m'} e^{j (β_{m} - β_{m'}) z'} d z' \int_{z}^{z + d z} C_{m m'}^{*} e^{- j (β_{m} - β_{m'}) z''} d z'' 〉 \end{array}

The first term on the RHS in Eq. (19) comes from the second order differential in Eq. (18), while the second term on the RHS comes from the first order differentials. The term inside the bracket can be further calculated as

〈 \int_{z}^{z + d z} C_{m m'} e^{j (β_{m} - β_{m'}) z'} d z' \int_{z}^{z + d z} C_{m m'}^{*} e^{- j (β_{m} - β_{m'}) z''} d z'' 〉 = d z κ_{m m'} \int_{- d z}^{d z} η (z') e^{j (β_{m} - β_{m'}) z'} d z'

where

κ_{m m'} = \frac{k_{0}^{4}}{16 ω^{2} μ^{2}} \iint_{x y} \iint_{x' y'} ϕ_{m}^{*} (x, y) ϕ_{m'} (x', y') R (x - x', y - y') d S d S'

If the differential dz is much longer than the correlation length of the random index variation in the z direction (in this case, the propagation distance should be much longer than the correlation length), which is the assumption adopted by Marcuse [4–6 ] and his followers [1], the results in Eq. (20) can be formulated as,

\int_{- d z}^{d z} η (z') e^{j (β_{m} - β_{m'}) z'} d z' \approx F T (η) (β_{m} - β_{m'})

where FT stands for the Fourier transform. And we have

\frac{d P_{m}}{d z} = \sum κ_{m m'} F T (η) (β_{m} - β_{m'}) (P_{m'} (z) - P_{m} (z))

When only the surface roughness is considered, the randomness in the x-y plane vanishes, and Eq. (22) becomes the coupled power equation proposed by Marcuse [4–6 ], which is derived under another theoretical framework [4–6 ]. Therefore, the validity of the proposed SDE is demonstrated. Furthermore, Eq. (22) considers the randomness in the x-y plane, which is not included in the Marcuse' model. Therefore, Eq. (22) cannot only be used for the analysis of the surface roughness impact for the optical waveguides, but also used for the analysis of the beam propagation inside random media, e. g. beam propagation in the free space with the presence of the atmosphere turbulence.

2.5 Coupled equation for the higher order statistics

Other than the second order statistics, i. e. the modal powers, higher order statistics of the amplitudes becomes more and more important. With the assistance of the SDE proposed in this paper, it is possible to derive the coupled ordinary differential equations for the power variations as well as the power correlation coefficients.

The second order statistics of the optical power is

〈 P_{m}^{2} (z) 〉 = 〈 a_{m}^{2} (z) a_{m} {^{*}}^{2} (z) 〉

Its differential is

\begin{array}{l} d P_{m}^{2} (z) = 2 a_{m} (z) a_{m} {^{*}}^{2} (z) d a_{m} (z) + 2 a_{m}^{*} (z) a_{m}^{2} (z) d a_{m}^{*} (z) \\ + 2 a_{m} (z) 2 a_{m}^{*} (z) d a_{m} (z) d a_{m}^{*} (z) + a_{m} {^{*}}^{2} (z) d a_{m} (z) d a_{m} (z) \\ + a_{m}^{2} (z) d a_{m}^{*} (z) d a_{m}^{*} (z) \end{array}

Substituting Eq. (14) into Eq. (24) and averaging on both sides of the Eq., one has

\begin{array}{l} \frac{d 〈 P_{m}^{2} (z) 〉}{d z} = - 2 〈 P_{m} {(z)}^{2} 〉 \sum_{m'} κ_{m m'} F T (η) (β_{m} - β_{m'}) \\ + 4 \sum_{m'} 〈 P_{m} (z) P_{m'} (z) 〉 κ_{m m'} F T (η) (β_{m} - β_{m'}) - 2 κ_{m m} F T (η) (0) 〈 P_{m} {(z)}^{2} 〉 \end{array}

The self-coupling term (3rd term on the RHS on Eq. (22)) is included in Eq. (22) while the Marcuse' approximated coupled power fluctuation equation neglects this term [12]. It can be not neglected because it is a random coupling term and cannot be absorbed by the deterministic propagation constant. However, this term will automatically vanish when it combines with the second term on the RHS. Furthermore, the power correlation term has appeared in Eq. (25), since it is not necessary that P_m and P_m' are independent to each other, the power correlation term cannot be evaluated by Eq. (22) (In Marcuse' approximated solution, P_m and P_m' are assumed to be independent). Therefore, the coupled equation for the power correlation term should be derived from the following differential equation.

\begin{array}{l} d P_{m} (z) P_{n} (z) = d a_{m} (z) a_{m}^{*} (z) a_{n} (z) a_{n}^{*} (z) + a_{m} (z) d a_{m}^{*} (z) a_{n} (z) a_{n}^{*} (z) \\ + a_{m} (z) a_{m}^{*} (z) d a_{n} (z) a_{n}^{*} (z) + a_{m} (z) a_{m}^{*} (z) a_{n} (z) d a_{n}^{*} (z) \\ + d a_{m} (z) d a_{m}^{*} (z) a_{n} (z) a_{n}^{*} (z) + a_{m} (z) a_{m}^{*} (z) d a_{n} (z) d a_{n}^{*} (z) \\ + d a_{m} (z) a_{m}^{*} (z) d a_{n} (z) a_{n}^{*} (z) + a_{m} (z) d a_{m}^{*} (z) a_{n} (z) d a_{n}^{*} (z) \\ + d a_{m} (z) a_{m}^{*} (z) a_{n} (z) d a_{n}^{*} (z) + a_{m} (z) d a_{m}^{*} (z) d a_{n} (z) a_{n}^{*} (z) \end{array}

The ODE for average value of P_m P_n (m≠n) is

\begin{array}{l} \frac{d 〈 P_{m} P_{n} 〉}{d z} = - 〈 P_{m} P_{n} 〉 \sum_{m'} κ_{m m'} F T (η) (β_{m} - β_{m'}) - 〈 P_{m} P_{n} 〉 \sum_{n'} κ_{n n'} F T (η) (β_{n} - β_{n'}) \\ + \sum_{m'} 〈 P_{m} (z) P_{m'} (z) 〉 κ_{m m'} F T (η) (β_{m} - β_{m'}) + \sum_{n'} 〈 P_{n} (z) P_{n'} (z) 〉 κ_{n n'} F T (η) (β_{n} - β_{n'}) \\ - 2 κ_{m n} F T (η) (0) 〈 P_{m} P_{n} 〉 \end{array}

The power standard deviations and the power correlation coefficients can be calculated as

\begin{array}{l} \sqrt{δ P_{m}^{2}} = \sqrt{〈 P_{m}^{2} 〉 - {〈 P_{m} 〉}^{2}} \\ C (P_{m}, P_{n}) = 〈 P_{m} P_{n} 〉 - 〈 P_{m} 〉 〈 P_{n} 〉 \end{array}

2.6 Analytical solutions for the ODEs of power and power variations

Equation (22), Eq. (25) and Eq. (27) are the ODEs with constant coefficients. They can be solved analytically. As indicated in [4–6 ], Eq. (22) can be reformulated as

\frac{d P}{d z} = K P

where vector P is composed by the average power of each mode, and the coupling matrix K is composed by the coupling coefficients indicated in Eq. (22). Equation (29) can be solved analytically as

P (z) = \exp (K z) P (0)

Equation (25) and Eq. (27) can be reformulated as

\frac{d Q (z)}{d z} = M Q (z)

where Q is a vector composed by the higher order statistics of the amplitudes,

Q = (\begin{matrix} 〈 P_{1} P_{1} 〉 \\ ⋮ \\ 〈 P_{1} P_{n} 〉 \\ \begin{matrix} 〈 P_{2} P_{1} 〉 \\ ⋮ \\ 〈 P_{2} P_{n} 〉 \\ \begin{matrix} ⋮ \\ 〈 P_{n} P_{1} 〉 \\ \begin{matrix} ⋮ \\ 〈 P_{n} P_{n} 〉 \end{matrix} \end{matrix} \end{matrix} \end{matrix})

Since M is a matrix with constant elements, Eq. (31) can be solved analytically as

Q (z) = \exp (M z) Q (0)

It should be noted that the above derivations have taken the leaky modes into account. When only the guided modes are considered, the impact of the leaky modes can be regarded as the radiation loss [4]. And this can be seen from Eq. (22), Eq. (25) and Eq. (27), by considering mode m' or mode n' as the radiation modes. Due to the fact that radiation modes dissipate rapidly, their power can be regarded as 0 in the coupling equation, and their contribution is only reflected in the radiation loss [4]. Not surprisingly, it can be interpreted from Eq. (22), Eq. (25) and Eq. (27) that the radiation loss can be replaced by a deterministic loss instead of a stochastic parameter. The radiation loss α_m of each mode can be evaluated according to the modal distribution and the waveguide parameters [4], and can be incorporated into Eq. (22), Eq. (25) and Eq. (27) easily. The details of loss evaluation will be presented in section 2.8.

2.7 Extension of formulas in the vector theory

Although the formulas have been derived based on the scalar Helmholtz equation, the results derived here can be easily extended for the vector case. A coupled equation can be derived based on the Maxwell's equations. It is in the form of Eq. (6), however, with the coupling coefficient modified as [19]

C_{m m'} (z) = \frac{- j ω ε_{0} \iint_{x y} e_{m}^{*} \cdot e_{m'} δ n^{2} d S}{\iint ((e_{m} \times h_{m}^{*}) + (e_{m}^{*} \times h_{m})) \cdot \vec{z} d S}

where e _m and h _m denote the electric field and the magnetizing field of the mth modes. if we define [19]

P_{m} = \frac{1}{4} \iint ((e_{m} \times h_{m}^{*}) + (e_{m}^{*} \times h_{m})) \cdot \vec{z} d S = 1

we have

C_{m m'} (z) = \frac{k_{0}^{2}}{j 4 ω μ_{0}} \iint_{x y} e_{m}^{*} \cdot e_{m'} δ n^{2} d S

which is almost in the same form as the scalar case. During the derivation of Eq. (34), the orthogonality between the modes is used [20]. The readers can find the derivations in [19] and in Marcuse' paper [21]. Obviously, the conclusion derived from the scalar coupled equation remains valid since the form of the coupled equation has not changed.

2.8 Evaluation of the modal losses in slab waveguides and planar waveguides

As suggested in the previous sections, the radiation modes will contribute to the loss of each guided mode. The loss can be evaluated by the coupling between the guided modes and the radiation modes.

For waveguides with high index contrast, e. g. lithographically etched ridge waveguides, the guided modes as well as the radiation modes should be studied based on the vector theory. While Eq. (10) stands, Eq. (11) should be modified as

κ_{m m'} = \frac{ω^{2} ε_{0}^{2}}{16} \iint_{x y} \iint_{x' y'} e_{m}^{*} (x, y) \cdot e_{m'} (x, y) e_{m} (x', y') \cdot e_{m'}^{*} (x', y') R (x - x', y - y') d S d S'

By replacing e_m' with the plane waves approximate the radiation modes and integrating over k _x and k _y, one may evaluate the mode loss.

In lithographically etched ridge waveguides, even in waveguides with high index contrast, the modes can be classified as quasi TE and quasi TM modes. For the quasi TE/TM modes, the z component of E/H is very small and can be neglected. Without loss of generality, we consider the quasi TE modes here, just like D. Lenz did in [9]. In such modes, the principal field components will be E_x and H_y [9], and E_x>>E_y. Therefore, Eq. (37) can be simplified as by only considering the x components of e_m and e_m' [9]. The waveguides usually have four different side walls, i.e. the top, the bottom, the left and the right side walls. For simplicity, we present the mode loss caused by the wall roughness on the top. The loss caused by other three walls can be evaluated similarly. The overall loss will be the summation of the losses caused by the roughness on the four walls, because the roughness on the walls can be regarded as independent random processes. We assume the top wall is located at the position of the y = b. In this case, Eq. (37) should be modified as

κ_{m m'} = \frac{ω^{2} ε_{0}^{2}}{16} \int_{x} \int_{x'} e_{m}^{*} (x, b) e_{m'} (x, b) e_{m} (x', b) e_{m'}^{*} (x', b) R_{x} (x - x') d x d x'

Here we directly use e_m and e_m' to denote the x components of the electrical fields so that the symbols do not mix with each other. As suggested above, we can use the plane waves to approximate the radiation modes [9], and therefore, e_m' should be replaced by

\begin{array}{l} e_{p} (k_{x}, k_{y}) = ρ_{p} (k_{x}, k_{y}) \exp (j k_{x} x + j k_{y} y) \\ ρ_{p} (k_{x}, k_{y}) = \frac{1}{2 π} \sqrt{\frac{2 ω μ_{0}}{β_{p} (k_{x}, k_{y})}} \\ β_{p} (k_{x}, k_{y}) = \sqrt{n_{2}^{2} k_{0}^{2} - k_{x}^{2} - k_{y}^{2}} \end{array}

where k _x and k _y are the x and y direction wave vectors of the plane wave, n₂ the cladding refractive index, β_p the propagation constant of the plane wave in the z direction. The plane wave in Eq. (39) has been normalized so that [9]

\iint_{x y} e_{p} (k_{x}, k_{y}) e_{p} (k_{x}', k_{y}') d x d y = \frac{2 ω μ_{0}}{β_{p} (k_{x}, k_{y})} δ (k_{x} - k_{x}') δ (k_{y} - k_{y}')

Substituting Eq. (39) into Eq. (38), one has the loss coefficient of the mth mode as

\begin{array}{l} α_{m} = \int_{k_{x}} \int_{k_{y}} κ_{m p} F T (η) (β_{m} - β_{p}) d k_{x} d k_{y} \\ κ_{m p} = \frac{ω^{2} ε_{0}^{2} {| ρ_{p} |}^{2}}{16} \int_{x} \int_{x'} e_{m}^{*} (x, b) e_{m} (x', b) \exp (j k_{x} (x - x')) R_{x} (x - x') d x d x' \end{array}

In the case of slab waveguide, we have e_m remain unchanged with respect to x and R(x-x') will be a constant with the value of (n₁ ²- n₂ ²) ², where n₁ is the core index. Therefore the integration over (x-x') will result in 2πδ(k _x). In this case, Eq. (42) is reduced to

α_{m} = \frac{ω^{2} ε_{0}^{2} ω μ_{0}}{16 π} {(n_{1}^{2} - n_{2}^{2})}^{2} e_{m}^{*} (b) e_{m} (b) \int_{k_{y}} F T (η) (β_{m} - β_{p}) \frac{1}{\sqrt{n_{2}^{2} k_{0}^{2} - k_{y}^{2}}} d k_{y}

Here we use e_m(b) to stand for the renormalized e_m(x,b) (It needs to be renormalized because the two dimensional field e_m(x,y) is reduced to the one dimensional field e_m(y)). By changing the variable,

k_{y} = n_{2} k_{0} \sin θ

we have the loss as

α_{m} = \frac{ω^{2} ε_{0}^{2} ω μ_{0}}{16 π} {(n_{1}^{2} - n_{2}^{2})}^{2} e_{m}^{*} (b) e_{m} (b) \int_{0}^{π} F T (η) (β_{m} - n_{2} k_{0} \cos θ) d θ

which is in exact agreement with the results derived by [8] if the normalization factor is properly chosen. For the waveguides with 3D structure, e. g. the rectangular waveguides, random index perturbation in the z direction is usually assumed in the published literatures [9]. Under such an assumption, for the top wall of the rectangular waveguide, i. e. the wall with the coordinates of y = b, the random index changes remain identical in the x direction [9], and therefore we have R(x-x') = (n₁ ²- n₂ ²) ². Substituting this into Eq. (41), we have the mode loss as

α_{m} = \int_{k_{x}} \int_{k_{y}} \frac{ω^{2} ε_{0}^{2} {| ρ_{p} |}^{2} {(n_{1}^{2} - n_{2}^{2})}^{2}}{16} {| \int_{x} e_{m}^{*} (x, b) \exp (j k_{x} x) d x |}^{2} F T (η) (β_{m} - β_{p}) d k_{x} d k_{y}

which is exactly the same as the mode loss derived in [9]. For a more general case, i. e. the index perturbation changes randomly both in the x and z directions on the top wall, Eq. (41) should be used.

In practical waveguides fabricated by the lithography process, the random variation of the side wall should be a two dimensional random process rather than the one dimensional random process in the z direction. The correlation length in the x/y direction should be in the same order as the correlation length in the z direction. Therefore, for the waveguides with short correlation length in the x/y direction, the radiation loss might be overestimated by the model in [9].

For example, the rectangular waveguide discussed in [9] is with the core index of 1.01 and cladding index of 1. The width and the height of the waveguide are both 30μm [9]. The autocorrelation function in the z direction is Gaussian function with the correlation length as 0.75μm [9]. The standard deviation of the roughness $\bar{σ}$ is 4 × 10⁻⁹ [9]. The signal wavelength is 1550nm. Using Eq. (45), the radiation loss for the $E_{x}^{11}$ mode can be evaluated as 0.0006 Np/m, or 0.005 dB/m, which is in agreement with the results in [9]. However, if we assume that the side wall roughness is a two dimensional process and the autocorrelation function in the x/y direction is Gaussian function, the radiation loss will be different. Figure 1 depicts the variation of the radiation loss with respect to the change of the correlation length in the x/y direction while other parameters are fixed.

Fig. 1 The radiation loss of $E_{x}^{11}$ mode versus correlation length in the x/y direction, the z directional correlation length is 0.75μm.

Download Full Size | PDF

It can be seen from the Fig. 1 that the radiation loss increases as the correlation length increases. Comparing the results when the x/y directional correlation length is infinitely long with the results when the x/y directional correlation length is the same as the z directional correlation length (i. e. 0.75μm), 0.0007 dB/m (15%) discrepancy is observed. However, when the z directional correlation length is shorter, e. g. 0.05 μm [10], the discrepancy between the two results can be as much as 80% (shown in Fig. 2 , comparing the results when the correlation length is 0.05 μm with the results when the correlation length is infinity). Therefore, the results in [9] are reliable only when the waveguides have large correlation lengths of roughness. For the waveguides fabricated by the advanced lithography process, which has the correlation length of roughness as short as 50 nm, Eq. (41) should be used to get a better estimation of the radiation loss.

Fig. 2 The radiation loss of $E_{x}^{11}$ mode versus correlation length in the x/y direction, the z directional correlation length is 50 nm.

Download Full Size | PDF

3. Simulation results and discussion

Numerical simulations are conducted to verify the derived analytical formulas. Without loss generality, a slab waveguide with random surfaces is simulated [4]. The slab waveguide is identical to the one in Marcuse' paper [4], which has the core index as 1.5 and the core and cladding index ratio as 1.01. The width of the waveguide is 2d and k₀d = 82 as indicated in [4], with the signal operation wavelength as 1550nm. The optical waveguide is able to hold ten modes [4]. The effective index of each mode is calculated and listed in Table 1 .

Table 1. Effective Indexes of the Ten Modes of the Waveguide

View Table | View all tables in this article

The optical waveguide has two random surfaces located at x = d and x = -d. Two independent processes f(z) and h(z) with the same correlation lengths are used to model the random index variations on the two surfaces [4]. The random index variations are usually caused by the fabrication process, but the random change of the temperature or the strain can also be a part of the reason.

According to Eq. (9) and Eq. (13), the coupling coefficient C_mm' can be modeled as [4]

\begin{array}{l} C_{m m'} (z) = \frac{(n_{1}^{2} - n_{2}^{2}) k_{0}^{2} ϕ_{m} ϕ_{m'} \sqrt{γ_{m} γ_{m'}}}{2 j \sqrt{β_{m} β_{m'} (1 + γ_{m} d) (1 + γ_{m'} d)}} (f (z) - {(- 1)}^{m + m'} h (z)) \\ ϕ_{m} = {\begin{matrix} \cos (κ_{m} d) m = 0, 2, 4 \dots \\ \sin (κ_{m} d) m = 1, 3, 5 \dots \end{matrix} \\ γ_{m} = \sqrt{β_{m}^{2} - n_{2}^{2} k_{0}^{2}} \\ κ_{m} = \sqrt{n_{1}^{2} k_{0}^{2} - β_{m}^{2}} \end{array}

where n ₁ and n ₂ are the core and cladding indexes, β _m the propagation constant of the m ^th mode, ϕ_m the modal distribution of the m ^th mode at the position of x = d.

The two random processes describing the index variations are assumed to have the following Gaussian correlation function [4]

〈 f (z) f (z') 〉 = 〈 h (z) h (z') 〉 = η (z - z') = {\bar{σ}}^{2} e^{- \frac{{(z - z')}^{2}}{D^{2}}}

where

{\bar{σ}}^{2}

is the strength of the roughness, D the correlation length.

The Fourier transform of the correlation function is [4]

F T ({\bar{σ}}^{2} e^{- \frac{{(z - z')}^{2}}{D^{2}}}) = {\bar{σ}}^{2} \sqrt{π} D e^{- \frac{D^{2} ω^{2}}{4}}

The correlation function of the coupling coefficient is [4]

〈 C_{m m'} (z) C_{m m'}^{*} (z') 〉 = \frac{{(n_{1}^{2} - n_{2}^{2})}^{2} k_{0}^{4} ϕ_{m}^{2} ϕ_{m'}^{2} γ_{m} γ_{m'}}{2 β_{m} β_{m'} (1 + γ_{m} d) (1 + γ_{m'} d)} {\bar{σ}}^{2} e^{- \frac{{(z - z')}^{2}}{D^{2}}}

As explained in the previous section, the radiation loss is the result of the random index variation, but it is equivalent to a deterministic loss in the SDEs. For each guided mode, its radiation loss can be evaluated as by section 2.8. Here, we directly use the expressions derived by Marcuse [4]

\begin{array}{l} α_{m} = \int_{- n_{2} k_{0}}^{n_{2} k_{0}} F T (η) (β_{m} - β) I_{m} (β) d β \\ I_{m} (β) = \frac{(n_{1}^{2} - n_{2}^{2}) k_{0}^{3} n_{1} \sin^{2} θ_{m}}{2 π d \cos θ_{m} (1 + γ_{m} d)} (\frac{ρ \cos^{2} σ d}{ρ^{2} \cos^{2} σ d + σ^{2} \sin^{2} σ d} + \frac{ρ \sin^{2} σ d}{ρ^{2} \sin^{2} σ d + σ^{2} \cos^{2} σ d}) \\ ρ = \sqrt{n_{2}^{2} k_{0}^{2} - β^{2}} \\ σ = \sqrt{n_{1}^{2} k_{0}^{2} - β^{2}} \\ \sin θ_{m} = \frac{(m + 1) π}{2 n_{1} k_{0} d} \end{array}

When the correlation length D >> d or D<<d, the analytical approximation of the radiation loss can be obtained. When D<<d, we have [4]

α_{m} = \frac{\sqrt{π} (n_{1}^{2} - n_{2}^{2}) k_{0}^{3} n_{1} \sin^{2} θ_{m}}{2 \cos θ_{m} (1 + \frac{1}{γ_{m} d})} {\bar{σ}}^{2} \frac{D}{d}

When D>>d, we have

\begin{array}{l} α_{m} = {\bar{σ}}^{2} e^{- \frac{D^{2}}{4} {(β_{m} - n_{2} k_{0})}^{2}} \frac{n_{1} k_{0}^{3} \sin^{2} θ_{m}}{2 d \sqrt{n_{2} k_{0} (β_{m} - n_{2} k_{0})} (1 + \frac{1}{γ_{m} d}) \cos θ_{m}} \\ \cdot {\begin{matrix} \begin{array}{l} \frac{2 n_{2}}{k_{0} D^{2} (β_{m} - n_{2} k_{0})} (\cot {}^{2}σ d + \tan^{2} σ d) a n d \begin{matrix} \tan σ d \neq 0 \\ \cot σ d \neq 0 \end{matrix} \end{array} \\ (n_{1}^{2} - n_{2}^{2}) o r \begin{matrix} \tan σ d = 0 \\ \cot σ d = 0 \end{matrix} \end{matrix} \end{array}

The analytical model derived in section 2 is used to calculate the power, the power variation and the power correlation coefficients. Meanwhile, Monte-Carlo simulation is performed based on the SDE in order to provide the comparisons for the analytical model.

Firstly, the correlation length D is assumed to be 35d [4]. Without loss generality, only the fundamental mode is injected into the waveguide with the power of 1. Normalized propagation distance and normalized power are used in the simulation. The distance is normalized as [4]

z \frac{{\bar{σ}}^{2} k_{0}^{2}}{d^{2}}

while the power at each mode is normalized with respect to the input power of the fundamental mode. In order to maintain the accuracy during the numerical integration, the propagation step cannot be too large. Therefore, it is quite time-consuming to conduct Monte Carlo simulation. In the simulation, the normalized propagation distance is limited to 16000, and totally 4000 rounds of simulations on the SDE are conducted to obtain the average values.

In Fig. 3 and Fig. 4 , the evolution of the powers of mode 1 to mode 4 with respect to the normalize propagation distance is illustrated. Due to the page limit, not all of the modes are illustrated in the Figs.. The output powers of all of the modes at the end the propagation are shown in Table 2 . The Monte Carlo simulation is accomplished by integrating the SDE (Eq. (14), while the analytical results are from the analytical solution of the ODEs (Eq. (29)). The curves from the Monte Carlo simulation agree well with the analytical predictions. It can be seen from the Figs. that the coupling between the fundamental mode and mode 2 is the strongest. With the increase of the mode order, the coupling decreases. This is because the coupling coefficients depend on the difference between the propagation constants. It gets smaller as the difference becomes larger. Furthermore, it can be seen that with the increase of the propagation distance, the average power of each mode tends to have an equal value, and this is in agreement with the Marcuse' theory [4]. The accuracy of the analytical coupled power model (which is also derived by Marcuse) with respect to the Mont Carlo simulation is evaluated. The relative error between the results is below 5%, and this demonstrates the validity of the proposed SDEs (i. e. Equation (14)) in this work.

Fig. 3 The power evolution of the fundamental mode (mode 1).

Download Full Size | PDF

Fig. 4 The power evolution of mode 2 to mode 4.

Download Full Size | PDF

Table 2. Output Power and Output Power Standard Deviations at the End of the Propagation

View Table | View all tables in this article

Afterwards, the higher order statistics of the amplitudes is calculated. The power standard derivations of mode 1 to mode 4 versus the normalized propagation distance are shown in Figs. 5-6 . The power standard deviations of each mode at the end of the propagation are shown in Table 2. It can be concluded that the analytical solution of the ODEs shown in Eq. (33) gives a very accurate prediction for the higher order statistics. The relative error is below 5%.

Fig. 5 The power standard deviation evolution of the fundamental mode.

Download Full Size | PDF

Fig. 6 The power standard deviation evolution of mode 2 to mode 4.

Download Full Size | PDF

Since the optical signal is deterministic at the beginning of the propagation, i. e. with the constant input optical power, the standard deviation of the mode power will be 0. Therefore, it is not possible to evaluate the normalized correlation coefficient (i.e. Equation (28) normalized by the standard deviations). The un-normalized power correlation coefficients between the powers of the fundamental mode and mode 2 to mode 4 are calculated according to Eq. (28) and are plotted in Fig. 7 . It can be seen from the Fig. that as the propagation distance increases, the correlation between the powers of the modes increases (regarding the absolution value of the correlation coefficient).

Fig. 7 The power correlation coefficients between mode 1 and mode 2 to mode 4 along the propagation distance.

Download Full Size | PDF

The information of the higher order statistics of the modal amplitudes, which is not studied before, is very useful for the system performance estimation as discussed in the previous sections. For example, from Fig. 5-6 , it can be seen that the power variation increases as the propagation distance increases, which brings sensitivity penalty during the receiver design. Also, as shown in Fig. 7, the correlation (indicated by the absolute value of the correlation coefficients) increases with respect to the increase of the propagation distance. If this correlation coefficient variation is caused by the random index variations induced by the temperature or strain in a single waveguide, it will result in the decrease of the system capacity.

The robustness of the analytical results is examined by changing the correlation length D, which varies from 0.01d to 100d [4]. The discrepancy between the analytical solution of the ODEs and the Monte Carlo simulation remains below 5%. Hence, the model proposed in this paper is quite robust to simulate random optical waveguides with different correlation lengths.

One thing worth mentioning is that the computational efficiency of the analytical model is greatly improved in comparison with Monte-Carlo simulation. The computational time has been reduced by a factor of 10000.

4. Conclusion

Based on the coupled mode theory, a set of SDEs have been proposed to analyze the random index variations inside the optical waveguides. ODEs are derived to characterize the power evolution as well as the higher order statistics of the modal amplitudes. Analytical solution to these ODEs can be provided because they are ODEs with constant coefficients. These ODEs not only save the computational efforts by avoiding the time-consuming Monte-Carlo simulations, but also bring a deeper understanding of the impact of the random index variations. Detailed Monte-Carlo simulations are provided to verify the derived ODEs. Although the simulation has been conducted on a slab waveguide suggested by Marcuse [4], the ODEs are expected to be functional for optical waveguides with arbitrary structures.

Appendix

In this appendix, we provide the detailed derivation of Eq. (14).

As discussed before, Eq. (6) should be interpreted in the Stratonovich sense, which indicates that the RHS of the Eq. is evaluated in the middle of the integration step, i.e.

d a_{m} (z) = \sum_{m'} a_{m'} (z + \frac{d z}{2}) \int_{z}^{z + d z} C_{m m'} e^{j (β_{m} - β_{m'}) z'} d z'

It is worth noting that C_mm' is actually a filtered Gaussian random process [22–24 ], i.e. the Gaussian process filtered by a linear filter with the impulse response of

F T^{- 1} (\sqrt{F T (〈 C_{m m'} (z) C_{m m'}^{*} (z') 〉)})

where FT denotes the Fourier transform with respect to the variable (z-z') and FT⁻¹ denotes the inverse Fourier transform. This linear filter can also be described by a linear ordinary differential equation (ODE) [22–24 ]. For example, when

〈 C_{m m'} (z) C_{m m'}^{*} (z') 〉 = \frac{1}{2 α} \exp (- α | z - z' |)

the corresponding linear ODE will be

\begin{array}{l} \frac{d C_{m m'}}{d z} + α C_{m m'} = n (z) \\ 〈 n (z) n (z') 〉 = δ (z - z') \end{array}

According to the SDE theory, the Stratonovich sense SDE should be converted to Ito sense, i.e. the RHS of the Eq. should be evaluated at the beginning of the integration step. Using the Taylor expansion, we have [25]

\begin{array}{l} a_{m'} (z + \frac{d z}{2}) = a_{m'} (z) + \frac{d a_{m'} (z)}{d z} \frac{d z}{2} + \dots \\ = a_{m'} (z) + \frac{d a_{m'} (z)}{2} + \dots \end{array}

Substituting Eq. (47) into Eq. (43), we have

d a_{m} (z) = \sum_{m'} a_{m'} (z) \int_{z}^{z + d z} C_{m m'} e^{j (β_{m} - β_{m'}) z'} d z' + \sum_{m'} \frac{d a_{m'} (z)}{2} \int_{z}^{z + d z} C_{m m'} e^{j (β_{m} - β_{m'}) z'} d z'

which is Eq. (14). It can be further simplified by submitting the expression of da_m' into Eq. (48), and we have

\begin{array}{l} d a_{m} (z) = \sum_{m'} a_{m'} (z) \int_{z}^{z + d z} C_{m m'} e^{j (β_{m} - β_{m'}) z'} d z' \\ + \sum_{m'} \frac{\sum_{n} a_{n} (z) \int_{z}^{z + d z} C_{m' n} e^{j (β_{m'} - β_{n}) z'} d z'}{2} \int_{z}^{z + d z} C_{m m'} e^{j (β_{m} - β_{m'}) z'} d z' + O (d z^{3 / 2}) \end{array}

It is shown in Marcuse's book [12] that the second term on the RHS of Eq. (49) can be simplified by letting n=m, because only under this condition, will the fast oscillation term in the z direction vanish, and will the value of the expression become non zero after averaging over z.

Hence, we have

\begin{array}{l} d a_{m} (z) = \sum_{m'} a_{m'} (z) \int_{z}^{z + d z} C_{m m'} e^{j (β_{m} - β_{m'}) z'} d z' \\ + \frac{1}{2} (\sum_{m'} \int_{z}^{z + d z} C_{m' m} e^{j (β_{m'} - β_{m}) z'} d z' \int_{z}^{z + d z} C_{m m'} e^{j (β_{m} - β_{m'}) z'} d z') a_{m} (z) + O (d z^{3 / 2}) \end{array}

Since C_m'm=C_mm'*, we have

\begin{array}{l} d a_{m} (z) = \sum_{m'} a_{m'} (z) \int_{z}^{z + d z} C_{m m'} e^{j (β_{m} - β_{m'}) z'} d z' \\ + \frac{1}{2} (\sum_{m'} \int_{z}^{z + d z} C_{m m'} e^{j (β_{m} - β_{m'}) z'} d z' \int_{z}^{z + d z} C_{m m'}^{*} e^{- j (β_{m} - β_{m'}) z'} d z') a_{m} (z) + O (d z^{3 / 2}) \end{array}

It should be noted that in most of the text books, the SDE is usually driven by the Brownian motion, which indicates that the correlation function of the driving force is the delta function. However, the filtered Brownian motion, such as the Ornstein–Uhlenbeck process, can also be regarded as the driving force of the SDE and the results have been verified by the Mont Carlo simulations [22,23 ].

Acknowledgment

This work is partially supported by the National science foundation of China (NSFC) (Grant No. 61201068) the Fundamental Research Funds for the Central Universities of China. The authors would like to thank the two anonymous reviewers for their helpful suggestions during the review process.

References and links

1. A. Di Donato, M. Farina, D. Mencarelli, A. Lucesoli, S. Fabiani, T. Rozzi, G. M. Di Gregorio, and G. Angeloni, “Stationary Mode Distribution and Sidewall Roughness Effects in Overmoded Optical Waveguides,” J. Lightwave Technol. 28(10), 1510–1520 (2010). [CrossRef]

2. T. Barwicz and H. Haus, “Three-dimensional analysis of scattering losses due to sidewall roughness in micro-photonic waveguides,” J. Lightwave Technol. 23(9), 2719–2732 (2005). [CrossRef]

3. A. W. Snyder and J. Love, Optical Waveguide Theory, 1st ed., Springer, 1983.

4. D. Marcuse, “Mode conversion caused by surface imperfections of a dielectric slab waveguide,” Bell Syst. Tech. J. 48(10), 3187–3215 (1969). [CrossRef]

5. D. Marcuse, “Derivation of coupled power equations,” Bell Syst. Tech. J. 51(1), 229–237 (1972). [CrossRef]

6. D. Marcuse, “Power distribution and radiation losses in multimode dielectric waveguides,” Bell Syst. Tech. J. 51(2), 429–454 (1972). [CrossRef]

7. J. D. Love, T. J. Senden, and F. Ladouceur, “Effect of side wall roughness in buried channel waveguides,” IEE Proc., Optoelectron. 141(4), 242–248 (1994). [CrossRef]

8. J. Lacey and F. Payne, “Radiation loss from planar waveguides with random wall imperfections,” IEE Proc. 137, 282–288, 1990. [CrossRef]

9. D. Lenz, D. Erni, and W. Bächtold, “Modal power loss coefficients for highly overmoded rectangular dielectric waveguides based on free space modes,” Opt. Express 12(6), 1150–1156 (2004). [CrossRef] [PubMed]

10. F. Grillot, L. Vivien, S. Laval, D. Pascal, and E. Cassan, “Size Influence on the Propagation Loss Induced by Sidewall Roughness in Ultrasmall SOI Waveguides,” IEEE Photonics Technol. Lett. 16(7), 1661–1663 (2004). [CrossRef]

11. A. D. Simard, N. Ayotte, Y. Painchaud, S. Bédard, and S. LaRochelle, “Impact of Sidewall Roughness on Integrated Bragg Gratings,” J. Lightwave Technol. 29(24), 3693–3704 (2011). [CrossRef]

12. D. Marcuse, Theory of Dielectric Optical Waveguides (Academic Press, 1974).

13. K.-P. Ho and J. M. Kahn, “Mode coupling and its impact on spatially multiplexed systems,” in Optical Fiber Telecommunications (OFC) VI, B, I. P. Kaminow, T. Li, and A. E. Willner, eds. (2013), pp. 491–568.

14. K.-P. Ho and J. M. Kahn, “Frequency diversity in mode-division multiplexing systems,” J. Lightwave Technol. 29(24), 3719–3726 (2011). [CrossRef]

15. S. O. Arik, D. Askarov, and J. M. Kahn, “Effect of mode coupling on signal processing complexity in mode-division multiplexing,” J. Lightwave Technol. 31(3), 423–431 (2013). [CrossRef]

16. P. J. Winzer and G. J. Foschini, “MIMO capacities and outage probabilities in spatially multiplexed optical transport systems,” Opt. Express 19(17), 16680–16696 (2011). [CrossRef] [PubMed]

17. M. Huang and X. Yan, “Thermal-stress effects on the temperature sensitivity of optical waveguides,” J. Opt. Soc. Am. B 20(6), 1326–1333 (2003). [CrossRef]

18. J. A. Anguita, M. A. Neifeld, and B. V. Vasic, “Spatial correlation and irradiance statistics in a multiple-beam terrestrial free-space optical communication link,” Appl. Opt. 46(26), 6561–6571 (2007). [CrossRef] [PubMed]

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

20. J. Seligson, “The orthogonality relation for TE- and TM-modes in guided-wave optics,” J. Lightwave Technol. 6(8), 1260–1264 (1988). [CrossRef]

21. D. Marcuse, “Coupled-mode theory for anisotropic optical waveguides,” Bell Syst. Tech. J. 54(6), 985–995 (1975). [CrossRef]

22. P. K. A. Wai and C. R. Menyuk, “Polarization mode dispersion, decorrelation and diffusion in optical fibers with randomly varying birefringence,” J. Lightwave Technol. 14(2), 148–157 (1996). [CrossRef]

23. A. Galtarossa, L. Palmieri, A. Pizzinat, B. S. Marks, and C. R. Menyuk, “An analytical formula for the mean differential group delay of randomly-birefringent spun fibers,” J. Lightwave Technol. 21(7), 1635–1643 (2003). [CrossRef]

24. A. Mecozzi, “A Theory of Polarization-Mode Dispersion of Spun Fibers,” J. Lightwave Technol. 27(7), 938–943 (2009). [CrossRef]

25. Q. Lin and G. P. Agrawal, “Vector theory of stimulated Raman scattering and its application to fiber-based Raman amplifiers,” J. Opt. Soc. Am. B 20(8), 1616–1631 (2003). [CrossRef]

Mode order	effective index
1	1.499891
2	1.499563
3	1.499017
4	1.498254
5	1.497275
6	1.496082
7	1.494678
8	1.493068
9	1.49126
10	1.489269

	power after propagation of 16000 normalized units		power standard deviation after propagation
Mode order	Analytical	Monte-Carlo	Analytical	Monte-Carlo
1	0.8943	0.8907	0.0689	0.0720
2	0.0588	0.0613	0.0565	0.0593
3	0.0323	0.0327	0.0326	0.0339
4	0.0106	0.0110	0.0128	0.0126
5	0.0032	0.0034	0.0049	0.0052
6	0.0007	0.0008	0.0013	0.0013
7	0.0001	0.0001	0.0002	0.0002
8	0.0000	0.0000	0.0000	0.0000
9	0.0000	0.0000	0.0000	0.0000
10	0.0000	0.0000	0.0000	0.0000

Mode order	effective index
1	1.499891
2	1.499563
3	1.499017
4	1.498254
5	1.497275
6	1.496082
7	1.494678
8	1.493068
9	1.49126
10	1.489269

	power after propagation of 16000 normalized units		power standard deviation after propagation
Mode order	Analytical	Monte-Carlo	Analytical	Monte-Carlo
1	0.8943	0.8907	0.0689	0.0720
2	0.0588	0.0613	0.0565	0.0593
3	0.0323	0.0327	0.0326	0.0339
4	0.0106	0.0110	0.0128	0.0126
5	0.0032	0.0034	0.0049	0.0052
6	0.0007	0.0008	0.0013	0.0013
7	0.0001	0.0001	0.0002	0.0002
8	0.0000	0.0000	0.0000	0.0000
9	0.0000	0.0000	0.0000	0.0000
10	0.0000	0.0000	0.0000	0.0000

Comprehensive analytical model to characterize randomness in optical waveguides

Abstract

1. Introduction

2. Mathematical model

2.1 Coupled mode theory inside optical waveguides

2.2 The Stochastic differential equation model

2.3 Energy conservation

2.4 Coupled power equation

2.5 Coupled equation for the higher order statistics

2.6 Analytical solutions for the ODEs of power and power variations

2.7 Extension of formulas in the vector theory

2.8 Evaluation of the modal losses in slab waveguides and planar waveguides

3. Simulation results and discussion

4. Conclusion

Appendix

Acknowledgment

References and links

Cited By

Figures (7)

Tables (2)

Equations (63)

Optics Express