A fourth-order Runge-Kutta in the interaction picture method for coupled nonlinear Schrodinger equation

Abstract

A fourth-order Runge-Kutta in the interaction picture (RK4IP) method is presented for solving the coupled nonlinear Schr $¨ o$ 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. For the systems with realistic parameters, the RK4IP results are consistent with those using Manakov-PMD approximation [1, 2, 3]. However, increased interaction between the linear and nonlinear terms in CNLSE leads to increased discrepancy between RK4IP and Manakov-PMD approximation.

A fourth-order Runge-Kutta in the interaction picture method for coupled nonlinear Schr $¨ o$ dinger equation

Zhongxi Zhang, Liang Chen, Xiaoyi Bao

Department of Physics, university of Ottawa, 150 Louis Pasteur, Ottawa, K1N 6N5, Ontario, Canada

zzhan3@uottawa.ca

© 2019 Optical Society of America OCIS codes: (060.2330) Fiber optics communications; (260.2030) Dispersion; (190.3270) Kerr effect; (060.4370) Nonlinear optics, fiber; (060.5530) Pulse propagation and temporal solitons

References and links

[1] P.K.A. Wai, C.R. Menyak, “Polarization mode dispersion, decorrelation, and diffusion in optical fibers with randomly varying birefringence,” J. Lightwave Technol. 14, 148-157 (1996).
[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).
[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).
[4] C. Xie, M. Karlsson, P.A. Andrekson, H. Sunnerud, J. Li, “Influences of polarization-mode dispersion on soliton transmission systems,” J. Sel. Top. Quantum Electron. 8, 575-590 (2002).
[5] E. Alperovich, A. Mecozzi, M. Shtaif, “PMD penalties in long nonsoliton systems and the effect of inline filtering,” IEEE Photon. Technol. Lett. 18, 1179-1181 (2006).
[6] C. R. Menyuk, “Nonlinear pulse propagation in birefringent optical fibers,” IEEE J. Quantum Electron. QE-23, 174-176 (1987).
[7] M. Karlsson, “Modulational instability in lossy optical fibers,” J. Opt. Soc. Amer. 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. QE-31, 190-200 (1995).
[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, no.5, pp. 213-222 (1999).
[10] M. J. Ablowitz, T. Hirooka, “Nonlinear effects in quasi-linear dispersion-managed pulse transmission,” IEEE Photon. Technol. Lett. 13 , 1082-1084 (2001).
[11] A. Vannucci, P. Serena, A. Bononi, “The RP method: a new tool for the iterative solution of the nonlinear Schrodinger equation,” J. Lightwave Technol. 20, 1102- 1112 (2002).
[12] O.V. Sinkin, R. Holzlohner, J. Zweck, C.R. Menyuk, “Optimization of the split-step Fourier method in modeling optical-fiber communications systems,” J. Lightwave Technol. 21, 61-68 (2003).
[13] E. Ciaramella and E. Forestieri, “Analytical approximation of nonlinear distortions,” IEEE Photon. Technol. Lett. 17, 91-93 (2005).
[14] M. Secondini, E. Forestieri, C. R. Menyuk, “A Combined Regular-Logarithmic Perturbation Method for Signal-Noise Interaction in Amplified Optical Systems,” J. Lightwave Technol. 27, 3358-3369 (2009).
[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).
[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).
[17] J. N. Damask, “Polarization Optics in Telecommunications,” (Springer 2005).
[18] S. G. Evangelides Jr., L.F. Mollenauer, J. P. Gordon, N. S. Bergano, “Polarization multiplexing with solitons,” J. Lightwave Technol. 10, 28-35 (1992).
[19] 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).
[20] J. Wang and J. M. Kahn, “Impact of chromatic and polarization-mode dispersions on DPSK systems using interferometric demodulation and direct detection,” J. Lightwave Technol. 22, 362-371 (2004).
[21] 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.
[22] 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.

Appendix A 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 $¨ o$ dinger equation (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_{c o r r}$ ) (the length within which birefringence axes are randomly reoriented), beating length $Λ_{b e a t}$ (the length determined by birefringence strength), and the PMD parameter $D_{P M D}$ (ps/ $\sqrt{km}$ ) [1]-[3]. Ignoring the birefringence effect, the two dimensional (2D) CNLSE can be reduced to one dimensional nonlinear Schr $¨ o$ dinger equation (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 $Λ_{b e a t}$ (10 $\sim$ 100 m) and $L_{c o r r}$ (0.3 $\sim$ 300 m) are much smaller than $L_{D}$ and $L_{N}$ (up to hundreds of kilometers) [1, 2, 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_{c o r r}$ and $Λ_{b e a t}$ 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 the approach of Ref. [2] can be increased to a scale significantly larger than $Λ_{b e a t}$ and $L_{c o r r}$ . (Detailed value of its step size depends on the interaction between linear and nonlinear effects.) Requirements for this approach are that the circular component of the local birefringence in the fiber is negligible and that the nonlinear Kerr effect can be approximated as its statistical average over the Poincar $´ e$ sphere (named Manakov-PMD approximation or M-PMD approx in this work). 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]

e^{A} e^{B} = e^{A + B + \frac{1}{2} [A, B] + \frac{1}{12} [A - B, [A, B]] + . . .},

(1)

where $[A, B] \equiv A B - B A$ . Because of this, even in the case of zero birefringence, the step size using approach [2] is restricted to $10 % \sim 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 D.2).

The concept of interaction picture (IP) was originally used in quantum mechanics. Combined with the fourth-order Runge-Kutta (RK4) algorithm, the method of RK4 in IP (RK4IP) was recently proposed to study the supercontinuum generation in optical fibers (for 1D field cases) [15]. Since 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 the Kerr nonlinearity without approximation. With the help of the local error method proposed in [12], the local error of our RK4IP can be improved to $\sim O (h^{6})$ , rather than $\sim O (h^{5})$ 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_{c o r r}$ and $Λ_{b e a t}$ . Without birefringence, the step size of RK4IP can be the same order of magnitude as $L_{D}$ and/or $L_{N}$ , or, $6 \sim 10$ times the step size using the approach of Ref. [2].

In the parameter regime where communication fibers operate, our results are consistent with those using M-PMD approx, which was used in many publications (cf. e.g., [1, 2, 3]). Increased interaction between the linear and nonlinear terms in CNLSE will lead to increased discrepancy between RK4IP and M-PMD approx.

Appendix B From CNLSE to M-PMD approx

b.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 [17]. Namely, given a field with its $x - y$ components being $u_{x}$ and $u_{y}$ , it can be denoted as $| u ⟩ \equiv (u_{x}, u_{y})^{T}$ [or $⟨ u | \equiv (u_{x}^{*}, u_{y}^{*})$ ]. Thus, the CNLSE discussed in [1, 2, 3] can be written in the following form with retarded time $t = t_{l a b} - β_{ω} z$ :

j | u ⟩_{z} + Σ [Δ β | u ⟩ + j Δ β_{ω} | u ⟩_{t}] - \frac{β_{ω ω}}{2} | u ⟩_{t t} + γ [\frac{5}{6} ⟨ u | u ⟩ | u ⟩ + \frac{}{1} 6 σ_{3} | u ⟩ ⟨ u | σ_{3} | u ⟩ + \frac{1}{3} | v ⟩] = 0,

(2)

j | u ⟩_{z} + Σ [Δ β | u ⟩ + j Δ β_{ω} | u ⟩_{t}] - \frac{β_{ω ω}}{2} | u ⟩_{t t} + γ [⟨ u | u ⟩ | u ⟩ - \frac{1}{3} σ_{31} | u ⟩ ⟨ u | σ_{13} | u ⟩] = 0,

(3)

where $| u ⟩_{z} \equiv \partial | u ⟩ / \partial z$ , $| v ⟩ \equiv (u_{x}^{*} u_{y}^{2}, u_{x}^{2} u_{y}^{*})^{T}$ , and $Δ β$ ( $Δ β_{ω} \equiv \partial Δ β / \partial ω$ ) is related to $Λ_{b e a t}$ ( $L_{c o r r}$ ) with $Δ β = π / Λ_{b e a t}$ ( $Δ β_{ω} = D_{P M D} / \sqrt{8 L_{c o r r}}$ ), respectively [2]. Here $D_{P M D}$ is the PMD coefficient (ps/ $\sqrt{km}$ ). The average DGD of a $L$ -km-long fiber can be obtained using DGD $_{a v g} = D_{P M D} \sqrt{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 this work, the three components of the Pauli spin matrices are denoted as

σ_{1} = (\begin{matrix} 0 & 1 1 & 0 \end{matrix}), σ_{2} = (\begin{matrix} 0 & - j j & 0 \end{matrix}), σ_{3} = (\begin{matrix} 1 & 0 0 & - 1 \end{matrix}) .

(4)

Generally, the birefringence-induced matrix $Σ$ in (2) and (3) has the form of [3]

(5)

with $β_{i} \equiv ⟨ β_{0} (z) | σ_{i} | β_{0} (z) ⟩$ ( $i = 1, 2, 3$ ) and $| β_{0} (z) ⟩$ [ ${\to β}_{0} (z)$ ] being, respectively, the unit vector representing the fiber birefringence in 2D Jones (3D Stokes) space [17]. For the fiber with circular birefringence (e.g., spun optical fiber), $β_{2} \neq 0$ in (5).

Parameter $β_{ω ω} = \partial^{2} β / \partial ω^{2}$ corresponds to the CD parameter at wavelength $λ$ by $β_{ω ω}$ =-CD $(λ) λ^{2} / (2 π c)$ (ps/nm $\cdot$ km) with c=3 $\times$ 10 $^{8}$ m/s.

In Eq. (3), $σ_{31} = - σ_{13} = σ_{3} σ_{1} = (\begin{matrix} 0 & 1 - 1 & 0 \end{matrix}) = j σ_{2}$ . They are introduced to show directly that for a unitary transformation, e.g.,

(6)

we have $T^{†} σ_{13} T = σ_{13}$ , provided $T = T^{*}$ . This means the nonlinear term in the CNLSE is covariant for any real rotation. In the following sections, further discussions on CNLSE are based on the form of (3).

When the amplitude of $| u ⟩$ is viewed as the square root of optical pulse power, the nonlinear coefficient $γ$ in Eqs. (2)-(3) relates to Kerr coefficient $n_{2}$ , effective mode area $A_{e f f}$ , and wavenumber $k = 2 π / λ$ by $γ = n_{2} k / A_{e f f}$ . Typical values of $λ = 1550$ nm and $n_{2} = 2.6 \times 10^{- 20}$ m $^{2}$ /W are used throughout of the paper, unless otherwise noted.

b.2 M-PMD approx

Eq. (3) can also be expressed as

	$j \| u ⟩_{z} + Σ [Δ β \| u ⟩ + j Δ β_{ω} \| u ⟩_{t}] - \frac{β_{ω ω}}{2} \| u ⟩_{t t} + γ \frac{8}{9} ⟨ u \| u ⟩ \| u ⟩ = - \frac{γ}{3} [\frac{1}{3} ⟨ u \| u ⟩ \| u ⟩ - σ_{31} \| u ⟩ ⟨ u \| σ_{13} \| u ⟩],$
			(7)

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

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

j | u ⟩_{z} + Σ [Δ β | u ⟩ + j Δ β_{ω} | u ⟩_{t}] - \frac{β_{ω ω}}{2} | u ⟩_{t t} + γ \frac{8}{9} ⟨ u | u ⟩ | u ⟩ = 0.

(8)

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

The form of Eq. (8) is different from, but equivalent to, the M-PMD approx proposed in Refs. [1, 2, 3]. In fact one can obtain the published M-PMD approx, e.g., Eq. (68) of Ref. [3], by substituting the unitary transformations $| u ⟩ = R T | ~ Ψ ⟩$ into (8), with $T$ being given by (6) and $R = (\begin{matrix} cos α & - sin α sin α & cos α \end{matrix})$ , yielding $R_{z} = α_{z} (\begin{matrix} - sin α & - cos α cos α & - sin α \end{matrix})$ , $j R^{†} R_{z} = α_{z} σ_{2}$ , and $[σ_{2} α_{z} + σ_{3} Δ β] T + j T_{z} = 0$ .

Appendix C RK4IP: extension from 1D to 2D

c.1 RK4IP solution for 1D NLSE

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

j \frac{\partial u}{\partial z} - \frac{β_{ω ω}}{2} \frac{\partial^{2} u}{\partial t^{2}} + γ | u |^{2} u = 0,

(9)

which can be formally viewed as

j \frac{\partial u}{\partial z} + (^D +^N) u = 0,^D = - \frac{β_{ω ω}}{2} \frac{\partial^{2}}{\partial t^{2}},^N = γ | u |^{2} .

(10)

Introducing transformation $u = exp [j (z - z_{0})^D] u^{I}$ , NLSE in the IP has the form [15]

j \frac{\partial u^{I}}{\partial z} + {^N}^{I} u^{I} = 0

(11)

\frac{\partial u^{I}}{\partial z} = j {^N}^{I} u^{I} \equiv f (z, u^{I})

(12)

with

{^N}^{I} = exp [- j (z - z_{0})^D]^N exp [j (z - z_{0})^D] \equiv {^I}_{D}^{†}^N {^I}_{D}

(13)

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. (11) can be obtained using RK4, with the local error of $O (h^{5})$ . Choosing $z_{0} = z_{n} + h / 2$ can reduce the required FFTs by $50 %$ (compared to the choice of $z_{0} = z_{n}$ ), yielding [15]

	$u (z_{n} + h, t) = exp [j \frac{h}{2}^D] [u_{n}^{I} + k_{1} h / 6 + k_{2} h / 3 + k_{3} h / 3] + k_{4} h / 6$
	$u_{n}^{I} = exp [j \frac{h}{2}^D] u (z_{n}, t)$
	$k_{1} = f (z_{n}, u_{n}^{I}) = j exp [j \frac{h}{2}^D]^N (u (z_{n}, t)) u (z_{n}, t)$
	$k_{2} = f (z_{n} + \frac{h}{2}, u_{n}^{I} + \frac{h}{2} k_{1}) = j^N (u_{n}^{I} + \frac{h}{2} k_{1}) [u_{n}^{I} + \frac{h}{} 2 k_{1}]$
	$k_{3} = f (z_{n} + \frac{h}{2}, u_{n}^{I} + \frac{h}{2} k_{2}) = j^N (u_{n}^{I} + \frac{h}{2} k_{2}) [u_{n}^{I} + \frac{h}{} 2 k_{2}]$
	$k_{4} = exp [- j \frac{h}{2}^D] f (z_{n} + h, u_{n}^{I} + h k_{3}) = j^N (exp [j \frac{h}{2}^D] [u_{n}^{I} + h k_{3}]) exp [j \frac{h}{2}^D] [u_{n}^{I} + h k_{3}]$		(14)

Note that in (13) and (14), the linear operator $^D$ for 1D case was assumed to be unchanged in $z$ direction. For a fiber with random birefringence, the transformation ${^I}_{D} (z, z_{0}) \equiv exp [j (z - z_{0})^D] = exp [- j (z - z_{0}) β_{ω ω} (\partial^{2} / \partial t^{2}) / 2]$ in (13) needs to be modified correspondly, which is the key point to extend the RK4IP of [15] to 2D case.

c.2 RK4IP solution for 2D CNLSE with random birefringence

Eq. (3) can also be viewed as

	$j \frac{\partial \| u ⟩}{\partial z} + ({^D}_{2} + {^N}_{2}) \| u ⟩ = 0,$
	${^D}_{2} = {\to β}_{0} (z) \cdot \to σ (Δ β + j Δ β_{ω} \frac{\partial}{\partial t}) - \frac{β_{ω ω}}{2} \frac{\partial^{2}}{\partial t^{2}}, {^N}_{2} = γ [⟨ u \| u ⟩ + \frac{1}{3} σ_{13} ⟨ u \| σ_{13} \| u ⟩] .$		(15)

Introducing the transformation $| u_{n} ⟩ =^I (z_{n}, z_{0}) | u_{n}^{I} ⟩$ , with the operator $^I (z, z_{0})$ given by (27), the CNLSE (15) can be expressed in IP as

\frac{\partial | u^{I} (z) ⟩}{\partial z} = j {^N}_{2}^{I} (z) | u^{I} (z) ⟩ \equiv | f (z, u^{I}) ⟩, {^N}_{2}^{I} (z) = {^I}^{†} (z, z_{0}) {^N}_{2} (z)^I (z, z_{0}) .

(16)

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_{0} = z_{n} + h / 2$ to minimize the required number of FFTs, which leads to

	$\| u_{n + 1} ⟩ =^I (z_{n + 1}, z_{0}) \| u_{n + 1}^{I} ⟩ =^d (\frac{h}{2}) {^M}_{2} [\| u_{n}^{I} ⟩ + \| k_{1} ⟩ \frac{h}{6} + \| k_{2} ⟩ \frac{h}{3} + \| k_{3} ⟩ \frac{h}{3}] + \| k_{4} ⟩ \frac{h}{6}$
	$\| u_{n}^{I} ⟩ = {^I}^{†} (z_{n}, z_{0}) \| u_{n} ⟩ = {^d}^{†} (- \frac{h}{2}) {_{L}}^{†} \| u_{n} ⟩ =^d (\frac{h}{2})_{1} \| u_{n} ⟩$
	$\| k_{1} ⟩ = \| f (z_{n}, u_{n}^{I}) ⟩$
	$= j γ^d (\frac{h}{2})_{1} [⟨ u_{n} \| u_{n} ⟩ + \frac{σ_{13}}{3} ⟨ u_{n} \| σ_{13} \| u_{n} ⟩] \| u_{n} ⟩$
	$\| k_{2} ⟩ = \| f (z_{n} + \frac{h}{2}, u_{n}^{I} + \frac{h}{2} k_{1}) ⟩$
	$= j γ [⟨ u_{n} + \frac{h k_{1}}{2} \| u_{n} + \frac{h k_{1}}{2} ⟩ + \frac{σ_{13}}{3} ⟨ u_{n}^{I} + \frac{h k_{1}}{2} \| σ_{13} \| u_{n}^{I} + \frac{h k_{1}}{2} ⟩] \| u_{n}^{I} + \frac{h k_{1}}{2} ⟩$
	$\| k_{3} ⟩ = \| f (z_{n} + \frac{h}{2}, u_{n}^{I} + \frac{h}{2} k_{2}) ⟩$
	$= j γ [⟨ u_{n} + \frac{h k_{2}}{2} \| u_{n} + \frac{h k_{2}}{2} ⟩ + \frac{σ_{13}}{3} ⟨ u_{n}^{I} + \frac{h k_{2}}{2} \| σ_{13} \| u_{n}^{I} + \frac{h k_{2}}{2} ⟩] \| u_{n}^{I} + \frac{h k_{2}}{2} ⟩$

	$\| k_{4} ⟩ =^d (\frac{h}{2})^{†} {^M}_{2}^{†} \| f (z_{n} + h, u_{n}^{I} + h k_{3}) ⟩$
	$= j γ [⟨ v_{n} \| v_{n} ⟩ + \frac{σ_{13}}{3} ⟨ v_{n} \| σ_{13} \| v_{n} ⟩] \| v_{n} ⟩, \| v_{n} ⟩ =^d (\frac{h}{2}) {^M}_{2} \| u_{n}^{I} + h k_{3} ⟩,$		(17)

where $| u_{n}^{I} + c k_{i} ⟩ \equiv | u_{n}^{I} ⟩ + c | k_{i} ⟩$ $(i = 1, 2, 3, c = \frac{h}{2}, h)$ ,

	${^M}_{1} \equiv {^m}_{N_{h / 2}} {^m}_{N_{h / 2} - 1} \dots {^m}_{2} \cdot {^m}_{1}, {^M}_{2} \equiv {^m}_{N_{h}} \cdot {^m}_{N_{h} - 1} \dots {^m}_{N_{h / 2} + 2} \cdot {^m}_{N_{h / 2} + 1},$
	${^M}_{L} \equiv {^m}_{1} (-) \cdot {^m}_{2} (-) \dots {^m}_{N_{h / 2}} (-) = {^m}_{1}^{†} \cdot {^m}_{2}^{†} \dots {^m}_{N_{h / 2}}^{†},$		(18)

and $^d (h / 2)$ can be obtained according to (27). Introduced in (27), $N_{h}$ in (18) is the index of the last plate, whereas the unitary matrix ${^m}_{i}$ in (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 (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 $i$ th plate discussed in the Appendix. Thus we have ${^m}_{i} (-) = {^m}_{i}^{†}$ . Obviously, when $γ = 0$ , (17) yields $| u_{n + 1} ⟩ =^d (\frac{h}{2}) {^M}_{2}^d (\frac{h}{2}) {^M}_{1} | u_{n} ⟩$ , which is consistent with (27) [or (22) and (24)] given in the Appendix.

Appendix D Applications and discussions

In the following numerical calculations, the input optical field $u_{i n} (t) \equiv | {\to u}_{i n} (t) |$ is assumed to be a periodic repetition of $N$ -bit (N=16) de Bruijn sequence, i.e., $u_{i n} (t) = \sum_{n = - \infty}^{\infty} d_{B} (t - n N T_{b})$ , where $d_{B} (t) = \sum_{i = 0}^{N - 1} a_{i} p (t - i T_{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 $i$ th bit. Within the time interval [0, $T_{b}$ ], the elementary forms of RZ and NRZ pulses (or Marks) are assumed to be $p (t) = \sqrt{2 E_{b} / T_{b}} cos [\frac{π}{2} {cos}^{2} (\frac{}{π t} T_{b})]$ and $p (t) = \sqrt{E_{b} / T_{b}}$ , respectively, with $E_{b}$ the optical energy per transmitted bit [19, 20]. 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 all simulations, the fiber loss has been neglected, implying that there are no optical amplifiers in the system and, consequently, no amplifier (ASE) noise.

Only OOK format is considered in this section. Before photodetection, the optical signal is assumed to be filtered by a Fabry-P $´ e$ rot type channel filter with bandwidth $B_{o} = 6.9 / T_{b}$ [20, 21]. The square-law-detected signal is then electrically filtered by a fifth-order Bessel filter with bandwidth $B_{e} = 0.8 / T_{b}$ [3, 21, 22]. In the following figures, the electrically filtered photoelectric pulses are represented in units of W, since detected current corresponds to optical power [2].

In our computation, the approach of adaptive step size is used. Namely, given step size $h$ , one can use RK4IP (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$ ) from $| u (x_{n}) ⟩_{o b}$ (the obtained solution at $x_{n}$ ). 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}$ . According to Ref.[12], the local error of RK4IP can be improved from $O (h^{5})$ [15] to $\sim O (h^{6})$ by using $| u (x_{n} + h) ⟩_{o b} = [16 | u (x_{n} + h) ⟩_{f} - | u (x_{n} + h) ⟩_{c}] / 15$ .

d.1 RK4IP and SSFM comparison: zero birefringence

Fig. 1: Electrically filtered pulse as a function of time ( $T_{b} = 200$ ps), with CD=17 ps/nm $\cdot$ km and fiber birefringence being neglected. The nonlinear coefficient $γ = 1.26$ /(W $\cdot$ km) is obtained using effective mode area $A_{e f f} = 80$ $μ$ m $^{2}$ . 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 (17) and SSFM solution (19).

Without birefringence, 2D CNLSE (3) can be reduced to 1D NLSE (9) with its split-step solution being ( $h = z_{n + 1} - z_{n}$ )

u (z_{n + 1}, t) \approx F^{- 1} [e^{j \frac{β_{ω ω}}{2} ω^{2} h} ~ U], ~ u = u (z_{n}, t) e^{j γ | u (z_{n}) |^{2} h},

(19)

where $~ U$ is the Fourier transformation of $~ u$ , while $F^{- 1}$ denotes inverse Fourier transformation. Note that, in the case of zero birefringence, the SSFM result of Ref. [2] can be reduced to (19).

Obviously, (19) should yield the same result as the RK4IP solution (17) (without birefringence). This is numerically confirmed by considering a NRZ-OOK pulse train launched into the fiber with CD=17 ps/(nm $\cdot$ km) and $A_{e f f}$ =80 $μ$ m $^{2}$ [ $γ$ =1.26 (W $\cdot$ km) $^{- 1}$ ]. The bandwidths of the electrical filter and the input optical filter are $B_{e} = 4$ G and $B = 1.75 B_{e}$ . The input mark power is 20 mW for L=100 km (a) and 2 mW for L=500 km (b). As shown in Fig. 1, the results of RK4IP without birefringence agrees very well with those using SSFM (19).

d.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} \sim 240$ km and the computational accuracy $< 1 %$ of the pulse peak.

	random birefringence	no	birefringence
	RK4IP	RK4IP	SSFM of [2]
$Δ z$ (L $_{N} \sim 40$ km, L=100 km)	7 km	30 km	3 km
$Δ z$ (L $_{N} \sim 400$ km, L=500 km)	45 km	250km	40 km

Table 1 shows the step sizes using RK4IP for the fibers of Fig. 1 (a) and (b) with random birefringence [ $Λ_{b e a t}$ =50m, $L_{c o r r}$ =10m, and $D_{P M D}$ =1.0ps/(km) $^{1 / 2}$ , yielding average DGD $\approx$ 22ps], compared with the step size using RK4IP (without birefringence) and the step size using (19), which is the SSFM of Ref. [2] in the case of zero birefringence. Notice that the fiber dispersion length $L_{D} \sim 1 / (β_{ω ω} Δ f^{2})$ ( $Δ f$ is the signal bandwidth) for the two cases of Fig. 1 is $\sim 240$ km, while the nonlinear length $L_{N} \sim 1 / (γ P)$ is $\sim$ 40 km for the case of Fig. 1 (a) and 400 km for (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 $\sim$ 10 times of the step size of SSFM of Ref. [2], which is around $10 % \sim 15 %$ of L $_{D}$ and/or L $_{N}$ .

Taking into account random birefringence, the RK4IP step size is decreased from 30km to 7 km for the case of L $_{N} \sim 40$ km (L=100km) and from 250 km to 45km for L $_{N} \sim 400$ km (L=500km). In general, the stronger the interaction between the linear and nonlinear parts in CNLSE, the more impact of $L_{c o r r}$ and $Λ_{b e a t}$ 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.)

d.3 RK4IP solutions and M-PMD approximations with random birefringence

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 $_{100}$ =17ps/(nm $\cdot$ km), $A_{e f f}$ =80 $μ$ m $^{2}$ ]. Each span is followed by a 13-km dispersion compensation fiber [CD $_{13}$ =-120ps/(nm $\cdot$ km), $A_{e f f}$ =30 $μ$ m $^{2}$ ]. Before the 1st 100-km fiber, there is a precompensation of -120 $\times$ 6 ps/nm, whereas the 5th 100-km fiber is followed by the compensation of -120 $\times$ 5 ps/nm. The input mark power is 10 mW. Birefringence parameters are $D_{P M D}$ =2.0 ps/(km) $^{1 / 2}$ , $L_{c o r r} = 10$ m, and $Λ_{b e a t} = 50$ m. There is very little difference between the RK4IP solution (17) and the M-PMD approx using (8).

Fig. 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. 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 $Λ_{b e a t} = 50$ m, $L_{c o r r}$ =10 m, and $D_{P M D}$ =2.0 ps/(km) $^{1 / 2}$ , which means average DGD is around 2.0 $\sqrt{5 (100 + 13) - 2} \approx$ 47 ps with the birefringence direction being randomly changed every 10 m. As in D.1, the input optical filter with bandwidth of $B = 1.75 B_{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 solution (17) and the M-PMD approx (8). To confirm that such agreement is not because of the weak enough nonlinearity, one can increase the nonlinear coefficient in (8) from $\frac{8}{9} γ = 1.12$ /(W $\cdot$ km) to $\frac{8}{9} γ = 1.26$ /(W $\cdot$ km) for the transmission fibers and from $\frac{8}{9} γ = 2.99$ /(W $\cdot$ km) to $\frac{8}{9} γ = 3.36$ /(W $\cdot$ km) for the compensation fibers. Our result shows that the discrepancy between the two approximations can be more than 20 $%$ (not plotted in Fig. 2).

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 $_{100}$ =4.5ps/(nm $\cdot$ km), $A_{e f f}$ =72 $μ$ m $^{2}$ ], precompensation of -117 ps/nm before the first 100-km fiber, -429 ps/nm compensation per span, and -234 ps/nm compensation after the 50th 100-km fiber. The input mark power is 2 mW and the PMD coefficient is $D_{P M D}$ =0.40 ps/(km) $^{1 / 2}$ .

Fig. 4: Same as Fig. 3, except that CD $_{100}$ =6.0ps/(nm $\cdot$ km). Because of the change of CD, each span is followed by a -585 ps/nm dispersion compensation. The last 100-km fiber is followed by the compensation of -390 ps/nm.

Fig. 3 shows the RK4IP solution and M-PMD approx for a RZ-OOK system using 50 100-km spans of eLEAF fiber [21], with effective mode area being $A_{e f f}$ =72 $μ$ m $^{2}$ and CD in the range of $3 \sim 6$ ps/(km $\cdot$ nm). To get the received photoelectric pulse shown in Fig. 3, each 100-km fiber with CD $_{100}$ =4.5ps/(km $\cdot$ 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_{c o r r} = 100$ m and $Λ_{b e a t} = 50$ m. As shown in Fig. 3, the agreement between the RK4IP solution and M-PMD approx is not as good as that of Fig. 2. 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 of the CNLSE solution obtained using RK4IP (dashed).

Fig. 4 shows that, when the CD parameter of eLEAF fiber is increased to its worst case [6.0ps/(km $\cdot$ nm)], the discrepancy between RK4IP and M-PMD approx becomes larger than that of Fig. 3, due to the increased interaction between linear and nonlinear terms in the CNLSE.

Appendix E 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. Besides, the Kerr nonlinearity in the CNLSE is treated without approximation. 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 $Λ_{b e a t}$ and $L_{c o r r}$ , 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 \sim 10$ times the step size using the approach of Ref. [2].

In the parameter regime where communication fibers normally operate, our results are consistent well with the results using M-PMD approx (8) [1, 2, 3]. Increased interaction between the linear and nonlinear terms in the CNLSE will lead to increased discrepancy between RK4IP solution and M-PMD approx.

Acknowledgment

The authors acknowledge the financial support from Canadian funding agencies: NSERC.

Appendix: Optical field in the fiber without nonlinearity

For a system with $γ = 0$ , Eq. (3) can be reduced to

j \frac{\partial}{\partial z} | u (z, t) ⟩ + {\to β}_{0} (z) \cdot \to σ (Δ β + j Δ β_{ω} \frac{\partial}{\partial t}) | u (z, t) ⟩ - \frac{β_{ω ω}}{2} \frac{\partial^{2}}{\partial t^{2}} | u (z, t) ⟩ = 0

(20)

In frequency domain, Eq. (20) can be simplified as

j \frac{\partial}{\partial z} | U (z, ω) ⟩ + [{\to β}_{0} (z) \cdot \to σ (Δ β - Δ β_{ω} ω) + \frac{β_{ω ω}}{2} ω^{2}] | U (z, ω) ⟩ = 0,

(21)

where $| U (z, ω) ⟩ = \int | u (z, t) ⟩ e^{- j ω t} d t$ . (Note that in Ref. [2], the Fourier transformation was defined to be $| U (z, ω) ⟩ = \int | u (z, t) ⟩ e^{j ω t} d t$ .)

Assuming

| U (z, ω) ⟩ = M (z, ω) e^{j \frac{β_{ω ω}}{2} ω^{2} (z - z_{0})} | U (z_{0}, ω) ⟩ \equiv I (z, z_{0}) | U (z_{0}, ω) ⟩,

(22)

Eq. (21) is equivalent to the following differential equation for the $2 \times 2$ matrix $M$ [2]:

j \frac{\partial M (z, ω)}{\partial z} + {\to β}_{0} (z) \cdot \to σ (Δ β - Δ β_{ω} ω) M (z, ω) = 0.

(23)

The correct solution of (23) should be in the form of [2]

	$M (z, ω) = m_{N_{h}} \cdot m_{N_{h} - 1} \dots m_{2} \cdot m_{1}$		(24)
	$m_{i} = e^{j δ_{i} η ({\to β}_{0} (z_{i}) \cdot \to σ)} (i = 1, \dots, N_{h}, η \equiv Δ β - Δ β_{ω} ω)$		(25)

where $δ_{i} = z_{i - 1} - z_{i}$ . When $z_{i} > z_{i - 1}$ , $δ_{i}$ is the length of the $i$ th plate with its birefringence direction being represented by ${\to β}_{0} (z_{i})$ . In (24), $N_{h}$ is the number of the plates between $z_{0}$ and $z = z_{0} + h$ and the index $N_{h}$ denotes the $N_{h}$ th plate. The expression of $M (z, ω)$ given by (24) can also be extended to the case of $z - z_{0} < 0$ , provided to notice that 1) for each plate, $δ_{i} < 0$ ; 2) $m_{i}$ ( $i = 1, \dots, N_{h}$ ) are ordered from $z_{0}$ to $z$ , i.e., $m_{1}$ is the first left plate of $z_{0}$ . [ In (25) $η$ is determined by the strengths of the phase birefringence ( $Δ β$ ) and the group birefringence ( $Δ β_{ω}$ ), which is assumed to be independent of $z$ in this work.]

Using a fixed coordinate system in our approach, the result of (24) can be simply obtained by algebra operations such as [17]

	$m_{1} = cos (δ_{1} η) + j sin (δ_{1} η) {\to β}_{0} (z_{1}) \cdot \to σ \equiv μ_{0}^{(1)} + j {\to μ}^{(1)} \cdot σ,$
	$m_{1} m_{2} = [μ_{0}^{(1)} μ_{0}^{(2)} - {\to μ}^{(1)} \cdot {\to μ}^{(2)}] + j \to σ \cdot [μ_{0}^{(1)} {\to μ}^{(2)} + μ_{0}^{(2)} {\to μ}^{(1)} - ({\to μ}^{(1)} \times {\to μ}^{(2)})] \equiv μ_{0}^{(12)} + j \to σ \cdot {\to μ}^{(12)} .$		(26)

Obviously, in time domain, the transform matrices in (22) have the form of ( $i = 1, \dots, N_{h}$ )

	$^I (z, z_{0}) =^M (z, t)^d (z - z_{0})$
	$^d (z - z_{0}) \|_{z = z_{0} + h} = e^{- j h \frac{β_{ω ω}}{2} \frac{\partial^{2}}{\partial t^{2}}}$
	$^M (z, t) = {^m}_{N_{h}} \cdot {^m}_{N_{h} - 1} \dots {^m}_{2} \cdot {^m}_{1}, {^m}_{i} = e^{j δ_{i} (Δ β + j Δ β_{ω} \frac{\partial}{\partial t}) (\to β_{0} (z_{i}) \cdot \to σ)} .$		(27)

Comments 0

Request Comment

You are adding the first comment!

How to quickly get a good reply:

Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.

The feedback must be of minimum 40 characters and the title a minimum of 5 characters

Cancel

Loading ...

	$\| u_{n + 1} ⟩ =^I (z_{n + 1}, z_{0}) \| u_{n + 1}^{I} ⟩ =^d (\frac{h}{2}) {^M}_{2} [\| u_{n}^{I} ⟩ + \| k_{1} ⟩ \frac{h}{6} + \| k_{2} ⟩ \frac{h}{3} + \| k_{3} ⟩ \frac{h}{3}] + \| k_{4} ⟩ \frac{h}{6}$
	$\| u_{n}^{I} ⟩ = {^I}^{†} (z_{n}, z_{0}) \| u_{n} ⟩ = {^d}^{†} (- \frac{h}{2}) {_{L}}^{†} \| u_{n} ⟩ =^d (\frac{h}{2})_{1} \| u_{n} ⟩$
	$\| k_{1} ⟩ = \| f (z_{n}, u_{n}^{I}) ⟩$
	$= j γ^d (\frac{h}{2})_{1} [⟨ u_{n} \| u_{n} ⟩ + \frac{σ_{13}}{3} ⟨ u_{n} \| σ_{13} \| u_{n} ⟩] \| u_{n} ⟩$
	$\| k_{2} ⟩ = \| f (z_{n} + \frac{h}{2}, u_{n}^{I} + \frac{h}{2} k_{1}) ⟩$
	$= j γ [⟨ u_{n} + \frac{h k_{1}}{2} \| u_{n} + \frac{h k_{1}}{2} ⟩ + \frac{σ_{13}}{3} ⟨ u_{n}^{I} + \frac{h k_{1}}{2} \| σ_{13} \| u_{n}^{I} + \frac{h k_{1}}{2} ⟩] \| u_{n}^{I} + \frac{h k_{1}}{2} ⟩$
	$\| k_{3} ⟩ = \| f (z_{n} + \frac{h}{2}, u_{n}^{I} + \frac{h}{2} k_{2}) ⟩$
	$= j γ [⟨ u_{n} + \frac{h k_{2}}{2} \| u_{n} + \frac{h k_{2}}{2} ⟩ + \frac{σ_{13}}{3} ⟨ u_{n}^{I} + \frac{h k_{2}}{2} \| σ_{13} \| u_{n}^{I} + \frac{h k_{2}}{2} ⟩] \| u_{n}^{I} + \frac{h k_{2}}{2} ⟩$

	$\| k_{4} ⟩ =^d (\frac{h}{2})^{†} {^M}_{2}^{†} \| f (z_{n} + h, u_{n}^{I} + h k_{3}) ⟩$
	$= j γ [⟨ v_{n} \| v_{n} ⟩ + \frac{σ_{13}}{3} ⟨ v_{n} \| σ_{13} \| v_{n} ⟩] \| v_{n} ⟩, \| v_{n} ⟩ =^d (\frac{h}{2}) {^M}_{2} \| u_{n}^{I} + h k_{3} ⟩,$		(17)