A Generalized Raman Scattering Model for Real-Time SNR Estimation of Multi-Band Systems

We investigate the impact of stimulated Raman scattering (SRS) in multi-band transmissions. We first derive and test against numerical results an extended Christodoulides-Zirngibl (ECZ) SRS model based on a triangular approximation of the Raman gain. Then, we exploit the ECZ SRS model in the context of the Gaussian noise (GN) model. We generalize some known closed-form expressions of the self- and cross-phase modulation variances, as well as the amplified spontaneous emission noise variance, over multi-band systems of bandwidth wider than the SRS bandwidth, i.e., beyond the C+L band. Finally, the proposed theoretical model is used to estimate the performance of multi-band systems in terms of signal-to-noise ratio (SNR) and achievable information rate (AIR).

still be used to amplify the L band (1565-1625 nm) at the expense of a higher noise figure, the amplification of the S (1460-1530 nm) and the E (1360-1460 nm) bands requires different rare-earth elements, such as Thulium and Neodymium [2]. Other possible strategies are represented by distributed Raman amplifiers (DRA) [3], [4] and semiconductor optical amplifiers (SOAs) [5], which inherently offer wide amplification bandwidths.
However, the response of the optical fiber to such wide bandwidths becomes non-trivial. First of all, the fiber parameters such as attenuation and dispersion have a non-negligible dependency on frequency, such that WDM channels belonging to different bands undergo linear effects with different strengths [6]. While the linear effects can be compensated, their frequency dependence yields complex nonlinear interactions among spectral components through the fiber Kerr effect. The different bands interact nonlinearly also due to stimulated Raman scattering (SRS), through which low frequencies are amplified at the expense of higher-frequency components [6].
The design of multi-band systems accounting for the impact of all such aspects on the quality of transmission (QoT) represents a challenging task. An accurate QoT estimation can be achieved by numerically simulating the optical system by solving the optical field propagation along the fiber through split-step Fourier method (SSFM) simulations. However, the computational effort and time required by such a method escalate very quickly with the signal bandwidth [7].
In this context, semi-analytical perturbative models offer a significant complexity reduction. Of particular interest is the family of models that treats the nonlinear impairments as additive nonlinear interference (NLI). Among these models, the Gaussian noise (GN) model [8] stands out for its simplicity providing an estimate of the NLI power in the presence of Gaussian distributed symbols. The accuracy of the estimation is further improved by accounting for the dependency of the NLI on the modulation format of the transmitted symbols, as done by the enhanced GN model [9]. Ultra-wideband extensions of the GN model to include SRS were proposed in [10], [11], [12], and later exploited to analyze multi-band optical systems in [1], [13], [14], [15]. In particular, the work in [12] exploited the Christodoulides-Zirngibl formula of the signal power evolution along distance in the presence of SRS [16], [17], which was derived under the assumption of linear Raman gain vs frequency.
Along the same lines, the modulation format was included in [18].
The search for closed-form expressions for the SRS-aware NLI power has gained interest since it can bring the ultimate complexity reduction by allowing real-time estimations. Notably, the work in [19] proposed closed-form expressions of the semi-analytical estimate of the variance of self-phase modulation (SPM) and cross-phase modulation (XPM) of [12], based on [16], [17]. However, the linear approximation of the Raman gain of [16], [17] is justified when the WDM bandwidth is below ≈15 THz, i.e., up to the C+L band, while it over-estimates the Raman gain for larger bandwidths. To overcome this limitation, the works in [20], [21], [22], [23] fitted the numerical solution of the power evolution along distance to extend the GN closed-form formulas.
In this work, which is an extended version of [24], we propose to leverage a truncated triangular approximation of the Raman gain to extend the Christodoulides-Zirngibl (CZ) SRS formula for the signal power evolution derived in [16], [17]. We show that the proposed extended CZ (ECZ) SRS model can be used to extend the GN model formulas in [19] to multi-band transmissions even beyond the C+L band and we compare it against a numerical implementation of the GN model. Moreover, the SRS formula can be used to estimate even the SRS-inflated amplified spontaneous emission (ASE) power, thus yielding an approximated real-time estimation of the signal-to-noise ratio (SNR).
The remainder of this paper is organized as follows. In Section II, we tackle the topic of signal power evolution and derive the ECZ SRS model based on the truncated triangular Raman gain approximation. In Section III we generalize the SPM and XPM variance formulas based on the novel ECZ SRS model. Section IV focuses on the estimation of QoT metrics. Finally, in Section V we draw our conclusions.

II. POWER EVOLUTION WITH SRS
The evolution along the fiber-optic propagation distance z of the ith WDM channel power P i is governed by the following ordinary differential equation (ODE) [6], for i = 1, . . . , N ch : where λ i is the carrier wavelength of the ith channel, f i is its low-pass frequency in a reference system centered at the WDM central frequency, N ch is the number of channels, α is the fiber loss coefficient, andg R the polarization-averaged Raman gain coefficient. The first term on the right-hand side of (1) accounts for the depletion of channel i to amplify lower frequency channels, while the second term accounts for the power acquired from higher frequency channels. Fig. 1(a) shows the fiber loss coefficient, as a function of the wavelength λ, obtained through a quadratic polynomial fitting,  [20] vs. wavelength. (b) Raman gain coefficientg R profile (interpolated from [6]) vs frequency shift, with the two approximations investigated in this work, i.e., linear and triangular. in a dB scale, of the experimental attenuation profile available in [20]: with λ 0 = 1550 nm and coefficients: α 2 = 3.7685 × 10 −6 [dB/(km · nm 2 )], α 1 = −7.3764 × 10 −5 [dB/(km · nm)], and α 0 = 0.162 [dB/km]. We observe a negligible variation over the C or the C+L band since the maximum excursion is 0.016 dB/km. On the other hand, such an excursion becomes significant for wider bandwidths, reaching 0.15 dB/km over the E+S+C+L bands. Fig. 1(b) shows the Raman gain coefficient [6], as a function of the frequency shift Δf = f i − f j between channel i and j.
Since the power evolution in the general form of (1) does not have an analytical solution, it needs to be solved numerically for each WDM channel. We next review the approximated analytical solution of (1) proposed in [16], [17], which is valid for frequency shifts within ≈ 15 THz, and then we propose a generalization for wider bandwidths.

A. The Christodoulides-Zirngibl SRS Model
An analytical solution of (1) has been derived in [16], [17] under the following assumptions: i) a linear dependency of g R (Δf ) on the frequency shift, namelyg R (Δf ) ≈ C r Δf , ii) a frequency-independent fiber attenuation coefficientᾱ, and iii) λ i λ j ≈ 1. The power evolution equation then simplifies as follows: where the first term in the right-hand side accounts both for amplification and depletion, depending on the sign of f j − f i . The coefficient C r represents the slope of the linear approximation of the Raman gain coefficient, as sketched in Fig. 1, and can be found through a linear regression on the experimental data.
Along the lines of [17], the system of coupled ODEs in (2) can be written as a single equation in terms of the WDM power spectral density (PSD) G(z, f ) instead of channel powers, namely: where the summation over the channels was replaced by integration in frequency. A key property of (3) is that the total WDM power at any coordinate z, defined as P t (z) ∞ −∞ G(z, μ)dμ, is set by the fiber loss only. In fact, by integrating both sides of (3) with respect to frequency, the first term on the right-hand side disappears resulting in P t (z) = P t e −ᾱz , where P t ≡ P t (0) is the WDM total power injected in the fiber. Therefore, dividing (3) by G(z, f ) and taking the derivative with respect to frequency yields: Following [17], (4) can be integrated first with respect to distance and then to frequency, thus obtaining: where L eff (z) = (1 − e −αz )/α is the fiber effective length at coordinate z, and F (z) an integration constant. By imposing the conservation of the total power we get: From (6) it can be seen that the SRS induces a frequencydependent tilt in the signal power profile without affecting the total power, which decreases only through the fiber loss. Fig. 1(b) suggests that a linear approximation of the Raman gain is justified when the maximum frequency shift is less than ≈15 THz, as the Raman gain coefficient quickly vanishes at wider spacings. For larger bandwidths, the Raman gain coefficient is better described through the so-called triangular approximation [25], where the gain is treated as a linear function up to a cut-off shift Δf co and forced to zero at higher frequencies, as sketched in Fig. 1(b). Under this approximation, a generic frequency f interacts with the neighboring ones only within an effective window, experiencing amplification/depletion from the higher/smaller frequency channels, respectively. We define the effective window as the intersection between the WDM bandwidth and the two-sided Raman triangular window of bandwidth 2Δf co . The idea is sketched in Fig. 2, where ν 1 and ν 2 represent the edges of such an effective window, while f m and f M are the minimum and maximum frequencies in the WDM bandwidth, respectively. It is worth noting that the special case 1 coincides with the scenario analyzed in the previous section since the whole WDM bandwidth B t is covered by the Raman window.

B. The Extended CZ SRS Model
In this more general context, (3) can be generalized to: where the frequency integral now spans only the frequencies within the effective window, whose limits are: The difference with respect to the case analyzed in the previous section is thus the frequency dependency of the integral limits in (7). We thus exploit the Leibniz integral rule obtaining: where we introduced the definition of integrated power in the effective window: which coincides with the total power at coordinate z only when the Raman window is wider than the WDM bandwidth. Since G(z, ν 1,2 ) = 0 only for f m ≤ ν 1,2 ≤ f M , we can remove the min/max of (8)-(9) by ν 1,2 ≡ f ± Δf co . This allows us to express (10) in the compact form: where we introduced the following definition: It is worth noting that P I,T (z; f, Δf co ) coincides with the PSD integrated in the 2Δf co Raman window through the trapezoidal rule. Fig. 2 shows the relation between P I and P I,T in the four cases under investigation. In particular, the two quantities are equal when the Raman window is narrower than the WDM bandwidth. Equation (12) thus generalizes (4), and coincides with it in the special case 1 depicted in Fig. 2, where the integrated power coincides with the total power and P I,T (z; f, Δf co ) = 0. Equation (12) does not have an analytical solution. We thus proceed by introducing the following key approximations: f m and f M are the minimum and maximum WDM frequencies, respectively. The frequencies ν 1 and ν 2 represent the edges of the effective window, i..e, the intersection between the WDM spectrum and the Raman triangular window, which is identified by the shaded region. Top part of the figure: relationship between the integrated powers P I and P I,T defined in (11) and (13), respectively.
which assume that the power is undepleted by SRS within a window (effective or Raman, respectively). Such an assumption holds when the window coincides with the WDM bandwidth and the attenuation coefficient is constant with respect to frequency, as in Section II-A. Along the lines of Section II-A, we then integrate with respect to distance and frequency, obtaining: which generalizes (5) through the Raman shaping profile r(f ). Such a term results from the evaluation of the following indefinite integral: we show in the Appendix that the Raman shaping profile can be expressed as: for the four cases sketched in Fig. 2. In the simplest scenario, case 1 , the Raman window encompasses the whole WDM bandwidth B t . Therefore, P I ≡ P t and P I,T = 0 yielding r(f ) = P t f as in the original theory [16], [17]. On the other hand, if the Raman window covers only an inner portion of the WDM bandwidth, r(f ) balances out to zero (case 2 ) since P I = P I,T . For all the other scenarios, the term r(f ) takes one of the expressions labeled with 3 or 4 in (18), depending on the position of the frequency of interest within the WDM. The integration constant F (z) is found by imposing the conservation of the total power, similarly to the previous section, and thus yielding the following solution: which is formally identical to (6) with the replacement of P t f → r(f ) in the exponential terms. Since r(f ) is slowly varying over a channel bandwidth, we can substitute G with a train of Dirac's delta obtaining the following equation for P i : where P is the launch power per channel and we also arbitrarily chose as the common fiber attenuationᾱ within a window its value at the channel under test, i.e., as this choice provides the best fit with simulations.

C. Numerical Results
We test the proposed analytical solution in (20) against a numerical solution of the ODEs in (1), which is based on the true attenuation profile and Raman gain coefficient in Fig. 1, through the MATLAB function ode45. We chose the slope C r and the cutoff frequency shift Δf co by minimizing the least-squares error over the WDM bandwidth between the solution of (1) with the trueg R and the triangular approximation.
In a first test, we considered the transmission of 64 Gbaud channels spaced 75 GHz, covering either the C+L (11.4 THz), or the S+C+L (20.8 THz), or the E+S+C+L (35.9 THz) band without any guard band intervals. The total power was set to 21 dBm, while the link consisted of L = 100 km of an SMF. For the generic ith WDM channel, we define the span SRS gain as the ratio between the channel power at fiber-end and its value without the SRS, namely: Fig. 3 shows the SRS gain computed both with the numerical solution of the ODE (markers) and by using the proposed formula (solid line) in (20). For completeness, we show with a dashed line the CZ SRS solution [16], [17] in (6). For the sake of comparison, in (6) we chose the constant attenuation coefficient as in (20). The figure shows a good agreement between the numerical results and the proposed ECZ model, with a root mean square error (RMSE) of 0.07 dB, 0.18 dB, and 0.29 dB, for the C+L, S+C+L, and E+S+C+L transmission, respectively. For the same transmissions, the CZ model (linear approximation of the Raman gain) yields an RMSE of 0.07 dB, 0.6 dB, and 2.9 dB. From the figure, we note that over the C+L band the triangular and the linear SRS approximations are overlapping. This is not surprising, since the Raman window encompasses the entire WDM band. Differences start appearing in the S+C+L case, mainly in the S and L band, where there are edge channels that fall in the cases 3 and 4 , since their Raman window does not cover the whole bandwidth. It can be seen from Fig. 3, that the SRS impact on these regions is smaller than the one predicted by [16], [17] since only a subset of the overall channels participates in it. Thanks to the triangular shape of the Raman gain adopted in the proposed model, (20) is able to capture such a phenomenon. On the other hand, in the central WDM region, the Raman window encompasses the whole WDM bandwidth (case 1 of Fig. 2) and thus the SRS gain is well predicted by the CZ formula. In the E+S+C+L band, the linear SRS approximation overestimates the Raman effect all over the WDM bandwidth by artificially increasing the gain tilt. It is worth noting that in such a case an inner portion of channels meets the conditions of case 2 in Fig. 2, which, in the proposed theory, have zero Raman shaping profile in (18). The numerical results reported in Fig. 3 are in agreement with such a prediction, exhibiting an SRS gain ≈ 0 dB in the central part of the bandwidth.
A key assumption of the proposed model is represented by the undepleted integrated power in (14)- (15). This assumption becomes critical in the presence of strong SRS, for instance when the total WDM power P t is high. We thus investigated selected cases of P t focusing on the S+C+L transmission case. The results are shown in Fig. 4 for P t = 21 dBm, 23 dBm, and 25 dBm corresponding to a channel power of −3.4 dBm, −1.4 dBm, and 1.4 dBm, respectively. We note that the ECZ SRS model is able to capture the qualitative three-region shape of the SRS gain. Nevertheless, the accuracy with respect to the Runge-Kutta results decreases in the edge regions as the power increases, with an RMSE of 0.18 dB, 0.3 dB, and 0.8 dB, respectively.
We then investigated the simple yet relevant case of signal power pre-emphasis in the form of an opposite-sign SRS gain. Namely, this power allocation aims to pre-compensate a portion of the span SRS gain by applying a pre-emphasis to an uniform launched power. As an approximation, we still use (20) with: withk the pre-emphasis factor [18].    Fig. 3 with the different power allocations reported in the inset. In particular, we fixed the total power to 21 dBm and then applied a pre-emphasis with a factork ∈ (0, 0.5, 1). The numerical results (markers) show that the pre-emphasis effectively reduces the SRS gain experienced by the signal. We then compared the results with the formula in (20) combined with (22), which is reported with solid lines. The RMSE is within 0.2 dB in all cases, justifying the use of (20) even with a non-flat power allocation.

III. CLOSED-FORM GN MODEL
In this Section, we exploit the novel ECZ SRS model in the context of the GN model. In this framework, the NLI at channel f i is treated as an additive noise [26] whose variance, after compensation of the linear effects and SRS on the signal, is [8], [19]: Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
with B i the channel bandwidth, γ(f, f 1 , f 2 ) the fiber nonlinear coefficient 1 , and η(f, f 1 , f 2 ) the link kernel (or link function) weighting the four-wave mixing interaction among signal frequencies, which can be written as: for a dispersion uncompensated link composed of N identical spans, each of length L, with lumped amplification in each span followed by a gain-flattening filter (GFF) or a dynamic gain equalizer (DGE) to restore the signal PSD. The link kernel accounts for the fiber loss and the SRS experienced by the signals along the fiber through the following function: with R(z, f ) G(z, f )/G(0, f) the loss/gain profile at coordinate z and frequency f . Equation (24) also includes the fiber chromatic dispersion through the phase matching coefficient , with β 2 the dispersion and β 3 its slope. It is worth noting that the integral in (24) accounts for the nonlinearity arising locally within the span, while the outer summation keeps track of the dispersive effects accumulated in the previous spans. The NLI variance in (23) can be evaluated by means of numerical integration in frequency and in distance. The integral over the distance can be approximated by following similar steps of [19]: where η 1 (ᾱ) is the single-span kernel without SRS [8]: while the factor S accounts for the presence of SRS through the Raman shaping profile r(f ) in (18) as follows: The approximation sign in (26) accounts for the first-order For a homogeneous link, the external summation in (24) can be closed in the so-called phased-array term [8]. It is worth noting that (26) coincides with the solution in [19] in case 1 of Fig. 2.
Since the power profile R(z, f ) can be approximated as flat within a channel band, we can generalize the closed-form expressions of [19] by: 1 Motivated by the analysis of [22] we included the frequency dependency of the fiber nonlinear coefficient, although not originally present in [8], [19].
The novelty of (29)-(30) relies in the SRSdependent term T i , which is here defined as: and generalizes the term (2α i − C r P t f i ) 2 in [19]. Moreover, in (29)-(30) we included a pre-emphasis factork by exploiting (20), (22) and the wavelength dependency of the fiber nonlinear coefficient [27]: with A eff (λ i ) the fiber effective area at the wavelength λ i and n 2 the fiber nonlinear index. The channel index discriminates between the nonlinear coefficient relative to SPM ( = i) or XPM ( = i). While the XPM variance accumulates incoherently with the number of spans, the factor ε in (29) accounts for a coherent accumulation [8], [19], [28] of the SPM variance along the distance.

A. Numerical Results
We tested the generalized closed-form expressions of the SPM and XPM variance in (29)-(30) against a fully-numerical implementation of the GN model. Such a benchmark is obtained by using the numerical solution of the ODE in (1) with the true profiles of Fig. 1 to obtain the function ρ(ζ, f, f 1 , f 2 ) for the link kernel in (24). The NLI variance in (23) is then evaluated numerically by exploiting Filon's method [18] for the integral along the distance in (24) and the Monte Carlo method [29] for the double frequency integral. The link under test was composed of 10 × 100 km of single-mode fibers, with full equalization of the SRS gain at the end of each span. The fibers had SPM nonlinear coefficient and chromatic dispersion coefficient reported in Fig. 6. It is worth noting that the nonlinear coefficient has an opposite trend with respect to wavelength compared to the chromatic dispersion. Comparing Fig. 6 with the results reported in Figs. 3-5 we also note that the nonlinear coefficient trend is Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. opposite to the one of the SRS gain. Throughout this section, we consider the transmission of a dual polarization WDM multiband signal, where each channel was modulated with complex Gaussian distributed symbols. The channel symbol rate was 64 Gbaud and the spacing 75 GHz, with guard bands of 500 GHz between neighboring multi-bands.
We first tested the accuracy of the SPM and XPM variance closed-form expressions with respect to the channel power. Fig. 7 shows such a variance for a flat WDM launch power over the E+S+C+L band, normalized to the cube of the channel power. The results are reported for four channels of interest placed in the middle of the E, S, C, and L bands, and plotted versus the channel power. The markers represent the benchmark while the solid lines are the formulas in (29)- (30). We can see from the figures that the variance does not scale with P 3 , contrary to what happens in absence of SRS [8]. This is caused by the SRS gain which makes the link kernel in (24) dependent on the WDM total power. In particular, we note that the SPM and XPM variance of the CUT in the central E band scales slower than P 3 , while the scaling is faster in the L band. Moreover, the inset of Fig. 7 shows that, due to SRS, the balance between the two nonlinear effects depends both on the channel power and the band. The worst RMSE between the XPM variance obtained with the proposed theory and the numerical results is 0.5 dB for the CUT in the S-band, while in the other bands we obtained an RMSE of 0.2 dB. On the other hand, the estimation of SPM is more critical, with a worst-case RMSE of 1.3 dB in the L band. However, since XPM is the dominant impairment in the transmission under test, combining the two nonlinear effects yields an RMSE of 0.5 dB for the CUT in the S-band and 0.3 dB for the others.
Then, we considered the transmission of signals spanning the C+L, S+C+L, and E+S+C+L bands with fixed channel power equal to -1 dBm. The results are shown in Fig. 8, where the markers indicate the benchmark GN model, and the solid lines are the generalized closed-form expressions in (29)- (30), while the dashed lines indicate the closed-form expressions derived in [19] for the C+L band assuming a linear Raman gain. For a fair comparison, in the latter expressions we used the same per-channel attenuation and wavelength-dependent fiber nonlinear coefficient γ as in (29)-(30) of this work. Fig. 8 (left) shows that in the C+L band the novel variance expressions coincide with those used in [19], and accurately estimate the NLI variance since the condition of case 1 (see Fig. 2) is always met. We note that the NLI variance is almost flat with respect to wavelength, contrary to what was observed in [24]. This is due to the wavelength dependency of the nonlinear coefficient γ reported in Fig. 6 (which was not accounted for in [24]) that partially counterbalances the SRS-induced tilt.
On the other hand, when the WDM signal covers the S+C+L bands as in Fig. 8 (center), the effect of the distributed interaction between the SRS and the Kerr effect is more visible, inducing a tilt on the NLI variance at the expense of the channels populating the L band. In this scenario, the generalized formulas of this work are in better agreement with the numerical estimation, with an RMSE of 0.4 dB against the 0.7 dB with a linear Raman gain coefficient (the dashed line). This result is consistent with comments drawn in Section II-C on the overestimation of the SRS gain in the WDM edges when a linear approximation is used. The linear approximation of the Raman gain becomes highly detrimental over wider bandwidths. In particular, for the E+S+C+L transmission in Fig. 8 (right), the proposed formulas allowed us to reduce the RMSE from 2 dB to 0.6 dB. It is worth noting that, in this last scenario, the NLI variance exhibits a tilted profile vs wavelength in the opposite direction compared to the one reported for the S+C+L band transmission. Such an NLI variance inflation in the E band is caused by the low values of the chromatic dispersion coefficient in this band [30], whose impact on the NLI out-weights the one of SRS.
For the sake of completeness, in extra simulations, we computed the numerical GN model imposing the same triangular profile and attenuation value used in the theoretical formulas, and we obtained RMSE values comparable to those reported for Fig. 8. Such results suggest that the main source of model inaccuracy results from the assumption that the integrated power within a sub-band is undepleted by SRS, i.e., (14)- (15).
Finally, we tested the proposed formulas in the presence of a signal power pre-emphasis, hence with the extra approximation of (22). We focused on an S+C+L band transmission with a fixed total power of 23 dBm. We varied the amount of pre-compensated SRS gain by consideringk = 0, 0.2, and 0.5. All he other parameters were the same as Figs. 7 and 8. Fig. 9 shows the benchmark results through markers and the theoretical results with solid lines. The different power allocations are reported in the inset. Such a pre-emphasis is designed to pre-compensate all the accumulated SRS gain on the signal power whenk = 1, which does not coincide with the SRS gain experienced by NLI variance since the NLI arises in a distributed way along propagation. In fact, the figure shows that the signal power pre-emphasis almost counterbalances the SRS-induced tilt on the NLI variance whenk = 0.2 while it induces an opposite-sign tilt whenk = 0.5. Although using the power profile derived in Section II-B with a non-uniform power  allocation represents an approximation, the estimated RMSE between the NLI variance numerical results and the proposed formulas is within 0.5 dB in all three cases.

IV. SNR ESTIMATION
By treating the ASE and the NLI as statistically independent additive noises, the SNR of a generic WDM channel, after SRS gain equalization on the signal power at the receiver, can be written as: where σ 2 NLI is the channel NLI variance and σ 2 ASE is the ASE noise variance. The variance of the ASE noise at the end of a link composed of N identical spans can be generalized as: where h is the Planck's constant, G(f i ) and F (f i ) are the gain and noise figure of each line amplifier + DGE. The gain of the equalized amplifier restores the target signal PSD by equalizing the SRS and the fiber loss. An ultra-fast theoretical estimation of the SNR in (33) can be obtained by exploiting the novel ECZ SRS model for the evaluation of both the ASE and the NLI variance. To test the accuracy of such an estimation, we compared the theoretical results with a fully-numerical SNR estimation relying on a Runge-Kutta solution of the ODEs in (1) for both the NLI and the ASE noise.
At first, we focused on an E+S+C+L transmission with total power equal to 25 dBm. All the system parameters were the same as in Section III. Each span had identical amplifiers with a per-band noise figure [30] equal to 6 dB (E band), 7 dB (S band), 5.5 dB (C band), and 6 dB (L band). Fig. 10 shows the SNR versus wavelength, with a good agreement between the theoretical SNR (line) and the exact SNR (markers). In the inset, we report the variances of the NLI and the ASE noise, showing that the E band is the most affected by both noise sources with particularly strong ASE noise, hence justifying the poor performance of such channels in terms of SNR.
To further explore the reliability of the proposed model in predicting the system SNR, we ran a campaign of simulations for many system setups. Namely, we randomly picked some of the link and WDM parameters, and we estimated the SNR for 1000 setup realizations. In particular, we picked the number of spans between 1 and 20, the chromatic dispersion between 8 and 18 ps/(nm · km), the symbol rate R between 32 and 96 Gbaud, and the channel spacing between 1.2×R and 1.7×R, the power P level between −4 and 0 dBm, and the pre-emphasis factork between 0 and 1. All the parameters had uniform Fig. 11. Wavelength RMSE on the SNR vs. the total WDM power and preemphasis factork, for 500 random setups in the E+S+C+L band (blue dots) and 500 in the S+C+L band (red circles).
distributions. We collected 500 random seeds for an E+S+C+L band transmission and 500 for S+C+L band and evaluated the wavelength RMSE of the theoretical SNR. The analysis revealed that the most challenging scenarios for the model accuracy are those associated with high power. Fig. 11 shows the RMSE versus total power (left) and versus pre-emphasis (right). The figure shows an increasing trend of the RMSE with the total power, both in the S+C+L band (red circles) and in the E+S+C+L band (blue dots). We note that, for the E+S+C+L band setups, the high RMSE points observed in the left figure are associated with a smallk. This is because a smallk increases the SRS gain tilt (see Fig. 5), such that some channels experience too high power over the E+S+C+L bands, thus outside the perturbative assumptions of the model. The pre-emphasis helps mitigate such a problem.
We next exploited the estimation of the SNR to evaluate the achievable information rate (AIR) for various transmissions. Treating both the signal and the NLI as Gaussian distributed, the AIR is given by the Shannon formula: Despite the outer summation over the number of channels, the presence of fiber nonlinear effects makes the growth of the AIR with the WDM bandwidth slower than linear. In addition, the presence of SRS contributes to setting a strong dependency of the SNR on the frequency. To investigate the scaling of the AIR with the transmitted bandwidth, we considered a gradual upgrade from the C band to the E+S+C+L band. For a fair comparison, in each case we adopted a power allocation (in the form of signal power preemphasis) maximizing the AIR. All the other parameters were identical to those used in Section II-C and Section III. Fig. 12 shows the spectral efficiency SE = 2log 2 (1 + SNR(f )) for the four transmissions under test and, in the inset, the corresponding AIR vs the WDM bandwidth. It can be seen that upgrading from C to C+L band, corresponding to a 65% increase in the bandwidth, yields a similar growth in terms of AIR (63%). Such a bandwidth upgrade penalizes the C band, which acts as a pump to the L band through SRS. Moving towards the S-band, we note that while its channels boost the performance of the C band and, partially, the L band, they offer lower SE. Overall, the 45% bandwidth increase results in 39% of AIR increase. The same phenomenon is observed for the upgrade toward the E band, which offers a limited contribution of the E-band per se, however with positive implications on the other bands thanks to its SRS pumping effect, as shown in Fig. 12. In this case, increasing the bandwidth by 42% with respect to the S+C+L transmission yields an AIR that is 31% higher.

V. CONCLUSION
Fast theoretical models for QoT estimation help alleviate the complex task of the design and the study of multi-band optical systems. In this work, we generalized the Christodoulides-Zirngibl SRS model of the signal power evolution along a fiber optic by leveraging a truncated triangular approximation of the Raman gain coefficient. The numerical results show that the proposed ECZ SRS formula is a viable tool to estimate in real-time the SRS-affected signal power, thus effectively mitigating the SRS overestimation of the CZ SRS formula in transmissions covering the S+C+L and E+S+C+L bands.
We showed that, by using the ECZ SRS formula, the GN model closed-form expressions of [19], originally derived with a linear Raman gain coefficient, based on the CZ model, can be extended to estimate the variance of SPM and XPM in the context of multi-band transmissions. We reported a significant reduction of the RMSE compared to the original GN model in [19] over the S+C+L and E+S+C+L bands.
The proposed analytical expression for the SRS gain can be used to quickly estimate even the SRS-affected ASE variance, and thus both the SNR and the AIR under the GN model assumptions. A random simulative campaign over 1000 different setups showed that, except for a few outliers, the proposed model estimates the SNR with an RMSE smaller than 0.5 dB up to a total power of 23 dBm over the E+S+C+L band. It is worth noting that using the GN model of [19] based on the linear approximation of the Raman gain, the RMSE diverges more rapidly with the total power, yielding an RMSE greater than 1 dB already with 19 dBm in the E+S+C+L band.
Finally, we exploited the analytical model to investigate the benefits in terms of AIR in upgrading the system bandwidth from the C to the E+S+C+L band. The analysis showed that the optical channels in the E band offer limited contributions to the AIR while boosting the performance on the S, C, and L bands through SRS, in agreement with the findings of [22].

APPENDIX A EVALUATION OF THE RAMAN SHAPING PROFILE
In this section we report the steps adopted in the evaluation of the following expression: P I (0; ν 1 (f ), ν 2 (f )) − P I,T (0; f, Δf co ) df (36) appearing in (17). Assuming that the transmitted PSD is flat, the integrated power can be expressed as: where the effective window edges ν 1 (f ) and ν 2 (f ), defined in (9)-(8), identify the four scenarios sketched in Fig. 2. Similarly, the integrated power in the Raman window with trapezoidal quadrature, defined in (13), can be written as: where we used the Heaviside function H(x) to set to zero the PSD outside the effective window. Note that in case 2 P I,T ≡ P I in the case of flat PSD.
We then evaluate the indefinite integral of the integrated power expressions in (37)-(38) by exploiting the properties: i) H(x)dx = xH(x), and ii) H(−x)dx = x(1 − H(x)). Finally, gluing the solutions at the domain edges we find the Raman shaping factor r(f ) in (18). For instance, the solution of case 3 has to coincide with 1 when the effective window upper edge ν 2 (f )=f + Δf co reaches the WDM spectrum highest frequency f M , i.e., f = f M − Δf co .