Analysis of Nonlinear Fiber Kerr Effects for Arbitrary Modulation Formats

Coherent optical transmission systems can be modeled as a four-dimensional (4D) signal space resulting from the two polarization states, each with two quadratures. Recently, nonlinear analytical models have been proposed capable of capturing the impact of Kerr nonlinearity on 4D constellations. None of these addresses the inter-channel nonlinear interference (NLI) imposed by arbitrary modulation formats in multi-channel wavelength division multiplexed (WDM) systems. In this paper, we introduce a general nonlinear model for multi-channel WDM systems that is valid for arbitrary modulation formats, even asymmetric ones. The proposed model converges to the previous models, including the EGN model, in the special case of polarization multiplexed systems. The model focuses on the cross-phase modulation (XPM) nonlinear term that lies at the heart of the NLI in multi-channel WDM systems operating on standard high dispersion single-mode fiber. We show that strategic mappings of the modulation format's coordinates to the polarization states can reduce the NLI undergone by these formats.


I. INTRODUCTION
A NALYTICAL nonlinear channel models in optical fiber communications have been developed that provide a powerful tool to estimate the nonlinear interference (NLI) caused by the Kerr nonlinearity. Although the literature provides many such models, most are restricted to polarization multiplexed (PM) modulation formats. This paper presents an analytical nonlinear model that has the power to predict the NLI in systems using an arbitrary modulation format, including those using asymmetric four dimensional (4D) constellations. The first nonlinear model was introduced in 1993 [1]. Then, an analytical solution to the nonlinear Schrödinger equation employing the Volterra series method was presented both in the time and frequency domains in [2]. Following years of neglect, greater efforts have been made recently to achieve more accurate nonlinear models. The Gaussian noise (GN) model was derived based on the assumption that the transmitted signal in a link follows a Gaussian distribution [3], [4], leading to an overestimate of the NLI. The first 4D GN-like nonlinear model was introduced in [5]. The GN model does not contain any modulation-formatdependent terms. A second-order perturbation technique for the self-phase modulation and cross-phase modulation (XPM) effects was developed in [6]. A modulation-format-dependent time-domain model was proposed for the first time in [7] by resorting to an asymptotic approximation reminiscent of the far-field approximation in paraxial optics. The authors of [8] found that there is a discrepancy between the time domain model in [7] and the GN model [3], and they attributed this deviation to the Gaussianity assumption of the signal in the GN model. To settle this discrepancy, [8] added a modulation-format-dependent correction term to the XPM term. Following the same approach as [8], the authors of [9] added correction terms to the GN model, taking the self-channel interference (SCI), cross-channel interference (XCI), and multi-channel interference (MCI) terms into account, giving rise to an enhanced Gaussian noise (EGN) model.
Other versions of the GN model have emerged to improve its accuracy under different scenarios. Modifications to the GN model to account for the presence of stimulated Raman scattering (SRS) were presented in [10], [11], which are capable of taking into account an arbitrary frequency-dependent signal power profile. These models are valid for Gaussian-modulated signals such as probabilistically-shaped high-order modulation signals. Very recently, [12] proposed an approximate GN model for SCI and XPM in the presence of SRS. The authors of [13], [14] added a modulation format correction term to the XPM, derived in [12,Eq. (8)], while the SCI was computed under a Gaussian assumption. Modulation-format-dependent models in the presence of SRS were proposed in [15], [16], accounting for all the NLI terms, the SCI, XCI, and MCI. The model in [16] introduced a general link function for heterogeneous fiber spans where the span loss and SRS gain/loss are not fully compensated by the amplifier at the end of each span.
Some nonlinear channel models have targeted space-division multiplexing (SDM) systems. An extended version of the GN model for SDM was proposed in [17], irrespective of modulation format dependence and modal dispersion. A modulation-formatdependent nonlinear model was derived in [18] for SDM fibers taking the variance of XPM into account. A comprehensive nonlinear model for SDM as an extension of [18] was introduced This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ in [19], including all the NLI terms such as the SCI, XCI, and MCI terms.
All of the analytical models described above are valid for PM systems but do not apply to 4D modulated systems. In [20], we derived a nonlinear model which quantifies the impact of Kerr nonlinearity on 4D symmetric constellations, accounting for the SCI and XPM nonlinear terms. We extended the model in [20] to a general nonlinear model [21] that accounts for all the NLI terms, including the SCI, XCI, and MCI terms. Although the model given in [21] has the power to take all the NLI terms into account, it still lacks the ability to predict the NLI of constellations that lack symmetry. The contribution of the SCI variance for arbitrary modulation formats was derived in [22], but it does not address the NLI terms that disturb multi-channel wavelength division multiplexed (WDM) systems.
In this paper, we derive a nonlinear model for arbitrary 4D constellations capable of capturing the predominant nonlinear term in multi-channel WDM systems, namely the XPM term. The SCI nonlinear term is caused by the interaction between symbols transmitted within the channel of interest (COI) and is much easier to compensate for [8]. The SCI can, for instance, be reduced by a signal processing technique known as backpropagation [23], [24]. For high-speed optical transmission, where multiple WDM channels occupy the C-band spectrum, the predominant nonlinear term in the high dispersion regime comes from the XPM terms [25]. For these reasons, we concentrate on deriving the XPM in this paper, as in [8]. To derive this model, we remove the restricting assumptions made in [20,Sec. III]. The emphasis in this paper is on the high dispersion regime where, for instance, high symbol rates of around 32 Gbaud are used in a single mode fiber (SMF); the proposed model should not be used for predicting the NLI in the low dispersion regime, where the impact of MCI nonlinear terms is significant. The derivation of a general analytical model that accounts for all NLI terms, such as SCI, XCI, and MCI, for an arbitrary modulation format falls outside the scope of this paper and is left to future work.
The significance of our model is its ability to capture the NLI on any 4D signal space. Not only does this paper give expressions to evaluate the NLI of asymmetric constellations, which has not been done before for multichannel systems, it also tackles the shortcoming in [20], [21] brought about by the assumptions made therein. There exist 4D symmetric constellations whose NLI cannot be predicted by [20], [21] such as BPSK, SP-QAM4_32 [26], SP-QAM4_512 [26], etc. However, we are able to compute the NLI of such constellation through the model proposed in this paper. Furthermore, we provide final expressions that the reader can directly apply to a wide range of purposes, such as geometric and probabilistic shaping. We benchmark the proposed model against the EGN model, and find that the EGN model may inaccurately predict the NLI by about 1.4 dB for a system with 80 WDM channels. Unique to this work, the model can be used to understand how the alignment of the signal constellation with the light polarization affects the severity of the NLI. In particular, we show that a good mapping of the constellation's coordinates to the polarization states can decrease the NLI up to 0.5 dB in certain scenarios.
The structure of the paper is as follows. In Sec. II, we review the first-order solution to the Manakov equation. Sec. III presents the key result of this work, which is an expression for the XCI power. Sec. IV is devoted to simulation results; a wide range of 4D constellations are compared in terms of the experienced XCI in this section. Sec. V provides the conclusion. Finally, the appendix provides a detailed derivation of the NLI model. Notation: We use (·) x and (·) y in this paper to refer to variables related to polarizations x and y, respectively. Two dimensional complex functions are designated by boldface symbols. The conjugate transpose is denoted by (·) † , whereas expectations are indicated by E{·}.

II. PRELIMINARIES
To describe the signal transmission, we start from the Manakov equation 1 in which u(t, z) is associated with the electrical field E(t, z), given in [20, Eq. (1)], by rescaling it to cancel out the attenuation contribution. The function f (z) is responsible for the link's loss/gain profile, which is equal to 1 in the case of perfectly distributed amplification, and equal to exp{−αmod(z, L)} in the case of lumped amplification where α is the loss coefficient, L is the span length, and mod(z, L) indicates the distance between the point z and the nearest preceding amplifier. In (1), β 2 is the group velocity dispersion and γ is the fiber nonlinearity coefficient. Throughout this work, our main focus is on interference due to the XPM effect. This effect involves only two-channel interactions, and as a result, the NLI contributions of multiple WDM channels add up independently. We thus carry out our initial analysis with only two channels, of which one is the COI, whose central frequency is arbitrarily set to zero, and the other is an interfering channel with central frequency Ω. The linear solution to (1) for two channels is then expressed as where a k and b k represent the k-th symbol transmitted in the COI and the interfering channel, respectively. The dispersed pulse waveform at point z along the fiber is is the injected waveform and ∂ 2 t is the time derivative operator. The symbol durations of the COI and interfering channel are denoted by T a and T b , respectively. The COI is matched filtered with a filter whose impulse response is proportional to g * a (L, t). Without loss of generality, we aim to detect the zeroth symbol a 0 .
The extracted symbol at the receiver may be expressed as a 0 + Δa 0 , where Δa 0 accounts for the NLI. By resorting to a perturbation approach, we can write the first-order solution of the Manakov equation as where I is the 2 × 2 identity matrix, and S h,k,l and X h,k,l are expressed as and respectively. The SCI and XPM terms are obtained via the first and second summations on the right-hand side of (3), respectively. We concentrate exclusively on the XPM term because the impact of the SCI was already described in [22]. Using the fact that g(t, z) = dwg(w)exp(−iwt + iw 2 β 2 z/2)/(2π), whereg(w) is the Fourier transform of g(t, 0), (5) can be written in the frequency domain as where d 3 w signifies dw 1 dw 2 dw 3 , and The x-polarized and y-polarized components of the second term of (3) can be assembled into a vector denoted as Δa XPM,0 with components and Δa XPM,0,y (Ω) = i The reader is referred to [20,Eq. 2], [28,Appendix], [8], and [7, Sec. II] to find the origin of (1)-(9).

III. KEY RESULT AT A GLANCE: THE NLI POWER
In this section, we give the final result of the paper, making the resulting expressions easily accessible to the reader. Detailed derivations are relegated to the Appendix.
To obtain the key result, we remove some of the simplifying assumptions made in [20, Sec. III] except the following. We first assume that the data symbols in the x-and y-polarization are correlated with each other. Channels across the spectrum can have different 4D modulation formats. The modulations are assumed to be zero mean, i.e., E{a x } = E{a y } = 0. In our expressions, channels within the spectrum have the same launch power, an assumption which can be easily removed to generalize the results. The key result is obtained for Nyquist rectangular spectral shape channels (sinc pulses) [3], [8], [9], yet we note that our model has the ability to compute the nonlinear disturbance resulting from near rectangular signal spectral shapes, such as a root raised cosine with a small roll-off factor [3]. Note that we no longer assume that the modulation's constellation is symmetric with respect to the two polarizations.
Given (3), the NLI covariance of the zeroth symbol of the COI is given by in which we suppressed the indication of the Ω-dependence in Δa XPM,0 for notational convenience. Because E{Δa XPM,0 } is equal to zero, (10) can be written as Although the non-diagonal terms of the covariance matrix in (11) are non-zero, the power of the NLI on the COI caused by the second term of (3) depends only on the diagonal terms and can be written as where σ 2 XPM,x and σ 2 XPM,y are the XPM variances on polarizations x and y, respectively. The term σ 2 XPM,x , given in (13), results in the final expression where P x is the launch power in polarization x so that the total optical transmit power becomes The terms χ 1 (Ω), Z(Ω), and χ 2 (Ω) in Table I depend on the spectral properties of the signal. The terms Φ 1 , Φ 2 and Φ 3 in Table II, on the other hand, depend on the modulation format.   (14) The term σ 2 XPM,y , given in (13), can be obtained from (14) by swapping x and y in these equations and the terms given in Table II. A detailed derivation of (14) is given in the Appendix. In the special case of independent polarizations where the same format is used in both polarizations, Table II reduces to Φ 1 = 5E{|b x | 4 }/E 2 {|b x | 2 } − 10, Φ 2 = 6, and Φ 3 = 0. These values used in combination with the integral expressions in Table I can be shown to coincide with the EGN model.
As mentioned above, the NLI contributions stemming from multiple channels in a multi-channel WDM system add up independently. The total NLI power on the n-th channel in the spectrum resulting from the XPM contributions of N WDM channels can be expressed as where the function P XPM (·) is provided in (12). In (16), ν i is the central frequency of channel i.

IV. NUMERICAL RESULTS
In this section, we first numerically validate our model using the split-step Fourier method (SSFM), described in detail in [20, Sec. III], for an optical fiber communication system accommodating 80 WDM channels, using the parameters listed in Table III. We then compare a wide range of 4D modulation formats regarding the NLI experienced in a fully-loaded C band transmission system. Lastly, we investigate how different mappings of the constellation's coordinates to the polarization states may affect the NLI experienced.

A. SSFM Validation
In this section, we compare our model with the classical EGN model (the XPM term given in [9, Eqs. (14)-(17)]) as a benchmark in a fully-loaded C-band transmission. Fig. 1 shows the SNR of the COI, channel n = 40, where σ 2 ASE is the amplified spontaneous emission noise generated by the EDFA amplifiers along the link, as a function of launch power for 16-, 256-, and 4096-point constellations. In order to validate (17), SSFM numerical simulations were conducted. Specifically, (17) can be estimated by the simulated SNR of the 40th channel; for a constellation with M symbols, the SNR was estimated through where X and Y are the random variables representing the transmitted and received symbols, respectively, x i is the i-th constellation point, andȳ i = E{Y |X = x i }. A total of 2 15 symbols were simulated per data point, of which the first and last 1500 symbols were removed from the transmitted and received sequences. All channels used the same launch power. In Fig. 1, we mark the SSFM results as filled circles. The modulation formats considered in our simulations are PM-QPSK, subset optimized PM-QPSK (SO-PM-QPSK) [29],  [30], c4_16 [31], [32], PM-16QAM, voronoi4_256 [33], [34], a4_256 [35], w4_256 [36], PM-64QAM, and a4_4096 [32].
Our proposed model, labeled as the 4D model, closely follows the results obtained via SSFM for all types of 4D formats shown in Fig. 1. By contrast, the gap between the conventional EGN model results (marked as dashed curves in Fig. 1) and the SSFM results shows the obvious shortcoming of the EGN model in predicting the NLI of 4D formats. From these results, we expect our model to be precise enough to predict the NLI in systems using any modulation type and operating in the high dispersion regime in which the chief NLI terms are the SCI and XPM. We found in our simulations that the discrepancy between the results of the model proposed and the SSFM results becomes greater as the channel spacing decreases. For systems operating at low symbol rates, this gap also increases. The main reason for this deviation is the impact of MCI nonlinear terms ignored in this paper. Note that for PM modulation formats that do not violate the assumptions made in [20], [21], our model yields the same results as the EGN model.

B. Analysis of the NLI Undergone by Channels Across the Spectrum
In this section we further analyze the experienced NLI of all channels across the spectrum using the parameters listed in Table III. The 4D modulation formats tested, selected from [30], are compared using the NLI noise experienced by channel n normalized by P −3 , where P NLI,n , defined in (16), is the NLI power at the COI. We assume that the whole spectrum can accommodate 80 WDM channels in a system with 50 GHz channel spacing and a symbol rate of 32 Gbaud. We use (13)- (14), and Tables I and II to evaluate the NLI of the 4D constellations. The figures in this section show η n as a function of the spectral location of the COI, identified by the channel number n. In Fig. 2, we compare 16-point, 256-point, and 4096-point constellations in terms of the NLI experienced. We first note that the symmetry of the results in Fig. 2 indicates that the channels located at either edge of the spectrum experience lower NLI than the ones in the middle. The EGN model results are marked as dashed lines. As can be seen in Fig. 2(a), the EGN model results are inaccurate for estimating the NLI for the 4D 16-point constellations considered. The EGN model overestimates the NLI encountered by c4_16 format. In contrast, the EGN model underestimates the NLI of SO-PM-QPSK. The difference between the NLI obtained from the EGN model and our proposed model is more pronounced for the c4_16 format, with a gap of about 1.4 dB as shown in Fig. 2(a). Fig. 2(a) also shows that all three 4D formats experience higher NLI than PM-QPSK, meaning that PM-QPSK outperforms 4D peers. SO-PM-QPSK is the most nonlinearity-prone constellation; the difference between the experienced NLI for SO-PM-QPSK and PM-QPSK is about 1.70 dB. The c4_16 format proposed in [31], although exquisitely evolved to increase the power efficiency, is more vulnerable to the destructive effect  of Kerr nonlinearity than its 2D counterpart (PM-QPSK). Although the mean of the c4_16 constellation is not exactly zero, is it close enough so that the 4D model is able to accurately approximate the NLI experienced. The experienced NLI disparities between different 4D modulation formats is attributed to Φ 1 , Φ 2 , and Φ 3 , defined in Table II and used in (14). Table IV quantifies the influence of Φ 1 , Φ 2 , and Φ 3 , and as a result, χ 1 , χ 2 , and Z, given in (14), on the NLI. Among the terms shown in this table, Φ 3 has the lowest value, implying that χ 2 in (14) has the least impact on the NLI. The terms Φ 1 and Φ 2 , on the other hand, are the dominant factors affecting the NLI. Fig. 2(b) shows the NLI undergone by PM-16QAM, voroni4-256, a4_256, and w4_256 formats. In this case, the EGN model overestimates the NLI of the 4D modulations, voronoi4_256 and a4_256 (rotated w4_256), by around 0.90 dB and 0.80 dB, respectively. Unlike the 16-point constellations, the 256-point 2D modulation, PM-16QAM, is at a disadvantage in comparison with its 4D peers. The difference in NLI between the PM-16QAM and the a4_256/w4_256 formats is about 0.3 dB. This deviation may be rooted in the value of Φ 1 , shown in Table IV, which is smaller for a4_256/w4_256 Recall that a4_256 is equivalent to w4_256 by rotation [30]. As can be seen in Fig. 2(b), there is no change in the NLI between a4_256 and w4_256, meaning that rotations of the constellation do not affect the NLI. Fig. 2(c) finally compares PM-64QAM and a4_4096 [30] regarding the experienced NLI. An inaccuracy of around 0.91 dB is seen in the NLI predicted by the EGN model, leading to an overestimate of the NLI. Looking at this figure, we can see that a4_4096 is more resistant to NLI than PM-64QAM and experiences roughly 0.50 dB lower NLI than PM-64QAM.

C. Influence of a Constellation's Coordinates on the NLI
In this section, we investigate the effects that a constellation's coordinates may have on the NLI. The symbol alphabet, or constellation, of a 4D modulation format with M symbols is given by the set of vectors where c k = (c k,1 , c k,2 , c k,3 , c k,4 ), of which two are mapped onto the x polarization and the other two onto the y. Different mappings have different tolerance to the NLI depending on the modulation format; we noted a strong impact of this mapping particularly on l4_13, c4_13, and l4_18 formats studied. The reader interested in visualizing how 4D symbols map to each polarization is invited to see an exquisite demonstration in [37,Fig. 3] and [38, Fig. 2(a)].
Two different mappings are shown in Fig. 3. The first one consists of mapping c k,1 , c k,2 to the x polarization and the other two to the y. The second one maps c k,1 , c k, 4 to the x polarization and the other two to the y. The term Φ 1 , given in Table II, for the second mapping is higher than for the first one, which explains why this mapping generates higher NLI. As can be seen in Fig. 3(a), the NLI that disturbs l4_13 increases from around 37 dB for the first mapping to 37.55 dB for the second mapping, a difference of more than half a dB, surprisingly. There is also an increase of 0.20 dB in the NLI obtained via the second mapping compared to the first mapping in the case of c4_13, as indicated in Fig. 3(b). This discrepancy for l4_18 is less, only about 0.10 dB, as shown in Fig. 3(c). These results indicate that a good mapping of the constellation's coordinates to the polarization states might curb the impact of the NLI. Interestingly, the experienced NLI of the constellations presented in [37], [38] remains unchanged under different mappings.

V. CONCLUSION
A detailed derivation of a general analytical nonlinear model for 4D formats is given in this paper. The derived model has the ability to quantify the impact of Kerr nonlinearity on a 4D signal space, irrespective of symmetries. The interpolarization dependency had to be taken into account to derive this model. Numerical results show that the EGN model overestimates the NLI by around 1.4 dB in the case of c4_16 for a system with 80 WDM channels. This erroneous prediction is ascribed to polarization dependency, ignored in the EGN model. We also show that l4_16, c4_16 and SO-PM-QPSK modulations have higher NLI than PM-QPSK; the SO-PM-QPSK format experiences the highest NLI amongst 16-point constellations. Because the model can capture the nonlinear disturbance of an arbitrary 4D format, it it can uniquely be used to quantify the influence of a constellation's coordinates on the NLI. The model presented in this paper is valid for high dispersion regimes where the majority of the NLI stems from the XPM terms. Extending this model to transmissions over low-dispersion fiber is the subject of future research.

APPENDIX
This appendix is devoted to evaluating σ 2 XPM,x (Ω) in (14) for a single pair of channels with fixed separation Ω. Thorough this section, the dependence on Ω is left out for notational convenience.
The variance of the perturbative term given in (8) can be written as For the sake of brevity, we only give the procedure to calculate the second term of (21), and the same approach can be followed for the other terms. We focus on calculating this term because it is more general to compute than the first term.
where the second order moment E{a l,x a * l ,y } = E{a x a * y }δ l,l (see [39,Appendix A]). To compute the fourth order moment, the following cases should be considered: Combining (6), (23), and (22) gives h =k e i(w 1 −w 1 )hT −(w 2 −w 2 )kT + 2E{b * Using [39,Eqs. (15) and (29) The term involving δ(w 1 − w 2 )δ(w 1 − w 2 ) is a bias term resulting in a constant phase shift, and should be ignored [ where χ 1 , Z and χ 2 are expressed in Table I. The same approach can be employed for the other terms in (21). Using the fact that E{|a x | 2 } = E{|b x | 2 } = P x , (21) is expressed as (14).