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

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 lowerand higherorder mode chromatic dispersion, differential mode attenuation, and mode mixing and power diffusion occurring due to real fiber irregularity and microand 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.


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][2][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][2][3][4].Because emission from the conventional VCSEL usually consists of about 5 or 6 transversal modes   with maximal azimuthal order  not more than 3, and single-mode LD injects just fundamental and lower-order modes  01 and  11 , only several guided modes are excited in multimode fiber link [5][6][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][9][10][11][12][13][14][15]-should also be noted: their authors are affiliated to IEEE 802.3z, 802.3ae, 802.3aq, and 802.3baTask Forces.Other works are also concerned with the design and simulation of multi-Gigabit multimode optical communication systems.Thus model [16][17][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][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][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 gradedindex 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 wellknown 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, patchpanels, 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 gradedindex 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.

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-section  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 .
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][29][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][26][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 .
In [32,33], mode components of few-mode signal are combined in mode groups with roughly equal propagation constants.In works [34][35][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][38][39].Also in [40][41][42][43][44] a special curvilinear coordinates are applied to describe mode coupling via fiber bends.In [45][46][47], mode mixing and power diffusion is simulated by combination of coupled mode theory and raytracing analysis.
The system of incoherently coupled nonlinear Schrodinger equations is applied for mathematical description of pulse propagation over multimode fiber in [48][49][50][51][52][53][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.
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 (10 3 , . . ., 10 4 ) ⋅ -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.
From this point of view, method described in [96][97][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 higherorder 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.

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 Δ.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.
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 by  mode components   over regular multimode fiber with length  [4, 8-15, 18, 99-101]: where  is direct Fourier transform; ℎ  () is initial pulse at the transmitter end;  (0)  and   are starting amplitude and mode attenuation of th guided mode   ( = 1, . . ., );  () 1 and  () 2 are first-and second-order dispersion parameters.These dispersion parameters are elements of wellknown Taylor series expansion approximation of propagation constant frequency dependence () [2,55,99]: where Here  () 1 =  () is mode delay and  () 2 is group velocity dispersion (GVD) associated with chromatic dispersion of th guided mode   .
Therefore, according to introduced piecewise regular representation of irregular multimode fiber link, the frequency response of a few-mode optical signal transferred by  guided modes over irregular multimode fiber with length  under given particular length of regular span Δ by taking into account expression (1) can be written in the following form: where   = (/Δ); () is the integer part of real number .
The resulting pulse response at the receiver end of irregular multimode link is computed by the following simple expression: where 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 axialsymmetric 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.

Extension of Modified Gaussian Approximation
EMGA is based on classical Gaussian approximation [25] of radial mode field distribution  ()  () 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: where  () −1 is Laguerre polynomial.EMGA leads to equivalent normalized mode field radius  0 estimation by solving a characteristic equation, which is derived from propagation constant variational expression under following passage to square core mode parameter  2 variational expression, written for analyzed weakly guiding optical waveguide with given refractive index profile.Parameter  0 is basic for this method and completely defines mode transmission parameters.
Unlike known methods and their modifications [25,[104][105][106][107][108][109][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  layers in which the refractive stays a constant [102,103]: and any profile function () related to refractive index profile () as where Δ = ( 2 max −  2  )/2 2 max is profile height parameter, can be written in terms of profile parameter ℎ  : where where   is refractive index of  layer ( = 0, . . ., );  max is the maximal core refractive index;   is cladding refractive index.This refractive index profile representation permits writing the variational expression for core mode parameter  2 and characteristic equation for normalized equivalent mode field radius  2 / 0 = 0 in the form of finite nested sums as follows [15,16]: where and  (,)  is coefficient of polynomial representation in the form of power series [111,112]: =  0  max √ 2Δ is normalized frequency;  0 = 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  layers.Fiber parameters and mode orders  and  are substituted to the characteristic equation (11).By means of numerical solution (11), the normalized equivalent mode field radius  0 will be obtained.Then  0 is substituted to the expression (10), and mode core parameter  is estimated; that permits evaluating the propagation constant  for guided mode   by the wellknown expression [4, 25-27, 55, 100]: Solution of the characteristic equation ( 11) is correct under normalized frequency  > 1, and it should satisfy the guided mode cutoff condition [4, 25-27, 55, 100]: Optical confinement factor  core can be considered as the second criterion for the identification of the ghost solutions: If we take into account Gaussian approximation, parameter  core is defined by analytical expression derived in [15] from the generalized integral form for weakly guiding optical fibers presented in [25]: where

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 velocity V  is related to propagation constant  by the well-known formula [4, 25-27, 55, 100]: Advances in Optical Technologies 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: After substitution of ( 20) to (19), an expression for the mode group velocity can be rewritten in the following form: Then by differentiating expressions (10) and ( 11) with respect to wavelength, first derivatives of square core mode parameter  2 and normalized equivalent mode field radius  0 will be obtained: ℎ  ( (1)  1 −  (1)  2 )} ,

Mode Chromatic Dispersion Parameter
The same technique is applied to derive the analytical expressions for mode chromatic dispersion parameter  associated with GVD.Chromatic dispersion parameter  is related to propagation constant  by the well-known formula [4, 25-27, 55, 100]: By differentiating (20), a second partial derivative of the square of propagation constant will be obtained: Therefore expression for  leads to following form: By differentiating expressions (22) with respect to wavelength, required second derivatives of square core mode parameter  2 and normalized equivalent mode field radius  0 will be obtained: 0 1 2 0 1 2

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]: where   and   are Sellmeier's coefficients (  is also denoted as the resonance wavelength) which has been empirically measured for GeO 2 -SiO 2 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: Therefore the first and second derivatives of the profile height parameter Δ can be expressed as follows: Then derivatives of profile parameter ℎ  defined by formula (9b) are determined by the following expressions: Finally the first and second derivatives of the normalized frequency can be written with a help of formulas ( 28)- (30) in the following form

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) where  = 2 +  + 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);  0 is total number of modes satisfying the cutoff condition for analyzed fiber: where  is gradient factor of smoothed -profile.

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][118][119][120][121], modeling, and research of the influence of launch conditions on optical waveguide mode composition excitation [25,[122][123][124] including the modeling of a fewmode optical signal propagation over multimode fibers [4,8,10].Generally, the coupling coefficient   of injected mode  with excited mode  is determined by the overlap integral of interacting mode fields [4,8,10,[117][118][119][120][121][122][123][124]: 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  01 .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][118][119][120][121][122][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 Advances in Optical Technologies    formula for arbitrary order mode coupling coefficient at the central splice is written as follows: Here mode coupling occurs only for modes with the same azimuthal order ; Γ is Gamma function;   and   are injected   and excited   mode field radiuses.
The final expression for mode coupling coefficients for 0th azimuth order modes  0 and  0 at the splice with transverse offset  has the form [125]: where the lower and upper limits of second and forth sums are  = max(0,  −  + 1) . . .min(,  − 1) and  = max(0,  −  + ) . . .min(, ) and the "odd" factorial is (2 − 1)!! = 1 ⋅ 3 . . .(2 − 1): The analytic formula mode coupling coefficients for the 1st azimuth order modes  1 and  1 at the splice with transverse offset  is written in the following form [125]: Finally the expression for mode coupling coefficients for the arbitrary higher-order ( ≥ 2) modes   and   at the splice with transverse offset  leads to the following form: where

Results and Discussion
Commercial GeO 2 -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.
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].
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.
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.
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).Microand macrobends were modeled by equivalent low  = 2 . . . 3 ∘ angular misalignment at the boundaries of regular spans with length Δ = 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%.
Another series of tests were produced under the control of launching conditions.Simplified block-scheme of the experimental setup is shown in Figure 10.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.

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.

Figure 1 :
Figure 1: Piecewise regular representation of 200 m irregular multimode optical fiber with varying core diameter.

Figure 4 :
Figure 4: Simplified block scheme if pulse response workstation R2D2 and experimental setup.

Figure 9 :
Figure 9: Measured and simulated pulse responses: (a) multimode fiber of the OM2+/OM3 category (profile is shown in Figure 3(a)); (b) (a) multimode fiber of the OM2 category (profile is shown in Figure 3(b)).

Figure 10 :Figure 11 :
Figure 10: Experimental setup modified for pulse response measurements under control of launching conditions.