A fourth-order Runge-Kutta in the interaction picture (RK4IP) method is presented for solving the coupled nonlinear Schrödinger equation (CNLSE) that governs the light propagation in optical fibers with randomly varying birefringence. The computational error of RK4IP is caused by the fourth-order Runge-Kutta algorithm, better than the split-step approximation limited by the step size. As a result, the step size of RK4IP can have the same order of magnitude as the dispersion length and/or the nonlinear length of the fiber, provided the birefringence effect is small. For communication fibers with random birefringence, the step size of RK4IP can be orders of magnitude larger than the correlation length and the beating length of the fibers, depending on the interaction between linear and nonlinear effects. Our approach can be applied to the fibers having the general form of local birefringence and treat the Kerr nonlinearity without approximation. Our RK4IP results agree well with those obtained from Manakov-PMD approximation, provided the polarization state can be mixed enough on the Poincaré sphere.
©2010 Optical Society of America
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  and further studied in many publications.
For a non-soliton system, the CNLSE is solved with the step size determined by dispersion length LD, nonlinear length LN, as well as the birefringence related parameters such as fiber correlation length (Lcorr) (the length within which birefringence axes are randomly reoriented), beating length Λbeat (the length determined by birefringence strength), and the PMD parameter DPMD [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 LD and LN. 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 Lcorr (0.3~300 m) are much smaller than LD and LN (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 Lcorr 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. , 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 Lcorr. (Detailed value of its step size depends on the interaction between linear and nonlinear effects.) Assumptions required by Ref.  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 , 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.  is that it needs split-step Fourier method (SSFM), which is based on the first order approximation of Baker-Hausdorff formula 
In Eq. (1), [A, B] ≡ AB − BA. Because of this, even in the case of zero birefringence, the step size using approach  is restricted to 10% ~ 15% of LD and/or LN, 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) , the RK4 in IP method (RK4IP) was recently proposed to study the supercontinuum generation in optical fibers (for 1D field cases) . 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.  cannot be applied directly to the cases with random birefringence, because in Ref. , 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.  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 , the local error of our RK4IP can be improved to ~ O(h6), rather than ~ O(h5) of 1D RK4IP . As results, for the communication fibers with random birefringence, the step size using RK4IP can be orders of magnitude larger than Lcorr and Λbeat. Without birefringence, the step size of RK4IP can be the same order of magnitude as LD and/or LN, or, 6 ~ 10 times the step size using the approach of Ref. .
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 . Namely, given a field with its x − y components being ux and uy, it can be denoted as ∣u〉 ≡ (ux,uy)T [or 〈u∣ ≡ (ux*,uy*)]. Thus, the CNLSE discussed in [1–3] can be written in the following form with retarded time t = tlab − βωz:
where ∣u〉z ≡ ∂ ∣u〉/∂z and Δβ (Δβω ≡ ∂Δβ/∂ω) is related to Λbeat (Lcorr) with , respectively (cf. Eq. (28) of  or P. 6 of  and others). Here DPMD is the PMD coefficient . When L >> Lcorr, the average DGD of a L-km-long fiber can be obtained using DGDavg = DPMD√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〉 ≡ (ux*uy2,ux2uy*)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., ).
In this work, the three components of the Pauli spin matrices are denoted as
with βi ≡ 〈β0(z)∣σi∣β0(z)〉 (i = 1,2,3) and ∣β0(z)〉 [β⃗0(z)] being, respectively, the unit vector representing the fiber birefringence in 2D Jones (3D Stokes) space . For the fiber with circular birefringence (e.g., spun optical fiber), β2 ≠ 0 in Eq. (5).
Parameter βωω = ∂2β/∂ω2 corresponds to the CD parameter at wavelength λ by βωω=-CD(λ)λ2/(2πc) (ps/nm·km) with c=3 × 108m/s.
we have T†σ13T = σ13, 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 n2, effective mode area Aeff, and wavenumber k = 2π/λ by γ = n2k/Aeff. Typical values of λ = 1550 nm and n2 = 2.6 × 10-20m2/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 .
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
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. , by substituting the unitary transformations ∣u〉 = RT∣ 〉 into Eq. (8), with T being given by Eq. (6)and , yielding and [σ2αz + σ3Δβ]T + jTz = 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
Equation (11) is equivalent to
N̂I in Eq. (13) is the nonlinear operator represented in the IP. Given u(zn,t), the next step field u(zn+1,t) (zn+1 = zn + h) governed by Eq. (12) can be obtained using RK4, with the local error of O(h5). According to RK4 , we have uI(zn+1,t) = uI(zn,t) + h[k1 + 2(k2 + k3) + k4]/6 with k1 = f[zn,uI(zn,t)], k2 = f[zn + h/2,uI(zn,t) + k1h/2], k3 = f[zn + h/2,uI(zn,t) + k2h/2], and k4 = f[zn + h,uI(zn,t) + k3h]. Setting z0 = zn, one needs 16 FFTs to get u(zn+1,t) from u(zn,t): 12 FFTs for k2, k3, k4, and 4 FFTs for u(zn,t) → uI(zn,t) and uI(zn+1,t) → u(zn+1,t).
Choosing z0 = zn + h/2, one only needs 8 FFTs, yielding 
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,z0) ≡ exp[j(z - z0)D̂] = exp[-j(z - z0)βωω(∂2/∂t2)/2] in Eq. (13) needs to be modified correspondly, which is the key point to extend the RK4IP of  to 2D case.
3.2. RK4IP solution for 2D CNLSE with random birefringence
Equation (3) can also be viewed as
Given ∣un〉 (the 2D field at zn), the next step field ∣un+1〉 can be obtained from Eq. (16) by RK4. As in 1D case of Ref. , one can choose z0 = zn + h/2 to minimize the required number of FFTs, which leads to
and d̂(h/2) can be obtained according to Eq. (27). Introduced in Eq. (27), Nh 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 , 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 uin(t) ≡ ∣u⃗in(t)∣ is assumed to be a periodic repetition of N-bit (N=16) de Bruijn sequence, i.e., uin(t) = ∑∞n=−∞dB(t-nNTb), where dB(t) = ∑N-1i=0aiP(t − iTb) and Tb is the time interval of each bit. Here p(t) determines the elementary input pulse shape and ai is the logic value of the ith bit. Within the time interval [0, Tb], the elementary forms of RZ and NRZ pulses (or Marks) are assumed to be and , respectively, with Eb the optical energy per transmitted bit [22, 23]. Outside this time interval, p(t) is zero. Obviously, Eb/Tb 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 . 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 . The square-law-detected signal is electrically filtered by a fifth-order Bessel filter . In this section, the bandwidths of optical (Bo) and electrical (Be) filters are Bo = 6.9/Tb  and Be = 0.8/Tb [3, 27, 28], unless otherwise noted. The electrically filtered photoelectric pulses are represented in units of W, since detected current corresponds to optical power .
In our computation, the approach of adaptive step size is used. As discussed in Ref. , given a step size h and the obtained solution at xn (∣u(xn)〉ob), one can use RK4IP [Eq. (17)] to obtain ∣u(xn + h/2)〉f (the fine solution at xn + h/2) and ∣u(xn + h)〉c (the coarse solution at xn + h). Similarly, based on ∣u(xn + h/2)〉f, the second fine solution ∣u(xn + h)〉f can be obtained. Then the next step size can be adjusted, depending on the difference between ∣u(xn + h)〉f and ∣u(xn + h)〉c (cf. Ref.  for details). According to the local error method discussed in Ref. , the local error of RK4IP can be improved from O(h5)  to ~ O(h6) by using ∣u(xn + h)〉ob = [16∣u(xn + h)〉f − ∣u(xn + h)〉c]/15.
4.1. RK4IP and SSFM comparison: zero birefringence
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 Aeff=80μm2 [γ=1.26 (W·km)-1]. The bitrate is 5Gb/s and the bandwidths of the electrical filter and the input optical filter are Be = 4GHz and B = 1.75Be. 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. , where (Fig. 2, p.1740) the detailed dispersion map was plotted. Fig. 1(d) shows our numerical results with Be = 1.0/Tb (Tb = 200 ps). Like Ref. , no optical filter is used here. Related parameters shown in Fig. 1(d) are given by Ref. . Indeed, the two curves plotted in Fig. 1(d) are almost the same as those shown in the Fig. 4(d) of Ref. . [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 shows the step sizes using RK4IP for the fibers of Figs. 1(a) and 1(b) with random birefringence [Λbeat =50m, Lcorr=10m, and DPMD=1.0ps/(km)1/2, yielding average DGD=DPMD√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.  in the case of zero birefringence. Notice that the fiber dispersion length LD ~ 1/(βωωΔf2) (Δf is the signal bandwidth) for the two cases of Fig. 1 is ~ 240 km, while the nonlinear length LN ~ 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 LD and LN, or, about 6~10 times of the step size of SSFM of Ref. , which is around 10% ~ 15% of LD and/or LN.
Taking into account random birefringence, the RK4IP step size is decreased from 30km to 7 km for the case of LN ~ 40 km (L=100km) and from 250 km to 45km for LN ~ 400 km (L=500km). In general, the stronger the interaction between the linear and nonlinear parts in CNLSE, the more impact of Lcorr and Λbeat on the step size, assuming that LD and LN are much larger than the former. (For the same reason, the step size of approach  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, Lcorr=10 m, and DPMD=2.0 ps/(km)1/2, which means average DGD is around with the birefringence direction being randomly changed every 10 m. As in subsection 4.1, the input optical filter with bandwidth of B = 1.75Be 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. . As plotted in Fig. 2, there is no significant difference between the RK4IP result [Eq. (17)] and the M-PMD approx [Eq. (8)]. Taking Lcorr=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 differs from the original M-PMD approx 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.
Figure 3 shows the RK4IP result and M-PMD approx for a RZ-OOK system using 50 100-km spans of eLEAF fiber , with effective mode area being Aeff=72 μm2 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 CD100=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 Lcorr = 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 Lcorr 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 (Lcorr=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 Lcorr = 100 m to Lcorr = 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 Lcorr.
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 ∣uI(zn,t)〉 and N̂2I 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(h6), provided the local error method of Ref.  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.  accurately solved the linear and nonlinear parts in Eq. (8) (cf. Ref. ). 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  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 . 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. ) 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 Lcorr = 10 m (a) and Lcorr = 100 m (b). The curves in Fig. 5 indicates that, in this case, the value of the Lcorr does not significantly affect the agreement between these two approaches.
However, Fig. 6 shows a different story. When Lcorr = 100 m, the RK4IP does not agree well with the M-PMD approx, e.g., the optical power near 9.0Tb 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 Lcorr = 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 Lcorr should be).
Based on the RK4IP method for 1D NLSE  as well as the analytical solution for the linear terms in CNLSE , 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/LN ~ 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.  is used, for normal fibers with random birefringence, the step size using RK4IP can be orders of magnitude larger than Λbeat and Lcorr, 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 LD and/or LN, or, around 6 ~ 10 times the step size using the approach of Ref. . 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.
Appendix: Optical field in the fiber without nonlinearity
For a system with γ = 0, Eq. (3) can be reduced to
In frequency domain, Eq. (20) can be simplified as
where ∣U(z,ω)〉 = ∫∣u(z, t)∣e-jωt dt. (Note that in Ref. , the Fourier transformation was defined to be ∣U(z,ω)〉 = ∫∣u(z, t)〉ejωt dt.)
where δi = zi - zi-1. When zi > zi-1, δi is the length of the ith plate. The birefringence direction in the ith plate is β⃗0 (zi). In Eq. (24) and Eq. (25), Nh is the number of the plates between z0 and z = z0 + h and the index Nh denotes the Nhth plate. The expression of M(z, ι) given by Eq. (24) can also be extended to the case of z - z0 < 0, provided to notice that 1) for each plate, δi < 0; 2) mi (i = 1, …, Nh) are ordered from z0 to z, i.e., m1 is the first left plate from z0.
In general, Eq. (24) and Eq. (25) can be applied to a h-long fiber that is viewed as the concatenation of a number of thin birefringent plates with its ith plate (length δi) having a constant Δβi [meaning constant (Δβω)i = ∂Δβi/∂ω], provided the distributions of Δβi and (Δβω)i (i = 1,…, Nh) in this h-long fiber are known. Like Ref. , in 4.2 and 4.3 we simply assume that, for a given fiber, its phase birefringence strength (Δβ) and group birefringence strength (Δβω) are constants in z direction. Thus η in Eq. (25) is also a constant. The birefringent directions [β⃗0(zi)] are randomly oriented.
where I is the unitary matrix.
In time domain, the transform matrices in Eq. (22) have the form of (i = 1, …,Nh)
With the help of Eq. (26), the interaction operator in Eq. (27) [or Eq. (16)] can be treated in a very efficient way. Within each step, although there are a lot of randomly oriented plates (e.g., 700~4500 or more plates in each step for the cases of Table 1), one can treat them using the algebra recursion given by Eq. (26), rather than using the “full split-step calculation” method illustrated in Ref.  (p. 1740), where the CNLSE needs to be solved many times for each randomly rotated plate.
Acknowledgment The authors acknowledge the financial support from Canadian funding agencies: NSERC. They thank the anonymous reviewers for their insightful comments and suggestions.
References and links
1. P. K. A. Wai and C. R. Menyuk, “Polarization mode dispersion, decorrelation, and diffusion in optical fibers with randomly varying birefringence,” J. Lightwave Technol. 14, 148–157 (1996). [CrossRef]
2. D. Marcuse, C. R. Manyuk, and P. K. A. Wai, “Application of the Manakov-PMD equation to studies of signal propagation in optical fibers with randomly varying birefringence,” J. Lightwave Technol. 15, 1735–1746 (1997). [CrossRef]
3. C. R. Menyuk and B. S. Marks, “Interaction of polarization mode dispersion and nonlinearity in optical fiber transmission systems,” J. Lightwave Technol. 24, 2806–2826 (2006). [CrossRef]
4. C. Xie, M. Karlsson, P. A. Andrekson, H. Sunnerud, and J. Li, “Influences of polarization-mode dispersion on soliton transmission systems,” J. Sel. Top. Quantum Electron. 8, 575–590 (2002). [CrossRef]
5. E. Alperovich, A. Mecozzi, and M. Shtaif, “PMD penalties in long nonsoliton systems and the effect of inline filtering,” IEEE Photon. Technol. Lett. 18, 1179–1181 (2006). [CrossRef]
6. C. R. Menyuk, “Nonlinear pulse propagation in birefringent optical fibers,” IEEE J. Quantum Electron. 23, 174–176 (1987). [CrossRef]
7. M. Karlsson, “Modulational instability in lossy optical fibers,” J. Opt. Soc. Am. B 12, 2071–2077 (1995).
8. H. Ghafouri-Shiraz, P. Shum, and M. Nagata, “A novel method for analysis of soliton propagation in optical fibers,” IEEE J. Quantum Electron. 31, 190–200 (1995). [CrossRef]
9. A. V. T. Cartaxo, “Small-signal analysis for nonlinear and dispersive optical fibres, and its application to design of dispersion supported transmission systems with optical dispersion compensation,” Proc. Inst. Elect. Eng., Optoelectron. 146(5), 213–222 (1999). [CrossRef]
10. M. J. Ablowitz and T. Hirooka, “Nonlinear effects in quasi-linear dispersion-managed pulse transmission,” IEEE Photon. Technol. Lett. 13, 1082–1084 (2001). [CrossRef]
11. A. Vannucci, P. Serena, and A. Bononi, “The RP method: a new tool for the iterative solution of the nonlinear Schrodinger equation,” J. Lightwave Technol. 20, 1102–1112 (2002). [CrossRef]
12. O. V. Sinkin, R. Holzlohner, J. Zweck, and C. R. Menyuk, “Optimization of the split-step Fourier method in modeling optical-fiber communications systems,” J. Lightwave Technol. 21, 61–68 (2003). [CrossRef]
13. E. Ciaramella and E. Forestieri, “Analytical approximation of nonlinear distortions,” IEEE Photon. Technol. Lett. 17, 91–93 (2005). [CrossRef]
14. M. Secondini, E. Forestieri, and C. R. Menyuk, “A combined regular-logarithmic perturbation method for signal-noise interaction in amplified optical systems,” J. Lightwave Technol. 27, 3358–3369 (2009). [CrossRef]
15. J. Hult, “A fourth-order Runge-Kutta in the interaction picture method for simulating supercontinuum generation in optical fibers,” J. Lightwave Technol. 25, 3770–3775 (2007). [CrossRef]
16. G. H. Weiss and A. A. Maradudin, “The Baker-Hausdorff formula and a problem in crystal physics,” J. Math. Phys. 3, 771–777 (1962). [CrossRef]
17. J. Butcher, Numerical Methods for Ordinary Differential Equations, (Wiley, 2003). [CrossRef]
18. J. N. Damask, Polarization Optics in Telecommunications, (Springer2005).
19. A. Galtarossa and C. R. Menyuk, Polarization Mode Dispersion, (Springer2005). [CrossRef]
20. K. W. Chow, K. K. Y. Wong, and K. Lam, “Modulation instabilities in a system of four coupled, nonlinear Schrodinger equations,” Phys. Lett. A 372, 4596–4600 (2008). [CrossRef]
21. S. G. Evangelides Jr., L. F. Mollenauer, J. P. Gordon, and N. S. Bergano, “Polarization multiplexing with solitons,” J. Lightwave Technol. 10, 28–35 (1992). [CrossRef]
22. E. Forestieri, “Evaluating the error probability in lightwave systems with chromatic dispersion, arbitrary pulse shape and pre-and postdetection filtering,” J. Lightwave Technol. 18, 1493–1503 (2000). [CrossRef]
23. J. Wang and J. M. Kahn, “Impact of chromatic and polarization-mode dispersions on DPSK systems using inter-ferometric demodulation and direct detection,” J. Lightwave Technol. 22, 362–371 (2004). [CrossRef]
24. F. Matera, A. Mecozzi, M. Romagnoli, and M. Settembre, “Sideband instability induced by periodic power variation in long-distance fiber links,” Opt. Lett. 18, 1499–1501 (1993). [CrossRef] [PubMed]
25. A. Carena, V. Curri, R. Gaudino, P. Poggiolini, and S. Benedetto, “A time-domain optical transmission system simulation package accounting for nonlinear and polarization-related effects in fiber,” IEEE J. Select. Areas Commun. 15, 751–765 (1997). [CrossRef]
26. G. Bosco, A. Carena, V. Curri, R. Gaudino, P. Poggiolini, and S. Benedetto, “A novel analytical approach to the evaluation of the impact of fiber Parametric Gain on the Bit Error Rate,” IEEE Trans. Commun. 49, 2154–2163 (2001). [CrossRef]
27. J. M. Wiesenfeld, L. D. Garrett, M. Shtaif, M. H. Eiselt, and R. W. Tkach, “Effects of DGE bandwidth on nonlinear ULH systems,” in Proc. Optical Fiber Communications Conf., Anaheim, CA, 2005, paper OWA2.
28. E. Ibragimov, C. Menyuk, and W. Kath, “PMD-induced reduction of nonlinear penalties in terrestrial optical fiber transmission,” in Proc. OFC, 2000, pp. 195–197, Paper WL3.