A Closed-Form Expression for the Gaussian Noise Model in the Presence of Inter-Channel Stimulated Raman Scattering Extended for Arbitrary Loss and Fibre Length

A closed-form formula for the nonlinear interference (NLI) estimation using the Gaussian noise (GN) model in the presence of inter-channel stimulated Raman scattering (ISRS) is derived. The formula enables accurate estimation of the NLI evolution along any portion of the fibre span together with arbitrary values of optical fibre losses. The formula also accounts for wavelength-dependent fibre parameters, variable modulation formats and launch power profiles. The formula is suitable for ultra-wideband (UWB) optical transmission systems and its accuracy is assessed for a system with 20 THz optical bandwidth over the entire S-, C-, and L- band through comparison with numerical integration of the ISRS GN model and split-step Fourier method (SSFM) simulations in point-to-point transmission and inline NLI estimation scenarios.

. Although bandwidth expansion can lead to higher data throughputs, several challenges arise in modelling UWB transmission. Among these challenges, the wavelength-dependent optical fibre parameters [5], [6], together with the ISRS effect [7], [8], [9], play a significant role. In contrast to conventional C-band systems, they must be taken into account. The combination of these effects, together with their Kerr-induced nonlinear interaction leads to additional signal degradation and variations in performance between channels. Some studies and strategies to compensate for these effects have recently been proposed in [10], [11], [12], [13].
Associated with this, research in adaptive network planning tools aims to introduce intelligence in the network and deliver capacity when and where it is needed [14]. This is an essential step to achieve efficient resource-utilization [15] and to build a self-controlled network infrastructure. To cope with this, one requirement is to bring physical layer awareness to the control plane level [16] enabling it to account for inline signal impairments, to predict failures, and to avoid wasting resources. To achieve this, an efficient, fast and accurate model to estimate NLI at any portion of the optical fibre link is essential.
To enable real-time prediction of the UWB system performance, formulations in closed form are needed. These formulations must offer a fast, yet accurate, evaluation of the network characteristics, so that they can be widely used for network optimisation purposes [1], [2]. The closed-form expressions derived using the ISRS Gaussian noise (GN) model [17] are a starting point due to their simplicity and efficiency in estimating NLI in UWB systems, and numerous closed-form equations have been developed to date [18], [19], [20], [21].
However, these studies can provide models for NLI estimation for a subset of scenarios only. The proposed formulas do not account for the cases of short-span lengths and extremely low losses, due to the approximations made in their derivation. The first case is essential for estimating the NLI in every portion of the fibre spans, while the second case is essential when modelling, for instance, Raman amplified links, in which the effective attenuation is much lower than the intrinsic fibre attenuation.
In our recent work published in [22], these limitations were overcome, allowing the derivation of a closed-form expression capable of accurately estimating the NLI in the presence of ISRS for any fibre span length and for fibres with extremely low losses (∼ 0.04 dB/km). This was enabled by removing one of the main approximations used in deriving the formulas in [18], [19], [23]. The proposed closed-form expression in [22] accounts for all modulation formats, wavelength-dependent attenuation and dispersion, and its accuracy was compared with the ISRS GN model in integral form [17], [19] in this same work. Note that the closed-form formula derived in [24], [25] actually account for short span lengths and extremely low losses but do not include the ISRS effect, and hence are not suitable for UWB modelling in conventional optical fibres, which is rapidly gaining interest now.
This paper presents all the assumptions and the mathematical derivations used to obtain the closed-form expression in [22]. We have also simplified this closed-form expression, validated its accuracy with SSFM simulations, and present a discussion on the validity range of this closed-form formula compared with those in [17], [19], [23]. Finally, the proposed closed-form expression is applied to estimate the evolution of NLI along a fibre link, demonstrating one of the multiple applications of the proposed closed-form expression.

II. THE DERIVATION OF THE CLOSED-FORM EXPRESSION FOR NLI-INDUCED SNR
This section describes the closed-form expression used to analytically estimate the NLI. The integral model is presented and the steps used to derive a closed-form expression of this model are detailed.

A. Preliminaries
After coherent detection and electronic dispersion compensation, the total received SNR for the i-th WDM channel (SNR i ) after n spans can be estimated as where SNR TRX , SNR ASE and SNR NLI are the SNR from the transceiver subsystem, the amplified spontaneous emission (ASE) from the optical amplifiers used to compensate for the fibre loss and the accumulated NLI, respectively. n is the number of spans, i is the channel of interest (COI), P i is its launch power, κ i = 1/SNR TRX i , P ASE i is the ASE noise power at the i-th channel frequency, and P NLI i = η n (f i )P 3 i is the NLI noise power after n spans. All three impairments are assumed to be statistically independent of one another. In this paper, we will focus on the SNR NLI calculation. The impact of ASE noise power and TRX noise power in the total SNR are extensively discussed in [2], [13].
To calculated the power spectral density (PSD) P NLI = η n (f i )P 3 i , the ISRS GN model approach is considered. This model is an extension of the GN model [26] accounting for the ISRS effect and was proposed in [17], [19]. This model also accounts for the modulation dependence of the NLI in the input symbol distribution [27], [28], [29]. This dependence is accounted for by calculating one of the NLI correction terms in [29] in the presence ISRS.

B. The ISRS GN Model in Integral Form
In this section, the integral expressions used to derive the proposed closed-form expression are presented. The nonlinear coefficient obtained at the end of the n-th span, η n (f i ), can be written as [29,Eq. (1)], where η GN,n (f i ) and η corr,n (f i ) are respectively the nonlinear coefficient contributions accounting for Gaussian modulated symbols [17] and its correction term accounting for the dependence of the NLI on the modulation format [19]. For Gaussian modulated signals, the correction term η corr,n (f i ) vanishes and one obtains the model in [17].
Following the assumptions described in [18], η GN,n (f i ) in (1) can be approximated as [18, (5)] where η SPM j (f i ) is the SPM contribution and η XPM j (f i ) is the total XPM contribution to the NLI, both generated in the j-th span. P i,j is the power of channel i launched into the j-th span, is the coherent factor [26,Eq. (22)]. For notation convenience, the j dependence of the SPM and XPM contribution is suppressed below.
The XPM contribution in (3), η XPM (f i ), is obtained by summing over all COI-interfering pairs present in the transmitted signal, i.e., where N ch is the number of WDM channels and η and where Π(x) denotes the rectangular function and B k is the bandwidth of the channel k. μ(f 1 , f 2 , f i ) is the so-called link function or FWM efficiency [26], which is given by [17,Eq. (4)] Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
, and ρ(z, f i ) is the normalized signal power profile. β 2 is the group velocity dispersion (GVD) parameter, β 3 is the linear slope of the GVD parameter. Now, the correction term contribution to the NLI in (8) η corr,n (f i ) is considered. Similar to (4) and following the assumptions described in [19], η corr,n (f i ) is obtained by summing over all COI-interfering-channel pairs present in the transmitted signal, The correction term η corr,n (f i ) is the XPM contribution of a single interfering channel k on channel i, which is given by [19,Eq. (4)] where L and L is the span length. Φ stands for the excess kurtosis of the given constellation, providing statistical characteristics of the signal, and reflecting how the constellation deviates from the Gaussian one. As shown in detail in [19], (9) can be accurately approximate as corr,1 (f i ) is a correction term originating in the first span and η (k) corr,a (f i ) is an asymptotic correction term originating in the limit of large span number. η corr,a (f i ) are given respectively by [19,Eqs. (9), (12)] and

C. The Closed-Form Expression
This section is divided into two parts. Firstly, we present the semi-analytical solution of the Raman differential equations, which is used to represent the signal power profile. Secondly, we present the novelty of this work, i.e., the closed-form expressions of (5), (6), (11) and (12), which are then used to calculate (2) in closed form.
1) The first step in deriving a closed-form expression for (2) is to derive a closed-form expression for the link function in (7). To that end, the normalised power evolution along the fibre ρ(z, f i ) is considered as the semi-analytical solution of the Raman differential equations [30], [31], approximated by a Taylor series to the first order as in [18], which is given by [18,Eq. (17)], [23, Eq. (2)]: where , P tot is the total launch power, α i and α i model the fibre loss, C r,i is the slope of the Raman gain spectrum. Eq. (13) is obtained by considering several assumptions, we describe them in the following. Firstly, as by [30], [31], a constant attenuation profile, a triangular approximation of the Raman gain spectrum and the approximation f k f i ≈ 1 are assumed. Afterwards, in the equations derived in [30], [31], a spectrally uniform launch power profile is assumed, and a firstorder Taylor expansion is used in [18] leading to (13) (see [18], Section II-E for further discussion of these assumptions).
To overcome the restrictive assumptions described above, a fitting strategy is used, whereas in (13), as in [18], two different loss coefficients are considered (α i andα i ) and together with C r,i are treated as channel-dependent fitting parameters and matched using a fitting algorithm to reproduce the true power profile, which is obtained by numerically solving the Raman differential equations [31] using the Raman profile shown in Fig. 1(b). Also, the utilisation of two separate loss coefficients (α i andα i ) enables an increase in the dimension of optimisation space.
Because of the fitting optimisation routine, we call (13) as a semi-analytical solution. Moreover, in this work, the fitting algorithm was found to overcome all the assumptions used to derive (13) making this equation valid for a wide range of simulation scenarios, such as nonuniform launch power profiles, wavelength-dependent attenuation, non-triangular Raman gain spectrum, etc. Note that, in its current form, for each new link configuration, the fitting optimisation needs to be done again, however, scaling rules of such fitting coefficients with the physical parameters might also be exploited.
2) We now introduce the novel equations of this work w.r.t. the ones in [18], [19], [23]. The novelty is a new approximation, shown in (30), which enables to extend the equations in [18], [19], [23] to short-span lengths and arbitrary fibre loss. This approximation is shown in Appendix A and discussed in Appendix D and is reflected in all mathematical derivations to obtain a closed-form expression for (5), (6), (11) and (12), as shown below. LetT By assuming the normalised power evolution along the fibre ρ(z, f i ) as (13), the link function in (7) can be approximated in closed form as whereα l,i and κ l,i are given respectively bỹ and The proof of (14) is given in Appendix A. The coefficientsα l ,i and κ l ,i are respectively the same as those in (15) and (16) with the indices l replaced by l . The same is valid for the variable α l ,i . The next step is to derive closed-form expressions for the XPM and SPM NLI contributions given by (5) and (6), respectively. Using (14) as an analytical solution of the link function, a closed-form expression for the XPM and SPM are given respectively by and The proof of (17) and (18) are given respectively in Appendix B and C. The final step is to derive closed-form expressions for the correction terms, which account for the dependence of the NLI on the modulation format, i.e, a closed-form expression for the terms given in (11) and (12). For (11), a similar integral was already solved in Appendix B for the XPM contribution, the solution is given by (35). Thus, (11) is given in closed form as For (12), it is enough to note that μ( (14), i.e. when substituting identical arguments. Thus, (12) can be written in closed form as We have now presented the complete set of equations to calculate SNR NLI in (1), i.e. (17), (18), (19) and (20). In order to write the complete equation for SNR NLI in (1) and simplify the notation, we will assume that the optical link under study is made up of identical spans in terms of fibre parameters (the homogeneous link assumption). This is equivalent of assuming that η SPM (f i ), η XPM (f i ), and η (k) corr,a (f i ) are independent of the fibre span j, and n j=1 (P i,j /P i ) 2 = n in (3) and (10). Under this condition, the SNR NLI can be written as (21) shown at the bottom of the next page, where in this equationñ = 0 for a single span orñ = n otherwise. Also, we include the indices i and k in all the variables to explicitly show their channel dependence. Moreover, if the homogeneous and transparent link assumption is removed, those variables will also be span-dependent. Finally, note that, in the limit α l,i L → ∞, (21) converges to that in [23].

III. RESULTS
This section describes the numerical validation of the closedform expression shown in (21). A comparison with previously reported closed-form expressions and the application of the new formula in a system scenario are also carried out.

A. Transmission Setup
The baseline transmission system, over which the derived expressions are validated, consists of a WDM transmission with N ch = 181 channels spaced by 100 GHz and centred at 1540 nm. Each channel was modulated at the symbol rate of 96 GBd. This resulted in a total bandwidth of 20 THz (158 nm), ranging from 1470 nm to 1615 nm, corresponding to the transmission over the S-(1470 nm-1530 nm), C-(1530 nm-1565 nm) and L-(1565 nm-1615 nm) bands. Spectral gaps of 10 nm and 5 nm were assumed between the S-/C-and C-/L-bands, respectively. The channels were transmitted over 5 spans using a single-mode fibre (SMF) where the span length is varied as described in the next sections. It is assumed that each amplifier fully compensates for the span losses (the transparent link assumption). A spectrally uniform input launch power profile was used, with each channel having a launch power of 1 dBm. Realistic wavelengthdependent attenuation profile and Raman gain spectrum were used, these profiles being measured for an ultra-low-loss (ULL), ITU-T G.652 fibre, and shown in Fig. 1(a) and (b), respectively. Dispersion and nonlinearity parameters were D = 16.5 ps nm·km , S = 0.067 ps nm 2 ·km and γ = 1.03 1 W·km . To verify the accuracy of the proposed-closed form expression, a variety of span lengths and losses were also considered; these values are described in detail in Section III-B. Additionally, Gaussian modulated and 64-QAM symbols were also used. For the latter, this is achieved by setting the excess kurtosis in (21) to Φ = −0.6190, against Φ = 0 for the former.

B. Numerical Validation
The transmission system performance estimation using the proposed closed-form formula, i.e., (21), was carried out for two different scenarios using the transmission setup described in Section III-A. The scenarios were chosen to assess the formula for short-span lengths and low losses. For the first scenario, a variety of span lengths were chosen. These results are shown in Fig. 2(a) for Gaussian constellations and in Fig. 2(b) for 64-QAM constellations. In the second scenario the span length was fixed to a value of 80 km and different spectrally uniform loss profiles were used. These results are also shown for Gaussian constellations in Fig. 3(a) and for 64-QAM constellations in Fig. 3(b).
The interaction between fibre attenuation, dispersion and normalised ISRS-power evolution profile, leads to the SNR NLI profile as shown in Figs. 2 and 3. The high dispersion towards the L-band reduces the NLI for the long-wavelength channels. This reduction however is counter-balanced by the ISRS-transferred power from short to long wavelength channels, increasing the NLI for these channels, and thus reducing the SNR NLI . In the case of Fig. 2, the effect of ISRS is increasingly seen as the span length increases. This is because the longer the span length, the greater the power which is transferred due to ISRS effect. In the case of Fig. 3, the same trend is observed among the curves, with a larger tilt due to the ISRS in the case of links with lower losses.
To verify the accuracy of the proposed closed-form expression, (21) was compared with the ISRS GN model in integral form, for Gaussian constellation [17] and with SSFM simulations for both Gaussian and 64-QAM constellations. For the SSFM, a local-error method [32] was used to ensure optimal step size with global local error δ G = 10 −10 . This value was found to be sufficient throughout all simulations by comparing simulations with smaller δ G values and observing insignificant changes in the signal output. The simulations were performed with each channel having 2 16 random symbols. Due to the relatively short symbol sequence, each result represents an average of eight simulations. The results of these validations are also shown in Figs. 2 and 3.
Additionally, for all the scenarios described by each figure, the maximum errors among all curves were computed. For the scenario shown in Fig. 2(a), in which the span length is varied, using Gaussian symbols, the closed-form formula in (21) shows maximum errors of 0.93 dB and 0.58 dB, from integral model and SSFM simulations respectively. For Fig. 2(b), using 64-QAM symbols 1.2 dB maximum error is obtained between (21) and SSFM. The same analysis is carried out for the scenario where the loss profile is varied. In that case, using Gaussian symbols, Fig. 3(a) shows maximum errors of 1.05 dB and 1.48 dB, respectively, while for 64-QAM symbols Fig. 3(b) shows maximum error of 1.61 dB. Moreover, a mathematical justification of the validity of the proposed closed-form expression can be found in Appendix D. Note that, as this closed-form expression is an extension of the ones in [18], [19],  its accuracy for different values of symbol rates for the different NLI contributions, namely, SPM and XPM can be found in [18], Figs. 3, 4.

C. Comparison to Previously Reported Closed-Form Expressions
In this section, we compare the accuracy of the closed-form expression proposed in this paper, i.e., (21), with that of the closed-form expressions reported in [18], [19], [23]. To that end, the simulation scenario was varied in two different ways: (a) the span length was swept from 1 km to 80 km and (b) the span length was fixed at 80 km and a spectrally-uniform loss profile ranging from 0.02 dB/km to 0.2 dB/km was considered. The results were obtained considering Gaussian constellations and transmission over 5 spans. Fig. 4 shows the maximum per-channel SNR NLI difference, i.e., the maximum per-channel error in terms of SNR NLI between the integral ISRS GN model and the proposed closed-form formula presented in this paper. For comparison, the same analysis is also carried out with the closed-form expression reported in [23]. As shown in Fig. 4, the new closed-form formula proposed in this paper can accurately account for any span length and fibre loss; among all the scenarios considered in Fig. 4, maximum errors of 0.93 dB and 1.27 dB were found respectively when considering different span lengths and losses. A mathematical justification of the inaccuracy of the closed-form expression reported in [18], [19], [23] and the validity range of the one shown in (21) are given in Appendix D.

D. NLI Evolution During Propagation in the Fibre
One of the motivations and importance of the closed-form formula derived in (21) is the possibility to perform an accurate estimate of the NLI in every portion of the fibre link, enabling the calculation of the NLI evolution during propagation in the fibre. The importance of such calculation was mentioned in Section I. In order to illustrate it, we consider the transmission Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.  setup described in Section III-A and a 400 km optical fibre link formed of 5 × 80 km spans. Equation (21) is then applied to estimate the NLI for the first 10 km of each span and at the end of the link. Fig. 5 shows the nonlinear coefficient η n , defined by (1), as a function of wavelength for different distances using (21) (continuous lines). The results were also matched using SSFM simulations (dotted lines); among all the curves a maximum per-channel error of 0.58 dB was found after a propagation distance of 400 km. Note that, the discussion of these results is similar to the ones described in Section III-B and is omitted in this section.

IV. CONCLUSION
A closed-form formula that can accurately evaluate the nonlinear interference (NLI) in the presence of ISRS at any step of the fibre span and in extremely low loss regimes (∼ 0.04 dB/km) was proposed. The formula was applied in modelling an S+C+L band (20 THz) transmission system and its accuracy was verified through comparisons with results obtained using integral model and split-step Fourier method simulations. Using the proposed closed-form formula, the NLI can be calculated in a few microseconds, enabling rapid performance evaluations (e.g., SNR, maximum reach, optimum launch power) in ultrawideband transmission systems.
The proposed formula enables accurate inline nonlinear interference estimation of any portion of the optical fibre link, representing an essential step towards the development of intelligent and dynamic optical fibre networks. The latter, together with the computational speed of the proposed closed-form formula, will enable effective network planning tools allowing an online assessment of the data rates, modulation formats, the number of channels and launch power profile, given the fibre and amplifier characteristics and allocated lightpaths.

IV. DATA AVAILABILITY STATEMENT
The data that supports the figures in this paper are available from the UCL Research Data Repository (DOI: 10.5522/04/21630251), hosted by FigShare.

APPENDIX A DERIVATION OF THE LINK FUNCTION
This section describes the derivation of (14). Let (13) can then be written as Inserting (22) in (7) yields The term x(ζ) can be written as (25) can be conveniently rewritten in terms of a summation using identity (43), which will facilitate all the mathematical derivations, Now, inserting (26) in (23) Solving the integral in (27) yields (28) Now, let us define α l,i = α i + lα i . Eq (28) can then be written as (29) Despite (29) being in closed form, it needs to be further integrated in frequency to obtain the XPM and SPM contributions to the NLI as shown in (5), (6), (11) and (12). In order to solve these integrals in closed-form, the approach in [24] is used, i.e., we approximate the fraction with exponential terms in (29) as where κ l,i andα l,i are chosen such that the first-order Taylor approximations of both the left and the right side of (30) around the variable φ = 0 become equal. This yields (15) and (16). Inserting the approximation shown in (30) into (29) and using identity (47), yields Finally, performing the multiplication in (31) together with the identity (45) yields (14), concluding the proof.

APPENDIX B DERIVATION OF THE XPM CONTRIBUTION
This section presents the derivation of (17). We start by approximating the phase mismatch term in (7). For the XPM contribution, let Δf = f k − f i be the frequency separation between channels k and i. Assuming that frequency separation is much larger than half of the bandwidth of channel k, i.e, |Δf | B k 2 , we can make the assumption that f 2 + Δf ≈ Δf . Also, we assume that the dispersion slope β 3 is constant over the channel bandwidth. Thus, the phase mismatch term can be approximate as [18, (15)], The channels most impacted by this approximation are those near the COI.
The error relative to this approximation is given by [18], 25. Now, we consider (5). For notation brevity, we will omit the factor 32 27 is neglected -this is equivalent to approximating the integration domain of the GN model to a rectangle [26]. Inserting (14) in (5) η (k) .
Because of the approximation in (32), φ no longer depends on f 2 . Thus, (33) can be written as .
The integral in (34) can be solved using identity (46), yielding By inserting again the pre-factor 32 27 (17) is obtained, concluding the proof.

APPENDIX C DERIVATION OF THE SPM CONTRIBUTION
This section presents the derivation of (18). We start by approximating the phase mismatch term. We assume that the dispersion slope β 3 is constant over the channel bandwidth. Thus, the phase mismatch term can be approximated as Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
with φ i = −4π 2 (β 2 + 2πβ 3 f i ). Now, using (6) together with (5) with k = i, and omitting the pre-factor of 16 27 Note that, similarly to Appendix B, the term Π( f 1 +f 2 B i ) is neglected. The integral in (37) can be rewritten in polar coordinates (r, ϕ) as where the relations f 1 = r cos (ϕ/2), f 2 = r sin (ϕ/2) and sin (ϕ/2) cos (ϕ/2) = sin (ϕ) 2 were used. Also the integration domain of (6) was approximated by a circular domain such that the area of both domains are equal [18], Fig. 3. This yields the variation of the radius in the outer integral as shown in (38). The inner integral in (38) can be solved using identity (47) This integral can be rewritten as:

APPENDIX D VALIDITY RANGE
This section shows the validity range of the closed-form formula proposed in (21) and the limitations of that derived in [18], [19], [23]. To that end we must consider (29). In the work [18], [19], [23], the following approximation is used where it is easily seen that its accuracy relies on the condition α l,i L → ∞. The breaking of this condition is the source of the inaccuracy shown in Fig. 4 for the closed-form expression published in [18], [19], [23].
On the other hand, (21) relies on the approximation shown in (30) in Appendix A, in which a Taylor expansion is performed around φ = 0. The following approach is inaccurate when α l,i L → 0. This is because in this condition, the oscillator function in the numerator 1 − e jφL ≈ 1, as required by the Taylor expansion around φ = 0. Thus, in order to satisfy the requirement 1 − e jφL ≈ 1, the condition α l,i L 0 should be satisfied.
Finally, by comparing the requirement of the accuracy of the closed-form formulas published in [18], [19], [23] (α l,i L → ∞) with that of (21) (α l,i L 0), it is noted that this last requirement is much less restrictive than that of α l,i L → ∞, justifying the accuracy of (21) for all the simulations and scenarios considered in this paper. Additionally, in the limit α l,i L → ∞, (21) converges to that in [23].