Theory of Intermodal Four-wave Mixing with Random Linear Mode Coupling in Few-mode Fibers References and Links

We study intermodal four-wave mixing (FWM) in few-mode fibers in the presence of birefringence fluctuations and random linear mode coupling. Two different intermodal FWM processes are investigated by including all nonlinear contributions to the phase-matching condition and FWM bandwidth. We find that one of the FWM processes has a much larger bandwidth than the other. We include random linear mode coupling among fiber modes using three different models based on an analysis of the impact of random coupling on differences of propagation constants between modes. We find that random coupling always reduces the FWM efficiency relative to its vale in the absence of linear coupling. The reduction factor is relatively small (about 3 dB) when only a few modes are linearly coupled but can become very large (> 40 dB) when all modes couple strongly. In the limit of a coupling length much shorter than the nonlinear length, intermodal FWM efficiency becomes vanishingly small. These results should prove useful in the context of space-division multiplexing with few-mode and multimode fibers.


Introduction
Four-wave mixing (FWM) is an important nonlinear process in the context of silica fibers.Although it was studied in few-mode fibers (FMFs) as early as 1975 [1], until recently most studies on FWM focused on single-mode fibers (SMFs) exclusively [2].From the recent advances in space-division multiplexing (SDM), few-mode, multimode and multicore fibers are seen as candidates for increasing the capacity of telecommunication systems [3].In the last few years, considerable work has been performed, both experimental and theoretical, to study the impact of fiber nonlinearity on SDM systems [4][5][6][7][8][9][10].In particular, Essiambre et  types of FWM processes [8] in an experimental investigation of non-degenerate intermodal FWM (IM-FWM) in FMFs many kilometers in length.Many mechanisms exist through which random coupling among fiber modes can occur [11,12].Long FMFs typically exhibit random linear coupling between modes [13].To fully understand the IM-FWM processes in such fibers, one should consider a model that incorporates the impact of random linear mode coupling between different groups of modes on IM-FWM.
In this paper, we present a theoretical model and perform detailed numerical simulations to understand the impact of random linear mode coupling on the IM-FWM process.In Sec. 2, the two non-degenerate IM-FWM configurations are discussed and the corresponding phasematching conditions are introduced.Section 3 presents the nonlinear propagation equations used in the numerical simulations.More specifically, the effects of random linear mode coupling are incorporated through a random coupling matrix that includes both polarization rotation and linear coupling between different fiber spatial modes.The results obtained in the absence of linear mode coupling are presented in Sec. 4.Under proper phase-matching conditions, both non-degenerate IM-FWM processes exhibit high efficiency.In this uncoupled mode limit, a simple analytical expression is used to estimate the power of the generated idler wave.The bandwidth of the two IM-FWM processes is analyzed in Sec. 5, again in the absence of linear mode coupling.The impact of nonlinear phase shifts of the pumps and the probe on the phase-matching condition of IM-FWM processes is also considered.Numerical results are compared with theoretical predictions and it is found that a simple theoretical model provides a good estimate of the idler power for both IM-FWM processes.Section 6 focuses on the impact of random linear mode coupling and introduces three different coupling models ranging from weak coupling to very strong coupling among all fiber modes.These models are justified by an analysis presented in the appendix that considers the impact of differences of propagation constants between modes on linear random coupling in the limit of infinite fiber lengths.Changes in the FWM efficiency occurring in three linear coupling regimes can be related to the Manakov equations derived for multimode fibers [5][6][7]14].

IM-FWM in FMFs and phase-matching condition
FWM is a process involving up to three optical fields at different frequencies to produce a fourth wave through the third-order nonlinearity [15].In this paper we focus on the so-called non-degenerate IM-FWM and study the generation of an idler wave using three waves distinct in frequencies (see Fig. 1).More specifically, we consider two pumps (labeled pump1 and pump2) at frequencies ω 1 and ω 3 that copropagate inside a FMF with a signal at frequency ω 2 .The idler wave is generated at a frequency ω 4 determined by the general energy conservation relation: where the indices i, j, k, l label the waves.As is well known [2], for an IM-FWM process to become efficient, the following phase-matching condition also needs to be satisfied: where β (r) (ω q ) is the propagation constant of the qth wave propagating in the rth fiber mode.
In practice, the frequencies of the four waves participating in the FWM process are close enough to each other that we can expand all propagation constants in a Taylor series around a common reference frequency ω 0 , often chosen somewhere near the middle of the frequency range involved.For a given mode r, expansion to third order provides where ∆ω = ω − ω 0 and β (r) 0 = β (r) (ω 0 ).The first, second and third derivatives, β 1 , β 2 and β 3 , evaluated at ω 0 , represent the inverse of the group velocity, its dispersion and third-order dispersion, respectively [2].For simplicity, we only consider terms up to second order in this section.For a SMF, all waves involved in the FWM process share the same spatial mode and hence have the same values of β n to all orders.It is easy to deduce that ∆β cannot vanish in SMFs for any non-degenerate FWM process unless one operates at wavelengths near the zero group-velocity dispersion (β 2 ∼ 0).Indeed, FWM in standard SMFs typically requires other techniques to realize phase matching, such as the use of nonlinear phase shifts [16], or the use of birefringent fibers [17].
The situation is different in FMFs whose spatial modes have in general different propagation constants, making it possible to make ∆β = 0 by simply propagating different waves in different fiber modes [18].Here we consider the specific fiber used in the 2013 experiment reported in Ref. [8].It is a graded-index fiber and supports 3 spatial modes, LP01, LP11a, and LP11b in the linearly-polarized (LP) mode approximation [19]; the last two among these are degenerate and share the same propagation constant and are referred together as LP11 modes.Such a FMF is capable of transmitting 6 sets of wavelength-division multiplexing (WDM) channels, because each spatial fiber mode has two orthogonal states of polarization (again in the LP-mode approximation).In order to generate an idler wave through IM-FWM, we consider the configuration shown in Fig. 1 where two waves (probe and pump2) are injected into the LP01 mode, and another wave (pump1) into the LP11a mode.In the experiment, frequencies of the probe and pump2 waves were swept over the entire C-band, and power of the idler wave generated in the LP11 mode was measured.We are interested in two IM-FWM processes, denoted as PROC1 and PROC2 following Ref.[8,20] and shown schematically in Fig. 1.The two processes differ in whether the probe or pump2 is conjugated in the IM-FWM process.The energy relations for PROC1 and PROC2 are described by Eq. ( 1) and, after identifying the waves from Fig. 1, become ) The phase matching conditions have been discussed in Ref. [8], and they turn out to be identical for both processes: where the superscripts 01 and 11 correspond to LP01 and LP11 modes, respectively.The contribution of β 3 is ignored here to simplify the discussion but is included later in (19).Equation ( 6) has a simple interpretation: phase matching of IM-FWM in a FMF is achieved when the group velocities at the average frequencies of the two waves present in each spatial mode are equal.
For the 6-mode FMF considered, based on the measured group velocities [8], equation ( 6) can be satisfied when the waves in the LP11 spatial modes are separated from the two waves in the LP01 spatial mode by about 16 nm.In contrast to SMFs, where effective phase matching can only be achieved at certain pump power levels, full phase matching in a FMF can be achieved at arbitrarily low power levels for a given wavelength configuration.

Nonlinear modal propagation equations
To predict the power of the idler wave generated through IM-FWM, we need to solve a set of coupled nonlinear propagation equations, developed previously by several authors [5,21].The set that we use for this purpose was given earlier as Eq. ( 29) in Ref. [14]: In this matrix equation, A is a column vector with N = 2M elements, where M is the number of spatial modes.The N modes are arranged in a descending order with respect to their propagation constants, i.e., The superscripts T and H denote the transpose and the Hermitian conjugate, respectively.B 0 , B 1 , B 2 , and B 3 are N × N diagonal matrices containing the values of β 0 , β 1 , β 2 , and β 3 for each mode, and all modes are assumed to have the same linear absorption coefficient α.
The nonlinear effects are included through the last term containing the nonlinear parameter γ, defined relative to the fundamental spatial mode [14].Their contributions depend on the mode profiles entering in Eq. ( 7) through two matrices G (1) and G (2) , whose elements are given by G (1) The transverse profile F i (x, y) of the ith mode is normalized such that |F i (x, y)| 2 dx dy = 1, and Γ i j is defined as Here the Kronecker delta function and the modular arithmetic operation lead to Γ i j = 1 when both i and j are odd or even numbers, but Γ i j = 0 in all remaining cases.As can be seen from the form of the nonlinear term in Eq. ( 7), the nonlinear effects depend on the integration of the product of two matrices G (1) and G (2) .Therefore, all nonlinear effects, including IM-FWM, involve spatial profiles of up to 4 modes participating in a specific nonlinear process [14], and are governed by the nonlinear overlap factors: 1 for 4 waves in LP01 mode 0.747 for 4 waves in LP11a or LP11b 0.496 for 2 in LP01 and 2 in LP11a or LP11b 0.249 for 2 in LP11a and 2 in LP11b 0 for all other cases (12) Random linear mode coupling [6,22,23] is included in Eq. ( 7) through a random matrix Q. Considering fluctuations in the refractive index profile, represented in the form of a scalar field ∆ε(x, y, z) as the source of random coupling, the elements q i j (z) of Q, with i, j = 1, 2, . . ., N, can be calculated from the coupled-mode theory using (see [12,Sec. 4.7] and [14, Eq. (40)]) For an ideal fiber, ∆ε(x, y, z) = 0, and there is no linear coupling between different spatial modes since they represent eigenmodes of an ideal fiber.However, random linear mode coupling may need to be included in practice, and its impact is discussed in Sec. 6.

IM-FWM in the absence of linear mode coupling
In this section, we ignore linear mode coupling by setting Q = 0 in Eq. ( 7) and solve this matrix equation numerically for the two IM-FWM processes shown in Fig. 1.In all cases, we assume that both pumps and the probe are launched using continuous-wave (CW) lasers operating in the wavelength range of 1530 to 1570 nm (C-band amplification).We first focus on the PROC1 process and assume that the pump 1 at λ pump1 = 1530 nm is injected into the LP11a mode with 16 dBm (40 mW) power, the pump 2 at λ pump2 = 1554 nm is launched into the LP01 mode with 18-dBm power, and a probe at λ probe = 1546 nm is injected into LP01 mode with 6-dBm power.We choose the reference frequency ω 0 = 2πc/λ 0 to correspond to λ 0 = 1540 nm.The dispersion parameters of our fiber are deduced from measurements using Eq. ( 3) and are summarized in Tab. 1.The LP11 values apply to both LP11a and LP11b modes.When expressed in terms of β 2 = −λ 2 D/(2πc), we obtain β 01 2 = −24.30ps 2 /km and β 11 2 = −23.04ps 2 /km.The differential group delay (DGD) between LP01 and LP11a (or LP11b) spatial modes is d = β 01 1 − β 11 1 ≈ 300 ps/km.All spatial modes are assumed to share the same value of dispersion slope S = 0.055 ps/(nm 2 -km).The nonlinear parameter for our FMF is γ = 1.77W −1 /km, and we include losses in our calculations using α = 0.226 dB/km.We assume that all three waves are polarized in the same direction, along the x axis.
We solve Eq. ( 7) numerically over a distance of 4.7 km and plot the spectrum of the total field at the fiber output.Figure 2 compares the input (solid green lines) and output (dotted red lines) spectra.As seen in the figure, an idler wave is generated in the LP11a mode.Its wavelength of 1538 nm corresponds to what is expected from the PROC1 IM-FWM process.The powers  of pump1 and pump2 are reduced by 1.18 and 1.09 dBm, respectively.Most of this power reduction is due to linear losses since the total loss for our 4.7-km-long fiber is 1.06 dB.A small part of the pump powers is transferred to the probe and idler waves through IM-FWM.Indeed, the 0.51 dB reduction in the probe power is less than the fiber loss of 1.06 dB.The difference of 0.55 dB corresponds to parametric amplification of the probe from the IM-FWM process.
We can estimate the idler power using a simple model.For the configurations shown in Fig. 1 and in the absence of linear coupling (Q = 0), Eq. ( 7) is reduced to the following two equations for the fields in LP01x and LP11ax modes: where Dr is the dispersion operator for the rth mode in the form The nonlinear term that generates the idler wave for both PROC1 and PROC2 is the last term governing intermodal cross-phase modulation (IM-XPM).Since power reduction of the two pumps is mainly due to linear fiber losses, we can employ the undepleted-pump approximation [15].Also, with the assumption of perfect phase matching and CW waves, the idler field A I within the LP11ax mode satisfies 1.53  where A p1 , A p2 , and A B are the amplitudes for pump1, pump2, and probe, respectively.As the preceding equation is linear, it can be easily solved to obtain the idler power For the parameter values used in our numerical simulations, we obtain P I (L) = −4 dBm, which is quite close to the numerical value of −3.78 dB obtained by solving Eq. ( 7).The reason our estimated value is slightly lower than the numerical value is that Eq. ( 17) does not consider parametric amplification of the probe through IM-FWM.
We also performed simulations for the PROC2 IM-FWM process by switching the wavelengths of the probe and the second pump.More specifically, we chose λ probe = 1552 nm and λ pump2 = 1546 nm, following the experimental work [8].We adjusted the value of ∆β 1 slightly (about 0.1%, β 11  1 −β 01 1 = 299.7 ps/km) to ensure that PROC2 is perfectly phase-matched (phase matching in PROC2 is very sensitive to the wavelength and fiber parameters as will be seen in Sec. 5).The results are shown in Fig. 3.As expected, an idler wave is now created in the LP11a mode at a wavelength of 1536 nm through PROC2.The power of this idler wave is quite close to the idler in Fig. 2.This is expected since both processes are perfectly phase-matched in our simulations.Notice that another much weaker idler is also generated in the LP01 mode at 1540 nm.Its origin lies in the cascaded IM-FWM process (of type PROC1) in which the idler wave generated from PROC2 plays the role of the probe.The main reason for the low intensity is that this process is not fully phase-matched.As discussed in Sec. 2, in order to have efficient FWM, the averaged wavelength for the LP01 mode should be about 16 nm longer than that for the LP11 mode.For the cascaded IM-FWM process, the difference in the averaged wavelength is about 10 nm, and this condition is not satisfied.

Bandwidth analysis for PROC1 and PROC2
Since a FWM process requires phase matching, any detuning of the two pumps or the probe from the exact phase-matching conditions will reduce IM-FWM efficiency.This section focuses on the tuning bandwidth of the two IM-FWM processes discussed in Sec. 4. We assume for simplicity that the power levels are sufficiently low that the SPM and XPM do not affect phase matching significantly.Then, the phase-matching bandwidth can be estimated from the linear phase mismatch given in Eq. ( 2).Its use predicts that the idler power in the undepleted-pump approximation is proportional to [15, Eq. ( 2 where ∆β is the linear phase mismatch and L is the fiber length.We stress that η is the phasematching efficiency relative to the ideal phase-matching conditions and should not be confused with the overall efficiency or gain of the underlying FWM process, which relates to the ratio of the generated idler power to the pumps or probe [24].Before we can use Eq. ( 18), we need to calculate ∆β for the two IM-FWM processes considered.Using the numbering of the waves as shown in Fig. 1 and substituting Eq. (3) into Eq.( 2), we obtain an identical expression for both PROC1 and PROC2: where ∆ω i = ω i −ω 0 for i = 1 . . . 4. Setting β 11 3 = β 01 3 = 0 leads to the phase-matching condition given in Eq. ( 6).
Let us assume that PROC1 is phase-matched for a specific probe frequency ω 2 so that ∆β = 0 at that frequency.If we now detune the probe by an amount δ ω, the new probe and idler frequencies become ω 2 = ω 2 +δ ω, and ω 4 = ω 4 −δ ω.Using them, we obtain the phase mismatch due to detuning of the probe frequency by δ ω for PROC1: Although PROC2 shares the same expression for ∆β with PROC1, the impact of probe detuning is quite different for PROC2 because relative locations for pump2 and probe are switched (see Fig. 1).For a detuning δ ω in the probe frequency, the idler frequency now changes to ω 4 = ω 4 + δ ω.As a result, the phase mismatch for PROC2 is given by The two expressions for the phase mismatch can be simplified further by noting that δ ω is a small quantity in practice compared to any ∆ω i (i = 1-3).With this simplification, the phase mismatch for both processes is proportional to the probe detuning δ ω as where the two constants M 1 and M 2 are given by   22).Numerical result (input power P p1 = 4, P p2 = 6, P B = -6 dBm, dotted red curve) agrees well with the analytic approach.
The main difference between the two constant M 1 and M 2 is the sign change of the last term, indicating that phase mismatch depends on the difference in β 2 and β 3 values for the LP01 and LP11 modes for PROC2, while it depends on their sum for PROC1.The differences in β 2 and β 3 values for the two modes are quite small in practice (typically less than 10%).As a result, PROC2 has generally a much larger bandwidth compared to PROC1 because M 1 is much larger than M 2 .
To verify the accuracy of the analytical approximation, we use Eq. ( 18) to calculate the phase-matching efficiency η as a function of λ probe and λ pump2 .The results for PROC1 using Eq. ( 22) are shown in Fig. 4(a).As can be seen, the ridge follows a straight line in the twodimensional plane formed by λ probe and λ pump2 .The slope of the ridge depends on the values of β 2 and β 3 for the two spatial modes, LP01 and LP11a, in which the four waves are propagating.The bandwidth with respect to the probe detuning is quite small (∼0.01 nm) for PROC1.On the other hand, λ pump2 can vary over more than 2 nm. Figure 4(b) shows the efficiency curve for the specific value λ pump2 = 1554 nm (full green curve).We also show the numerical results with a dotted red curve obtained using Eq. ( 7) using low power levels for all waves to rule out the impact of nonlinearity on phase matching.
We performed similar efficiency calculations for PROC2, and the results using Eq. ( 23) are shown in Fig. 5(a).Notice that the x and y axes have been switched in this figure compared to Fig. 4(a).The efficiency map for PROC2 is similar to that of PROC1 in the sense that the peak position follows a straight line in the two-dimensional plane formed by λ probe and λ pump2 .However there is a practical difference related to the switching of the axes.Now, for PROC2,  23).Notice the difference in scale of the λ probe axis and difference in bandwidth of the process.The analytical formula agrees well with the numerical results.
the bandwidth with respect to pump detuning is ∼0.01 nm but the probe wavelength can be detuned more than 2 nm.This feature is evident in Fig. 4(b) where we show the efficiency curve at a specific value λ pump2 = 1546 nm (full green curve).We also show the numerical results with a dotted red curve obtained using Eq. ( 7), using once again low power levels for all waves.The excellent agreement between the two curves for both PROC1 and PROC2 supports well the analytic approximations using Eqs.( 18), ( 22) and (23).
One may ask why PROC2 has a much larger bandwidth (∼2 nm) than that of PROC1 (∼0.01 nm) with respect to probe detuning.A physical reason behind this difference is related to the energy conservation conditions that determine the probe frequency.For PROC1, the idler and the probe frequencies move in the opposite directions, as evident from Eq. ( 1).In contrast, the two waves move in the same direction in the case of PROC2.Clearly the frequency difference ω 2 − ω 4 is preserved in the case of PROC2.It is this property that enhances the bandwidth of PROC2.As the probe wavelength is detuned, the idler wavelength changes in such a way that the phase-matching condition in Eq. ( 6) is maintained over a much wider range compared to PROC1.Using the parameter values for our FMF, we find from Eqs. ( 22) and ( 23) that (∆β ) 1 ≈ 1.4963 × 10 −10 δ ω and (∆β ) 2 ≈ −8.0439 × 10 −13 δ ω.We can estimate the bandwidth from η in Eq. ( 18) by noting that η = 0 for ∆β L/2 = π.Taking L = 4.7 km, it follows that δ ν = δ ω/(2π) ≈ 1.5 GHz for PROC1 and 280 GHz for PROC2.In wavelength units, these numbers correspond to 12 pm and 2.2 nm, respectively.These values agree with the numerical curves in Figs. 4 and 5.
We now discuss briefly the impact of nonlinearity on the phase-matching conditions that increases in importance as the power of the input waves increases.For this purpose, we need to find how the phase of each wave is affected by self-phase modulation (SPM) and crossphase modulation (XPM).In our case, the powers of the two pumps (16 and 18 dBm) are much higher than those of the probe and idler waves (6 and -4 dBm).If we only retain the nonlinear contributions coming from the two pumps, the nonlinear contribution to the phase mismatch for PROC1 and PROC2 is found to be It is interesting to note that nonlinearity affects the two IM-FWM processes differently.The nonlinear contributions to the phase matching increase with increasing pump powers for both pumps for PROC1 but go in opposite directions for pump2 in PROC2.Clearly, PROC2 is generally less sensitive to pump-power-dependent phase matching.The phase-matching condition become insensitive to the pump powers in PROC2 when P p1 /P p2 = f 1111 / f 3333 .
To help assess how nonlinearity impacts phase matching in the experiments, we performed a series of simulations for both IM-FWM processes where we vary the pump powers while keeping other parameters unchanged.We plotted the efficiency curves for PROC1 and PROC2 in Figs.6(a) and (b), respectively.As shown in these two figures, the impact of nonlinearity is mainly to shift slightly the efficiency curves with minimum impact on the bandwidth.As seen from Fig. 6(a), the efficiency curve shifts towards the same direction (shorter wavelength) when either P p1 or P p2 increases.On the other hand, the efficiency curve shifts in the opposite directions when P p1 or P p2 increases for PROC2 (b).These observed behaviors agree with Eqs. ( 26) and (27).By carefully looking at the shift of the efficiency curves in Fig. 6, one can observe that increasing P p2 produces a larger shift than that from increasing P p1 .This behavior can be understood directly from Eqs. ( 26) and ( 27), since the nonlinear overlapping factor f 1111 is greater than f 3333 (For our FMF, f 1111 = 1, f 3333 ≈ 0.75).
Looking from the efficiency curve for PROC2, the amount of wavelength shift between the dashed-green curve [P = (10,10,6) dBm] and solid-purple curve with "+" marker [P = (10,16,6) dBm] is about 0.09 nm.This wavelength shift can be estimated using Eq. ( 27).The nonlinear phase difference between these two curves is γ f 1111 (16 dBm -10dBm) = 0.053 km −1 , which is compensated by an amount of frequency detuning about 65.9 GHz from Eq. ( 23).In wavelength units, this number corresponds to 0.085 nm, quite close to the value 0.09 nm from numerical simulation of Eq. (7).Similar estimation can also be done for PROC1.

Impact of linear mode coupling on IM-FWM
In this section, we consider the impact of random linear mode coupling on the two IM-FWM processes.For SMFs, linear coupling between the two polarization modes originates from the presence of random birefringence fluctuations that lead to random rotations of the polarization state of the field along the fiber length [11,12].The length of fiber over which the state of polarization at the output becomes uncorrelated with its input value (i.e., the correlation length) depends on details of fiber design and fabrication.Experimental measurements on specific fibers have shown values ranging from a few meters to few tens of meters [25,26].In contrast, for fibers supporting multiple spatial modes, linear coupling between two different spatial modes varies greatly depending on the difference in their propagation constants [12,27,28].Many physical effects lead to random linear coupling between modes.Examples of such physical mechanisms are bending, twisting, tension, kinks, pressure, and fluctuations in the fiber core shape and refractive index [11,12].Fluctuations in the refractive index can be modeled by a random quantity ∆ε(x, y, z) that is used to find the elements of the matrix Q(z) in Eq. ( 13).The Fig. 6.Efficiency curves all shift to shorter wavelengths as P p1 or P p2 increase for (a) PROC1 while they shift in opposite directions for (b) PROC2 depending on whether P p1 or P p2 is increased.Powers vector is defined as P = (P p1 , P p2 , P B ).
way polarization and spatial modes couple for a particular fiber depends on the strength and spatial frequency of each coupling mechanism for that fiber.
To evaluate the impact of random linear coupling on nonlinear propagation, we first assume that the linear coupling and phase evolution from difference in refractive indices occur on a length scale shorter than that associated with the nonlinear effects.Under this condition, we can treat the terms B 0 and Q(z) separately from the other terms in Eq. ( 7) and write, The exact way various fiber imperfections impact the coupling matrix Q(z) of multimode fibers still remains an active topic of investigation [28].For our numerical simulations, we assume that the elements of the coupling matrix Q(z) follow Eq. ( 13).Because ∆ε(x, y, z) is a scalar field in Eq. ( 13), the matrix Q(z) does not introduce random coupling between the two polarization states of the same spatial mode.However, as mentioned earlier, we expect several physical mechanisms to randomly couple the two polarization states of the same spatial modes.To take this into account in a model that resembles what occurs in SMFs, we introduce the following matrix R(z) that randomly couples the two polarization states of each spatial mode [29]: where the subscripts p, l, m indicate that all 2 × 2 unitary matrices are independently random.R(z) is also unitary.Note that we omitted writing explicitly the z-dependence of each matrix element for clarity.
In the following numerical simulations, we solve Eq. ( 7) with the split-step Fourier method and incorporate the combined effects of B 0 and Q by multiplying the field after each numerical step ∆z by the matrix Note that in this model, we assume full scrambling at each integration step, and therefore the rotation matrix R is not "proportional" to ∆z.Even though not explicitly studied in this paper, this formulation of coupling allows to separately consider the impact of random linear coupling within each spatial mode from linear coupling between spatial modes.For instance, one can find an analytic derivation of the impact of random linear coupling between the two polarization modes of the same spatial mode alone on the nonlinear propagation equation in [5,14].
The extent of linear coupling between modes depends on many fiber parameters and length of the system considered.In the Appendix, we present a simplified analysis that presents evidence of the existence of three different regimes of coupling.The analysis is based on Eq. ( 28) and considers a 6 × 6 coupling matrix.The three different regimes of propagation are found based on differences in propagation constants between spatial modes, the magnitude of the variances of the elements of Q(z) (all elements are assumed to have identical variances), and the factor by which the fiber length exceeds the correlation length L c of random mode coupling.In the following, we therefore consider three distinct regimes of random linear coupling for T(z).
Model 1: No linear coupling among spatial modes In this case, none of the spatial modes linearly couple, but random linear coupling between the two polarization components of each spatial mode occurs.This simplified model is intended to address the case where all spatial modes have sufficiently distinct propagation constants.In this model, the 6 × 6 transfer matrix T(z) for a 6-mode fiber is identical to R(z) of Eq. ( 29) and assumes the following blockdiagonal form In Eq. ( 31) and equations below, we omit the z-dependence of the elements of the transfer matrices to simplify the notation.
Model 2: Strong random linear coupling between LP11a and LP11b In this case, the degenerate spatial modes LP11a and LP11b experience a strong random linear coupling, while their coupling to the LP01 mode is neglected by assuming a sufficiently large difference in their where the 4 × 4 random unitary matrix takes into account random linear coupling among the LP11ax, LP11ay, LP11bx, and LP11by modes.
Model 3: Strong random linear coupling among all modes In this case, all modes are strongly linearly coupled, and the transfer matrix T 3 (z) is a full 6 × 6 random unitary matrix: We performed numerical simulations for PROC2-type IM-FWM using these three random linear coupling models.The transfer matrix T i (z) (i = 1-3) was applied to the fields at each step of the split-step Fourier method.We propagated the six field components over 4.7 km of the FMF described in Tab. 1.A step size of 7.8 m was used as smaller step sizes generates similar results within the statistical uncertainty of the three models.Numerical simulations performed using Eq. ( 7) are repeated with different random seeds for the generation of the coupling matrices T i (z), with i being the model number.Table 2 shows how the idler power changes over five realizations of the random linear coupling process for the three linear coupling models.As in Sec. 4, the launched pump powers are 16 and 18 dBm, while the probe power is 6 dBm.As seen in Table 2, variations of the idler power for different random matrices realizations are relatively small (< 0.1 dBm) for Models 1 and 2, with an average value close to −7.4 dBm.These results suggest that strong linear coupling between the LP11a and LP11b modes does not affect much the IM-FWM process.
In contrast, in the case of Model 3 where the LP01 and LP11 modes are also linearly coupled, the idler power is drastically reduced and is in the range of −50 ± 2 dBm over five realizations of the random coupling process.Such larger changes indicate that fluctuations in phase  Linear coupling between modes can be reduced in practice by designing fibers with large differences in their modal propagation constants [12], such as between the LP01 and LP11 modes for our FMF.Model 2 would be most realistic in this case.In Sec. 4, the power of the idler was about −3.93 dBm for PROC2 in the absence of any linear mode coupling.The average value of the idler power for Model 2 is about −7.4 dBm.Since random linear mode coupling decreases the idler power by 3.47 dB, it will help in reducing the penalty associated to the interference of the idler with other channels in a WDM system.To estimate the validity of using Model 1 in this section, we used the generalized Manakov equation [Eq.( 26)] of Ref. [14] that includes the coupling model described by the matrix T 1 .The numerical value obtained for the idler power was −7.36 dBm, close to the values obtained in Table 2.
Analytic considerations can also be used to estimate the linear mode-coupling penalty.As discussed in Sec. 4, the nonlinear term responsible for generating the idler wave is proportional to |A 3 | 2 A 1 , which is a "XPM-type" term as described in Refs.[14] and [20].It was shown in Ref. [14] that the factor of 2 associated with the XPM terms in the absence of linear coupling is reduced to a factor 4/3 in the presence of random polarization rotations.As a result, the idler field is reduced by a factor of 2/3, which results in reduction by "(3/2) 2 " or 3.5 dB of the idler power.
In contrast, any linear coupling between the LP01 and LP11 (a or b) reduces the FWM efficiency drastically.The reason for this decrease in idler power is that the phase matching condition can no longer be satisfied.A generalized Manakov equation has also been derived for the strong coupling case (Model 3) [5,14].In fact, as can be seen from this Manakov equation, only "SPM" terms survive after the random averaging.Therefore, it is expected that no efficient IM-FWM happens in such a strong coupling region.Indeed, one can confirm this by performing simulations of Eq. ( 7) with increased number of steps.We found that the idler power did not change for Model 1 and Model 2, but it kept decreasing with the increasing number of steps in the case of Model 3.For a sufficient large step number (i.e., a sufficient small coupling length), the idler power became vanishingly small.Clearly, the results for Model 3 depend on the coupling length of the fiber.For a fixed coupling length, the idler power converges to a stable fixed value that depends on the fiber length.Note also that, in the strong coupling regime, a FMF approaches the behavior of a SMF with the noticeable difference that many modes can occupy the same frequencies.
To verify the power dependence of the IM-FWM processes studied, we reduced the power of all input waves by 3 dB, i.e., P p1 = 13 dBm, P p2 = 15 dBm, and P probe = 3 dBm.All other parameters remain unchanged.The results for PROC2 are summarized in Table 3. Calculations based on the generalized Manakov equation [14, Eq. ( 26)] give -16.31dBm for the idler power.by Q j the random coupling matrix associated to the i th fiber segment of length L c , the total evolution matrix over the fiber length T can then be written as, where we defined p = 1/(∆ L c ) and made use of z = ẑ / ∆ to return to physical units of distance.The ranges of the three variables that impact T in Eq. ( 39) are chosen to be σ ∈ The asymptotic behavior in N of T in Eq. ( 39) for different values of p and σ is studied below.The parameter p can be interpreted as a ratio of a dephasing fiber length (L β = 1/∆ = 2/(β 01 0 − β 11 0 )) to the correlation length of random fluctuations (L c ).The parameter σ is the amplitude of the random coupling coefficient.The ratio σ /p can therefore be loosely interpreted as the ratio of the amplitude of the random coupling to the dephasing length measured in units of correlation length.In the following we prove that, based on the relative values of σ and p, the product in Eq. (39) converges for N → ∞ to a proper subgroup of the group of unitary matrices.To this end we recall the following well-known theorem: Theorem 8.1 ( [30,31]) Consider the group of unitary matrices U(d) of dimension d and a probability measure µ whose support is not contained in a proper subgroup of U(d).Then the measure µ n , obtained by doing the convolution of µ with itself n times, converges in distribution as n → ∞ to the uniform probability measure on U(d).
Specifically we study the following nine different regimes: • p is fixed and p is fixed: In this case, it is straightforward to prove, using 8.1, that the distribution of the product will converge to the uniform measure on the unitary group as N → ∞.
σ p diverges: In this case also, one can prove using Theorem 8.1 that the distribution of the product will converge to the uniform measure on the unitary group as N → ∞ -σ p converges to zero: It is straightforward to prove that the matrix B + σ Q j converges to a diagonal matrix.Hence exp[(i/p)(B + σ Q j )] will also be diagonal.
• p diverges σ p is fixed: In this case, it is straightforward to prove using Theorem 8.1 that the distribution of the product will converge to the uniform measure on the unitary group as N → ∞ -σ p diverges: In this case the matrix B + σ Q j will behave more and more as a GUE matrix, and we can use Theorem 8.1 to conclude convergence to a full Haar matrix as N → ∞.
σ p converges to zero: In this case the argument of the exponential goes to 0, hence the solution converges to the identity matrix.
• p converges to zero σ p is fixed: This case will be analyzed in the rest of this appendix.We prove that exp[ i p (B + σ Q j )] converges to a block diagonal Haar matrix.

Fig. 2 .
Fig.2.Generation of the idler through PROC1 in a 4.7-km-long FMF supporting three spatial modes.Comparison of input (solid green) and output (dotted red) spectra indicates that the idler is created at the expected wavelength of 1538 nm.The idler is generated in the LP11a mode.

Fig. 3 .
Fig.3.Same as in Fig.2but for PROC2.The wavelength of the idler is λ idler = 1536 nm.A second idler is generated in the LP01 mode at λ = 1540nm through cascaded IM-FWM (see main text).

#Fig. 5 .
Fig.5.Same as Fig.4but for PROC2 and with λ pump2 = 1546 nm from Eq. (23).Notice the difference in scale of the λ probe axis and difference in bandwidth of the process.The analytical formula agrees well with the numerical results.

Table 1 .
Dispersion parameters for the LP01 and LP11 spatial modes for the 6 mode FMF.

Table 2 .
Power of the idler (in dBm) generated by the IM-FWM PROC2 for five different realizations of the random linear mode coupling for the three models considered.

Table 3 .
Same as Table2but for 3 dB lower powers for the two pumps and the probe.