# Modeling and Simulation of Piecewise Regular Multimode Fiber Links Operating in a Few-Mode Regime

Authors: Anton Bourdine
Journal: Advances in Optical Technologies
 (2013)

---

## Abstract

This work presents an alternative model of multimode fiber links with conventional silica weakly-guiding graded-index irregular multimode fibers under a few-mode optical signal propagation generated by laser source. The proposed model is based on the piecewise regular representation. It takes into account launch conditions, differential mode delay, both lower- and higher-order mode chromatic dispersion, differential mode attenuation, and mode mixing and power diffusion occurring due to real fiber irregularity and micro- and macrobends. We present some results of introduced model approbation with following pulse propagation simulations. A close matching with measured pulse responses at the output of test fibers is noticed.


---

## Body

## 1. Introduction

Silica multimode graded-index fibers are used in a wide variety of applications ranging from on-board to in-premises networks links with length not more than 1–2 km. Since IEEE 802.3z standard was ratified and commercial available SFP transceivers with Vertical Cavity Surface Emitting Lasers (VCSELs) appeared on telecommunication market, multimode fibers became very popular both for in-building structural cabling systems (SCS) and typical distributed networks with backbone/vertical fiber cabling systems, private networks in premises, and campus environments. Nowadays the most top applications of multimode fibers are associated with data center SCS, high bit-rate storage area networks, and radio-over-fiber (RoF) techniques over already installed multimode fiber infrastructure inside building [1–3].Modern commercial multi-Gigabit transmission transceivers are realized on Vertical Cavity Surface Emitting Lasers (VCSELs) or single-mode Fabry-Perot laser diodes (LD) [1–4]. Because emission from the conventional VCSEL usually consists of about 5 or 6 transversal modes LPlm with maximal azimuthal order l not more than 3, and single-mode LD injects just fundamental and lower-order modes LP01 and LP11, only several guided modes are excited in multimode fiber link [5–7]. Therefore optical signal propagates over multimode fiber link in a so-called few-mode regime, and passage to the simulation of a few-mode pulse transmission over multimode fiber requires taking into account both “individual” dispersion parameters of mode component with particular order (amplitude, attenuation, delay, chromatic dispersion, etc.) and mode coupling.A variety of methods have been developed for modeling and simulation of laser-based multi-Gigabit data transmission over multimode fibers. Monograph [4] can be considered as fundamental complete basis work. The following groups of publications—[8–15]—should also be noted: their authors are affiliated to IEEE 802.3z, 802.3ae, 802.3aq, and 802.3ba Task Forces. Other works are also concerned with the design and simulation of multi-Gigabit multimode optical communication systems. Thus model [16–18] was realized in commercial RSoft tool ModeSYS, [19, 20] are based on vector methods, and [21] applies ray-tracing analysis. All mentioned models [4–20] are primarily assigned for the estimation of differential mode delay (DMD), which is the main issue of pulse distortion for laser-based multi-Gigabit data transmission over multimode fiber links and is concerned with theoretical research of DMD dependence on source parameters, launch conditions, and real fiber refractive index profile defects and deviations from the optimum graded-index profile form. They also provide accurate simulation technique of signal generation, conversion and processing in transceivers of multi-Gigabit communication systems.A generalized chromatic dispersion parameter is applied in some models mentioned above. It is based on well-known empiric formula, which is cited in the most of fiber specifications and data sheets [22]. However this approach does not differ of lower- and higher-order modes and does not solve the problem of “individual” chromatic dispersion estimation for guided mode with particular order. Also mode coupling is neglected in most of the listed above. It is justified by the short length of multimode fiber link. In some works correction factors or modal noise penalty are applied, but only mode power redistribution at fiber interconnections in fiber optic distribution systems, consolidation points, patch-panels, and so forth, is recalculated without taking into account both real fiber irregularities and fiber bends.Here we introduce an alternative model of multimode fiber links with conventional silica weakly-guiding graded-index irregular multimode fibers under a few-mode optical signal propagation generated by laser source. Proposed model is based on the piecewise regular representation. It takes into account launch conditions, differential mode delay, both lower- and higher-order mode chromatic dispersion, differential mode attenuation, and mode mixing and power diffusion occurring due to real fiber irregularity and micro and macro-bends. We also present some results of introduced model approbation with following pulse propagation simulations and matching with measured pulse responses at the output of test fibers.


## 2. Brief Overview of Methods for Simulation of Pulse Propagation over Irregular Optical Fibers

A several methods have been developed for computing electromagnetic wave fields during propagation over irregular waveguides. Monograph [23] is one of the fundamental works of this direction. Effects of wave interaction in irregular waveguides are exact described by the coupled mode theory. Complete overview, description, classification, and list of applications of coupled mode theory methods can be found in monograph [24].Generally coupled mode theory is based on representation of mode fields along longitudinal cross-sectionz of irregular optical waveguide in the kind of mode field superposition of complete system of guided and leaky modes in regular optical waveguide. This expansion varies along the fiber longitudinal axis and it is described by system of coupled mode equations. The system of equations is formulated by substitution of mentioned irregular optical waveguide mode field representation to Maxwell equations [25] or Helmholtz equation [26, 27]. Result of substitution is infinite system of differential equations of the first or second order (it depends on initial equations—Maxwell or Helmholtz) for each mode amplitude depending on longitudinal coordinate z.The system of equations exactly describes mode fields in longitudinal cross-section of irregular optical waveguide, but it has not general analytical solution. Problem can be much simplified for the limited number of modes—just one or two [28–30]. For example coupled mode theory is effectively applied for analysis of the fundamental mode propagation over irregular optical fiber under taking into account the birefringence effect [25–27, 31]. However few-mode regime of signal transmission over multimode fibers requires increase of mode components that much complicates solutions of the system of equations. This fact limits coupled mode theory application in the case of random fluctuations of fiber core geometry parameters along the longitudinal axis z.In [32, 33], mode components of few-mode signal are combined in mode groups with roughly equal propagation constants. In works [34–36], a passage to infinite number of mode components (mode continuum) is proposed followed by replacing the sum of discrete mode component fields by integration. By neglecting core geometry variation, random fiber macrobends is considered as the main factor of mode coupling in [37–39]. Also in [40–44] a special curvilinear coordinates are applied to describe mode coupling via fiber bends. In [45–47], mode mixing and power diffusion is simulated by combination of coupled mode theory and ray-tracing analysis.The system of incoherently coupled nonlinear Schrodinger equations is applied for mathematical description of pulse propagation over multimode fiber in [48–54]. Although the generalized system of equations was written for infinite number of modes by taking into account fiber nonlinearity, the main application of this method is used just for simulation of coupling between two orthogonally polarized modes in single mode fibers occurring due to birefringence.Piecewise regular representation [25, 26, 55] of irregular optical waveguide is an alternative method to coupled mode theory. It is basic for direct methods like finite-difference time-domain (FDTD) method—one of the most accurate, universal, and “popular” numerical methods, which is widely used for analysis of electromagnetic wave propagation over waveguides with arbitrarily complicated configuration (e.g., [56–62], etc.). Also well-known direct rigorous methods are method of moments [58, 62, 63], finite-difference frequency-domain (FDFD) method [55, 58, 62], finite element method [55, 58, 62], wavelets [64–66] and their modifications, and combinations with elements of coupled mode theory [41–44, 67–73].Rigorous methods provide accurately computing of the transversal mode field distribution at the particular distance of emission source and time interval under the lowest error. The main issue of their application is complicated formulation of boundary and initial conditions according to individual problem and waveguide properties and configuration. They are very effective for the modeling of waveguide structures with quite short dimensions—not more than(103,…,104)·λ—like integral optics components. However simulation of pulse propagation over conventional fiber link (with length of hundreds or even tens meters) by the mentioned direct rigorous methods requires extra computing resources and is almost impossible.In terms of piecewise regular representation of both LAN, MAN (metropolitan area network) and long-haul fiber links, matrix methods [62, 74–80], recurrent methods [81–85] and quasi-analytic time-domain models based on quasianalytic methods [86–89] and spatial transform methods [55, 86, 90–95] are widely used. The most of them were developed for simulation of pulse propagation over single mode fibers.From this point of view, method described in [96–98] can be considered as an exception. It combined coupled mode theory with piecewise regular representation and was developed for modeling of signal transmission processes over long irregular multimode fiber links. Here a random fiber bends are preset and supposed as the main reason of mode coupling occurring. Also it takes into account birefringence. However the method neglects core geometry variation along longitudinal fiber axis and considers only multimode fibers with ideal infinite α-profiles.Therefore the most considered methods above for modeling and analysis of pulse propagation over optical fibers should be adapted for irregular multimode fiber link operating in a few-mode regime. Model has to take into account the whole total selected excited guided mode composition of transferred few-mode optical signal and their “individual” parameters during propagation over multimode fibers with real graded-index profiles differing form idealα-profiles by local defects and refractive index fluctuation. Moreover, in this case division of mode composition on lower- and higher-order mode groups with averaged propagation constant is unacceptable. Also model has to take into account launch conditions and mode coupling due to both variation of core geometry parameter along fiber length and micro- and macrobends of fiber occurring in real cable links.


## 3. Model of Piecewise Regular Multimode Fiber Link

In this work we also apply piecewise regular representation combined with general approach of split-step method [99] to simulate a processes of mode mixing and power diffusion due to mode coupling (Figure 1). Here single silica weakly guiding circular multimode fibers with arbitrary axial symmetric refractive index profile with single continue outer claddings are considered. According to piecewise regular representation the fiber is divided into regular spans with length Δz. Inside the span fiber geometry parameters are considered as constant, and modes propagate independently without interaction and mixing. It is supposed that each excited guided mode with propagation constant, varying from one regular span to another span, satisfies to cut-off condition for whole regular spans composing fiber. Also it is assumed that at the link transmitter end each excited guided mode starts transferring single optical pulse in particular form identical to input signal (e.g., Gaussian). During pulse propagation over regular span its amplitude decreases due to mode attenuation. The signal is mainly distorted due to difference between group velocities and amplitudes of modes, that is, DMD effect. Also transferred by each mode pulse, spreading due to chromatic dispersion for particular mode is taken into account.Figure 1
Piecewise regular representation of 200 m irregular multimode optical fiber with varying core diameter.Boundaries of regular spans can be represented generally as ideal axially alignment splice of two almost similar optical fibers with mismatching parameters. However it is correct only for “straight” fibers. Therefore it is propose to simulate fiber bends by representation of boundaries in the form of splice of two mismatched fibers with some low angular misalignment. Mode power redistribution between the amplitudes of signal components as a result of mode interaction is estimated by mode coupling coefficient matrix computing at the joints of regular spans. Here only guided modes are considered as the main issue under pulse dynamics research during propagation over multimode fiber link in a few-mode regime. However, in fact, power loss due lo component transformation from guided to leaky mode, and reflections are also indirectly taken into account.At the receiver end the resulting pulse envelope is considered as superposition of all mode components of signal. Here it is propose to apply a well-known expression for the frequency response of the signal transferred byM mode components LPlm over regular multimode fiber with length z [4, 8–15, 18, 99–101]:
(1)HRxt(ω,z)=F[hTx(t)]∑pMAp(0)exp(-αpz)×exp[-j(ω-ω0)β1(p)z]×exp[-j12(ω-ω0)2β2(p)z],
where F is direct Fourier transform; hTx(t) is initial pulse at the transmitter end; Ap(0) and αp are starting amplitude and mode attenuation of pth guided mode LPlm(p=1,…,M); β1(p) and β2(p) are first- and second-order dispersion parameters. These dispersion parameters are elements of well-known Taylor series expansion approximation of propagation constant frequency dependence β(ω) [2, 55, 99]:
(2)β(ω)=(ω-ω0)β1+12(ω-ω0)2β2+16(ω-ω0)3β3,…,
where
(3)βn=[∂nβ∂ωn]ω=ω0.
Here β1(p)=τ(p) is mode delay and β2(p) is group velocity dispersion (GVD) associated with chromatic dispersion of pth guided mode LPlm.Therefore, according to introduced piecewise regular representation of irregular multimode fiber link, the frequency response of a few-mode optical signal transferred byM guided modes over irregular multimode fiber with length z under given particular length of regular span Δz by taking into account expression (1) can be written in the following form:
(4)HRx(ω,z)=F[hTx(t)]∑pMA(Nz+1)(p)exp[-α(Nz+1)(p)(z-NzΔz)]×exp[-j(ω-ω0)τ(p;Nz+1)(z-NzΔz)]×exp[-j12(ω-ω0)2β2(p;Nz+1)(z-NzΔz)]×∏q=1NzAq(p)exp(-αq(p)Δz)×exp(-j(ω-ω0)τ(p;q)Δz)×exp(-j12(ω-ω0)2β2(p;q)Δz),
where Nz=E(z/Δz); E(x) is the integer part of real number x.The resulting pulse response at the receiver end of irregular multimode link is computed by the following simple expression:(5)hRx(t)=F-1(HRx(ω))·[F-1(HRx(ω))]*,
where F-1 is inverse Fourier transform; [x]* is the complex conjugate of x.It is obvious that the proposed model requires fast and simple method for the analysis of the optical fibers with complicated refractive index profile form close to profiles of real fibers, which provides evaluation of parameters not only for the fundamental but also higher-order guided modes. For this purpose we propose to apply an extension of modified Gaussian approximation (EMGA) generalized for estimation of any order guided mode parameters, propagating over silica weakly guiding optical fibers with an arbitrary axial-symmetric refractive index profile [102]. EMGA is based on combination of Gaussian approximation and stratification method. It permits a much reduce computing time, especially for the calculations of higher-order mode parameters under the low error [103] by making a passage to analytical expressions both for the variational expression and characteristic equation. Below a detailed description of EMGA and its application for the proposed model is presented.


## 4. Extension of Modified Gaussian Approximation

EMGA is based on classical Gaussian approximation [25] of radial mode field distribution Fm(l)(R) in the weakly guiding optical waveguide with an arbitrary refractive index profile by the well-known Laguerre-Gauss function expression [4, 25, 100], describing a mode field distribution in weakly guiding optical waveguide with ideal infinite parabolic index profile:
(6)Ψm(l)(R)=(RR0)lLm-1(l)(R2R02)exp(-R22R02),
where Lm-1(l)is Laguerre polynomial.EMGA leads to equivalent normalized mode field radiusR0 estimation by solving a characteristic equation, which is derived from propagation constant variational expression under following passage to square core mode parameter U2 variational expression, written for analyzed weakly guiding optical waveguide with given refractive index profile. Parameter R0 is basic for this method and completely defines mode transmission parameters.Unlike known methods and their modifications [25, 104–110], based on conventional Gaussian approximation, in EMGA optical fiber with an arbitrary graded axial-symmetric index profile is considered as multicladding optical fiber. Therefore, refractive index profile inside core region can be represented in the form of the set of N layers in which the refractive stays a constant [102, 103]:
(7)n(R)={nk,Rk=kN,0≤k≤N-1,nN,1<R≤+∞,
and any profile function f(R) related to refractive index profile n(R) as
(8)n2(R)=nmax2[1-2Δ·f(R)],
where Δ=(nmax2-nN2)/2nmax2 is profile height parameter, can be written in terms of profile parameter hk:(9a)f(R)={hk,Rk=kN,0≤k≤N-1,1,1<R≤+∞,
where
(9b)hk=nmax2-nk2nmax2-nN2,where nk is refractive index of k layer (k=0,…,N); nmax is the maximal core refractive index; nN is cladding refractive index.This refractive index profile representation permits writing the variational expression for core mode parameterU2 and characteristic equation for normalized equivalent mode field radius ∂U2/∂R0=0 in the form of finite nested sums as follows [15, 16]:
(10)U2=(m-1)!(l+m-1)!×{QR02+V2[X0+∑k=0N-1hk(X1-X2)]},X0=exp(-1R02)∑q=02m-2Dq∑p=0l+q(l+q)!p!R02p,X1=exp(-k2N2R02)∑q=02m-2Dq∑p=0l+q(l+q)!p!R02pk2pN2p,X2=exp(-(k+1)2N2R02)∑q=02m-2Dq∑p=0l+q(l+q)!p!R02p(k+1)2pN2p,-Q+V2[S0+∑k=0N-1hk(S1-S2)]=0,(11)S0=exp(-1R02)∑q=02m-2Dq∑p=0l+q(l+q)!p!R02p(1-pR02),S1=exp(-k2N2R02)∑q=02m-2Dq×∑p=0l+q(l+q)!p!R02pk2pN2p(k2N2-pR02),S2=exp(-(k+1)2N2R02)∑q=02m-2Dq×∑p=0l+q(l+q)!p!R02p(k+1)2pN2p((k+1)2N2-pR02),Q=(l+m-1)!(3l+2m-1)(m-1)!+2l2∑q=02m-2Dq(q+l-1)!-4l×∑q=02m-2Cq(q+l)!,
where(12a)Dq=∑p=max(0,q-m+1)min(q,m-1)bp(l,m-1)bq-p(l,m-1),(12b)Lm-1(l)(x)Lm-1(l+1)(x)=∑q=02m-2Cqxq,and bp(l,m) is coefficient of polynomial representation in the form of power series [111, 112]: (13a)Lm(l)(x)=∑q=0mbq(l,m)xq,(13b)bq(l,m)=(-1)q(l+m)!(l+q)!(m-q)!q!,V=k0anmax2Δ is normalized frequency; k0=2π/λ is wavenumber; λ is wavelength.Therefore analysis of weakly guiding single-cladding optical fiber with an arbitrary profile leads to the following. Refractive index profile is represented by the profile function (9a) and (9b) in the form of N layers. Fiber parameters and mode orders l and m are substituted to the characteristic equation (11). By means of numerical solution (11), the normalized equivalent mode field radius R0 will be obtained. Then R0 is substituted to the expression (10), and mode core parameter U is estimated; that permits evaluating the propagation constant β for guided mode LPlm by the well-known expression [4, 25–27, 55, 100]:
(14)β2=k02nmax2-U2a2.Solution of the characteristic equation (11) is correct under normalized frequency V>1, and it should satisfy the guided mode cutoff condition [4, 25–27, 55, 100]:
(15)k0nN<β≤k0nmax.Optical confinement factorPcore can be considered as the second criterion for the identification of the ghost solutions:
(16)Pco(lm)≥0.5.
If we take into account Gaussian approximation, parameter Pcore is defined by analytical expression derived in [15] from the generalized integral form for weakly guiding optical fibers presented in [25]:
(17)Pco(lm)=(m-1)!(l+m-1)!∑q=02m-2Dq(l+q)!·Wq,
where
(18)Wq=1-exp(-1R02)∑p=0l+q1p!R02p.


## 5. Mode Group Velocity

According to what is mentioned above, realization of introduced model of piecewise regular multimode fiber link requires a passage to analytical expressions for mode group velocity. Guided mode group velocityvg is related to propagation constant β by the well-known formula [4, 25–27, 55, 100]:
(19)vg=-2πcλ2∂λ∂β.By squaring and following differentiation of the expression (14) with respect to wavelength, the first partial derivative of the square of propagation constant will be obtained:
(20)∂β2∂λ=-2k02nmax2λ+k02∂nmax2∂λ-1a2∂U2∂λ.After substitution of (20) to (19), an expression for the mode group velocity can be rewritten in the following form:
(21)vg=-4πβcλ2∂λ∂β2=-k02βcπ∂λ∂β2.Then by differentiating expressions (10) and (11) with respect to wavelength, first derivatives of square core mode parameter U2 and normalized equivalent mode field radius R0 will be obtained:
(22)∂U2∂λ=(m-1)!(l+m-1)!×{-2WR03∂R0∂λ+∂V2∂λ[X0+∑k=0N-1hk(X1-X2)]000+2V2X0(1)∂R0∂λ+V2∑k=0N-1∂hk∂λ(X1-X2)000+2V2∂R0∂λ∑k=0N-1hk(X1(1)-X2(1))},X0(1)=exp(-1R02)∑q=02m-2Dq∑p=0l+q(l+q)!p!R02p+3(1-pR02),X1(1)=exp(-k2N2R02)∑q=02m-2Dq×∑p=0l+q(l+q)!p!R02p+3k2pN2p(k2N2-pR02),X2(1)=exp(-(k+1)2N2R02)∑q=02m-2Dq×∑p=0l+q(l+q)!p!R02p+3(k+1)2pN2p((k+1)2N2-pR02),∂R0∂λ=(∑k=0N-1[(S0-S1)(hk∂V2∂λ+V2∂hk∂λ)]S0∂V2∂λ0+∑k=0N-1[(S0-S1)(hk∂V2∂λ+V2∂hk∂λ)])×(-2V2·[S0(1)+∑k=0N-1hk(S1(1)-S2(1))])-1,S0(1)=exp(-1R02)∑q=02m-2Dq∑p=0l+q(l+q)!p!R02p+3[(1-pR02)2-pR04],S1(1)=exp(-k2N2R02)∑q=02m-2Dq×∑p=0l+q(l+q)!k2pp!N2pR02p+3[(k2N2-pR02)2-pR04],S2(1)=exp(-(k+1)2N2R02)∑q=02m-2Dq×∑p=0l+q(l+q)!(k+1)2pp!N2pR02p+3×[((k+1)2N2-pR02)2-pR04].


## 6. Mode Chromatic Dispersion Parameter

The same technique is applied to derive the analytical expressions for mode chromatic dispersion parameterD associated with GVD. Chromatic dispersion parameter D is related to propagation constant β by the well-known formula [4, 25–27, 55, 100]:
(23)D=-λ2πc(2∂β∂λ+λ∂2β∂λ2).By differentiating (20), a second partial derivative of the square of propagation constant will be obtained:
(24)∂2β2∂λ2=k02(∂2nmax2∂λ2+6nmax2λ2-4λ∂nmax2∂λ)-1a2∂2U2∂λ2.Therefore expression forD leads to following form:
(25)D=-πk02βc[2λ∂β2∂λ+∂2β2∂λ2-12β2(∂β2∂λ)2].By differentiating expressions (22) with respect to wavelength, required second derivatives of square core mode parameter U2 and normalized equivalent mode field radius R0 will be obtained:
(26)∂2U2∂λ200=(m-1)!(l+m-1)!0000×{2MR03[3R0(∂R0∂λ)2-∂2R0∂λ2]0000000+∑k=0N-1(X1-X2)0000000×(hk∂2V2∂λ2+2∂hk∂λ∂V2∂λ+V2∂2hk∂λ2)0000000+2∑k=0N-1(X1(1)-X2(1))0000000×(2hk∂V2∂λ∂R0∂λ+2V2∂hk∂λ∂R0∂λ+V2hk∂2R0∂λ2)0000000+4V2(∂R0∂λ)2[X0(2)+∑k=0N-1hk(X1(2)-X2(2))]0000000+X0∂2V2∂λ2+2X0(1)(2∂V2∂λ∂R0∂λ+V2∂2R0∂λ2)},X0(2)=exp(-1R02)∑q=02m-2Dq×∑p=0l+q(l+q)!p!R02p+6[(1-pR02)2+R022(pR02-3)],X1(2)=exp(-k2N2R02)∑q=02m-2Dq×∑p=0l+q(l+q)!p!R02p+6k2pN2p×[(k2N2-pR02)2+R022(pR02-3k2N2)],X2(2)=exp(-(k+1)2N2R02)∑q=02m-2Dq×∑p=0l+q(l+q)!p!R02p+6(k+1)2pN2p×{[(k+1)2N2-pR02]200+R022[pR02-3(k+1)2N2]},∂X0(1)∂λ=2X0(2)∂R0∂λ,∂X1(1)∂λ=2X1(2)∂R0∂λ,∂X2(1)∂λ=2X2(2)∂R0∂λ,∂2R0∂λ2=-[V2S0(1)+∑k=0N-1hk(S1(1)-S2(1))]-1×{2∂R0∂λ∑k=0N-1(S1(1)-S2(1))000×(V2∂hk∂λ+hk∂V2∂λ)000+∑k=0N-1(S1-S2)000×(hk2∂2V2∂λ2+∂hk∂λ∂V2∂λ+V22∂2hk∂λ2)000+2V2(∂R0∂λ)2000×[S0(2)+∑k=0N-1hk(S1(2)-S2(2))]000+S02∂2V2∂λ2+2S0(1)∂R0∂λ∂V2∂λ},S0(2)=exp(-1R02)∑q=02m-2Dq×∑p=0l+q(l+q)!p!R02p+6×{(1-pR02)3+R022[pR04(3p-1)-3]},S1(2)=exp(-k2N2R02)∑q=02m-2Dq×∑p=0l+q(l+q)!k2pp!N2pR02p+6×{(k2N2-pR02)300+R022[pR04(3p-1)-3k4N4]},S2(2)=exp(-(k+1)2N2R02)∑q=02m-2Dq×∑p=0l+q(l+q)!(k+1)2pp!N2pR02p+6×{((k+1)2N2-pR02)300+R022[pR04(3p-1)-3(k+1)4N4]((k+1)2N2-pR02)3},∂S0(1)∂λ=2S0(2)∂R0∂λ,∂S1(1)∂λ=2S1(2)∂R0∂λ,∂S2(1)∂λ=2S2(2)∂R0∂λ.


## 7. Material Dispersion and Profile Parameters

We shall apply a well-known Sellmeier equation to take into account the refractive index dependence on wavelength [4, 25–27, 55, 100]:
(27)n(λ)=1+∑i=13Aiλ2λ2-Bi2,
where Ai and Bi are Sellmeier’s coefficients (Bi is also denoted as the resonance wavelength) which has been empirically measured for GeO2-SiO2 glasses [113, 114] under the several particular dopant concentration. Here we shall apply the method described in [115] to estimate Sellmeier coefficients at the graded-index profile points.A passage from the differentiation operator to its square is used by analogously to propagation constant derivatives. This passage will much simplify expressions for the first and second derivatives of the refractive index:(28)∂n∂λ=12n∂n2∂λ,∂n2∂λ=-2λ∑i=13AiBi2(λ2-Bi2)2,(29a)∂2n∂λ2=12n[∂2n2∂λ2-12n2(∂n2∂λ)2],(29b)∂2n2∂λ2=∑i=132AiBi2(Bi2+3λ2)(λ2-Bi2)3.Therefore the first and second derivatives of the profile height parameterΔcan be expressed as follows:
(30)∂Δ∂λ=12nmax2[(1-2Δ)∂nmax2∂λ-∂nN2∂λ],∂2Δ∂λ2=12nmax2[(1-2Δ)∂2nmax2∂λ2000000-∂2nN2∂λ2-4∂Δ∂λ∂nmax2∂λ].Then derivatives of profile parameterhk defined by formula (9b) are determined by the  following expressions:
(31)∂hk∂λ=1nmax2-nN2[(1-hk)∂nmax2∂λ-∂nk2∂λ+hk∂nN2∂λ],∂2hk∂λ2=1nmax2-nN2[(1-hk)∂2nmax2∂λ20000000000-∂2nk2∂λ2+hk∂2nN2∂λ20000000000-2∂hk∂λ(∂nmax2∂λ-∂nN2∂λ)].Finally the first and second derivatives of the normalized frequency can be written with a help of formulas (28)–(30) in the following form
(32)∂V2∂λ=a2k02[∂nmax2∂λ-∂nN2∂λ-2λ(nmax2-nN2)],∂2V2∂λ2=-2λ∂V2∂λ+a2k02×{∂2nmax2∂λ2-∂2nN2∂λ2000+2λ[1λ(nmax2-nN2)-∂nmax2∂λ+∂nN2∂λ]}.


## 8. Differential Mode Attenuation

We shall apply a simple empirical relation proposed in work [116] using experimental data from [117] to estimate differential mode attenuation depending on mode order: (33a)αμ(λ)=α0(λ)+α0(λ)I9[7,35(μ-1M0)2g/(g+2)],
where μ=2m+l+1 is principal mode number; α0(λ) is the attenuation of lower-order modes (it is supposed equal to the attenuation at the correspondence wavelength mentioned in fiber specification); M0 is total number of modes satisfying the cutoff condition for analyzed fiber:
(33b)M0(λ)=2πnmax(λ)λgΔ(λ)g+2,where g is gradient factor of smoothed α-profile.


## 9. Mode Coupling

Nowadays calculation of mode coupling coefficients is still one of the most well-known and simple methods for optical waveguide junction analysis. It is widely applied for estimation of insertion loss and reflection at the splices of the same type optical fibers with nonidentical technologic parameters [117–121], modeling, and research of the influence of launch conditions on optical waveguide mode composition excitation [25, 122–124] including the modeling of a few-mode optical signal propagation over multimode fibers [4, 8, 10].Generally, the coupling coefficientηmn of injected mode m with excited mode n is determined by the overlap integral of interacting mode fields [4, 8, 10, 117–124]:
(34)ηmn=|∫0∞∫02πΨm(lm)Ψn(ln)rdrdφ|2∫0∞∫02π|Ψm(lm)|2rdrdφ∫0∞∫02π|Ψn(ln)|2rdrdφ.Therefore substitution of radial mode field distribution expressions to integrals (34) is required for estimation of mode coupling coefficients. The simplest mode field structure corresponds to the fundamental mode LP01. As a result, it permits obtaining analytical formulas for the overlap integral (34) by applying simple approximation expressions or by a particular function series expansion for the fundamental mode field representation. These algorithms are widely used for analysis of the single mode fiber splices [25, 117–123].Contrariwise, higher-order mode field structure is more complicated. For this reason in certain papers concerning multimode fiber analysis [4, 8, 10] it is propose to integrate expression (34) numerically or to substitute simple approximation expressions, which are solutions of wave equation for the fibers with ideal step or infinite parabolic refractive index profile. However this substitution is not correct entirely: as it was mentioned above, real refractive index profiles differ much from the model ideal smoothed α-profiles [100] even for the last generation of commercial silica optical fibers both single mode and multimode.Here we propose to apply the overlap integral method combined with introduced EMGA. Under the EMPG application, the use of the algorithm mentioned above for the following mode coupling estimation becomes acceptable: here local features of real silica optical fiber refractive index profile are taken into account that provides decreasing the error of calculations. Thereby the fiber splice analysis will led to substitution of simple approximation expression (6) to the overlap integral (34) under preliminary estimation of equivalent mode field radiuses of injected and excited modes with given azimuthal and radial orders by EMPG. In terms of introduced model of piecewise regular multimode fiber link, this approach provides both simulation of mode mixing and power diffusion by recalculation of the matrix of mode coupling coefficients at the junctions of regular span boundaries and modeling of initial mode composition excitation in multimode fiber by laser source at the transmitter end under taking into account launch condition (single mode pigtail, central launching, offset, air gap, special matching fiber, etc.).Another advantage of proposed method combining the overlap integral method and EMGA is the ability of passage to analytical expressions both for mode coupling coefficients at the ideal centralized splice and splice with offset or tilt for modes with arbitrary order. In recent work [125] the derivation of them is described in detail. Therefore here only final heavy analytical expressions are presented.First of all let us consider the idealized model of diverse optical fiber splice: jointed mismatched type fibers have ideal smooth ends cleaved under 90° to the core axes, and there are no any longitudinal, axial, or angular misalignments between core centers. On the one hand it is very simplified fiber splice description; on the other hand this ideal splice representation corresponds to central launching during the modeling of fiber excitation by laser source.A passage from the generalized form of the overlap integral (1) to analytical expression for the arbitrary order mode coupling coefficient estimation at the central splice of the optical fibers with mismatched parameters and without any misalignments was proposed in work [121]. Analytical expression deriving is described in detail in [125]. Finally the formula for arbitrary order mode coupling coefficient at the central splice is written as follows:
(35)ηmn=Γ(m)Γ(n)Γ(l+m)Γ(l+n)(2ρmρn)2l+2[(ρm-ρn)m+n-2(ρm2+ρn2)m+n+l-1]2×{∑k=0min(m-1n-1)(-1)kΓ(m+n+l-k-1)Γ(m-k)Γ(n-k)k!0000×(ρm2+ρn2ρm2-ρn2)k∑k=0min(m-1n-1)(-1)kΓ(m+n+l-k-1)Γ(m-k)Γ(n-k)k!}.
Here mode coupling occurs only for modes with the same azimuthal order l; Γ is Gamma function; ρm and ρn are injected LPlm and excited LPln mode field radiuses.The analytical expression for arbitrary order mode coupling coefficient at the optical fiber splice under low angle misalignmentθ<10° has the following form [125]:
(36)ηmn(θ)=4(m-1)!(lm+m-1)!(n-1)!(ln+n-1)!×{∑p=0m+n-2∑q=max(0,p-n+1)min(p,m-1)bq(lm,m-1)bp-q(ln,n-1)0000×2(lm/2)+p(k0nθθ)lnΓ((lm/2)+ln+p+1)(ρm2+ρn2)(lm/2)+ln+p+1(ln)!0000×ρm2(ln+p-q)+1ρnlm+ln+2q+10000·F11[ρm2ρn2(ρm2+ρn2)(lm2+ln+p+1),0000000000(ln+1);-(k0nθθ)2ρm2ρn2(ρm2+ρn2)]∑p=0m+n-2∑q=max(0,p-n+1)min(p,m-1)bq(lm,m-1)bp-q(ln,n-1)}2,
where F11 is confluent hypergeometric function of the 1st kind [111, 112]:
(37)F11(a,b;x)=1+∑p=1∞[∏q=0p-1(a+q)x(1+q)(b+q)],
wherebq(lm,m-1)and bq(ln,n-1) are Laguerre polynomial expansion factors of (13a), (13b); nθ is refractive index of launching medium (air gap, core of adjusting/exiting fiber, etc.).The final expression for mode coupling coefficients for 0th azimuth order modesLP0m and LP0n at the splice with transverse offset d has the form [125]:
(38)η0mn(d)=M0·{∑p=0m+n-2∑qbq(0,m-1)bp-q(0,n-1)00000×∑s=0p∑tQt(q,p-q)00000×∑k=02(q-t)(-1)u(2s-1)!![2(q-t)]!u![2(q-t)-k]!k!00000×∑w=0E(u/2)12w(u-2w)!w!·d2(p-s-w)(ρm2+ρn2)s+u-w00000×ρm2(s-q+w)ρn2(s-p+q+u-w)∑p=0m+n-2∑qbq(0,m-1)bp-q(0,n-1)}2,
where the lower and upper limits of second and forth sums are q=max(0,p-n+1)…min(p,m-1) and t=max(0,s-p+q)…min(s,q) and the “odd” factorial is (2s-1)!!=1·3…(2s-1):
(39)M0=4ρm2ρn2(ρm2+ρn2)2exp[-d2ρm2(1-ρn2ρm2+ρn2)],Qt(q,p-q)=q!(q-t)!t!·(p-q)!(p-q-s+t)!(s-t)!,u=k+2(p-q-s+t).The analytic formula mode coupling coefficients for the 1st azimuth order modesLP1m and LP1n at the splice with transverse offset d is written in the following form [125]:
(40)η1mn(d)=4M0mn·{∑p=0m+n-2∑qbq(1,m-1)bp-q(1,n-1)∑s=0p∑tQt(q,p-q)0000000×∑k=02(q-t)+1((-1)u+1(2s-1)!![2(q-t)+1]!0000000000000000×(u+1)!(-1)u+1(2s-1)‼[2(q-t)+1]!)000000000000000×([2(q-t)+1-k]!k!)-10000000×∑w=0E((u+1)/2)12w(u+1-2w)!w!0000000·d2(p-s-w+1)(ρm2+ρn2)s+u+1-w0000000×ρm2(s-q+w)-1ρn2(s-p+q+u-w)+1∑p=0m+n-2∑qbq(1,m-1)bp-q(1,n-1)}2.Finally the expression for mode coupling coefficients for the arbitrary higher-order (l≥2) modes LPlm and LPln at the splice with transverse offset d leads to the following form:
(41)η1mn(d)=M0(m-1)!(l+m-1)!(n-1)!(l+n-1)!×{22l[M1(l)+M2(l)+M3(l)+M4(l)]}2,
where
(42)M1(l)=12∑p=0m+n-2∑q=max(0,p-n+1)min(p,m-1)bq(l,m-1)bp-q(l,n-1)×∑s=0p∑t=max(0,s-p+q)min(s,q)Qt(q,p-q)∑k=0l+2(q-t)Ωk(l)∑w=0E((u+l)/2)Ψw(l),M2(l)=∑μ=0E((l-2)/2)Pμ(l)∑p=0m+n-2∑q=max(0,p-n+1)min(p,m-1)bq(l,m-1)bp-q(l,n-1)×∑s=0p+μ+1∑t=max(0,s-p+q)min(s,q+μ+1)Qt(q+μ+1,p-q)×∑k=0l+2(q-t)Ωk(l)∑w=0E((u+l)/2)Ψw(l),M3(l)=∑μ=0E((l-2)/2)Pμ(l)∑p=0m+n-2∑q=max(0,p-n+1)min(p,m-1)bq(l,m-1)bp-q(l,n-1)×∑s=0p+μ+1∑t=max(0,s-p+q-μ-1)min(s,q)Qt(q,p-q+μ+1)×∑k=0l+2(q-t)Ωk(l)∑w=0E((u+l)/2)Ψw(l),M4(l)=2∑i=02E((l-2)/2)∑j=max[0,E((l-2)/2)]min[i,E((l-2)/2)]Pj(l)Pi-j(l)×∑p=0m+n-2∑q=max(0,p-n+1)min(p,m-1)bq(l,m-1)bp-q(l,n-1)×∑s=0p+i+2∑t=max(0,s-p+q-i+j-1)min(s,q+j+1)Qt(q+j+1,p-q+i-j+1)×∑k=0l+2(q-t)Ωk(l)∑w=0E((u+l)/2)Ψw(l),Pμ(l)=(-1)μ+122μ+3lμ+1(l-μ-2)![l-2(μ+1)]!μ!,Ωk(l)=(-1)u+l(2s-1)!![l+2(q-t)]!(u+l)![l+2(q-t)-k]!k!,Ψw(l)=12w(u+l-2w)!w!·d2(l+p-s-t-w)(ρm2+ρn2)l+s+u-w×ρm2(s-q+w)-lρn2(s-p+q+u-w)+l.


## 10. Results and Discussion

Commercial GeO2-doped graded-index multimode fiber samples of different generations of categories OM2 and OM2+/OM3 with length about 300–500 m packed in coils and conventional spools (Figure 2) were selected for experimental approbations of introduced model.Multimode fiber samples: (a) spools of modern multimode fibers Corning of cat. OM2+/OM3; (b) coil of the first-generation multimode fiber of cat. OM2.
(a)
(b)The refractive index profiles of selected samples were preliminary measured [126] by optical fiber analyzer EXFO NR-9200 via refracted near field scan technique [127]. Profiles reconstructed from the protocol data are presented in Figure 3. Also each fiber sample passed an OTDR test with following fiber length measurement by backscattering method. The OTDR Hewlett Packard E6000A with single mode plug-in unit 1310/1550 was used, and test was produced at the wavelength λ=1310 nm [128, 129].Refractive index profiles measured by EXFO NR-9200: (a) sample no. 01: multimode fiber Corning of cat. OM2+/OM3, length is 373 m; sample no. 02: first-generation multimode fiber of cat. OM2, length is 456 m.
(a)
(b)We applied a typical optical pulse response workstation R2D2 for DMD diagnostics [128, 129]. Simplified block scheme of R2D2 is presented in Figure 4. The experimental setup and software R2D2 screenshot are presented in Figures 5 and 6.Figure 4
Simplified block scheme if pulse response workstation R2D2 and experimental setup.Figure 5
Experimental setup.Figure 6
Software screenshot.The FP LD emits light with a wavelengthλ=1310 nm. The input pulse is shown in Figure 7. It has the full width at a half maximum (FWHM) of about 340 ps. The first series of measurements did not apply special launching conditions. Here only conventional FC/PC adaptors were applied. Therefore an angular misalignment was taking into account θ=4.2° according to [8] during the modeling of transmitter end of link.Figure 7
Input pulse.We simulated fiber irregularity by core diameter variation. It was set from the data of fiber diameter measurements been produced during fiber drawing [130] (Figure 8). Micro- and macrobends were modeled by equivalent low θ=2…3° angular misalignment at the boundaries of regular spans with length Δz=4.79 m [130]. Measured and simulated pulse responses of fibers OM2+/OM3 and OM2 with refractive index profiles shown in Figure 3 are presented in Figures 9(a) and 9(b). Here a close matching can be noticed. Thus for OM2+/OM3 fiber the error of the FWHM is 1.16%, and output pulse dispersion error is 7.42%. For OM2 fiber the errors of second peak amplitude and time scale position is 1.63% and 0.43%, and DMD error is 1.39%.Figure 8
Core diameter variations over fiber length 340 m.Measured and simulated pulse responses: (a) multimode fiber of the OM2+/OM3 category (profile is shown in Figure3(a)); (b) (a) multimode fiber of the OM2 category (profile is shown in Figure 3(b)).
(a)
(b)Another series of tests were produced under the control of launching conditions. Simplified block-scheme of the experimental setup is shown in Figure10. Here emission from the LD launched to tested multimode fiber via single mode pigtail (ITU-T Rec. G.652) which was spliced with multimode fiber with precision offset. Refractive index profile of the single mode fiber is presented in Figure 11. The offset was realized by Ericsson FSU-975 program number 8 “Attenuator” [131]. Here parameter “attenuation” was set to 0, “offset”—to particular value, and “ECF factor” (surface tension compensation factor) also was set to 0 to avoid deformations in the heat zone of the offset fusion splice (Figure 12). Some measured and simulated pulse responses of fibers OM2+/OM3 and OM2 under particular offset launching conditions are presented in Figure 13. A close matching also can be noticed.Figure 10
Experimental setup modified for pulse response measurements under control of launching conditions.Figure 11
Refractive index profile of the single mode fiber.Figure 12
Screenshot of the Ericsson FSU-975 display with offset fusion splice.(a)
Central launching
(b)
Offsetd=18μm
(c)
Offsetd=23,5μm
(d)
Central launching
(e)
Offsetd=5μm
(f)
Offsetd=23,5μm


## 11. Conclusion

An alternative model of multimode fiber links with conventional silica weakly-guiding graded-index irregular multimode fibers under a few-mode optical signal propagation generated by laser source is introduced. The model is based on the piecewise regular representation. It takes into account launch conditions, differential mode delay, both lower- and higher-order mode chromatic dispersion, differential mode attenuation, and mode mixing and power diffusion occurring due to real fiber irregularity and micro- and macrobends. Results of theoretical approbation are presented. A good agreement between developed model and measurements of pulse response were obtained.



---

*Source: 469389-2013-11-12.xml*