1. 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–10]. In particular, Essiambre et al. have observed two 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 phase-matching 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–7, 14].

2. 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 ω₁ and ω₃ that copropagate inside a FMF with a signal at frequency ω₂. The idler wave is generated at a frequency ω₄ determined by the general energy conservation relation:

 

Fig. 1 Wave configurations and notations for two non-degenerate IM-FWM processes: (a) PROC1 and (b) PROC2.

Download Full Size | PDF

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 ω₀, often chosen somewhere near the middle of the frequency range involved. For a given mode r, expansion to third order provides

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:

3. 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]:

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⁽¹⁾ and G⁽²⁾, whose elements are given by

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⁽¹⁾ and G⁽²⁾. 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:

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_ij(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)])

4. 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 ω₀ = 2πc/λ₀ to correspond to λ₀ = 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 β₂ = −λ²D/(2πc), we obtain ${\beta }_2^{01}=-24.30\hspace{0.17em}{\text{ps}}^2/\text{km}$ and ${\beta }_2^{11}=-23.04\hspace{0.17em}{\text{ps}}^2/\text{km}$. The differential group delay (DGD) between LP01 and LP11a (or LP11b) spatial modes is $d={\beta }_1^{01}-{\beta }_1^{11}\approx 300\hspace{0.17em}\text{ps}/\text{km}$. All spatial modes are assumed to share the same value of dispersion slope S = 0.055 ps/(nm²-km). The nonlinear parameter for our FMF is γ = 1.77 W⁻¹/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.

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

View Table | View all tables in this article

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.

 

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.

Download Full Size | PDF

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:

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 Δβ₁ slightly (about 0.1%, ${\beta }_1^{11}-{\beta }_1^{01}=299.7\hspace{0.17em}\text{ps}/\text{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.

 

Fig. 3 Same as in Fig. 2 but 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).

Download Full Size | PDF

5. 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.3.2)]

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:

Let us assume that PROC1 is phase-matched for a specific probe frequency ω₂ so that Δβ = 0 at that frequency. If we now detune the probe by an amount δω, the new probe and idler frequencies become ω′₂ = ω₂ + δω, and ω′₄ = ω₄ − δω. 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 ω′₄ = ω₄ + δω. 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

The main difference between the two constant M₁ and M₂ is the sign change of the last term, indicating that phase mismatch depends on the difference in β₂ and β₃ values for the LP01 and LP11 modes for PROC2, while it depends on their sum for PROC1. The differences in β₂ and β₃ 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₁ is much larger than M₂.

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 two-dimensional plane formed by λ_probe and λ_pump2. The slope of the ridge depends on the values of β₂ and β₃ 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.

 

Fig. 4 (a) Efficiency as a function of λ_probe and λ_pump2 for PROC1 when λ_pump1 = 1530 nm. (b) Efficiency curve for λ_pump2 = 1554 nm (full green curve) from Eq. (22). Numerical result (input power P_p1 = 4, P_p2 = 6, P_B = − 6 dBm, dotted red curve) agrees well with the analytic approach.

Download Full Size | PDF

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, 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).

 

Fig. 5 Same as Fig. 4 but 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.

Download Full Size | PDF

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 ω₂ − ω₄ 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.4963 × 10⁻¹⁰δω and (Δβ)₂ ≈ −8.0439 × 10⁻¹³δω. 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 cross-phase 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

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₁₁₁₁ is greater than f₃₃₃₃ (For our FMF, f₁₁₁₁ = 1, f₃₃₃₃ ≈ 0.75).

 

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).

Download Full Size | PDF

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₁₁₁₁ (16 dBm − 10dBm) = 0.053 km⁻¹, 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.

6. 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 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₀ 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]:

In the following numerical simulations, we solve Eq. (7) with the split-step Fourier method and incorporate the combined effects of B₀ and Q by multiplying the field after each numerical step Δz by the matrix

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 block-diagonal form

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 propagation constants. As a result, the 6 × 6 transfer matrix for this case takes the following block-diagonal form

Model 3: Strong random linear coupling among all modes

In this case, all modes are strongly linearly coupled, and the transfer matrix T₃(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.

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.

View Table | View all tables in this article

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 matching conditions resulting from linear coupling between the LP01 and LP11 modes are so important that IM-FWM is highly suppressed. These results indicate that phase matching cannot be achieves when three nondegenerate waves are launched in a FMF with nearly identical propagation constants.

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₁. 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₃|²A₁, 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)²” 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.31 dBm for the idler power. As can be seen by comparing Tables 2 and 3, the power of the idler generated for all cases is reduced by about 9 dB, a scaling in agreement with Eq. (17).

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

View Table | View all tables in this article

7. Conclusions

In this paper, we have studied a simple analytical model and performed numerical studies of intermodal four-wave mixing (IM-FWM) in few-mode fibers, including birefringence fluctuations and random linear mode coupling. Two different intermodal IM-FWM processes have been investigated. We found that one of these IM-FWM processes has a phase-matching bandwidth about 100 times larger than the other. The impact of nonlinearity on the phase-matching condition and the bandwidth was also studied. We included random linear mode coupling between fiber modes using three different models and found that random linear coupling always reduces the IM-FWM efficiency relative to the absence of any linear mode coupling. The reduction in efficiency is relatively small (on the order of 3 dB) for weakly coupled spatial modes, but becomes very large when all modes are strongly coupled. These results may prove useful in the context of SDM systems designed using wavelength-division multiplexing in multimode and multicore fibers.

8. Appendix: Asymptotic Behavior of Concatenated Random Mode Coupling Matrices

We are interested in the asymptotic behavior of the solution of Eq. (28), where B₀ is a diagonal matrix of the form ${\mathbf{B}}_0=\text{diag}\left({\beta }_0^{01},{\beta }_0^{01},{\beta }_0^{11},{\beta }_0^{11},{\beta }_0^{11},{\beta }_0^{11}\right)$ and Q(z) is a Hermitian random matrix with all entries assumed to have the same standard deviation σ. Note that, in contrast to Eq. (13) of the main text, the model considered for Q(z) includes more general fiber deformations that randomly couple all modes, including those with orthogonal polarizations.

First, we introduce the following changes of variables,

(34)$$\begin{array}{cc}{\mathbf{B}}_0=\mathrm{\Delta }\widehat{\mathbf{B}}-\mathrm{\Sigma }I,& \mathbf{A}=\widehat{\mathbf{A}}\text{exp}(-i\mathrm{\Sigma }z),\end{array}$$

(35)$$\begin{array}{cc}z=\widehat{z}/\mathrm{\Delta },& \mathbf{Q}=\mathrm{\Delta }\sigma \widehat{\mathbf{Q}},\end{array}$$

where $\mathrm{\Delta }=\left({\beta }_0^{01}-{\beta }_0^{11}\right)/2$ and $\mathrm{\Sigma }=\left({\beta }_0^{01}+{\beta }_0^{11}\right)/2$. In terms of the new normalized “hat” variables, Eq. (28) becomes

(36)$$\frac{\partial \widehat{\mathbf{A}}(\widehat{z})}{\partial \widehat{z}}=i\left[\widehat{\mathbf{B}}+\sigma \widehat{\mathbf{Q}}(\widehat{z})\right]\widehat{\mathbf{A}}(\widehat{z}),$$

where B̂ = diag(1, 1, −1, −1, −1, −1) is a normalized diagonal matrix and Q̂ is a normalized Hermitian random matrix with all elements having a standard deviation of 1. Mathematically, this model of Q̂ corresponds to the well known Gaussian Unitary Ensemble (GUE).

We define the evolution matrix in the normalized space T̂(ẑ) as

(37)$$\widehat{\mathbf{A}}(\widehat{z})=\widehat{\mathbf{T}}(\widehat{z}){\widehat{\mathbf{A}}}^{\prime },$$

where Â′ is stationary in ẑ. Using Eqs. (36) and (37), the solution for T̂(ẑ) from ẑ₁ to ẑ₂ is given by,

(38)$$\widehat{\mathbf{T}}({\widehat{z}}_2)=\text{exp}\left[i\left(\widehat{\mathbf{B}}+\sigma 〈\widehat{\mathbf{Q}}(\widehat{z})〉\right)\hspace{0.17em}\hspace{0.17em}\left({\widehat{z}}_2-{\widehat{z}}_1\right)\right]\hspace{0.17em}\hspace{0.17em}\widehat{\mathbf{T}}({\widehat{z}}_1).$$

where 〈Q̂(ẑ)〉 is the z-independent value of Q̂(ẑ), assumed to be nearly constant over the interval ẑ₂ − ẑ₁ ≡ L_c. Let’s separate the fiber of length L into N = L/L_c segments of equal lengths L_c, where L_c is the correlation length over which the entries of 〈Q̂(ẑ)〉 can be considered as constant. Assuming a fiber that is statistically uniform along the length, the entries of 〈Q̂(ẑ)〉 can be considered as independent and identically distributed between distinct segments. Denoting 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,

(39)$$\mathbf{T}=\prod _{j=1}^N\text{exp}\left[(i/p)\hspace{0.17em}\hspace{0.17em}\left(\widehat{\mathbf{B}}+\sigma {\widehat{\mathbf{Q}}}_j\right)\right],$$

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 σ ∈ [10⁻⁵, 10²], N ∈ [10¹, 10⁴] and p ∈ [10⁻⁵, 10⁶].

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_{\beta }=1/\mathrm{\Delta }=2/\left({\beta }_0^{01}-{\beta }_0^{11}\right)$) 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 μⁿ, 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

  • – $\frac{\sigma }{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 → ∞.
  • – $\frac{\sigma }{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 → ∞
  • – $\frac{\sigma }{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

  • – $\frac{\sigma }{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 → ∞
  • – $\frac{\sigma }{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 → ∞.
  • – $\frac{\sigma }{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

  • – $\frac{\sigma }{p}$ is fixed: This case will be analyzed in the rest of this appendix. We prove that $\text{exp}\left[\frac{i}{p}(\mathbf{B}+\sigma {\mathbf{Q}}_j)\right]$ converges to a block diagonal Haar matrix.
  • – $\frac{\sigma }{p}$ diverges: In this case the matrix B + σQ_j will behave increasingly as a GUE matrix as $\frac{\sigma }{p}$ increases, and we can use Theorem 8.1 to conclude convergence to a full Haar matrix as N → ∞.
  • – $\frac{\sigma }{p}$ converges to zero: We expect several different behaviors to be exhibited depending on the rate of convergence of σ to 0.

We now prove that when p = σ and σ → 0, the matrix $\text{exp}\left[\frac{i}{p}(\mathbf{B}+\sigma {\mathbf{Q}}_j)\right]$ converges to a block diagonal matrix. To this end, we divide the individual entries of this matrix into three groups: diagonal entries, off-diagonal entries within a block, and off-block entries. The size of the blocks is dictated by the structure of the matrix B.

Off-block entries

Let us suppose σ ≪ 1. We can divide the eigenvalues of $\frac{1}{p}(\mathbf{B}+\sigma {\mathbf{Q}}_j)$ into two subsets: the eigenvalues S_Δ close to Δ and the eigenvalues S_−Δ close to −Δ. Let μ ∈ S_Δ (respectively S_−Δ)) and v_μ the associated eigenvector. By writing down the associated eigenvalue eigenvector decomposition, we see that the last four (the first two) coordinates are of the order of σ. Since the off block entries of $\text{exp}\left[\frac{i}{p}(\mathbf{B}+\sigma {\mathbf{Q}}_j)\right]$ are of the form ${\sum }_{\lambda }{e}^{i{\lambda }_j}{v}_{\mu }\overline{v_{{\mu }^{\prime }}}$ with μ and μ′ in subsets attached to different eigenvalues, we conclude that the magnitude of an off block entry is upper bounded by nr where n is the dimension of the matrix. Hence, what remains to be shown is that the magnitude of the block entries does not vanish when σ = p → 0.

Block entries

We will show that the eigenvector matrix of $A=\frac{1}{p}((\mathbf{B}+\sigma {\mathbf{Q}}_j)$ is block diagonal in the limit σ = p → 0, with non vanishing block entries. We prove this for $C={(1/p)}^2\left[\sigma (\mathbf{B}{\mathbf{Q}}_j+{\mathbf{Q}}_j\mathbf{B})+{\sigma }^2{\mathbf{Q}}_j^2\right]$ instead of A, as the two statements are equivalent. Indeed, C is A² minus a scalar matrix, which does not alter the eigenvectors, i.e., C, A and e^iA all have the same eigenvectors. Here, we use the fact the matrix A is Hermitian to be able to apply functional calculus to A.

Observe that the matrix BQ_j + Q_jB is block diagonal with 2 blocks of the desired dimensions. This is because of the particular structure of B. We will show that in the limit as σ = p → 0, the eigenvectors of C converge to those of BQ_j + Q_jB. This is the content of the following theorem, which bounds the Frobenius norm of the eigenvector matrix.

Theorem 8.2 ([32, 33]) Let M ∈ ℂ^n×n be a Hermitian matrix with distinct eigenvalues and have spectral decomposition M = UΛU^*. Let M̃ = M + δM be a Hermitian matrix and

$$\begin{array}{cc}\epsilon =\frac{{\Vert \delta M\Vert }_F}{{\text{min}}_{i\ne j}\left|{\lambda }_i-{\lambda }_j\right|},& \alpha =\frac{\left(1+\sqrt{1-1/n}\right){\Vert M\Vert }_2}{{\text{min}}_{i\ne j}\left|{\lambda }_i-{\lambda }_j\right|}.\end{array}$$

If $\epsilon \left(2\alpha +\sqrt{1+4{\alpha }^2}\right)\le 1$, then M̃ has a spectral decomposition such that

(40)$${\Vert \tilde{U}-U\Vert }_F\le \frac{\sqrt{2}\epsilon }{{\left[1-2\alpha \epsilon +\sqrt{1-4\alpha \epsilon -{\epsilon }^2}\right]}^{1/2}}.$$

We apply this theorem with $M=\frac{1}{p}(\mathbf{B}{\mathbf{Q}}_j+{\mathbf{Q}}_j\mathbf{B})$ and $\delta M={\mathbf{Q}}_j^2$. Hence U is the eigenvector matrix of $\frac{1}{p}(\mathbf{B}{\mathbf{Q}}_j+{\mathbf{Q}}_j\mathbf{B})$ and Ũ the eigenvector matrix of C. Fix p and observe that the eigenvalues of $\frac{1}{p}\left(\mathbf{B}{\mathbf{Q}}_j+\mathbf{B}{\mathbf{Q}}_j\right)$ are distinct with probability 1 since the density of the eigenvalue distribution is continuous. As p → 0, min_i≠j|λ_i − λ_j| and ‖A‖₂ grow as 1/p. Using this, a straightforward computation shows that ε → 0 and αε ≤ 1 for p small enough, and the theorem is satisfied. Since the right side of Eq. (40) decreases to 0, the Frobenius norm also converges to 0. This implies that the entries of Ũ converge to those of U. Since Q_j is a GUE matrix, we can say that the blocks in the eigenvector matrix of BQ_j + Q_jB (and therefore also those of C) will be Haar distributed.

Having this result in hand, we can show that e^iA is block diagonal, which is the result we seek. All off-block entries converge to 0 as stated earlier. Consider a block entry (e^iA)_lm. The eigenvalue/eigenvector decomposition of A gives that this entry is of the form ${\sum }_k{e}^{i{\lambda }_k}{v}_k{w}_k$ where v and w are respectively the l^th and m^th eigenvectors of A. Observe that the only nonzero entries of v_k and w_k are those corresponding to the block under consideration. By a standard result in perturbation theory (B and Q_j being hermitian, see [34]), the eigenvalues of $\frac{1}{p}(\mathbf{B}+\sigma {\mathbf{Q}}_j)$ are close to those of B/p, and the distance between the corresponding eigenvalues is of the order of σ/p. Since σ/p = 1, we conclude that the eigenvalues λ_k are distinct modulo 2π with probability 1. Hence, the sum ${\left(e^{iA}\right)}_{lm}={\sum }_k{e}^{i{\lambda }_k}{v}_k{w}_k$ is non zero. Putting all these results together, we conclude that in the limit p = σ → 0, e^iA is block diagonal. Since multiplication respects the block structure, we use Theorem 8.1 to conclude that the blocks of the matrix will be Haar distributed in the limit as N → ∞.

We show examples of the shape of the matrix T in Eq. (39) by varying parameters p, σ and N. In all numerical experiments, we consider an ensemble with 1000 members, each representing a random draw of the matrix T, and compute the empirical expected value of the absolute value of each entry. Figure 7 shows the gray-coded values of the matrix elements of T using N = 1 and σ = p for values of σ ranging from 10⁻¹ to 10⁻⁴. We observe two matrix blocks of sizes 2 × 2 and a 4 × 4, and a strong diagonal over the entire 6 × 6 matrix. The bottom row shows the situation for N = 10. We clearly see two blocks; the distribution of the elements in each block become much more uniform for N = 10 compared to the N = 1 case. These numerical results show behavior predicted by the asymptotic analysis in this appendix.

 

Fig. 7 Gray-coded values of the elements of the matrix T with white and black representing 0 and 1, respectively. The values of p with p = σ are shown above each plot. Top row N = 1; bottom row N = 10.

Download Full Size | PDF

Finally, Figure 8 shows an example where σ is fixed and p decreases using N = 10. In this situation, we do not obtain a block structure, and the matrix has its full dimension. As N → ∞, the product in this case converged to a uniform distribution over the whole matrix.

 

Fig. 8 Gray-coded values of the elements of the matrix T for N = 10 when σ = 1 is fixed and p decreases from 0.1 to 10⁻⁴.

Download Full Size | PDF

Acknowledgments

This work was supported in part by Alcatel–Lucent in the framework of GreenTouch. We acknowledge an anonymous reviewer for the insightful comments.

References and links

1. R. H. Stolen, “Phase-matched-stimulated four-photon mixing in silica-fiber waveguides,” IEEE J. Quantum Electron. 11, 100–103 (1975). [CrossRef]  

2. G. P. Agrawal, Nonlinear Fiber Optics, 5th ed. (Academic Press, New York, 2012).

3. R. Ryf, S. Randel, A. H. Gnauck, C. Bolle, A. Sierra, S. Mumtaz, M. Esmaeelpour, E. C. Burrows, R.-J. Essiambre, P. J. Winzer, D. W. Peckham, A. McCurdy, and R. Lingle Jr., “Mode-division multiplexing over 96 km of few-mode fiber using coherent 6 × 6 MIMO processing,” J. Lightwave Technol. 30, 521–531 (2012). [CrossRef]  

4. C. Koebele, M. Salsi, G. Charlet, and S. Bigo, “Nonlinear effects in mode division multiplexed transmission over few-mode optical fiber,” IEEE Photon. Technol. Lett. 23, 1316–1318 (2011). [CrossRef]  

5. A. Mecozzi, C. Antonelli, and M. Shtaif, “Nonlinear propagation in multi-mode fibers in the strong coupling regime,” Opt. Express 20, 11673–11678 (2012). [CrossRef]   [PubMed]  

6. S. Mumtaz, R.-J. Essiambre, and G. P. Agrawal, “Reduction of nonlinear penalties due to linear coupling in multicore optical fibers,” IEEE Photon. Technol. Lett. 24, 1574–1576 (2012). [CrossRef]  

7. A. Mecozzi, C. Antonelli, and M. Shtaif, “Coupled Manakov equations in multimode fibers with strongly coupled groups of modes,” Opt. Express 20, 23436–23441 (2012). [CrossRef]   [PubMed]  

8. R.-J. Essiambre, R. Ryf, M. A. Mestre, A. H. Gnauck, R. W. Tkach, A. R. Chraplyvy, Y. Sun, X. Jiang, and R. Lingle Jr, “Experimental investigation of inter-modal four-wave mixing in few-mode fibers,” IEEE Photon. Technol. Lett. 25, 539–541 (2013). [CrossRef]  

9. J. Demas, P. Steinvurzel, B. Tai, Y. Chen, and S. Ramachandran, “Two octaves of frequency generation by cascaded intermodal nonlinear mixing in solid optical fiber,” Proc. Conf. on Lasers and Electro-optics (2013). Paper CTu2E.5.

10. A. D. Ellis, N. Mac Suibhne, F. C. Gunning, and S. Sygletos, “Expressions for the nonlinear transmission performance of multi-mode optical fiber,” Opt. Express 21, 22834–22846 (2013). [CrossRef]   [PubMed]  

11. S. Rashleigh, “Origins and control of polarization effects in single-mode fibers,” J. Lightwave Technol. 1, 312–331 (1983). [CrossRef]  

12. D. Marcuse, Theory of dielectric optical waveguides (Academic Press, 1991).

13. J. M. Fini, B. Zhu, T. F. Taunay, M. F. Yan, and K. S. Abedin, “Crosstalk in multicore fibers with randomness: gradual drift vs. short-length variations,” Opt. Express 20, 949–959 (2012). [CrossRef]   [PubMed]  

14. S. Mumtaz, R.-J. Essiambre, and G. P. Agrawal, “Nonlinear propagation in multimode and multicore fibers: generalization of the Manakov equations,” J. Lightwave Technol. 31, 398–406 (2013). [CrossRef]  

15. R. W. Boyd, Nonlinear optics (Academic Press, 1992).

16. Z. Su, X. Zhu, and W. Sibbett, “Conversion of femtosecond pulses from the 1.5-to the 1.3-μm region by self-phase-modulation-mediated four-wave mixing,” J. Opt. Soc. Am. B 10, 1050–1053 (1993). [CrossRef]  

17. R. Stolen, M. Bösch, and C. Lin, “Phase matching in birefringent fibers,” Opt. Lett. 6, 213–215 (1981). [CrossRef]   [PubMed]  

18. R. Stolen, J. Bjorkholm, and A. Ashkin, “Phase-matched three-wave mixing in silica fiber optical waveguides,” Appl. Phys. Lett. 24, 308–310 (1974). [CrossRef]  

19. D. Gloge, et al., “Weakly guiding fibers,” Appl. Opt 10, 2252–2258 (1971). [CrossRef]   [PubMed]  

20. R.-J. Essiambre, R. W. Tkach, and R. Ryf, “Fiber nonlinearity and capacity: Single-mode and multimode fibers,” in “Optical Fiber Telecommunications VI B,”, I. Kaminow, T. Li, and A. E. Willner, eds. (Academic Press, 2013), chap. 1, pp. 1–37. [CrossRef]  

21. F. Poletti and P. Horak, “Description of ultrashort pulse propagation in multimode optical fibers,” J. Opt. Soc. Am. B 25, 1645–1654 (2008). [CrossRef]  

22. J. M. Fini, B. Zhu, T. F. Taunay, M. F. Yan, and K. S. Abedin, “Statistical models of multicore fiber crosstalk including time delays,” J. Lightwave Technol. 30, 2003–2010 (2012). [CrossRef]  

23. C. Antonelli, A. Mecozzi, M. Shtaif, and P. J. Winzer, “Random coupling between groups of degenerate fiber modes in mode multiplexed transmission,” Opt. Express 21, 9484–9490 (2013). [CrossRef]   [PubMed]  

24. M. E. Marhic, Fiber optical parametric amplifiers, oscillators and related devices (Cambridge university press, 2008).

25. M. Wuilpart, P. Mégret, M. Blondel, A. J. Rogers, and Y. Defosse, “Measurement of the spatial distribution of birefringence in optical fibers,” IEEE Photon. Technol. Lett. 13, 836–838 (2001). [CrossRef]  

26. A. Galtarossa and L. Palmieri, “Spatially resolved PMD measurements,” J. Lightwave Technol. 22, 1103–1115 (2004). [CrossRef]  

27. R. Olshansky, “Mode coupling effects in graded-index optical fibers,” Appl. Opt. 14, 935–945 (1975). [CrossRef]   [PubMed]  

28. L. Palmieri, “Coupling mechanism in multimode fibers,” in Proc. SPIE 9009, Next-Generation Optical Communication: Components, Sub-Systems, and Systems III (2013). Paper 90090G.

29. C. Poole and R. Wagner, “Phenomenological approach to polarisation dispersion in long single-mode fibres,” Electronics Letters 22, 1029–1030 (1986). [CrossRef]  

30. P. Harremoës, “Maximum entropy on compact groups,” Entropy 11, 222–237 (2009). [CrossRef]  

31. E. Breuillard, “Random walks on lie groups,” Survey, http://www.math.u-psud.fr/breuilla/part0gb.pdf.

32. G. W. Stewart and J. guang Sun, Matrix Perturbation Theory, Computer Science and Scientific computing (Academic Press, 1990).

33. X. S. Chen, W. Li, and W. W. Xu, “Perturbation analysis of the eigenvector matrix and singular vector matrices,” Taiwanese J. Math. 16, 179–194 (2012).

34. T. Kato, Perturbation theory for linear operators, vol. 132 of Grundlehren der mathematischen Wissenschaften (Springer, 1995).