Polynomial Closed Form Model for Ultra-Wideband Transmission Systems

TL;DR

This paper introduces a polynomial closed-form model for ultra-wideband optical transmission systems that offers improved accuracy and generality, validated against numerical models in complex scenarios, and accompanied by reconfigurable software.

Contribution

The paper presents a novel Polynomial Closed-Form Model (PCFM) that analytically computes the GN/EGN model for ultra-wideband systems using polynomial spatial power profiles, enhancing reliability and applicability.

Findings

High accuracy in multiband scenarios including Raman effects

Effective in systems with multiple lumped losses

Broad applicability demonstrated through validation

Abstract

Ultrafast and accurate physical layer models are essential for designing, optimizing and managing ultra-wideband optical transmission systems. We present a closed-form GN/EGN model, named Polynomial Closed-Form Model (PCFM), improving reliability, accuracy, and generality. The key to deriving PCFM is expressing the spatial power profile of each channel along a span as a polynomial. Then, under reasonable approximations, the integral calculation can be carried out analytically, for any chosen degree of the polynomial. We present a full detailed derivation of the model. We then validate it vs. the numerically integrated GN-model in a challenging multiband (C+L+S) scenario, including Raman amplification and inter-channel Raman scattering. We then show that the approach works well also in the special case of the presence of multiple lumped loss along the fiber. Overall, the approach shows very good accuracy and broad applicability. A software implementing the model, fully reconfigurable to any type of system layout, is available for download under the Creative Commons 4.0 License.

Equations

$$G_{\text{NLI}}^{\text{end}}\left(f\right)\approx \sum _{n_s=1}^{N_s}{G}_{\text{NLI}}^{(n_s),\text{end}}\left(f\right)$$

$$G_{\text{NLI}}^{\text{end}}\left(f\right)=\sum _{n_s=1}^{N_s}{G}_{\text{NLI}}^{(n_s)}\left(f\right)\cdot {\left|\text{H}\left(f;n_s+1,N_s\right)\right|}^2$$

$${\left|\text{H}\left(f;n_s,m_s\right)\right|}^2=\prod _{n_l=n_s}^{m_s}{\left|\text{H}\left(f;n_l,n_l\right)\right|}^2$$

$$\text{H}\left(f;n_l,n_l\right)=\sqrt{{\Gamma }^{(n_l)}\left(f\right)}\cdot e^{{\int }_0^{L_s^{\left(n_l\right)}}{\kappa }^{(n_l)}\left(f,z\right)dz}$$

$${\kappa }^{(n_l)}\left(f,z\right)=-j{\beta }^{(n_l)}\left(f\right)-{{\alpha }_T}^{(n_l)}\left(f,z\right)$$

$$\text{H}(f;n_l,n_l)$$

$$\cdot e^{-{\int }_0^{L_s^{(n_l)}}{\alpha }_T^{(n_l)}(f,z)dz}$$

$${\left|\text{H}\left(f;n_s,m_s\right)\right|}^2=\prod _{n_l=n_s}^{m_s}{\Gamma }^{(n_l)}\left(f\right)e^{-2{\int }_0^{L_s^{\left(n_l\right)}}{{\alpha }_T}^{(n_l)}\left(f,z\right)dz}$$

$$G_{\text{NLI}}^{\text{end}}\left(f\right)$$

$$\cdot \prod _{n_l=n_s+1}^{N_s}{\Gamma }^{(n_l)}\left(f\right)e^{-2{\int }_0^{L_s^{\left(n_l\right)}}{{\alpha }_T}^{(n_l)}\left(f,z\right)dz}$$

$$G_{\text{NLI}}^{\text{end}}\left(f\right)=\sum _{n_s=1}^{N_s}{G}_{\text{NLI}}^{(n_s)}\left(f\right)$$

$${\left\{f_{n_{\text{ch}}}^{(n_s)}\right\}}_{n_{\text{ch}}=1}^{N_{\text{ch}}^{(n_s)}},{\left\{B_{n_{\text{ch}}}^{(n_s)}\right\}}_{n_{\text{ch}}=1}^{N_{\text{ch}}^{(n_s)}}$$

$$G_{\text{NLI}}^{(n_s)}\left(f\right)$$

$$\cdot {\left({\gamma }^{(n_s)}\left(f_1,f_2,f\right)\right)}^2\cdot {\left|{\rho }^{(n_s)}\left(f_1,f_2,f\right)\right|}^{2}d{f}_{1}d{f}_2$$

$${\left|{\rho }^{(n_s)}\left(f_1,f_2,f\right)\right|}^2$$

$$G_{\text{WDM}}^{(n_s)}\left(f\right)=\sum _{n_{\text{ch}}=1}^{N_{\text{ch}}^{(n_s)}}{G}_{\text{WDM,}{n}_{\text{ch}}}^{(n_s)}\cdot {\Pi }_{B_{n_{\text{ch}}}^{(n_s)}}\left(f-f_{n_{\text{ch}}}^{(n_s)}\right)$$

$${\Pi }_B(f)=\{\begin{array}{ll}1,& f\in \left[-\frac{B}{2},\frac{B}{2}\right]\\ 0,& \text{otherwise}\end{array}$$

$$G_{\text{NLI}}^{(n_s)}(f)=\frac{16}{27}{\Gamma }^{(n_s)}(f)$$

$$\sum _{m_{\text{ch}}=1}^{N_{\text{ch}}^{(n_s)}}\sum _{k_{\text{ch}}=1}^{N_{\text{ch}}^{(n_s)}}\sum _{n_{\text{ch}}=1}^{N_{\text{ch}}^{(n_s)}}{G}_{\text{WDM},m_{\text{ch}}}^{(n_s)}{G}_{\text{WDM},k_{\text{ch}}}^{(n_s)}{G}_{\text{WDM},n_{\text{ch}}}^{(n_s)}$$

$${\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\Pi }_{B_{m_{\text{ch}}}^{(n_s)}}(f_1-f_{m_{\text{ch}}}^{(n_s)}){\Pi }_{B_{k_{\text{ch}}}^{(n_s)}}(f_2-f_{k_{\text{ch}}}^{(n_s)})$$

$$\cdot {\Pi }_{B_{n_{\text{ch}}}^{(n_s)}}(f_1+f_2-f-f_{n_{\text{ch}}}^{(n_s)}){\left({\gamma }^{(n_s)}(f_1,f_2,f)\right)}^2$$

$$\cdot {\left|{\rho }^{(n_s)}(f_1,f_2,f)\right|}^{2}d{f}_{1}d{f}_2$$

$$f_1$$

$$\text{which is where }{\Pi }_{B_{m_{\text{ch}}}^{(n_s)}}\left(f_1-f_{m_{\text{ch}}}^{(n_s)}\right)\ne 0,$$

$$f_2$$

$$\text{which is where }{\Pi }_{B_{k_{\text{ch}}}^{(n_s)}}\left(f_2-f_{k_{\text{ch}}}^{(n_s)}\right)\ne 0,$$

$$f_1+f_2$$

$$\text{which is where }{\Pi }_{B_{n_{\text{ch}}}^{(n_s)}}\left(f_1+f_2-f_{\text{CUT}}-f_{n_{\text{ch}}}^{(n_s)}\right)\ne 0.$$

$$x$$

$$\mathcal{X}$$

$$1\le m_{\text{ch}}^{(n_s)},k_{\text{ch}}^{(n_s)},n_{\text{ch}}^{(n_s)}\le N_{\text{ch}}^{(n_s)}$$

$$\&f_{\text{CUT}}=f_{m_{\text{ch}}}^{(n_s)}+f_{k_{\text{ch}}}^{(n_s)}-f_{n_{\text{ch}}}^{(n_s)}\}$$

$$G_{\text{NLI}}^{(n_s)}\left(f_{\text{CUT}}\right)=\sum _{\begin{array}{c}x\end{array}\in \mathcal{X}}{G}_{\text{NLI,}x}^{(n_s)}\left(f_{\text{CUT}}\right)$$

$$G_{\text{NLI,}x}^{(n_s)}\left(f_{\text{CUT}}\right)=\frac{16}{27}{\Gamma }^{(n_s)}\left(f_{\text{CUT}}\right)G_{\text{WDM,}{m}_{\text{ch}}}^{(n_s)}{G}_{\text{WDM,}{k}_{\text{ch}}}^{(n_s)}$$

$$\cdot G_{\text{WDM,}{n}_{\text{ch}}}^{(n_s)}{\int }_{f_{k_{\text{ch}}}^{(n_s)}-B_{k_{\text{ch}}}^{(n_s)}/2}^{f_{k_{\text{ch}}}^{(n_s)}+B_{k_{\text{ch}}}^{(n_s)}/2}{\int }_{f_{m_{\text{ch}}}^{(n_s)}-B_{m_{\text{ch}}}^{(n_s)}/2}^{f_{m_{\text{ch}}}^{(n_s)}+B_{m_{\text{ch}}}^{(n_s)}/2}$$

$${\left({\gamma }^{(n_s)}(f_1,f_2,f_{\text{CUT}})\right)}^2{\left|{\rho }^{(n_s)}(f_1,f_2,f_{\text{CUT}})\right|}^{2}d{f}_{1}d{f}_2$$

$${\gamma }_x^{(n_s)}\stackrel{△}{=}{\gamma }^{(n_s)}(f_1,f_2,f_{\text{CUT}})\approx {\gamma }^{(n_s)}(f_{m_{\text{ch}}}^{(n_s)},f_{k_{\text{ch}}}^{(n_s)},f_{\text{CUT}})$$

$$=\frac{2\pi f_{\text{CUT}}}{c}\cdot \frac{n_2}{A_{\text{eff}}(f_{\text{CUT}},f_{m_{\text{ch}}}^{(n_s)},f_{k_{\text{ch}}}^{(n_s)},f_{n_{\text{ch}}}^{(n_s)})}$$

$$A_{\text{eff}}(f_{\text{CUT}},f_{m_{\text{ch}}}^{(n_s)},f_{k_{\text{ch}}}^{(n_s)},f_{n_{\text{ch}}}^{(n_s)})$$

$$\approx \frac{A_{\text{eff}}(f_{\text{CUT}})+A_{\text{eff}}(f_{m_{\text{ch}}}^{(n_s)})+A_{\text{eff}}(f_{k_{\text{ch}}}^{(n_s)})+A_{\text{eff}}(f_{n_{\text{ch}}}^{(n_s)})}{4}$$

Full text

I Introduction

Ultra-wideband (UWB) optical transmission systems (also called ‘multiband’) serve as a promising technology to meet the ever-increasing capacity demands of modern communication networks. Systems using C+L bands are being widely deployed. Research is exploring adding further bands, such as S and E, and others [1]-[5]. However, UWB systems are affected by very significant inter-channel Raman scattering (ISRS) which together with the deterioration of most propagation parameters towards higher optical frequencies, makes the design of UWB systems challenging. Careful optimization of launch power and other system parameters is then required. In addition, UWB systems greatly benefit from backward Raman amplification, but this in turn requires the optimization of the number, frequency and power of the pumps [2], [6].

Optimization is carried out by means of iterative algorithms that require ultra-fast system performance assessment. Non-linear-interference (NLI) estimation is a key component of performance assessment and typically relies on GN/EGN-type models. Since solving the related integrals numerically is not a practicable option due to the large computational effort, these models are implemented as closed-form approximate formulas (CFMs). Over the decade, mainly two groups have pursued the derivation of CFMs, one from University College London (UCL, see [7], [8]), and one from Politecnico di Torino (PoliTo) in collaboration with CISCO (see [9], [10]). Recently, other groups have started looking into similar investigations too [11]-[13].

This paper reports on a new approach at handling the core integrals of the GN model to obtain a CFM. It consists of expressing the spatial power-profile (SPP) of each channel along a span as a polynomial. Once this is done, some of the core GN model integrals admit closed-form near-exact solutions, involving less drastic approximations than previous CFMs needed to take the SPP into account, especially in the presence of Raman amplification and ISRS. We call the newly derived result ‘PCFM’, where ‘P’ stands for ‘polynomial’.

The PCFM can model UWB systems with ISRS, forward and backward Raman amplification, short spans, lumped loss and other possible SPPs perturbations. It enjoys better accuracy and reliability than previous CFMs, because it does not hinge on either neglecting intra-span NLI coherence or resorting to infinite series expansions, or other problematic approximations. The coherent beat of NLI noise along a single span, especially in the presence of Raman amplification, is fully captured, as opposed for instance to our previous most advanced model, CFM6 [10].

In addition, PCFM makes it possible to remove some of the other approximations that CFMs generally make. In particular, the PCFM approach potentially allows to retain the Multi-Channel Interference (MCI) contribution to NLI, making it possible to handle low-dispersion and low-symbol-rate scenarios where the MCI ‘islands’ in the integration domain of the GN and EGN model [15] core integrals need to be considered. In principle, it can also allow to remove other approximations that are typically made by CFMs, such as that channel spectra are rectangular and that the power spectral density of NLI is flat over each channel. However, these extensions mentioned in this paragraph are left for a specifically devoted forthcoming submission and will not be dealt with here.

This paper is a follow-up to the OFC 2025 invited paper [16]. Compared to [16], we provide a full and detailed derivation of the PCFM. Starting with the GNRF (GN-model reference formula), we work out the complete set of integrals for all types of NLI contributions, which can be solved in closed form using the SPP polynomial representation. In particular, the Self-Channel Interference (SCI) and Cross-Channel Interference (XCI) contributions are explicitly expressed in closed form. We discuss in depth the polynomial fitting of the channel SPPs. We expand on the testing of the model, including the case of multiple lumped losses in a 100 km fiber with significant ISRS and backward Raman amplification.

This paper is organized as follows. In Sect. II, we provide the complete theoretical background involved in the derivation of the PCFM. In Sect. III, we present the PCFM derivation in detail, focusing on the use of a polynomial to express the SPP of each channel. In Sect. IV, we show accuracy tests of SPP representation and of NLI estimation in UWB systems. We provide at the end some hints on what directions the development of the approach is taking and the link to a software implementation which is available for download. Comments and conclusions follow.

II Premises

The NLI noise is produced by the Wavelength Division Multiplexing (WDM) signal in each span. We assume the approximation of incoherent NLI accumulation, that is, the NLI produced in each span sums up in power at the end of the link. This approximation can be removed a posteriori, similar to what was done in [17], but it is necessary in this phase. So, the power spectral density (PSD) of NLI at the end of the link, $G_{\text{NLI}}^{\text{end}}\left(f\right)$ is:

[TABLE]

where ${N}_{s}$ is the number of spans in the link and the ‘approximately equal’ symbol is used here to point out that incoherent NLI accumulation is an approximation. However, henceforth we will use an ‘equal-to’ symbol for ease of readability. Note that $G_{\text{NLI}}^{(n_{s}),\text{end}}(f)$ represents the NLI noise produced within the $n_{s}$-th span, then propagated linearly till the end of the link.

For an illustration of the definition of ‘span’, please see Fig. 1. A span is comprised first of a stretch of fiber, which we assume to be of the same type within the span. It may of course differ from span to span. The span fiber is characterized through its dispersion $\beta^{(n_{s})}(f)$, attenuation $\alpha^{(n_{s})}(f)$ and non-linearity coefficient $\gamma^{(n_{s})}({\bar{f}})$. Note that $\gamma^{(n_{s})}$ depends on multiple frequencies, as it is discussed later, hence the vector symbol above $f$. Both forward and backward-pumped optional Raman amplification may be present in the span, as indicated by the pump injection in Fig. 1.

We then make a distinction between fiber attenuation $\alpha^{(n_{s})}(f)$, which is assumed $z$-independent, and total distributed span attenuation/gain, which also includes possible lumped loss due to splicing and connectors along the span length, as well as forward/backward pumped Raman amplification and ISRS, which is $z$-dependent: $\alpha_{T}^{(n_{s})}(f,z)$.

The span then includes a lumped element at the end, shown as a red triangle. Such lumped-element can account for the following: a Doped-Fiber-Amplifier (DFA), a linear filtering element as well as lumped loss before and after amplification. Overall, the lumped-element is modeled using a frequency-dependent lumped power gain/loss $\Gamma^{(n_{s})}(f)$.

At the input of the span, there is the WDM overall signal PSD $G_{\text{WDM}}^{(n_{s})}(f)$. The output of the $n_{s}\text{-th}$ span is the input of the $(n_{s}+1)\text{-th}$ span $G_{\text{WDM}}^{(n_{s}+1)}(f)$.

We call the transfer function accounting for all linear elements and effects from the input of the $n_{s}\text{-th}$ span to the output of the ${{m}_{s}}\text{-th}$ span as ${\left|\text{H}(f;n_{s},m_{s})\right|}^{2}$. We can then re-write Eq. (1) as follows:

[TABLE]

where $G_{\text{NLI}}^{(n_{s})}(f)$ is the NLI noise produced in the $n_{s}\text{-th}$ span, as it shows up at the end of the $n_{s}\text{-th}$ span, i.e., after the lumped element.

Note that the last term in the summation has index ${{n}_{s}}={{N}_{s}}$, which generates a transfer function: ${{\left|\text{H}\left(f;{{N}_{s}}+1,{{N}_{s}}\right)\right|}^{2}}$. The fact that the first index is larger than the second may appear to generate an inconsistency. We could easily remove the problem by simply imposing ${{\left|\text{H}\left(f;{{N}_{s}}+1,{{N}_{s}}\right)\right|}^{2}}=1$. There is however a justification for doing so. The writing ${{\left|\text{H}\left(f;{{N}_{s}}+1,{{N}_{s}}\right)\right|}^{2}}$ means that we are looking at the transfer function between two points that coincide: the input of the $\left({{N}_{s}}+1\right)$-th span and the output of the ${{N}_{s}}\text{-th}$ span are the same location. Since they spatially coincide, it makes sense that ${{\left|\text{H}\left(f;{{N}_{s}}+1,{{N}_{s}}\right)\right|}^{2}}=1$ .

II-A The linear transfer function

We first focus on expressing ${{\left|\text{H}\left(f;n_{s},m_{s}\right)\right|}^{2}}$ in a general closed form. To express it, we first point out that:

[TABLE]

So, we can concentrate on the simpler problem of finding ${{\left|\text{H}\left(f;n_{l},n_{l}\right)\right|}^{2}}$, which is the square modulus of the linear transfer function from input to output of the $n_{l}$-th span. We have:

[TABLE]

where ${{\kappa}^{(n_{l})}}\left(f,z\right)$ is the complex propagation constant defined as:

[TABLE]

where, for the $n_{l}$-th span, ${{\beta}^{(n_{l})}}\left(f\right)$ is the propagation constant at frequency $f$ and ${{\alpha_{T}}^{(n_{l})}}\left(f,z\right)$ is the fiber attenuation, at frequency $f$ and location $z$ in the span. Note that in Eq. (4) it is assumed that the spatial variable $z$ is re-initialized to zero at the start of each span, that is, it is local to the span. It runs from 0 to the $n_{l}$-th span length $L_{s}^{\left(n_{l}\right)}$.

Also, we remind the reader that distributed gain due to ISRS or Raman amplification is incorporated within ${{\alpha_{T}}^{(n_{l})}}\left(f,z\right)$. This means that there may be stretches of fiber whose ${{\alpha_{T}}^{(n_{l})}}\left(f,z\right)<0.$

Then Eq. (4) simplifies to:

[TABLE]

Eq. (3) then reads:

[TABLE]

and note the important aspect that, due to the absolute value squared on the left-hand side, dispersion disappears completely from Eq. (7). As a consequence, it also disappears from Eq. (2) which becomes:

[TABLE]

This equation embodies our fundamental premises, including the incoherent NLI accumulation assumption. If all spans were transparent, i.e., with gain exactly compensating for loss, at all frequencies, the linear propagation factors would all be 1 and NLI at the receiver would simply be:

[TABLE]

Of course, if transparency is not verified, Eq. (8) accounts for the general case.

II-B NLI PSD

According to the assumptions and derivations of the incoherent GN model [20], Eq. (10) provides the NLI PSD produced by the generic $n_{s}\text{-th}$ span, $G_{\text{NLI}}^{({{n}_{s}})}\left(f\right)$.

The integral contains the factors $(\gamma^{(n_{s})})^{2}$ and $\left|\rho^{(n_{s})}\right|^{2}$ which will be discussed later. We first focus instead on the role of the WDM signal $G_{\text{WDM}}^{(n_{s})}(f)$ at the input of the $n_{s}\text{-th}$ span, which appears three times in the integrand function of Eq. (10). An example of a possible $G_{\text{WDM}}^{(n_{s})}(f)$ is shown in Fig. 2. Note that the superscript $n_{s}$ is present because we assume that in each span there can be a different set of WDM channels, that is, a different $G_{\text{WDM}}^{(n_{s})}(f)$. We also assume that the channel under test (CUT) is present in each span, at the same frequency. All other channels can change, in number, frequency and bandwidth. Specifically, the number of WDM channels in the ${{n}_{s}}\text{-th}$ span is $N_{\text{ch}}^{({{n}_{s}})}$. The set of center frequencies and the set of channel bandwidths in the ${{n}_{s}}\text{-th}$ span are:

[TABLE]

We explicitly write the CUT channel index in the ${{n}_{s}}\text{-th}$ span as $n_{\text{CUT}}^{({{n}_{s}})}$. Of course, $n_{\text{CUT}}^{({{n}_{s}})}$ is one of the channel indices of that span, that is $n_{\text{CUT}}^{({{n}_{s}})}\in\left\{1,\ldots,N_{\text{ch}}^{({{n}_{s}})}\right\}$. The frequency and bandwidth of the CUT in the ${{n}_{s}}\text{-th}$ span should then be written $f_{n_{\text{CUT}}^{({{n}_{s}})}}^{({{n}_{s}})}$, $B_{n_{\text{CUT}}^{({{n}_{s}})}}^{({{n}_{s}})}$. However, the frequency and bandwidth of the CUT do not change span-by-span, so for ease of notation we will write them as: ${{f}_{\text{CUT}}}$, ${{B}_{\text{CUT}}}$, dropping all indices.

We then approximate the channel spectra as rectangles, whose bandwidth is $B_{{{n}_{\text{ch}}}}^{({{n}_{s}})}$. As a result, we can write the WDM comb at the input of the ${{n}_{s}}\text{-th}$ span as:

[TABLE]

where ${{\Pi}_{B}}\left(f\right)$ is a Heaviside ‘pi’ function (i.e., a rectangle) of bandwidth $B$, defined as:

[TABLE]

and $G_{\text{WDM},n_{\text{ch}}}^{({{n}_{s}})}$ is the value of the channel PSD, setting the height of the rectangular spectrum.

Regarding the approximation of a channel spectrum with a rectangle in the context of the GN-model, this topic has been discussed in depth in the past, and possibly for the first time in [14]. As time goes by, this approximation seems to become increasingly immaterial, as modern transceivers operate at very small roll-offs, typically on the order of 0.1 or less. Throughout this paper we assume rectangular spectra.

Introducing Eq. (11) into Eq. (10) we get:

[TABLE]

We then concentrate on calculating $G_{\text{NLI}}^{({{n}_{s}})}\left(f\right)$ related to a single channel, specifically the CUT. Note that the CUT can be any of the WDM comb channels, so there is actually no loss of generality in this assumption.

We also limit our interest to finding the PSD of NLI at the center frequency of the CUT, that is $G_{\text{NLI}}^{({{n}_{s}})}\left({{f}_{\text{CUT}}}\right)$. Given the formalism that we are introducing in this paper, it would be possible to calculate in closed-form the PSD of NLI at frequencies other than the center frequency of the CUT. On the other hand, it has been consistently found (for the first time in [14]) that, while not perfectly, $G_{\text{NLI}}^{({{n}_{s}})}\left(f\right)$ is rather flat over the bandwidth occupied by each channel. As a result, the assumption of ‘locally white’ NLI noise over each channel, with a PSD value equal to the value that it assumes at the center frequency of that channel, is an accepted approximation. We do however mean to investigate the possibility of removing such approximation, in a forthcoming study.

Each of the terms of the triple summation then creates a specific integration domain in ${{f}_{1}}$ and ${{f}_{2}}$. Specifically, considering a generic term of the triple summation, the integration domain for ${{f}_{1}}$ and ${{f}_{2}}$ results from the intersection of three conditions:

[TABLE]

The intersection of the above three conditions creates distinct sub-domains, sometimes called ‘islands’. Looking for instance at an example of a 7-channel WDM comb, made up of equally-spaced, identical-bandwidth channels, where the CUT is the center channel, then the integration islands that contribute to the NLI PSD at the center of the CUT are those shown in Fig. 3-top.

Since each island accounts for the NLI which is produced by the ‘beating’ of three specific channels of the WDM comb, We can therefore label each island by means of an identifier $x$ made up of a triple of integers, where each integer is the index of one of the three WDM channels that produce the NLI contribution related to that specific island. Not all triples of channels create an island, though, because the intersection of the conditions Eq. (14) may turn up empty. The set $\mathcal{X}$ of the triples of channel indices $x$ that do create an islands for a given CUT is formally:

[TABLE]

Degeneracy is possible, in the sense that the three indices need not be all different. With reference to the example of Fig. 3, the red island has indeed all three indices identical, corresponding to the channel index of the CUT, $n^{(n_{s})}_{\rm CUT}$. This island accounts for the beating of the CUT with itself, i.e., SCI. The yellow islands have one index that is the CUT and the other two indices identical, spanning all WDM channels other than the CUT. These islands account for the NLI caused by each single channel onto the CUT, called XCI. The light blue islands have all three indices different from that of the CUT and account for MCI.

The NLI PSD on the CUT in Eq. (13) is calculated by summing up the contribution $G_{\text{NLI,}x}^{(n_{s})}\left(f_{\text{CUT}}\right)$ of each individual island $x$, as follows:

[TABLE]

Each contribution therefore requires integrating over a lozenge-shaped island like the ones in Fig. 3. However, this is quite challenging to do analytically. The problem would be much simpler if we could approximate each island with the square or rectangle that inscribes it, such as shown in Fig. 4.

It can be shown that this approximation does not cause excessive error, for various reasons. One reason is the interesting aspect that the lozenge-shaped integration islands tend to become rectangle-shaped when the channel spacing $\Delta f$ becomes close to the channel bandwidth $B_{\rm ch}$. Specifically, this manifests itself with the appearance of smaller triangular sub-islands, when $\Delta f<3/2\cdot B_{\rm ch}$, which then grow and merge with the lozenge-shaped islands to form a rectangle ([14], Fig. 20). Given the trend towards ultra-high Baud-rate channels, with tight spacing, this condition is increasingly met.

Besides this, there are also other reasons that justify this approximation, even when the spacing is not so tight, based on the decay features of the integrand function $\rho^{(n_{s})}\left(f_{1},f_{2},f_{\text{CUT}}\right)$, originally pointed out in [14]. We will come back to this aspect later.

Converting the islands into the inscribing rectangles, we get:

[TABLE]

We will focus on this integral to derive its closed-form solution for any generic island $x$ in the integration domain.

III PCFM derivation

Focusing on the integrand in Eq. (16), we find the fiber non-linearity coefficient $\gamma^{(n_{s})}(f_{1},f_{2},{f_{\text{CUT}})}$. It is indicated as a function of three frequencies but it actually depends on four frequencies. However, in this context, the fourth frequency would be $(f_{1}+f_{2}-f_{\rm CUT})$. Since it is a function of the other three, we omit to indicate it.

As a reasonable approximation, we assume that $\gamma^{(n_{s})}$ is a constant over each island $x$. Then, according to [18], $\gamma^{(n_{s})}_{x}$ can be written as:

[TABLE]

where $c$ denotes the speed of light in vacuum and $n_{2}$ is the nonlinear refractive index of the fiber. The quantity ${A_{\text{eff}}(f_{\text{CUT}},f_{m_{\text{ch}}}^{(n_{s})},f_{k_{\text{ch}}}^{(n_{s})},f_{n_{\text{ch}}}^{(n_{s})})}$ is the cross-effective area of the four frequency components involved in the non-linear Kerr interaction, which can be approximately expressed as [18], [19]:

[TABLE]

The above formulas show that $\gamma^{(n_{s})}$ depends on the type of NLI contribution. For the SCI island, for instance, all four effective areas at the numerator of Eq. (18) are equal to $A_{\text{eff}}(f_{\text{CUT}})$ and the resulting value of $\gamma^{(n_{s})}$ is the one that is generally used in models:

[TABLE]

For the other islands, somewhat different values are found, that relate to XCI and MCI. To compute the $A_{\text{eff}}$ terms in Eq. (18), we used [12] Eq. (4) (derived from [21]). Note that $A_{\text{eff}}$ depends on the type of fiber and hence on the span index $n_{s}$, but we omitted it to avoid clutter.

Since $\gamma^{(n_{s})}$ is now a constant in each island $x$, we can pull it out of the integrals:

[TABLE]

We then focus on the link function $\rho^{(n_{s})}(f_{1},f_{2},f_{\text{CUT}})$:

[TABLE]

In the following we discuss the various quantities appearing in the above formula.

III-A Dispersion

To accurately capture dispersion in UWB systems, we assume that it is expressed through its fourth-order series expansion as:

[TABLE]

where $f_{c}^{({{n}_{s}})}$ is the arbitrary frequency where the expansion is taken in the ${{n}_{s}}\text{-th}$ span. Substituting Eq. (22) into the expression of $\Delta{{\beta}^{({{n}_{s}})}}$ in Eq. (21), we then get the expression:

[TABLE]

For reasons of convenience, we now change the integration variables of Eq. (20) as follows:

[TABLE]

This is equivalent to placing the origin of the new integration variables $f_{1}^{{}^{\prime}}$ at the center frequency of the CUT. As each channel becomes the CUT, such origin is moved onto its frequency. This way, we obtain:

[TABLE]

This important quantity depends on the value of $\beta_{2}^{({{n}_{s}})}$ and on a ‘correction’ that in turn depends on $\beta_{3}^{({{n}_{s}})}$, $\beta_{4}^{({{n}_{s}})}$. Then, as an approximation, only in the inner bracket we replace $f_{1}^{{}^{\prime}}$ and $f_{2}^{{}^{\prime}}$ with their center values over the integration island $x$, that is:

[TABLE]

We get:

[TABLE]

To better understand the meaning of this approximation, we point out that the result would be exact if $\beta_{2}^{({{n}_{s}})}$ was piecewise constant over the bandwidth of each channel. So, the above approximation does account for the change over frequency of $\Delta{{\beta}^{({{n}_{s}})}}$, due to the presence of a non-zero $\beta_{3}^{({{n}_{s}})}$, $\beta_{4}^{({{n}_{s}})}$. We are simply assuming that dispersion is locally constant over each channel, a reasonable assumption given the typical bandwidth of a channel and the small values of $\beta_{3}^{({{n}_{s}})}$, $\beta_{4}^{({{n}_{s}})}$. This approximation however turns out to be crucial in achieving a closed-form solution of the integrals. We can now define:

[TABLE]

which turns out to be a sort of ‘effective’ value of $\beta_{2}^{({{n}_{s}})}$, constant over the island $x$. Eq. (27) can then be drastically simplified as:

[TABLE]

III-B Loss

We can then apply the same strategy to the loss term:

[TABLE]

Once more, the result would be exact if loss was frequency-flat over the bandwidth of each of the channels involved, a very reasonable assumption.

III-C *Spatial power-profile formalism *

Substituting Eqs. (24), (29), (30) into Eq. (20) yields:

[TABLE]

with the integrand function:

[TABLE]

Note that Eq. (32) is indicated as an approximation of $\mid\rho^{(n_{s})}\mid^{2}$ because two approximated expressions, Eq. (29) and Eq. (30), are substituted into the defining equation for $\mid\rho^{(n_{s})}\mid^{2}$ in Eq. (21). The same is true for Eq. (35) below.

We now define the ‘normalized Spatial Power-Profile’, or SPP, of the ${i}_{\text{ch}}\text{-th}$ channel along the ${n}_{s}\text{-th}$ span as:

[TABLE]

Note that $p_{{{i}_{\text{ch}}}}^{({{n}_{s}})}\left(0\right)=1$. We then point out that such SPP can be expressed in terms of the generalized loss/gain coefficient $\alpha_{T}$ as follows:

[TABLE]

which takes into account the accumulated loss/gain along the propagation fiber, experienced up to length $z$ into the span, due to fiber intrinsic loss, Raman amplification, ISRS and possible lumped loss along the fiber.

Using Eq. (34) in Eq. (32), we get:

[TABLE]

This form of the integrand function departs substantially from what has been so far used or proposed for all other CFMs in the literature. Typically, from Eq. (32) onward, a completely different approach at performing calculations was used, focusing on the loss/gain coefficient $\alpha_{T}$. Here instead, we focus on the SPP of each channel.

We then express the combination of the SPPs under square root in Eq. (35) as a single function which depends on the island $x$ being addressed, $p_{x}^{(n_{s})}(z)$, as follows:

[TABLE]

Using Eqs. (35), (36) in Eq. (31) we obtain the expression for the NLI PSD contribution of any generic island $x$ (approximated as a rectangle, see Fig. 4) in any generic span $n_{s}$:

[TABLE]

This important formula shows that the key analytical hurdle to achieving a CFM is solving for the ‘core’ integrals of Eq. (37), which we call $K_{x}^{(n_{s})}(f_{\text{CUT}})$:

[TABLE]

Henceforth we will focus on analytically dealing with these core integrals.

III-D Integrals for SCI and XCI

The SCI and XCI islands are the red and yellow islands, respectively, in Fig. 3. As discussed in Sect. (II-B), we have approximated each lozenge-shaped island with the rectangle or square that inscribes it, as shown in Fig. 4. This approximation is already reflected in the integration domain of the ‘core’ frequency integrals Eqs. (37) and (38).

Note that, if channels had non-uniform spacing and non-identical bandwidths, the corresponding island layout would be much less regular than shown in Fig. 3. See [22] for a few examples of islands in these cases. However, it is always possible to inscribe an island into a square and therefore the following results are completely general. They can deal with non-uniform combs too, with no change.

III-D1 XCI contribution

We first focus on the yellow XCI islands straddling the horizontal $f^{\prime}_{1}$ axis, i.e., the axis with ${{f}^{\prime}_{2}}=0$. They are found by imposing $k_{\text{ch}}^{(n_{s})}=n_{\text{CUT}}^{(n_{s})}$ and $m_{\text{ch}}^{(n_{s})}=n_{\text{ch}}^{(n_{s})}$, resulting in island identifiers of the form:

[TABLE]

with ${n}_{\text{ch}}=1\ldots N_{\rm ch}$, and ${n}_{\text{ch}}\neq{n}_{\text{CUT}}$ because ${n}_{\text{ch}}={n}_{\text{CUT}}$ would be SCI. We call $\mathcal{X}_{\rm{\rm XCI}_{h}}$ the set of these ‘horizontal axis’ XCI islands.

Focusing on these islands in Eqs. (37) and (38), we obtain:

[TABLE]

[TABLE]

where the nonlinear coefficient, specializing the general expression Eq. (17) to the XCI case, can be written:

[TABLE]

The quantity $A_{\rm eff}(f)$ is the effective area of the fiber at the frequency $f$, for which we used [12] Eq. (4) (derived from [21]). Importantly, the power-profile factor $p_{x}^{(n_{s})}(z)$ becomes simply the SPP of the channel with index $n_{\rm ch}$, once the islands indices $x\in\mathcal{X}_{\rm{\rm XCI}_{h}}$ are plugged into Eq. (36) :

[TABLE]

This way, the core integrals simply become:

[TABLE]

We then focus on the yellow XCI islands straddling the vertical $f^{\prime}_{2}$ axis, i.e., the axis with ${{f}^{\prime}_{1}}=0$. They are found by imposing $m_{\text{ch}}^{(n_{s})}=n_{\text{CUT}}^{(n_{s})}$ and $k_{\text{ch}}^{(n_{s})}=n_{\text{ch}}^{(n_{s})}$, resulting in island identifiers of the form:

[TABLE]

again with ${n}_{\text{ch}}=1\ldots N_{\rm ch}$ and ${n}_{\text{ch}}\neq{n}_{\text{CUT}}$. We call $\mathcal{X}_{\rm{\rm XCI}_{v}}$ the set of these ‘vertical axis’ XCI islands.

Going through the same procedure used in the case of the horizontal axis XCI islands, we get:

[TABLE]

The only difference between Eq. (46) and Eq. (44) is the inverted order of the differentials $df^{\prime}_{1}$ and $df^{\prime}_{2}$. However, a simple change of variables, i.e., swapping $f^{\prime}_{1}$ with $f^{\prime}_{2}$, shows that, given $n_{\text{CUT}}^{(n_{s})}$ and $n_{\text{ch}}^{(n_{s})}$, the result of Eq. (46) and Eq. (44) are identical. This implies that, to account for all of XCI, we only need to carry out the integrals on, for instance, the horizontal axis XCI islands and then multiply the result by two, as depicted in Fig. 5 with reference to the example of Fig. 4.

This also means that we can merge together Eqs. (44) and (46) and simplify the notation for the XCI equations as follows:

[TABLE]

[TABLE]

where we recognize that the only index needed to identify the XCI islands is $n_{\text{ch}}^{(n_{s})}$. Each contribution physically represents the XCI produced by the WDM channel $n_{\text{ch}}^{(n_{s})}$ onto the CUT. Also, $\gamma_{\text{XCI,}n_{\text{ch}}}^{(n_{s})}$ coincides with Eq. (42).

III-E SPP polynomial representation

The next step towards reaching a CFM is to express the SPP $p_{n_{\text{ch}}}^{(n_{s})}(z)$ of the $n_{\text{ch}}^{(n_{s})}$-th channel as a polynomial in the distance variable $z$:

[TABLE]

This step, to the best of our knowledge, has not been previously tried in the context of CFMs. The key aspect of this representation is that it allows fully closed-form solutions for the core NLI integrals to be found, for all types of islands.

These closed-form solutions may however become rather complex, specifically as the polynomial degree $N_{p}$ is increased. In the next section, we present instead a remarkably simple general closed-form solution for XCI, of limited complexity for any polynomial degree $N_{p}$, which, however, requires one further approximation. Such approximation causes negligible error in systems with highly dispersive fibers and high symbol rates (see below) and therefore it covers the majority of modern systems. We will also discuss in detail accuracy vs. $N_{p}$ in Sects. IV-B and IV-C.

III-E1 The XCI contribution

To achieve a simple general closed-form solution for XCI, a further approximation is needed. The approximation consists of stretching the XCI islands to $\left[-\infty,\infty\right]$ in the dimension where their size would be $\left[-{{{B}_{\text{CUT}}}}/{2}\;,{{{B}_{\text{CUT}}}}/{2}\;\right]$. An example is shown in Fig. 6. This may seem a quite substantial approximation but it is very common in CFMs since it typically leads to negligible error in modern high-capacity, long-haul SMF systems.

In Fig. 6-bottom, the contour lines of the integrand (resulting from the inner integral in $z$ in Eq. (44)) are shown, in dB, normalized to 0 dB at their maximum, which occurs all along the axes, assuming 100 km of SMF. Note that it is essentially dispersion, and specifically $\beta_{2,\text{eff},{n_{\text{ch}}}}^{(n_{s})}$, that sets the decay of the integrand moving away from the axes. The more dispersion, the faster the decay. In the figure, the channels are relatively low-rate: 28 GBaud, with spacing 50 GHz. The plot clearly shows that at the edge of the stretched island, a decay of $-30$ dB has already occurred. If we consider the next island closer to the origin, even though less decay has occurred at its edges, it is still a decay greater than $-20$ dB. So, stretching upward that island does not cause major error either. The argument gets stronger as one moves away from the origin.

Note that this general feature of the integrand was originally pointed out already in [14], and has since been exploited in many papers to achieve CFMs of various types. Whether this approximation is acceptable or not, it depends on the system parameters. With low dispersion fibers, the decay of the integrand is slower. But also the vertical size of the island before stretching is important, which is set by the CUT symbol rate. As a reasonable guideline, it must be that the product $\beta_{2,\text{eff},{n_{\text{ch}}}}^{(n_{s})}\cdot R_{\rm CUT}^{2}>0.01$ 1/km. This is achieved for instance for $\beta_{2,\text{eff},{n_{\text{ch}}}}^{(n_{s})}=20$ ps2/km (typical for G.652 SMF) and $R_{\rm CUT}=0.028$ TBaud (or 28 GBaud). With modern systems trending towards 100+ GBaud, even fibers with dispersion as low as 2 ps2/km could be dealt with. However, if dispersion is even lower, or otherwise symbol rates are very low such as in multi-subcarrier systems, then these scenarios fall outside of the scope of this approximation. In that case, the island-stretching approximation must be avoided and PCFM allows that, resulting in more complex but still fully closed-form formulas. A presentation of these results is however left for a separate submission.

Quite remarkably, once the stretching of the islands has been performed, the XCI kernel integrals admit a compact fully closed-form exact solution, with no further approximation:

[TABLE]

Notice that the summations are finite and run over the degree of the polynomial that is needed to describe the SPP. We investigated what degrees are typically needed. The results are shown in Sects. IV-B and IV-C.

III-E2 The SCI contribution

The SCI integration takes place over the red islands represented in the examples of Fig. 3. Formally, it corresponds to the triple:

[TABLE]

As in the case of XCI, we first approximate the SCI lozenge with the square that inscribes it, as shown in Fig. 4. Unlike the XCI case, here it is not possible to stretch the integration domain to infinity, along either one of the axes. The reason is that integration to infinity straddling either one of the axes inevitably results in divergence. The SCI core integrals therefore need to be integrated without any further approximations.

In the SCI island, we have:

[TABLE]

Therefore, taking the above into account in Eqs. (37) and (38), we get for SCI:

[TABLE]

[TABLE]

Similarly to the case of XCI, the SPP-related factor $p_{x}^{n_{s}}(z)$ of Eq. (36) gets drastically simplified to just the SPP of the CUT, which we then express here too as a polynomial:

[TABLE]

Quite remarkably, the core integral $K_{\text{SCI}}^{(n_{s})}$ can be fully integrated analytically and exactly. However, without the possibility to resort to the stretching of the island towards infinity, the expressions for SCI turn out to be rather complex. So far we have been able to derive the full expressions for ${{N}_{p}}$ ranging from 0 to 9. A general parametric form vs. $N_{p}$ is not yet available. You think it is likely possible and we are actively researching it.

In Eq. (55) we show as examples the individual formulas for ${{N}_{p}}=0,1,2,3$. Note that these formulas could be optimized for computational efficiency, such as by appropriate factorization. Different strategies are available to do so, but we will not investigate them here. Also notice the presence of the hypergeometric function ${}_{2}F_{3}$ in the formulas. Computationally, it is not a problem since a single evaluation of it is needed for each island, independently of the value of $N_{p}$. Also, it is implemented in Matlab, which typically takes less than 1 ms to compute it on a conventional PC.

Higher-${{N}_{p}}$ formulas, such as for ${{N}_{p}}=5,9$, which would be too long to reproduce here, are reported in [22]. For a detailed discussion of the appropriate value of ${{N}_{p}}$ please see Sects. IV-B and IV-C.

IV PCFM applications and tests

The polynomial representation of the SPPs in PCFM opens up several avenues of application that have been difficult to deal with otherwise. The problem of accounting for backward-pumped Raman amplification, which required cumbersome and unreliable approximations, is overcome. Low loss fibers, ISRS, forward Raman, but also lumped loss and other critical conditions can be dealt with using these formulas. In the following, we focus on testing PCFM in single-span systems that exhibit the above SPPs behaviors, including lumped loss.

The general methodology of using polynomials has great potential even in situations such as low dispersion or small symbol rates, but this requires using XCI formulas obtained without resorting to stretching the islands to infinity, and the inclusion of MCI islands too. Thanks to the SPP polynomial approach, this is possible in closed-form too, but we leave these special applications for a forthcoming submission.

To carry out the testing of the PCFM presented here, we first define a test system which will be our reference throughout.

IV-A Test system description

The system schematic is shown in Fig. 7. It consists of one span of 100 km, which was the first span in the 10-span system in [29]. It was characterized from the experimental set-up used for CFM5 validation in [28]. Loss and dispersion were measured in the C and L bands and then extrapolated to the S band (Fig. 8(a)) using well-known formulas [23]. The Raman gain spectrum $C_{R}(f,f_{p})$ was experimentally characterized using a pump at $f_{p}$ = 206.5 THz (the green curve in Fig. 8(b)). It was shifted and scaled as a function of $f$ and $f_{p}$, according to [19]. The black line in Fig. 8(b) shows the SCI Kerr non-linearity coefficient $\gamma_{\rm SCI}$ of Eq.(19) vs. frequency. Its values and those of $\gamma_{\rm XCI}$ were obtained as indicated next to Eq. (42).

The WDM band boundaries set here were: L-band 184.50 to 190.35; C-band 190.75 to 196.60; S-band 197.00 to 202.85. Each band contained 50 equally-spaced channels, with a symbol rate of 100 GBaud, roll-off 0.1, and spacing 118.75 GHz. The modulation was taken to be Gaussian-shaped.

Backward Raman pumps were included, consisting of three pumps at 205.1 THz, 21.5 dBm; 211.5 THz, 27.7 dBm; 214.0 THz, 26.6 dBm. Their frequency and power were jointly optimized with the launch power of each channel, to maximize the overall system throughput, please see Fig.4 (a) of [29]. The resulting optimum signal launch power is shown in Fig. 9. The figure clearly shows significant Raman amplification, especially in the C and S bands, as well as obvious signs of ISRS, especially in the L-band, where the low-frequency channels initially propagate at zero loss for a few km.

IV-B Accuracy of the SPP polynomial representation

In multiband systems, the computation of the SPPs is typically based on numerically solving the coupled differential Raman equations, which account for both ISRS and Raman amplification [24]. Conventional solvers, such as the bvp4c algorithm in MATLAB, provide accurate results but they are not optimized for this specific use. We developed a faster SPP estimation algorithm, which we recently reported in [29]. The SPPs in Fig. 9 are computed using such algorithm.

Starting from the SPPs, a 9th-degree polynomial is then used to bestfit the SPP of each channel. Fig. 10 displays the polynomial fit for three channels in each band.

Note that the L-band channels experience very low loss, or no loss, at the start, while the S-band channels experience high loss at the start. This is only partially due to the fiber having lower intrinsic loss in the L-band than in the S-band. Much of the behavior is due to the S-band channels losing power to the benefit of the L-band channels, because of ISRS. C-band and S-band channels also experience substantial gain at the end of the fiber due to backward Raman amplification.

The $N_{p}$ = 9 polynomial fit of Fig. 10 appears excellent across all the different SPPs. To provide an idea of what type of error can otherwise be incurred in SPP fitting when using lower degree polynomials, we focus on one of the most challenging channels of the system. It is channel 100 of Fig. 10(f), which exhibits a very steep power increase towards the end of the span. Its polynomial fitting with $N_{p}$=1, 3, 5 and 7 is displayed in Fig. 11. Expectedly, degrees 1, 3 do not look good. As for 5, it still does not appear fully satisfactory.

However, when the actual error on NLI estimation is computed, it turns out that rather accurate results can be obtained even for relatively low values of $N_{p}$. This seem to indicate that some deviation of the polynomial fitting from the SPP, such as visible for $N_{p}$=5, may be ‘averaged out’ during integration and eventually have limited impact on NLI assessment. We will come back to this in detail in the next section.

IV-C Accuracy of NLI estimation

After SPP polynomial fitting, we use the PCFM to compute a key NLI indicator for all 150 channels in the system. We focus on $\text{GSNR}_{\text{NLI}}$, defined as:

[TABLE]

We then compare the result of $\text{GSNR}_{\text{NLI}}$ obtained through the PCFM with that of the numerically-integrated full EGN model [15], which is our benchmark. We call $\Delta\text{GSNR}_{\text{NLI}}$ the error, in dB scale:

[TABLE]

Before looking at the result, we would like to remind the Readers that an error of $\delta$ dB in $\rm GSNR_{\rm NLI}$ estimation causes an error in total $\rm GSNR$ of only $\delta/3$ dB, if operating at the optimum launch power. This is a well-known result from the early days of NLI modeling [25], [26], which has helped make the use of CFMs more viable over the years, by substantially quelling the error coming from the approximations needed to achieve closed-forms. The factor 1/3 should be kept in mind while looking at any NLI assessment accuracy, and we are going to call it the ‘1/3 rule’ henceforth.

Fig. 12(a) shows the violin plot of $\Delta\text{GSNR}_{\text{NLI}}$ across all channels, where different polynomial degrees are used to fit the SPPs. As $N_{p}$ increases, the error distribution stabilizes already for $N_{p}$= 4. This is a somewhat unexpected result given that the plots of Fig. 11 show significant SPP fitting error still for $N_{p}$= 5. We conjecture that the repeated integration in the core integrals $K_{x}^{(n_{s})}(f_{\text{CUT}})$ averages out the fitting errors in the SPPs. Notwithstanding, the violin plots of Fig. 12(a) suggest that, even in complex systems like the one being probed, relatively low values of $N_{p}$ are potentially usable to get fast preliminary NLI estimation, if speed is critical. The final results could then be validated using higher degrees, such as $N_{p}$= 9.

The violin plots for $N_{p}\geq 4$ have a very small standard deviation of 0.1 dB, an excellent result. However, the picture shows a systematic bias of approximately $-0.5$ dB, meaning that the PCFM overestimates NLI by 0.5 dB, leading to a lower $\rm GSNR_{\rm NLI}$ than the EGN benchmark. This systematic overestimation arises from the two main approximations described in detail in the previous sections, namely: assuming that the NLI spectrum over the CUT is flat at the value that it assumes at the center frequency of the CUT (whereas it tapers off somewhat); approximating the lozenge-shaped integration islands with bigger shapes. On the other hand, the 1/3 rule means that the bias on the total GSNR at optimum launch power would be only about $-0.2$ dB. Note also that, since negative, it would be conservative.

Nonetheless, it would be nice to remove such bias. In previous CFM versions we confronted a similar problem and resorted to an analytical machine-trained correction, which we called the ‘machine-learning coefficient’ $\rho$, introduced in Eq. (14) of [17]. We decided to apply such coefficient to the PCFM too, with no change vs. [17]. As the the violin diagrams of Fig. 12(b) show, the bias disappears almost completely. We are currently re-training $\rho$ specifically for the PCFM so that it is possible that even better results, for both residual bias and standard deviation, could be achieved.

The plots of Fig. 12 also show the violin diagram for CFM6. The result for most channels is comparable to those of the PCFM but there are several outliers whose error is as large as -2.5 dB. This is due to a problem intrinsic to the analytical derivation of CFM6 where a series expansion followed by integration produces a denominator that can go to zero. When calculating NLI values near these points, we conjecture that loss of machine precision occurs, leading to errors. While such singularities are removable by means of ad-hoc workarounds, and the outliers can be brought within the violin main body, there remains the need to check for their occurrence during the calculations and then apply the workarounds, a procedure which is cumbersome. This is one of the reasons that compelled us to pursue a new approach, resulting in the PCFM, which is free from any singularity.

There were other reasons too, as mentioned before, such as the difficulty or impossibility of upgrading CFM6 to situations such as low loss, lumped loss, short spans, multi-subcarrier multiplexing and ultra-low dispersion, where PCFM provides a solution or promises to make one possible. One more important reason had to do with modeling intra-span NLI coherence. This aspect is clarified in the the following.

We explore a new test scenario where we reduce the fiber length to 60 km. All SPPs are shown in Fig. 13 (a). Some are shown again in detail, together with their 9-th degree polynomial fit, in Fig. 13 (b). From the plots, it can be seen that with such shorter fiber length, in the presence of ISRS and Raman amplification, some SPPs experience almost ideal distributed amplification. For instance, channel 25 in the L-band sees a maximum loss of less than 3 dB, and its overall profile along the fiber is remarkably flat.

The violin plot of $\Delta\text{GSNR}_{\text{NLI}}$ is depicted in Fig. 13 (c). The distribution gets stable at $N_{p}$ = 5, again with a quite contained standard deviation of about 0.1 dB. As in the 100 km example, a $-0.5$ dB bias is present. Once again, applying the machine-learning coefficient $\rho$ from [17] the bias is eliminated (Fig. 13 (d)).

What is different between the violin plots of Fig. 13 and Fig. 12 is that CFM6 has a much wider distribution than PCFM, which remains such also after applying the machine-learning coefficient. Note that this greater error is not due to the removable singularity problem of CFM6 previously mentioned. Here it is a fundamental shortcoming, which we explain in the following.

CFM6 resorts to breaking up the span into two sub-spans, approximately in the middle of the span. It then separately calculates NLI for both sub-spans. The two contributions are incoherently summed at the end of the span. This is done to allow CFM6 to model the NLI produced by the power growth at the end of the fiber, because of backward-pumped Raman amplification. As Fig. 12 shows, the approach works well in the 100 km case, where mid-span loss is substantial (8-10 dB, see Fig. 10). When this is not the case, such as in the 60 km system, where mid-span loss is as low as 3 dB, a quite substantial error is incurred by CFM6, as visible in Fig. 13 (c) and (d), because of the neglected coherence between the NLI contributions of the two sub-spans. PCFM has no such problem, since the span is dealt with in its entirety, faithfully accounting for intra-span NLI coherence.

IV-D SPP polynomial representation and lumped loss

Polynomials are clearly effective in accurately describing SPPs subject to ‘smooth’ influences, such as ISRS or Raman gain. It is not immediately evident that they could be equally effective when abrupt or ‘lumped’ perturbations are experienced, such as the case of lumped loss.

To assess this capability, we revisit the 100 km configuration and introduce lumped losses of varying magnitudes at different locations along the fiber. Fig. 14 (a) displays the SPPs of the center channels of each of the three bands, with a 1 dB lumped loss introduced at 10 km, such as due to a bad splice or a lossy connector.

Even though an intrinsically continuous curve cannot perfectly fit a discontinuous profile, for $N_{p}=9$ it appears that the polynomial fit does a reasonable job of closely interpolating the SPP. The polynomial curve then goes through repeated integration and it is conceivable that SPP fitting inaccuracy can be averaged out. The resulting GSNR prediction by the PCFM indeed aligns closely with the GN-model markers, as shown by the blue curve in Fig. 14 (c), with no apparent loss of accuracy with respect to the case of no lumped loss.

A more extreme test includes two lumped losses along the span: 2 dB at 5 km and and 0.5 dB at 97 km. Both losses have a significant impact on the SPPs, the first directly and the second because it affects both the signal and the Raman pumps power. Here too the SPPs are not perfectly fitted by the polynomials. However, the interpolation they provide results again in a final GSNR prediction by PCFM (green curve in Fig. 14 (c)) which appears as accurate as the no-loss case.

Overall, considering the complexity of the SPP shapes of Fig. 14 (a) and (b), characterized by substantial ISRS, strong backward-pumped Raman amplification and multiple lumped loss, resulting in profiles that are very far from the simple decreasing exponential of a conventional fiber, we consider the accuracy displayed by the NLI prediction of Fig. 14 (c) as remarkable. Note that the plot displays $\rm GSNR_{NLI}$, so there is no 1/3-rule at play here to attenuate the error.

IV-E Non-Gaussian formats and inter-span coherence

So far we have focused on the discussion of the PCFM performance in the context of Gaussian-shaped modulation.

As debated previously, the PCFM, like CFM6 and essentially all CFMs in the literature, is derived from the GN-model, which is most accurate when Gaussian-shaped modulation is used. One could use PCFM ‘as is’ to model QAM systems, but there would be some error. Such error is relatively modest in links with high dispersion, long spans and high symbol rates, but it is present and it is the reason why the EGN-model was developed [27], [15].

To achieve accurate format-dependent, EGN-like prediction with a GN-model derived CFM, we developed in the past suitable corrections which account for the effect of different formats, as well as the coherence of NLI among multiple spans (inter-span coherence), please see [17] and the experimental validations [10], [28]. We are also in the process of developing and training an improved set of corrections, in the quest for greater accuracy, which we leave for a forthcoming paper. Due also to space constraints, we do not pursue this topic further here.

V Downloadable software and developments

We have implemented in software the SCI+XCI PCFM model described in this paper, with SPP fitting polynomial degree $N_{p}\!=\!9$. It includes the corrections mentioned above for multi-span and multi-format (EGN-model) support.

The software is quite flexible and allows to estimate the $\rm GSNR_{NLI}$ and total $\rm GSNR$ of rather general WDM UWB systems, where each span can be different from any other. It also handles forward and backward-pumped Raman amplification, as well as ISRS. ASE is computed for lumped amplification as well as Raman. Fiber data is fully customizable, span by span. The WDM comb does not need to be uniform: channels can have all different: spacing, symbol rate, modulation format and launch power. The launch powers can change span by span.

Currently, the main file and parameter file are in editable Matlab code, whereas the model is provided as executable files for Windows and Linux [30]. It is available for download under the Creative Commons 4.0 License. Editing the main file and parameter file provides full system customization flexibility.

We will be upgrading it to provide more pre-programmed examples and to add further features. Future releases will include Double Rayleigh Backscattering, XCI without the island-stretching approximation, as well as MCI with rectangular islands, to support low dispersion and multi-subcarrier systems.

VI Conclusion

We have presented a closed-form GN/EGN model, called Polynomial-CFM or PCFM, based on a new analytical approach consisting of expressing the spatial power profile, or SPP, of each channel as a polynomial. Once this is done, remarkably, closed-form analytical solutions of the core integrals of the GN-model can be found with few approximations. This leads to improved CFM reliability and generality. The polynomial SPP representation allows to model quite general system conditions: we have shown results and examples of the accuracy and generality of the PCFM, for UWB systems with ISRS and Raman amplification, as well as lumped loss along the span.

The polynomial SPP approach promises to allow tackling even more general system configurations, that were difficult to deal with by previous CFMs, such as low symbol rates, multi-subcarrier transmission, very low dispersion and others.

We provide a software implementation of the model described in this paper, available for download and fully customizable.