Theoretical study of gain-induced mode coupling and mode beating in few-mode optical fiber amplifiers.

We developed a generalized field-propagating model for active optical fibers that takes into account mode beating and mode coupling through the amplifying medium. We applied the model to the particular case of a few-mode erbium doped fiber amplifier. Results from the model predict that mode coupling mediated by the amplifying medium is very low. Furthermore, we applied the model to a typical amplifier configuration. In this particular configuration, the new model predicts much lower differential modal gain than that predicted by a classical intensity model.


Introduction
In addition to Dense Wavelength Division Multiplexing (DWDM), Polarization Division Multiplexing (PDM) and multi-level modulation format, the most promising solution, for enhancing data transmission rates beyond those of SMF technology capacity crunch, seems to take benefit of an unused degree of freedom, namely space [1]. Two technical approaches based on Spatial Division Multiplexing (SDM) are now largely studied, namely multicore Fibers (MCF) and Mode Division Multiplexing (MDM) using Few-Mode Fibers (FMF). As a ground-breaking technology, MDM prompts a need to re-adapt each element of an optical transmission line, among which the amplifiers which are essential in long-haul lines where they are located at regular intervals [2].
In order to adapt Erbium Doped Fiber Amplifiers (EDFAs) to MDM context and so develop efficient Few-Mode EDFAs (FM-EDFAs), numerical simulations are essential tools for design and optimization. To do so, many different models have been developed, and most of them are based on the pioneer model proposed by Giles and Desurvire [3] with generalization to the FMF case [4][5][6][7]. These models consider a two-level system to describe the energy level configuration of Er 3+ ions, Linearly Polarized (LP) modes and have shown a good agreement with experimental measurements. The models describing amplification in single-mode or strongly multimode fibers are largely based on intensity as the evolving z-dependent quantity, which is sufficient in most cases. Situation can be different in few-mode systems where modal competition involving a well-defined number of modes suggests that field models be developed so as to get a real and accurate knowledge of the problem. Such a model is particularly relevant to describe modal competition in fiber laser systems that are rarely single-mode or to evaluate quantities such as Differential Modal Gain (DMG) in few-mode optical amplifiers that are now developed for SDM. However, as it has been previously explained, it is important to have a better understanding of phenomena which take place during amplification and especially in multimode systems, so as to design a FM-EDFA. Some fine phenomena such as global transverse intensity evolution (induced by mode beating) and mode coupling are known to occur in such a system. To our knowledge, two models neglecting noise and counter-propagating beams, considering only mode beating effects and using vector mode basis have been proposed by Lim et al. Another model proposed by Nasiri Mahalati et al. [12], similar to the one presented in this paper, considers both mode beating and mode coupling in FM-EDFA, but doesn't take into account counter-propagating beams and noise. In this paper, by using equations based on electric fields in a complete amplifier model, both mode beating and coupling through gain medium are taken into account, together with co and counter-propagating beams and noise. Such a model gives, for example, the opportunity to predict cross-talk induced by the amplification media.
The purpose of the paper is: • Develop field model (including mode beating), • Estimate gain-induced cross-talk, • Highlight difference between field model and classical intensity model in a typical FM-EDFA configuration.
In Sec. 2 intensity model will be briefly recalled. Then, in Sec. 3 fundamentals and approximations used for field model will be reviewed in order to describe the architecture of the model itself in Sec. 4. In the last section, some results obtained with this model will be compared to results coming from an intensity FM-EDFA algorithm.

Reminders on intensity model
So as to describe the amplifier, an intensity model generally considers a classical two-level system for Erbium ions for which population in the two levels is driven by absorption, stimulated emission and spontaneous emission. For classical intensity models, propagation equations (Eq. (1)) and population equations (Eq. (2)) in the steady-state approximation are presented as a reminder here below [7]: where N 2 (x, y, z) is the density of ions in the upper energy level ( 4 I 13/2 ), N (x, y, z) is the erbium doping profile, σ a,k and σ e,k are respectively the absorption and emission cross sections at wavelength λ k , I k,i (x, y) stands for the normalized transverse intensity repartition at λ k for spatial mode i, φ k ,i ,± (z) is the photon flux (at λ k ) which pass through the transverse section at z in co-propagating (+) or counter-propagating (−) direction, and A 21 is the relaxation rate associated to spontaneous emission. The above model is insufficient for taking field effects into account. In the next section, we present a field model that accommodates field effects.

Fundamental and approximations of the field model
The generalized model presented here can be applied to many different optical fiber amplifiers and adapted to fiber lasers. Moreover, it is compatible with any eigenvectors basis. Since our research activity is focused on MDM and FM-EDFA, the model is here used in this case in order to evaluate the impact of the amplification mechanism on the degradation of the signal quality in MDM systems as accurately as possible.
In the following, a low refractive index contrast will be considered and spatial component of the electric field ( − → E (x, y, z)) will be written with the formalism : where i is the index of the eigenmode, k is the index of the considered wavelength λ k , z is the spatial coordinate along propagation axis (both co-and counter-propagating beams are considered), β k ,i , A k,i (z) and − → F k,i (x, y) are respectively the propagation constant of mode i at wavelength λ k , the amplitude of the electric field related to the power carried by this mode and the normalized transverse field distribution for this mode (real in the low refractive index contrast assumption). The initial phase shift is set as zero for simplicity, which means that at z=0 all modes are launched in-phase.

Mode beating
Since the different modes do not propagate at the same velocity and so accumulate different phase shifts, the instantaneous global transverse intensity profile at position z is given by the result of their mutual interferences, which is the definition of mode beating. The total transverse intensity distribution (I k (x, y, z)) for the coherent combination of modes of arbitrary amplitudes, at the same wavelength λ k , is given by : Before presenting an example of mode beating, some reminders about the propagation of modes in optical fibers will be made. Conventionally, approximation of LP mode basis (well adapted for low index contrast) is used to model propagating modes. However, the exact solutions of Maxwell's equations for modes propagating in an optical fiber, thus called vector modes, are divided into three families : Transverse Electric (TE), Transverse Magnetic (TM) and hybrids (HE and EH). These modes can be combined (in-phase, i.e. propagation phase shift ∆ βz = 0 [2 π]) to reconstruct the well-known LP l,m modes, in the low refractive index contrast assumption: This procedure will be used in Sec. 5 to compare the solutions of two algorithms based on LP mode basis (intensity model) and vector mode basis (field model). To simplify notation in the following, the propagation phase shift ∆ βz will be noted Θ(z). The resulting expression of I k (x, y, z) for two vector modes constituting a LP mode can be used as an example to illustrate this phenomenon. For example, if interferences between arbitrary mode 1 and mode 2 are considered, Eq. (4) becomes then : where Θ k,1,2 (z) = ∆ β k ,1,2 z = ( β k ,2 − β k,1 )z. This equation is similar, for example, to the one obtained for Spatially and Spectrally resolved imaging (S 2 ), since it considers the same phenomenon [13]. An example of mode beating within a pseudo-LP mode between its vector mode constituents, considering a step index weakly guiding fiber, is presented with HE even 31 ( Fig. 1(a)) and EH even 11 ( Fig. 1(b)) modes summed to reconstruct the LP a, x 21 mode when Θ(z) = 0 [2 π] ( Fig. 1(c)). We denote as pseudo-LP modes the coherent addition of the vector modes constituent fields of a LP mode when Θ(z) 0 [2 π]. When Θ(z) = 0.5 π [2 π], the resulting transverse intensity exhibits a doughnut shape ( Fig. 1(d)). Then, the modes pursue their propagations and reach Θ(z) = π [2 π]: the transverse intensity then became that of the LP b,y 21 mode ( Fig. 1(e)). At Θ(z) = 1.5 π [2 π] ( Fig. 1(f)), conclusions similar to that at Θ(z) = 0.5 π [2 π] can be drawn. Finally, when the accumulated phase shift reaches Θ(z) = 2 π [2 π] (z is an integral multiple of a beat length) the transverse intensity and electric field become again identical to those at the fiber input ( Fig. 1 (c)). This phenomenon is not limited to pairs of vector modes constituting a LP mode, and occurs between all modes at a specific wavelength. This can result in many different transverse intensity distributions, no matter the mode basis. As can be guessed, mode beating can greatly alter photons repartition over propagation and so have a significant impact on FM-EDFA behavior. One aim of this study is to evaluate the impact of this phenomena on FM-EDFA modeling.

Mode coupling, propagation and population equations
In the field model, equations presented in Sec. 2 and used in the intensity model are modified. Light at pump wavelength (λ p ) is considered to be monochromatic and to propagate on n(λ p ) different modes. For wavelength range used for signal, light is considered as Ω s independent optical beams that can propagate on n(λ s ) different modes. Over this wavelength range, signals are attributed to particular modes at specific wavelengths λ s (some discrete λ s values in the so-called C and L bands). For pump beam and signal beams (in this second case, at the same wavelength λ s ), since light is coherent, effects of mode beating and coupling are taken into account. Regarding propagation equations, the model uses coupled differential equations for signal and pump beams, by supposing a slowly varying envelope A. It has to be noticed that only the impact of amplification on the imaginary part of the susceptibility is considered in our theoretical model. For λ k corresponding to λ p or λ s , n being the number of guided modes at λ p or λ s , field propagation equations for the field model are written: Mode coupling then appears in the equation through the terms g k,i , j in Eq. (6b), that represents, to simplify, the overlap integrals between electric fields distributions of two modes through the population inversion "seen" by these modes. Amplified Spontaneous Emission (ASE) is considered as a sum of Ω e incoherent light beams at different wavelengths (λ e ∈ [1500 1600] nm) propagating on n(λ e ) different modes, that are processed independently. In this context, local intensity for ASE at λ k is found by summing incoherently the n(λ e ) different electrical fields : Then, propagation equation used to model ASE is the same as for intensity model and can be written as follows : Regarding population equations, for pump and signal beams, in Eq. (2), the term x, y, z) (including scale factor), in accordance with Eq. (4), and is dependent of position z now that intensity profile depends on phase shifts between all the different modes existing at λ p or λ s and their respective amplitude. This leads to the general population equations for field model : Now that both population and propagation equations for signal, pump and ASE have been introduced, the architecture of the model will be presented in the next section.

FM-EDFA modeling
In order to solve Eqs. (6), (8) and (9), numerically, we need to divide the three spatial dimensions into elements of small size. In other words, the optical fiber is divided in N x × N y × N z elements of volume. So, the algorithm makes a longitudinal, transverse, spectral and modal resolution of the system of equations. Note that longitudinal resolution dz is shorter than the smallest beat length.
For each step z, the different incoming photons flux φ k ,i ,± (z) are used to calculate the transverse distribution of ground-and excited-state populations at z position. Then after these populations are used to calculate photon flux at z + dz. Eqs. (9) are used to calculate population inversion for each coordinate (x, y). Once populations at z are calculated, propagation equations are used to determine the photon flux φ k ,i ,± (z + dz) at z + dz with a fourth-order Runge Kutta method. To do so, Eq. (6) is used for beams corresponding to signal and pump, whereas Eq. (8) is used for ASE.
The algorithm is divided into two parts : initial guess and convergence. Each part is further subdivided into two steps : co-propagating and counter-propagating beams.
Firstly, calculation of populations and flux are done in co-propagating direction with initial conditions at z = 0 (co-propagating pump power coupled on each mode, co-propagating signal power carried by each mode at each wavelength λ s , noise term for ASE). Then, these approximate populations and intensities of co-propagating beams (pump, signal and ASE) are considered fixed in order to evaluate counter-propagating beams with initial conditions at z = z ma x . This constitutes the initial guess, the results of which having to be refined with convergence loop.
Secondly, before entering in convergence loop, populations are evaluated with both copropagating and counter-propagating beams from z = 0 to z = z ma x . The convergence loop evaluates alternatively co-propagating and counter-propagating beams with population re-evaluation between each step. Convergence is reached when photon flux for iteration j (φ j k ,i ,± ) are close to photon flux for iteration j − 1 (φ j −1 k ,i ,± ), the criterion of convergence being define by the user. This criterion also lays down dx, dy, dz values, and so acts on computation time-accuracy ratio.
Then, numerical results provide the longitudinal evolution of : • the density of ions in the upper energy level : N 2 (x, y, z), • the photon flux for signal, pump or ASE co-propagating and counter-propagating beams : These results are used to calculate gain for each signal modes i at each wavelength λ k . Specific parameters can be analyzed such as average gain (G ave ), Differential Modal Gain (DMG(k)) that represents the maximum gain difference between modes at the same wavelength λ k , Differential Spectral Gain (DSG(m)) that represents the spectral gain excursion for one mode m and differential gain (∆G) that merges both DMG and DSG. These numerical results also permit to calculate Signal to Noise Ratio (SNR) and Noise Figure (NF).

6-modes FM-EDFA
In this section, the amplification properties of a 6-mode Few-Mode Erbium Doped Fiber are presented. In this example, refractive index (∆n = 9.7.10 −3 ) and Erbium doping (concentration of 14 × 10 24 ions/m 3 ) profiles are step profiles. This fiber has a circular core radius of 7.5 µm, supports 34 vector modes at 980 nm (pump wavelength) and 12 vector modes at 1550 nm. The amplifier length is chosen as 2.6 m. Signal is multiplexed at wavelengths from 1530 to 1560 nm with 1 nm channel spacing, and 20 µW input power at every wavelength for each signal vector mode, corresponding to 40 µW for each equivalent LP mode. Signal modes are chosen to be HE odd 11 and HE even 11 , HE odd 21 and TE 01 , HE even 21 and TM 01 equivalent respectively to LP 01 , LP 11a and LP 11b . The FM-EDFA is co-propagating pumped with 600 mW total pump power, 300 mW in each of the HE odd 21 and TE 01 modes, equivalent to LP 11a mode. This configuration has been chosen to highlight, with a single example: • effects of mode beating on the calculated gain for the 3(6) LP(vector) modes, • effects of mode coupling induced by the gain medium on the 3(6) LP(vector) modes unused as signal channels (i.e. with zero input power).
First of all, results extracted from the model previously introduced are presented. Firstly, the longitudinal evolutions of terms g k ,i , j (Eq. (6b)).
Looking at the term g k ,i , j , only some coupling coefficients of modes i with all modes j are presented, these modes i being the one with nonzero input power. These coefficients g k,i , j have  Fig. 2. Longitudinal evolution of the term g p,i , j presented in Eq. (6b) at pump wavelength (λ p = 980 nm) for modes TE 01 (a) and HE odd 21 (b). Only modes n corresponding to higher values of this term are identified. It has to be noted that some curves are superposed and are identified together with the symbol "&".
higher values for modes with the same orientation (odd or even) and the same parity l. TE 0m and TM 0m modes are two particular cases that couples mainly to HE odd/even 2m . This is illustrated in Fig. 2   wavelengths, the longitudinal evolutions of terms g k ,i ,i (coupling coefficient for one mode with itself) is mainly driven by longitudinal evolution of population inversion. As can be observed in Eq. (6b), the g p,i ,i term is negative at λ p , reflecting pump absorption (σ e , p = 0). This is depicted in Figs. 2(a) and 2(b)). By contrast, the g s ,i ,i term is positive (at λ s ) as can be seen in Figs. 3(a)-(f), since signal is amplified if population inversion is nonzero. Concerning coupling with other modes (g s,i ,n i ), coefficients oscillate around zero due to periodic phase matching conditions between modes. This coupling is alternatively power loss or noise coming from another channel.
Looking more in details at modal powers evolutions at pump wavelength (Fig. 4), it can be noticed that the most important coupling occurs with HE even 21 and TM 01 modes. This would appear to contradict conclusions extracted from Fig. 2, since these two modes do not correspond to higher values of coefficient g k ,i , j . However, this can be explained by the values of β for these two modes which are close to the values of modes chosen for pump beam (introduced in Eq. (6a) through phase term exp( j ( β k,n − β k ,i ) z)). It is also due to the fact that the modal power evolution presented in Fig. 4 is the result of the contributions from all the modes while Fig. 2 presents evolutions of term g k ,i , j for one particular mode. It must be pointed out that coupling coefficient g k ,i , j is different from coupling phenomenon and represents the overlap between electric fields of two different modes and population inversion as written in Eq. (6b), it could be then simplified as a coupling tendency coefficient. In order to consider coupling phenomenon, phase matching and powers carried by each modes must be considered as in Eq. (6a). Similar observations can be made for the different modes at signal wavelengths (Fig. 5). However, the most important point to underline is that mode coupling due to amplification generated by the FM-EDFA is low for both pump and signal beams (ten order of magnitude lower).
One can note that mode beating has an effect on both g k ,i , j and modal power AND manifests itself as oscillations in the evolutions of these quantities. This effect has been presented in Fig. 1 and will be further studied in the following.
The results are compared to those obtained with a classical intensity model. This comparison is possible by "reconstructing" the LP modes with the procedure described in Sec. 3.2. The longitudinal evolutions of normalized average N 2 population is presented on Fig.6, whereas longitudinal evolutions of mode powers at pump wavelength and signal wavelengths are presented on Fig. 7 and 8. For these two last figures, in order to simplify observation only modes with   nonzero input power are presented since mode coupling is low. The purpose of this comparison is to highlight that the predictions of the intensity model and those of the field model can differ significantly for some specific configurations. As such, the results present the outcome of both models using the same parameters and the same initial conditions, which incidentally corresponds to a realistic FM-EDFA configuration. The normalized inversion profile on Fig. 6 exhibits oscillations along the first meter of the fiber due to beating between HE odd 21 and TE 01 pump modes. The period of these oscillations corresponds to the beat length of these modes (L beat = 0.36 m). Beyond this length, oscillations are more complex because signal modes become powerful enough to have a significant impact on population inversion while pump power decreases. This effect is more and more pronounced when the fiber length increases. The longitudinal evolution of pump modes presented in Fig. 7 permits to observe a slight difference between the two models. This difference is a result of both beating and coupling that modify the interaction between pump beam and erbium ions. The homogeneous distribution of ions dopant explains the small change in the absorption of pump beam. Regarding the propagation of signals beams represented in Fig. 8, results are very different between the two simulations. Indeed, it can be easily observed that with the intensity model, LP 11a mode is the most amplified mode and LP 11b is the less amplified one, with an important gap between them. With the field model, curves representing modal power evolution for these two modes are almost superposed. Mode beating is mainly responsible for these differences as was the case for N 2 population. Indeed, for the intensity model the overlap integrals between the pump beam and signal modes remain constant during propagation. In this example, the lower value of this integral being with the signal mode LP 11b , the highest with LP 11a mode. With the field model presented in this paper, mode beating at pump and signals wavelengths changes the values of these integrals along propagation, generating the observed differences. This leads to the situation that signal modes corresponding to LP 11 modes have very similar average overlap integrals with the pump beam. The LP 01 mode is the mode the less affected by these changes, since its transverse intensity profile has a rotational symmetry. The differences for the fundamental mode, contrary to LP 11 modes, are not due to the modification of the overlap integrals for this mode but are the consequences of changes for LP 11 modes explained above. The consequence is the shift in the balance of modal gain within the amplifier, that results in higher gain for the fundamental mode as will be seen below.
It is important to note that the amplification is performed via erbium ions and therefore the overlap integrals should also consider their transverse spatial distribution. In this example, erbium ions are distributed evenly, thereby simplifying these integrals in order to facilitate the understanding of the phenomenon and its consequences. The effect of mode beating can also be observed in the evolution of signal power beams presented in Fig. 8. The previous differences are also reflected in the gain values. In the case of the intensity model presented in Fig. 9(a), the modal gain is, in ascending order : G L P 11b < G L P 01 < G L P 11a , which is consistent with previous observations made about overlap integrals. For this simulation, average gain, G ave , is 17.1 dB, ∆G is 16 dB. DSG values for four representative wavelengths and DMG values calculated with both models are reported respectively in Table 1 and Table 2. With the new code described in this paper, the results are very different ( Fig. 9(b)), as suggested by previous  observations. Indeed, average gain G ave is 18.9 dB, ∆G is 3.0 dB. So, 13 dB ∆G difference between the two models is reported. In this case, mode beating effect acts as an intrinsic gain equalizer. It appears that taking into account the effects of mode beating and mode coupling can greatly modifies the results of the amplifier modeling. In the case presented here, the simulated FM-EDFA would be declared to have very poor performances if evaluated with the intensity code. Considering these effects with the field model presented here, its performances are declared good for the different reasons explained previously. Hence, in the amplifier configuration tested in this paper, the intensity model gives misleading results when applied to FM-EDFAs. It is important, for each FM-EDFA configuration, to check that effects usually neglected by intensity models do not have a significant impact on amplifier performances.

Conclusion
In conclusion, a generalized field model for few-mode fiber amplifiers including co-and counterpropagating beams and noise has been developed. This model is compatible with any eigenvector basis (vector modes, LP modes,...), and takes into account mode beating and mode coupling through the gain medium. It has been applied to a particular case of FM-EDFA. On one hand, the comparison with the usual intensity model shows that the gain medium is not an important source of coupling in FM-EDFA. The results also suggest that the change on the real part of the refractive index due to amplification will also have low impact on mode coupling since it is of the same order of magnitude than the modifications considered in our model for the imaginary part.
On the other hand, the effect of mode beating could have an important impact on the gains (13 dB difference on ∆G) and the exact solutions of Maxwell's equations must be used in order to take into account all the possible mode beatings.To determine if mode beating effects should be considered or not, the global amplifier configuration should be taken into account. The amplification is driven by the transverse overlap integral between pump intensity profile, erbium doping profile and signal modes transverse intensity profiles. Every effect that modifies these integrals along the propagation axis will have an effect on gain. It is difficult to determine precisely the parameters that have to be considered to evaluate if the field model should be used or if the intensity model is sufficient. However, one can imagine that if the transverse population inversion induced by the pump beam(s) does not have a rotational symmetry, mode beating effects should be considered. It is the case in the example presented in this paper with LP 11a (or HE odd 21 and TE 01 modes) pump mode. Phase modulated signals will strongly decrease the impact of this effect between signal modes but will still occur between vector modes constituting a LP mode, if LP mode basis is used for signal modulations. This would decrease the impact of mode beating, but not enough for it to be negligible. Mode beating is obviously more pronounced at the pump wavelength. In many cases, the coherence length of the pump laser diode could be as small as few millimeters. In this case, the effect of mode beating will only occur at the beginning of the amplifier where population inversion is higher and thus can have a significant impact on gain.
The example chosen here depicts a core pumping scheme that is probably the case for which the mode beating is exacerbated. However, mode beating still occurs at the pump wavelength in the cladding pumping configuration but the large number of modes in this case and their low overlap integrals with the doped area would probably reduce its impact on gain.
Moreover, in practice, a very low core ellipticity is sufficient for vector modes to become linearly polarized. Whatever, as it has been previously written, the model presented here is compatible with any eigenmode basis and effects such as mode beating or coupling are still present [10,11].