Moment-generating function method used to accurately evaluate the impact of the linearized optical noise amplified by EDFAs

In a nonlinear optical fiber communication (OFC) system with signal power much stronger than noise power, the noise field in the fiber can be described by linearized noise equation (LNE). In this case, the noise impact on the system performance can be evaluated by moment-generating function (MGF) method. Many published MGF calculations were based on the LNE using continuous wave (CW) approximation, where the modulated signal needs to be artificially simplified as an unmodulated signal. Results thus obtained should be treated carefully. More reliable results can be obtained by improving the CW-based LNE with the accurate LNE proposed by Holzlohner et al in Ref. [1]. In this work we show that, for the case of linearized noise amplified by EDFAs, its MGF can be calculated by obtaining the noise propagation information directly from the accurate LNE. Our results agree well with the experimental data of multi-span DPSK systems.


Introduction
The amplified spontaneous emission (ASE) noise from optical amplifiers, e.g., Erbium Doped Fiber Amplifiers (EDFAs), is one of the fundamental reasons for the bit-error-rate (BER) in an optical fiber communication (OFC) system. For an OFC system with non-negligible Kerr nonlinearity, the ASE impact evaluation is complicated due to the nonlinear interaction between signal and ASE noise. In the case of signal power much stronger than the noise power, the noise-noise beating is relatively small so that the noise field in the fiber can be approximately described by linearized noise equation (LNE), which was proposed separately in Ref. [1] and Refs. [2]- [8].
Noise propagator is a matrix used to show noise field propagation in the fiber. In the case of linear perturbation, it is independent of the specific noise realizations, which makes it possible to calculate the moment-generating function (MGF) of the filtered photoelectric current at the receiver [1]- [8].
MGF method is an approach making use of MGF to evaluate the noise impact on the system performance. Now, it is well known that this method is accurate for various linear OFC systems. For nonlinear OFC systems, the MGF method can also be applied, provided their noise fields obey LNE. One can see this from Doob's theorem, which means that, in a linearizable system driven by Gaussian-distributed noise, each of the independent random variable keeps Gaussian (cf. P. 35 of Ref. [9]). Thus the MGF of the received photoelectric current can be calculated as those for linear OFC systems.
The common form of LNE [2]- [8] is based on the continuous wave (CW) assumption, i.e., the noise-free signal in the LNE is artificially simplified as a CW wave. As a result, a semianalytical form of noise propagator and the noise power spectral density (PSD), so called parametric gain (PG), can be obtained. The drawback of this simplification is that the noise-free signal in this LNE neglects chromatic dispersion (CD) effect. As a result, the couplings between noise components (in frequency domain) cannot be taken into account.
A LNE beyond CW, named accurate LNE in this work, was first proposed and discussed in Ref. [1]. Dynamically taking into account the local CD and Kerr nonlinearity along the fiber, this LNE provides accurate noise information, with its computational cost being much higher than the CW approach. For example, given a noise-free signal obtained from nonlinear Schrödinger equation (NLSE), the computation required to update the accurate LNE has cubic complexity in the number of Fourier components [1]. To reduce the computational complexity, covariance matrix method (CMM) was proposed in Ref. [1], where the noise covariance matrix was obtained by processing large noise realizations. In Refs. [10,11], the computational cost of CMM was further reduced by a deterministic approach using perturbation solution. Since the raw covariance matrix obtained from NLSE via Monte Carlo noise realizations [1] or perturbation solution [10,11] may contain nonlinear noise contribution, it needs to be separated from the nonlinearity-induced phase and timing jitter. Thus, the obtained pdfs of the receiver voltage agrees well with Monte Carlo simulation [1,10,12].
With the help of linear perturbation, the noise covariance matrix can also be solved from its ordinary differential equation (ODE) proposed by [1,13]. The covariance matrix obtained by solving such linear ODE does not need jitter separation, although this ODE is more complicated than the accurate LNE [13]. So far there is little comparison between the approaches of Ref. [13] and Refs. [1,10,11].
In this work, we simplify the CMM by showing that the noise propagator matrix can be obtained directly from the accurate LNE. Therefore, there is no nonlinearity-induced jitter. To effectively reduce the computational complexity in updating the LNE, one can decompose the Kerr effect related matrix into a symmetric and an antisymmetric matrices [cf. the discussion after Eq. (29) in Appendix A]. Making use of the fourth-order Runge-Kutta in the interaction picture (RK4IP) method [14,15], the accurate LNE can be solved with large step size, as detailed in Sec. 2. We evaluate the impacts of noise propagator on moment-generating function (MGF) and BER in Sec. 3. The accuracy of this new approach depends on how far the linearized noise deviates from the actual noise. To numerically verify this new approach, we consider the BERs in a 20-span DPSK system with nonlinear phase ofΦ N = 0.2π [5] in Sec. 5. Our BER calculations agree well with the published CMM results. In Sec. 5, we also simulate the experiments of the multi-span DPSK systems discussed in Ref. [7] and show that, to fit the experimental data, one needs to take into account the nonlinearity induced phase difference between noise and noise-free signal, which will affect the signal-noise beating significantly.

Noise propagator obtained from the accurate LNE
In an OFC system amplified by EDFAs, the noise propagator is a fundamental matrix that determines the noise impacts on MGF and BER. In this section, we show that the noise propagator matrix in a fiber of length L can be obtained from the accurate LNE given by Eq. (31). For a multi-span OFC system, one needs to introduce an equivalent noise propagator which can be obtained from PG.

Noise propagator in a fiber of length L
The noise propagator in a fiber of length L can be obtained by extending the RK4IP in Refs. [14,15] to the accurate LNE given by Eq. (31) in Appendix A, where the linear operatorL is associated with CD effect, whereas the nonlinear operatorN is caused by Kerr nonlinearity. By introducingã = eL (z−z 0 )ãI andN I = e −L(z−z 0 )N eL (z−z 0 ) , Eq. (31) in the interaction picture (IP) has the form Taking z 0 = z n + h/2 with step size h = z n+1 − z n and denotingã n =ã(z n ),ã n+1 =ã(z n+1 ), [14,15] to to solve Eq. (1) with which means the noise propagator for the fiber of length h = z n+1 − z n can be calculated as For the fiber of length L, the noise propagator has the form Note that the RK4IP used here is different from the RK4IP in Ref. [15], where what to be solved was the noise-free signal (a 1D matrix), while here what we want is the noise propagator (a 2D matrix). The computational complexity for this 2D matrix is O(N 3 w ), due to that eachk i (i = 2, 3, 4) in Eq. (2) needs one dense matrix multiplication. Here N w is the number of Fourier components used for signal representation, as mentioned between Eqs. (29) and (30).

Equivalent noise propagator of a multi-span system
As discussed in Appendix A, the ASE from an EDFA can be modeled as additive white Gaussian noise (AWGN) with its variance given by Eq. (32). Given the ASE injected at the input of a fiber and the noise propagator obtained from Eqs. (2), (4), and (5), the noise PSD (or PG) at the output of a fiber of length L can be written as [1,5,7,8] where I is a unit matrix and Eq. (32) has been used. In Eq. (6), p T n is the transpose of p n .
For a K-span system consisting of (K + 1) EDFAs (cf. Fig. 4), its PG has the form where σ 2 in = N in /(2T 0 ) and σ 2 is given by Eqs. (21) and (32). In Eq. (7), the real symmetric matrix PG is positive definite. It can be factorized as where the equivalent noise propagator P n,eq can be obtained either by using Cholesky decomposition or symmetric (square root) decomposition [16]. The latter yields P n,eq = P T n,eq .

MGF calculation
With the noise propagator matrix, one can evaluate the BER in the OFC system by calculating the MGF of the electrically filtered current I(t s ) expressed using Karhunen-Loève series expansion (KLSE). Due to the noise in the OFC system, the received (or filtered photoelectric) current fluctuates around its expectation value. The MGF of such current is a useful form of its probability distribution. To get a simple form of MGF, the received current needs to be expressed using KLSE. For a nonlinear OFC system with its noise being linearizable, the KLSE form of the received current can be formulated as Eq. (37) and the related formulas in Appendix B. All the parameters for the nonlinear case are generalized from those for the linear case, which has been well discussed in Refs. [17,18] and other publications. In Eqs. (36) and (37), the Dirac notation |Z is used to represent the normalized noise (in Re-Im form) expressed by Karhunen-Loève bases (in Re-Im form). Averaging |Z i (i = 1, · · · , 4M n + 2) with formula [17,18] the MGF of the filtered current, denoted as Ψ t s (s) here, can be written as where I(t s ) is the filtered photoelectric current at time t s . It consists of signal-signal beating (y ss ), noise-noise beating (y nn ), and signal-noise beating (y ns ), which are detailed in Eqs. (33) and (37)  With the help of Eq. (10) as well as Eqs. (7) and (8) in Ref. [18], the BER can be calculated.

OSNR at the receiver
For an optical system with ASE power being much larger than other noise sources, the OSNR with reference bandwidth B r (0.1nm) can be calculated as In Eq. (11),P s is the time-averaged (noise free) signal power, while P ASE (B r ) is the noise power within B r . To obtainP s and P ASE (B r ), one needs to notice that the measurement bandwidth B m [e.g., the bandwidth of the transfer function of an optical spectrum analyzer (OSA)] may not be the same as B r . Thus theP s in Eq. (11) becomes the power of the signal filtered by B m , while P ASE (B r ) becomes the ASE filtered by B m and weighted by a factor B r /B m [19]. In a linear optical system, the ASE noise along the fiber can be treated as AWGN. Thus, for the system of Fig. 4, its OSNR can be simply calculated as where P ASE (B m ) is the ASE power within B m and N 0 is given by Eq. (21). In Eq. (12), the filter (B m ) effect on the ASE has been neglected. In a nonlinear optical system, the ASE noise "amplified by" PG cannot be treated as a white noise. Similar to the noise-noise beating given in Eq. (37), the measured ASE power will only relate with the self-beating terms of the noise, which yields Here O m is the low-pass transfer function of the bandpass filter (bandwidth B m ). In Eq. (13), σ 2 and PG are given by Eq. (32) and Eq. (7), respectively.
In the case of traditional OSA-based out-of-band OSNR monitoring, the ASE power can be interpolated using where O m (±∆λ ) is the filter function centered at ±∆λ . When ∆λ = 0, Eq. (14) returns to the ASE power in Eq. (13), where P ASE (B m ) ≡ P ASE (B m , 0).

Applications to DPSK systems
To show that the new approach to get the noise propagator is numerically applicable, we will compare our RK4IP results with the CMM results given by Ref. [5] and with the experimental data given by Ref. [7]. Both consider systems with R b =20 Gb/s, using RZ-50% DPSK modulation. In the receiver, the optical filter is Gaussian type, while the electric filter is the fifth-order Bessel type. In the following calculations, we set T 0 = NT b by changing µ in Eq. (35). This means, given the noise propagator matrix, the computational cost for BER is much higher than that in the linear case. For BER calculations, the length of the de Bruijn sequence is N = 2 5 [17]. Based on the relation between the RK4IP step and the fiber dispersion length (L D ) or nonlinear length (L N ) discussed in Ref. [15] as well as the detailed fiber parameters in the following discussion, we let the RK4IP step for the transmission fiber (h tr ) and that for the DCF fiber (h DCF ) be related with h tr : h DCF = (5 ∼ 6) : 1. The value of h tr and the required computational time will be detailed below.

Comparison with CMM results
We consider the 20-span DPSK system discussed in Ref. [5], where BERs using the CMM and the (improved) CW approaches were plotted against the received OSNR in the Fig. 8 of Ref. [5]. In fact this system is basically the same as the one shown in Fig. 4 in Appendix A, provided that one removes the pre-and postcompensating fibers and their amplifiers in the Fig. 2 of Ref. [5] and removes the first EDFA and N in in our Fig. 4. Thus the first term of PG [in Eq. (7)] needs to be ignored. Like Refs. [5,20], where OSNR was calculated in the absence of PG, we obtain OSNR from Eq. (12) with N in = 0 and (K + 1) being replaced by K. According to Ref.
[5], we change OSNR by changing the n sp in Eq. (21). As plotted in Fig. 4, each span contains a transmission fiber followed by a dispersion-compensating fiber (DCF). The transmission fiber is l =100 km long with its CD parameter D tx = 8 ps/nm/km. Each span is fully compensated. The nonlinear phase accumulated in the fiber, defined asΦ NL = z 0 γ(ξ )P in e −α(ξ )ξ dξ with P in being the time averaged signal power (at the input of the fiber), is 0.2π. The bandwidth of the optical (electrical) filter in the receiver is To let our results be reproducible, we provide, as detailed as possible, other related parameters below. The DCF in each span is 8 km long with D DCF = −100 ps/nm/km. Transmission fiber and DCF are assumed to have same fiber loss (α =0.2 dB/km) and same nonlinear coefficient ( γ =2.0 /W/km). The EDFA in each span is used to compensate the total loss in the fiber of L = (100 + 8) km. Therefore the signal power at the input of each span (P in ) keeps constant. Ignoring the nonlinear phase contribution of DCF [21], we set P in =0.7307 mW, obtained from Φ NL = KγP in (1 − e −αl )/α = 0.2π with K = 20 and l = 100 km. The OSNR is obtained using Eq. (12) with B m /B r = 1.35.
As shown in Fig. 1, the curve using the proposed RK4IP approach (dash-dotted) agrees very well with the CMM curve (solid) given by Ref. [5]. In Fig. 1, the curves using improved CW approach [5] with CW power being transmitted peak power (dashed) and average power (dotted) are plotted for comparison.
For each calculated point in Fig. 1, the CPU time for the noise propagator calculation is ∼1.8 hr with RK4IP step for transmission fiber (20×100 km long) being h tr =3.5 km. The CPU time for each BER calculation is ∼ 0.5 hr. In fact, h tr ranged within 0.3 ∼ 5 km yields almost the same curve.   2. BER versus received OSNR for the multi-span RZ-50% DPSK systems withΦ NL = 0.9. Each span consists of a SMF fiber (42 km long) and a DCF file (7 km long). Other fiber parameters were detailed in Table 1 of Ref. [7]. According to Fig. 7 (b) in Ref. [7], where the experimental curves for 5-, 10-, and 25-span systems were almost the same, here we replot these three curves using a thick solid curve.

Comparison with experimental data
The optical system discussed in Ref. [7] can be modelled by Fig. 4, except that each EDFA should be replaced by an EDFA followed by an optical filter (Gaussian) O lk with bandwidth B lk =5 nm. Also, the input noise N in needs to be filtered by an optical filter O in (Gaussian, B in =3 nm). As a result, in Eq. (7), the noise propagator p n (k + 1)L, kL (k = 0, · · · , K − 1) should be replaced with O lk p n (k + 1)L, kL , P n (K, K) = 1 with P n (K, K) = O lk , and Gσ 2 in P n (K, 0)P T n (K, 0) with Gσ 2 in P n (K, 0)O in O T in P T n (K, 0). Since we only consider the curves plotted in Fig. 7 (b) of Ref. [7], each fiber (L km long) in Fig. 4 contains a SMF (42 km long) followed by a DCF (7 km long). In our calculation, all the related fiber parameters are same as those given in Table 1 of Ref. [7]. In the receiver, the bandwidth of the optical filter is 1.87R b , while the bandwidth of the electrical filter is 0.75R b .
We first consider the back-to-back case. Similar to Ref. [7], we modify the 20Gb/s RZ-50% signal at the transmitter by comparing its calculated spectrum [22] and its measured spectrum [7]. Due to that the input noise N in is filtered by O in , the OSNR is calculated using Eq. (13) with B m /B r = 0.95, yielding the back-to-back RK4IP curve shown in Fig. 2.
For the 5-, 10-, 25-span systems, their accumulated nonlinear phases are calculated according to Eq. (48) of Ref. [7]. Because of the spectral modification of the input signal, the optical power at the input of each span P in = P SMF is smaller than E b /T b , where E b is the energy per bit before the spectral modification. For example, to get nonlinear phaseΦ N = 0.9 for the 25-span system, the fiber input power P in = P SMF should be 1.516 mW, which means E b /T b =0.127  Fig.4, where the delay is T b = 1/R b = 50 ps, the delay in the receiver of Ref. [7] was T ′ b =(24.84 GHz) −1 =40.26 ps. Thus, the DPSK phase factors given in Eq. (34) should be modified as where < δ φ >, given by Eq. (40) in Appendix C, is the nonlinear phase difference between noise and noise-free signal. As shown in Fig. 2, all RK4IP curves (φ 0 = 0.31) agree very well with the experiment results. The ASE power is calculated using Eq. (14) with ∆λ = 2B m . In Eq. (16), φ 0 is a calibration constant that basically shifts the RK4IP curves in the OSNR direction, while < δ φ > determines the slope of the RK4IP curves. To show this, we plot in Fig. 3 the RK4IP results for the 25-span system with ∆ = 0.31− < δ φ > and ∆ = 0. Also, we consider the RK4IP curves using Eq. (13) to calculate ASE power. Our results for the 5-span, 10-span, and 25-span systems confirm that there is almost no difference between the curve using Eq. (14) with φ 0 = 0.31 and the curve using Eq. (13) with φ 0 = 0.57. In Eq. (16), the calibration constant φ 0 can be temporally considered as a fitting parameter. As mentioned above, it basically affects the BER vs OSNR curve in the OSNR direction and is related with the detailed OSNR monitoring technique used in Ref. [7]. As there is few information about its OSNR measurement, evaluation of φ 0 is expected to be discussed elsewhere.
In the above calculation, our CPU time for the noise propagator calculation is ∼0.8 hr with RK4IP step for transmission fibers (25×40 km long) being h tr =6.0 km. The CPU time for each BER calculation is ∼ 0.5 hr. The step size h tr within 0.3 ∼ 10 km will result in almost the same curve.

Summary
For linear OFC systems, MGF method is useful for one to evaluate the noise impacts on BERs. This is true not only because it is computationally efficient for the cases with low BER (e.g., < 10 −9 ) but also it can provide reliable information for the cases using coherent detection, which is now widely used in modern OFC systems. It is now well recognized that traditional Gaussian fitting Q-factor approximation is accurate for OOK detection, while MGF method is accurate for various linear OFC systems.
To extend the MGF method to a nonlinear OFC system, one needs to make sure its noise propagator varies within linear regime. This means the noise-noise interaction needs to be neglected and the noise field should follow LNE.
In CMM of Refs. [1,11,12], the noise propagator was obtained from NLSE. It may contain the nonlinearity-induced jitter, which is beyond the linear regime and should be removed. In this work we simplify the CMM by directly solving the accurate LNE [Eq. (1)] with RK4IP. Like the CW approximations (cf. Refs. [2]- [8] and many others) as well as the approach of Ref. [13], where the noise propagation information obtained from the related linearized noise equation is automatically free from nonlinearity-induced jitter, the noise propagator obtained from accurate LNE also varies within linear regime.
To numerically verify this new approach, we consider a 20-span RZ-DPSK system discussed in Ref. [5]. The BERs obtained using this new RK4IP agree well with those using CMM in Ref. [5]. Taking account the phase difference between noise and noise-free signal leads to quantitative matching between numerical evaluation and experimental result [7].