Introduction

Optical fibre propagation can be modelled by the nonlinear Schrödinger equation (NLSE)1. The NLSE is a partial differential equation that has three main effects: attenuation, second order dispersion, and Kerr nonlinearity. When considering the three effects together, the NLSE has no analytical solution for an arbitrary input pulse. The solution is obtained via a numerical method known as split step Fourier method (SSFM)1,2. The SSFM divides the fibre in small segments, known as steps, and increasing the number of steps increases the method’s accuracy2,3. For a high number of steps, the method’s computational complexity is a limiting factor4. To overcome this limitation, analytical models are generally used.

An analytical model that is easily mathematically manipulable is highly desirable as it can be used for improving fibre-optical transmission. As models are approximations of the NLSE, they can be used to build improved receivers and mitigate fibre effects5,6,7,8, design signal shaping and coding9,10, to estimate channel capacity10,11,12,13,14, and even to predict system performance15,16. An extensive review on optical channel models can be found in ref. 10. Each model will have a regime of operation based on the approximation used to derive it. The regimes classify models with respect to the group-velocity dispersion parameter β2 (which in this paper we refer to as the linear coefficient), the Kerr nonlinear coefficient γ, and the input power.

Some of the main regimes of operation and models present in the literature are schematically represented in Fig. 1. If both the linear and nonlinear coefficients are zero, the fibre degenerates to an additive Gaussian noise (AWGN) channel (in the presence of noise), and no interesting effects of the fibre propagation appear. One of the simplest regimes that accounts for fibre propagation effects is when the linear coefficient is zero (β2 = 0), represented by the green region of Fig. 1. In this case, the NLSE is modelled by a nonlinear phase shift1, known as the nonlinear phase noise (NLPN) model11,17. The assumption of a zero linear coefficient results in a memoryless channel. This assumption was used in the literature, for example, when chromatic dispersion is completely compensated18. The same premise was considered when analysing the highly nonlinear regime, i.e., when the nonlinearities are predominant and may be an important effect to take into account12,19. As will be seen in the Results section, even in that regime the dispersion can deteriorate the performance of the model.

Fig. 1: Validity region for different models for the nonlinear Schrödinger equation.
figure 1

Each region is characterized by a combination of β2 and γ values. The models are derived by using approximations based on the magnitude of these two parameters.

Another simple model arises when the nonlinear coefficient is zero (γ = 0), represented by the blue region of Fig. 1. In this case, the NLSE again admits an exact solution given by the so-called dispersion-only model1. This model considers the fibre propagation as an all-pass filter, whose phase response grows with the square of the frequency. As the model considers zero nonlinearities, it is ideal for low power regimes, where the dispersion is the major effect.

For high power regimes, if the nonlinear coefficient is low but nonzero, the dispersion-only model becomes inaccurate, as will be shown in the Results section. In such scenario, regular perturbation (RP) theory on the nonlinear coefficient20,21,22 becomes a more suitable model, represented by the yellow region in Fig. 1. The nonlinearities depend on the signal times the square of the absolute value of the signal, as will be seen in the nonlinear term of Eq. (3). This dependence makes the nonlinearities grow with the cube of the signal power, thus compromising the accuracy of the RP for high powers20,21. As will be shown in the Results section, RP on the nonlinear coefficient is accurate for a wider range of powers than the dispersion-only model. This wider range allows the RP model to model many communication systems8,22. The assumption of low nonlinearity is also used for other models, such as logarithmic perturbation21 and Volterra series23. With respect to the Volterra series model, it was proved in ref. 20 that its (2n + 1)-th order is equivalent to the n-th order RP on the nonlinear coefficient model.

The mentioned models only cover regions , , and of Fig. 1. Models for the red region and the brown region do not exist in the published literature. The latter represents regimes where both the linear and nonlinear coefficient are high, and might not be achievable by perturbation models as, by the RP’s definition, one of the coefficients should be low. In region , the absolute value of the linear coefficient is low; therefore, performing a RP is a reasonable approach.

In this paper, we propose a perturbation on the linear coefficient of the NLSE, providing a model for the weak-dispersion regime represented by region . The proposed RP on the linear coefficient covers regimes where the nonlinear coefficient is high, and in contrast to the NLPN model, small amounts of dispersion are allowed. The RP on the linear coefficient is a model in closed mathematical form, and depends on the input field and the fibre parameters. A closed-form equation is derived for the continuous-time fibre output with single polarization. The proposed model is compared with the RP on the nonlinear coefficient, the dispersion only, and the NLPN models. For comparison purposes, the fibre and simulation parameters (such as bandwidth and span length) are varied to identify regimes where each model is accurate. As will be shown in this paper, the RP on the linear coefficient is accurate for a wider range of powers than the RP on the nonlinear coefficient in low accumulated dispersion systems.

Results

Fibre propagation and metrics

The noiseless propagation of the optical field E at the retarded time frame t and distance z for a single polarization in a single-mode fibre can be represented by the NLSE1

$$\frac{\partial E(t,z)}{\partial z}=-\frac{\alpha }{2}E(t,z)-\frac{j{\beta }_{2}}{2}\frac{{\partial }^{2}E(t,z)}{\partial {t}^{2}}+j\gamma | E(t,z){| }^{2}E(t,z),$$
(1)

where α is the attenuation coefficient, β2 the group-velocity dispersion parameter, and γ the nonlinear coefficient. When normalizing the field E via

$$E(t,z)=A(t,z){{\rm{e}}}^{-\frac{\alpha }{2}z},$$
(2)

Eq. (1) is simplified to

$$\dfrac{\partial A(t,z)}{\partial z} = {\underbrace{-\dfrac{j\beta_2}{2}\dfrac{\partial^2 A(t,z)}{\partial t^2 }}_{{\text{linear}} \, \, {\text{term}}}} + \underbrace{j\gamma {\mathrm{e}}^{-\alpha z}|A(t,z)|^2 A(t,z)}_{{\text{nonlinear}} \, \, {\text{term}}}.$$
(3)

By writing the NLSE as Eq. (3), the RP method can be more easily applied.

This version of the NLSE disregards dual-polarization effects, such as polarization-mode dispersion, as well as other effects like stimulated Raman scattering. The latter can be reasonably ignored given the low-bandwidth scenarios studied in this paper1. We chose to analyse the single-polarization equation as a first step toward RP on β2. Nevertheless, single-polarization transmission is still attractive for low cost optical systems24.

The first and second terms on the right-hand side of Eq. (3) are linear and nonlinear terms in the NLSE. In this paper, for simplicity, we will refer to them as linear and nonlinear terms, respectively, even though they refer to a normalized version of Eq. (1). When both the linear and nonlinear terms are considered together, there is no analytical solution for Eq. (1) for an arbitrary input pulse A(, 0). However, by setting β2 = 0 or γ = 0, it is possible to obtain simple analytical solutions1. These solutions are the basis of the models described below, and each of them has a regime where they can predict well the NLSE solution of Eq. (3).

In this paper, to quantify how well a model approximates the solution of the NLSE in Eq. (3), we will use a metric that relates the error between two waveforms: the normalized square deviation (NSD) (previously used in ref. 20). The output of the SSFM algorithm will be considered as the solution A of Eq. (3). The NSD calculation between A(, z) and its approximation made by a certain model AM(, z) is illustrated in Fig. 2. For a certain propagation distance z > 0, the error ξ between the model and the fibre output A is

$$\xi (t,z)={A}_{M}(t,z)-A(t,z).$$
(4)

Based on Eq. (4), we define the NSD as

$${\text{NSD}}\,\triangleq \frac{{\int }_{\!\! -\infty }^{\infty }| \xi (t,z){| }^{2}\,{\text{d}}t}{{\int }_{\!\! -\infty }^{\infty }| A(t,z){| }^{2}\,{\text{d}}t}.$$
(5)

The NSD captures the average of the squared absolute error over the time dimension, which corresponds to the error energy. To enable a fair comparison between different input powers, the NSD is normalized by the energy of the fibre output A(, z). Following ref. 21, we will use a threshold of 0.1% for comparing models. We will say that a model is precise if it gives an NSD below 0.1%.

Fig. 2: Schematic of NSD and SNR calculations.
figure 2

The transmitted symbols (\({x_1},\ldots,{x_K}\)) are filtered by a root-raised-cosine (RRC) filter, originating the input signal A(, 0). This input signal is processed either by the SSFM or by an analytical model (AM), resulting in the outputs A(, z) and AM(, z), respectively. The NSD is obtained from the SSFM output A(, z) and a model output AM(, z), whereas the SNR is obtained from the transmitted symbols (\({x_1},\ldots,{x_K}\)) and with the received symbols (\({y_1},\ldots,{y_K}\)) from the SSFM or from an analytical model. To obtain the received symbols, the output waveforms are submitted to an optional chromatic dispersion compensation (CDC) block, followed by a matched filter (MF) and sampling.

In addition to the NSD, which characterizes the continuous-time performance, we also observe the discrete-time output, for which the signal-to-noise ratio (SNR) is a suitable metric. As shown in Fig. 2, we transmit a sequence of K symbols (x1, , xK), shaped by a root-raised-cosine (RRC) filter. The receiver consists of an optional chromatic dispersion compensation (CDC) block, followed by a matched filter and sampling. Even though no noise is added to the system, the received symbols are not exactly the transmitted ones. This difference is owing to limitations on the linear receiver we assume, which cannot undo the fibre propagation effects on the signal. Therefore, the SNR only accounts for signal–signal interactions in this case.

The SNR for a constellation with M symbols is defined as

$${\text{SNR}}\,\triangleq\frac{{\sum }_{m=1}^{M}| {\overline{y}}_{m}{| }^{2}}{{\sum }_{m=1}^{M}\frac{1}{{N}_{m}}{\sum }_{k=1}^{{N}_{m}}| {y}_{km}-{\overline{y}}_{m}{| }^{2}},$$
(6)

where \({\overline{y}}_{m}=\frac{1}{{N}_{m}}{\sum }_{k=1}^{{N}_{m}}{y}_{km}\)is the average received symbol corresponding to the m-th constellation point, Nm is the number of times the m-th constellation point was transmitted, and ykm is the k-th received symbol given a fixed transmitted m-th constellation point. This SNR calculation assumes that we know the corresponding transmitted symbol for a given received symbol, and that the mean values \({\overline{y}}_{m}\) would be the signal components with ideal reception.

Notation convention: throughout this paper, we use A(, z) to represent the evaluation of a two-variable function A that depends on the retarded time frame, evaluated at distance z. In other words, we use this notation to emphasize that A(, z) is still a function of the retarded time frame. The complex conjugate of A is denoted by A*. {} and {} give the real and imaginary parts of a complex number, respectively. Operators are denoted by calligraphic letters.

The numerical examples we will present investigate the limits and the operational regimes of each model. To vary the group-velocity dispersion, two types of fibre were considered: standard single-mode fibre (SSMF) and nonzero dispersion-shifted fibre (NZDSF). The SSMF has α = 0.2 dB km−1, β2 = − 21.67 ps2 km−1, and γ = 1.2 W−1 km−1, whereas the considered NZDSF has α = 0.22 dB km−1, β2 = − 5.42 ps2 km−1, and γ = 1.46 W−1 km−1. Except in the symbol rate variation analysis, all the simulations consider a symbol rate of 10 Gbaud. The modulation format is 64-ary quadrature amplitude modulation (64-QAM) unless otherwise stated.

In what follows, we review 3 models available in the literature and then, present the RP on β2.

Dispersion-only model

When considering γ = 0, Eq. (3) simplifies to1

$$\frac{\partial {A}_{M}(t,z)}{\partial z} \, = \,-\frac{j{\beta }_{2}}{2}\frac{{\partial }^{2}{A}_{M}(t,z)}{\partial {t}^{2}},$$
(7)

which has the exact solution25

$${A}_{M}(t,z)=\left(A(\cdot ,0)* h(\cdot ,z)\right)(t)={{\mathcal{D}}}_{z}\{A(\cdot ,0)\}(t).$$
(8)

In Eq. (8), * represents convolution, h is given by

$$h(t,z)=\frac{1}{\sqrt{j2\pi {\beta }_{2}z}}{{\rm{e}}}^{-\frac{j}{2{\beta }_{2}z}{t}^{2}},$$
(9)

and \({{\mathcal{D}}}_{z}\) is the dispersion operator defined as

$${{\mathcal{D}}}_{z}\{f\}(t) \, \triangleq \, \left(f* h(\cdot ,z)\right)(t),$$
(10)

where f is a function of t. The solution in Eq. (8) is called the dispersion-only model, and is a linear, time-invariant all-pass filter. It corresponds to in Fig. 1.

Example 1. Figure 3a shows the NSD vs. input power for the dispersion-only model for the SSMF and an NZDSF. As shown in Fig. 3a, the model is an accurate approximation of the fibre output in this system for powers lower than  + 2 dBm for the 10 km fibre. This is as a regime where the nonlinearities are not predominant. As the input power increases, the NSD increases by ~ 2 dB/dBm. This figure also shows that a change in β2 of approximately four times does not considerably change the NSD. By increasing the distance from 10 to 80 km, the NSD grows almost one order of magnitude, which can be justified by the extended interaction between the nonlinearities and the chromatic dispersion, not modelled in this solution.

Fig. 3: NSD for the four described models as a function of the input power.
figure 3

Two fibre lengths (10 and 80 km) and two types of fibre (NZDSF and SSMF). a Dispersion-only model given by (8). b NLPN model given by (12). c RP on γ model given by (15). d RP on β2 model given by (19). Each model has a different accuracy, measured in NSD, for a specific set of input powers, β2 coefficient, and distance. The black horizontal line represents a constant 0.1% NSD, which is used for comparison between models.

NLPN model

For β2 = 0, Eq. (3) simplifies to1

$$\frac{\partial {A}_{M}(t,z)}{\partial z}=j\gamma {{\rm{e}}}^{-\alpha z}| {A}_{M}(t,z){| }^{2}{A}_{M}(t,z),$$
(11)

which has the exact solution

$${A}_{M}(t,z)=A(t,0){{\rm{e}}}^{j\gamma | A(t,0){| }^{2}G(z)},$$
(12)

where

$$G(z)={\int }_{\!\!\!\! 0}^{z}{{\rm{e}}}^{-\alpha u}\,{\text{d}}u=\frac{1-{{\rm{e}}}^{-\alpha z}}{\alpha }$$
(13)

is the effective length. The solution in Eq. (12) is called the NLPN model, and is a memoryless, signal-dependent phase shift. This model is in Fig. 1.

Example 2. Figure 3b shows the NSD of the NLPN model for the same scenario as Example 1. For the NLPN, increasing the value of β2 (changing from NZDSF to SSMF) can deteriorate the NSD by a factor of almost 10 times for the same distance at 0 dBm. This effect can be justified by the increased dispersion contribution to the signal, not modelled by this solution. Furthermore, the accumulated dispersion is also increased when extending the fibre length from 10 to 80 km, which increases the NSD by two orders of magnitude. This last increase can also be seen by the fact that, at long distances, \(G(z) \approx 1/ \alpha\) and Eq. (12) barely changes with respect to z. However, this result is not consistent with the actual fibre output, which still depends on the distance z in the presence of chromatic dispersion. Although only the 10 km NZDSF presented an NSD below the threshold of 0.1%, the NLPN model has an NSD almost constant for powers lower than 10 dBm.

RP on γ

The RP method consists of representing the solution of an equation with an expansion in terms of simplified solutions and a coefficient. In the RP on γ, the general solution can be written as20,22

$$A(t,z)=\mathop{\sum }_{k=0}^{\infty }{\gamma }^{k}{A}_{k}(t,z),$$
(14)

where the functions Ak are functions that depend on the initial field A(, 0). In the first order RP on γ, Eq. (14) is truncated at k = 1, and the functions A0 and A1 are given in the next theorem. The result of Theorem 1 is well known in the literature (see ref. 20,21,22). However, we include its proof in Supplementary Note 2 for consistency with the notation.

Theorem 1. Let A be the solution of the NLSE in Eq. (3) with initial condition A(, 0). Then, A can be approximated by AM, the first order RP on the nonlinear coefficient γ of Eq. (3), written as

$${A}_{M}(t,z)={A}_{0}(t,z)+\gamma {A}_{1}(t,z),$$
(15)

where

$${A}_{0}(t,z)={{\mathcal{D}}}_{z}\{A(\cdot ,0)\}(t)$$
(16)

is the dispersion-only solution in Eq. (8) and

$${A}_{1}(t,z)=j{\int }_{\!\!\!\! 0}^{z}{{\rm{e}}}^{-\alpha u}{{\mathcal{D}}}_{z-u}\left\{| {A}_{0}(\cdot ,u){| }^{2}{A}_{0}(\cdot ,u)\right\}(t)\,{\text{d}}u.$$
(17)

The proof is given in Supplementary Note 2.

Example 3. In Fig. 3c, the NSD results for SSMF and NZDSF fibres are presented. This figure shows that the NSD scales as 4 dB/dBm, which is twice as much as the value found for the dispersion-only model example (see Fig. 3a). This behaviour might be explained by the cubic signal power dependence of the nonlinear term and A1 (see nonlinear term in Eq. (3)). Despite this faster growth of the NSD, the crossing point of the curves with the 0.1% threshold happens at higher powers compared with the dispersion-only case. This shows that the RP on γ model has a wider range of validity than the dispersion-only model. Figure 3c also shows that, for the RP on γ model, just like for the dispersion-only case, increasing β2 slightly improves the NSD. The increase in the fibre length also deteriorates the performance of the RP on γ. As expected, the NSD grows when the power of the input signal is increased. To reduce the approximation error, more terms of the expansion in Eq. (14) can be considered. This comes at the cost of a higher model complexity, as done in ref. 20.

Proposed model: RP on β 2

We now present the proposed model, which was derived based on the RP method to accurately represent the NLSE in the highly nonlinear regime, illustrated by region in Fig. 1.

In analogy with the RP on γ, the RP method can be applied to β2. The only difference is that now the expansion of A is written in terms of β2 as

$$A(t,z)=\mathop{\sum }\limits_{k=0}^{\infty }{\beta }_{2}^{k}{A}_{k}(t,z).$$
(18)

The terms A0 and A1 for the first order RP on β2 are given in the next theorem, which is the main contribution of this paper.

Theorem 2. Let A be the solution of the NLSE in Eq. (3) with initial condition A(, 0). Then, A can be approximated by AM, the first order RP on the linear coefficient β2 of Eq. (3), written as

$${A}_{M}(t,z)={A}_{0}(t,z)+{\beta }_{2}{A}_{1}(t,z),$$
(19)

where

$${A}_{0}(t,z)=A(t,0){{\rm{e}}}^{j\gamma | A(t,0){| }^{2}G(z)},$$
(20)

and

$${A}_{1}(t,z)=B(t,z){{\rm{e}}}^{j\gamma | A(t,0){| }^{2}G(z)},$$
(21)

with B given by

$$ B(t,z) = -M(t)z+ G_{1}(z) R(t)+G_{2}(z)P(t) \\ -2j\gamma A(t,0) \Re\{A^*(t,0) V(t,z)\} ,\qquad\qquad$$
(22)
$$V(t,z)= \ G(z) \left[M(t)z- G_{1}(z) R(t)-G_{2}(z)P(t) \right] \\ -G_{1}(z)M(t) +G_{2}(z)R(t)+G_{3}(z)P(t) ,$$
(23)
$$M(t)=\frac{j}{2}\frac{{\partial }^{2}A(t,0)}{\partial {t}^{2}}, \qquad \qquad \qquad$$
(24)
$$R(t)=\frac{\gamma }{2}A(t,0)\frac{{\partial }^{2}| A(t,0){| }^{2}}{\partial {t}^{2}}+\gamma \frac{\partial A(t,0)}{\partial t}\frac{\partial | A(t,0){| }^{2}}{\partial t},$$
(25)
$$P(t)=\frac{j{\gamma }^{2}}{2}A(t,0){\left(\frac{\partial | A(t,0){| }^{2}}{\partial t}\right)}^{2},$$
(26)
$${G}_{1}(z)=\frac{\alpha z+{{\rm{e}}}^{-\alpha z}-1}{{\alpha }^{2}}, \qquad \qquad \, \, \,\,$$
(27)
$${G}_{2}(z)=\frac{2\alpha z+4{{\rm{e}}}^{-\alpha z}-{{\rm{e}}}^{-2\alpha z}-3}{2{\alpha }^{3}},\, \, \,$$
(28)
$${G}_{3}(z)=\frac{6\alpha z+18{{\rm{e}}}^{-\alpha z}-9{{\rm{e}}}^{-2\alpha z}+2{{\rm{e}}}^{-3\alpha z}-11}{6{\alpha }^{4}}. \quad \qquad$$
(29)

The proof is postponed to Supplementary Note 3.

By comparing Eq. (19) and Eq. (20) with Eq. (12), it is clear that the RP on β2 corresponds to the NLPN solution perturbed by the dispersion. The perturbation term A1 in Eq. (21) is in closed form and it depends on the derivatives of the input field A(, 0). Note also that both A1 and A0 are multiplied by the same phase rotation, as shown in Eq. (20) and Eq. (21).

The functions A0 and A1 depend on elementary functions of z (see equations (13), (27), (28), and (29)). This results in the same calculation time for AM at any distance z, in contrast to the SSFM, for which the number of necessary steps increases with the distance. For this reason, RP on β2 calculations are significantly faster than SSFM calculations for the parameters simulated in this paper.

Example 4. Figure 3d presents the NSD for the RP on β2 for the SSMF and the NZDSF. As shown in Fig. 3d, the NSD is approximately constant below a certain input power (0 dBm in this case). In addition, changes in the amount of dispersion severely impact the NSD. For example, at 0 dBm, increasing β2 (going from NZDSF to SSMF) in the 80 km system, the NSD rises more than two orders of magnitude. This model is also very sensitive to the distance. Increasing the fibre length from 10 to 80 km in the SSMF system, the NSD rises more than three orders of magnitude for 0 dBm.

As discussed above, the RP on β2 in Eq. (19) is the sum of the NLPN solution (A0) and a term accounting for the dispersion (β2A1). Therefore, some similarities in the NSD curves of these two models are expected (in analogy to the RP on γ and the dispersion-only model, discussed in Example 3). By comparing Fig. 3d with Fig. 3b, the RP on β2 is an improved version of the NLPN model, just like the RP on γ is an improved version of the dispersion-only model.

Example 5. Figure 4a presents a comparison of the four models for an 80 km NZDSF. As shown in Fig. 4a, only the dispersion-only, the RP on γ, and the RP on β2 models have NSDs below the threshold of 0.1% for some powers. The dispersion-only model crosses the threshold at  − 2 dBm, whereas the RP on γ crosses it at 6.2 dBm, presenting a gain of 8.2 dB. The RP on β2 crosses the threshold at even higher powers (9.2 dBm), with a gain of 3 dB with respect to the RP on γ.

Fig. 4: NSD for the four models in 80 km of an NZDSF.
figure 4

a With attenuation (α = 0.22 dB km−1). b Without attenuation (α = 0 dB km−1). The power gap in dB between the NSD curves is related to different input powers that allow each model to achieve 0.1% NSD. In the system without attenuation b, this gap is increased when comparing RP on γ with RP on β2.

In systems with distributed Raman amplification, the power profile is approximately flat26 and the attenuation term in (1) is often neglected. In this case, a simpler analytical form for B(tz) when compared with Eq. (22) is achieved. This simplification is given in the next theorem.

Theorem 3. With ideal distributed amplification, the functions A0 and A1 in Eq. (20) and Eq. (21) can be written as

$${A}_{0}(t,z)=A(t,0){{\rm{e}}}^{j\gamma | A(t,0){| }^{2}z},$$
(30)
$${A}_{1}(t,z)=B(t,z){{\rm{e}}}^{j\gamma | A(t,0){| }^{2}z},$$
(31)

where

$$B(t,z) = -M(t)z + \dfrac{z^2}{2} R(t) +\dfrac{z^3}{3}P(t) \\ -2j\gamma A(t,0) \Re\left\{ A^*(t,0) \left[\dfrac{z^2}{2}M(t)-\dfrac{z^3}{6}R(t)-\dfrac{z^4}{12}P(t)\right]\right\}$$
(32)

and where M(t), R(t), and P(t) are given by Eq. (24), Eq. (25), and Eq. (26), respectively. The proof is postponed to Supplementary Note 4.

Example 6. Figure 4b shows the NSD for the four models with the same parameters as in Example 5, except that in this case there is no attenuation. For powers below  − 5 dBm, the NLPN model and the RP on β2 present almost the same NSD compared with Fig. 4a. This behaviour could be justified by the small impact of the nonlinearities for low powers. The NSD values in that regime become close to each other, since the attenuation is mostly connected to the nonlinear effect (as can be seen in Eq. (3)). The curves cross the threshold at lower powers compared with Fig. 4a, excluding the NLPN model, which remains above the threshold for all analysed powers. We believe that the lower threshold crossing happens owing to the additional interactions between nonlinearities and dispersion present in the lack of attenuation. Although the RP on γ crosses the threshold at 0 dBm, the RP on β2 crosses at 5 dBm, representing a gain of 5 dB (2 dB more compared with the attenuation case).

In the previous examples, the parameters γ and β2 were fixed, along with the simulation bandwidth and the fibre length. These parameters are further investigated in the next examples, where the four models are compared with each other in systems with attenuation. For the next simulations, most of the parameters given in the previous examples are still considered; however, some of them will be changed. We will discuss fibre length variation, symbol rate variation, and β2 and γ variation.

Variation of β 2 and γ

To analyse the impact of changes in β2 and γ, we consider a fibre length of 80 km with fixed input power of 5 dBm, symbol rate of 10 Gbaud, and α = 0.2 dB km−1. In this analysis, only the RP on γ and β2 will be considered. Figure 5a shows the NSD for the RP on γ for different values of β2 (negative β2) and γ. As depicted in Fig. 5a, variations in β2 practically do not affect the accuracy of the model, whereas changes in γ have a major impact. In analogy with Fig. 1, the lower region of Fig. 5a is equivalent to region , where the model is accurate (NSD < 0.1%). This region covers values of γ of up to ~ 1.78 W−1 km−1 for this system. Making the same analysis for RP on β2 leads to Fig. 5b, which shows the NSD for the same range of γ and β2 values. In this case, the area where the model is accurate is vertical, in analogy to region in Fig. 1. The NSD for the RP on β2 changes mostly with the value of β2; however, changes in γ can also significantly affect the performance, specially for high β2 values. The intersection of the areas that are not accurate for any of the models is related to region in Fig. 1.

Fig. 5: NSD for the RP on γ and the RP on β2 for different values of β2 and γ.
figure 5

The system consists of an 80 km fibre with fixed input power 5 dBm, symbol rate of 10 Gbaud, α = 0.2 dB km−1, and negative β2. a Regular Perturbation on γ. b Regular Perturbation on β2. Region represents low γ and large β2 values (see Fig. 1). Region represents low β2 and large γ values. Region represents large γ and large β2 values. The thick black line of constant 0.1% NSD is used as a threshold for accuracy comparison. The results indicate, as expected, that RP on γ is accurate for region and RP on β2 is accurate for region . Both models yield low accuracy in region .

Fibre length variation

For the fibre length variation analysis, we consider an NZDSF with fixed input power of 5 dBm and a symbol rate of 10 Gbaud. Figure 6a shows the NSD versus fibre length for the four different models. All the four models increase the NSD when increasing the fibre length; however, the dispersion-only and the RP on γ seem to converge to a constant NSD value. This convergence is owing to the attenuation on the nonlinear term in Eq. (3): for high distances, the nonlinearities do not considerably affect the signal, and the major effect is the dispersion. As these two models fully predict the dispersion effect, they do not lose accuracy in that regime. The RP on β2, for this system, can reach ~ 120 km within the NSD threshold of 0.1%, and for fibre lengths lower than 90 km, the model outperforms the RP on γ.

Fig. 6: NSD versus fibre length and symbol rate for the four models.
figure 6

The system consists of an NZDSF with fixed input power 5 dBm. a NSD versus fibre length at a symbol rate of 10 Gbaud. b NSD versus symbol rate for a fibre length of 80 km. Changes in fibre length or symbol rate result in a wider range of NSD variation for RP on β2 than for RP on γ.

Symbol rate variation

For the symbol rate variation analysis, we consider 80 km of an NZDSF with fixed input power of 5 dBm. Figure 6b depicts the NSD variation with respect to the symbol rate. As shown in Fig. 6b, the dispersion-only model and the RP on γ do not change their accuracy when varying the symbol rate. On the other hand, the NLPN model and the RP on β2 drastically drop the NSD when decreasing the symbol rate. For bandwidths lower than 4 Gbaud, even the NLPN can outperform the RP on γ. This behaviour may be justified by the decreasing influence of the dispersion on the signal, since according to Eq. (7), higher frequencies are more affected by dispersion owing to their high second derivative.

Fibre length versus symbol rate

The previous two sections analysed the NSD by varying the fibre length and symbol rate separately. Both of these parameters influence the accumulated dispersion. Therefore, given a fibre length or a symbol rate, we can find the values for the other parameter in which RP on β2 is accurate27.

 Figure 7 depicts the NSD given a fibre length and a symbol rate for an NZDSF with fixed input power 5 dBm. As shown in Fig. 7, the model can accurately handle arbitrarily large fibre lengths if the symbol rate is small enough and vice versa. By fixing a fibre length of 80 km, the maximum symbol rate in which RP on β2 is still accurate is 12.55 Gbaud (see triangular marker in Fig. 7). If the distance is reduced to 20 km, the symbol rate can be increased until 27.38 Gbaud (see square marker in Fig. 7). These values show that by reducing the fibre length by a factor of 4, the symbol rate can be increased by a factor of ~ 2.18. This difference in scaling factors was already expected considering that the accumulated chromatic dispersion increases linearly with the distance and with the square of the signal bandwidth.

Fig. 7: NSD for RP on β2 for different values of fibre length and symbol rate.
figure 7

The system consists of a 80 km NZDSF with fixed input power 5 dBm. The triangular marker on the top of the figure is associated with the distance of 80 km and symbol rate of 12.55 Gbaud. The square marker on the bottom right is associated with the distance of 20 km and symbol rate of 27.38 Gbaud. Both markers lie on the thick black line representing a constant 0.1% NSD.

The thick solid line in Fig. 7 at an NSD of 0.1% can be seen as a conservative threshold. Choosing other metrics, such as SNR, might motivate different conclusions. For example, we will show in the discrete-time performance section (ahead), RP on β2 can have an SNR close to that of the SSFM, even though RP on β2 yields an NSD in the order of 30% in that scenario.

Modulation format impact

The previous simulations were based on 64-QAM. For different modulation formats, the transmitted signal statistics change. As signal statistics impact the nonlinear effect, the NSD curves can be different for other modulation formats27.

 Figure 8 illustrates the NSD versus input power for three modulation formats on 80 km of an NZDSF with a symbol rate of 10 Gbaud. As shown in Fig. 8, quadrature phase-shift keying (QPSK) yields lower NSD than the other two modulation formats. We believe this behaviour is justified by the QPSK’s high tolerance to nonlinearities, which reduces the error in the first-order RP approximation. The performance of star 8-QAM and 64-QAM are almost the same for the considered input powers. The gap between QPSK and 64-QAM for the 0.1% threshold crossing is 0.8 dB for RP on β2 and 1.1 dB for RP on γ. The higher gap for RP on γ might indicate that this model is more sensitive to changes in the modulation format than RP on β2.

Fig. 8: Impact of QPSK, star 8-QAM, and 64-QAM modulation formats on the NSD versus input power.
figure 8

The system consists of an 80 km NZDSF with 10 Gbaud symbol rate. The modulation format impact is measured by the power gap in dB between the crossing points of the NSD curves for each modulation format with the horizontal black curve representing a constant 0.1% NSD.

Discrete-time performance

As discussed in the “Fibre propagation and metrics” section, the analysis of the received symbols might bring insightful conclusions about the models. In order to clearly visualize the fibre effects on the received constellations, the SNR simulations were based on QPSK modulation over a 20 km NZDSF. At the receiver, we considered two cases: with and without CDC, followed by matched filtering and sampling (see Fig. 2). We emphasize that these operations at the receiver are applied to the continuous-time output of the models and SSFM, and we are not using discrete-time analytical models.

Figure 9a shows the SNR for SSFM, RP on γ, and RP on β2 at 10 Gbaud for both receivers. As depicted in Figure 9a, for input powers lower than 6 dBm, the SNR for the receiver with CDC is higher than the one without CDC. The latter converges to ~35.9 dB. This convergence could be explained by the uncompensated dispersion effect, which does not depend on the signal power (a linear effect). For input powers >8 dBm, the systems with and without CDC are equivalent in SNR performance. This behaviour might be explained by the predominance of the nonlinear effect at these input powers. For input powers >11 dBm, the SNR is higher for RP on γ than for SSFM and RP on β2. In this power regime, RP on β2 can be seen as a more accurate model for SNR calculations.

Fig. 9: SNR versus input power, received constellations, and their respective PDFs.
figure 9

On the SNR results, receivers with (w/) and without (w/o) CDC were considered for a QPSK modulation over a 20 km NZDSF. The constellations were obtained for 40 Gbaud with CDC at 16 dBm. a SNR results for the three models at a symbol rate of 10 Gbaud. b SNR results for the three models at 40 Gbaud. c Received RP on β2 constellation and respective PDF. d Received SSFM constellation and respective PDF. e Received RP on γ constellation and respective PDF. The arrows pointing to the markers in a and b represent which curves are overlapped in those parts of the figure. In c, d, and e, the respective SNR values are shown on top of each of them. For c and e, the NSD values are also displayed.

 Figure 9b illustrates the SNR for SSFM, RP on γ, and RP on β2 at 40 Gbaud for both receivers. At this symbol rate and input powers lower than 15 dBm, RP on β2 with CDC does not follow the SNR of SSFM and RP on γ. This difference might be related to the precarious estimation of the dispersion effect by the RP on β2 for this higher bandwidth scenario (see Fig. 7). With a poor estimation of the chromatic dispersion, the CDC will not compensate the exact dispersion effect predicted by the RP on β2 model. For the system without CDC, the performance of RP on β2 is close to the SSFM performance for all displayed input powers. The SNR for RP on γ for systems with and without CDC diverges again from the SSFM performance at input powers higher than 11 and 16 dBm, respectively. Meanwhile, at input powers higher than 16 dBm, the performance of RP on β2 for the CDC system approaches the SSFM performance.

For an input power of 16 dBm in the 40 Gbaud system with CDC, the received constellations for RP on β2, SSFM, and RP on γ are shown in Figs. 9c, d, and e, respectively, with their corresponding probability density functions (PDFs). As illustrated in these figures, the constellation for RP on β2 approximately preserves the circular shape of the SSFM’s constellation, which is not observed in the constellation for RP on γ. This preservation could be explained by the better prediction of nonlinear effects at this high input power by the RP on β2. The received symbols for RP on γ are mostly outside the unitary square given by y ≤ 1, which shows a high gain of energy when using this model. In addition, the PDF for RP on γ shows that its symbols are highly concentrated toward a single point, which was expected by its high SNR. Even though the SNR for RP on β2 is close to the SNR of SSFM (8.49 and 8.29 dB, respectively), the received constellation shapes are slightly different. For example, as observed in the constellation PDFs, the RP on β2 received symbols are more spread than the SSFM symbols. This contrast means that the SNR alone might not indicate precisely the accuracy of the model. On the other hand, a high NSD for RP on β2 (37.13%) may not indicate that the received signal in discrete time is severely different from the reference given by the SSFM, showing that our proposed model is accurate in scenarios where nonlinearities are the dominant effect.

Discussion

This paper presented a new closed form analytical approximation for the solution of the NLSE: the RP on β2. The derived approximation is a suitable model for low symbol rates, low fibre lengths, and/or high input powers. The RP on β2 was compared with three other models with respect to variations in the bandwidth, fibre length, input power, and fibre parameters.

The main comparison was with the RP on γ, a well-known model in the literature that is accurate in the regime of high dispersion and low nonlinearities. In a NZDSF of 80 km with attenuation, the RP on β2 can be used as an accurate model until input powers of 9 dBm, whereas the RP on γ is accurate only to input powers lower than 6 dBm. In addition, the RP on β2 is accurate for high γ values where the RP on γ is not. Thus, the new proposed model is convenient for the opposite regime, where the nonlinearities are predominant, and the dispersion has a minor effect.

As all simplified models, the proposed model has a limited range of validity. At the moment, the main applicability of the model is for applications that rely on low bandwidths (below 11 GHz) and short distances (below 80 km). This includes, for example, passive optical networks28. Another short-distance low-bandwidth application is hybrid fibre coax systems29. The model is in its present form not intended for long-haul or wavelength-division-multiplexed systems.

Another drawback of the proposed model is that it neglects the noise. This effect has been considered in the literature for the RP on γ in ref. 22, where noise was added in the zeroth order linear equation, followed by a Karhunen–Loève expansion to account for nonlinear signal–noise interactions. The same approach for the RP on β2 would lead to cumbersome equations, as the zeroth order equation is nonlinear in this case.

This paper represents the first steps in the theory of the proposed model. Possible extensions of this work are designing a receiver based on the RP on β2, deriving higher-order perturbations of the model, and adding noise within the RP analysis. The derivations were conducted for a single-polarization system, and the equations for dual-polarization are still a subject of further investigation. Separating the contribution of individual pulses and finding a discrete model are also left as future work. Although the focus of this paper was on optical-fibre communications, we believe that our model can be applied to other fields where the NLSE is applicable.

Methods

Simulation specifications

The simulations were conducted in Matlab® and considered 215 symbols randomly chosen from different constellations during the paper. The constellations were generated with unitary energy. To generate the constellation figures and PDFs in Fig. 9c, d, and e, we used 220 symbols for a smoother plot. The colour for each received symbol was attributed from a colormap according to the PDF values. The symbols were oversampled by 16 samples per symbol. After oversampling, the signal was shaped by an ideal RRC filter, with roll-off factor of 0.1, implemented in the frequency domain. After filtering, the signal was scaled to adjust to the desired input power. The resulting waveform was used as the input of either the SSFM or one of the models. For the SSFM case, we considered a symmetric SSFM implementation, with step-size 0.1 km. The step-size and the simulation bandwidth substantially impact the SSFM accuracy2,3. Further reducing the step-size and increasing the simulated bandwidth by increasing the numbers of samples per symbol did not impact the displayed results, which indicates 0.1 km step-size and 16 samples per symbol are accurate enough for the systems in this paper. For the RP on γ model, the integral of Eq. (17) was evaluated numerically using an integration step of also 0.1 km. For the RP on β2, the derivatives were obtained in the frequency domain. The dispersion operator was ideally implemented in the SSFM and in all simulated equations in the frequency domain. The CDC was implemented in the frequency domain with the dispersion operator applied with negative fibre length. Before the matched filter, the signal was scaled with the inverse scaling factor used in the transmission to adjust the input power. The matched filter was the same RRC filter as used in the transmission.