1. Introduction

In optical communication fibers, polarization-mode dispersion (PMD, i.e., group birefringence), chromatic dispersion (CD), and Kerr nonlinearity are the well known effects causing signal distortion. In linear region (i. e., without Kerr nonlinearity), the signal distortion due to PMD and CD can be predicted and compensated. On the other hand, in the limit of pure nonlinear effect (i.e., without linear effects such as PMD and CD), the signal amplitude will not be affected and the pulse shape will be unchanged along the fiber. Unfortunately both linear and nonlinear effects in real fibers are not negligible. To accurately evaluate the signal distortion, one must treat them simultaneously [1–5] by solving the coupled nonlinear Schrödinger Eq. (CNLSE), which was proposed in 1980’s [6] and further studied in many publications.

For a non-soliton system, the CNLSE is solved with the step size determined by dispersion length L_D, nonlinear length L_N, as well as the birefringence related parameters such as fiber correlation length (L_corr) (the length within which birefringence axes are randomly reoriented), beating length Λ_beat (the length determined by birefringence strength), and the PMD parameter D_PMD $\left(\mathrm{ps}/\sqrt{\mathrm{km}}\right)$ [1–3]. Ignoring the birefringence effect, the two dimensional (2D) CNLSE can be reduced to one dimensional nonlinear Schrödinger Eq. (1D NLSE), with its step size being only related with L_D and L_N. To efficiently solve the 2D CNLSE and 1D NLSE, various approaches have been proposed (cf. Refs. [1–15] and the references therein).

In optical communication fibers, the values of Λ_beat (10~100 m) and L_corr (0.3~300 m) are much smaller than L_D and L_N (up to hundreds of kilometers) [1–3]. Dealing with the rapidly and randomly varying birefringence in the multispan communication fibers, decreasing the computational step size to a very small percentage of L_corr and Λ_beat not only needs too much computation but also cannot ensure good enough accuracy because of the computational error accumulated in all steps. To efficiently solve the CNLSE, one needs an accurate approach (or algorithm) with large enough step size.

In Ref. [2], a coordinate system rotating with the principal axes in each wave plate was introduced to get the analytical solutions for the linear and nonlinear effects in CNLSE. As a result, the computational step size using this approach can be increased to a scale significantly larger than Λ_beat and L_corr. (Detailed value of its step size depends on the interaction between linear and nonlinear effects.) Assumptions required by Ref. [2] are 1) the circular component of the local birefringence in the fiber is negligible and 2) the nonlinear term in the CNLSE can be approximated as the SPM (self-phase modulation) modified by a factor of 8/9. For telecom fibers, the linear local birefringence approximation is good enough [1], while the SPM approximation, also called Manakov-PMD approximation (M-PMD approx) in this work, needs the condition that the polarization state on the Poincaré sphere should be fully mixed [1, 2].

Like many approaches used to solve 2D CNLSE or 1D NLSE, another essential feature of Ref. [2] is that it needs split-step Fourier method (SSFM), which is based on the first order approximation of Baker-Hausdorff formula [16]

In Eq. (1), [A, B] ≡ AB − BA. Because of this, even in the case of zero birefringence, the step size using approach [2] is restricted to 10% ~ 15% of L_D and/or L_N, assuming the computational error tolerance is less than 1% of the pulse peak (cf. the results in Sec. 4.2).

The concept of interaction picture (IP) was originally used in quantum mechanics. Combined with the fourth-order Runge-Kutta (RK4) [17], the RK4 in IP method (RK4IP) was recently proposed to study the supercontinuum generation in optical fibers (for 1D field cases) [15]. As the IP method does not introduce the error inherent in various SSFMs, the computational error of RK4IP is determined by the accuracy of RK4.

The approach of Ref. [15] cannot be applied directly to the cases with random birefringence, because in Ref. [15], the linear dispersion operator in 1D NLSE was required to be constant along the fiber. Motivated by the analytical solution for the linear terms in CNLSE, which was proposed in Ref. [2] and is further discussed in the Appendix for our calculation, we extend the RK4IP approach from 1D NLSE to 2D CNLSE.

As there is no need to rotate the coordinate system, our approach can be applied to the fibers having the general form of local birefringence. Moreover, it treats Kerr nonlinearity without approximation. So, 2D RK4IP does not need to concern the mixing condition required by M-PMD approx. With the help of the local error method of [12], the local error of our RK4IP can be improved to ~ O(h⁶), rather than ~ O(h⁵) of 1D RK4IP [15]. As results, for the communication fibers with random birefringence, the step size using RK4IP can be orders of magnitude larger than L_corr and Λ_beat. Without birefringence, the step size of RK4IP can be the same order of magnitude as L_D and/or L_N, or, 6 ~ 10 times the step size using the approach of Ref. [2].

Under the condition that the polarization state is mixed enough on the Poincaré sphere, our RK4IP result agrees well the M-PMD approx, which was used in many publications (cf. e.g., [1–3]). In the case of strong linear-nonlinear interaction, as the mixing on the Poincaré sphere may not be enough, the discrepancy between RK4IP and M-PMD approx can be observed (cf. Fig. 4 and Fig. 6).

2. From CNLSE to Manakov-PMD approximation

2.1. CNLSE expressed in different forms

To deal with the birefringence related problem, it is convenient to represent a 2D optical field with Dirac bra or ket notation [18]. Namely, given a field with its x − y components being u_x and u_y, it can be denoted as ∣u〉 ≡ (u_x,u_y)^T [or 〈u∣ ≡ (u_x^*,u_y^*)]. Thus, the CNLSE discussed in [1–3] can be written in the following form with retarded time t = t_lab − β_ωz:

or

where ∣u〉_z ≡ ∂ ∣u〉/∂z and Δβ (Δβ_ω ≡ ∂Δβ/∂ω) is related to Λ_beat (L_corr) with $\mathrm{\Delta }\beta =\pi /{\mathrm{\Lambda }}_{\mathit{beat}}\left(\mid \mathrm{\Delta }{\beta }_{\omega }\mid =D_{\mathit{PMD}}/\sqrt{8L_{\mathit{corr}}}\right)$, respectively (cf. Eq. (28) of [2] or P. 6 of [19] and others). Here D_PMD is the PMD coefficient $\left(\mathrm{ps}/\sqrt{\mathrm{km}}\right)$. When L >> L_corr, the average DGD of a L-km-long fiber can be obtained using DGD_avg = D_PMD√L (ps). Obviously, the second term on the left-hand side of Eq. (3) represents the phase birefringence, while the third term relates to the group birefringence (or linear PMD). In Eq. (2), ∣v〉 ≡ (u_x^*u_y²,u_x²u_y^*)^T is often called DFWM (degenerate four-wave mixing), while the first term in the nonlinear part of Eq. (3) is the SPM (self-phase modulation) term (cf. e.g., [20]).

In this work, the three components of the Pauli spin matrices are denoted as

Generally, the birefringence-induced matrix ∑ in Eq. (2) and Eq. (3) has the form of [3]

with β_i ≡ 〈β₀(z)∣σ_i∣β₀(z)〉 (i = 1,2,3) and ∣β₀(z)〉 [β⃗₀(z)] being, respectively, the unit vector representing the fiber birefringence in 2D Jones (3D Stokes) space [18]. For the fiber with circular birefringence (e.g., spun optical fiber), β₂ ≠ 0 in Eq. (5).

Parameter β_ωω = ∂²β/∂ω² corresponds to the CD parameter at wavelength λ by β_ωω=-CD(λ)λ²/(2πc) (ps/nm·km) with c=3 × 10⁸m/s.

In Eq. (3), ${\sigma }_{31}=-{\sigma }_{13}={\sigma }_3{\sigma }_1=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)=j{\sigma }_2$, with σ_i (i = 1, 2, 3) being given by Eq. (4). Obviously, for a unitary transformation, e.g.,

we have T^†σ₁₃T = σ₁₃, provided T = T^*. This means the nonlinear term in the CNLSE is covariant for any real rotation. Or, the nonlinear effect in the CNLSE is independent of the artificial reference frame we choose. In the following sections, further discussions on CNLSE are based on the form of Eq. (3).

When the amplitude of ∣u〉 is viewed as the square root of optical pulse power, the nonlinear coefficient γ in Eq. (2) and Eq. (3) relates to Kerr coefficient n₂, effective mode area A_eff, and wavenumber k = 2π/λ by γ = n₂k/A_eff. Typical values of λ = 1550 nm and n₂ = 2.6 × 10^-20m²/W are used throughout of the paper, unless otherwise noted.

2.2. Manakov-PMD approximation

Equation (3) can also be expressed as

which was named Manakov-PMD equation in Refs. [2, 3]. The last term on the left-hand side of Eq. (7) is the Kerr nonlinearity averaged over the Poincaré sphere with the 8/9 factor [1, 2, 21]. The right-hand side of Eq. (7) is the nonlinear PMD due to incomplete mixing on the Poincaré sphere [2, 3]. Physically it corresponds the rapidly varying fluctuations as the polarization state changes [2].

In this work, the M-PMD approx means that the variation of the second term on the right-hand side of Eq. (7) is approximated by its statistical average over the Poincaré sphere [the first term on the right-hand side of Eq. (7)]. Thus, the right-hand side of Eq. (7) is approximated as zero, yielding

Without linear birefringence, M-PMD approx [Eq. (8)] reduces to Manakov Eq. [1, 2, 21].

The form of Eq. (8) is different from, but equivalent to, the M-PMD approx proposed in Refs. [1–3]. In fact one can obtain the published M-PMD approx, e.g., Eq. (68) of Ref. [3], by substituting the unitary transformations ∣u〉 = RT∣$\tilde{\Psi }$ 〉 into Eq. (8), with T being given by Eq. (6)and $R=\left(\begin{array}{cc}\cos \alpha & -\sin \alpha \\ \sin \alpha & \cos \alpha \end{array}\right)$, yielding $R_z={\alpha }_z\left(\begin{array}{cc}-\sin \alpha & -\cos \alpha \\ \cos \alpha & -\sin \alpha \end{array}\right),j{R}^{\mathrm{†}}{R}_z={\alpha }_z{\sigma }_2,$ and [σ₂α_z + σ₃Δβ]T + jT_z = 0.

3. RK4IP: extension from 1D to 2D

3.1. RK4IP solution for 1D NLSE

Without birefringence, Eq. (3) can be reduced to 1D NLSE

Equation (9) can be formally viewed as

Introducing transformation u = e^j(z-z₀)D̂u^I, NLSE in the IP [i.e., Eq. (10)] has the form [15]

Equation (11) is equivalent to

with

N̂^I in Eq. (13) is the nonlinear operator represented in the IP. Given u(z_n,t), the next step field u(z_n+1,t) (z_n+1 = z_n + h) governed by Eq. (12) can be obtained using RK4, with the local error of O(h⁵). According to RK4 [17], we have u^I(z_n+1,t) = u^I(z_n,t) + h[k₁ + 2(k₂ + k₃) + k₄]/6 with k₁ = f[z_n,u^I(z_n,t)], k₂ = f[z_n + h/2,u^I(z_n,t) + k₁h/2], k₃ = f[z_n + h/2,u^I(z_n,t) + k₂h/2], and k₄ = f[z_n + h,u^I(z_n,t) + k₃h]. Setting z₀ = z_n, one needs 16 FFTs to get u(z_n+1,t) from u(z_n,t): 12 FFTs for k₂, k₃, k₄, and 4 FFTs for u(z_n,t) → u^I(z_n,t) and u^I(z_n+1,t) → u(z_n+1,t).

Choosing z₀ = z_n + h/2, one only needs 8 FFTs, yielding [15]

In Eq. (14), there is no FFT for k₂ and k₃. Also, unlike k₁, k₂, and k₃, k₄ in Eq. (14) has been transformed back from the IP.

Note that in Eq. (13) and Eq. (14), the linear operator D̂ for 1D case was assumed to be unchanged in z direction. For a fiber with random birefringence, the transformation Î_D(z,z₀) ≡ exp[j(z - z₀)D̂] = exp[-j(z - z₀)β_ωω(∂²/∂t²)/2] in Eq. (13) needs to be modified correspondly, which is the key point to extend the RK4IP of [15] to 2D case.

3.2. RK4IP solution for 2D CNLSE with random birefringence

Equation (3) can also be viewed as

Introducing the transformation ∣u_n〉 = Î(z_n, z₀)∣u^I_n〉, with the operator Î(z, z₀) given by Eq. (27), the CNLSE [Eq. (15)] can be expressed in IP as

Given ∣u_n〉 (the 2D field at z_n), the next step field ∣u_n+1〉 can be obtained from Eq. (16) by RK4. As in 1D case of Ref. [15], one can choose z₀ = z_n + h/2 to minimize the required number of FFTs, which leads to

where $|u_n^I+{\mathit{ck}}_i〉\equiv |u_n^I〉+c|k_i〉(i=1,2,3,c=\frac{h}{2},h$ ,

and d̂(h/2) can be obtained according to Eq. (27). Introduced in Eq. (27), N_h in Eq. (18) is the index of the last plate, whereas the unitary matrix m̂_i in Eq. (18) and its Hermitian (i.e., its self-adjoint matrix) m̂_i^† satisfy m̂_i^† = m̂_i^-1. To order all m̂_i in one direction, the indexes in M̂_L, [given in Eq. (18)] does not follow the ordering rule introduced in the Appendix. The minus sign in m̂_i(−) (i = 1,2…) is used to denote δ_i < 0, where ∣δ_i∣ is the length of the ith plate discussed in the Appendix. Thus we have m̂_i(−) =m̂_i^†. Obviously, when γ = 0, Eq. (17) yields $\mid u_{n+1}〉=\hat{d}\left(\frac{h}{2}\right){\hat{M}}_2\hat{d}\left(\frac{h}{2}\right){\hat{M}}_1\mid u_n〉$, which is consistent with Eq. (27) [or Eq. (22)and Eq. (24)] given in the Appendix.

4. RK4IP and other methods: numerical comparisons for various transmissions

In this work, the input optical field u_in(t) ≡ ∣u⃗_in(t)∣ is assumed to be a periodic repetition of N-bit (N=16) de Bruijn sequence, i.e., u_in(t) = ∑^∞_n=−∞d_B(t-nNT_b), where d_B(t) = ∑^N-1_i=0a_iP(t − iT_b) and T_b is the time interval of each bit. Here p(t) determines the elementary input pulse shape and a_i is the logic value of the ith bit. Within the time interval [0, T_b], the elementary forms of RZ and NRZ pulses (or Marks) are assumed to be $p\left(t\right)=\sqrt{2E_b/T_b}\cos \left[\frac{\pi }{2}{\cos }^2\left(\frac{\mathit{\pi t}}{T_b}\right)\right]$ and $p\left(t\right)=\sqrt{E_{\mathrm{b}}/T_{\mathrm{b}}}$, respectively, with E_b the optical energy per transmitted bit [22, 23]. Outside this time interval, p(t) is zero. Obviously, E_b/T_b is the mark power. To give the NRZ optical pulses slightly rounded edges, the input pulses are generated by passing through an input optical filter [2]. In this work, it is assumed to be fifth-order Bessel type with bandwidth B.

In the following numerical calculations, RK4IP is applied to various transmission fibers, where the parameter values can be quite different. The numerical results obtained from RK4IP are verified by comparing them with those using other approaches. For simplicity, we neglect the fiber loss, implying that there are no optical amplifiers in the system and, consequently, no amplifier (ASE) noise. In principle, the current RK4IP can be extended to the real systems where the fiber loss, amplifier gain, ASE from amplifiers, and the related parametric gain also need to be taken into account.

Only OOK format is considered in this section. Also, we insert an optical channel filter before photodetection, not only because it is used in a real optical system to suppress the ASE noise, but also it can suppress the spurious spectral lines discussed in Refs. [24–26] and others, from where one can see that such spurious components can be basically removed by limiting the simulation bandwidth and using varied step size. This optical filter can also suppress the linear effects (CD and PMD), which is small near the carrier frequency (ω = 0). The optical channel filter is Fabry-Pérot type [23]. The square-law-detected signal is electrically filtered by a fifth-order Bessel filter [23]. In this section, the bandwidths of optical (B_o) and electrical (B_e) filters are B_o = 6.9/T_b [27] and B_e = 0.8/T_b [3, 27, 28], unless otherwise noted. The electrically filtered photoelectric pulses are represented in units of W, since detected current corresponds to optical power [2].

 

Fig. 1. (a) and (b): Electrically filtered pulse as a function of time (T_b = 200ps), with CD=17 ps/nm·km and fiber birefringence being neglected. The nonlinear coefficient γ = 1.26/(W·km) is obtained using effective mode area A_eff = 80 μm². The input mark power E_b/T_b is 20 mW for L=100 km (a) and 2 mW for L=500 km (b). There is no visible difference between RK4IP solution [Eq. (17)] and SSFM solution for Eq. (19). (c): Optical power as a function of time (before the optical filter). Other parameters are same as those in (b). (d): Electrically filtered pulse as a function of time (no optical filter). Related parameters shown in (d) were given by Ref. [2] [p. 1740, Fig. 4(c)].

Download Full Size | PDF

In our computation, the approach of adaptive step size is used. As discussed in Ref. [12], given a step size h and the obtained solution at x_n (∣u(x_n)〉_ob), one can use RK4IP [Eq. (17)] to obtain ∣u(x_n + h/2)〉_f (the fine solution at x_n + h/2) and ∣u(x_n + h)〉_c (the coarse solution at x_n + h). Similarly, based on ∣u(x_n + h/2)〉_f, the second fine solution ∣u(x_n + h)〉_f can be obtained. Then the next step size can be adjusted, depending on the difference between ∣u(x_n + h)〉_f and ∣u(x_n + h)〉_c (cf. Ref. [12] for details). According to the local error method discussed in Ref. [12], the local error of RK4IP can be improved from O(h⁵) [15] to ~ O(h⁶) by using ∣u(x_n + h)〉_ob = [16∣u(x_n + h)〉_f − ∣u(x_n + h)〉_c]/15.

4.1. RK4IP and SSFM comparison: zero birefringence

Without birefringence, 2D CNLSE [Eq. (3)] can be reduced to 1D NLSE [Eq. (9)] with its split-step solution being (h = z_n+1 - z_n)

where Ũ is the Fourier transformation of ũ, while F^-1 denotes inverse Fourier transformation. In the case of zero birefringence, the SSFM result of Ref. [2] can be reduced to Eq. (19).

Obviously, Eq. (19) should yield the same result as the RK4IP solution [Eq. (17), without birefringence]. This can be numerically confirmed by considering a NRZ-OOK pulse train launched into the fiber with CD=17 ps/(nm·km) and A_eff=80μm² [γ=1.26 (W·km)^-1]. The bitrate is 5Gb/s and the bandwidths of the electrical filter and the input optical filter are B_e = 4GHz and B = 1.75B_e. The input mark power is 20 mW for L=100 km (a) and 2 mW for L=500 km (b). As shown in Figs. 1(a) and 1(b), the electrically filtered currents obtained from RK4IP (without birefringence) agrees very well with those from SSFM [Eq. (19)].

To show that the good agreement between the two filtered currents is due to the good agreement between the RK4IP and SSFM solutions, in Fig. 1(c), we plot the optical power distributions obtained from these two approaches for the case of Fig. 1(b). Obviously, similar optical distributions lead to similar electrically filtered currents.

Our RK4IP and SSFM results not only agree well with each other, but also agree with the published results. As an example, one can simulate the 5 (80+20)-km spans system discussed in Ref. [2], where (Fig. 2, p.1740) the detailed dispersion map was plotted. Fig. 1(d) shows our numerical results with B_e = 1.0/T_b (T_b = 200 ps). Like Ref. [2], no optical filter is used here. Related parameters shown in Fig. 1(d) are given by Ref. [2]. Indeed, the two curves plotted in Fig. 1(d) are almost the same as those shown in the Fig. 4(d) of Ref. [2]. [Decreasing γ from 2.02/W·km to 8/9γ=1.80/W·km will reduced the current shown in Fig. 1(d) by ~ 1%.]

4.2. Step size of RK4IP

Here, the step size means that, within given computational error tolerance (< 1% of the pulse peak), the maximum allowable step size at the end of the fiber.

Table 1. Step size Δz using RK4IP and SSFM of Ref. [2] obtained with given L_N and fiber length L. The dispersion length L_D ~ 240 km and the computational accuracy < 1% of the pulse peak.

View Table

Table 1 shows the step sizes using RK4IP for the fibers of Figs. 1(a) and 1(b) with random birefringence [Λ_beat =50m, L_corr=10m, and D_PMD=1.0ps/(km)^1/2, yielding average DGD=D_PMD√L ≈ 10 ps for L = 100 km and 22 ps for L = 500 km], compared with the step size using RK4IP (without birefringence) and the step size using Eq. (19), which is the SSFM of Ref. [2] in the case of zero birefringence. Notice that the fiber dispersion length L_D ~ 1/(β_ωωΔf²) (Δf is the signal bandwidth) for the two cases of Fig. 1 is ~ 240 km, while the nonlinear length L_N ~ 1/(γ^P) is ~ 40 km for the case of Fig. 1(a) and 400 km for Fig. 1(b).

In each case of zero birefringence, the RK4IP step size is around the smaller one of L_D and L_N, or, about 6~10 times of the step size of SSFM of Ref. [2], which is around 10% ~ 15% of L_D and/or L_N.

 

Fig. 2. The filtered photoelectric current (W) vs time (T_b = 100 ps) in a NRZ-OOK system consisting of 5 100-km spans of transmission fiber [CD₁₀₀=17ps/(nm·km), A_eff =80 μm²]. Each span is followed by a 13-km dispersion compensation fiber [CD₁₃=-120ps/(nm·km), A_eff=30 μm²]. Before the 1st 100-km fiber, there is a precompensation of -120×6 ps/nm, whereas the 5th 100-km fiber is followed by the compensation of -120×5 ps/nm. The input mark power is 10 mW. Birefringence parameters are D_PMD=2.0 ps/(km)^1/2, L_corr = 10 m, and Λ_beat = 50 m. There is very little difference between the RK4IP solution [Eq. (17)] and the M-PMD approx [Eq. (8)].

Download Full Size | PDF

Taking into account random birefringence, the RK4IP step size is decreased from 30km to 7 km for the case of L_N ~ 40 km (L=100km) and from 250 km to 45km for L_N ~ 400 km (L=500km). In general, the stronger the interaction between the linear and nonlinear parts in CNLSE, the more impact of L_corr and Λ_beat on the step size, assuming that L_D and L_N are much larger than the former. (For the same reason, the step size of approach [2] also needs to be decreased correspondingly.)

4.3. RK4IP solutions and Manakov-PMD approximations with random birefringence

Figure 2 shows good agreement between RK4IP result and M-PMD approx for a NRZ-OOK system with 5 100-km spans of transmission fiber and 5 13-km spans of dispersion compensation fiber. The bit-rate is 10Gb/s. To get the received photoelectric pulses shown in Fig. 2, a precompensation fiber (6 km) is inserted before the first 100-km fiber, whereas the last compensation fiber is 5 km long. The birefringence related parameters are Λ_beat = 50 m, L_corr=10 m, and D_PMD=2.0 ps/(km)^1/2, which means average DGD is around $2.0\sqrt{5\left(100+13\right)-2}\approx 47\mathrm{ps}$ with the birefringence direction being randomly changed every 10 m. As in subsection 4.1, the input optical filter with bandwidth of B = 1.75B_e is used to generate rounded edges for NRZ signal. Other paratemeters given in the caption of Fig. 2 are based on the discussion of Ref. [3]. As plotted in Fig. 2, there is no significant difference between the RK4IP result [Eq. (17)] and the M-PMD approx [Eq. (8)]. Taking L_corr=100 m yields almost the same results. To confirm that such agreement is not because of the weak nonlinearity, one can artificially increase the nonlinear coefficient in M-PMD approx [Eq. (8)] by 12.5%. Indeed, the M-PMD approx thus obtained $[\frac{8}{9}\gamma =3.36/(\mathrm{W}·\mathrm{km})]$ differs from the original M-PMD approx $[\frac{8}{9}\gamma =2.99/(\mathrm{W}·\mathrm{km})]$ by more than 20%. (not plotted in Fig. 2). The good agreement shown in Fig. 2 verifies both the the RK4IP result and M-PMD approx for such 5 100-km spans of NRZ transmission.

 

Fig. 3. The filtered photoelectric current (W) vs time (T_b = 100 ps) in a RZ-OOK system using 50 100-km eLEAF fibers [average CD₁₀₀=4.5ps/(nm·km), A_eff =72 μm²], precom-pensation of -(3 × 39) ps/nm before the first 100-km fiber, -(11 × 39) ps/nm compensation per span, and -(6 × 39) ps/nm compensation after the 50th 100-km fiber. The input mark power is 2 mW and the PMD coefficient is D_PMD=0.40 ps/(km)^1/2.

Download Full Size | PDF

Figure 3 shows the RK4IP result and M-PMD approx for a RZ-OOK system using 50 100-km spans of eLEAF fiber [27], with effective mode area being A_eff=72 μm² and CD in the range of 3 ~ 6 ps/(km·nm). To get the received photoelectric pulse shown in Fig. 3, each 100-km fiber with CD₁₀₀=4.5ps/(km·nm) is followed by -429 ps/nm compensation. Also, precompensation of -117 ps/nm before the first 100-km fiber and last compensation of -234 ps/nm after the 50th 100-km fiber are used. The birefringence related parameters are L_corr = 100 m and Λ_beat = 50 m. As shown in Fig. 3, the agreement between the RK4IP result and M-PMD approx is not as good as that of Fig. 2. Decreasing L_corr to 10 m, we get a similar agreement. Considering that, due to the linear-nonlinear interaction accumulated in such 50 100-km fibers, any small change in the related parameters will lead to significant change in received pulse shape, the M-PMD approx (dot-dashed) in Fig. 3 can be viewed as a reasonably good approximation.

In Fig. 4, the thin dotted curve (L_corr=100 m) is used to show that, when the CD parameter of eLEAF fiber is increased to its worst case [6.0ps/(km·nm)], the discrepancy between RK4IP result and M-PMD approx becomes larger than that shown in Fig. 3, due to the increased interaction between linear and nonlinear effects in the CNLSE. Other related parameters for this dotted curve are same as those used in Fig. 3, except that each span is followed by a -585 ps/nm dispersion compensation and the 50th 100-km fiber is followed by a -390 ps/nm compensation. To improve the mixing on the Poincaré sphere, one can reduce the fiber correlation length from L_corr = 100 m to L_corr = 10 m. The dashed curve thus obtained is very close to the filtered current using M-PMD approx (dot-dashed). The M-PMD approx is not sensitive to the value of L_corr.

 

Fig. 4. Same as Fig. 3, except that CD₁₀₀=6.0ps/(nm·km) and that each span is followed by a -(15 × 39) ps/nm dispersion compensation and a -(10 × 39) ps/nm compensation after the 50th 100-km fiber.

Download Full Size | PDF

The difference between the RK4IP result and the M-PMD approx shown in Fig. 4 is further discussed in Subsection 5.3.

5. Further discussions based on the obtained results

5.1. Key points for RK4IP method being able to capture the whole randomness existing in large step size

Why RK4IP method can capture the whole randomness existing in large step size? One key point is the interaction picture (IP). For each step, the role of the interaction operator [actually the Î^† in Eq. (16), caused by the linear part in CNLSE] is to rotate the reference system, while ∣u^I(z_n,t)〉 and N̂₂^I are the 2D optical field and the nonlinear effect observed in this rotating system. When one wants to return from the rotating system to the laboratory coordinate system, the inverse transformation [actually the Î in Eq. (16) and given by Eq. (27)] is needed. If there is no random birefringence, the system rotates with z constantly (for a given frequency). It will be disturbed when random birefringence exists. With the help of Eq. (26) in the Appendix, the random birefringence in Î can be treated efficiently. The final result obtained from algebra recursion [Eq. (26)] is the lumped birefringence (for a given frequency). Note that no matter how large the step size is and how randomly the birefringence varies, the CNLSE in IP is correct. No approximation is introduced. (This is different from the various SSFM methods, based on which the larger the step-size the less the accuracy.) The only reason to limit the step size of RK4IP is the RK algorithm, which has local error of O(h⁶), provided the local error method of Ref. [12] is used. A RK4IP solution with too large step size will deviate from the true value, like any numerical method used to solve a given differential equation.

Different from RK4IP, the SSFM of Ref. [2] accurately solved the linear and nonlinear parts in Eq. (8) (cf. Ref. [2]). However it combined them using split-step approximation. Thus, the step size of this approach should be small enough.

5.2. Computational time and the complexity of RK4IP

To reduce the complexity of RK4IP, an efficient way to calculate the interaction operator is not enough. One still needs to concern the computation of the RK4 part (in IP).

In the case of no birefringence, the computational times of RK4IP and SSFM for solving the NLSE are dominated by the number of FFT per step. It is true that SSFM of [2] needs 2 FFTs per step, while RK4IP needs 8 FFTs per step (see the related discussion in 3.1). However, the step size of RK4IP is much larger than the SSFM one. Taking the case with L=100 km in Table 1 as an example, for the simulation length of 30 km, the RK4IP needs 8 FFTs, while the SSFM needs 2(FFTs/step)×10(steps)=20 FFTs. Thus, for a given fiber length, the required FFTs using RK4IP is around 40% of the FFTs required by SSFM of [2]. In fact, for such 100-km long fiber, the simulation time using each method is less than 0.5 minute. In the cases with random birefringence, it takes RK4IP about a few minutes to get each curve in Figs. 1–4.

Again, we indicate that a good algorithm with large step is better than a simple algorithm with small step size, even if, sometimes, the former needs longer time than the later. This is because a simple algorithm with too small step size (e.g., the “full split-step” illustrated in p. 1740 of Ref. [2]) not only needs too much computation, but also may lead to a questionable result due to the computational error accumulated.

5.3. Accuracies of RK4IP, SSFM, and M-PMD approx

1) RK4IP and SSFM

In principle, there should be no reason to see any discrepancy between RK4IP and SSFM results, provided the numerical error of each approach is negligible. Therefore, even in the cases of strong linear-nonlinear interaction [e.g., high bit-rate (meaning strong CD and PMD effects) and/or high optical power] the results obtained from these two approaches still should be the same.

However, due to the numerical error per step, there is a step-size range for each approach, out of which one cannot get the result with good enough accuracy. The left-side of such step-size range (lower boundary) is determined by the computational error per step (e.g., caused by the numerical error in FFT). A step size smaller than this range does not mean any improvement in accuracy. The right-side of such range (upper boundary) is determined by the strength of linear and nonlinear effects. The stronger the linear and nonlinear effects, the smaller the value of the upper boundary. For a given computational error per step, this step-size range will be smaller and smaller when the linear-nonlinear interaction becomes stronger and stronger. Thus, RK4IP should be more tolerant to the linear-nonlinear interaction than SSFM, because its step size is larger than SSFM’s, as illustrated in Table 1.

Taking the 500-km SMF in Fig. 1(b) as an example, one can increase the input power from 2 mW to 40 mW to get strong enough linear and nonlinear effects. In fact, there is no visible difference between the RK4IP current and the SSFM current. The obtained RK4IP current keeps stable when we decrease the step size from a few kilometers to 100 m, while the obtained SSFM current is stable with the step size around 50~30 m. Decreasing the optical power from 40 mW, the step-size range for the SSFM becomes larger and larger.

2) RK4IP and M-PMD approx

Considering that the currents shown in Figs. 2–4 are filtered by the optical and electrical filters, which may suppress the discrepancy between RK4IP result and M-PMD approx, in this subsection, we consider the optical distribution before the optical filter.

Figure 5 shows the optical power distributions obtained using RK4IP and M-PMD approx with L_corr = 10 m (a) and L_corr = 100 m (b). The curves in Fig. 5 indicates that, in this case, the value of the L_corr does not significantly affect the agreement between these two approaches.

 

Fig. 5. The optical power (W) vs time (before the optical filter, T_b = 100 ps), with L_corr=10 m (a) and L_corr=100 m (b). Other related parameters are same as those in Fig. 2. Here, the agreement between the two approaches is not very sensitive to the value of L_corr.

Download Full Size | PDF

 

Fig. 6. The optical power (W) vs time (before the optical filter, T_b = 100 ps), with L_corr=10 m (a) and L_corr=100 m (b). Other related parameters are the same as those in Fig. 4. The agreement between the two approaches is affected by the value of L_corr because of the strong linear-nonlinear interaction.

Download Full Size | PDF

However, Fig. 6 shows a different story. When L_corr = 100 m, the RK4IP does not agree well with the M-PMD approx, e.g., the optical power near 9.0T_b in Fig. 6(b). It is such discrepancy in optical powers that results in the discrepancy in filtered currents (the dotted and dot-dashed curves) shown in Fig. 4. When one improves the mixing of the polarization state on the Poincaré sphere by decreasing the fiber correlation length to L_corr = 10 m, the good agreement shown in Fig. 6(a) can be obtained. This leads to the good agreement between the dashed and dot-dashed curves in Fig. 4.

Results given in Figs. 5–6 means that the mixing condition required by M-PMD approx is also related with the strength of the linear-nonlinear interaction. The stronger the linear-nonlinear interaction the harder to meet the mixing requirement for M-PMD approx (e.g, the shorter the L_corr should be).

6. Summary

Based on the RK4IP method for 1D NLSE [15] as well as the analytical solution for the linear terms in CNLSE [2], RK4IP method for 2D CNLSE is presented in this work. Without rotating the coordinate system for each step, our approach can be applied to a fiber with general form of birefringence. Moreover, the Kerr nonlinearity in the CNLSE is treated without approximation. Thus, there is no need for RK4IP to concern the mixing condition required by the M-PMD approx. In the cases of strong linear-nonlinear interaction (e.g., bit-rate > 40 Gb/s and/or the nonlinear phase of the fiber (ϕ_NL = L/L_N ~ LγP > 5π), the fully mixing requirement for M-PMD approx may not be well satisfied.

Since there is no split-step approximation for each step and the local error method of Ref. [12] is used, for normal fibers with random birefringence, the step size using RK4IP can be orders of magnitude larger than Λ_beat and L_corr, depending on the intensity of the linear-nonlinear interaction. Without birefringence effect, the RK4IP step size can be increased to the same order of magnitude as L_D and/or L_N, or, around 6 ~ 10 times the step size using the approach of Ref. [2]. Our numerical RK4IP results agree well with those obtained from the M-PMD approx, provided the polarization state can be mixed enough on the Poincaré sphere.